36 (u,v,w,p,t,uhx,uhy,uhz,vhx,vhy,vhz,whx,why,ff,ff0,gg,gg0,hh,hh0,qu,qv,qw)
40 real(mytype),
dimension(:,:,:),
intent(in) ::u,v,w,p,t,ff0,gg0, &
41 hh0,uhx,uhy,uhz,vhx,vhy,vhz,whx,why
42 real(mytype),
dimension(:,:,:),
intent(out)::ff,gg,hh,qu,qv,qw
46 real(mytype) :: uux, uvy, uwz, uxx, uyy, uzz, dpdx, ux
47 real(mytype) :: vux, vvy, vwz, vxx, vyy, vzz, dpdy, vx
48 real(mytype) :: wux, wvy, wwz, wxx, wyy, wzz, dpdz, wx
62 uux= ( uhx(i+1,j,k)**2.d0 - &
63 uhx(i,j,k)**2.d0 ) /
dx 65 uvy= ( uhy(i,j,k)*vhx(i,j,k) - &
66 uhy(i,j-1,k)*vhx(i,j-1,k) ) /
dy 68 ux = ( uhx(i+1,j,k)-uhx(i,j,k) ) /
dx 70 elseif (
cds .eq. 2)
then 72 uux= ( -uhx(i+2,j,k)**2.d0 + &
73 27.d0*uhx(i+1,j,k)**2.d0 - &
74 27.d0*uhx(i,j,k)**2.d0 + &
75 uhx(i-1,j,k)**2.d0 ) / 24.d0 /
dx 77 uvy= ( -uhy(i,j+1,k)*vhx(i,j+1,k) + &
78 27.d0*uhy(i,j,k) *vhx(i,j,k) - &
79 27.d0*uhy(i,j-1,k)*vhx(i,j-1,k) + &
80 uhy(i,j-2,k)*vhx(i,j-2,k) ) / 24.d0 /
dy 82 ux = ( -uhx(i+2,j,k) + &
83 27.d0*uhx(i+1,j,k) - &
85 uhx(i-1,j,k) ) / 24.d0 /
dx 90 uwz= ( uhz(i,j,k) *whx(i,j,k) - &
91 uhz(i,j,k-1)*whx(i,j,k-1) ) /
dz(k)
101 uyy= ( u(i,j+1,k) - &
105 elseif (
cds .eq. 2)
then 107 uxx= ( -u(i+2,j,k) + &
111 u(i-2,j,k) ) / 12.d0 /
dx2 113 uyy= ( -u(i,j+2,k) + &
117 u(i,j-2,k) ) / 12.d0 /
dy2 122 uzz= ( (u(i,j,k+1)-u(i,j,k))/
dz_t(k) - &
123 (u(i,j,k)-u(i,j,k-1))/
dz_b(k) ) /
dz(k)
130 ff(i,j,k) = -uux-uvy-uwz -
u_mrf*ux + &
131 fc*( 0.5d0*(vhx(i,j,k)+vhx(i,j-1,k)) -
vg )
141 if (
cds .eq. 1)
then 143 dpdx= ( p(i+1,j,k)-p(i,j,k) ) /
dx 145 elseif (
cds .eq. 2)
then 147 dpdx= ( -p(i+2,j,k) + &
150 p(i-1,j,k) ) / 24.d0 /
dx 155 qu(i,j,k)= u(i,j,k) +
dt*(
rk_a*ff(i,j,k)+
rk_b*ff0(i,j,k) ) + &
171 if (
cds .eq. 1)
then 173 vux= ( vhx(i,j,k) *uhy(i,j,k) - &
174 vhx(i-1,j,k)*uhy(i-1,j,k) ) /
dx 176 vvy= ( vhy(i,j+1,k)**2.d0 - &
177 vhy(i,j,k)**2.d0 ) /
dy 179 vx = ( vhx(i,j,k)-vhx(i-1,j,k) ) /
dx 181 elseif (
cds .eq. 2)
then 183 vux= ( -vhx(i+1,j,k)*uhy(i+1,j,k) + &
184 27.d0*vhx(i,j,k) *uhy(i,j,k) - &
185 27.d0*vhx(i-1,j,k)*uhy(i-1,j,k) + &
186 vhx(i-2,j,k)*uhy(i-2,j,k) ) / 24.d0 /
dx 188 vvy= ( -vhy(i,j+2,k)**2.d0 + &
189 27.d0*vhy(i,j+1,k)**2.d0 - &
190 27.d0*vhy(i,j,k)**2.d0 + &
191 vhy(i,j-1,k)**2.d0 ) / 24.d0 /
dy 193 vx = ( -vhx(i+1,j,k) + &
195 27.d0*vhx(i-1,j,k) + &
196 vhx(i-2,j,k) ) / 24.d0 /
dx 201 vwz= ( vhz(i,j,k) *why(i,j,k) - &
202 vhz(i,j,k-1)*why(i,j,k-1) ) /
dz(k)
205 if (
cds .eq. 1)
then 207 vxx= ( v(i+1,j,k) - &
211 vyy= ( v(i,j+1,k) - &
215 elseif (
cds .eq. 2)
then 217 vxx= ( -v(i+2,j,k) + &
221 v(i-2,j,k) ) / 12.d0 /
dx2 223 vyy= ( -v(i,j+2,k) + &
227 v(i,j-2,k) ) / 12.d0 /
dy2 232 vzz= ( (v(i,j,k+1)-v(i,j,k))/
dz_t(k) - &
233 (v(i,j,k)-v(i,j,k-1))/
dz_b(k) ) /
dz(k)
240 gg(i,j,k) = -vux-vvy-vwz -
u_mrf*vx + &
241 fc*(
ug - 0.5d0*(uhy(i,j,k)+uhy(i-1,j,k)) -
u_mrf )
245 gg(i,j,k) = -vux-vvy-vwz -
u_mrf*vx
250 if (
cds .eq. 1)
then 252 dpdy = ( p(i,j+1,k)-p(i,j,k) ) /
dy 254 elseif (
cds .eq. 2)
then 256 dpdy= ( -p(i,j+2,k) + &
259 p(i,j-1,k) ) / 24.d0 /
dy 264 qv(i,j,k)=v(i,j,k) +
dt*(
rk_a*gg(i,j,k)+
rk_b*gg0(i,j,k) ) + &
279 if (
cds .eq. 1)
then 281 wux= ( whx(i,j,k) *uhz(i,j,k) - &
282 whx(i-1,j,k)*uhz(i-1,j,k) ) /
dx 284 wvy= ( why(i,j,k) *vhz(i,j,k) - &
285 why(i,j-1,k)*vhz(i,j-1,k) ) /
dy 287 wx = ( whx(i,j,k)-whx(i-1,j,k) ) /
dx 289 elseif (
cds .eq. 2)
then 291 wux= ( -whx(i+1,j,k)*uhz(i+1,j,k) + &
292 27.d0*whx(i,j,k) *uhz(i,j,k) - &
293 27.d0*whx(i-1,j,k)*uhz(i-1,j,k) + &
294 whx(i-2,j,k)*uhz(i-2,j,k) ) / 24.d0 /
dx 296 wvy= ( -why(i,j+1,k)*vhz(i,j+1,k) + &
297 27.d0*why(i,j,k) *vhz(i,j,k) - &
298 27.d0*why(i,j-1,k)*vhz(i,j-1,k) + &
299 why(i,j-2,k)*vhz(i,j-2,k) ) / 24.d0 /
dy 301 wx = ( -whx(i+1,j,k) + &
303 27.d0*whx(i-1,j,k) + &
304 whx(i-2,j,k) ) / 24.d0 /
dx 309 a1 = ( 0.5d0*(w(i,j,k)+w(i,j,k+1)) )**2.d0
310 a2 = ( 0.5d0*(w(i,j,k)+w(i,j,k-1)) )**2.d0
311 wwz= (a1-a2) /
dz_t(k)
314 if (
cds .eq. 1)
then 316 wxx= ( w(i+1,j,k) - &
320 wyy= ( w(i,j+1,k) - &
324 elseif (
cds .eq. 2)
then 326 wxx= ( -w(i+2,j,k) + &
330 w(i-2,j,k) ) / 12.d0 /
dx2 332 wyy= ( -w(i,j+2,k) + &
336 w(i,j-2,k) ) / 12.d0 /
dy2 341 wzz= ( (w(i,j,k+1)-w(i,j,k))/
dz(k+1) - &
342 (w(i,j,k)-w(i,j,k-1))/
dz(k) ) /
dz_t(k)
348 hh(i,j,k) = -wux-wvy-wwz -
u_mrf*wx + &
353 hh(i,j,k) = -wux-wvy-wwz -
u_mrf*wx
356 print *,
'isscalar error!' 361 dpdz = ( p(i,j,k+1)-p(i,j,k) ) /
dz_t(k)
364 qw(i,j,k)=w(i,j,k) +
dt*(
rk_a*hh(i,j,k)+
rk_b*hh0(i,j,k) ) + &
385 real(mytype),
dimension(:,:,:),
intent(in) :: u,v,w,t,ss0,thx,thy,thz
386 real(mytype),
dimension(:,:,:),
intent(out) :: ss,qt
389 real(mytype) :: tux, tvy, twz, txx, tyy, tzz, tx
396 if (
cds .eq. 1)
then 398 tux= ( thx(i,j,k) *u(i,j,k) - &
399 thx(i-1,j,k)*u(i-1,j,k) ) /
dx 401 tvy= ( thy(i,j,k) *v(i,j,k) - &
402 thy(i,j-1,k)*v(i,j-1,k) ) /
dy 404 tx = ( thx(i,j,k)-thx(i-1,j,k) ) /
dx 406 elseif (
cds .eq. 2)
then 408 tux= ( -thx(i+1,j,k)*u(i+1,j,k) + &
409 27.d0*thx(i,j,k) *u(i,j,k) - &
410 27.d0*thx(i-1,j,k)*u(i-1,j,k) + &
411 thx(i-2,j,k)*u(i-2,j,k) ) / 24.d0 /
dx 413 tvy= ( -thy(i,j+1,k)*v(i,j+1,k) + &
414 27.d0*thy(i,j,k) *v(i,j,k) - &
415 27.d0*thy(i,j-1,k)*v(i,j-1,k) + &
416 thy(i,j-2,k)*v(i,j-2,k) ) / 24.d0 /
dy 418 tx = ( -thx(i+1,j,k) + &
420 27.d0*thx(i-1,j,k) + &
421 thx(i-2,j,k) ) / 24.d0 /
dx 426 twz= ( thz(i,j,k)*w(i,j,k) - thz(i,j,k-1)*w(i,j,k-1) ) /
dz(k)
429 if (
cds .eq. 1)
then 431 txx= ( t(i+1,j,k) - &
435 tyy= ( t(i,j+1,k) - &
439 elseif (
cds .eq. 2)
then 441 txx= ( -t(i+2,j,k) + &
445 t(i-2,j,k) ) / 12.d0 /
dx2 447 tyy= ( -t(i,j+2,k) + &
451 t(i,j-2,k) ) / 12.d0 /
dy2 456 tzz= ( (t(i,j,k+1)-t(i,j,k))/
dz_t(k) - &
457 (t(i,j,k)-t(i,j,k-1))/
dz_b(k) ) /
dz(k)
460 ss(i,j,k) = -tux-tvy-twz -
u_mrf*tx
463 qt(i,j,k)= t(i,j,k) +
dt*(
rk_a*ss(i,j,k)+
rk_b*ss0(i,j,k) ) + &
real(mytype), save, protected rk_b
real(mytype), save, protected fc
real(mytype), save, protected dt
real(mytype), save, protected rk_a
integer, save, protected jend3
real(mytype), save, protected nu
real(mytype), save, protected ug
real(mytype), dimension(:), allocatable, save, protected dz
real(mytype), save, protected dy2
integer, save, protected isscalar
real(mytype), save, protected alpha
integer, save, protected istr3
real(mytype), save, protected vg
real(mytype), save, protected t_ref
real(mytype), save, protected dx
integer, save, protected jstr3
real(mytype), dimension(:), allocatable, save, protected dz_b
real(mytype), save, protected dx2
subroutine get_temperature_cnvdiff(u, v, w, t, thx, thy, thz, ss, ss0, qt)
real(mytype), save, protected u_mrf
integer, save, protected cds
real(mytype), save, protected dpx_drive
real(mytype), dimension(:), allocatable, save, protected dz_t
real(mytype), save, protected dy
subroutine update_dpx_drive
real(mytype), save, protected got
real(mytype), save, protected rk_c
integer, save, protected iend3
real(mytype), dimension(:), allocatable, save, protected l_t
subroutine get_momentum_cnvdiff(u, v, w, p, t, uhx, uhy, uhz, vhx, vhy, vhz, whx, why, ff, ff0, gg, gg0, hh, hh0, qu, qv, qw)
integer, save, protected dp_opt
integer, save, protected kstr3
integer, save, protected kend3