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