SCALE-RM
Functions/Subroutines
scale_atmos_phy_tb_smg Module Reference

module ATMOSPHERE / Physics Turbulence More...

Functions/Subroutines

subroutine, public atmos_phy_tb_smg_setup (FZ, CZ, CDX, CDY, MAPF, horizontal)
 Setup. More...
 
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)
 
real(rp) function mixlen (dz, dx, dy, filter_fact)
 

Detailed Description

module ATMOSPHERE / Physics Turbulence

Description
Sub-grid scale turbulence process Smagolinsky-type
Author
Team SCALE
Reference
  • Brown et al., 1994: Large-eddy simulaition of stable atmospheric boundary layers with a revised stochastic subgrid model. Roy. Meteor. Soc., 120, 1485-1512
  • Scotti et al., 1993: Generalized Smagorinsky model for anisotropic grids. Phys. Fluids A, 5, 2306-2308
NAMELIST
  • PARAM_ATMOS_PHY_TB_SMG
    nametypedefault valuecomment
    ATMOS_PHY_TB_SMG_CS real(RP)
    ATMOS_PHY_TB_SMG_NU_MAX real(RP) 10000.0_RP
    ATMOS_PHY_TB_SMG_FILTER_FACT real(RP) 2.0_RP
    ATMOS_PHY_TB_SMG_BOTTOM logical .true.
    ATMOS_PHY_TB_SMG_BACKSCATTER logical .false.
    ATMOS_PHY_TB_SMG_IMPLICIT logical .false.
    ATMOS_PHY_TB_SMG_CONSISTENT_TKE logical .true.
    ATMOS_PHY_TB_SMG_HORIZONTAL logical .false.

History Output
namedescriptionunitvariable
TKE_SMG turbulent kinetic energy (Smagorinsky) m2/s2 TKE

Function/Subroutine Documentation

◆ atmos_phy_tb_smg_setup()

subroutine, public scale_atmos_phy_tb_smg::atmos_phy_tb_smg_setup ( real(rp), dimension (0:ka,ia,ja), intent(in)  FZ,
real(rp), dimension ( ka,ia,ja), intent(in)  CZ,
real(rp), dimension (ia), intent(in)  CDX,
real(rp), dimension (ja), intent(in)  CDY,
real(rp), dimension(ia,ja,2), intent(in)  MAPF,
logical, intent(in), optional  horizontal 
)

Setup.

Definition at line 100 of file scale_atmos_phy_tb_smg.F90.

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

References scale_const::const_karman, scale_atmos_grid_cartesc_index::ia, scale_atmos_grid_cartesc_index::ie, scale_io::io_fid_conf, scale_atmos_grid_cartesc_index::is, scale_atmos_grid_cartesc_index::ja, scale_atmos_grid_cartesc_index::je, scale_atmos_grid_cartesc_index::js, scale_tracer::k, scale_atmos_grid_cartesc_index::ka, scale_atmos_grid_cartesc_index::ke, scale_atmos_grid_cartesc_index::ks, mixlen(), and scale_prc::prc_abort().

Referenced by mod_atmos_phy_tb_driver::atmos_phy_tb_driver_setup().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ atmos_phy_tb_smg()

subroutine, public scale_atmos_phy_tb_smg::atmos_phy_tb_smg ( real(rp), dimension(ka,ia,ja,3), intent(out)  qflx_sgs_momz,
real(rp), dimension(ka,ia,ja,3), intent(out)  qflx_sgs_momx,
real(rp), dimension(ka,ia,ja,3), intent(out)  qflx_sgs_momy,
real(rp), dimension(ka,ia,ja,3), intent(out)  qflx_sgs_rhot,
real(rp), dimension(ka,ia,ja,3,qa), intent(out)  qflx_sgs_rhoq,
real(rp), dimension (ka,ia,ja), intent(out)  MOMZ_t,
real(rp), dimension (ka,ia,ja), intent(out)  MOMX_t,
real(rp), dimension (ka,ia,ja), intent(out)  MOMY_t,
real(rp), dimension (ka,ia,ja), intent(out)  RHOT_t,
real(rp), dimension (ka,ia,ja,qa), intent(out)  RHOQ_t,
real(rp), dimension (ka,ia,ja), intent(out)  Nu,
real(rp), dimension (ka,ia,ja), intent(out)  Ri,
real(rp), dimension (ka,ia,ja), intent(out)  Pr,
real(rp), dimension (ka,ia,ja), intent(in)  MOMZ,
real(rp), dimension (ka,ia,ja), intent(in)  MOMX,
real(rp), dimension (ka,ia,ja), intent(in)  MOMY,
real(rp), dimension (ka,ia,ja), intent(in)  POTT,
real(rp), dimension (ka,ia,ja), intent(in)  DENS,
real(rp), dimension (ka,ia,ja,qa), intent(in)  QTRC,
real(rp), dimension (ka,ia,ja), intent(in)  N2,
real(rp), dimension (0:ka,ia,ja), intent(in)  FZ,
real(rp), dimension (ka-1), intent(in)  FDZ,
real(rp), dimension (ka), intent(in)  RCDZ,
real(rp), dimension (ka-1), intent(in)  RFDZ,
real(rp), dimension (ia), intent(in)  CDX,
real(rp), dimension (ia-1), intent(in)  FDX,
real(rp), dimension (ja), intent(in)  CDY,
real(rp), dimension (ja-1), intent(in)  FDY,
real(rp), dimension (ka,ia,ja,7), intent(in)  GSQRT,
real(rp), dimension (ka,ia,ja,7), intent(in)  J13G,
real(rp), dimension (ka,ia,ja,7), intent(in)  J23G,
real(rp), intent(in)  J33G,
real(rp), dimension(ia,ja,2,4), intent(in)  MAPF,
real(dp), intent(in)  dt 
)
Parameters
[in]gsqrtvertical metrics {G}^1/2
[in]j13g(1,3) element of Jacobian matrix
[in]j23g(1,3) element of Jacobian matrix
[in]j33g(3,3) element of Jacobian matrix
[in]mapfmap factor

