SCALE-RM
scale_atmos_phy_tb_common.F90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
13 !-------------------------------------------------------------------------------
14 module scale_atmos_phy_tb_common
15  !-----------------------------------------------------------------------------
16  !
17  !++ used modules
18  !
19  use scale_precision
20  use scale_stdio
21  use scale_prof
22  use scale_grid_index
23  use scale_tracer
24 
25 #ifdef DEBUG
26  use scale_debug, only: &
27  check
28  use scale_const, only: &
29  undef => const_undef, &
30  iundef => const_undef2
31 #endif
32  !-----------------------------------------------------------------------------
33  implicit none
34  private
35  !-----------------------------------------------------------------------------
36  !
37  !++ Public procedure
38  !
39  public :: atmos_phy_tb_calc_strain_tensor
40  public :: atmos_phy_tb_diffusion_solver
41  public :: atmos_phy_tb_calc_tend_momz
42  public :: atmos_phy_tb_calc_tend_momx
43  public :: atmos_phy_tb_calc_tend_momy
44  public :: atmos_phy_tb_calc_tend_phi
45  public :: atmos_phy_tb_calc_flux_phi
46 
47  !-----------------------------------------------------------------------------
48  !
49  !++ Public parameters & variables
50  !
51  !-----------------------------------------------------------------------------
52  !
53  !++ Private procedure
54  !
55  !-----------------------------------------------------------------------------
56  !
57  !++ Private parameters & variables
58  !
59  !-----------------------------------------------------------------------------
60 contains
61  !-----------------------------------------------------------------------------
62  subroutine atmos_phy_tb_calc_strain_tensor( &
63  S33_C, S11_C, S22_C, &
64  S31_C, S12_C, S23_C, &
65  S12_Z, S23_X, S31_Y, &
66  S2, &
67  DENS, MOMZ, MOMX, MOMY, &
68  GSQRT, J13G, J23G, J33G, MAPF )
69  use scale_gridtrans, only: &
70  i_xyz, &
71  i_xyw, &
72  i_uyw, &
73  i_xvw, &
74  i_uyz, &
75  i_xvz, &
76  i_uvz, &
77  i_xy, &
78  i_uy, &
79  i_xv, &
80  i_uv
81  use scale_grid, only: &
82  fdz => grid_fdz, &
83  fdx => grid_fdx, &
84  fdy => grid_fdy, &
85  rcdz => grid_rcdz, &
86  rcdx => grid_rcdx, &
87  rcdy => grid_rcdy, &
88  rfdz => grid_rfdz, &
89  rfdx => grid_rfdx, &
90  rfdy => grid_rfdy
91  implicit none
92 
93  real(RP), intent(out) :: S33_C (ka,ia,ja) ! (cell center)
94  real(RP), intent(out) :: S11_C (ka,ia,ja)
95  real(RP), intent(out) :: S22_C (ka,ia,ja)
96  real(RP), intent(out) :: S31_C (ka,ia,ja)
97  real(RP), intent(out) :: S12_C (ka,ia,ja)
98  real(RP), intent(out) :: S23_C (ka,ia,ja)
99  real(RP), intent(out) :: S12_Z (ka,ia,ja) ! (z edge or x-y plane)
100  real(RP), intent(out) :: S23_X (ka,ia,ja) ! (x edge or y-z plane)
101  real(RP), intent(out) :: S31_Y (ka,ia,ja) ! (y edge or z-x plane)
102  real(RP), intent(out) :: S2 (ka,ia,ja) ! |S|^2
103 
104  real(RP), intent(in) :: DENS (ka,ia,ja)
105  real(RP), intent(in) :: MOMZ (ka,ia,ja)
106  real(RP), intent(in) :: MOMX (ka,ia,ja)
107  real(RP), intent(in) :: MOMY (ka,ia,ja)
108 
109  real(RP), intent(in) :: GSQRT (ka,ia,ja,7)
110  real(RP), intent(in) :: J13G (ka,ia,ja,7)
111  real(RP), intent(in) :: J23G (ka,ia,ja,7)
112  real(RP), intent(in) :: J33G
113  real(RP), intent(in) :: MAPF (ia,ja,2,4)
114 
115  ! velocity
116  real(RP) :: VELZ_C (ka,ia,ja)
117  real(RP) :: VELZ_XY(ka,ia,ja)
118  real(RP) :: VELX_C (ka,ia,ja)
119  real(RP) :: VELX_YZ(ka,ia,ja)
120  real(RP) :: VELY_C (ka,ia,ja)
121  real(RP) :: VELY_ZX(ka,ia,ja)
122 
123  ! work space
124  real(RP) :: WORK_V(ka,ia,ja) ! work space (vertex)
125  real(RP) :: WORK_Z(ka,ia,ja) ! (z edge or x-y plane)
126  real(RP) :: WORK_X(ka,ia,ja) ! (x edge or y-z plane)
127  real(RP) :: WORK_Y(ka,ia,ja) ! (y edge or z-x plane)
128 
129  integer :: IIS, IIE, JJS, JJE
130  integer :: k, i, j
131 
132 #ifdef DEBUG
133  s33_c(:,:,:) = undef
134  s11_c(:,:,:) = undef
135  s22_c(:,:,:) = undef
136  s31_c(:,:,:) = undef
137  s12_c(:,:,:) = undef
138  s23_c(:,:,:) = undef
139  s12_z(:,:,:) = undef
140  s23_x(:,:,:) = undef
141  s31_y(:,:,:) = undef
142  s2(:,:,:) = undef
143 
144  velz_c(:,:,:) = undef
145  velz_xy(:,:,:) = undef
146  velx_c(:,:,:) = undef
147  velx_yz(:,:,:) = undef
148  vely_c(:,:,:) = undef
149  vely_zx(:,:,:) = undef
150 
151  work_v(:,:,:) = undef
152  work_z(:,:,:) = undef
153  work_x(:,:,:) = undef
154  work_y(:,:,:) = undef
155 #endif
156 
157  ! momentum -> velocity
158  do j = js-1, je+1
159  do i = is-1, ie+1
160  do k = ks, ke-1
161 #ifdef DEBUG
162  call check( __line__, momz(k,i,j) )
163  call check( __line__, dens(k+1,i,j) )
164  call check( __line__, dens(k,i,j) )
165 #endif
166  velz_xy(k,i,j) = 2.0_rp * momz(k,i,j) / ( dens(k+1,i,j)+dens(k,i,j) )
167  enddo
168  enddo
169  enddo
170 #ifdef DEBUG
171  i = iundef; j = iundef; k = iundef
172 #endif
173  do j = js-1, je+1
174  do i = is-1, ie+1
175  velz_xy(ke,i,j) = 0.0_rp
176  enddo
177  enddo
178 #ifdef DEBUG
179  i = iundef; j = iundef; k = iundef
180 #endif
181  do j = js-2, je+2
182  do i = is-2, ie+2
183  do k = ks+1, ke
184 #ifdef DEBUG
185  call check( __line__, momz(k,i,j) )
186  call check( __line__, momz(k-1,i,j) )
187  call check( __line__, dens(k,i,j) )
188 #endif
189  velz_c(k,i,j) = 0.5_rp * ( momz(k,i,j) + momz(k-1,i,j) ) / dens(k,i,j)
190  enddo
191  enddo
192  enddo
193 #ifdef DEBUG
194  i = iundef; j = iundef; k = iundef
195 #endif
196  do j = js-2, je+2
197  do i = is-2, ie+2
198 #ifdef DEBUG
199  call check( __line__, momz(ks,i,j) )
200  call check( __line__, dens(ks,i,j) )
201 #endif
202  velz_c(ks,i,j) = 0.5_rp * momz(ks,i,j) / dens(ks,i,j) ! MOMZ(KS-1,i,j) = 0
203  enddo
204 
205  enddo
206 #ifdef DEBUG
207  i = iundef; j = iundef; k = iundef
208 #endif
209 
210  do j = js-1, je+1
211  do i = is-2, ie+1
212  do k = ks, ke
213 #ifdef DEBUG
214  call check( __line__, momx(k,i,j) )
215  call check( __line__, dens(k,i+1,j) )
216  call check( __line__, dens(k,i,j) )
217 #endif
218  velx_yz(k,i,j) = 2.0_rp * momx(k,i,j) / ( dens(k,i+1,j)+dens(k,i,j) )
219  enddo
220  enddo
221  enddo
222 #ifdef DEBUG
223  i = iundef; j = iundef; k = iundef
224 #endif
225  do j = js-1, je+1
226  do i = is-2, ie+1
227  velx_yz(ke+1,i,j) = 0.0_rp
228  enddo
229  enddo
230 #ifdef DEBUG
231  i = iundef; j = iundef; k = iundef
232 #endif
233  do j = js-2, je+2
234  do i = is-1, ie+2
235  do k = ks, ke
236 #ifdef DEBUG
237  call check( __line__, momx(k,i,j) )
238  call check( __line__, momx(k,i-1,j) )
239  call check( __line__, dens(k,i,j) )
240 #endif
241  velx_c(k,i,j) = 0.5_rp * ( momx(k,i,j) + momx(k,i-1,j) ) / dens(k,i,j)
242  enddo
243  enddo
244  enddo
245 #ifdef DEBUG
246  i = iundef; j = iundef; k = iundef
247 #endif
248 
249  do j = js-2, je+1
250  do i = is-1, ie+1
251  do k = ks, ke
252 #ifdef DEBUG
253  call check( __line__, momy(k,i,j) )
254  call check( __line__, dens(k,i,j+1) )
255  call check( __line__, dens(k,i,j) )
256 #endif
257  vely_zx(k,i,j) = 2.0_rp * momy(k,i,j) / ( dens(k,i,j+1)+dens(k,i,j) )
258  enddo
259  enddo
260  enddo
261 #ifdef DEBUG
262  i = iundef; j = iundef; k = iundef
263 #endif
264  do j = js-2, je+1
265  do i = is-1, ie+1
266  vely_zx(ke+1,i,j) = 0.0_rp
267  enddo
268  enddo
269 #ifdef DEBUG
270  i = iundef; j = iundef; k = iundef
271 #endif
272  do j = js-1, je+2
273  do i = is-2, ie+2
274  do k = ks, ke
275 #ifdef DEBUG
276  call check( __line__, momy(k,i,j) )
277  call check( __line__, momy(k,i,j-1) )
278  call check( __line__, dens(k,i,j) )
279 #endif
280  vely_c(k,i,j) = 0.5_rp * ( momy(k,i,j) + momy(k,i,j-1) ) / dens(k,i,j)
281  enddo
282  enddo
283  enddo
284 #ifdef DEBUG
285  i = iundef; j = iundef; k = iundef
286 #endif
287 
288  do jjs = js, je, jblock
289  jje = jjs+jblock-1
290  do iis = is, ie, iblock
291  iie = iis+iblock-1
292 
293 #ifdef DEBUG
294  work_z(:,:,:) = undef; work_x(:,:,:) = undef; work_y(:,:,:) = undef; work_v(:,:,:) = undef
295 #endif
296  ! w
297  ! (x-y plane; x,y,w)
298  ! WORK_Z = VELZ_XY
299  ! (y-z plane; u,y,z)
300  do j = jjs-1, jje+1
301  do i = iis-1, iie+1
302  do k = ks, ke
303 #ifdef DEBUG
304  call check( __line__, velz_c(k,i+1,j) )
305  call check( __line__, velz_c(k,i,j) )
306 #endif
307  work_x(k,i,j) = 0.5_rp * ( velz_c(k,i+1,j) + velz_c(k,i,j) )
308  enddo
309  enddo
310  enddo
311 #ifdef DEBUG
312  i = iundef; j = iundef; k = iundef
313 #endif
314  do j = jjs-1, jje+1
315  do i = iis-1, iie+1
316  work_x(ke+1,i,j) = 0.0_rp
317  enddo
318  enddo
319 #ifdef DEBUG
320  i = iundef; j = iundef; k = iundef
321 #endif
322  ! (z-x plane; x,v,z)
323  do j = jjs-1, jje+1
324  do i = iis-1, iie+1
325  do k = ks, ke
326 #ifdef DEBUG
327  call check( __line__, velz_c(k,i,j+1) )
328  call check( __line__, velz_c(k,i,j) )
329 #endif
330  work_y(k,i,j) = 0.5_rp * ( velz_c(k,i,j+1) + velz_c(k,i,j) )
331  enddo
332  enddo
333  enddo
334 #ifdef DEBUG
335  i = iundef; j = iundef; k = iundef
336 #endif
337  do j = jjs-1, jje+1
338  do i = iis-1, iie+1
339  work_y(ke+1,i,j) = 0.0_rp
340  enddo
341  enddo
342 #ifdef DEBUG
343  i = iundef; j = iundef; k = iundef
344 #endif
345 
346  ! dw/dz
347  ! (cell center; x,y,z)
348  do j = jjs-1, jje+1
349  do i = iis-1, iie+1
350  do k = ks+1, ke
351 #ifdef DEBUG
352  call check( __line__, velz_xy(k,i,j) )
353  call check( __line__, velz_xy(k-1,i,j) )
354  call check( __line__, rcdz(k) )
355 #endif
356  s33_c(k,i,j) = ( velz_xy(k,i,j) - velz_xy(k-1,i,j) ) * rcdz(k) &
357  * j33g / gsqrt(k,i,j,i_xyz)
358  enddo
359  enddo
360  enddo
361 #ifdef DEBUG
362  i = iundef; j = iundef; k = iundef
363 #endif
364  do j = jjs-1, jje+1
365  do i = iis-1, iie+1
366 #ifdef DEBUG
367  call check( __line__, velz_xy(ks,i,j) )
368  call check( __line__, gsqrt(ks,i,j,i_xyz) )
369  call check( __line__, rcdz(ks) )
370 #endif
371  s33_c(ks,i,j) = velz_xy(ks,i,j) * rcdz(ks) & ! VELZ_XY(KS-1,i,j) == 0
372  * j33g / gsqrt(ks,i,j,i_xyz)
373  enddo
374  enddo
375 #ifdef DEBUG
376  i = iundef; j = iundef; k = iundef
377 #endif
378 
379  ! 1/2 * dw/dx
380  ! (cell center; x,y,z)
381  do j = jjs-1, jje+1
382  do i = iis-1, iie+1
383  do k = ks+1, ke-1
384 #ifdef DEBUG
385  call check( __line__, velz_c(k,i+1,j) )
386  call check( __line__, velz_c(k,i-1,j) )
387  call check( __line__, gsqrt(k,i+1,j,i_xyz) )
388  call check( __line__, gsqrt(k,i-1,j,i_xyz) )
389  call check( __line__, velz_xy(k,i,j) )
390  call check( __line__, velz_xy(k-1,i,j) )
391  call check( __line__, j13g(k,i,j,i_xyw) )
392  call check( __line__, j13g(k-1,i,j,i_xyw) )
393  call check( __line__, fdx(i) )
394  call check( __line__, fdx(i-1) )
395 #endif
396  s31_c(k,i,j) = 0.5_rp * ( &
397  ( gsqrt(k,i+1,j,i_xyz)*velz_c(k,i+1,j) - gsqrt(k,i-1,j,i_xyz)*velz_c(k,i-1,j) ) / ( fdx(i) + fdx(i-1) ) &
398  + ( j13g(k,i,j,i_xyw)*velz_xy(k,i,j) - j13g(k-1,i,j,i_xyw)*velz_xy(k-1,i,j) ) * rcdz(k) &
399  ) * mapf(i,j,1,i_xy)
400 
401  enddo
402  enddo
403  enddo
404 #ifdef DEBUG
405  i = iundef; j = iundef; k = iundef
406 #endif
407  do j = jjs-1, jje+1
408  do i = iis-1, iie+1
409 #ifdef DEBUG
410  call check( __line__, velz_c(ks,i+1,j) )
411  call check( __line__, velz_c(ks,i-1,j) )
412  call check( __line__, gsqrt(ks,i+1,j,i_xyz) )
413  call check( __line__, gsqrt(ks,i-1,j,i_xyz) )
414  call check( __line__, velz_xy(ks,i,j) )
415  call check( __line__, j13g(ks,i,j,i_xyw) )
416  call check( __line__, velz_c(ke,i+1,j) )
417  call check( __line__, velz_c(ke,i-1,j) )
418  call check( __line__, gsqrt(ke,i+1,j,i_xyz) )
419  call check( __line__, gsqrt(ke,i-1,j,i_xyz) )
420  call check( __line__, velz_xy(ke,i,j) )
421  call check( __line__, j13g(ke,i,j,i_xyw) )
422  call check( __line__, fdx(i) )
423  call check( __line__, fdx(i-1) )
424 #endif
425  s31_c(ks,i,j) = 0.5_rp * ( &
426  ( gsqrt(ks,i+1,j,i_xyz)*velz_c(ks,i+1,j) - gsqrt(ks,i-1,j,i_xyz)*velz_c(ks,i-1,j) ) / ( fdx(i) + fdx(i-1) ) &
427  + ( j13g(ks,i,j,i_xyw)*velz_xy(ks,i,j) ) * rcdz(ks) &
428  ) * mapf(i,j,1,i_xy)
429  s31_c(ke,i,j) = 0.5_rp * ( &
430  ( gsqrt(ke,i+1,j,i_xyz)*velz_c(ke,i+1,j) - gsqrt(ke,i-1,j,i_xyz)*velz_c(ke,i-1,j) ) / ( fdx(i) + fdx(i-1) ) &
431  - ( j13g(ke-1,i,j,i_xyw)*velz_xy(ke-1,i,j) ) * rcdz(ke) &
432  ) * mapf(i,j,1,i_xy)
433  enddo
434  enddo
435 #ifdef DEBUG
436  i = iundef; j = iundef; k = iundef
437 #endif
438 
439  ! (y edge, u,y,w)
440  do j = jjs , jje
441  do i = iis-1, iie
442  do k = ks, ke-1
443 #ifdef DEBUG
444  call check( __line__, velz_xy(k,i+1,j) )
445  call check( __line__, velz_xy(k,i,j) )
446  call check( __line__, rfdx(i) )
447 #endif
448  s31_y(k,i,j) = 0.5_rp * ( &
449  ( gsqrt(k,i+1,j,i_xyw)*velz_xy(k,i+1,j) - gsqrt(k,i,j,i_xyw)*velz_xy(k,i,j) ) * rfdx(i) &
450  + ( j13g(k+1,i,j,i_uyz)*work_x(k+1,i,j) - j13g(k,i,j,i_uyz)*work_x(k,i,j)) * rfdz(k) &
451  ) * mapf(i,j,1,i_uy)
452  enddo
453  enddo
454  enddo
455 #ifdef DEBUG
456  i = iundef; j = iundef; k = iundef
457 #endif
458 
459  ! 1/2 * dw/dy
460  ! (cell center; x,y,z)
461  do j = jjs-1, jje+1
462  do i = iis-1, iie+1
463  do k = ks+1, ke-1
464 #ifdef DEBUG
465  call check( __line__, velz_c(k,i,j+1) )
466  call check( __line__, velz_c(k,i,j-1) )
467  call check( __line__, gsqrt(k,i,j+1,i_xyz) )
468  call check( __line__, gsqrt(k,i,j-1,i_xyz) )
469  call check( __line__, velz_xy(k,i,j) )
470  call check( __line__, velz_xy(k-1,i,j) )
471  call check( __line__, j23g(k,i,j,i_xyw) )
472  call check( __line__, j23g(k-1,i,j,i_xyw) )
473  call check( __line__, fdy(j) )
474  call check( __line__, fdy(j-1) )
475 #endif
476  s23_c(k,i,j) = 0.5_rp * ( &
477  ( gsqrt(k,i,j+1,i_xyz)*velz_c(k,i,j+1) - gsqrt(k,i,j-1,i_xyz)*velz_c(k,i,j-1) ) / ( fdy(j) + fdy(j-1) ) &
478  + ( j23g(k,i,j,i_xyw)*velz_xy(k,i,j) - j23g(k-1,i,j,i_xyw)*velz_xy(k-1,i,j) ) * rcdz(k) &
479  ) * mapf(i,j,2,i_xy)
480  enddo
481  enddo
482  enddo
483 #ifdef DEBUG
484  i = iundef; j = iundef; k = iundef
485 #endif
486  do j = jjs-1, jje+1
487  do i = iis-1, iie+1
488 #ifdef DEBUG
489  call check( __line__, velz_c(ks,i,j+1) )
490  call check( __line__, velz_c(ks,i,j-1) )
491  call check( __line__, gsqrt(ks,i,j+1,i_xyz) )
492  call check( __line__, gsqrt(ks,i,j-1,i_xyz) )
493  call check( __line__, velz_xy(ks,i,j) )
494  call check( __line__, j23g(ks,i,j,i_xyw) )
495  call check( __line__, velz_c(ke,i,j+1) )
496  call check( __line__, velz_c(ke,i,j-1) )
497  call check( __line__, gsqrt(ke,i,j+1,i_xyz) )
498  call check( __line__, gsqrt(ke,i,j-1,i_xyz) )
499  call check( __line__, velz_xy(ke,i,j) )
500  call check( __line__, j23g(ke,i,j,i_xyw) )
501  call check( __line__, fdy(j) )
502  call check( __line__, fdy(j-1) )
503 #endif
504  s23_c(ks,i,j) = 0.5_rp * ( &
505  ( gsqrt(ks,i,j+1,i_xyz)*velz_c(ks,i,j+1) - gsqrt(ks,i,j-1,i_xyz)*velz_c(ks,i,j-1) ) / ( fdy(j) + fdy(j-1) ) &
506  + ( j23g(ks,i,j,i_xyw)*velz_xy(ks,i,j) ) * rcdz(ks) &
507  ) * mapf(i,j,2,i_xy)
508  s23_c(ke,i,j) = 0.5_rp * ( &
509  ( gsqrt(ke,i,j+1,i_xyz)*velz_c(ke,i,j+1) - gsqrt(ke,i,j-1,i_xyz)*velz_c(ke,i,j-1) ) / ( fdy(j) + fdy(j-1) ) &
510  - ( j23g(ke-1,i,j,i_xyw)*velz_xy(ke-1,i,j) ) * rcdz(ke) &
511  ) * mapf(i,j,2,i_xy)
512  enddo
513  enddo
514 #ifdef DEBUG
515  i = iundef; j = iundef; k = iundef
516 #endif
517 
518  ! (x edge; x,v,w)
519  do j = jjs-1, jje
520  do i = iis , iie
521  do k = ks, ke-1
522 #ifdef DEBUG
523  call check( __line__, velz_xy(k,i,j+1) )
524  call check( __line__, velz_xy(k,i,j) )
525  call check( __line__, rfdy(j) )
526 #endif
527  s23_x(k,i,j) = 0.5_rp * ( &
528  ( gsqrt(k,i,j+1,i_xyw)*velz_xy(k,i,j+1) - gsqrt(k,i,j,i_xyw)*velz_xy(k,i,j) ) * rfdy(j) &
529  + ( j23g(k+1,i,j,i_xvz)*work_y(k+1,i,j) - j23g(k,i,j,i_xvz)*work_y(k,i,j) ) * rfdz(k) &
530  ) * mapf(i,j,2,i_xv)
531  enddo
532  enddo
533  enddo
534 #ifdef DEBUG
535  i = iundef; j = iundef; k = iundef
536 #endif
537 
538 #ifdef DEBUG
539  work_z(:,:,:) = undef; work_x(:,:,:) = undef; work_y(:,:,:) = undef; work_v(:,:,:) = undef
540 #endif
541  ! u
542  ! (x-y plane; x,y,w)
543  do j = jjs-1, jje+1
544  do i = iis-1, iie+1
545  do k = ks, ke-1
546 #ifdef DEBUG
547  call check( __line__, velx_c(k+1,i,j) )
548  call check( __line__, velx_c(k,i,j) )
549 #endif
550  work_z(k,i,j) = 0.5_rp * ( velx_c(k+1,i,j) + velx_c(k,i,j) )
551  enddo
552  enddo
553  enddo
554 #ifdef DEBUG
555  i = iundef; j = iundef; k = iundef
556 #endif
557  ! (y-z plane; u,y,z)
558  ! WORK_X = VELX_YZ
559  ! (z-x plane; x,v,z)
560  do j = jjs-1, jje
561  do i = iis-1, iie+1
562  do k = ks, ke
563 #ifdef DEBUG
564  call check( __line__, velx_c(k,i,j+1) )
565  call check( __line__, velx_c(k,i,j) )
566 #endif
567  work_y(k,i,j) = 0.5_rp * ( velx_c(k,i,j+1) + velx_c(k,i,j) )
568  enddo
569  enddo
570  enddo
571 #ifdef DEBUG
572  i = iundef; j = iundef; k = iundef
573 #endif
574  ! (vertex; u,v,w)
575  do j = jjs-1, jje
576  do i = iis-1, iie
577  do k = ks, ke-1
578 #ifdef DEBUG
579  call check( __line__, velx_yz(k,i,j) )
580  call check( __line__, velx_yz(k,i,j+1) )
581  call check( __line__, velx_yz(k+1,i,j) )
582  call check( __line__, velx_yz(k+1,i,j+1) )
583  call check( __line__, j23g(k ,i,j ,i_uvz) )
584  call check( __line__, j23g(k+1,i,j ,i_uvz) )
585  call check( __line__, j23g(k ,i,j+1,i_uvz) )
586  call check( __line__, j23g(k+1,i,j+1,i_uvz) )
587 #endif
588  work_v(k,i,j) = 0.25_rp &
589  * ( j23g(k ,i,j ,i_uyz)*velx_yz(k ,i,j ) &
590  + j23g(k+1,i,j ,i_uyz)*velx_yz(k+1,i,j ) &
591  + j23g(k ,i,j+1,i_uyz)*velx_yz(k ,i,j+1) &
592  + j23g(k+1,i,j+1,i_uyz)*velx_yz(k+1,i,j+1) )
593  enddo
594  enddo
595  enddo
596 #ifdef DEBUG
597  i = iundef; j = iundef; k = iundef
598 #endif
599 
600  ! du/dx
601  ! (cell center; x,y,z)
602  do j = jjs-1, jje+1
603  do i = iis-1, iie+1
604  do k = ks+1, ke-1
605 #ifdef DEBUG
606  call check( __line__, velx_yz(k,i,j) )
607  call check( __line__, velx_yz(k,i-1,j) )
608  call check( __line__, gsqrt(k,i,j,i_uyz) )
609  call check( __line__, gsqrt(k,i-1,j,i_uyz) )
610  call check( __line__, work_z(k,i,j) )
611  call check( __line__, work_z(k-1,i,j) )
612  call check( __line__, j13g(k,i,j,i_xyw) )
613  call check( __line__, j13g(k-1,i,j,i_xyw) )
614  call check( __line__, gsqrt(k,i,j,i_xyz) )
615  call check( __line__, rcdx(i) )
616 #endif
617  s11_c(k,i,j) = ( &
618  ( gsqrt(k,i,j,i_uyz)*velx_yz(k,i,j) - gsqrt(k,i-1,j,i_uyz)*velx_yz(k,i-1,j) ) * rcdx(i) &
619  + ( j13g(k,i,j,i_xyw)*work_z(k,i,j) - j13g(k-1,i,j,i_xyw)*work_z(k-1,i,j) ) * rcdz(k) &
620  ) * mapf(i,j,1,i_xy) / gsqrt(k,i,j,i_xyz)
621  enddo
622  enddo
623  enddo
624 #ifdef DEBUG
625  i = iundef; j = iundef; k = iundef
626 #endif
627  do j = jjs-1, jje+1
628  do i = iis-1, iie+1
629 #ifdef DEBUG
630  call check( __line__, velx_yz(ks,i,j) )
631  call check( __line__, velx_yz(ks,i-1,j) )
632  call check( __line__, gsqrt(ks,i,j,i_uyz) )
633  call check( __line__, gsqrt(ks,i-1,j,i_uyz) )
634  call check( __line__, velx_c(ks+1,i,j) )
635  call check( __line__, velx_c(ks,i,j) )
636  call check( __line__, j13g(ks+1,i,j,i_xyz) )
637  call check( __line__, j13g(ks,i,j,i_xyz) )
638  call check( __line__, gsqrt(ks,i,j,i_xyz) )
639  call check( __line__, velx_yz(ke,i,j) )
640  call check( __line__, velx_yz(ke,i-1,j) )
641  call check( __line__, gsqrt(ke,i,j,i_uyz) )
642  call check( __line__, gsqrt(ke,i-1,j,i_uyz) )
643  call check( __line__, velx_c(ke,i,j) )
644  call check( __line__, velx_c(ke-1,i,j) )
645  call check( __line__, j13g(ke,i,j,i_xyz) )
646  call check( __line__, j13g(ke-1,i,j,i_xyz) )
647  call check( __line__, gsqrt(ke,i,j,i_xyz) )
648  call check( __line__, rcdx(i) )
649 #endif
650  s11_c(ks,i,j) = ( &
651  ( gsqrt(ks,i,j,i_uyz)*velx_yz(ks,i,j) - gsqrt(ks,i-1,j,i_uyz)*velx_yz(ks,i-1,j) ) * rcdx(i) &
652  + ( j13g(ks+1,i,j,i_xyz)*velx_c(ks+1,i,j) - j13g(ks,i,j,i_xyz)*velx_c(ks,i,j) ) * rfdz(ks) &
653  ) * mapf(i,j,1,i_xy) / gsqrt(ks,i,j,i_xyz)
654  s11_c(ke,i,j) = ( &
655  ( gsqrt(ke,i,j,i_uyz)*velx_yz(ke,i,j) - gsqrt(ke,i-1,j,i_uyz)*velx_yz(ke,i-1,j) ) * rcdx(i) &
656  + ( j13g(ke,i,j,i_xyz)*velx_c(ke,i,j) - j13g(ke-1,i,j,i_xyz)*velx_c(ke-1,i,j) ) * rfdz(ke-1) &
657  ) * mapf(i,j,1,i_xy) / gsqrt(ke,i,j,i_xyz)
658  enddo
659  enddo
660 #ifdef DEBUG
661  i = iundef; j = iundef; k = iundef
662 #endif
663 
664  ! 1/2 * du/dz
665  ! (cell center; x,y,z)
666  do j = jjs-1, jje+1
667  do i = iis-1, iie+1
668  do k = ks+1, ke-1
669 #ifdef DEBUG
670  call check( __line__, s31_c(k,i,j) )
671  call check( __line__, velx_c(k+1,i,j) )
672  call check( __line__, velx_c(k-1,i,j) )
673  call check( __line__, fdz(k) )
674  call check( __line__, fdz(k-1) )
675 #endif
676  s31_c(k,i,j) = ( s31_c(k,i,j) & ! dw/dx
677  + 0.5_rp * ( velx_c(k+1,i,j) - velx_c(k-1,i,j) ) * j33g / ( fdz(k) + fdz(k-1) ) &
678  ) / gsqrt(k,i,j,i_xyz)
679  enddo
680  enddo
681  enddo
682 #ifdef DEBUG
683  i = iundef; j = iundef; k = iundef
684 #endif
685  do j = jjs-1, jje+1
686  do i = iis-1, iie+1
687 #ifdef DEBUG
688  call check( __line__, s31_c(ks,i,j) )
689  call check( __line__, velx_c(ks+1,i,j) )
690  call check( __line__, velx_c(ks,i,j) )
691  call check( __line__, rfdz(ks) )
692  call check( __line__, s31_c(ke,i,j) )
693  call check( __line__, velx_c(ke,i,j) )
694  call check( __line__, velx_c(ke-1,i,j) )
695  call check( __line__, rfdz(ke-1) )
696 #endif
697  s31_c(ks,i,j) = ( s31_c(ks,i,j) &
698  + 0.5_rp * ( velx_c(ks+1,i,j) - velx_c(ks,i,j) ) * j33g * rfdz(ks) &
699  ) / gsqrt(ks,i,j,i_xyz)
700  s31_c(ke,i,j) = ( s31_c(ke,i,j) &
701  + 0.5_rp * ( velx_c(ke,i,j) - velx_c(ke-1,i,j) ) * j33g * rfdz(ke-1) &
702  ) / gsqrt(ke,i,j,i_xyz)
703  enddo
704  enddo
705 #ifdef DEBUG
706  i = iundef; j = iundef; k = iundef
707 #endif
708  ! (y edge; u,y,w)
709  do j = jjs , jje
710  do i = iis-1, iie
711  do k = ks, ke-1
712 #ifdef DEBUG
713  call check( __line__, s31_y(k,i,j) )
714  call check( __line__, velx_yz(k+1,i,j) )
715  call check( __line__, velx_yz(k,i,j) )
716  call check( __line__, rfdz(k) )
717 #endif
718  s31_y(k,i,j) = ( s31_y(k,i,j) & ! dw/dx
719  + 0.5_rp * ( velx_yz(k+1,i,j) - velx_yz(k,i,j) ) * j33g * rfdz(k) &
720  ) / gsqrt(k,i,j,i_uyw)
721  enddo
722  enddo
723  enddo
724 #ifdef DEBUG
725  i = iundef; j = iundef; k = iundef
726 #endif
727 
728  ! 1/2 * du/dy
729  ! (cell center; x,y,z)
730  do j = jjs-1, jje+1
731  do i = iis-1, iie+1
732  do k = ks+1, ke-1
733 #ifdef DEBUG
734  call check( __line__, velx_c(k,i,j+1) )
735  call check( __line__, velx_c(k,i,j-1) )
736  call check( __line__, gsqrt(k,i,j+1,i_xyz) )
737  call check( __line__, gsqrt(k,i,j-1,i_xyz) )
738  call check( __line__, work_z(k,i,j) )
739  call check( __line__, work_z(k-1,i,j) )
740  call check( __line__, j23g(k,i,j,i_xyw) )
741  call check( __line__, j23g(k-1,i,j,i_xyw) )
742 #endif
743  s12_c(k,i,j) = 0.5_rp * ( &
744  ( gsqrt(k,i,j+1,i_xyz)*velx_c(k,i,j+1) - gsqrt(k,i,j-1,i_xyz)*velx_c(k,i,j-1) ) / ( fdy(j) + fdy(j-1) ) &
745  + ( j23g(k,i,j,i_xyw)*work_z(k,i,j) - j23g(k-1,i,j,i_xyw)*work_z(k-1,i,j) ) * rcdz(k) &
746  ) * mapf(i,j,2,i_xy) / gsqrt(k,i,j,i_xyz)
747  enddo
748  enddo
749  enddo
750 #ifdef DEBUG
751  i = iundef; j = iundef; k = iundef
752 #endif
753  do j = jjs-1, jje+1
754  do i = iis-1, iie+1
755 #ifdef DEBUG
756  call check( __line__, velx_c(ks,i,j+1) )
757  call check( __line__, velx_c(ks,i,j-1) )
758  call check( __line__, gsqrt(ks,i,j+1,i_xyz) )
759  call check( __line__, gsqrt(ks,i,j-1,i_xyz) )
760  call check( __line__, velx_c(ks+1,i,j) )
761  call check( __line__, velx_c(ks,i,j) )
762  call check( __line__, j23g(ks+1,i,j,i_xyz) )
763  call check( __line__, j23g(ks,i,j,i_xyz) )
764  call check( __line__, gsqrt(ks,i,j,i_xyz) )
765  call check( __line__, velx_c(ke,i,j+1) )
766  call check( __line__, velx_c(ke,i,j-1) )
767  call check( __line__, gsqrt(ke,i,j+1,i_xyz) )
768  call check( __line__, gsqrt(ke,i,j-1,i_xyz) )
769  call check( __line__, velx_c(ke,i,j) )
770  call check( __line__, velx_c(ke-1,i,j) )
771  call check( __line__, j23g(ke,i,j,i_xyz) )
772  call check( __line__, j23g(ke-1,i,j,i_xyz) )
773  call check( __line__, gsqrt(ke,i,j,i_xyz) )
774  call check( __line__, fdy(j) )
775  call check( __line__, fdy(j-1) )
776 #endif
777  s12_c(ks,i,j) = 0.5_rp * ( &
778  ( gsqrt(ks,i,j+1,i_xyz)*velx_c(ks,i,j+1) - gsqrt(ks,i,j-1,i_xyz)*velx_c(ks,i,j-1) ) / ( fdy(j) + fdy(j-1) ) &
779  + ( j23g(ks+1,i,j,i_xyz)*velx_c(ks+1,i,j) - j23g(ks,i,j,i_xyz)*velx_c(ks,i,j) ) * rfdz(ks) &
780  ) * mapf(i,j,2,i_xy) / gsqrt(ks,i,j,i_xyz)
781  s12_c(ke,i,j) = 0.5_rp * ( &
782  ( gsqrt(ke,i,j+1,i_xyz)*velx_c(ke,i,j+1) - gsqrt(ke,i,j-1,i_xyz)*velx_c(ke,i,j-1) ) / ( fdy(j) + fdy(j-1) ) &
783  + ( j23g(ke,i,j,i_xyz)*velx_c(ke,i,j) - j23g(ke-1,i,j,i_xyz)*velx_c(ke-1,i,j) ) * rfdz(ke-1) &
784  ) * mapf(i,j,2,i_xy) / gsqrt(ke,i,j,i_xyz)
785  enddo
786  enddo
787 #ifdef DEBUG
788  i = iundef; j = iundef; k = iundef
789 #endif
790 
791  ! (z edge; u,v,z)
792  do j = jjs-1, jje
793  do i = iis-1, iie
794  do k = ks+1, ke-1
795 #ifdef DEBUG
796  call check( __line__, velx_yz(k,i,j+1) )
797  call check( __line__, velx_yz(k,i,j) )
798  call check( __line__, work_v(k,i,j) )
799  call check( __line__, work_v(k-1,i,j) )
800  call check( __line__, rfdy(j) )
801 #endif
802  s12_z(k,i,j) = 0.5_rp * ( &
803  ( gsqrt(k,i,j+1,i_uyz)*velx_yz(k,i,j+1) - gsqrt(k,i,j,i_uyz)*velx_yz(k,i,j) ) * rfdy(j) &
804  + ( work_v(k,i,j) - work_v(k-1,i,j) ) * rcdz(k) &
805  ) * mapf(i,j,2,i_uv)
806  enddo
807  enddo
808  enddo
809 #ifdef DEBUG
810  i = iundef; j = iundef; k = iundef
811 #endif
812  do j = jjs-1, jje
813  do i = iis-1, iie
814 #ifdef DEBUG
815  call check( __line__, velx_yz(ks,i,j+1) )
816  call check( __line__, velx_yz(ks,i,j) )
817  call check( __line__, velx_yz(ks+1,i,j) )
818  call check( __line__, velx_yz(ks+1,i,j+1) )
819  call check( __line__, j23g(ks+1,i,j,i_uvz) )
820  call check( __line__, j23g(ks ,i,j,i_uvz) )
821  call check( __line__, velx_yz(ke,i,j+1) )
822  call check( __line__, velx_yz(ke,i,j) )
823  call check( __line__, velx_yz(ke-1,i,j) )
824  call check( __line__, velx_yz(ke-1,i,j+1) )
825  call check( __line__, j23g(ke ,i,j,i_uvz) )
826  call check( __line__, j23g(ke-1,i,j,i_uvz) )
827 #endif
828  s12_z(ks,i,j) = 0.25_rp * ( &
829  ( gsqrt(ks,i,j+1,i_uyz)*velx_yz(ks,i,j+1) - gsqrt(ks,i,j,i_uyz)*velx_yz(ks,i,j) ) * rfdy(j) &
830  + ( j23g(ks+1,i,j,i_uvz) * ( velx_yz(ks+1,i,j) + velx_yz(ks+1,i,j+1) ) &
831  - j23g(ks ,i,j,i_uvz) * ( velx_yz(ks ,i,j) + velx_yz(ks ,i,j+1) ) ) * rfdz(ks) &
832  ) * mapf(i,j,2,i_uv)
833  s12_z(ke,i,j) = 0.25_rp * ( &
834  ( gsqrt(ke,i,j+1,i_uyz)*velx_yz(ke,i,j+1) - gsqrt(ke,i,j,i_uyz)*velx_yz(ke,i,j) ) * rfdy(j) &
835  + ( j23g(ke ,i,j,i_uvz) * ( velx_yz(ke ,i,j) + velx_yz(ke ,i,j+1) ) &
836  - j23g(ke-1,i,j,i_uvz) * ( velx_yz(ke-1,i,j) + velx_yz(ke-1,i,j+1) ) ) * rfdz(ke-1) &
837  ) * mapf(i,j,2,i_uv)
838  enddo
839  enddo
840 #ifdef DEBUG
841  i = iundef; j = iundef; k = iundef
842 #endif
843 
844 #ifdef DEBUG
845  work_z(:,:,:) = undef; work_x(:,:,:) = undef; work_y(:,:,:) = undef; work_v(:,:,:) = undef
846 #endif
847  ! v
848  ! (x-y plane; x,y,w)
849  do j = jjs-1, jje+1
850  do i = iis-1, iie+1
851  do k = ks, ke-1
852 #ifdef DEBUG
853  call check( __line__, vely_c(k+1,i,j) )
854  call check( __line__, vely_c(k,i,j) )
855 #endif
856  work_z(k,i,j) = 0.5_rp * ( vely_c(k+1,i,j) + vely_c(k,i,j) )
857  enddo
858  enddo
859  enddo
860 #ifdef DEBUG
861  i = iundef; j = iundef; k = iundef
862 #endif
863  ! (y-z plane; u,y,z)
864  do j = jjs-1, jje+1
865  do i = iis-1, iie
866  do k = ks, ke
867 #ifdef DEBUG
868  call check( __line__, vely_c(k,i+1,j) )
869  call check( __line__, vely_c(k,i,j) )
870 #endif
871  work_x(k,i,j) = 0.5_rp * ( velz_c(k,i+1,j) + velz_c(k,i,j) )
872  enddo
873  enddo
874  enddo
875 #ifdef DEBUG
876  i = iundef; j = iundef; k = iundef
877 #endif
878  ! (z-x plane; x,v,z)
879  ! WORK_Y = VELY_ZX
880  ! (vertex; u,v,w)
881  do j = jjs-1, jje
882  do i = iis-1, iie
883  do k = ks, ke-1
884 #ifdef DEBUG
885  call check( __line__, vely_zx(k,i,j) )
886  call check( __line__, vely_zx(k+1,i,j) )
887  call check( __line__, vely_zx(k,i+1,j) )
888  call check( __line__, vely_zx(k+1,i+1,j) )
889 #endif
890  work_v(k,i,j) = 0.25_rp &
891  * ( j13g(k ,i ,j,i_xvz)*vely_zx(k ,i ,j) &
892  + j13g(k+1,i ,j,i_xvz)*vely_zx(k+1,i ,j) &
893  + j13g(k ,i+1,j,i_xvz)*vely_zx(k ,i+1,j) &
894  + j13g(k+1,i+1,j,i_xvz)*vely_zx(k+1,i+1,j) )
895  enddo
896  enddo
897  enddo
898 #ifdef DEBUG
899  i = iundef; j = iundef; k = iundef
900 #endif
901 
902  ! dv/dy
903  ! (cell center; x,y,z)
904  do j = jjs-1, jje+1
905  do i = iis-1, iie+1
906  do k = ks+1, ke-1
907 #ifdef DEBUG
908  call check( __line__, vely_zx(k,i,j) )
909  call check( __line__, vely_zx(k,i,j-1) )
910  call check( __line__, gsqrt(k,i,j,i_xvz) )
911  call check( __line__, gsqrt(k,i,j-1,i_xvz) )
912  call check( __line__, work_z(k,i,j) )
913  call check( __line__, work_z(k-1,i,j) )
914  call check( __line__, j23g(k,i,j,i_xyw) )
915  call check( __line__, j23g(k-1,i,j,i_xyw) )
916  call check( __line__, rcdy(j) )
917 #endif
918  s22_c(k,i,j) = ( &
919  ( gsqrt(k,i,j,i_xvz)*vely_zx(k,i,j) - gsqrt(k,i,j-1,i_xvz)*vely_zx(k,i,j-1) ) * rcdy(j) &
920  + ( j23g(k,i,j,i_xyw)*work_z(k,i,j) - j23g(k-1,i,j,i_xyw)*work_z(k-1,i,j) ) * rcdz(k) &
921  ) * mapf(i,j,2,i_xy) / gsqrt(k,i,j,i_xyz)
922  enddo
923  enddo
924  enddo
925 #ifdef DEBUG
926  i = iundef; j = iundef; k = iundef
927 #endif
928  do j = jjs-1, jje+1
929  do i = iis-1, iie+1
930 #ifdef DEBUG
931  call check( __line__, vely_zx(ks,i,j) )
932  call check( __line__, vely_zx(ks,i,j-1) )
933  call check( __line__, gsqrt(ks,i,j,i_xvz) )
934  call check( __line__, gsqrt(ks,i,j-1,i_xvz) )
935  call check( __line__, vely_c(ks+1,i,j) )
936  call check( __line__, vely_c(ks,i,j) )
937  call check( __line__, j23g(ks+1,i,j,i_xyz) )
938  call check( __line__, j23g(ks,i,j,i_xyz) )
939  call check( __line__, rcdy(j) )
940  call check( __line__, vely_zx(ke,i,j) )
941  call check( __line__, vely_zx(ke,i,j-1) )
942  call check( __line__, gsqrt(ke,i,j,i_xvz) )
943  call check( __line__, gsqrt(ke,i,j-1,i_xvz) )
944  call check( __line__, vely_c(ke,i,j) )
945  call check( __line__, vely_c(ke-1,i,j) )
946  call check( __line__, j23g(ke,i,j,i_xyz) )
947  call check( __line__, j23g(ke-1,i,j,i_xyz) )
948 #endif
949  s22_c(ks,i,j) = ( &
950  ( gsqrt(ks,i,j,i_xvz)*vely_zx(ks,i,j) - gsqrt(ks,i,j-1,i_xvz)*vely_zx(ks,i,j-1) ) * rcdy(j) &
951  + ( j23g(ks+1,i,j,i_xyz)*vely_c(ks+1,i,j) - j23g(ks,i,j,i_xyz)*vely_c(ks,i,j) ) * rfdz(ks) &
952  ) * mapf(i,j,2,i_xy) / gsqrt(ks,i,j,i_xyz)
953  s22_c(ke,i,j) = ( &
954  ( gsqrt(ke,i,j,i_xvz)*vely_zx(ke,i,j) - gsqrt(ke,i,j-1,i_xvz)*vely_zx(ke,i,j-1) ) * rcdy(j) &
955  + ( j23g(ke,i,j,i_xyz)*vely_c(ke,i,j) - j23g(ke-1,i,j,i_xyz)*vely_c(ke-1,i,j) ) * rfdz(ke-1) &
956  ) * mapf(i,j,2,i_xy) / gsqrt(ke,i,j,i_xyz)
957  enddo
958  enddo
959 #ifdef DEBUG
960  i = iundef; j = iundef; k = iundef
961 #endif
962 
963  ! 1/2 * dv/dx
964  ! (cell center; x,y,z)
965  do j = jjs-1, jje+1
966  do i = iis-1, iie+1
967  do k = ks+1, ke-1
968 #ifdef DEBUG
969  call check( __line__, s12_c(k,i,j) )
970  call check( __line__, vely_c(k,i+1,j) )
971  call check( __line__, vely_c(k,i-1,j) )
972  call check( __line__, gsqrt(k,i+1,j,i_xyz) )
973  call check( __line__, gsqrt(k,i-1,j,i_xyz) )
974  call check( __line__, work_z(k,i,j) )
975  call check( __line__, work_z(k-1,i,j) )
976  call check( __line__, j13g(k,i,j,i_xyw) )
977  call check( __line__, j13g(k-1,i,j,i_xyw) )
978  call check( __line__, gsqrt(k,i,j,i_xyz) )
979  call check( __line__, fdx(i) )
980  call check( __line__, fdx(i-1) )
981 #endif
982  s12_c(k,i,j) = ( s12_c(k,i,j) & ! du/dy
983  + 0.5_rp * ( &
984  ( gsqrt(k,i+1,j,i_xyz)*vely_c(k,i+1,j) - gsqrt(k,i-1,j,i_xyz)*vely_c(k,i-1,j) ) / ( fdx(i) + fdx(i-1) ) &
985  + ( j13g(k,i,j,i_xyw)*work_z(k,i,j) - j13g(k-1,i,j,i_xyw)*work_z(k-1,i,j) ) * rcdz(k) ) * mapf(i,j,1,i_xy) &
986  ) / gsqrt(k,i,j,i_xyz)
987  enddo
988  enddo
989  enddo
990 #ifdef DEBUG
991  i = iundef; j = iundef; k = iundef
992 #endif
993  do j = jjs-1, jje+1
994  do i = iis-1, iie+1
995 #ifdef DEBUG
996  call check( __line__, s12_c(ks,i,j) )
997  call check( __line__, vely_c(ks,i+1,j) )
998  call check( __line__, vely_c(ks,i-1,j) )
999  call check( __line__, gsqrt(ks,i+1,j,i_xyz) )
1000  call check( __line__, gsqrt(ks,i-1,j,i_xyz) )
1001  call check( __line__, vely_c(ks+1,i,j) )
1002  call check( __line__, vely_c(ks,i,j) )
1003  call check( __line__, j13g(ks+1,i,j,i_xyz) )
1004  call check( __line__, j13g(ks,i,j,i_xyz) )
1005  call check( __line__, gsqrt(ks,i,j,i_xyz) )
1006  call check( __line__, s12_c(ke,i,j) )
1007  call check( __line__, vely_c(ke,i+1,j) )
1008  call check( __line__, vely_c(ke,i-1,j) )
1009  call check( __line__, gsqrt(ke,i+1,j,i_xyz) )
1010  call check( __line__, gsqrt(ke,i-1,j,i_xyz) )
1011  call check( __line__, vely_c(ke,i,j) )
1012  call check( __line__, vely_c(ke-1,i,j) )
1013  call check( __line__, j13g(ke,i,j,i_xyz) )
1014  call check( __line__, j13g(ke-1,i,j,i_xyz) )
1015  call check( __line__, gsqrt(ke,i,j,i_xyz) )
1016  call check( __line__, fdx(i) )
1017  call check( __line__, fdx(i-1) )
1018 #endif
1019  s12_c(ks,i,j) = ( s12_c(ks,i,j) & ! du/dy
1020  + 0.5_rp * ( &
1021  ( gsqrt(ks,i+1,j,i_xyz)*vely_c(ks,i+1,j) - gsqrt(ks,i-1,j,i_xyz)*vely_c(ks,i-1,j) ) / ( fdx(i) + fdx(i-1) ) &
1022  + ( j13g(ks+1,i,j,i_xyz)*vely_c(ks+1,i,j) - j13g(ks,i,j,i_xyz)*vely_c(ks,i,j) ) * rfdz(ks) ) &
1023  * mapf(i,j,1,i_xy) &
1024  ) / gsqrt(ks,i,j,i_xyz)
1025  s12_c(ke,i,j) = ( s12_c(ke,i,j) & ! du/dy
1026  + 0.5_rp * ( &
1027  ( gsqrt(ke,i+1,j,i_xyz)*vely_c(ke,i+1,j) - gsqrt(ke,i-1,j,i_xyz)*vely_c(ke,i-1,j) ) / ( fdx(i) + fdx(i-1) ) &
1028  + ( j13g(ke,i,j,i_xyz)*vely_c(ke,i,j) - j13g(ke-1,i,j,i_xyz)*vely_c(ke-1,i,j) ) * rfdz(ke-1) ) &
1029  * mapf(i,j,1,i_xy) &
1030  ) / gsqrt(ke,i,j,i_xyz)
1031  enddo
1032  enddo
1033 #ifdef DEBUG
1034  i = iundef; j = iundef; k = iundef
1035 #endif
1036  ! (z edge; u,v,z)
1037  do j = jjs-1, jje
1038  do i = iis-1, iie
1039  do k = ks+1, ke-1
1040 #ifdef DEBUG
1041  call check( __line__, s12_z(k,i,j) )
1042  call check( __line__, vely_zx(k,i+1,j) )
1043  call check( __line__, vely_zx(k,i,j) )
1044  call check( __line__, work_v(k,i,j) )
1045  call check( __line__, work_v(k-1,i,j) )
1046  call check( __line__, rfdx(i) )
1047 #endif
1048  s12_z(k,i,j) = ( s12_z(k,i,j) &
1049  + 0.5_rp * ( &
1050  ( gsqrt(k,i+1,j,i_xvz)*vely_zx(k,i+1,j) - gsqrt(k,i,j,i_xvz)*vely_zx(k,i,j) ) * rfdx(i) &
1051  + ( work_v(k,i,j) - work_v(k-1,i,j) ) * rcdz(k) ) * mapf(i,j,1,i_uv) &
1052  ) / gsqrt(k,i,j,i_uvz)
1053  enddo
1054  enddo
1055  enddo
1056 #ifdef DEBUG
1057  i = iundef; j = iundef; k = iundef
1058 #endif
1059  do j = jjs-1, jje
1060  do i = iis-1, iie
1061 #ifdef DEBUG
1062  call check( __line__, s12_z(ks,i,j) )
1063  call check( __line__, vely_zx(ks,i+1,j) )
1064  call check( __line__, vely_zx(ks,i,j) )
1065  call check( __line__, vely_zx(ks+1,i,j) )
1066  call check( __line__, vely_zx(ks+1,i+1,j) )
1067  call check( __line__, s12_z(ke,i,j) )
1068  call check( __line__, vely_zx(ke,i+1,j) )
1069  call check( __line__, vely_zx(ke,i,j) )
1070  call check( __line__, vely_zx(ke-1,i,j) )
1071  call check( __line__, vely_zx(ke-1,i+1,j) )
1072  call check( __line__, rfdx(i) )
1073 #endif
1074  s12_z(ks,i,j) = ( s12_z(ks,i,j) &
1075  + 0.5_rp * ( &
1076  ( gsqrt(ks,i+1,j,i_xvz)*vely_zx(ks,i+1,j) - gsqrt(ks,i,j,i_xvz)*vely_zx(ks,i,j) ) * rfdx(i) &
1077  + ( j13g(ks+1,i,j,i_uvz) * ( vely_zx(ks+1,i,j) + vely_zx(ks+1,i+1,j) ) &
1078  - j13g(ks ,i,j,i_uvz) * ( vely_zx(ks ,i,j) + vely_zx(ks ,i+1,j) ) ) * rfdz(ks) ) * mapf(i,j,1,i_uv) &
1079  ) / gsqrt(ks,i,j,i_uvz)
1080  s12_z(ke,i,j) = ( s12_z(ke,i,j) &
1081  + 0.5_rp * ( &
1082  ( gsqrt(ke,i+1,j,i_xvz)*vely_zx(ke,i+1,j) - gsqrt(ke,i,j,i_xvz)*vely_zx(ke,i,j) ) * rfdx(i) &
1083  + ( j13g(ke ,i,j,i_uvz) * ( vely_zx(ke ,i,j) + vely_zx(ke ,i+1,j) ) &
1084  - j13g(ke-1,i,j,i_uvz) * ( vely_zx(ke-1,i,j) + vely_zx(ke-1,i+1,j) ) ) * rfdz(ke-1) ) * mapf(i,j,1,i_uv) &
1085  ) / gsqrt(ke,i,j,i_uvz)
1086  enddo
1087  enddo
1088 #ifdef DEBUG
1089  i = iundef; j = iundef; k = iundef
1090 #endif
1091 
1092  ! 1/2 * dv/dz
1093  ! (cell center; x,y,z)
1094  do j = jjs-1, jje+1
1095  do i = iis-1, iie+1
1096  do k = ks+1, ke-1
1097 #ifdef DEBUG
1098  call check( __line__, s23_c(k,i,j) )
1099  call check( __line__, vely_c(k+1,i,j) )
1100  call check( __line__, vely_c(k-1,i,j) )
1101  call check( __line__, fdz(k) )
1102  call check( __line__, fdz(k-1) )
1103 #endif
1104  s23_c(k,i,j) = ( s23_c(k,i,j) & ! dw/dy
1105  + 0.5_rp * ( vely_c(k+1,i,j) - vely_c(k-1,i,j) ) * j33g / ( fdz(k) + fdz(k-1) ) &
1106  ) / gsqrt(k,i,j,i_xyz)
1107  enddo
1108  enddo
1109  enddo
1110 #ifdef DEBUG
1111  i = iundef; j = iundef; k = iundef
1112 #endif
1113  do j = jjs-1, jje+1
1114  do i = iis-1, iie+1
1115 #ifdef DEBUG
1116  call check( __line__, s23_c(ks,i,j) )
1117  call check( __line__, vely_c(ks+1,i,j) )
1118  call check( __line__, vely_c(ks,i,j) )
1119  call check( __line__, rfdz(ks) )
1120  call check( __line__, s23_c(ke,i,j) )
1121  call check( __line__, vely_c(ke,i,j) )
1122  call check( __line__, vely_c(ke-1,i,j) )
1123  call check( __line__, rfdz(ke-1) )
1124 #endif
1125  s23_c(ks,i,j) = ( s23_c(ks,i,j) &
1126  + 0.5_rp * ( vely_c(ks+1,i,j) - vely_c(ks,i,j) ) * j33g * rfdz(ks) &
1127  ) / gsqrt(ks,i,j,i_xyz)
1128  s23_c(ke,i,j) = ( s23_c(ke,i,j) &
1129  + 0.5_rp * ( vely_c(ke,i,j) - vely_c(ke-1,i,j) ) * j33g * rfdz(ke-1) &
1130  ) / gsqrt(ke,i,j,i_xyz)
1131  enddo
1132  enddo
1133 #ifdef DEBUG
1134  i = iundef; j = iundef; k = iundef
1135 #endif
1136 
1137  ! (x edge; x,v,w)
1138  do j = jjs-1, jje
1139  do i = iis , iie
1140  do k = ks, ke-1
1141 #ifdef DEBUG
1142  call check( __line__, s23_x(k,i,j) )
1143  call check( __line__, vely_zx(k+1,i,j) )
1144  call check( __line__, vely_zx(k,i,j) )
1145  call check( __line__, rfdz(k) )
1146 #endif
1147  s23_x(k,i,j) = ( s23_x(k,i,j) &
1148  + 0.5_rp * ( vely_zx(k+1,i,j) - vely_zx(k,i,j) ) * j33g * rfdz(k) &
1149  ) / gsqrt(k,i,j,i_xvw)
1150  enddo
1151  enddo
1152  enddo
1153 #ifdef DEBUG
1154  i = iundef; j = iundef; k = iundef
1155 #endif
1156 
1157 
1158  ! |S|^2 = 2*Sij*Sij
1159 #ifdef DEBUG
1160  work_z(:,:,:) = undef; work_x(:,:,:) = undef; work_y(:,:,:) = undef
1161 #endif
1162  ! (cell center)
1163  do j = jjs-1, jje+1
1164  do i = iis-1, iie+1
1165  do k = ks, ke
1166 #ifdef DEBUG
1167  call check( __line__, s11_c(k,i,j) )
1168  call check( __line__, s22_c(k,i,j) )
1169  call check( __line__, s33_c(k,i,j) )
1170  call check( __line__, s31_c(k,i,j) )
1171  call check( __line__, s12_c(k,i,j) )
1172  call check( __line__, s23_c(k,i,j) )
1173 #endif
1174  s2(k,i,j) = max( 1e-10_rp, &
1175  2.0_rp * ( s11_c(k,i,j)**2 + s22_c(k,i,j)**2 + s33_c(k,i,j)**2 ) &
1176  + 4.0_rp * ( s31_c(k,i,j)**2 + s12_c(k,i,j)**2 + s23_c(k,i,j)**2 ) )
1177  enddo
1178  enddo
1179  enddo
1180 #ifdef DEBUG
1181  i = iundef; j = iundef; k = iundef
1182 #endif
1183 
1184  enddo
1185  enddo
1186 
1187  return
1188  end subroutine atmos_phy_tb_calc_strain_tensor
1189 
1190  !-----------------------------------------------------------------------------
1191  subroutine atmos_phy_tb_calc_flux_phi( &
1192  qflx_phi, &
1193  DENS, PHI, Kh, FACT, &
1194  GSQRT, J13G, J23G, J33G, MAPF, &
1195  a, b, c, dt, &
1196  implicit, &
1197  IIS, IIE, JJS, JJE )
1198  use scale_gridtrans, only: &
1199  i_xyz, &
1200  i_xyw, &
1201  i_uyw, &
1202  i_xvw, &
1203  i_uyz, &
1204  i_xvz, &
1205  i_uvz, &
1206  i_xy, &
1207  i_uy, &
1208  i_xv, &
1209  i_uv
1210  use scale_grid, only: &
1211  fdz => grid_fdz, &
1212  fdx => grid_fdx, &
1213  fdy => grid_fdy, &
1214  rcdz => grid_rcdz, &
1215  rcdx => grid_rcdx, &
1216  rcdy => grid_rcdy, &
1217  rfdz => grid_rfdz, &
1218  rfdx => grid_rfdx, &
1219  rfdy => grid_rfdy
1220  implicit none
1221 
1222  real(RP), intent(out) :: qflx_phi(ka,ia,ja,3)
1223  real(RP), intent(in) :: DENS(ka,ia,ja)
1224  real(RP), intent(in) :: PHI(ka,ia,ja)
1225  real(RP), intent(in) :: Kh(ka,ia,ja)
1226  real(RP), intent(in) :: FACT
1227  real(RP), intent(in) :: GSQRT(ka,ia,ja,7)
1228  real(RP), intent(in) :: J13G(ka,ia,ja,7)
1229  real(RP), intent(in) :: J23G(ka,ia,ja,7)
1230  real(RP), intent(in) :: J33G
1231  real(RP), intent(in) :: MAPF(ia,ja,2,4)
1232  real(RP), intent(in) :: a(ka,ia,ja)
1233  real(RP), intent(in) :: b(ka,ia,ja)
1234  real(RP), intent(in) :: c(ka,ia,ja)
1235  real(DP), intent(in) :: dt
1236  logical, intent(in) :: implicit
1237  integer, intent(in) :: IIS
1238  integer, intent(in) :: IIE
1239  integer, intent(in) :: JJS
1240  integer, intent(in) :: JJE
1241 
1242  real(RP) :: TEND(ka,ia,ja)
1243  real(RP) :: d(ka)
1244 
1245  integer :: k, i, j
1246 
1247  ! (x-y plane; x,y,w)
1248  do j = jjs, jje
1249  do i = iis, iie
1250  do k = ks, ke-1
1251 #ifdef DEBUG
1252  call check( __line__, dens(k,i,j) )
1253  call check( __line__, dens(k+1,i,j) )
1254  call check( __line__, kh(k,i,j) )
1255  call check( __line__, kh(k+1,i,j) )
1256  call check( __line__, phi(k+1,i,j) )
1257  call check( __line__, phi(k,i,j) )
1258  call check( __line__, rfdz(k) )
1259 #endif
1260  qflx_phi(k,i,j,zdir) = - 0.25_rp & ! 1/2/2
1261  * ( dens(k,i,j)+dens(k+1,i,j) ) &
1262  * ( kh(k,i,j) + kh(k+1,i,j) ) &
1263  * ( phi(k+1,i,j)-phi(k,i,j) ) * rfdz(k) * j33g &
1264  / gsqrt(k,i,j,i_xyw)
1265  enddo
1266  enddo
1267  enddo
1268 #ifdef DEBUG
1269  i = iundef; j = iundef; k = iundef
1270 #endif
1271  do j = jjs, jje
1272  do i = iis, iie
1273  qflx_phi(ks-1,i,j,zdir) = 0.0_rp
1274  qflx_phi(ke ,i,j,zdir) = 0.0_rp
1275  enddo
1276  enddo
1277 #ifdef DEBUG
1278  i = iundef; j = iundef; k = iundef
1279 #endif
1280 
1281  ! (y-z plane; u,y,z)
1282  do j = jjs, jje
1283  do i = iis-1, iie
1284  do k = ks+1, ke-1
1285 #ifdef DEBUG
1286  call check( __line__, dens(k,i,j) )
1287  call check( __line__, dens(k,i+1,j) )
1288  call check( __line__, kh(k,i,j) )
1289  call check( __line__, kh(k,i+1,j) )
1290  call check( __line__, phi(k,i+1,j) )
1291  call check( __line__, phi(k,i,j) )
1292  call check( __line__, rfdx(i) )
1293 #endif
1294  qflx_phi(k,i,j,xdir) = - 0.25_rp & ! 1/2/2
1295  * ( dens(k,i,j)+dens(k,i+1,j) ) &
1296  * ( kh(k,i,j) + kh(k,i+1,j) ) &
1297  * ( &
1298  ( gsqrt(k,i+1,j,i_xyz) * phi(k,i+1,j) &
1299  - gsqrt(k,i ,j,i_xyz) * phi(k,i ,j) ) * rfdx(i) &
1300  + ( j13g(k+1,i,j,i_uyz) * ( phi(k+1,i+1,j)+phi(k+1,i,j) ) &
1301  - j13g(k-1,i,j,i_uyz) * ( phi(k-1,i+1,j)+phi(k-1,i,j) ) &
1302  ) * 0.5_rp / ( fdz(k)+fdz(k-1) ) &
1303  ) / gsqrt(k,i,j,i_uyz)
1304  enddo
1305  enddo
1306  enddo
1307 #ifdef DEBUG
1308  i = iundef; j = iundef; k = iundef
1309 #endif
1310  do j = jjs, jje
1311  do i = iis-1, iie
1312 #ifdef DEBUG
1313  call check( __line__, dens(ks,i,j) )
1314  call check( __line__, dens(ks,i+1,j) )
1315  call check( __line__, kh(ks,i,j) )
1316  call check( __line__, kh(ks,i+1,j) )
1317  call check( __line__, phi(ks,i+1,j) )
1318  call check( __line__, phi(ks,i,j) )
1319  call check( __line__, rfdx(i) )
1320 #endif
1321  qflx_phi(ks,i,j,xdir) = - 0.25_rp & ! 1/2/2
1322  * ( dens(ks,i,j)+dens(ks,i+1,j) ) &
1323  * ( kh(ks,i,j) + kh(ks,i+1,j) ) &
1324  * ( &
1325  ( gsqrt(ks,i+1,j,i_xyz) * phi(ks,i+1,j) &
1326  - gsqrt(ks,i ,j,i_xyz) * phi(ks,i ,j) ) * rfdx(i) &
1327  + ( j13g(ks+1,i,j,i_uyz) * ( phi(ks+1,i+1,j)+phi(ks+1,i,j) ) &
1328  - j13g(ks ,i,j,i_uyz) * ( phi(ks ,i+1,j)+phi(ks ,i,j) ) &
1329  ) * 0.5_rp * rfdz(ks) &
1330  ) * mapf(i,j,1,i_uy) / gsqrt(ks,i,j,i_uyz)
1331  qflx_phi(ke,i,j,xdir) = - 0.25_rp & ! 1/2/2
1332  * ( dens(ke,i,j)+dens(ke,i+1,j) ) &
1333  * ( kh(ke,i,j) + kh(ke,i+1,j) ) &
1334  * ( &
1335  ( gsqrt(ke,i+1,j,i_xyz) * phi(ke,i+1,j) &
1336  - gsqrt(ke,i ,j,i_xyz) * phi(ke,i ,j) ) * rfdx(i) &
1337  + ( j13g(ke ,i,j,i_uyz) * ( phi(ke ,i+1,j)+phi(ke ,i,j) ) &
1338  - j13g(ke-1,i,j,i_uyz) * ( phi(ke-1,i+1,j)+phi(ke-1,i,j) ) &
1339  ) * 0.5_rp * rfdz(ke-1) &
1340  ) * mapf(i,j,1,i_uy) / gsqrt(ke,i,j,i_uyz)
1341  enddo
1342  enddo
1343 #ifdef DEBUG
1344  i = iundef; j = iundef; k = iundef
1345 #endif
1346  ! (z-x plane; x,v,z)
1347  do j = jjs-1, jje
1348  do i = iis, iie
1349  do k = ks+1, ke-1
1350 #ifdef DEBUG
1351  call check( __line__, dens(k,i,j) )
1352  call check( __line__, dens(k,i,j+1) )
1353  call check( __line__, kh(k,i,j) )
1354  call check( __line__, kh(k,i,j+1) )
1355  call check( __line__, phi(k,i,j+1) )
1356  call check( __line__, phi(k,i,j) )
1357  call check( __line__, rfdy(j) )
1358 #endif
1359  qflx_phi(k,i,j,ydir) = - 0.25_rp &
1360  * ( dens(k,i,j)+dens(k,i,j+1) ) &
1361  * ( kh(k,i,j) + kh(k,i,j+1) ) &
1362  * ( &
1363  ( gsqrt(k,i,j+1,i_xyz) * phi(k,i,j+1) &
1364  - gsqrt(k,i,j ,i_xyz) * phi(k,i,j ) ) * rfdy(j) &
1365  + ( j23g(k+1,i,j,i_xvz) * ( phi(k+1,i,j+1)+phi(k+1,i,j) ) &
1366  - j23g(k-1,i,j,i_xvz) * ( phi(k-1,i,j+1)+phi(k-1,i,j) ) &
1367  ) * 0.5_rp / ( fdz(k)+fdz(k-1) ) &
1368  ) * mapf(i,j,2,i_xv) / gsqrt(k,i,j,i_xvz)
1369  enddo
1370  enddo
1371  enddo
1372 #ifdef DEBUG
1373  i = iundef; j = iundef; k = iundef
1374 #endif
1375  do j = jjs-1, jje
1376  do i = iis, iie
1377 #ifdef DEBUG
1378  call check( __line__, dens(ks,i,j) )
1379  call check( __line__, dens(ks,i,j+1) )
1380  call check( __line__, kh(ks,i,j) )
1381  call check( __line__, kh(ks,i,j+1) )
1382  call check( __line__, phi(ks,i,j+1) )
1383  call check( __line__, phi(ks,i,j) )
1384  call check( __line__, rfdy(j) )
1385 #endif
1386  qflx_phi(ks,i,j,ydir) = - 0.25_rp &
1387  * ( dens(ks,i,j)+dens(ks,i,j+1) ) &
1388  * ( kh(ks,i,j) + kh(ks,i,j+1) ) &
1389  * ( &
1390  ( gsqrt(ks,i,j+1,i_xyz) * phi(ks,i,j+1) &
1391  - gsqrt(ks,i,j ,i_xyz) * phi(ks,i,j ) ) * rfdy(j) &
1392  + ( j23g(ks+1,i,j,i_xvz) * ( phi(ks+1,i,j+1)+phi(ks+1,i,j) ) &
1393  - j23g(ks ,i,j,i_xvz) * ( phi(ks ,i,j+1)+phi(ks ,i,j) ) &
1394  ) * 0.5_rp * rfdz(ks) &
1395  ) * mapf(i,j,2,i_xv) / gsqrt(ks,i,j,i_xvz)
1396  qflx_phi(ke,i,j,ydir) = - 0.25_rp &
1397  * ( dens(ke,i,j)+dens(ke,i,j+1) ) &
1398  * ( kh(ke,i,j) + kh(ke,i,j+1) ) &
1399  * ( &
1400  ( gsqrt(ke,i,j+1,i_xyz) * phi(ke,i,j+1) &
1401  - gsqrt(ke,i,j ,i_xyz) * phi(ke,i,j ) ) * rfdy(j) &
1402  + ( j23g(ke ,i,j,i_xvz) * ( phi(ke ,i,j+1)+phi(ke ,i,j) ) &
1403  - j23g(ke-1,i,j,i_xvz) * ( phi(ke-1,i,j+1)+phi(ke-1,i,j) ) &
1404  ) * 0.5_rp * rfdz(ke-1) &
1405  ) * mapf(i,j,2,i_xv) / gsqrt(ke,i,j,i_xvz)
1406  enddo
1407  enddo
1408 #ifdef DEBUG
1409  i = iundef; j = iundef; k = iundef
1410 #endif
1411 
1412  if ( implicit ) then
1413  call atmos_phy_tb_calc_tend_phi( tend, & ! (out)
1414  qflx_phi, & ! (in)
1415  gsqrt, j13g, j23g, j33g, mapf, & ! (in)
1416  iis, iie, jjs, jje ) ! (in)
1417 
1418  do j = jjs, jje
1419  do i = iis, iie
1420 
1421  do k = ks, ke
1422  d(k) = tend(k,i,j)
1423  end do
1424 
1425  call atmos_phy_tb_diffusion_solver( &
1426  tend(:,i,j), & ! (out)
1427  a(:,i,j), b(:,i,j), c(:,i,j), d, & ! (in)
1428  ke ) ! (in)
1429 
1430  do k = ks, ke-1
1431  qflx_phi(k,i,j,zdir) = qflx_phi(k,i,j,zdir) &
1432  - 0.25_rp & ! 1/2/2
1433  * ( dens(k,i,j)+dens(k+1,i,j) ) &
1434  * ( kh(k,i,j) + kh(k+1,i,j) ) &
1435  * dt * ( tend(k+1,i,j)-tend(k,i,j) ) * rfdz(k) * j33g &
1436  / gsqrt(k,i,j,i_xyw)
1437  end do
1438 
1439  end do
1440  end do
1441 
1442  end if
1443 
1444  return
1445  end subroutine atmos_phy_tb_calc_flux_phi
1446 
1447  !-----------------------------------------------------------------------------
1448 !OCL SERIAL
1449  subroutine atmos_phy_tb_diffusion_solver( &
1450  phi, &
1451  a, b, c, d, &
1452  KE_TB )
1453  implicit none
1454  real(RP), intent(out) :: phi(ka)
1455  real(RP), intent(in) :: a(ka)
1456  real(RP), intent(in) :: b(ka)
1457  real(RP), intent(in) :: c(ka)
1458  real(RP), intent(in) :: d(ka)
1459  integer, intent(in) :: KE_TB
1460  real(RP) :: e(ka)
1461  real(RP) :: f(ka)
1462  real(RP) :: denom
1463  integer :: k
1464 
1465  e(ks) = - a(ks) / b(ks)
1466  f(ks) = d(ks) / b(ks)
1467  do k = ks+1, ke_tb-1
1468  denom = b(k) + c(k)*e(k-1)
1469  e(k) = - a(k) / denom
1470  f(k) = ( d(k) - c(k)*f(k-1) ) / denom
1471  end do
1472 
1473  ! flux at the top boundary is zero
1474  phi(ke_tb) = ( d(ke_tb) - c(ke_tb)*f(ke_tb-1) ) / ( b(ke_tb) + c(ke_tb)*e(ke_tb-1) ) ! = f(KE_PBL)
1475 
1476  do k = ke_tb-1, ks, -1
1477  phi(k) = e(k) * phi(k+1) + f(k)
1478  end do
1479  do k = 1, ks-1
1480  phi(k) = 0.0_rp
1481  end do
1482  do k = ke_tb+1, ka
1483  phi(k) = 0.0_rp
1484  end do
1485 
1486  return
1487  end subroutine atmos_phy_tb_diffusion_solver
1488 
1489  !-----------------------------------------------------------------------------
1490  subroutine atmos_phy_tb_calc_tend_momz( &
1491  MOMZ_t_TB, &
1492  QFLX_MOMZ, &
1493  GSQRT, J13G, J23G, J33G, MAPF, &
1494  IIS, IIE, JJS, JJE )
1495  use scale_grid, only: &
1496  rcdz => grid_rcdz, &
1497  rcdx => grid_rcdx, &
1498  rcdy => grid_rcdy, &
1499  rfdz => grid_rfdz, &
1500  cdz => grid_cdz
1501  use scale_gridtrans, only: &
1502  i_xyz, &
1503  i_xyw, &
1504  i_uyw, &
1505  i_xvw, &
1506  i_xy
1507  implicit none
1508 
1509  real(RP), intent(out) :: MOMZ_t_TB(ka,ia,ja)
1510  real(RP), intent(in) :: QFLX_MOMZ(ka,ia,ja,3)
1511  real(RP), intent(in) :: GSQRT(ka,ia,ja,7)
1512  real(RP), intent(in) :: J13G(ka,ia,ja,7)
1513  real(RP), intent(in) :: J23G(ka,ia,ja,7)
1514  real(RP), intent(in) :: J33G
1515  real(RP), intent(in) :: MAPF(ia,ja,2,4)
1516  integer , intent(in) :: IIS
1517  integer , intent(in) :: IIE
1518  integer , intent(in) :: JJS
1519  integer , intent(in) :: JJE
1520 
1521  integer :: k, i, j
1522 
1523  do j = jjs, jje
1524  do i = iis, iie
1525  do k = ks+1, ke-2
1526  momz_t_tb(k,i,j) = &
1527  - ( ( gsqrt(k,i ,j,i_uyw) * qflx_momz(k,i ,j,xdir) &
1528  - gsqrt(k,i-1,j,i_uyw) * qflx_momz(k,i-1,j,xdir) ) * rcdx(i) * mapf(i,j,1,i_xy) &
1529  + ( gsqrt(k,i,j ,i_xvw) * qflx_momz(k,i,j ,ydir) &
1530  - gsqrt(k,i,j-1,i_xvw) * qflx_momz(k,i,j-1,ydir) ) * rcdy(j) &
1531  + ( ( j13g(k+1,i,j,i_xyz) * ( qflx_momz(k+1,i,j,xdir) + qflx_momz(k+1,i-1,j,xdir) ) &
1532  - j13g(k-1,i,j,i_xyz) * ( qflx_momz(k-1,i,j,xdir) + qflx_momz(k-1,i-1,j,xdir) ) &
1533  ) * mapf(i,j,1,i_xy) &
1534  + ( j23g(k+1,i,j,i_xyz) * ( qflx_momz(k+1,i,j,ydir) + qflx_momz(k+1,i,j-1,ydir) ) &
1535  - j23g(k-1,i,j,i_xyz) * ( qflx_momz(k-1,i,j,ydir) + qflx_momz(k-1,i,j-1,ydir) ) &
1536  ) * mapf(i,j,2,i_xy) &
1537  ) * 0.5_rp / ( cdz(k+1)+cdz(k) ) &
1538  + j33g * ( qflx_momz(k+1,i,j,zdir) - qflx_momz(k,i,j,zdir) ) * rfdz(k) ) &
1539  / gsqrt(k,i,j,i_xyw)
1540  enddo
1541  enddo
1542  enddo
1543 
1544  do j = jjs, jje
1545  do i = iis, iie
1546  momz_t_tb(ks,i,j) = &
1547  - ( ( gsqrt(ks,i ,j,i_uyw) * qflx_momz(ks,i ,j,xdir) &
1548  - gsqrt(ks,i-1,j,i_uyw) * qflx_momz(ks,i-1,j,xdir) ) * rcdx(i) * mapf(i,j,1,i_xy) &
1549  + ( gsqrt(ks,i,j ,i_xvw) * qflx_momz(ks,i,j ,ydir) &
1550  - gsqrt(ks,i,j-1,i_xvw) * qflx_momz(ks,i,j-1,ydir) ) * rcdy(j) * mapf(i,j,2,i_xy) &
1551  + ( ( j13g(ks+1,i,j,i_xyz) * ( qflx_momz(ks+1,i,j,xdir) + qflx_momz(ks+1,i-1,j,xdir) ) &
1552  - j13g(ks ,i,j,i_xyz) * ( qflx_momz(ks ,i,j,xdir) + qflx_momz(ks ,i-1,j,xdir) ) &
1553  ) * mapf(i,j,1,i_xy) &
1554  + ( j23g(ks+1,i,j,i_xyz) * ( qflx_momz(ks+1,i,j,ydir) + qflx_momz(ks+1,i,j-1,ydir) ) &
1555  - j23g(ks ,i,j,i_xyz) * ( qflx_momz(ks ,i,j,ydir) + qflx_momz(ks ,i,j-1,ydir) ) &
1556  ) * mapf(i,j,2,i_xy) &
1557  ) * 0.5_rp * rcdz(ks+1) &
1558  + j33g * ( qflx_momz(ks+1,i,j,zdir) - qflx_momz(ks,i,j,zdir) ) * rfdz(ks) ) &
1559  / gsqrt(ks,i,j,i_xyw)
1560  momz_t_tb(ke-1,i,j) = &
1561  - ( ( gsqrt(ke-1,i ,j,i_uyw) * qflx_momz(ke-1,i ,j,xdir) &
1562  - gsqrt(ke-1,i-1,j,i_uyw) * qflx_momz(ke-1,i-1,j,xdir) ) * rcdx(i) * mapf(i,j,1,i_xy) &
1563  + ( gsqrt(ke-1,i,j ,i_xvw) * qflx_momz(ke-1,i,j ,ydir) &
1564  - gsqrt(ke-1,i,j-1,i_xvw) * qflx_momz(ke-1,i,j-1,ydir) ) * rcdy(j) * mapf(i,j,2,i_xy) &
1565  + ( ( j13g(ke-1,i,j,i_xyz) * ( qflx_momz(ke-1,i,j,xdir) + qflx_momz(ke-1,i-1,j,xdir) ) &
1566  - j13g(ke-2,i,j,i_xyz) * ( qflx_momz(ke-2,i,j,xdir) + qflx_momz(ke-2,i-1,j,xdir) ) &
1567  ) * mapf(i,j,1,i_xy) &
1568  + ( j23g(ke-1,i,j,i_xyz) * ( qflx_momz(ke-1,i,j,ydir) + qflx_momz(ke-1,i,j-1,ydir) ) &
1569  - j23g(ke-2,i,j,i_xyz) * ( qflx_momz(ke-2,i,j,ydir) + qflx_momz(ke-2,i,j-1,ydir) ) &
1570  ) * mapf(i,j,2,i_xy) &
1571  ) * 0.5_rp * rcdz(ke-1) &
1572  + j33g * ( qflx_momz(ke,i,j,zdir) - qflx_momz(ke-1,i,j,zdir) ) * rfdz(ke-1) ) &
1573  / gsqrt(ke-1,i,j,i_xyw)
1574  enddo
1575  enddo
1576 
1577  return
1578  end subroutine atmos_phy_tb_calc_tend_momz
1579 
1580  subroutine atmos_phy_tb_calc_tend_momx( &
1581  MOMX_t_TB, &
1582  QFLX_MOMX, &
1583  GSQRT, J13G, J23G, J33G, MAPF, &
1584  IIS, IIE, JJS, JJE )
1585  use scale_grid, only: &
1586  rcdz => grid_rcdz, &
1587  rcdy => grid_rcdy, &
1588  rfdz => grid_rfdz, &
1589  rfdx => grid_rfdx, &
1590  fdz => grid_fdz
1591  use scale_gridtrans, only: &
1592  i_xyz, &
1593  i_uyw, &
1594  i_uyz, &
1595  i_uvz, &
1596  i_uy
1597  implicit none
1598  real(RP), intent(out) :: MOMX_t_TB(ka,ia,ja)
1599  real(RP), intent(in) :: QFLX_MOMX(ka,ia,ja,3)
1600  real(RP), intent(in) :: GSQRT(ka,ia,ja,7)
1601  real(RP), intent(in) :: MAPF(ia,ja,2,4)
1602  real(RP), intent(in) :: J13G(ka,ia,ja,7)
1603  real(RP), intent(in) :: J23G(ka,ia,ja,7)
1604  real(RP), intent(in) :: J33G
1605  integer , intent(in) :: IIS
1606  integer , intent(in) :: IIE
1607  integer , intent(in) :: JJS
1608  integer , intent(in) :: JJE
1609 
1610  integer :: k, i, j
1611 
1612  do j = jjs, jje
1613  do i = iis, iie
1614  do k = ks+1, ke-1
1615  momx_t_tb(k,i,j) = &
1616  - ( ( gsqrt(k,i+1,j,i_xyz) * qflx_momx(k,i+1,j,xdir) &
1617  - gsqrt(k,i ,j,i_xyz) * qflx_momx(k,i ,j,xdir) ) * rfdx(i) * mapf(i,j,1,i_uy) &
1618  + ( gsqrt(k,i,j ,i_uvz) * qflx_momx(k,i,j ,ydir) &
1619  - gsqrt(k,i,j-1,i_uvz) * qflx_momx(k,i,j-1,ydir) ) * rcdy(j) * mapf(i,j,2,i_uy) &
1620  + ( ( j13g(k+1,i,j,i_uyw) * ( qflx_momx(k+1,i+1,j,xdir) + qflx_momx(k+1,i,j ,xdir) ) &
1621  - j13g(k-1,i,j,i_uyw) * ( qflx_momx(k-1,i+1,j,xdir) + qflx_momx(k-1,i,j ,xdir) ) &
1622  ) * mapf(i,j,1,i_uy) &
1623  + ( j23g(k+1,i,j,i_uyw) * ( qflx_momx(k+1,i ,j,ydir) + qflx_momx(k+1,i,j-1,ydir) ) &
1624  - j23g(k-1,i,j,i_uyw) * ( qflx_momx(k-1,i ,j,ydir) + qflx_momx(k-1,i,j-1,ydir) ) &
1625  ) * mapf(i,j,2,i_uy) &
1626  ) * 0.5_rp / ( fdz(k)+fdz(k-1) ) &
1627  + j33g * ( qflx_momx(k,i,j,zdir) - qflx_momx(k-1,i,j,zdir) ) * rcdz(k) ) &
1628  / gsqrt(k,i,j,i_uyz)
1629  enddo
1630  enddo
1631  enddo
1632  do j = jjs, jje
1633  do i = iis, iie
1634  momx_t_tb(ks,i,j) = &
1635  - ( ( gsqrt(ks,i+1,j,i_xyz) * qflx_momx(ks,i+1,j,xdir) &
1636  - gsqrt(ks,i ,j,i_xyz) * qflx_momx(ks,i ,j,xdir) ) * rfdx(i) * mapf(i,j,1,i_uy) &
1637  + ( gsqrt(ks,i,j ,i_uvz) * qflx_momx(ks,i,j ,ydir) &
1638  - gsqrt(ks,i,j-1,i_uvz) * qflx_momx(ks,i,j-1,ydir) ) * rcdy(j) &
1639  + ( ( j13g(ks+1,i,j,i_uyw) * ( qflx_momx(ks+1,i+1,j,xdir) + qflx_momx(ks+1,i,j ,xdir) ) &
1640  - j13g(ks ,i,j,i_uyw) * ( qflx_momx(ks ,i+1,j,xdir) + qflx_momx(ks ,i,j ,xdir) ) &
1641  ) * mapf(i,j,1,i_uy) &
1642  + ( j23g(ks+1,i,j,i_uyw) * ( qflx_momx(ks+1,i ,j,ydir) + qflx_momx(ks+1,i,j-1,ydir) ) &
1643  - j23g(ks ,i,j,i_uyw) * ( qflx_momx(ks ,i ,j,ydir) + qflx_momx(ks ,i,j-1,ydir) ) &
1644  ) * mapf(i,j,2,i_uy) &
1645  ) * 0.5_rp * rcdz(ks) &
1646  + j33g * ( qflx_momx(ks,i,j,zdir) ) * rfdz(ks) ) &
1647  / gsqrt(ks,i,j,i_uyz)
1648  momx_t_tb(ke,i,j) = &
1649  - ( ( gsqrt(ke,i+1,j,i_xyz) * qflx_momx(ke,i+1,j,xdir) &
1650  - gsqrt(ke,i ,j,i_xyz) * qflx_momx(ke,i ,j,xdir) ) * rfdx(i) * mapf(i,j,1,i_uy) &
1651  + ( gsqrt(ke,i,j ,i_uvz) * qflx_momx(ke,i,j ,ydir) &
1652  - gsqrt(ke,i,j-1,i_uvz) * qflx_momx(ke,i,j-1,ydir) ) * rcdy(j) * mapf(i,j,2,i_uy)&
1653  + ( ( j13g(ke ,i,j,i_uyw) * ( qflx_momx(ke ,i+1,j,xdir) + qflx_momx(ke ,i,j ,xdir) ) &
1654  - j13g(ke-1,i,j,i_uyw) * ( qflx_momx(ke-1,i+1,j,xdir) + qflx_momx(ke-1 ,i,j ,xdir) ) &
1655  ) * mapf(i,j,1,i_uy) &
1656  + ( j23g(ke ,i,j,i_uyw) * ( qflx_momx(ke ,i ,j,ydir) + qflx_momx(ke-1+1,i,j-1,ydir) ) &
1657  - j23g(ke-1,i,j,i_uyw) * ( qflx_momx(ke-1,i ,j,ydir) + qflx_momx(ke-1 ,i,j-1,ydir) ) &
1658  ) * mapf(i,j,2,i_uy) &
1659  ) * 0.5_rp * rfdz(ke-1) &
1660  - j33g * ( qflx_momx(ke-1,i,j,zdir) ) * rcdz(ke) ) &
1661  / gsqrt(ke,i,j,i_uyz)
1662  enddo
1663  enddo
1664 
1665  return
1666  end subroutine atmos_phy_tb_calc_tend_momx
1667 
1668  subroutine atmos_phy_tb_calc_tend_momy( &
1669  MOMY_t_TB, &
1670  QFLX_MOMY, &
1671  GSQRT, J13G, J23G, J33G, MAPF, &
1672  IIS, IIE, JJS, JJE )
1673  use scale_grid, only: &
1674  rcdz => grid_rcdz, &
1675  rcdx => grid_rcdx, &
1676  rfdz => grid_rfdz, &
1677  rfdy => grid_rfdy, &
1678  fdz => grid_fdz
1679  use scale_gridtrans, only: &
1680  i_xyz, &
1681  i_xvw, &
1682  i_uvz, &
1683  i_xv
1684  implicit none
1685 
1686  real(RP), intent(out) :: MOMY_t_TB(ka,ia,ja)
1687  real(RP), intent(in) :: QFLX_MOMY(ka,ia,ja,3)
1688  real(RP), intent(in) :: GSQRT(ka,ia,ja,7)
1689  real(RP), intent(in) :: MAPF(ia,ja,2,4)
1690  real(RP), intent(in) :: J13G(ka,ia,ja,7)
1691  real(RP), intent(in) :: J23G(ka,ia,ja,7)
1692  real(RP), intent(in) :: J33G
1693  integer , intent(in) :: IIS
1694  integer , intent(in) :: IIE
1695  integer , intent(in) :: JJS
1696  integer , intent(in) :: JJE
1697 
1698  integer :: k, i, j
1699 
1700  do j = jjs, jje
1701  do i = iis, iie
1702  do k = ks+1, ke-1
1703  momy_t_tb(k,i,j) = &
1704  - ( ( gsqrt(k,i ,j ,i_uvz) * qflx_momy(k,i ,j,xdir) &
1705  - gsqrt(k,i-1,j ,i_uvz) * qflx_momy(k,i-1,j,xdir) ) * rcdx(i) * mapf(i,j,1,i_xv) &
1706  + ( gsqrt(k,i ,j+1,i_xyz) * qflx_momy(k,i,j+1,ydir) &
1707  - gsqrt(k,i ,j ,i_xyz) * qflx_momy(k,i,j ,ydir) ) * rfdy(j) * mapf(i,j,2,i_xv) &
1708  + ( ( j13g(k+1,i,j ,i_xvw) * ( qflx_momy(k+1,i,j ,xdir) + qflx_momy(k+1,i-1,j,xdir) ) &
1709  - j13g(k-1,i,j ,i_xvw) * ( qflx_momy(k-1,i,j ,xdir) + qflx_momy(k-1,i-1,j,xdir) ) &
1710  ) * mapf(i,j,1,i_xv) &
1711  + ( j23g(k+1,i,j+1,i_xvw) * ( qflx_momy(k+1,i,j+1,ydir) + qflx_momy(k+1,i ,j,ydir) ) &
1712  - j23g(k-1,i,j+1,i_xvw) * ( qflx_momy(k-1,i,j+1,ydir) + qflx_momy(k-1,i ,j,ydir) ) &
1713  ) * mapf(i,j,2,i_xv) &
1714  ) * 0.5_rp / ( fdz(k)+fdz(k-1) ) &
1715  + j33g * ( qflx_momy(k,i,j,zdir) - qflx_momy(k-1,i,j,zdir) ) * rcdz(k) ) &
1716  / gsqrt(k,i,j,i_xvw)
1717  enddo
1718  enddo
1719  enddo
1720  do j = jjs, jje
1721  do i = iis, iie
1722  momy_t_tb(ks,i,j) = &
1723  - ( ( gsqrt(ks,i ,j ,i_uvz) * qflx_momy(ks,i ,j,xdir) &
1724  - gsqrt(ks,i-1,j ,i_uvz) * qflx_momy(ks,i-1,j,xdir) ) * rcdx(i) * mapf(i,j,1,i_xv) &
1725  + ( gsqrt(ks,i ,j+1,i_xyz) * qflx_momy(ks,i,j+1,ydir) &
1726  - gsqrt(ks,i ,j ,i_xyz) * qflx_momy(ks,i,j ,ydir) ) * rfdy(j) * mapf(i,j,2,i_xv) &
1727  + ( ( j13g(ks+1,i,j ,i_xvw) * ( qflx_momy(ks+1,i,j ,xdir) + qflx_momy(ks+1,i-1,j,xdir) ) &
1728  - j13g(ks ,i,j ,i_xvw) * ( qflx_momy(ks ,i,j ,xdir) + qflx_momy(ks ,i-1,j,xdir) ) &
1729  ) * mapf(i,j,1,i_xv) &
1730  + ( j23g(ks+1,i,j+1,i_xvw) * ( qflx_momy(ks+1,i,j+1,ydir) + qflx_momy(ks+1,i ,j,ydir) ) &
1731  - j23g(ks ,i,j+1,i_xvw) * ( qflx_momy(ks ,i,j+1,ydir) + qflx_momy(ks ,i ,j,ydir) ) &
1732  ) * mapf(i,j,2,i_xv) &
1733  ) * 0.5_rp * rfdz(ks) &
1734  + j33g * ( qflx_momy(ks,i,j,zdir) ) * rcdz(ks) ) &
1735  / gsqrt(ks,i,j,i_xvw)
1736  momy_t_tb(ke,i,j) = &
1737  - ( ( gsqrt(ke,i ,j ,i_uvz) * qflx_momy(ke,i ,j,xdir) &
1738  - gsqrt(ke,i-1,j ,i_uvz) * qflx_momy(ke,i-1,j,xdir) ) * rcdx(i) * mapf(i,j,1,i_xv) &
1739  + ( gsqrt(ke,i ,j+1,i_xyz) * qflx_momy(ke,i,j+1,ydir) &
1740  - gsqrt(ke,i ,j ,i_xyz) * qflx_momy(ke,i,j ,ydir) ) * rfdy(j) * mapf(i,j,2,i_xv) &
1741  + ( ( j13g(ke ,i,j ,i_xvw) * ( qflx_momy(ke ,i,j ,xdir) + qflx_momy(ke ,i-1,j,xdir) ) &
1742  - j13g(ke-1,i,j ,i_xvw) * ( qflx_momy(ke-1,i,j ,xdir) + qflx_momy(ke-1,i-1,j,xdir) ) &
1743  ) * mapf(i,j,1,i_xv) &
1744  + ( j23g(ke ,i,j+1,i_xvw) * ( qflx_momy(ke ,i,j+1,ydir) + qflx_momy(ke ,i ,j,ydir) ) &
1745  - j23g(ke-1,i,j+1,i_xvw) * ( qflx_momy(ke-1,i,j+1,ydir) + qflx_momy(ke-1,i ,j,ydir) ) &
1746  ) * mapf(i,j,2,i_xv) &
1747  ) * 0.5_rp * rfdz(ke-1) &
1748  - j33g * ( qflx_momy(ke-1,i,j,zdir) ) * rcdz(ke) ) &
1749  / gsqrt(ke,i,j,i_xvw)
1750  end do
1751  end do
1752 
1753  return
1754  end subroutine atmos_phy_tb_calc_tend_momy
1755 
1756  subroutine atmos_phy_tb_calc_tend_phi( &
1757  phi_t_TB, &
1758  QFLX_phi, &
1759  GSQRT, J13G, J23G, J33G, MAPF, &
1760  IIS, IIE, JJS, JJE )
1761  use scale_grid, only: &
1762  rcdz => grid_rcdz, &
1763  rcdx => grid_rcdx, &
1764  rcdy => grid_rcdy, &
1765  rfdz => grid_rfdz, &
1766  fdz => grid_fdz
1767  use scale_gridtrans, only: &
1768  i_xyz, &
1769  i_xyw, &
1770  i_uyz, &
1771  i_xvz, &
1772  i_uvz, &
1773  i_xy
1774  implicit none
1775 
1776  real(RP), intent(out) :: phi_t_TB(ka,ia,ja)
1777  real(RP), intent(in) :: QFLX_phi(ka,ia,ja,3)
1778  real(RP), intent(in) :: GSQRT(ka,ia,ja,7)
1779  real(RP), intent(in) :: MAPF(ia,ja,2,4)
1780  real(RP), intent(in) :: J13G(ka,ia,ja,7)
1781  real(RP), intent(in) :: J23G(ka,ia,ja,7)
1782  real(RP), intent(in) :: J33G
1783  integer , intent(in) :: IIS
1784  integer , intent(in) :: IIE
1785  integer , intent(in) :: JJS
1786  integer , intent(in) :: JJE
1787 
1788  integer :: k, i, j
1789 
1790  do j = jjs, jje
1791  do i = iis, iie
1792  do k = ks+1, ke-1
1793  phi_t_tb(k,i,j) = &
1794  - ( ( gsqrt(k,i ,j,i_uyz) * qflx_phi(k,i ,j,xdir) &
1795  - gsqrt(k,i-1,j,i_uvz) * qflx_phi(k,i-1,j,xdir) ) * rcdx(i) * mapf(i,j,1,i_xy) &
1796  + ( gsqrt(k,i,j ,i_xvz) * qflx_phi(k,i,j ,ydir) &
1797  - gsqrt(k,i,j-1,i_xvz) * qflx_phi(k,i,j-1,ydir) ) * rcdy(j) * mapf(i,j,2,i_xy) &
1798  + ( ( j13g(k+1,i,j,i_xyw) * ( qflx_phi(k+1,i,j,xdir) + qflx_phi(k+1,i-1,j,xdir) ) &
1799  - j13g(k-1,i,j,i_xyw) * ( qflx_phi(k-1,i,j,xdir) + qflx_phi(k-1,i-1,j,xdir) ) &
1800  ) * mapf(i,j,1,i_xy) &
1801  + ( j23g(k+1,i,j,i_xyw) * ( qflx_phi(k+1,i,j,ydir) + qflx_phi(k+1,i,j-1,ydir) ) &
1802  - j23g(k-1,i,j,i_xyw) * ( qflx_phi(k-1,i,j,ydir) + qflx_phi(k-1,i,j-1,ydir) ) &
1803  ) * mapf(i,j,2,i_xy) &
1804  ) * 0.5_rp / ( fdz(k) + fdz(k-1) ) &
1805  + j33g * ( qflx_phi(k,i,j,zdir) - qflx_phi(k-1,i,j,zdir) ) * rcdz(k) ) &
1806  / gsqrt(k,i,j,i_xyz)
1807  enddo
1808  enddo
1809  enddo
1810  do j = jjs, jje
1811  do i = iis, iie
1812  phi_t_tb(ks,i,j) = &
1813  - ( ( gsqrt(ks,i ,j,i_uyz) * qflx_phi(ks,i ,j,xdir) &
1814  - gsqrt(ks,i-1,j,i_uvz) * qflx_phi(ks,i-1,j,xdir) ) * rcdx(i) * mapf(i,j,1,i_xy) &
1815  + ( gsqrt(ks,i,j ,i_xvz) * qflx_phi(ks,i,j ,ydir) &
1816  - gsqrt(ks,i,j-1,i_xvz) * qflx_phi(ks,i,j-1,ydir) ) * rcdy(j) * mapf(i,j,2,i_xy) &
1817  + ( ( j13g(ks+1,i,j,i_xyw) * ( qflx_phi(ks+1,i,j,xdir) + qflx_phi(ks+1,i-1,j,xdir) ) &
1818  - j13g(ks ,i,j,i_xyw) * ( qflx_phi(ks ,i,j,xdir) + qflx_phi(ks ,i-1,j,xdir) ) &
1819  ) * mapf(i,j,1,i_xy) &
1820  + ( j23g(ks+1,i,j,i_xyw) * ( qflx_phi(ks+1,i,j,ydir) + qflx_phi(ks+1,i,j-1,ydir) ) &
1821  - j23g(ks ,i,j,i_xyw) * ( qflx_phi(ks ,i,j,ydir) + qflx_phi(ks ,i,j-1,ydir) ) &
1822  ) * mapf(i,j,2,i_xy) &
1823  ) * 0.5_rp * rfdz(ks) &
1824  + j33g * ( qflx_phi(ks,i,j,zdir) ) * rcdz(ks) ) &
1825  / gsqrt(ks,i,j,i_xyz)
1826  phi_t_tb(ke,i,j) = &
1827  - ( ( gsqrt(ke,i ,j,i_uyz) * qflx_phi(ke,i ,j,xdir) &
1828  - gsqrt(ke,i-1,j,i_uvz) * qflx_phi(ke,i-1,j,xdir) ) * rcdx(i) * mapf(i,j,1,i_xy) &
1829  + ( gsqrt(ke,i,j ,i_xvz) * qflx_phi(ke,i,j ,ydir) &
1830  - gsqrt(ke,i,j-1,i_xvz) * qflx_phi(ke,i,j-1,ydir) ) * rcdy(j) * mapf(i,j,2,i_xy) &
1831  + ( ( j13g(ke ,i,j,i_xyw) * ( qflx_phi(ke ,i,j,xdir) + qflx_phi(ke ,i-1,j,xdir) ) &
1832  - j13g(ke-1,i,j,i_xyw) * ( qflx_phi(ke-1,i,j,xdir) + qflx_phi(ke-1,i-1,j,xdir) ) &
1833  ) * mapf(i,j,1,i_xy) &
1834  + ( j23g(ke ,i,j,i_xyw) * ( qflx_phi(ke ,i,j,ydir) + qflx_phi(ke ,i,j-1,ydir) ) &
1835  - j23g(ke-1,i,j,i_xyw) * ( qflx_phi(ke-1,i,j,ydir) + qflx_phi(ke-1,i,j-1,ydir) ) &
1836  ) * mapf(i,j,2,i_xy) &
1837  ) * 0.5_rp * rfdz(ke-1) &
1838  - j33g * ( qflx_phi(ke-1,i,j,zdir) ) * rcdz(ke) ) &
1839  / gsqrt(ke,i,j,i_xyz)
1840  end do
1841  end do
1842 
1843  return
1844  end subroutine atmos_phy_tb_calc_tend_phi
1845 
1846 end module scale_atmos_phy_tb_common
1847 
scale_grid_index::is
integer, public is
start point of inner domain: x, local
Definition: scale_grid_index.F90:60
scale_debug
module DEBUG
Definition: scale_debug.F90:13
scale_gridtrans::i_xvz
integer, public i_xvz
Definition: scale_gridtrans.F90:54
scale_atmos_phy_tb_common::atmos_phy_tb_calc_tend_momz
subroutine, public atmos_phy_tb_calc_tend_momz(MOMZ_t_TB, QFLX_MOMZ, GSQRT, J13G, J23G, J33G, MAPF, IIS, IIE, JJS, JJE)
Definition: scale_atmos_phy_tb_common.F90:1495
scale_grid_index::je
integer, public je
end point of inner domain: y, local
Definition: scale_grid_index.F90:63
scale_grid::grid_rcdy
real(rp), dimension(:), allocatable, public grid_rcdy
reciprocal of center-dy
Definition: scale_grid_cartesian.F90:61
scale_atmos_phy_tb_common::atmos_phy_tb_diffusion_solver
subroutine, public atmos_phy_tb_diffusion_solver(phi, a, b, c, d, KE_TB)
Definition: scale_atmos_phy_tb_common.F90:1453
scale_atmos_phy_tb_common::atmos_phy_tb_calc_flux_phi
subroutine, public atmos_phy_tb_calc_flux_phi(qflx_phi, DENS, PHI, Kh, FACT, GSQRT, J13G, J23G, J33G, MAPF, a, b, c, dt, implicit, IIS, IIE, JJS, JJE)
Definition: scale_atmos_phy_tb_common.F90:1198
scale_grid::grid_fdy
real(rp), dimension(:), allocatable, public grid_fdy
y-length of grid(j+1) to grid(j) [m]
Definition: scale_grid_cartesian.F90:68
scale_grid_index::iblock
integer, public iblock
block size for cache blocking: x
Definition: scale_grid_index.F90:47
scale_atmos_phy_tb_common::atmos_phy_tb_calc_tend_momy
subroutine, public atmos_phy_tb_calc_tend_momy(MOMY_t_TB, QFLX_MOMY, GSQRT, J13G, J23G, J33G, MAPF, IIS, IIE, JJS, JJE)
Definition: scale_atmos_phy_tb_common.F90:1673
scale_grid::grid_rcdx
real(rp), dimension(:), allocatable, public grid_rcdx
reciprocal of center-dx
Definition: scale_grid_cartesian.F90:60
scale_grid_index::zdir
integer, parameter, public zdir
Definition: scale_grid_index.F90:35
scale_stdio
module STDIO
Definition: scale_stdio.F90:12
scale_grid_index::ydir
integer, parameter, public ydir
Definition: scale_grid_index.F90:37
scale_grid_index::ke
integer, public ke
end point of inner domain: z, local
Definition: scale_grid_index.F90:59
scale_grid_index::xdir
integer, parameter, public xdir
Definition: scale_grid_index.F90:36
scale_gridtrans::i_xy
integer, public i_xy
Definition: scale_gridtrans.F90:57
scale_grid::grid_rfdy
real(rp), dimension(:), allocatable, public grid_rfdy
reciprocal of face-dy
Definition: scale_grid_cartesian.F90:71
scale_debug::check
subroutine, public check(current_line, v)
Undefined value checker.
Definition: scale_debug.F90:58
scale_grid::grid_rcdz
real(rp), dimension(:), allocatable, public grid_rcdz
reciprocal of center-dz
Definition: scale_grid_cartesian.F90:59
scale_gridtrans::i_xvw
integer, public i_xvw
Definition: scale_gridtrans.F90:52
scale_const::const_undef
real(rp), public const_undef
Definition: scale_const.F90:43
scale_grid_index
module grid index
Definition: scale_grid_index.F90:14
scale_atmos_phy_tb_common::atmos_phy_tb_calc_tend_phi
subroutine, public atmos_phy_tb_calc_tend_phi(phi_t_TB, QFLX_phi, GSQRT, J13G, J23G, J33G, MAPF, IIS, IIE, JJS, JJE)
Definition: scale_atmos_phy_tb_common.F90:1761
scale_tracer
module TRACER
Definition: scale_tracer.F90:18
scale_grid_index::ia
integer, public ia
of x whole cells (local, with HALO)
Definition: scale_grid_index.F90:55
scale_atmos_phy_tb_common::atmos_phy_tb_calc_tend_momx
subroutine, public atmos_phy_tb_calc_tend_momx(MOMX_t_TB, QFLX_MOMX, GSQRT, J13G, J23G, J33G, MAPF, IIS, IIE, JJS, JJE)
Definition: scale_atmos_phy_tb_common.F90:1585
scale_gridtrans
module GRIDTRANS
Definition: scale_gridtrans.F90:15
scale_grid::grid_fdz
real(rp), dimension(:), allocatable, public grid_fdz
z-length of grid(k+1) to grid(k) [m]
Definition: scale_grid_cartesian.F90:66
scale_grid_index::ka
integer, public ka
of z whole cells (local, with HALO)
Definition: scale_grid_index.F90:54
scale_gridtrans::i_uy
integer, public i_uy
Definition: scale_gridtrans.F90:58
scale_grid_index::jblock
integer, public jblock
block size for cache blocking: y
Definition: scale_grid_index.F90:48
scale_gridtrans::i_xyw
integer, public i_xyw
Definition: scale_gridtrans.F90:50
scale_grid_index::js
integer, public js
start point of inner domain: y, local
Definition: scale_grid_index.F90:62
scale_atmos_phy_tb_common::atmos_phy_tb_calc_strain_tensor
subroutine, public atmos_phy_tb_calc_strain_tensor(S33_C, S11_C, S22_C, S31_C, S12_C, S23_C, S12_Z, S23_X, S31_Y, S2, DENS, MOMZ, MOMX, MOMY, GSQRT, J13G, J23G, J33G, MAPF)
Definition: scale_atmos_phy_tb_common.F90:69
scale_const::const_undef2
integer, parameter, public const_undef2
undefined value (INT2)
Definition: scale_const.F90:40
scale_gridtrans::i_uyw
integer, public i_uyw
Definition: scale_gridtrans.F90:51
scale_const
module CONSTANT
Definition: scale_const.F90:14
scale_grid_index::ks
integer, public ks
start point of inner domain: z, local
Definition: scale_grid_index.F90:58
scale_atmos_phy_tb_common
module ATMOSPHERE / Physics Turbulence
Definition: scale_atmos_phy_tb_common.F90:14
scale_grid
module GRID (cartesian)
Definition: scale_grid_cartesian.F90:17
scale_gridtrans::i_xv
integer, public i_xv
Definition: scale_gridtrans.F90:59
scale_prof
module profiler
Definition: scale_prof.F90:10
scale_grid_index::ie
integer, public ie
end point of inner domain: x, local
Definition: scale_grid_index.F90:61
scale_gridtrans::i_uyz
integer, public i_uyz
Definition: scale_gridtrans.F90:53
scale_precision
module PRECISION
Definition: scale_precision.F90:13
scale_grid::grid_cdz
real(rp), dimension(:), allocatable, public grid_cdz
z-length of control volume [m]
Definition: scale_grid_cartesian.F90:56
scale_grid::grid_fdx
real(rp), dimension(:), allocatable, public grid_fdx
x-length of grid(i+1) to grid(i) [m]
Definition: scale_grid_cartesian.F90:67
scale_grid::grid_rfdx
real(rp), dimension(:), allocatable, public grid_rfdx
reciprocal of face-dx
Definition: scale_grid_cartesian.F90:70
scale_grid::grid_rfdz
real(rp), dimension(:), allocatable, public grid_rfdz
reciprocal of face-dz
Definition: scale_grid_cartesian.F90:69
scale_gridtrans::i_uvz
integer, public i_uvz
Definition: scale_gridtrans.F90:55
scale_gridtrans::i_xyz
integer, public i_xyz
Definition: scale_gridtrans.F90:49
scale_gridtrans::i_uv
integer, public i_uv
Definition: scale_gridtrans.F90:60
scale_grid_index::ja
integer, public ja
of y whole cells (local, with HALO)
Definition: scale_grid_index.F90:56
Generated by   doxygen 1.8.13