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