Definition at line 242 of file scale_atmos_phy_tb_smg.F90.

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

References scale_atmos_phy_tb_common::atmos_phy_tb_calc_flux_phi(), scale_atmos_phy_tb_common::atmos_phy_tb_calc_strain_tensor(), scale_atmos_phy_tb_common::atmos_phy_tb_calc_tend_momx(), scale_atmos_phy_tb_common::atmos_phy_tb_calc_tend_momy(), scale_atmos_phy_tb_common::atmos_phy_tb_calc_tend_momz(), scale_atmos_phy_tb_common::atmos_phy_tb_diffusion_solver(), scale_debug::check(), scale_const::const_eps, scale_const::const_grav, scale_atmos_hydrometeor::i_qv, scale_atmos_grid_cartesc_index::i_uy, scale_atmos_grid_cartesc_index::i_uyw, scale_atmos_grid_cartesc_index::i_uyz, scale_atmos_grid_cartesc_index::i_xv, scale_atmos_grid_cartesc_index::i_xvw, scale_atmos_grid_cartesc_index::i_xvz, scale_atmos_grid_cartesc_index::i_xy, scale_atmos_grid_cartesc_index::i_xyw, scale_atmos_grid_cartesc_index::i_xyz, scale_atmos_grid_cartesc_index::ia, scale_atmos_grid_cartesc_index::iblock, scale_atmos_grid_cartesc_index::ie, scale_atmos_grid_cartesc_index::is, scale_atmos_grid_cartesc_index::ja, scale_atmos_grid_cartesc_index::jblock, scale_atmos_grid_cartesc_index::je, scale_atmos_grid_cartesc_index::js, scale_tracer::k, scale_atmos_grid_cartesc_index::ke, scale_atmos_grid_cartesc_index::ks, scale_tracer::qa, scale_tracer::tracer_advc, scale_atmos_grid_cartesc_index::xdir, scale_atmos_grid_cartesc_index::ydir, and scale_atmos_grid_cartesc_index::zdir.

Referenced by mod_atmos_phy_tb_driver::atmos_phy_tb_driver_calc_tendency().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ mixlen()

real(rp) function scale_atmos_phy_tb_smg::mixlen ( real(rp), intent(in)  dz,
real(rp), intent(in)  dx,
real(rp), intent(in)  dy,
real(rp), intent(in)  filter_fact 
)

Definition at line 1396 of file scale_atmos_phy_tb_smg.F90.

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

Referenced by atmos_phy_tb_smg_setup().

Here is the caller graph for this function:
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::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_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::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_grid_cartesc_index::i_xy
integer, public i_xy
Definition: scale_atmos_grid_cartesC_index.F90:98
scale_atmos_grid_cartesc_index::i_uyz
integer, public i_uyz
Definition: scale_atmos_grid_cartesC_index.F90:94
scale_atmos_grid_cartesc_index::i_xvw
integer, public i_xvw
Definition: scale_atmos_grid_cartesC_index.F90:93
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::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_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_tracer
module TRACER
Definition: scale_tracer.F90:12
scale_atmos_grid_cartesc_index::xdir
integer, parameter, public xdir
Definition: scale_atmos_grid_cartesC_index.F90:33
scale_atmos_phy_tb_common
module ATMOSPHERE / Physics Turbulence
Definition: scale_atmos_phy_tb_common.F90:12
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_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_grid_cartesc_index::i_xyw
integer, public i_xyw
Definition: scale_atmos_grid_cartesC_index.F90:91
Generated by   doxygen 1.8.17