SCALE-RM
scale_atmos_phy_tb_smg.F90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
19 !-------------------------------------------------------------------------------
20 #include "scalelib.h"
21 module scale_atmos_phy_tb_smg
22  !-----------------------------------------------------------------------------
23  !
24  !++ used modules
25  !
26  use scale_precision
27  use scale_io
28  use scale_prof
29  use scale_atmos_grid_cartesc_index
30  use scale_tracer
31 
32 #if defined DEBUG || defined QUICKDEBUG
33  use scale_debug, only: &
34  check
35  use scale_const, only: &
36  undef => const_undef, &
37  iundef => const_undef2
38 #endif
39  !-----------------------------------------------------------------------------
40  implicit none
41  private
42  !-----------------------------------------------------------------------------
43  !
44  !++ Public procedure
45  !
46  public :: atmos_phy_tb_smg_setup
47  public :: atmos_phy_tb_smg_finalize
48  public :: atmos_phy_tb_smg
49 
50  !-----------------------------------------------------------------------------
51  !
52  !++ Public parameters & variables
53  !
54  !-----------------------------------------------------------------------------
55  !
56  !++ Private procedure
57  !
58  !-----------------------------------------------------------------------------
59  !
60  !++ Private parameters & variables
61  !
62  real(RP), private, parameter :: OneOverThree = 1.0_rp / 3.0_rp
63  real(RP), private, parameter :: twoOverThree = 2.0_rp / 3.0_rp
64  real(RP), private, parameter :: FourOverThree = 4.0_rp / 3.0_rp
65 
66  real(RP), private :: Cs = 0.13_rp ! Smagorinsky constant (Scotti et al. 1993)
67  real(RP), private, parameter :: PrN = 0.7_rp ! Prandtl number in neutral conditions
68  real(RP), private, parameter :: RiC = 0.25_rp ! critical Richardson number
69  real(RP), private, parameter :: FmC = 16.0_rp ! fum = sqrt(1 - c*Ri)
70  real(RP), private, parameter :: FhB = 40.0_rp ! fuh = sqrt(1 - b*Ri)/PrN
71  real(RP), private :: RPrN ! 1 / PrN
72  real(RP), private :: RRiC ! 1 / RiC
73  real(RP), private :: PrNovRiC ! PrN / RiC
74 
75  ! for backscatter
76  real(RP), private, parameter :: CB = 1.4_rp
77  real(RP), private, parameter :: CBt = 0.45_rp
78  real(RP), private, parameter :: aN = 0.47958315233127197_rp ! a_N = sqrt(0.23)
79  real(RP), private, parameter :: atN4 = 0.09_rp ! a_{\theta N}^4 = 0.3**2
80  real(RP), private, parameter :: C1o = an**3
81  real(RP), private, parameter :: D1o = prn * atn4 / an
82 
83  real(RP), private, allocatable :: lambda0(:,:,:)
84  real(RP), private, allocatable :: lambda (:,:,:)
85 
86  real(RP), private :: ATMOS_PHY_TB_SMG_NU_MAX = 10000.0_rp
87  logical, private :: ATMOS_PHY_TB_SMG_backscatter = .false.
88  logical, private :: ATMOS_PHY_TB_SMG_bottom = .true.
89  logical, private :: ATMOS_PHY_TB_SMG_implicit = .false.
90  logical, private :: ATMOS_PHY_TB_SMG_horizontal = .false.
91 
92  real(RP), private :: tke_fact
93 
94 
95  !-----------------------------------------------------------------------------
96 contains
97  !-----------------------------------------------------------------------------
99  subroutine atmos_phy_tb_smg_setup( &
100  FZ, CZ, CDX, CDY, MAPF, &
101  horizontal )
102  use scale_prc, only: &
103  prc_abort
104  use scale_const, only: &
105  karman => const_karman
106  implicit none
107 
108  real(rp), intent(in) :: fz (0:ka,ia,ja)
109  real(rp), intent(in) :: cz ( ka,ia,ja)
110  real(rp), intent(in) :: cdx (ia)
111  real(rp), intent(in) :: cdy (ja)
112  real(rp), intent(in) :: mapf(ia,ja,2)
113 
114  logical, intent(in), optional :: horizontal
115 
116  real(rp) :: atmos_phy_tb_smg_cs
117  real(rp) :: atmos_phy_tb_smg_filter_fact
118  logical :: atmos_phy_tb_smg_consistent_tke
119 
120  namelist / param_atmos_phy_tb_smg / &
121  atmos_phy_tb_smg_cs, &
122  atmos_phy_tb_smg_nu_max, &
123  atmos_phy_tb_smg_filter_fact, &
124  atmos_phy_tb_smg_bottom, &
125  atmos_phy_tb_smg_backscatter, &
126  atmos_phy_tb_smg_implicit, &
127  atmos_phy_tb_smg_consistent_tke, &
128  atmos_phy_tb_smg_horizontal
129 
130  integer :: ierr
131  integer :: k, i, j
132  !---------------------------------------------------------------------------
133 
134  !$acc data copyin(FZ, CZ, CDX, CDY, MAPF)
135 
136  log_newline
137  log_info("ATMOS_PHY_TB_smg_setup",*) 'Setup'
138  log_info("ATMOS_PHY_TB_smg_setup",*) 'Smagorinsky-type Eddy Viscocity Model'
139 
140  atmos_phy_tb_smg_cs = cs
141  atmos_phy_tb_smg_filter_fact = 2.0_rp
142  atmos_phy_tb_smg_consistent_tke = .true.
143 
144  if ( present(horizontal) ) atmos_phy_tb_smg_horizontal = horizontal
145  !--- read namelist
146  rewind(io_fid_conf)
147  read(io_fid_conf,nml=param_atmos_phy_tb_smg,iostat=ierr)
148  if( ierr < 0 ) then !--- missing
149  log_info("ATMOS_PHY_TB_smg_setup",*) 'Not found namelist. Default used.'
150  elseif( ierr > 0 ) then !--- fatal error
151  log_error("ATMOS_PHY_TB_smg_setup",*) 'Not appropriate names in namelist PARAM_ATMOS_PHY_TB_SMG. Check!'
152  call prc_abort
153  endif
154  log_nml(param_atmos_phy_tb_smg)
155 
156  cs = atmos_phy_tb_smg_cs
157 
158  rprn = 1.0_rp / prn
159  rric = 1.0_rp / ric
160  prnovric = ( 1.0_rp - prn ) * rric
161 
162  allocate( lambda0(ka,ia,ja) )
163  allocate( lambda(ka,ia,ja) )
164  !$acc enter data create(lambda0,lambda)
165 
166 #ifdef DEBUG
167  lambda0(:,:,:) = undef
168  lambda(:,:,:) = undef
169 #endif
170  if ( atmos_phy_tb_smg_horizontal ) then
171  !$omp parallel do collapse(2)
172  !$acc kernels
173  do j = js-1, je+1
174  do i = is-1, ie+1
175  do k = ks, ke
176  lambda0(k,i,j) = cs * sqrt( cdx(i) * cdy(j) / ( mapf(i,j,1) * mapf(i,j,2) ) )
177  lambda(k,i,j) = lambda0(k,i,j)
178  enddo
179  enddo
180  enddo
181  !$acc end kernels
182 
183 #ifdef DEBUG
184  i = iundef; j = iundef; k = iundef
185 #endif
186  atmos_phy_tb_smg_consistent_tke = .false.
187  atmos_phy_tb_smg_implicit = .false. ! flux in the z-direction is not necessary
188  atmos_phy_tb_smg_backscatter = .false.
189  else
190  !$omp parallel do collapse(2)
191  !$acc kernels
192  do j = js-1, je+1
193  do i = is-1, ie+1
194  do k = ks, ke
195  lambda0(k,i,j) = cs * mixlen( fz(k,i,j) - fz(k-1,i,j), &
196  cdx(i) / mapf(i,j,1), &
197  cdy(j) / mapf(i,j,2), &
198  atmos_phy_tb_smg_filter_fact )
199  enddo
200  enddo
201  enddo
202  !$acc end kernels
203 #ifdef DEBUG
204  i = iundef; j = iundef; k = iundef
205 #endif
206  if ( atmos_phy_tb_smg_bottom ) then
207  !$omp parallel do collapse(2)
208  !$acc kernels
209  do j = js-1, je+1
210  do i = is-1, ie+1
211  do k = ks, ke
212  lambda(k,i,j) = sqrt( 1.0_rp / ( 1.0_rp / lambda0(k,i,j)**2 + 1.0_rp / ( karman*(cz(k,i,j)-fz(ks-1,i,j)) )**2 ) )
213  enddo
214  enddo
215  enddo
216  !$acc end kernels
217 #ifdef DEBUG
218  i = iundef; j = iundef; k = iundef
219 #endif
220  else
221  !$omp parallel do collapse(2)
222  !$acc kernels
223  do j = js-1, je+1
224  do i = is-1, ie+1
225  do k = ks, ke
226  lambda(k,i,j) = lambda0(k,i,j)
227  enddo
228  enddo
229  enddo
230  !$acc end kernels
231 #ifdef DEBUG
232  i = iundef; j = iundef; k = iundef
233 #endif
234  end if
235  end if
236 
237  ! flag for isotropic stress tensor
238  if ( atmos_phy_tb_smg_consistent_tke ) then
239  tke_fact = 1.0_rp
240  else
241  tke_fact = 0.0_rp ! neglect
242  end if
243 
244  !$acc end data
245 
246  return
247  end subroutine atmos_phy_tb_smg_setup
248 
249  !-----------------------------------------------------------------------------
251  subroutine atmos_phy_tb_smg_finalize
252 
253  !$acc exit data delete(lambda0, lambda)
254  deallocate( lambda0 )
255  deallocate( lambda )
256 
257  return
258  end subroutine atmos_phy_tb_smg_finalize
259 
260  !-----------------------------------------------------------------------------
261  subroutine atmos_phy_tb_smg( &
262  qflx_sgs_momz, qflx_sgs_momx, qflx_sgs_momy, &
263  qflx_sgs_rhot, qflx_sgs_rhoq, &
264  MOMZ_t, MOMX_t, MOMY_t, RHOT_t, RHOQ_t, &
265  Nu, Ri, Pr, &
266  MOMZ, MOMX, MOMY, POTT, DENS, QTRC, N2, &
267  FZ, FDZ, RCDZ, RFDZ, CDX, FDX, CDY, FDY, &
268  GSQRT, J13G, J23G, J33G, MAPF, dt )
269  use scale_precision
270  use scale_atmos_grid_cartesc_index
271  use scale_tracer
272  use scale_const, only: &
273  eps => const_eps, &
274  grav => const_grav
275  use scale_file_history, only: &
276  file_history_in
277  use scale_atmos_phy_tb_common, only: &
278  calc_strain_tensor => atmos_phy_tb_calc_strain_tensor, &
279  calc_tend_momz => atmos_phy_tb_calc_tend_momz, &
280  calc_tend_momx => atmos_phy_tb_calc_tend_momx, &
281  calc_tend_momy => atmos_phy_tb_calc_tend_momy, &
282  calc_flux_phi => atmos_phy_tb_calc_flux_phi
283  use scale_atmos_hydrometeor, only: &
284  i_qv
285  use scale_random, only: &
286  random_normal
287  use scale_matrix, only: &
288  matrix_solver_tridiagonal
289  implicit none
290 
291  ! SGS flux
292  real(rp), intent(out) :: qflx_sgs_momz(ka,ia,ja,3)
293  real(rp), intent(out) :: qflx_sgs_momx(ka,ia,ja,3)
294  real(rp), intent(out) :: qflx_sgs_momy(ka,ia,ja,3)
295  real(rp), intent(out) :: qflx_sgs_rhot(ka,ia,ja,3)
296  real(rp), intent(out) :: qflx_sgs_rhoq(ka,ia,ja,3,qa)
297 
298  real(rp), intent(out) :: momz_t (ka,ia,ja)
299  real(rp), intent(out) :: momx_t (ka,ia,ja)
300  real(rp), intent(out) :: momy_t (ka,ia,ja)
301  real(rp), intent(out) :: rhot_t (ka,ia,ja)
302  real(rp), intent(out) :: rhoq_t (ka,ia,ja,qa) ! tendency of rho * QTRC
303 
304  real(rp), intent(out) :: nu (ka,ia,ja) ! eddy viscosity (center)
305  real(rp), intent(out) :: ri (ka,ia,ja) ! Richardson number
306  real(rp), intent(out) :: pr (ka,ia,ja) ! Prantle number
307 
308  real(rp), intent(in) :: momz (ka,ia,ja)
309  real(rp), intent(in) :: momx (ka,ia,ja)
310  real(rp), intent(in) :: momy (ka,ia,ja)
311  real(rp), intent(in) :: pott (ka,ia,ja)
312  real(rp), intent(in) :: dens (ka,ia,ja)
313  real(rp), intent(in) :: qtrc (ka,ia,ja,qa)
314  real(rp), intent(in) :: n2 (ka,ia,ja)
315 
316  real(rp), intent(in) :: fz (0:ka,ia,ja)
317  real(rp), intent(in) :: fdz (ka-1)
318  real(rp), intent(in) :: rcdz (ka)
319  real(rp), intent(in) :: rfdz (ka-1)
320  real(rp), intent(in) :: cdx (ia)
321  real(rp), intent(in) :: fdx (ia-1)
322  real(rp), intent(in) :: cdy (ja)
323  real(rp), intent(in) :: fdy (ja-1)
324 
325  real(rp), intent(in) :: gsqrt (ka,ia,ja,7)
326  real(rp), intent(in) :: j13g (ka,ia,ja,7)
327  real(rp), intent(in) :: j23g (ka,ia,ja,7)
328  real(rp), intent(in) :: j33g
329  real(rp), intent(in) :: mapf(ia,ja,2,4)
330  real(dp), intent(in) :: dt
331 
332  ! tke
333  real(rp) :: tke(ka,ia,ja)
334 
335  ! deformation rate tensor
336  real(rp) :: s33_c(ka,ia,ja) ! (cell center)
337  real(rp) :: s11_c(ka,ia,ja)
338  real(rp) :: s22_c(ka,ia,ja)
339  real(rp) :: s31_c(ka,ia,ja)
340  real(rp) :: s12_c(ka,ia,ja)
341  real(rp) :: s23_c(ka,ia,ja)
342  real(rp) :: s12_z(ka,ia,ja) ! (z edge or x-y plane)
343  real(rp) :: s23_x(ka,ia,ja) ! (x edge or y-z plane)
344  real(rp) :: s31_y(ka,ia,ja) ! (y edge or z-x plane)
345  real(rp) :: s2 (ka,ia,ja) ! |S|^2
346 
347  real(rp) :: kh(ka,ia,ja) ! eddy diffusion
348  real(rp) :: fm(ka)
349  real(rp) :: e(ka)
350  real(rp) :: et
351  real(rp) :: lambda_r(ka)
352  real(rp) :: rf
353  real(rp) :: c1(ka)
354  real(rp) :: c2
355  real(rp) :: d2
356 
357  ! implicit scheme
358  real(rp) :: tend(ka,ia,ja)
359  real(rp) :: tend2(ka,ia,ja)
360  real(rp) :: a (ka,ia,ja)
361  real(rp) :: b (ka,ia,ja)
362  real(rp) :: c (ka,ia,ja)
363  real(rp) :: ap
364 
365  ! backscatter
366  real(rp) :: random (ka,ia,ja)
367  real(rp) :: random_mz(ka,ia,ja)
368  real(rp) :: random_mx(ka,ia,ja)
369  real(rp) :: random_my(ka,ia,ja)
370  real(rp) :: random_qz(ka,ia,ja)
371  real(rp) :: random_qx(ka,ia,ja)
372  real(rp) :: random_qy(ka,ia,ja)
373  real(rp) :: dd (ka,ia,ja)
374  real(rp) :: leovleo5
375  real(rp) :: dz, dx, dy
376  real(rp) :: fact
377  real(rp) :: flxz(ka)
378 
379  integer :: iis, iie
380  integer :: jjs, jje
381 
382  integer :: k, i, j, iq
383  !---------------------------------------------------------------------------
384 
385  !$acc data copyout(qflx_sgs_momz, qflx_sgs_momx, qflx_sgs_momy, qflx_sgs_rhot, qflx_sgs_rhoq, &
386  !$acc MOMZ_t, MOMX_t, MOMY_t, RHOT_t, RHOQ_t, &
387  !$acc nu, Ri, Pr) &
388  !$acc copyin(MOMZ, MOMX, MOMY, POTT, DENS, QTRC, N2, &
389  !$acc FZ, FDZ, RCDZ, RFDZ, CDX, FDX, CDY, FDY, GSQRT, J13G, J23G, MAPF) &
390  !$acc create(TKE, S33_C, S11_C, S22_C, S31_C, S12_C, S23_C, S12_Z, S23_X, S31_Y, S2, &
391  !$acc Kh, TEND, TEND2, a, b, c)
392 
393  !$acc data create(random_mz, random_mx, random_my, random_qz, random_qx, random_qy, dd) if (ATMOS_PHY_TB_SMG_backscatter)
394 
395 
396  log_progress(*) 'atmosphere / physics / turbulence / Smagorinsky'
397 
398 #ifdef DEBUG
399  !$acc kernels
400  qflx_sgs_momz(:,:,:,:) = undef
401  qflx_sgs_momx(:,:,:,:) = undef
402  qflx_sgs_momy(:,:,:,:) = undef
403  qflx_sgs_rhot(:,:,:,:) = undef
404  qflx_sgs_rhoq(:,:,:,:,:) = undef
405 
406  nu(:,:,:) = undef
407  tke(:,:,:) = undef
408  pr(:,:,:) = undef
409  ri(:,:,:) = undef
410  kh(:,:,:) = undef
411  !$acc end kernels
412 #endif
413 
414 #ifdef QUICKDEBUG
415  !$acc kernels
416  qflx_sgs_momz(ks:ke, 1:is-1, : ,:) = undef
417  qflx_sgs_momz(ks:ke,ie+1:ia , : ,:) = undef
418  qflx_sgs_momz(ks:ke, : , 1:js-1,:) = undef
419  qflx_sgs_momz(ks:ke, : ,je+1:ja ,:) = undef
420  qflx_sgs_momx(ks:ke, 1:is-1, : ,:) = undef
421  qflx_sgs_momx(ks:ke,ie+1:ia , : ,:) = undef
422  qflx_sgs_momx(ks:ke, : , 1:js-1,:) = undef
423  qflx_sgs_momx(ks:ke, : ,je+1:ja ,:) = undef
424  qflx_sgs_momy(ks:ke, 1:is-1, : ,:) = undef
425  qflx_sgs_momy(ks:ke,ie+1:ia , : ,:) = undef
426  qflx_sgs_momy(ks:ke, : , 1:js-1,:) = undef
427  qflx_sgs_momy(ks:ke, : ,je+1:ja ,:) = undef
428  !$acc end kernels
429 #endif
430 
431 
432  !##### Start Upadate #####
433 
434  call calc_strain_tensor( &
435  s33_c, s11_c, s22_c, & ! (out)
436  s31_c, s12_c, s23_c, & ! (out)
437  s12_z, s23_x, s31_y, & ! (out)
438  s2 , & ! (out)
439  dens, momz, momx, momy, & ! (in)
440  gsqrt, j13g, j23g, j33g, mapf ) ! (in)
441 
442  if ( atmos_phy_tb_smg_backscatter ) then
443  call random_normal( random(:,:,:) )
444  ! 1:2:1 filter
445  !$omp parallel do collapse(2)
446  !$acc kernels copyin(random)
447  do j = js-1, je+1
448  do i = is-1, ie+1
449  do k = ks+1, ke-1
450  random_mz(k,i,j) = ( random(k,i,j) * 2.0_rp &
451  + random(k+1,i,j) + random(k-1,i,j) &
452  + random(k,i+1,j) + random(k,i-1,j) &
453  + random(k,i,j+1) + random(k,i,j-1) ) / 8.0_rp
454  end do
455  random_mz(ks,i,j) = ( random(ks,i,j) * 2.0_rp &
456  + random(ks+1,i,j) &
457  + random(ks,i+1,j) + random(ks,i-1,j) &
458  + random(ks,i,j+1) + random(ks,i,j-1) ) / 7.0_rp
459  random_mz(ke,i,j) = ( random(ke,i,j) * 2.0_rp &
460  + random(ke-1,i,j) &
461  + random(ke,i+1,j) + random(ke,i-1,j) &
462  + random(ke,i,j+1) + random(ke,i,j-1) ) / 7.0_rp
463  end do
464  end do
465  !$acc end kernels
466 
467  call random_normal( random(:,:,:) )
468  ! 1:2:1 filter
469  !$omp parallel do collapse(2)
470  !$acc kernels copyin(random)
471  do j = js-1, je+1
472  do i = is-1, ie+1
473  do k = ks+1, ke-1
474  random_mx(k,i,j) = ( random(k,i,j) * 2.0_rp &
475  + random(k+1,i,j) + random(k-1,i,j) &
476  + random(k,i+1,j) + random(k,i-1,j) &
477  + random(k,i,j+1) + random(k,i,j-1) ) / 8.0_rp
478  end do
479  random_mx(ks,i,j) = ( random(ks,i,j) * 2.0_rp &
480  + random(ks+1,i,j) &
481  + random(ks,i+1,j) + random(ks,i-1,j) &
482  + random(ks,i,j+1) + random(ks,i,j-1) ) / 7.0_rp
483  random_mx(ke,i,j) = ( random(ke,i,j) * 2.0_rp &
484  + random(ke-1,i,j) &
485  + random(ke,i+1,j) + random(ke,i-1,j) &
486  + random(ke,i,j+1) + random(ke,i,j-1) ) / 7.0_rp
487  end do
488  end do
489  !$acc end kernels
490 
491  call random_normal( random(:,:,:) )
492  ! 1:2:1 filter
493  !$omp parallel do collapse(2)
494  !$acc kernels copyin(random)
495  do j = js-1, je+1
496  do i = is-1, ie+1
497  do k = ks+1, ke-1
498  random_my(k,i,j) = ( random(k,i,j) * 2.0_rp &
499  + random(k+1,i,j) + random(k-1,i,j) &
500  + random(k,i+1,j) + random(k,i-1,j) &
501  + random(k,i,j+1) + random(k,i,j-1) ) / 8.0_rp
502  end do
503  random_my(ks,i,j) = ( random(ks,i,j) * 2.0_rp &
504  + random(ks+1,i,j) &
505  + random(ks,i+1,j) + random(ks,i-1,j) &
506  + random(ks,i,j+1) + random(ks,i,j-1) ) / 7.0_rp
507  random_my(ke,i,j) = ( random(ke,i,j) * 2.0_rp &
508  + random(ke-1,i,j) &
509  + random(ke,i+1,j) + random(ke,i-1,j) &
510  + random(ke,i,j+1) + random(ke,i,j-1) ) / 7.0_rp
511  end do
512  end do
513  !$acc end kernels
514 
515  call random_normal( random(:,:,:) )
516  ! 1:2:1 filter
517  !$omp parallel do collapse(2)
518  !$acc kernels copyin(random)
519  do j = js-1, je+1
520  do i = is-1, ie+1
521  do k = ks+1, ke-1
522  random_qz(k,i,j) = ( random(k,i,j) * 2.0_rp &
523  + random(k+1,i,j) + random(k-1,i,j) &
524  + random(k,i+1,j) + random(k,i-1,j) &
525  + random(k,i,j+1) + random(k,i,j-1) ) / 8.0_rp
526  end do
527  random_qz(ks,i,j) = ( random(ks,i,j) * 2.0_rp &
528  + random(ks+1,i,j) &
529  + random(ks,i+1,j) + random(ks,i-1,j) &
530  + random(ks,i,j+1) + random(ks,i,j-1) ) / 7.0_rp
531  random_qz(ke,i,j) = ( random(ke,i,j) * 2.0_rp &
532  + random(ke-1,i,j) &
533  + random(ke,i+1,j) + random(ke,i-1,j) &
534  + random(ke,i,j+1) + random(ke,i,j-1) ) / 7.0_rp
535  end do
536  end do
537  !$acc end kernels
538 
539  call random_normal( random(:,:,:) )
540  ! 1:2:1 filter
541  !$omp parallel do collapse(2)
542  !$acc kernels copyin(random)
543  do j = js-1, je+1
544  do i = is-1, ie+1
545  do k = ks+1, ke-1
546  random_qx(k,i,j) = ( random(k,i,j) * 2.0_rp &
547  + random(k+1,i,j) + random(k-1,i,j) &
548  + random(k,i+1,j) + random(k,i-1,j) &
549  + random(k,i,j+1) + random(k,i,j-1) ) / 8.0_rp
550  end do
551  random_qx(ks,i,j) = ( random(ks,i,j) * 2.0_rp &
552  + random(ks+1,i,j) &
553  + random(ks,i+1,j) + random(ks,i-1,j) &
554  + random(ks,i,j+1) + random(ks,i,j-1) ) / 7.0_rp
555  random_qx(ke,i,j) = ( random(ke,i,j) * 2.0_rp &
556  + random(ke-1,i,j) &
557  + random(ke,i+1,j) + random(ke,i-1,j) &
558  + random(ke,i,j+1) + random(ke,i,j-1) ) / 7.0_rp
559  end do
560  end do
561  !$acc end kernels
562 
563  call random_normal( random(:,:,:) )
564  ! 1:2:1 filter
565  !$omp parallel do collapse(2)
566  !$acc kernels copyin(random)
567  do j = js-1, je+1
568  do i = is-1, ie+1
569  do k = ks+1, ke-1
570  random_qy(k,i,j) = ( random(k,i,j) * 2.0_rp &
571  + random(k+1,i,j) + random(k-1,i,j) &
572  + random(k,i+1,j) + random(k,i-1,j) &
573  + random(k,i,j+1) + random(k,i,j-1) ) / 8.0_rp
574  end do
575  random_qy(ks,i,j) = ( random(ks,i,j) * 2.0_rp &
576  + random(ks+1,i,j) &
577  + random(ks,i+1,j) + random(ks,i-1,j) &
578  + random(ks,i,j+1) + random(ks,i,j-1) ) / 7.0_rp
579  random_qy(ke,i,j) = ( random(ke,i,j) * 2.0_rp &
580  + random(ke-1,i,j) &
581  + random(ke,i+1,j) + random(ke,i-1,j) &
582  + random(ke,i,j+1) + random(ke,i,j-1) ) / 7.0_rp
583  end do
584  end do
585  !$acc end kernels
586 
587  end if
588 
589 
590  !$omp parallel do collapse(2) &
591  !$omp private(fm,Rf,lambda_r,leOvleo5,C1,C2,D2,e,et,fact,dz,dx,dy)
592  !$acc kernels
593  do j = js-1, je+1
594  !$acc loop private(fm,lambda_r,C1,e,fact)
595  do i = is-1, ie+1
596  ! Ri = N^2 / |S|^2, N^2 = g / theta * dtheta/dz
597  do k = ks, ke
598  ri(k,i,j) = n2(k,i,j) / max(s2(k,i,j), eps)
599  enddo
600 
601  ! Nu
602  ! Pr = Nu / Kh = fm / fh
603  do k = ks, ke
604  if ( ri(k,i,j) < 0.0_rp ) then ! stable
605  fm(k) = sqrt( 1.0_rp - fmc * ri(k,i,j) )
606  nu(k,i,j) = lambda(k,i,j)**2 * sqrt( s2(k,i,j) ) * fm(k)
607  pr(k,i,j) = fm(k) / sqrt( 1.0_rp - fhb*ri(k,i,j) ) * prn
608  else if ( ri(k,i,j) < ric ) then ! weakly stable
609  fm(k) = ( 1.0_rp - ri(k,i,j)*rric )**4
610  nu(k,i,j) = lambda(k,i,j)**2 * sqrt( s2(k,i,j) ) * fm(k)
611  pr(k,i,j) = prn / ( 1.0_rp - prnovric * ri(k,i,j) )
612  else ! strongly stable
613  fm(k) = 0.0_rp
614  nu(k,i,j) = 0.0_rp
615  kh(k,i,j) = 0.0_rp
616  pr(k,i,j) = 1.0_rp
617  endif
618 
619  if ( ri(k,i,j) < ric ) then
620  kh(k,i,j) = max( min( nu(k,i,j) / pr(k,i,j), atmos_phy_tb_smg_nu_max ), eps )
621  nu(k,i,j) = max( min( nu(k,i,j), atmos_phy_tb_smg_nu_max ), eps )
622  pr(k,i,j) = nu(k,i,j) / kh(k,i,j)
623  rf = ri(k,i,j) / pr(k,i,j)
624  lambda_r(k) = lambda(k,i,j) * sqrt( fm(k) / sqrt( 1.0_rp - rf ) )
625  else
626  lambda_r(k) = 0.0_rp
627  end if
628 
629  enddo
630 
631 
632  if ( atmos_phy_tb_smg_backscatter ) then
633  do k = ks, ke
634  lambda_r(k) = min( 1.8_rp * lambda0(k,i,j), lambda_r(k) )
635  leovleo5 = ( lambda_r(k) / lambda0(k,i,j) )**5
636  c2 = cb * leovleo5 / ( 1.0_rp + cb * leovleo5 )
637  d2 = cbt * leovleo5 / ( 1.0_rp + cbt * leovleo5 )
638  c1(k) = c1o / sqrt( 1.0_rp - c2 )
639  ! D1 = D1o / sqrt( 1.0_RP - C2 )
640  e(k) = nu(k,i,j)**3 * ( 1.0_rp - c2 ) / ( lambda_r(k)**4 + eps )
641  et = kh(k,i,j) * ( 1.0_rp - d2 ) ! epsilon_t / D^2
642 
643  dz = fz(k,i,j) - fz(k-1,i,j)
644  dx = cdx(i) / mapf(i,j,1,i_xy)
645  dy = cdy(j) / mapf(i,j,2,i_xy)
646 
647  fact = sqrt( cb * leovleo5 * e(k) / ( 3.0_rp * dt ) ) * dens(k,i,j)
648  random_mz(k,i,j) = random_mz(k,i,j) * fact * dx * dy / sqrt( dx**2 + dy**2 )
649  random_mx(k,i,j) = random_mx(k,i,j) * fact * dy * dz / sqrt( dy**2 + dz**2 )
650  random_my(k,i,j) = random_my(k,i,j) * fact * dz * dx / sqrt( dz**2 + dx**2 )
651 
652  fact = sqrt( cbt * leovleo5 * et / ( 3.0_rp * dt ) ) * dens(k,i,j)
653  random_qz(k,i,j) = random_qz(k,i,j) * fact * dz
654  random_qx(k,i,j) = random_qx(k,i,j) * fact * dx
655  random_qy(k,i,j) = random_qy(k,i,j) * fact * dy
656  end do
657 
658  else
659 
660  do k = ks, ke
661  e(k) = nu(k,i,j)**3 / ( lambda_r(k)**4 + eps )
662  c1(k) = c1o
663  end do
664 
665  end if
666 
667  ! TKE
668  do k = ks, ke
669  tke(k,i,j) = ( e(k) * lambda_r(k) / c1(k) )**twooverthree
670  enddo
671 
672  enddo
673  enddo
674  !$acc end kernels
675 #ifdef DEBUG
676  i = iundef; j = iundef; k = iundef
677 #endif
678 
679  do jjs = js, je, jblock
680  jje = jjs+jblock-1
681  do iis = is, ie, iblock
682  iie = iis+iblock-1
683 
684  !##### momentum equation (z) #####
685 
686  !$omp parallel private(i,j,k)
687  ! (cell center)
688  if ( atmos_phy_tb_smg_horizontal ) then
689  !$omp workshare
690  !$acc kernels
691  qflx_sgs_momz(:,:,:,zdir) = 0.0_rp
692  !$acc end kernels
693  !$omp end workshare nowait
694  else
695  !$omp do collapse(2)
696  !$acc kernels
697  do j = jjs, jje
698  do i = iis, iie
699  do k = ks+1, ke-1
700 #ifdef DEBUG
701  call check( __line__, dens(k,i,j) )
702  call check( __line__, nu(k,i,j) )
703  call check( __line__, s33_c(k,i,j) )
704  call check( __line__, s11_c(k,i,j) )
705  call check( __line__, s22_c(k,i,j) )
706  call check( __line__, tke(k,i,j) )
707 #endif
708  qflx_sgs_momz(k,i,j,zdir) = dens(k,i,j) * ( &
709  - 2.0_rp * nu(k,i,j) &
710  * ( s33_c(k,i,j) - ( s11_c(k,i,j) + s22_c(k,i,j) + s33_c(k,i,j) ) * oneoverthree ) &
711  + twooverthree * tke(k,i,j) * tke_fact )
712  enddo
713  enddo
714  enddo
715  !$acc end kernels
716 #ifdef DEBUG
717  i = iundef; j = iundef; k = iundef
718 #endif
719  !$omp do
720  !$acc kernels
721  do j = jjs, jje
722  do i = iis, iie
723  ! momentum will not be conserved
724  qflx_sgs_momz(ks,i,j,zdir) = dens(ks,i,j) * twooverthree * tke(ks,i,j) * tke_fact
725  qflx_sgs_momz(ke,i,j,zdir) = dens(ke,i,j) * twooverthree * tke(ke,i,j) * tke_fact
726  ! anti-isotropic stress is calculated by the surface scheme
727  enddo
728  enddo
729  !$acc end kernels
730  !$omp end do nowait
731 #ifdef DEBUG
732  i = iundef; j = iundef; k = iundef
733 #endif
734  end if
735  ! (y edge)
736  !$omp do collapse(2)
737  !$acc kernels
738  do j = jjs, jje
739  do i = iis-1, iie
740  do k = ks, ke-1
741 #ifdef DEBUG
742  call check( __line__, dens(k,i,j) )
743  call check( __line__, dens(k,i+1,j) )
744  call check( __line__, dens(k+1,i,j) )
745  call check( __line__, dens(k+1,i+1,j) )
746  call check( __line__, nu(k,i,j) )
747  call check( __line__, nu(k,i+1,j) )
748  call check( __line__, nu(k+1,i,j) )
749  call check( __line__, nu(k+1,i+1,j) )
750  call check( __line__, s31_y(k,i,j) )
751 #endif
752  qflx_sgs_momz(k,i,j,xdir) = - 0.125_rp & ! 2.0 / 4 / 4
753  * ( dens(k,i,j)+dens(k+1,i,j)+dens(k,i+1,j)+dens(k+1,i+1,j) ) &
754  * ( nu(k,i,j)+nu(k+1,i,j)+nu(k,i+1,j)+nu(k+1,i+1,j)) &
755  * s31_y(k,i,j)
756  enddo
757  enddo
758  enddo
759  !$acc end kernels
760  !$omp end do nowait
761 #ifdef DEBUG
762  i = iundef; j = iundef; k = iundef
763 #endif
764  ! (x edge)
765  !$omp do collapse(2)
766  !$acc kernels
767  do j = jjs-1, jje
768  do i = iis, iie
769  do k = ks, ke-1
770 #ifdef DEBUG
771  call check( __line__, dens(k,i,j) )
772  call check( __line__, dens(k,i,j+1) )
773  call check( __line__, dens(k+1,i,j) )
774  call check( __line__, dens(k+1,i,j+1) )
775  call check( __line__, nu(k,i,j) )
776  call check( __line__, nu(k,i,j+1) )
777  call check( __line__, nu(k+1,i,j) )
778  call check( __line__, nu(k+1,i,j+1) )
779  call check( __line__, s23_x(k,i,j) )
780 #endif
781  qflx_sgs_momz(k,i,j,ydir) = - 0.125_rp & ! 2/4/4
782  * ( dens(k,i,j)+dens(k+1,i,j)+dens(k,i,j+1)+dens(k+1,i,j+1) ) &
783  * ( nu(k,i,j)+nu(k+1,i,j)+nu(k,i,j+1)+nu(k+1,i,j+1) ) &
784  * s23_x(k,i,j)
785  enddo
786  enddo
787  enddo
788  !$acc end kernels
789  !$omp end do nowait
790 #ifdef DEBUG
791  i = iundef; j = iundef; k = iundef
792 #endif
793 
794  !$omp end parallel
795 
796  if ( atmos_phy_tb_smg_implicit ) then
797 
798  call calc_tend_momz( tend, & ! (out)
799  qflx_sgs_momz, & ! (in)
800  gsqrt, j13g, j23g, j33g, mapf, & ! (in)
801  iis, iie, jjs, jje ) ! (in)
802 
803  !$omp parallel do collapse(2) &
804  !$omp private (ap)
805  !$acc kernels
806  do j = jjs, jje
807  do i = iis, iie
808  ap = - fouroverthree * dt &
809  * dens(ks+1,i,j)*nu(ks+1,i,j) &
810  * rcdz(ks+1) / gsqrt(ks+1,i,j,i_xyz)
811  a(ks,i,j) = ap * rfdz(ks) / gsqrt(ks,i,j,i_xyw)
812  c(ks,i,j) = 0.0_rp
813  b(ks,i,j) = - a(ks,i,j) + 0.5_rp * ( dens(ks,i,j)+dens(ks+1,i,j) )
814  do k = ks+1, ke-2
815 #ifdef _OPENACC
816  ap = - fouroverthree * dt &
817  * dens(k,i,j)*nu(k,i,j) &
818  * rcdz(k) / gsqrt(k,i,j,i_xyz)
819 #endif
820  c(k,i,j) = ap * rfdz(k+1) / gsqrt(k+1,i,j,i_xyw)
821  ap = - fouroverthree * dt &
822  * dens(k+1,i,j)*nu(k+1,i,j) &
823  * rcdz(k+1) / gsqrt(k+1,i,j,i_xyz)
824  a(k,i,j) = ap * rfdz(k) / gsqrt(k,i,j,i_xyw)
825  b(k,i,j) = - a(k,i,j) - c(k,i,j) + 0.5_rp * ( dens(k,i,j)+dens(k+1,i,j) )
826  end do
827 #ifdef _OPENACC
828  ap = - fouroverthree * dt &
829  * dens(ke-1,i,j)*nu(ke-1,i,j) &
830  * rcdz(ke-1) / gsqrt(ke-1,i,j,i_xyz)
831 #endif
832  a(ke-1,i,j) = 0.0_rp
833  c(ke-1,i,j) = ap * rfdz(ke) / gsqrt(ke,i,j,i_xyw)
834  b(ke-1,i,j) = - c(ke-1,i,j) + 0.5_rp * ( dens(ke-1,i,j)+dens(ke,i,j) )
835  end do
836  end do
837  !$acc end kernels
838 
839  call matrix_solver_tridiagonal( ka, ks, ke-1, ia, iis, iie, ja, jjs, jje, &
840  a(:,:,:), b(:,:,:), c(:,:,:), tend(:,:,:), &
841  tend2(:,:,:) )
842 
843  !$omp parallel do collapse(2)
844  !$acc kernels
845  do j = jjs, jje
846  do i = iis, iie
847  do k = ks+1, ke-1
848  qflx_sgs_momz(k,i,j,zdir) = qflx_sgs_momz(k,i,j,zdir) &
849  - fouroverthree * dens(k,i,j) * nu(k,i,j) * dt &
850  * ( tend2(k,i,j) - tend2(k-1,i,j) ) * rcdz(k) / gsqrt(k,i,j,i_xyz)
851  end do
852  end do
853  end do
854  !$acc end kernels
855 
856  end if
857 
858  !##### momentum equation (x) #####
859 
860  !$omp parallel private(i,j,k)
861 
862  ! (y edge)
863  if ( atmos_phy_tb_smg_horizontal ) then
864  !$omp workshare
865  !$acc kernels
866  qflx_sgs_momx(:,:,:,zdir) = 0.0_rp
867  !$acc end kernels
868  !$omp end workshare nowait
869  else
870  !$omp do collapse(2)
871  !$acc kernels
872  do j = jjs, jje
873  do i = iis, iie
874  do k = ks, ke-1
875 #ifdef DEBUG
876  call check( __line__, dens(k,i,j) )
877  call check( __line__, dens(k,i+1,j) )
878  call check( __line__, dens(k+1,i,j) )
879  call check( __line__, dens(k+1,i+1,j) )
880  call check( __line__, nu(k,i,j) )
881  call check( __line__, nu(k,i+1,j) )
882  call check( __line__, nu(k+1,i,j) )
883  call check( __line__, nu(k+1,i+1,j) )
884  call check( __line__, s31_y(k,i,j) )
885 #endif
886  qflx_sgs_momx(k,i,j,zdir) = - 0.125_rp & ! 2/4/4
887  * ( dens(k,i,j)+dens(k+1,i,j)+dens(k,i+1,j)+dens(k+1,i+1,j) ) &
888  * ( nu(k,i,j)+nu(k+1,i,j)+nu(k,i+1,j)+nu(k+1,i+1,j) ) &
889  * s31_y(k,i,j)
890  enddo
891  enddo
892  enddo
893  !$acc end kernels
894  !$omp end do nowait
895 #ifdef DEBUG
896  i = iundef; j = iundef; k = iundef
897 #endif
898  !$omp do
899  !$acc kernels
900  do j = jjs, jje
901  do i = iis, iie
902  qflx_sgs_momx(ks-1,i,j,zdir) = 0.0_rp ! bottom boundary
903  qflx_sgs_momx(ke ,i,j,zdir) = 0.0_rp ! top boundary
904  enddo
905  enddo
906  !$acc end kernels
907  !$omp end do nowait
908 #ifdef DEBUG
909  i = iundef; j = iundef; k = iundef
910 #endif
911  end if
912 
913  ! (cell center)
914  !$omp do collapse(2)
915  !$acc kernels
916  do j = jjs, jje
917  do i = iis, iie+1
918  do k = ks, ke
919 #ifdef DEBUG
920  call check( __line__, dens(k,i,j) )
921  call check( __line__, nu(k,i,j) )
922  call check( __line__, s11_c(k,i,j) )
923  call check( __line__, s22_c(k,i,j) )
924  call check( __line__, s33_c(k,i,j) )
925  call check( __line__, tke(k,i,j) )
926 #endif
927  qflx_sgs_momx(k,i,j,xdir) = dens(k,i,j) * ( &
928  - 2.0_rp * nu(k,i,j) &
929  * ( s11_c(k,i,j) - ( s11_c(k,i,j) + s22_c(k,i,j) + s33_c(k,i,j) ) * oneoverthree ) &
930  + twooverthree * tke(k,i,j) * tke_fact )
931  enddo
932  enddo
933  enddo
934  !$acc end kernels
935  !$omp end do nowait
936 #ifdef DEBUG
937  i = iundef; j = iundef; k = iundef
938 #endif
939 
940  ! (z edge)
941  !$omp do collapse(2)
942  !$acc kernels
943  do j = jjs-1, jje
944  do i = iis, iie
945  do k = ks, ke
946 #ifdef DEBUG
947  call check( __line__, dens(k,i,j) )
948  call check( __line__, dens(k,i+1,j) )
949  call check( __line__, dens(k,i,j+1) )
950  call check( __line__, dens(k,i+1,j+1) )
951  call check( __line__, nu(k,i,j) )
952  call check( __line__, nu(k,i+1,j) )
953  call check( __line__, nu(k,i,j+1) )
954  call check( __line__, nu(k,i+1,j+1) )
955  call check( __line__, s12_z(k,i,j) )
956 #endif
957  qflx_sgs_momx(k,i,j,ydir) = - 0.125_rp & ! 2/4/4
958  * ( dens(k,i,j)+dens(k,i+1,j)+dens(k,i,j+1)+dens(k,i+1,j+1) ) &
959  * ( nu(k,i,j)+nu(k,i+1,j)+nu(k,i,j+1)+nu(k,i+1,j+1) ) &
960  * s12_z(k,i,j)
961  enddo
962  enddo
963  enddo
964  !$acc end kernels
965  !$omp end do nowait
966 #ifdef DEBUG
967  i = iundef; j = iundef; k = iundef
968 #endif
969 
970  !$omp end parallel
971 
972  if ( atmos_phy_tb_smg_implicit ) then
973  call calc_tend_momx( tend, & ! (out)
974  qflx_sgs_momx, & ! (in)
975  gsqrt, j13g, j23g, j33g, mapf, & ! (in)
976  iis, iie, jjs, jje ) ! (in)
977 
978  !$omp parallel do collapse(2) &
979  !$omp private(ap)
980  !$acc kernels
981  do j = jjs, jje
982  do i = iis, iie
983 
984  ap = - dt * 0.25_rp * ( dens(ks ,i ,j)*nu(ks ,i ,j) &
985  + dens(ks+1,i ,j)*nu(ks+1,i ,j) &
986  + dens(ks ,i+1,j)*nu(ks ,i+1,j) &
987  + dens(ks+1,i+1,j)*nu(ks+1,i+1,j) ) &
988  * rfdz(ks) / gsqrt(ks,i,j,i_uyw)
989  a(ks,i,j) = ap * rcdz(ks) / gsqrt(ks,i,j,i_uyz)
990  c(ks,i,j) = 0.0_rp
991  b(ks,i,j) = - a(ks,i,j) + 0.5_rp * ( dens(ks,i,j)+dens(ks,i+1,j) )
992  do k = ks+1, ke-1
993 #ifdef _OPENACC
994  ap = - dt * 0.25_rp * ( dens(k-1,i ,j)*nu(k-1,i ,j) &
995  + dens(k ,i ,j)*nu(k ,i ,j) &
996  + dens(k-1,i+1,j)*nu(k-1,i+1,j) &
997  + dens(k ,i+1,j)*nu(k ,i+1,j) ) &
998  * rfdz(k-1) / gsqrt(k-1,i,j,i_uyw)
999 #endif
1000  c(k,i,j) = ap * rcdz(k) / gsqrt(k,i,j,i_uyz)
1001  ap = - dt * 0.25_rp * ( dens(k ,i ,j)*nu(k ,i ,j) &
1002  + dens(k+1,i ,j)*nu(k+1,i ,j) &
1003  + dens(k ,i+1,j)*nu(k ,i+1,j) &
1004  + dens(k+1,i+1,j)*nu(k+1,i+1,j) ) &
1005  * rfdz(k) / gsqrt(k,i,j,i_uyw)
1006  a(k,i,j) = ap * rcdz(k) / gsqrt(k,i,j,i_uyz)
1007  b(k,i,j) = - a(k,i,j) - c(k,i,j) + 0.5_rp * ( dens(k,i,j)+dens(k,i+1,j) )
1008  end do
1009 #ifdef _OPENACC
1010  ap = - dt * 0.25_rp * ( dens(ke-1,i ,j)*nu(ke-1,i ,j) &
1011  + dens(ke ,i ,j)*nu(ke ,i ,j) &
1012  + dens(ke-1,i+1,j)*nu(ke-1,i+1,j) &
1013  + dens(ke ,i+1,j)*nu(ke ,i+1,j) ) &
1014  * rfdz(ke-1) / gsqrt(ke-1,i,j,i_uyw)
1015 #endif
1016  a(ke,i,j) = 0.0_rp
1017  c(ke,i,j) = ap * rcdz(ke) / gsqrt(ke,i,j,i_uyz)
1018  b(ke,i,j) = - c(ke,i,j) + 0.5_rp * ( dens(ke,i,j)+dens(ke,i+1,j) )
1019  end do
1020  end do
1021  !$acc end kernels
1022 
1023  call matrix_solver_tridiagonal( ka, ks, ke, ia, iis, iie, ja, jjs, jje, &
1024  a(:,:,:), b(:,:,:), c(:,:,:), tend(:,:,:), &
1025  tend2(:,:,:) )
1026 
1027  !$omp parallel do collapse(2)
1028  !$acc kernels
1029  do j = jjs, jje
1030  do i = iis, iie
1031  do k = ks, ke-1
1032  qflx_sgs_momx(k,i,j,zdir) = qflx_sgs_momx(k,i,j,zdir) &
1033  - 0.25_rp * ( dens(k ,i ,j)*nu(k ,i ,j) &
1034  + dens(k+1,i ,j)*nu(k+1,i ,j) &
1035  + dens(k ,i+1,j)*nu(k ,i+1,j) &
1036  + dens(k+1,i+1,j)*nu(k+1,i+1,j) ) &
1037  * dt * ( tend2(k+1,i,j) - tend2(k,i,j) ) * rfdz(k) / gsqrt(k,i,j,i_uyw)
1038  end do
1039  end do
1040  end do
1041  !$acc end kernels
1042 
1043  end if
1044 
1045  !##### momentum equation (y) #####
1046  ! (x edge)
1047 
1048  !$omp parallel private(i,j,k)
1049 
1050  if ( atmos_phy_tb_smg_horizontal ) then
1051  !$omp workshare
1052  !$acc kernels
1053  qflx_sgs_momy(:,:,:,zdir) = 0.0_rp
1054  !$acc end kernels
1055  !$omp end workshare nowait
1056  else
1057  !$omp do collapse(2)
1058  !$acc kernels
1059  do j = jjs, jje
1060  do i = iis, iie
1061  do k = ks, ke-1
1062 #ifdef DEBUG
1063  call check( __line__, dens(k,i,j) )
1064  call check( __line__, dens(k,i,j+1) )
1065  call check( __line__, dens(k+1,i,j) )
1066  call check( __line__, dens(k+1,i,j+1) )
1067  call check( __line__, nu(k,i,j) )
1068  call check( __line__, nu(k,i,j+1) )
1069  call check( __line__, nu(k+1,i,j) )
1070  call check( __line__, nu(k+1,i,j+1) )
1071  call check( __line__, s23_x(k,i,j) )
1072 #endif
1073  qflx_sgs_momy(k,i,j,zdir) = - 0.125_rp & ! 2/4/4
1074  * ( dens(k,i,j)+dens(k+1,i,j)+dens(k,i,j+1)+dens(k+1,i,j+1) ) &
1075  * ( nu(k,i,j)+nu(k+1,i,j)+nu(k,i,j+1)+nu(k+1,i,j+1) ) &
1076  * s23_x(k,i,j)
1077  enddo
1078  enddo
1079  enddo
1080  !$acc end kernels
1081  !$omp end do nowait
1082 #ifdef DEBUG
1083  i = iundef; j = iundef; k = iundef
1084 #endif
1085  !$omp do
1086  !$acc kernels
1087  do j = jjs, jje
1088  do i = iis, iie
1089  qflx_sgs_momy(ks-1,i,j,zdir) = 0.0_rp ! bottom boundary
1090  qflx_sgs_momy(ke ,i,j,zdir) = 0.0_rp ! top boundary
1091  enddo
1092  enddo
1093  !$acc end kernels
1094  !$omp end do nowait
1095 #ifdef DEBUG
1096  i = iundef; j = iundef; k = iundef
1097 #endif
1098  end if
1099 
1100  ! (z edge)
1101  !$omp do collapse(2)
1102  !$acc kernels
1103  do j = jjs, jje
1104  do i = iis-1, iie
1105  do k = ks, ke
1106 #ifdef DEBUG
1107  call check( __line__, dens(k,i,j) )
1108  call check( __line__, dens(k,i+1,j) )
1109  call check( __line__, dens(k,i,j+1) )
1110  call check( __line__, dens(k,i+1,j+1) )
1111  call check( __line__, nu(k,i,j) )
1112  call check( __line__, nu(k,i+1,j) )
1113  call check( __line__, nu(k,i,j+1) )
1114  call check( __line__, nu(k,i+1,j+1) )
1115  call check( __line__, s12_z(k,i,j) )
1116 #endif
1117  qflx_sgs_momy(k,i,j,xdir) = - 0.125_rp & !
1118  * ( dens(k,i,j)+dens(k,i+1,j)+dens(k,i,j+1)+dens(k,i+1,j+1) ) &
1119  * ( nu(k,i,j)+nu(k,i+1,j)+nu(k,i,j+1)+nu(k,i+1,j+1) ) &
1120  * s12_z(k,i,j)
1121  enddo
1122  enddo
1123  enddo
1124  !$acc end kernels
1125  !$omp end do nowait
1126 #ifdef DEBUG
1127  i = iundef; j = iundef; k = iundef
1128 #endif
1129 
1130  ! (z-x plane)
1131  !$omp do collapse(2)
1132  !$acc kernels
1133  do j = jjs, jje+1
1134  do i = iis, iie
1135  do k = ks, ke
1136 #ifdef DEBUG
1137  call check( __line__, dens(k,i,j) )
1138  call check( __line__, nu(k,i,j) )
1139  call check( __line__, s11_c(k,i,j) )
1140  call check( __line__, s22_c(k,i,j) )
1141  call check( __line__, s33_c(k,i,j) )
1142  call check( __line__, tke(k,i,j) )
1143 #endif
1144  qflx_sgs_momy(k,i,j,ydir) = dens(k,i,j) * ( &
1145  - 2.0_rp * nu(k,i,j) &
1146  * ( s22_c(k,i,j) - ( s11_c(k,i,j) + s22_c(k,i,j) + s33_c(k,i,j) ) * oneoverthree ) &
1147  + twooverthree * tke(k,i,j) * tke_fact)
1148  enddo
1149  enddo
1150  enddo
1151  !$acc end kernels
1152  !$omp end do nowait
1153 #ifdef DEBUG
1154  i = iundef; j = iundef; k = iundef
1155 #endif
1156 
1157  !$omp end parallel
1158 
1159  if ( atmos_phy_tb_smg_implicit ) then
1160  call calc_tend_momy( tend, & ! (out)
1161  qflx_sgs_momy, & ! (in)
1162  gsqrt, j13g, j23g, j33g, mapf, & ! (in)
1163  iis, iie, jjs, jje ) ! (in)
1164 
1165  !$omp parallel do collapse(2) &
1166  !$omp private(ap)
1167  !$acc kernels
1168  do j = jjs, jje
1169  do i = iis, iie
1170 
1171  ap = - dt * 0.25_rp * ( dens(ks ,i,j )*nu(ks ,i,j ) &
1172  + dens(ks+1,i,j )*nu(ks+1,i,j ) &
1173  + dens(ks ,i,j+1)*nu(ks ,i,j+1) &
1174  + dens(ks+1,i,j+1)*nu(ks+1,i,j+1) ) &
1175  * rfdz(ks) / gsqrt(ks,i,j,i_xvw)
1176  a(ks,i,j) = ap * rcdz(ks) / gsqrt(ks,i,j,i_xvz)
1177  c(ks,i,j) = 0.0_rp
1178  b(ks,i,j) = - a(ks,i,j) + 0.5_rp * ( dens(ks,i,j)+dens(ks,i,j+1) )
1179  do k = ks+1, ke-1
1180 #ifdef _OPENACC
1181  ap = - dt * 0.25_rp * ( dens(k-1,i,j )*nu(k-1,i,j ) &
1182  + dens(k ,i,j )*nu(k ,i,j ) &
1183  + dens(k-1,i,j+1)*nu(k-1,i,j+1) &
1184  + dens(k ,i,j+1)*nu(k ,i,j+1) ) &
1185  * rfdz(k-1) / gsqrt(k-1,i,j,i_xvw)
1186 #endif
1187  c(k,i,j) = ap * rcdz(k) / gsqrt(k,i,j,i_xvz)
1188  ap = - dt * 0.25_rp * ( dens(k ,i,j )*nu(k ,i,j ) &
1189  + dens(k+1,i,j )*nu(k+1,i,j ) &
1190  + dens(k ,i,j+1)*nu(k ,i,j+1) &
1191  + dens(k+1,i,j+1)*nu(k+1,i,j+1) ) &
1192  * rfdz(k) / gsqrt(k,i,j,i_xvw)
1193  a(k,i,j) = ap * rcdz(k) / gsqrt(k,i,j,i_xvz)
1194  b(k,i,j) = - a(k,i,j) - c(k,i,j) + 0.5_rp * ( dens(k,i,j)+dens(k,i,j+1) )
1195  end do
1196 #ifdef _OPENACC
1197  ap = - dt * 0.25_rp * ( dens(ke-1,i,j )*nu(ke-1,i,j ) &
1198  + dens(ke ,i,j )*nu(ke ,i,j ) &
1199  + dens(ke-1,i,j+1)*nu(ke-1,i,j+1) &
1200  + dens(ke ,i,j+1)*nu(ke ,i,j+1) ) &
1201  * rfdz(ke-1) / gsqrt(ke-1,i,j,i_xvw)
1202 #endif
1203  a(ke,i,j) = 0.0_rp
1204  c(ke,i,j) = ap * rcdz(ke) / gsqrt(ke,i,j,i_xvz)
1205  b(ke,i,j) = - c(ke,i,j) + 0.5_rp * ( dens(ke,i,j)+dens(ke,i,j+1) )
1206  end do
1207  end do
1208  !$acc end kernels
1209 
1210  call matrix_solver_tridiagonal( ka, ks, ke, ia, iis, iie, ja, jjs, jje, &
1211  a(:,:,:), b(:,:,:), c(:,:,:), tend(:,:,:), &
1212  tend2(:,:,:) )
1213 
1214  !$omp parallel do collapse(2)
1215  !$acc kernels
1216  do j = jjs, jje
1217  do i = iis, iie
1218  do k = ks, ke-1
1219  qflx_sgs_momy(k,i,j,zdir) = qflx_sgs_momy(k,i,j,zdir) &
1220  - 0.25_rp * ( dens(k ,i,j )*nu(k ,i,j ) &
1221  + dens(k+1,i,j )*nu(k+1,i,j ) &
1222  + dens(k ,i,j+1)*nu(k ,i,j+1) &
1223  + dens(k+1,i,j+1)*nu(k+1,i,j+1) ) &
1224  * dt * ( tend2(k+1,i,j) - tend2(k,i,j) ) * rfdz(k) / gsqrt(k,i,j,i_xvw)
1225  end do
1226 
1227  end do
1228  end do
1229  !$acc end kernels
1230 
1231  end if
1232 
1233  if ( atmos_phy_tb_smg_backscatter ) then
1234 
1235 #define f2h(k,i,j,p) ( ( FZ(k+p-1,i,j) - FZ(k+p-2,i,j) ) / ( FZ(k+1,i,j) - FZ(k-1,i,j) ) )
1236 
1237  !$omp parallel private(flxz)
1238 
1239  ! MOMZ : dfy/dx - dfx/dy
1240  !$omp do collapse(2)
1241  !$acc kernels
1242  do j = jjs, jje
1243  !$acc loop private(flxz)
1244  do i = iis, iie
1245  do k = ks+1, ke-1
1246  flxz(k) = j13g(k,i,j,i_xyz) * random_my(k,i,j) &
1247  - j23g(k,i,j,i_xyz) * random_mx(k,i,j)
1248  end do
1249  flxz(ks) = 0.0_rp
1250  flxz(ke) = 0.0_rp
1251  do k = ks, ke-1
1252  momz_t(k,i,j) = ( ( gsqrt(k,i+1,j,i_xyw) * ( f2h(k,i+1,j,1) * random_my(k+1,i+1,j) &
1253  + f2h(k,i+1,j,2) * random_my(k,i+1,j) ) &
1254  - gsqrt(k,i-1,j,i_xyw) * ( f2h(k,i-1,j,1) * random_my(k+1,i-1,j) &
1255  + f2h(k,i-1,j,2) * random_my(k,i-1,j) ) &
1256  ) / ( fdx(i) + fdx(i-1) ) * mapf(i,j,1,i_xy) &
1257  - ( gsqrt(k,i,j+1,i_xyw) * ( f2h(k,i,j+1,1) * random_mx(k+1,i,j+1) &
1258  + f2h(k,i,j+1,2) * random_mx(k,i,j+1) ) &
1259  - gsqrt(k,i,j-1,i_xyw) * ( f2h(k,i,j-1,1) * random_mx(k+1,i,j-1) &
1260  + f2h(k,i,j+1,2) * random_mx(k,i,j-1) ) &
1261  ) / ( fdy(j) + fdy(j-1) ) * mapf(i,j,2,i_xy) &
1262  + ( flxz(k+1) - flxz(k) ) * rfdz(k) &
1263  ) / gsqrt(k,i,j,i_xyw)
1264  end do
1265  momz_t(ke,i,j) = 0.0_rp
1266  end do
1267  end do
1268  !$acc end kernels
1269  !$omp end do nowait
1270 
1271  ! MOMX : dfz/dy - dfy/dz
1272  !$omp do collapse(2)
1273  !$acc kernels
1274  do j = jjs, jje
1275  !$acc loop private(flxz)
1276  do i = iis, iie
1277  do k = ks, ke-1
1278  flxz(k) = j23g(k,i,j,i_uyw) * ( f2h(k,i+1,j,1) * random_mz(k+1,i+1,j) + f2h(k,i,j+1,2) * random_mz(k,i+1,j) &
1279  + f2h(k,i ,j,1) * random_mz(k+1,i ,j) + f2h(k,i ,j,2) * random_mz(k,i ,j) ) * 0.5_rp &
1280  - j33g * ( f2h(k,i+1,j,1) * random_my(k+1,i+1,j) + f2h(k,i+1,j,2) * random_my(k,i+1,j) &
1281  + f2h(k,i ,j,1) * random_my(k+1,i ,j) + f2h(k,i ,j,2) * random_my(k,i ,j) ) * 0.5_rp
1282  end do
1283  flxz(ks-1) = 0.0_rp
1284  flxz(ke ) = 0.0_rp
1285  do k = ks, ke
1286  momx_t(k,i,j) = ( ( gsqrt(k,i,j+1,i_uyz) * ( random_mz(k,i+1,j+1) + random_mz(k,i,j+1) ) * 0.5_rp &
1287  - gsqrt(k,i,j-1,i_uyz) * ( random_mz(k,i+1,j-1) + random_mz(k,i,j-1) ) * 0.5_rp ) / ( fdy(j) + fdy(j-1) ) * mapf(i,j,2,i_uy) &
1288  + ( flxz(k) - flxz(k-1) ) * rcdz(k) &
1289  ) / gsqrt(k,i,j,i_uyz)
1290  end do
1291  end do
1292  end do
1293  !$acc end kernels
1294  !$omp end do nowait
1295 
1296  ! MOMY : dfx/dz - dfz/dx
1297  !$omp do collapse(2)
1298  !$acc kernels
1299  do j = jjs, jje
1300  !$acc loop private(flxz)
1301  do i = iis, iie
1302  do k = ks, ke-1
1303  flxz(k) = j33g * ( f2h(k,i,j+1,1) * random_mx(k+1,i,j+1) + f2h(k,i,j+1,2) * random_mx(k,i,j+1) &
1304  + f2h(k,i,j ,1) * random_mx(k+1,i,j ) + f2h(k,i,j ,2) * random_mx(k,i,j ) ) * 0.5_rp &
1305  - j13g(k,i,j,i_xvw) * ( f2h(k,i,j+1,1) * random_mz(k+1,i,j+1) + f2h(k,i,j+1,2) * random_mz(k,i,j+1) &
1306  + f2h(k,i,j ,1) * random_mz(k+1,i,j ) + f2h(k,i,j ,2) * random_mz(k,i,j ) ) * 0.5_rp
1307  end do
1308  flxz(ks-1) = 0.0_rp
1309  flxz(ke ) = 0.0_rp
1310  do k = ks, ke
1311  momy_t(k,i,j) = ( ( flxz(k) - flxz(k-1) ) * rcdz(k) &
1312  - ( gsqrt(k,i+1,j,i_xvz) * ( random_mz(k,i+1,j+1) + random_mz(k,i+1,j) ) * 0.5_rp &
1313  - gsqrt(k,i-1,j,i_xyz) * ( random_mz(k,i-1,j+1) + random_mz(k,i-1,j) ) * 0.5_rp ) / ( fdx(i) + fdx(i-1) ) * mapf(i,j,1,i_xv) &
1314  ) / gsqrt(k,i,j,i_xvz)
1315  end do
1316  end do
1317  end do
1318  !$acc end kernels
1319  !$omp end do nowait
1320 
1321  !$omp end parallel
1322 
1323  else
1324 
1325  !$omp parallel
1326 
1327 !OCL XFILL
1328  !$omp do collapse(2)
1329  !$acc kernels
1330  do j = jjs, jje
1331  do i = iis, iie
1332  do k = ks, ke
1333  momz_t(k,i,j) = 0.0_rp
1334  end do
1335  end do
1336  end do
1337  !$acc end kernels
1338  !$omp end do nowait
1339 
1340 !OCL XFILL
1341  !$omp do collapse(2)
1342  !$acc kernels
1343  do j = jjs, jje
1344  do i = iis, iie
1345  do k = ks, ke
1346  momx_t(k,i,j) = 0.0_rp
1347  end do
1348  end do
1349  end do
1350  !$acc end kernels
1351  !$omp end do nowait
1352 
1353 !OCL XFILL
1354  !$omp do collapse(2)
1355  !$acc kernels
1356  do j = jjs, jje
1357  do i = iis, iie
1358  do k = ks, ke
1359  momy_t(k,i,j) = 0.0_rp
1360  end do
1361  end do
1362  end do
1363  !$acc end kernels
1364  !$omp end do nowait
1365 
1366  !$omp end parallel
1367 
1368  end if
1369 
1370  !##### Thermodynamic Equation #####
1371 
1372  if ( atmos_phy_tb_smg_implicit ) then
1373 
1374  !$omp parallel do collapse(2) &
1375  !$omp private (ap)
1376  !$acc kernels
1377  do j = jjs, jje
1378  do i = iis, iie
1379  ap = - dt * 0.25_rp * ( dens(ks,i,j)+dens(ks+1,i,j) ) &
1380  * ( kh(ks,i,j)+kh(ks+1,i,j) ) &
1381  * rfdz(ks) / gsqrt(ks,i,j,i_xyw)
1382  a(ks,i,j) = ap * rcdz(ks) / gsqrt(ks,i,j,i_xyz)
1383  c(ks,i,j) = 0.0_rp
1384  b(ks,i,j) = - a(ks,i,j) + dens(ks,i,j)
1385  do k = ks+1, ke-1
1386 #ifdef _OPENACC
1387  ap = - dt * 0.25_rp * ( dens(k-1,i,j)+dens(k,i,j) ) &
1388  * ( kh(k-1,i,j)+kh(k,i,j) ) &
1389  * rfdz(k-1) / gsqrt(k-1,i,j,i_xyw)
1390 #endif
1391  c(k,i,j) = ap * rcdz(k) / gsqrt(k,i,j,i_xyz)
1392  ap = - dt * 0.25_rp * ( dens(k,i,j)+dens(k+1,i,j) ) &
1393  * ( kh(k,i,j)+kh(k+1,i,j) ) &
1394  * rfdz(k) / gsqrt(k,i,j,i_xyw)
1395  a(k,i,j) = ap * rcdz(k) / gsqrt(k,i,j,i_xyz)
1396  b(k,i,j) = - a(k,i,j) - c(k,i,j) + dens(k,i,j)
1397  end do
1398 #ifdef _OPENACC
1399  ap = - dt * 0.25_rp * ( dens(ke-1,i,j)+dens(ke,i,j) ) &
1400  * ( kh(ke-1,i,j)+kh(ke,i,j) ) &
1401  * rfdz(ke-1) / gsqrt(ke-1,i,j,i_xyw)
1402 #endif
1403  a(ke,i,j) = 0.0_rp
1404  c(ke,i,j) = ap * rcdz(ke) / gsqrt(ke,i,j,i_xyz)
1405  b(ke,i,j) = - c(ke,i,j) + dens(ke,i,j)
1406 
1407  end do
1408  end do
1409  !$acc end kernels
1410 
1411  end if
1412 
1413  call calc_flux_phi( &
1414  qflx_sgs_rhot, &
1415  dens, pott, kh, 1.0_rp, &
1416  gsqrt, j13g, j23g, j33g, mapf, &
1417  atmos_phy_tb_smg_horizontal, &
1418  atmos_phy_tb_smg_implicit, &
1419  a, b, c, dt, &
1420  iis, iie, jjs, jje )
1421 
1422  if ( atmos_phy_tb_smg_backscatter ) then
1423 
1424  !$omp parallel do collapse(2)
1425  !$acc kernels
1426  do j = jjs-1, jje+1
1427  do i = iis-1, iie+1
1428  do k = ks+1, ke-1
1429  dd(k,i,j) = sqrt( ( ( pott(k+1,i,j) - pott(k-1,i,j) ) * j33g / ( fdz(k) + fdz(k-1) ) )**2 &
1430  + ( ( ( gsqrt(k,i+1,j,i_xyz)*pott(k,i+1,j) - gsqrt(k,i-1,j,i_xyz)*pott(k,i-1,j) ) / ( fdx(i) + fdx(i-1) ) &
1431  + ( j13g(k+1,i,j,i_xyz)*pott(k+1,i,j) - j13g(k-1,i,j,i_xyz)*pott(k-1,i,j) ) / ( fdz(k) + fdz(k-1) ) ) * mapf(i,j,1,i_xy) )**2 &
1432  + ( ( ( gsqrt(k,i,j+1,i_xyz)*pott(k,i,j+1) - gsqrt(k,i,j-1,i_xyz)*pott(k,i,j-1) ) / ( fdy(j) + fdy(j-1) ) &
1433  + ( j23g(k+1,i,j,i_xyz)*pott(k+1,i,j) - j23g(k-1,i,j,i_xyz)*pott(k-1,i,j) ) / ( fdz(k) + fdz(k-1) ) ) * mapf(i,j,2,i_xy) )**2 &
1434  ) / gsqrt(k,i,j,i_xyz)
1435  end do
1436  dd(ks,i,j) = sqrt( ( ( pott(ks+1,i,j) - pott(ks,i,j) ) * j33g * rfdz(ks) )**2 &
1437  + ( ( ( gsqrt(ks,i+1,j,i_xyz)*pott(ks,i+1,j) - gsqrt(ks,i-1,j,i_xyz)*pott(ks,i-1,j) ) / ( fdx(i) + fdx(i-1) ) &
1438  + ( j13g(ks+1,i,j,i_xyz)*pott(ks+1,i,j) - j13g(ks,i,j,i_xyz)*pott(ks,i,j) ) * rfdz(ks) ) * mapf(i,j,1,i_xy) )**2 &
1439  + ( ( ( gsqrt(ks,i,j+1,i_xyz)*pott(ks,i,j+1) - gsqrt(ks,i,j-1,i_xyz)*pott(ks,i,j-1) ) / ( fdy(j) + fdy(j-1) ) &
1440  + ( j23g(ks+1,i,j,i_xyz)*pott(ks+1,i,j) - j23g(ks,i,j,i_xyz)*pott(ks,i,j) ) * rfdz(ks) ) * mapf(i,j,2,i_xy) )**2 &
1441  ) / gsqrt(ks,i,j,i_xyz)
1442  dd(ke,i,j) = sqrt( ( ( pott(ke,i,j) - pott(ke-1,i,j) ) * j33g * rfdz(ke-1) )**2 &
1443  + ( ( ( gsqrt(ke,i+1,j,i_xyz)*pott(ke,i+1,j) - gsqrt(ke,i-1,j,i_xyz)*pott(ke,i-1,j) ) / ( fdx(i) + fdx(i-1) ) &
1444  + ( j13g(ke,i,j,i_xyz)*pott(ke,i,j) - j13g(ke-1,i,j,i_xyz)*pott(ke-1,i,j) ) * rfdz(ke-1) ) * mapf(i,j,1,i_xy) )**2 &
1445  + ( ( ( gsqrt(ke,i,j+1,i_xyz)*pott(ke,i,j+1) - gsqrt(ke,i,j-1,i_xyz)*pott(ke,i,j-1) ) / ( fdy(j) + fdy(j-1) ) &
1446  + ( j23g(ke,i,j,i_xyz)*pott(ke,i,j) - j23g(ke-1,i,j,i_xyz)*pott(ke-1,i,j) ) * rfdz(ke-1) ) * mapf(i,j,2,i_xy) )**2 &
1447  ) / gsqrt(ke,i,j,i_xyz)
1448 
1449  end do
1450  end do
1451  !$acc end kernels
1452 
1453  !$omp parallel do collapse(2) &
1454  !$omp private(flxz)
1455  !$acc kernels
1456  do j = jjs, jje
1457  !$acc loop private(flxz)
1458  do i = iis, iie
1459  do k = ks, ke-1
1460  flxz(k) = j33g * ( f2h(k,i,j,1) * random_qz(k+1,i,j) * dd(k+1,i,j) + f2h(k,i,j,2) * random_qz(k,i,j) * dd(k,i,j) ) &
1461  + j13g(k,i,j,i_xyw) * ( f2h(k,i,j,1) * random_qx(k+1,i,j) * dd(k+1,i,j) + f2h(k,i,j,2) * random_qx(k,i,j) * dd(k,i,j) ) &
1462  + j23g(k,i,j,i_xyw) * ( f2h(k,i,j,1) * random_qy(k+1,i,j) * dd(k+1,i,j) + f2h(k,i,j,2) * random_qy(k,i,j) * dd(k,i,j) )
1463  end do
1464  flxz(ks-1) = 0.0_rp
1465  flxz(ke ) = 0.0_rp
1466  do k = ks, ke
1467  rhot_t(k,i,j) = ( ( flxz(k) - flxz(k-1) ) * rcdz(k) &
1468  + ( gsqrt(k,i+1,j,i_xyz) * random_qx(k,i+1,j) * dd(k,i+1,j) - gsqrt(k,i-1,j,i_xyz) * random_qx(k,i-1,j) * dd(k,i-1,j) ) / ( fdx(i) + fdx(i-1) ) * mapf(i,j,1,i_xy) &
1469  + ( gsqrt(k,i,j+1,i_xyz) * random_qy(k,i,j+1) * dd(k,i,j+1) - gsqrt(k,i,j-1,i_xyz) * random_qy(k,i,j-1) * dd(k,i,j-1) ) / ( fdy(j) + fdy(j-1) ) * mapf(i,j,2,i_xy) &
1470  ) / gsqrt(k,i,j,i_xyz)
1471  end do
1472  end do
1473  end do
1474  !$acc end kernels
1475 
1476  else
1477 
1478  !$omp parallel do collapse(2)
1479  !$acc kernels
1480 !OCL XFILL
1481  do j = jjs, jje
1482  do i = iis, iie
1483  do k = ks, ke
1484  rhot_t(k,i,j) = 0.0_rp
1485  end do
1486  end do
1487  end do
1488  !$acc end kernels
1489 
1490  end if
1491 
1492  enddo
1493  enddo
1494 
1495 
1496  !##### Tracers #####
1497  do iq = 1, qa
1498 
1499  if ( .not. tracer_advc(iq) ) cycle
1500 
1501  do jjs = js, je, jblock
1502  jje = jjs+jblock-1
1503  do iis = is, ie, iblock
1504  iie = iis+iblock-1
1505 
1506  call calc_flux_phi( &
1507  qflx_sgs_rhoq(:,:,:,:,iq), &
1508  dens, qtrc(:,:,:,iq), kh, 1.0_rp, &
1509  gsqrt, j13g, j23g, j33g, mapf, &
1510  atmos_phy_tb_smg_horizontal, &
1511  atmos_phy_tb_smg_implicit, &
1512  a, b, c, dt, &
1513  iis, iie, jjs, jje )
1514 
1515 
1516  if ( atmos_phy_tb_smg_backscatter .and. iq == i_qv ) then
1517 
1518  !$omp parallel do collapse(2)
1519  !$acc kernels
1520  do j = jjs-1, jje+1
1521  do i = iis-1, iie+1
1522  do k = ks+1, ke-1
1523  dd(k,i,j) = sqrt( ( ( qtrc(k+1,i,j,iq) - qtrc(k-1,i,j,iq) ) * j33g / ( fdz(k) + fdz(k-1) ) )**2 &
1524  + ( ( ( gsqrt(k,i+1,j,i_xyz)*qtrc(k,i+1,j,iq) - gsqrt(k,i-1,j,i_xyz)*qtrc(k,i-1,j,iq) ) / ( fdx(i) + fdx(i-1) ) &
1525  + ( j13g(k+1,i,j,i_xyz)*qtrc(k+1,i,j,iq) - j13g(k-1,i,j,i_xyz)*qtrc(k-1,i,j,iq) ) / ( fdz(k) + fdz(k-1) ) ) * mapf(i,j,1,i_xy) )**2 &
1526  + ( ( ( gsqrt(k,i,j+1,i_xyz)*qtrc(k,i,j+1,iq) - gsqrt(k,i,j-1,i_xyz)*qtrc(k,i,j-1,iq) ) / ( fdy(j) + fdy(j-1) ) &
1527  + ( j23g(k+1,i,j,i_xyz)*qtrc(k+1,i,j,iq) - j23g(k-1,i,j,i_xyz)*qtrc(k-1,i,j,iq) ) / ( fdz(k) + fdz(k-1) ) ) * mapf(i,j,2,i_xy) )**2 &
1528  ) / gsqrt(k,i,j,i_xyz)
1529  end do
1530  dd(ks,i,j) = sqrt( ( ( qtrc(ks+1,i,j,iq) - qtrc(ks,i,j,iq) ) * j33g * rfdz(ks) )**2 &
1531  + ( ( ( gsqrt(ks,i+1,j,i_xyz)*qtrc(ks,i+1,j,iq) - gsqrt(ks,i-1,j,i_xyz)*qtrc(ks,i-1,j,iq) ) / ( fdx(i) + fdx(i-1) ) &
1532  + ( j13g(ks+1,i,j,i_xyz)*qtrc(ks+1,i,j,iq) - j13g(ks,i,j,i_xyz)*qtrc(ks,i,j,iq) ) * rfdz(ks) ) * mapf(i,j,1,i_xy) )**2 &
1533  + ( ( ( gsqrt(ks,i,j+1,i_xyz)*qtrc(ks,i,j+1,iq) - gsqrt(ks,i,j-1,i_xyz)*qtrc(ks,i,j-1,iq) ) / ( fdy(j) + fdy(j-1) ) &
1534  + ( j23g(ks+1,i,j,i_xyz)*qtrc(ks+1,i,j,iq) - j23g(ks,i,j,i_xyz)*qtrc(ks,i,j,iq) ) * rfdz(ks) ) * mapf(i,j,2,i_xy) )**2 &
1535  ) / gsqrt(ks,i,j,i_xyz)
1536  dd(ke,i,j) = sqrt( ( ( qtrc(ke,i,j,iq) - qtrc(ke-1,i,j,iq) ) * j33g * rfdz(ke-1) )**2 &
1537  + ( ( ( gsqrt(ke,i+1,j,i_xyz)*qtrc(ke,i+1,j,iq) - gsqrt(ke,i-1,j,i_xyz)*qtrc(ke,i-1,j,iq) ) / ( fdx(i) + fdx(i-1) ) &
1538  + ( j13g(ke,i,j,i_xyz)*qtrc(ke,i,j,iq) - j13g(ke-1,i,j,i_xyz)*qtrc(ke-1,i,j,iq) ) * rfdz(ke-1) ) * mapf(i,j,1,i_xy) )**2 &
1539  + ( ( ( gsqrt(ke,i,j+1,i_xyz)*qtrc(ke,i,j+1,iq) - gsqrt(ke,i,j-1,i_xyz)*qtrc(ke,i,j-1,iq) ) / ( fdy(j) + fdy(j-1) ) &
1540  + ( j23g(ke,i,j,i_xyz)*qtrc(ke,i,j,iq) - j23g(ke-1,i,j,i_xyz)*qtrc(ke-1,i,j,iq) ) * rfdz(ke-1) ) * mapf(i,j,2,i_xy) )**2 &
1541  ) / gsqrt(ke,i,j,i_xyz)
1542 
1543  end do
1544  end do
1545  !$acc end kernels
1546 
1547  !$omp parallel do collapse(2) &
1548  !$omp private (flxz)
1549  !$acc kernels
1550  do j = jjs, jje
1551  !$acc loop private(flxz)
1552  do i = iis, iie
1553  do k = ks, ke-1
1554  flxz(k) = j33g * ( f2h(k,i,j,1) * random_qz(k+1,i,j) * dd(k+1,i,j) + f2h(k,i,j,2) * random_qz(k,i,j) * dd(k,i,j) ) &
1555  + j13g(k,i,j,i_xyw) * ( f2h(k,i,j,1) * random_qx(k+1,i,j) * dd(k+1,i,j) + f2h(k,i,j,2) * random_qx(k,i,j) * dd(k,i,j) ) &
1556  + j23g(k,i,j,i_xyw) * ( f2h(k,i,j,1) * random_qy(k+1,i,j) * dd(k+1,i,j) + f2h(k,i,j,2) * random_qy(k,i,j) * dd(k,i,j) )
1557  end do
1558  flxz(ks-1) = 0.0_rp
1559  flxz(ke ) = 0.0_rp
1560  do k = ks, ke
1561  rhoq_t(k,i,j,iq) = ( ( flxz(k) - flxz(k-1) ) * rcdz(k) &
1562  + ( gsqrt(k,i+1,j,i_xyz) * random_qx(k,i+1,j) * dd(k,i+1,j) - gsqrt(k,i-1,j,i_xyz) * random_qx(k,i-1,j) * dd(k,i-1,j) ) / ( fdx(i) + fdx(i-1) ) * mapf(i,j,1,i_xy) &
1563  + ( gsqrt(k,i,j+1,i_xyz) * random_qy(k,i,j+1) * dd(k,i,j+1) - gsqrt(k,i,j-1,i_xyz) * random_qy(k,i,j-1) * dd(k,i,j-1) ) / ( fdy(j) + fdy(j-1) ) * mapf(i,j,2,i_xy) &
1564  ) / gsqrt(k,i,j,i_xyz)
1565  end do
1566  end do
1567  end do
1568  !$acc end kernels
1569 
1570  else
1571 
1572  !$omp parallel do collapse(2)
1573  !$acc kernels
1574 !OCL XFILL
1575  do j = jjs, jje
1576  do i = iis, iie
1577  do k = ks, ke
1578  rhoq_t(k,i,j,iq) = 0.0_rp
1579  end do
1580  end do
1581  end do
1582  !$acc end kernels
1583 
1584  end if
1585 
1586  enddo
1587  enddo
1588 
1589  enddo ! scalar quantities loop
1590 #ifdef DEBUG
1591  iq = iundef
1592 #endif
1593 
1594  call file_history_in( tke(:,:,:), 'TKE_SMG', 'turbulent kinetic energy (Smagorinsky)', 'm2/s2', fill_halo=.true. )
1595 
1596  !$acc end data
1597  !$acc end data
1598 
1599  return
1600  end subroutine atmos_phy_tb_smg
1601 
1602 
1603  function mixlen(dz, dx, dy, filter_fact)
1604  !$acc routine seq
1605  implicit none
1606  real(rp), intent(in), value :: dz
1607  real(rp), intent(in), value :: dx
1608  real(rp), intent(in), value :: dy
1609  real(rp), intent(in), value :: filter_fact
1610  real(rp) :: mixlen ! (out)
1611 
1612  mixlen = fact(dz, dx, dy) * filter_fact * ( dz * dx * dy )**oneoverthree ! Scotti et al. (1993)
1613 
1614  return
1615  end function mixlen
1616 
1617  function fact(dz, dx, dy)
1618  !$acc routine seq
1619  real(rp), intent(in) :: dz
1620  real(rp), intent(in) :: dx
1621  real(rp), intent(in) :: dy
1622  real(rp) :: fact ! (out)
1623 
1624  real(rp), parameter :: oot = -1.0_rp/3.0_rp
1625  real(rp), parameter :: fot = 5.0_rp/3.0_rp
1626  real(rp), parameter :: eot = 11.0_rp/3.0_rp
1627  real(rp), parameter :: tof = -3.0_rp/4.0_rp
1628  real(rp) :: a1, a2, b1, b2, dmax
1629 
1630  dmax = max(dz, dx, dy)
1631  if ( dz == dmax ) then
1632  a1 = dx / dmax
1633  a2 = dy / dmax
1634  else if ( dx == dmax ) then
1635  a1 = dz / dmax
1636  a2 = dy / dmax
1637  else ! dy == dmax
1638  a1 = dz / dmax
1639  a2 = dx / dmax
1640  end if
1641  b1 = atan( a1/a2 )
1642  b2 = atan( a2/a1 )
1643 
1644  fact = 1.736_rp * (a1*a2)**oot &
1645  * ( 4.0_rp*p1(b1)*a1**oot + 0.222_rp*p2(b1)*a1**fot + 0.077*p3(b1)*a1**eot - 3.0_rp*b1 &
1646  + 4.0_rp*p1(b2)*a2**oot + 0.222_rp*p2(b2)*a2**fot + 0.077*p3(b2)*a2**eot - 3.0_rp*b2 &
1647  )**tof
1648  return
1649  end function fact
1650  function p1(z)
1651  !$acc routine seq
1652  real(rp), intent(in) :: z
1653  real(rp) :: p1 ! (out)
1654 
1655  p1 = 2.5_rp * p2(z) - 1.5_rp * sin(z) * cos(z)**twooverthree
1656  return
1657  end function p1
1658  function p2(z)
1659  !$acc routine seq
1660  real(rp), intent(in) :: z
1661  real(rp) :: p2 ! (out)
1662 
1663  p2 = 0.986_rp * z + 0.073_rp * z**2 - 0.418_rp * z**3 + 0.120_rp * z**4
1664  return
1665  end function p2
1666  function p3(z)
1667  !$acc routine seq
1668  real(rp), intent(in) :: z
1669  real(rp) :: p3 ! (out)
1670 
1671  p3 = 0.976_rp * z + 0.188_rp * z**2 - 1.169_rp * z**3 + 0.755_rp * z**4 - 0.151_rp * z**5
1672  return
1673  end function p3
1674 
1675 end module scale_atmos_phy_tb_smg
scale_const::const_grav
real(rp), public const_grav
standard acceleration of gravity [m/s2]
Definition: scale_const.F90:49
scale_atmos_grid_cartesc_index::i_uy
integer, public i_uy
Definition: scale_atmos_grid_cartesC_index.F90:100
scale_atmos_grid_cartesc_index::ke
integer, public ke
end point of inner domain: z, local
Definition: scale_atmos_grid_cartesC_index.F90:52
scale_atmos_grid_cartesc_index::i_xv
integer, public i_xv
Definition: scale_atmos_grid_cartesC_index.F90:101
scale_prc::prc_abort
subroutine, public prc_abort
Abort Process.
Definition: scale_prc.F90:350
scale_tracer::qa
integer, public qa
Definition: scale_tracer.F90:35
scale_const::const_undef2
integer, parameter, public const_undef2
undefined value (INT2)
Definition: scale_const.F90:40
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:66
scale_precision
module PRECISION
Definition: scale_precision.F90:14
scale_atmos_grid_cartesc_index::ka
integer, public ka
Definition: scale_atmos_grid_cartesC_index.F90:47
scale_atmos_grid_cartesc_index::jblock
integer, public jblock
block size for cache blocking: y
Definition: scale_atmos_grid_cartesC_index.F90:41
scale_tracer::tracer_advc
logical, dimension(qa_max), public tracer_advc
Definition: scale_tracer.F90:46
scale_const::const_eps
real(rp), public const_eps
small number
Definition: scale_const.F90:35
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:1679
scale_atmos_grid_cartesc_index::iblock
integer, public iblock
block size for cache blocking: x
Definition: scale_atmos_grid_cartesC_index.F90:40
scale_atmos_hydrometeor
module atmosphere / hydrometeor
Definition: scale_atmos_hydrometeor.F90:12
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:1741
scale_atmos_phy_tb_smg::atmos_phy_tb_smg_finalize
subroutine, public atmos_phy_tb_smg_finalize
finalize
Definition: scale_atmos_phy_tb_smg.F90:252
scale_random
module RANDOM
Definition: scale_random.F90:11
scale_file_history
module file_history
Definition: scale_file_history.F90:15
scale_atmos_grid_cartesc_index::i_xyz
integer, public i_xyz
Definition: scale_atmos_grid_cartesC_index.F90:91
scale_prc
module PROCESS
Definition: scale_prc.F90:11
scale_atmos_phy_tb_smg::mixlen
real(rp) function mixlen(dz, dx, dy, filter_fact)
Definition: scale_atmos_phy_tb_smg.F90:1604
scale_atmos_grid_cartesc_index::i_xy
integer, public i_xy
Definition: scale_atmos_grid_cartesC_index.F90:99
scale_atmos_grid_cartesc_index::ie
integer, public ie
end point of inner domain: x, local
Definition: scale_atmos_grid_cartesC_index.F90:54
scale_atmos_grid_cartesc_index::i_uyz
integer, public i_uyz
Definition: scale_atmos_grid_cartesC_index.F90:95
scale_io
module STDIO
Definition: scale_io.F90:10
scale_atmos_grid_cartesc_index::i_xvw
integer, public i_xvw
Definition: scale_atmos_grid_cartesC_index.F90:94
scale_tracer::k
real(rp), public k
Definition: scale_tracer.F90:45
scale_atmos_grid_cartesc_index
module atmosphere / grid / cartesC index
Definition: scale_atmos_grid_cartesC_index.F90:12
scale_const
module CONSTANT
Definition: scale_const.F90:11
scale_atmos_grid_cartesc_index::ia
integer, public ia
Definition: scale_atmos_grid_cartesC_index.F90:48
scale_debug::check
subroutine, public check(current_line, v)
Undefined value checker.
Definition: scale_debug.F90:59
scale_atmos_grid_cartesc_index::zdir
integer, parameter, public zdir
Definition: scale_atmos_grid_cartesC_index.F90:32
scale_atmos_grid_cartesc_index::i_xvz
integer, public i_xvz
Definition: scale_atmos_grid_cartesC_index.F90:96
scale_prof
module profiler
Definition: scale_prof.F90:11
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, horizontal, implicit, a, b, c, dt, IIS, IIE, JJS, JJE)
Definition: scale_atmos_phy_tb_common.F90:1337
scale_atmos_grid_cartesc_index::i_uyw
integer, public i_uyw
Definition: scale_atmos_grid_cartesC_index.F90:93
scale_atmos_grid_cartesc_index::is
integer, public is
start point of inner domain: x, local
Definition: scale_atmos_grid_cartesC_index.F90:53
scale_precision::dp
integer, parameter, public dp
Definition: scale_precision.F90:32
scale_atmos_grid_cartesc_index::ja
integer, public ja
Definition: scale_atmos_grid_cartesC_index.F90:49
scale_tracer
module TRACER
Definition: scale_tracer.F90:12
scale_atmos_hydrometeor::i_qv
integer, public i_qv
Definition: scale_atmos_hydrometeor.F90:93
scale_atmos_grid_cartesc_index::xdir
integer, parameter, public xdir
Definition: scale_atmos_grid_cartesC_index.F90:33
scale_atmos_grid_cartesc_index::ks
integer, public ks
start point of inner domain: z, local
Definition: scale_atmos_grid_cartesC_index.F90:51
scale_debug
module DEBUG
Definition: scale_debug.F90:11
scale_atmos_phy_tb_common
module ATMOSPHERE / Physics Turbulence
Definition: scale_atmos_phy_tb_common.F90:12
scale_atmos_phy_tb_smg::atmos_phy_tb_smg_setup
subroutine, public atmos_phy_tb_smg_setup(FZ, CZ, CDX, CDY, MAPF, horizontal)
Setup.
Definition: scale_atmos_phy_tb_smg.F90:102
scale_atmos_grid_cartesc_index::js
integer, public js
start point of inner domain: y, local
Definition: scale_atmos_grid_cartesC_index.F90:55
scale_atmos_grid_cartesc_index::ydir
integer, parameter, public ydir
Definition: scale_atmos_grid_cartesC_index.F90:34
scale_const::const_karman
real(rp), parameter, public const_karman
von Karman constant
Definition: scale_const.F90:54
scale_matrix
module MATRIX
Definition: scale_matrix.F90:17
scale_const::const_undef
real(rp), public const_undef
Definition: scale_const.F90:43
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:1616
scale_atmos_phy_tb_smg
module ATMOSPHERE / Physics Turbulence
Definition: scale_atmos_phy_tb_smg.F90:21
scale_io::io_fid_conf
integer, public io_fid_conf
Config file ID.
Definition: scale_io.F90:57
scale_atmos_grid_cartesc_index::je
integer, public je
end point of inner domain: y, local
Definition: scale_atmos_grid_cartesC_index.F90:56
scale_atmos_phy_tb_smg::atmos_phy_tb_smg
subroutine, public atmos_phy_tb_smg(qflx_sgs_momz, qflx_sgs_momx, qflx_sgs_momy, qflx_sgs_rhot, qflx_sgs_rhoq, MOMZ_t, MOMX_t, MOMY_t, RHOT_t, RHOQ_t, Nu, Ri, Pr, MOMZ, MOMX, MOMY, POTT, DENS, QTRC, N2, FZ, FDZ, RCDZ, RFDZ, CDX, FDX, CDY, FDY, GSQRT, J13G, J23G, J33G, MAPF, dt)
Definition: scale_atmos_phy_tb_smg.F90:269
scale_atmos_grid_cartesc_index::i_xyw
integer, public i_xyw
Definition: scale_atmos_grid_cartesC_index.F90:92
Generated by   doxygen 1.8.17