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)
 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 
)

Setup.

Definition at line 99 of file scale_atmos_phy_tb_smg.F90.

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_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().

99  use scale_prc, only: &
100  prc_abort
101  use scale_const, only: &
102  karman => const_karman
103  implicit none
104 
105  real(RP), intent(in) :: fz (0:ka,ia,ja)
106  real(RP), intent(in) :: cz ( ka,ia,ja)
107  real(RP), intent(in) :: cdx (ia)
108  real(RP), intent(in) :: cdy (ja)
109  real(RP), intent(in) :: mapf(ia,ja,2)
110 
111  real(RP) :: atmos_phy_tb_smg_cs
112  real(RP) :: atmos_phy_tb_smg_filter_fact = 2.0_rp
113  logical :: atmos_phy_tb_smg_consistent_tke = .true.
114 
115  namelist / param_atmos_phy_tb_smg / &
116  atmos_phy_tb_smg_cs, &
117  atmos_phy_tb_smg_nu_max, &
118  atmos_phy_tb_smg_filter_fact, &
119  atmos_phy_tb_smg_bottom, &
120  atmos_phy_tb_smg_backscatter, &
121  atmos_phy_tb_smg_implicit, &
122  atmos_phy_tb_smg_consistent_tke, &
123  atmos_phy_tb_smg_horizontal
124 
125  integer :: ierr
126  integer :: k, i, j
127  !---------------------------------------------------------------------------
128 
129  log_newline
130  log_info("ATMOS_PHY_TB_smg_setup",*) 'Setup'
131  log_info("ATMOS_PHY_TB_smg_setup",*) 'Smagorinsky-type Eddy Viscocity Model'
132 
133  atmos_phy_tb_smg_cs = cs
134 
135  !--- read namelist
136  rewind(io_fid_conf)
137  read(io_fid_conf,nml=param_atmos_phy_tb_smg,iostat=ierr)
138  if( ierr < 0 ) then !--- missing
139  log_info("ATMOS_PHY_TB_smg_setup",*) 'Not found namelist. Default used.'
140  elseif( ierr > 0 ) then !--- fatal error
141  log_error("ATMOS_PHY_TB_smg_setup",*) 'Not appropriate names in namelist PARAM_ATMOS_PHY_TB_SMG. Check!'
142  call prc_abort
143  endif
144  log_nml(param_atmos_phy_tb_smg)
145 
146  cs = atmos_phy_tb_smg_cs
147 
148  rprn = 1.0_rp / prn
149  rric = 1.0_rp / ric
150  prnovric = ( 1.0_rp - prn ) * rric
151 
152  allocate( lambda0(ka,ia,ja) )
153  allocate( lambda(ka,ia,ja) )
154 
155 #ifdef DEBUG
156  lambda0(:,:,:) = undef
157  lambda(:,:,:) = undef
158 #endif
159  if ( atmos_phy_tb_smg_horizontal ) then
160  !$omp parallel do
161  do j = js-1, je+1
162  do i = is-1, ie+1
163  do k = ks, ke
164  lambda0(k,i,j) = cs * sqrt( cdx(i) * cdy(j) / ( mapf(i,j,1) * mapf(i,j,2) ) )
165  lambda(k,i,j) = lambda0(k,i,j)
166  enddo
167  enddo
168  enddo
169 #ifdef DEBUG
170  i = iundef; j = iundef; k = iundef
171 #endif
172  atmos_phy_tb_smg_consistent_tke = .false.
173  atmos_phy_tb_smg_implicit = .false. ! flux in the z-direction is not necessary
174  atmos_phy_tb_smg_backscatter = .false.
175  else
176  !$omp parallel do
177  do j = js-1, je+1
178  do i = is-1, ie+1
179  do k = ks, ke
180  lambda0(k,i,j) = cs * mixlen( fz(k,i,j) - fz(k-1,i,j), &
181  cdx(i) / mapf(i,j,1), &
182  cdy(j) / mapf(i,j,2), &
183  atmos_phy_tb_smg_filter_fact )
184  enddo
185  enddo
186  enddo
187 #ifdef DEBUG
188  i = iundef; j = iundef; k = iundef
189 #endif
190 
191  if ( atmos_phy_tb_smg_bottom ) then
192  !$omp parallel do
193  do j = js-1, je+1
194  do i = is-1, ie+1
195  do k = ks, ke
196  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 ) )
197  enddo
198  enddo
199  enddo
200 #ifdef DEBUG
201  i = iundef; j = iundef; k = iundef
202 #endif
203  else
204  !$omp parallel do
205  do j = js-1, je+1
206  do i = is-1, ie+1
207  do k = ks, ke
208  lambda(k,i,j) = lambda0(k,i,j)
209  enddo
210  enddo
211  enddo
212 #ifdef DEBUG
213  i = iundef; j = iundef; k = iundef
214 #endif
215  end if
216  end if
217 
218  ! flag for isotropic stress tensor
219  if ( atmos_phy_tb_smg_consistent_tke ) then
220  tke_fact = 1.0_rp
221  else
222  tke_fact = 0.0_rp ! neglect
223  end if
224 
225  return
scale_atmos_grid_cartesc_index::ia
integer, public ia
of whole cells: x, local, with HALO
Definition: scale_atmos_grid_cartesC_index.F90:48
scale_atmos_grid_cartesc_index::ja
integer, public ja
of whole cells: y, local, with HALO
Definition: scale_atmos_grid_cartesC_index.F90:49
scale_io::io_fid_conf
integer, public io_fid_conf
Config file ID.
Definition: scale_io.F90:55
scale_const::const_karman
real(rp), parameter, public const_karman
von Karman constant
Definition: scale_const.F90:50
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_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::ke
integer, public ke
end point of inner domain: z, local
Definition: scale_atmos_grid_cartesC_index.F90:52
scale_prc
module PROCESS
Definition: scale_prc.F90:11
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_grid_cartesc_index::ks
integer, public ks
start point of inner domain: z, local
Definition: scale_atmos_grid_cartesC_index.F90:51
scale_prc::prc_abort
subroutine, public prc_abort
Abort Process.
Definition: scale_prc.F90:338
scale_const
module CONSTANT
Definition: scale_const.F90:11
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::ka
integer, public ka
of whole cells: z, local, with HALO
Definition: scale_atmos_grid_cartesC_index.F90:47
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 237 of file scale_atmos_phy_tb_smg.F90.

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_uyw, scale_atmos_grid_cartesc_index::i_uyz, 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_atmos_grid_cartesc_index::ke, scale_atmos_grid_cartesc_index::ks, scale_tracer::qa, scale_random::random_normal(), 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().

237  use scale_precision
238  use scale_atmos_grid_cartesc_index
239  use scale_tracer
240  use scale_const, only: &
241  eps => const_eps, &
242  grav => const_grav
243  use scale_file_history, only: &
244  file_history_in
245  use scale_atmos_phy_tb_common, only: &
246  calc_strain_tensor => atmos_phy_tb_calc_strain_tensor, &
247  diffusion_solver => atmos_phy_tb_diffusion_solver, &
248  calc_tend_momz => atmos_phy_tb_calc_tend_momz, &
249  calc_tend_momx => atmos_phy_tb_calc_tend_momx, &
250  calc_tend_momy => atmos_phy_tb_calc_tend_momy, &
251  calc_flux_phi => atmos_phy_tb_calc_flux_phi
252  use scale_atmos_hydrometeor, only: &
253  i_qv
254  use scale_random, only: &
255  random_normal
256  implicit none
257 
258  ! SGS flux
259  real(RP), intent(out) :: qflx_sgs_momz(ka,ia,ja,3)
260  real(RP), intent(out) :: qflx_sgs_momx(ka,ia,ja,3)
261  real(RP), intent(out) :: qflx_sgs_momy(ka,ia,ja,3)
262  real(RP), intent(out) :: qflx_sgs_rhot(ka,ia,ja,3)
263  real(RP), intent(out) :: qflx_sgs_rhoq(ka,ia,ja,3,qa)
264 
265  real(RP), intent(out) :: momz_t (ka,ia,ja)
266  real(RP), intent(out) :: momx_t (ka,ia,ja)
267  real(RP), intent(out) :: momy_t (ka,ia,ja)
268  real(RP), intent(out) :: rhot_t (ka,ia,ja)
269  real(RP), intent(out) :: rhoq_t (ka,ia,ja,qa) ! tendency of rho * QTRC
270 
271  real(RP), intent(out) :: nu (ka,ia,ja) ! eddy viscosity (center)
272  real(RP), intent(out) :: ri (ka,ia,ja) ! Richardson number
273  real(RP), intent(out) :: pr (ka,ia,ja) ! Prantle number
274 
275  real(RP), intent(in) :: momz (ka,ia,ja)
276  real(RP), intent(in) :: momx (ka,ia,ja)
277  real(RP), intent(in) :: momy (ka,ia,ja)
278  real(RP), intent(in) :: pott (ka,ia,ja)
279  real(RP), intent(in) :: dens (ka,ia,ja)
280  real(RP), intent(in) :: qtrc (ka,ia,ja,qa)
281  real(RP), intent(in) :: n2 (ka,ia,ja)
282 
283  real(RP), intent(in) :: fz (0:ka,ia,ja)
284  real(RP), intent(in) :: fdz (ka-1)
285  real(RP), intent(in) :: rcdz (ka)
286  real(RP), intent(in) :: rfdz (ka-1)
287  real(RP), intent(in) :: cdx (ia)
288  real(RP), intent(in) :: fdx (ia-1)
289  real(RP), intent(in) :: cdy (ja)
290  real(RP), intent(in) :: fdy (ja-1)
291 
292  real(RP), intent(in) :: gsqrt (ka,ia,ja,7)
293  real(RP), intent(in) :: j13g (ka,ia,ja,7)
294  real(RP), intent(in) :: j23g (ka,ia,ja,7)
295  real(RP), intent(in) :: j33g
296  real(RP), intent(in) :: mapf(ia,ja,2,4)
297  real(DP), intent(in) :: dt
298 
299  ! tke
300  real(RP) :: tke(ka,ia,ja)
301 
302  ! deformation rate tensor
303  real(RP) :: s33_c(ka,ia,ja) ! (cell center)
304  real(RP) :: s11_c(ka,ia,ja)
305  real(RP) :: s22_c(ka,ia,ja)
306  real(RP) :: s31_c(ka,ia,ja)
307  real(RP) :: s12_c(ka,ia,ja)
308  real(RP) :: s23_c(ka,ia,ja)
309  real(RP) :: s12_z(ka,ia,ja) ! (z edge or x-y plane)
310  real(RP) :: s23_x(ka,ia,ja) ! (x edge or y-z plane)
311  real(RP) :: s31_y(ka,ia,ja) ! (y edge or z-x plane)
312  real(RP) :: s2 (ka,ia,ja) ! |S|^2
313 
314  real(RP) :: kh(ka,ia,ja) ! eddy diffusion
315  real(RP) :: fm(ka)
316  real(RP) :: e(ka)
317  real(RP) :: et
318  real(RP) :: lambda_r(ka)
319  real(RP) :: rf
320  real(RP) :: c1(ka)
321  real(RP) :: c2
322  real(RP) :: d2
323 
324  ! implicit scheme
325  real(RP) :: tend(ka,ia,ja)
326  real(RP) :: a (ka,ia,ja)
327  real(RP) :: b (ka,ia,ja)
328  real(RP) :: c (ka,ia,ja)
329  real(RP) :: d (ka)
330  real(RP) :: ap
331 
332  ! backscatter
333  real(RP) :: random (ka,ia,ja)
334  real(RP) :: random_mz(ka,ia,ja)
335  real(RP) :: random_mx(ka,ia,ja)
336  real(RP) :: random_my(ka,ia,ja)
337  real(RP) :: random_qz(ka,ia,ja)
338  real(RP) :: random_qx(ka,ia,ja)
339  real(RP) :: random_qy(ka,ia,ja)
340  real(RP) :: acc (ka,ia,ja)
341  real(RP) :: dd (ka,ia,ja)
342  real(RP) :: leovleo5
343  real(RP) :: dz, dx, dy
344  real(RP) :: fact
345 
346  integer :: iis, iie
347  integer :: jjs, jje
348 
349  integer :: k, i, j, iq
350  !---------------------------------------------------------------------------
351 
352  log_progress(*) 'atmosphere / physics / turbulence / Smagorinsky'
353 
354 #ifdef DEBUG
355  qflx_sgs_momz(:,:,:,:) = undef
356  qflx_sgs_momx(:,:,:,:) = undef
357  qflx_sgs_momy(:,:,:,:) = undef
358  qflx_sgs_rhot(:,:,:,:) = undef
359  qflx_sgs_rhoq(:,:,:,:,:) = undef
360 
361  nu(:,:,:) = undef
362  tke(:,:,:) = undef
363  pr(:,:,:) = undef
364  ri(:,:,:) = undef
365  kh(:,:,:) = undef
366 #endif
367 
368 #ifdef QUICKDEBUG
369  qflx_sgs_momz(ks:ke, 1:is-1, : ,:) = undef
370  qflx_sgs_momz(ks:ke,ie+1:ia , : ,:) = undef
371  qflx_sgs_momz(ks:ke, : , 1:js-1,:) = undef
372  qflx_sgs_momz(ks:ke, : ,je+1:ja ,:) = undef
373  qflx_sgs_momx(ks:ke, 1:is-1, : ,:) = undef
374  qflx_sgs_momx(ks:ke,ie+1:ia , : ,:) = undef
375  qflx_sgs_momx(ks:ke, : , 1:js-1,:) = undef
376  qflx_sgs_momx(ks:ke, : ,je+1:ja ,:) = undef
377  qflx_sgs_momy(ks:ke, 1:is-1, : ,:) = undef
378  qflx_sgs_momy(ks:ke,ie+1:ia , : ,:) = undef
379  qflx_sgs_momy(ks:ke, : , 1:js-1,:) = undef
380  qflx_sgs_momy(ks:ke, : ,je+1:ja ,:) = undef
381 #endif
382 
383 
384  !##### Start Upadate #####
385 
386  call calc_strain_tensor( &
387  s33_c, s11_c, s22_c, & ! (out)
388  s31_c, s12_c, s23_c, & ! (out)
389  s12_z, s23_x, s31_y, & ! (out)
390  s2 , & ! (out)
391  dens, momz, momx, momy, & ! (in)
392  gsqrt, j13g, j23g, j33g, mapf ) ! (in)
393 
394  if ( atmos_phy_tb_smg_backscatter ) then
395 
396  call random_normal( random(:,:,:) )
397  ! 1:2:1 filter
398  !$omp parallel do
399  do j = js-1, je+1
400  do i = is-1, ie+1
401  do k = ks+1, ke-1
402  random_mz(k,i,j) = ( random(k,i,j) * 2.0_rp &
403  + random(k+1,i,j) + random(k-1,i,j) &
404  + random(k,i+1,j) + random(k,i-1,j) &
405  + random(k,i,j+1) + random(k,i,j-1) ) / 8.0_rp
406  end do
407  random_mz(ks,i,j) = ( random(ks,i,j) * 2.0_rp &
408  + random(ks+1,i,j) &
409  + random(ks,i+1,j) + random(ks,i-1,j) &
410  + random(ks,i,j+1) + random(ks,i,j-1) ) / 7.0_rp
411  random_mz(ke,i,j) = ( random(ke,i,j) * 2.0_rp &
412  + random(ke-1,i,j) &
413  + random(ke,i+1,j) + random(ke,i-1,j) &
414  + random(ke,i,j+1) + random(ke,i,j-1) ) / 7.0_rp
415  end do
416  end do
417 
418  call random_normal( random(:,:,:) )
419  ! 1:2:1 filter
420  !$omp parallel do
421  do j = js-1, je+1
422  do i = is-1, ie+1
423  do k = ks+1, ke-1
424  random_mx(k,i,j) = ( random(k,i,j) * 2.0_rp &
425  + random(k+1,i,j) + random(k-1,i,j) &
426  + random(k,i+1,j) + random(k,i-1,j) &
427  + random(k,i,j+1) + random(k,i,j-1) ) / 8.0_rp
428  end do
429  random_mx(ks,i,j) = ( random(ks,i,j) * 2.0_rp &
430  + random(ks+1,i,j) &
431  + random(ks,i+1,j) + random(ks,i-1,j) &
432  + random(ks,i,j+1) + random(ks,i,j-1) ) / 7.0_rp
433  random_mx(ke,i,j) = ( random(ke,i,j) * 2.0_rp &
434  + random(ke-1,i,j) &
435  + random(ke,i+1,j) + random(ke,i-1,j) &
436  + random(ke,i,j+1) + random(ke,i,j-1) ) / 7.0_rp
437  end do
438  end do
439 
440  call random_normal( random(:,:,:) )
441  ! 1:2:1 filter
442  !$omp parallel do
443  do j = js-1, je+1
444  do i = is-1, ie+1
445  do k = ks+1, ke-1
446  random_my(k,i,j) = ( random(k,i,j) * 2.0_rp &
447  + random(k+1,i,j) + random(k-1,i,j) &
448  + random(k,i+1,j) + random(k,i-1,j) &
449  + random(k,i,j+1) + random(k,i,j-1) ) / 8.0_rp
450  end do
451  random_my(ks,i,j) = ( random(ks,i,j) * 2.0_rp &
452  + random(ks+1,i,j) &
453  + random(ks,i+1,j) + random(ks,i-1,j) &
454  + random(ks,i,j+1) + random(ks,i,j-1) ) / 7.0_rp
455  random_my(ke,i,j) = ( random(ke,i,j) * 2.0_rp &
456  + random(ke-1,i,j) &
457  + random(ke,i+1,j) + random(ke,i-1,j) &
458  + random(ke,i,j+1) + random(ke,i,j-1) ) / 7.0_rp
459  end do
460  end do
461 
462  call random_normal( random(:,:,:) )
463  ! 1:2:1 filter
464  !$omp parallel do
465  do j = js-1, je+1
466  do i = is-1, ie+1
467  do k = ks+1, ke-1
468  random_qz(k,i,j) = ( random(k,i,j) * 2.0_rp &
469  + random(k+1,i,j) + random(k-1,i,j) &
470  + random(k,i+1,j) + random(k,i-1,j) &
471  + random(k,i,j+1) + random(k,i,j-1) ) / 8.0_rp
472  end do
473  random_qz(ks,i,j) = ( random(ks,i,j) * 2.0_rp &
474  + random(ks+1,i,j) &
475  + random(ks,i+1,j) + random(ks,i-1,j) &
476  + random(ks,i,j+1) + random(ks,i,j-1) ) / 7.0_rp
477  random_qz(ke,i,j) = ( random(ke,i,j) * 2.0_rp &
478  + random(ke-1,i,j) &
479  + random(ke,i+1,j) + random(ke,i-1,j) &
480  + random(ke,i,j+1) + random(ke,i,j-1) ) / 7.0_rp
481  end do
482  end do
483 
484  call random_normal( random(:,:,:) )
485  ! 1:2:1 filter
486  !$omp parallel do
487  do j = js-1, je+1
488  do i = is-1, ie+1
489  do k = ks+1, ke-1
490  random_qx(k,i,j) = ( random(k,i,j) * 2.0_rp &
491  + random(k+1,i,j) + random(k-1,i,j) &
492  + random(k,i+1,j) + random(k,i-1,j) &
493  + random(k,i,j+1) + random(k,i,j-1) ) / 8.0_rp
494  end do
495  random_qx(ks,i,j) = ( random(ks,i,j) * 2.0_rp &
496  + random(ks+1,i,j) &
497  + random(ks,i+1,j) + random(ks,i-1,j) &
498  + random(ks,i,j+1) + random(ks,i,j-1) ) / 7.0_rp
499  random_qx(ke,i,j) = ( random(ke,i,j) * 2.0_rp &
500  + random(ke-1,i,j) &
501  + random(ke,i+1,j) + random(ke,i-1,j) &
502  + random(ke,i,j+1) + random(ke,i,j-1) ) / 7.0_rp
503  end do
504  end do
505 
506  call random_normal( random(:,:,:) )
507  ! 1:2:1 filter
508  !$omp parallel do
509  do j = js-1, je+1
510  do i = is-1, ie+1
511  do k = ks+1, ke-1
512  random_qy(k,i,j) = ( random(k,i,j) * 2.0_rp &
513  + random(k+1,i,j) + random(k-1,i,j) &
514  + random(k,i+1,j) + random(k,i-1,j) &
515  + random(k,i,j+1) + random(k,i,j-1) ) / 8.0_rp
516  end do
517  random_qy(ks,i,j) = ( random(ks,i,j) * 2.0_rp &
518  + random(ks+1,i,j) &
519  + random(ks,i+1,j) + random(ks,i-1,j) &
520  + random(ks,i,j+1) + random(ks,i,j-1) ) / 7.0_rp
521  random_qy(ke,i,j) = ( random(ke,i,j) * 2.0_rp &
522  + random(ke-1,i,j) &
523  + random(ke,i+1,j) + random(ke,i-1,j) &
524  + random(ke,i,j+1) + random(ke,i,j-1) ) / 7.0_rp
525  end do
526  end do
527 
528 
529  end if
530 
531 
532  !$omp parallel do &
533  !$omp private(fm,Rf,lambda_r,leOvleo5,C1,C2,D2,e,et,fact,dz,dx,dy)
534  do j = js-1, je+1
535  do i = is-1, ie+1
536 
537  ! Ri = N^2 / |S|^2, N^2 = g / theta * dtheta/dz
538  do k = ks, ke
539  ri(k,i,j) = n2(k,i,j) / max(s2(k,i,j), eps)
540  enddo
541 
542  ! Nu
543  ! Pr = Nu / Kh = fm / fh
544  do k = ks, ke
545  if ( ri(k,i,j) < 0.0_rp ) then ! stable
546  fm(k) = sqrt( 1.0_rp - fmc * ri(k,i,j) )
547  nu(k,i,j) = lambda(k,i,j)**2 * sqrt( s2(k,i,j) ) * fm(k)
548  pr(k,i,j) = fm(k) / sqrt( 1.0_rp - fhb*ri(k,i,j) ) * prn
549  else if ( ri(k,i,j) < ric ) then ! weakly stable
550  fm(k) = ( 1.0_rp - ri(k,i,j)*rric )**4
551  nu(k,i,j) = lambda(k,i,j)**2 * sqrt( s2(k,i,j) ) * fm(k)
552  pr(k,i,j) = prn / ( 1.0_rp - prnovric * ri(k,i,j) )
553  else ! strongly stable
554  fm(k) = 0.0_rp
555  nu(k,i,j) = 0.0_rp
556  kh(k,i,j) = 0.0_rp
557  pr(k,i,j) = 1.0_rp
558  endif
559 
560  if ( ri(k,i,j) < ric ) then
561  kh(k,i,j) = max( min( nu(k,i,j) / pr(k,i,j), atmos_phy_tb_smg_nu_max ), eps )
562  nu(k,i,j) = max( min( nu(k,i,j), atmos_phy_tb_smg_nu_max ), eps )
563  pr(k,i,j) = nu(k,i,j) / kh(k,i,j)
564  rf = ri(k,i,j) / pr(k,i,j)
565  lambda_r(k) = lambda(k,i,j) * sqrt( fm(k) / sqrt( 1.0_rp - rf ) )
566  else
567  lambda_r(k) = 0.0_rp
568  end if
569 
570  enddo
571 
572  if ( atmos_phy_tb_smg_backscatter ) then
573 
574  do k = ks, ke
575  lambda_r(k) = min( 1.8_rp * lambda0(k,i,j), lambda_r(k) )
576  leovleo5 = ( lambda_r(k) / lambda0(k,i,j) )**5
577  c2 = cb * leovleo5 / ( 1.0_rp + cb * leovleo5 )
578  d2 = c2 * cbt / cb
579  c1(k) = c1o / sqrt( 1.0_rp - c2 )
580  ! D1 = D1o / sqrt( 1.0_RP - C2 )
581  e(k) = nu(k,i,j)**3 * ( 1.0_rp - c2 ) / ( lambda_r(k)**4 + eps )
582  et = kh(k,i,j) * ( 1.0_rp - d2 )
583 
584  dz = fz(k,i,j) - fz(k-1,i,j)
585  dx = cdx(i) / mapf(i,j,1,i_xy)
586  dy = cdy(j) / mapf(i,j,2,i_xy)
587 
588  fact = sqrt( cb * leovleo5 * e(k) / ( dt * ( dz**2 + dx**2 + dy**2 ) ) ) * dens(k,i,j)
589  random_mz(k,i,j) = random_mz(k,i,j) * fact * dx * dy
590  random_mx(k,i,j) = random_mx(k,i,j) * fact * dy * dz
591  random_my(k,i,j) = random_my(k,i,j) * fact * dz * dx
592 
593  fact = sqrt( twooverthree * cbt * leovleo5 * et / dt ) * dens(k,i,j)
594  random_qz(k,i,j) = random_qz(k,i,j) * fact * dz
595  random_qx(k,i,j) = random_qx(k,i,j) * fact * dx
596  random_qy(k,i,j) = random_qy(k,i,j) * fact * dy
597  end do
598 
599  else
600 
601  do k = ks, ke
602  e(k) = nu(k,i,j)**3 / ( lambda_r(k)**4 + eps )
603  c1(k) = c1o
604  end do
605 
606  end if
607 
608  ! TKE
609  do k = ks, ke
610  tke(k,i,j) = ( e(k) * lambda_r(k) / c1(k) )**twooverthree
611  enddo
612 
613  enddo
614  enddo
615 #ifdef DEBUG
616  i = iundef; j = iundef; k = iundef
617 #endif
618 
619 
620 
621  do jjs = js, je, jblock
622  jje = jjs+jblock-1
623  do iis = is, ie, iblock
624  iie = iis+iblock-1
625 
626  !##### momentum equation (z) #####
627  ! (cell center)
628  if ( atmos_phy_tb_smg_horizontal ) then
629  qflx_sgs_momz(:,:,:,zdir) = 0.0_rp
630  else
631  !$omp parallel do
632  do j = jjs, jje
633  do i = iis, iie
634  do k = ks+1, ke-1
635 #ifdef DEBUG
636  call check( __line__, dens(k,i,j) )
637  call check( __line__, nu(k,i,j) )
638  call check( __line__, s33_c(k,i,j) )
639  call check( __line__, s11_c(k,i,j) )
640  call check( __line__, s22_c(k,i,j) )
641  call check( __line__, tke(k,i,j) )
642 #endif
643  qflx_sgs_momz(k,i,j,zdir) = dens(k,i,j) * ( &
644  - 2.0_rp * nu(k,i,j) &
645  * ( s33_c(k,i,j) - ( s11_c(k,i,j) + s22_c(k,i,j) + s33_c(k,i,j) ) * oneoverthree ) &
646  + twooverthree * tke(k,i,j) * tke_fact )
647  enddo
648  enddo
649  enddo
650 #ifdef DEBUG
651  i = iundef; j = iundef; k = iundef
652 #endif
653  !$omp parallel do
654  do j = jjs, jje
655  do i = iis, iie
656  ! momentum will not be conserved
657  qflx_sgs_momz(ks,i,j,zdir) = dens(ks,i,j) * twooverthree * tke(ks,i,j) * tke_fact
658  qflx_sgs_momz(ke,i,j,zdir) = dens(ke,i,j) * twooverthree * tke(ke,i,j) * tke_fact
659  ! anti-isotropic stress is calculated by the surface scheme
660  enddo
661  enddo
662 #ifdef DEBUG
663  i = iundef; j = iundef; k = iundef
664 #endif
665  end if
666  ! (y edge)
667  !$omp parallel do
668  do j = jjs, jje
669  do i = iis-1, iie
670  do k = ks, ke-1
671 #ifdef DEBUG
672  call check( __line__, dens(k,i,j) )
673  call check( __line__, dens(k,i+1,j) )
674  call check( __line__, dens(k+1,i,j) )
675  call check( __line__, dens(k+1,i+1,j) )
676  call check( __line__, nu(k,i,j) )
677  call check( __line__, nu(k,i+1,j) )
678  call check( __line__, nu(k+1,i,j) )
679  call check( __line__, nu(k+1,i+1,j) )
680  call check( __line__, s31_y(k,i,j) )
681 #endif
682  qflx_sgs_momz(k,i,j,xdir) = - 0.125_rp & ! 2.0 / 4 / 4
683  * ( dens(k,i,j)+dens(k+1,i,j)+dens(k,i+1,j)+dens(k+1,i+1,j) ) &
684  * ( nu(k,i,j)+nu(k+1,i,j)+nu(k,i+1,j)+nu(k+1,i+1,j)) &
685  * s31_y(k,i,j)
686  enddo
687  enddo
688  enddo
689 #ifdef DEBUG
690  i = iundef; j = iundef; k = iundef
691 #endif
692  ! (x edge)
693  !$omp parallel do default(none) &
694  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,DENS,S23_X,nu,qflx_sgs_momz) &
695  !$omp private(i,j,k) OMP_SCHEDULE_ collapse(2)
696  do j = jjs-1, jje
697  do i = iis, iie
698  do k = ks, ke-1
699 #ifdef DEBUG
700  call check( __line__, dens(k,i,j) )
701  call check( __line__, dens(k,i,j+1) )
702  call check( __line__, dens(k+1,i,j) )
703  call check( __line__, dens(k+1,i,j+1) )
704  call check( __line__, nu(k,i,j) )
705  call check( __line__, nu(k,i,j+1) )
706  call check( __line__, nu(k+1,i,j) )
707  call check( __line__, nu(k+1,i,j+1) )
708  call check( __line__, s23_x(k,i,j) )
709 #endif
710  qflx_sgs_momz(k,i,j,ydir) = - 0.125_rp & ! 2/4/4
711  * ( dens(k,i,j)+dens(k+1,i,j)+dens(k,i,j+1)+dens(k+1,i,j+1) ) &
712  * ( nu(k,i,j)+nu(k+1,i,j)+nu(k,i,j+1)+nu(k+1,i,j+1) ) &
713  * s23_x(k,i,j)
714  enddo
715  enddo
716  enddo
717 #ifdef DEBUG
718  i = iundef; j = iundef; k = iundef
719 #endif
720 
721  if ( atmos_phy_tb_smg_implicit ) then
722 
723  call calc_tend_momz( tend, & ! (out)
724  qflx_sgs_momz, & ! (in)
725  gsqrt, j13g, j23g, j33g, mapf, & ! (in)
726  iis, iie, jjs, jje ) ! (in)
727 
728  !$omp parallel do &
729  !$omp private (ap,d)
730  do j = jjs, jje
731  do i = iis, iie
732 
733  ap = - fouroverthree * dt &
734  * dens(ks+1,i,j)*nu(ks+1,i,j) &
735  * rcdz(ks+1) / gsqrt(ks+1,i,j,i_xyz)
736  a(ks,i,j) = ap * rfdz(ks) / gsqrt(ks,i,j,i_xyw)
737  c(ks,i,j) = 0.0_rp
738  b(ks,i,j) = - a(ks,i,j) + 0.5_rp * ( dens(ks,i,j)+dens(ks+1,i,j) )
739  do k = ks+1, ke-2
740  c(k,i,j) = ap * rfdz(k+1) / gsqrt(k+1,i,j,i_xyw)
741  ap = - fouroverthree * dt &
742  * dens(k+1,i,j)*nu(k+1,i,j) &
743  * rcdz(k+1) / gsqrt(k+1,i,j,i_xyz)
744  a(k,i,j) = ap * rfdz(k) / gsqrt(k,i,j,i_xyw)
745  b(k,i,j) = - a(k,i,j) - c(k,i,j) + 0.5_rp * ( dens(k,i,j)+dens(k+1,i,j) )
746  end do
747  a(ke-1,i,j) = 0.0_rp
748  c(ke-1,i,j) = ap * rfdz(ke) / gsqrt(ke,i,j,i_xyw)
749  b(ke-1,i,j) = - c(ke-1,i,j) + 0.5_rp * ( dens(ke-1,i,j)+dens(ke,i,j) )
750 
751  do k = ks, ke-1
752  d(k) = tend(k,i,j)
753  end do
754 
755  call diffusion_solver( &
756  tend(:,i,j), & ! (out)
757  a(:,i,j), b(:,i,j), c(:,i,j), d, & ! (in)
758  ke-1 ) ! (in)
759 
760  do k = ks+1, ke-1
761  qflx_sgs_momz(k,i,j,zdir) = qflx_sgs_momz(k,i,j,zdir) &
762  - fouroverthree * dens(k,i,j) * nu(k,i,j) * dt &
763  * ( tend(k,i,j) - tend(k-1,i,j) ) * rcdz(k) / gsqrt(k,i,j,i_xyz)
764  end do
765 
766  end do
767  end do
768 
769  end if
770 
771  !##### momentum equation (x) #####
772  ! (y edge)
773  if ( atmos_phy_tb_smg_horizontal ) then
774  qflx_sgs_momx(:,:,:,zdir) = 0.0_rp
775  else
776  !$omp parallel do
777  do j = jjs, jje
778  do i = iis, iie
779  do k = ks, ke-1
780 #ifdef DEBUG
781  call check( __line__, dens(k,i,j) )
782  call check( __line__, dens(k,i+1,j) )
783  call check( __line__, dens(k+1,i,j) )
784  call check( __line__, dens(k+1,i+1,j) )
785  call check( __line__, nu(k,i,j) )
786  call check( __line__, nu(k,i+1,j) )
787  call check( __line__, nu(k+1,i,j) )
788  call check( __line__, nu(k+1,i+1,j) )
789  call check( __line__, s31_y(k,i,j) )
790 #endif
791  qflx_sgs_momx(k,i,j,zdir) = - 0.125_rp & ! 2/4/4
792  * ( dens(k,i,j)+dens(k+1,i,j)+dens(k,i+1,j)+dens(k+1,i+1,j) ) &
793  * ( nu(k,i,j)+nu(k+1,i,j)+nu(k,i+1,j)+nu(k+1,i+1,j) ) &
794  * s31_y(k,i,j)
795  enddo
796  enddo
797  enddo
798 #ifdef DEBUG
799  i = iundef; j = iundef; k = iundef
800 #endif
801  !$omp parallel do
802  do j = jjs, jje
803  do i = iis, iie
804  qflx_sgs_momx(ks-1,i,j,zdir) = 0.0_rp ! bottom boundary
805  qflx_sgs_momx(ke ,i,j,zdir) = 0.0_rp ! top boundary
806  enddo
807  enddo
808 #ifdef DEBUG
809  i = iundef; j = iundef; k = iundef
810 #endif
811  end if
812  ! (cell center)
813  !$omp parallel do default(none) &
814  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,DENS,nu,S11_C,S22_C,S33_C,TKE,tke_fact,qflx_sgs_momx) &
815  !$omp private(i,j,k) OMP_SCHEDULE_ collapse(2)
816  do j = jjs, jje
817  do i = iis, iie+1
818  do k = ks, ke
819 #ifdef DEBUG
820  call check( __line__, dens(k,i,j) )
821  call check( __line__, nu(k,i,j) )
822  call check( __line__, s11_c(k,i,j) )
823  call check( __line__, s22_c(k,i,j) )
824  call check( __line__, s33_c(k,i,j) )
825  call check( __line__, tke(k,i,j) )
826 #endif
827  qflx_sgs_momx(k,i,j,xdir) = dens(k,i,j) * ( &
828  - 2.0_rp * nu(k,i,j) &
829  * ( s11_c(k,i,j) - ( s11_c(k,i,j) + s22_c(k,i,j) + s33_c(k,i,j) ) * oneoverthree ) &
830  + twooverthree * tke(k,i,j) * tke_fact )
831  enddo
832  enddo
833  enddo
834 #ifdef DEBUG
835  i = iundef; j = iundef; k = iundef
836 #endif
837  ! (z edge)
838  !$omp parallel do default(none) &
839  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,DENS,nu,S12_Z,qflx_sgs_momx) &
840  !$omp private(i,j,k) OMP_SCHEDULE_ collapse(2)
841  do j = jjs-1, jje
842  do i = iis, iie
843  do k = ks, ke
844 #ifdef DEBUG
845  call check( __line__, dens(k,i,j) )
846  call check( __line__, dens(k,i+1,j) )
847  call check( __line__, dens(k,i,j+1) )
848  call check( __line__, dens(k,i+1,j+1) )
849  call check( __line__, nu(k,i,j) )
850  call check( __line__, nu(k,i+1,j) )
851  call check( __line__, nu(k,i,j+1) )
852  call check( __line__, nu(k,i+1,j+1) )
853  call check( __line__, s12_z(k,i,j) )
854 #endif
855  qflx_sgs_momx(k,i,j,ydir) = - 0.125_rp & ! 2/4/4
856  * ( dens(k,i,j)+dens(k,i+1,j)+dens(k,i,j+1)+dens(k,i+1,j+1) ) &
857  * ( nu(k,i,j)+nu(k,i+1,j)+nu(k,i,j+1)+nu(k,i+1,j+1) ) &
858  * s12_z(k,i,j)
859  enddo
860  enddo
861  enddo
862 #ifdef DEBUG
863  i = iundef; j = iundef; k = iundef
864 #endif
865 
866  if ( atmos_phy_tb_smg_implicit ) then
867  call calc_tend_momx( tend, & ! (out)
868  qflx_sgs_momx, & ! (in)
869  gsqrt, j13g, j23g, j33g, mapf, & ! (in)
870  iis, iie, jjs, jje ) ! (in)
871 
872  !$omp parallel do &
873  !$omp private(ap,d)
874  do j = jjs, jje
875  do i = iis, iie
876 
877  ap = - dt * 0.25_rp * ( dens(ks ,i ,j)*nu(ks ,i ,j) &
878  + dens(ks+1,i ,j)*nu(ks+1,i ,j) &
879  + dens(ks ,i+1,j)*nu(ks ,i+1,j) &
880  + dens(ks+1,i+1,j)*nu(ks+1,i+1,j) ) &
881  * rfdz(ks) / gsqrt(ks,i,j,i_uyw)
882  a(ks,i,j) = ap * rcdz(ks) / gsqrt(ks,i,j,i_uyz)
883  c(ks,i,j) = 0.0_rp
884  b(ks,i,j) = - a(ks,i,j) + 0.5_rp * ( dens(ks,i,j)+dens(ks,i+1,j) )
885  do k = ks+1, ke-1
886  c(k,i,j) = ap * rcdz(k) / gsqrt(k,i,j,i_uyz)
887  ap = - dt * 0.25_rp * ( dens(k ,i ,j)*nu(k ,i ,j) &
888  + dens(k+1,i ,j)*nu(k+1,i ,j) &
889  + dens(k ,i+1,j)*nu(k ,i+1,j) &
890  + dens(k+1,i+1,j)*nu(k+1,i+1,j) ) &
891  * rfdz(k) / gsqrt(k,i,j,i_uyw)
892  a(k,i,j) = ap * rcdz(k) / gsqrt(k,i,j,i_uyz)
893  b(k,i,j) = - a(k,i,j) - c(k,i,j) + 0.5_rp * ( dens(k,i,j)+dens(k,i+1,j) )
894  end do
895  a(ke,i,j) = 0.0_rp
896  c(ke,i,j) = ap * rcdz(ke) / gsqrt(ke,i,j,i_uyz)
897  b(ke,i,j) = - c(ke,i,j) + 0.5_rp * ( dens(ke,i,j)+dens(ke,i+1,j) )
898 
899  do k = ks, ke
900  d(k) = tend(k,i,j)
901  end do
902 
903  call diffusion_solver( &
904  tend(:,i,j), & ! (out)
905  a(:,i,j), b(:,i,j), c(:,i,j), d, & ! (in)
906  ke ) ! (in)
907 
908  do k = ks, ke-1
909  qflx_sgs_momx(k,i,j,zdir) = qflx_sgs_momx(k,i,j,zdir) &
910  - 0.25_rp * ( dens(k ,i ,j)*nu(k ,i ,j) &
911  + dens(k+1,i ,j)*nu(k+1,i ,j) &
912  + dens(k ,i+1,j)*nu(k ,i+1,j) &
913  + dens(k+1,i+1,j)*nu(k+1,i+1,j) ) &
914  * dt * ( tend(k+1,i,j) - tend(k,i,j) ) * rfdz(k) / gsqrt(k,i,j,i_uyw)
915  end do
916 
917  end do
918  end do
919 
920  end if
921 
922  !##### momentum equation (y) #####
923  ! (x edge)
924  if ( atmos_phy_tb_smg_horizontal ) then
925  qflx_sgs_momy(:,:,:,zdir) = 0.0_rp
926  else
927  !$omp parallel do
928  do j = jjs, jje
929  do i = iis, iie
930  do k = ks, ke-1
931 #ifdef DEBUG
932  call check( __line__, dens(k,i,j) )
933  call check( __line__, dens(k,i,j+1) )
934  call check( __line__, dens(k+1,i,j) )
935  call check( __line__, dens(k+1,i,j+1) )
936  call check( __line__, nu(k,i,j) )
937  call check( __line__, nu(k,i,j+1) )
938  call check( __line__, nu(k+1,i,j) )
939  call check( __line__, nu(k+1,i,j+1) )
940  call check( __line__, s23_x(k,i,j) )
941 #endif
942  qflx_sgs_momy(k,i,j,zdir) = - 0.125_rp & ! 2/4/4
943  * ( dens(k,i,j)+dens(k+1,i,j)+dens(k,i,j+1)+dens(k+1,i,j+1) ) &
944  * ( nu(k,i,j)+nu(k+1,i,j)+nu(k,i,j+1)+nu(k+1,i,j+1) ) &
945  * s23_x(k,i,j)
946  enddo
947  enddo
948  enddo
949 #ifdef DEBUG
950  i = iundef; j = iundef; k = iundef
951 #endif
952  !$omp parallel do
953  do j = jjs, jje
954  do i = iis, iie
955  qflx_sgs_momy(ks-1,i,j,zdir) = 0.0_rp ! bottom boundary
956  qflx_sgs_momy(ke ,i,j,zdir) = 0.0_rp ! top boundary
957  enddo
958  enddo
959 #ifdef DEBUG
960  i = iundef; j = iundef; k = iundef
961 #endif
962  end if
963 
964  ! (z edge)
965  !$omp parallel do default(none) &
966  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,DENS,nu,S12_Z,qflx_sgs_momy) &
967  !$omp private(i,j,k) OMP_SCHEDULE_ collapse(2)
968  do j = jjs, jje
969  do i = iis-1, iie
970  do k = ks, ke
971 #ifdef DEBUG
972  call check( __line__, dens(k,i,j) )
973  call check( __line__, dens(k,i+1,j) )
974  call check( __line__, dens(k,i,j+1) )
975  call check( __line__, dens(k,i+1,j+1) )
976  call check( __line__, nu(k,i,j) )
977  call check( __line__, nu(k,i+1,j) )
978  call check( __line__, nu(k,i,j+1) )
979  call check( __line__, nu(k,i+1,j+1) )
980  call check( __line__, s12_z(k,i,j) )
981 #endif
982  qflx_sgs_momy(k,i,j,xdir) = - 0.125_rp & !
983  * ( dens(k,i,j)+dens(k,i+1,j)+dens(k,i,j+1)+dens(k,i+1,j+1) ) &
984  * ( nu(k,i,j)+nu(k,i+1,j)+nu(k,i,j+1)+nu(k,i+1,j+1) ) &
985  * s12_z(k,i,j)
986  enddo
987  enddo
988  enddo
989 #ifdef DEBUG
990  i = iundef; j = iundef; k = iundef
991 #endif
992 
993  ! (z-x plane)
994  !$omp parallel do default(none) &
995  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,DENS,nu,S11_C,S22_C,S33_C,tke_fact,TKE,qflx_sgs_momy) &
996  !$omp private(i,j,k) OMP_SCHEDULE_ collapse(2)
997  do j = jjs, jje+1
998  do i = iis, iie
999  do k = ks, ke
1000 #ifdef DEBUG
1001  call check( __line__, dens(k,i,j) )
1002  call check( __line__, nu(k,i,j) )
1003  call check( __line__, s11_c(k,i,j) )
1004  call check( __line__, s22_c(k,i,j) )
1005  call check( __line__, s33_c(k,i,j) )
1006  call check( __line__, tke(k,i,j) )
1007 #endif
1008  qflx_sgs_momy(k,i,j,ydir) = dens(k,i,j) * ( &
1009  - 2.0_rp * nu(k,i,j) &
1010  * ( s22_c(k,i,j) - ( s11_c(k,i,j) + s22_c(k,i,j) + s33_c(k,i,j) ) * oneoverthree ) &
1011  + twooverthree * tke(k,i,j) * tke_fact)
1012  enddo
1013  enddo
1014  enddo
1015 #ifdef DEBUG
1016  i = iundef; j = iundef; k = iundef
1017 #endif
1018 
1019  if ( atmos_phy_tb_smg_implicit ) then
1020  call calc_tend_momy( tend, & ! (out)
1021  qflx_sgs_momy, & ! (in)
1022  gsqrt, j13g, j23g, j33g, mapf, & ! (in)
1023  iis, iie, jjs, jje ) ! (in)
1024 
1025  !$omp parallel do &
1026  !$omp private(ap,d)
1027  do j = jjs, jje
1028  do i = iis, iie
1029 
1030  ap = - dt * 0.25_rp * ( dens(ks ,i,j )*nu(ks ,i,j ) &
1031  + dens(ks+1,i,j )*nu(ks+1,i,j ) &
1032  + dens(ks ,i,j+1)*nu(ks ,i,j+1) &
1033  + dens(ks+1,i,j+1)*nu(ks+1,i,j+1) ) &
1034  * rfdz(ks) / gsqrt(ks,i,j,i_xvw)
1035  a(ks,i,j) = ap * rcdz(ks) / gsqrt(ks,i,j,i_xvz)
1036  c(ks,i,j) = 0.0_rp
1037  b(ks,i,j) = - a(ks,i,j) + 0.5_rp * ( dens(ks,i,j)+dens(ks,i,j+1) )
1038  do k = ks+1, ke-1
1039  c(k,i,j) = ap * rcdz(k) / gsqrt(k,i,j,i_xvz)
1040  ap = - dt * 0.25_rp * ( dens(k ,i,j )*nu(k ,i,j ) &
1041  + dens(k+1,i,j )*nu(k+1,i,j ) &
1042  + dens(k ,i,j+1)*nu(k ,i,j+1) &
1043  + dens(k+1,i,j+1)*nu(k+1,i,j+1) ) &
1044  * rfdz(k) / gsqrt(k,i,j,i_xvw)
1045  a(k,i,j) = ap * rcdz(k) / gsqrt(k,i,j,i_xvz)
1046  b(k,i,j) = - a(k,i,j) - c(k,i,j) + 0.5_rp * ( dens(k,i,j)+dens(k,i,j+1) )
1047  end do
1048  a(ke,i,j) = 0.0_rp
1049  c(ke,i,j) = ap * rcdz(ke) / gsqrt(ke,i,j,i_xvz)
1050  b(ke,i,j) = - c(ke,i,j) + 0.5_rp * ( dens(ke,i,j)+dens(ke,i,j+1) )
1051 
1052  do k = ks, ke
1053  d(k) = tend(k,i,j)
1054  end do
1055 
1056  call diffusion_solver( &
1057  tend(:,i,j), & ! (out)
1058  a(:,i,j), b(:,i,j), c(:,i,j), d, & ! (in)
1059  ke ) ! (in)
1060 
1061  do k = ks, ke-1
1062  qflx_sgs_momy(k,i,j,zdir) = qflx_sgs_momy(k,i,j,zdir) &
1063  - 0.25_rp * ( dens(k ,i,j )*nu(k ,i,j ) &
1064  + dens(k+1,i,j )*nu(k+1,i,j ) &
1065  + dens(k ,i,j+1)*nu(k ,i,j+1) &
1066  + dens(k+1,i,j+1)*nu(k+1,i,j+1) ) &
1067  * dt * ( tend(k+1,i,j) - tend(k,i,j) ) * rfdz(k) / gsqrt(k,i,j,i_xvw)
1068  end do
1069 
1070  end do
1071  end do
1072 
1073  end if
1074 
1075  if ( atmos_phy_tb_smg_backscatter ) then
1076 
1077  ! MOMZ : dfy/dx - dfx/dy
1078  !$omp parallel do
1079  do j = jjs, jje
1080  do i = iis, iie
1081  do k = ks+1, ke-1
1082  acc(k,i,j) = ( ( gsqrt(k,i+1,j,i_xyz) * random_my(k,i+1,j) - gsqrt(k,i-1,j,i_xyz) * random_my(k,i-1,j) ) / ( fdx(i) + fdx(i-1) ) * mapf(i,j,1,i_xy) &
1083  + ( j13g(k+1,i,j,i_xyz) * random_my(k+1,i,j) - j13g(k-1,i,j,i_xyz) * random_my(k-1,i,j) ) / ( fdz(k) + fdz(k-1) ) &
1084  - ( gsqrt(k,i,j+1,i_xyz) * random_mx(k,i,j+1) - gsqrt(k,i,j-1,i_xyz) * random_mx(k,i,j-1) ) / ( fdy(j) + fdy(j-1) ) * mapf(i,j,2,i_xy) &
1085  + ( j23g(k+1,i,j,i_xyz) * random_mx(k+1,i,j) - j23g(k-1,i,j,i_xyz) * random_mx(k-1,i,j) ) / ( fdz(k) + fdz(k-1) ) &
1086  ) / gsqrt(k,i,j,i_xyz)
1087  end do
1088  acc(ks,i,j) = ( ( gsqrt(ks,i+1,j,i_xyz) * random_my(ks,i+1,j) - gsqrt(ks,i-1,j,i_xyz) * random_my(ks,i-1,j) ) / ( fdx(i) + fdx(i-1) ) * mapf(i,j,1,i_xy) &
1089  + ( j13g(ks+1,i,j,i_xyz) * random_my(ks+1,i,j) - j13g(ks,i,j,i_xyz) * random_my(ks,i,j) ) * rfdz(ks) * 0.5_rp &
1090  - ( gsqrt(ks,i,j+1,i_xyz) * random_mx(ks,i,j+1) - gsqrt(ks,i,j-1,i_xyz) * random_mx(ks,i,j-1) ) / ( fdy(j) + fdy(j-1) ) * mapf(i,j,2,i_xy) &
1091  + ( j23g(ks+1,i,j,i_xyz) * random_mx(ks+1,i,j) - j23g(ks,i,j,i_xyz) * random_mx(ks,i,j) ) * rfdz(ks) * 0.5_rp &
1092  ) / gsqrt(ks,i,j,i_xyz)
1093  acc(ke,i,j) = ( ( gsqrt(ke,i+1,j,i_xyz) * random_my(ke,i+1,j) - gsqrt(ke,i-1,j,i_xyz) * random_my(ke,i-1,j) ) / ( fdx(i) + fdx(i-1) ) * mapf(i,j,1,i_xy) &
1094  + ( j13g(ke,i,j,i_xyz) * random_my(ke,i,j) - j13g(ke-1,i,j,i_xyz) * random_my(ke-1,i,j) ) * rfdz(ke-1) * 0.5_rp &
1095  - ( gsqrt(ke,i,j+1,i_xyz) * random_mx(ke,i,j+1) - gsqrt(ke,i,j-1,i_xyz) * random_mx(ke,i,j-1) ) / ( fdy(j) + fdy(j-1) ) * mapf(i,j,2,i_xy) &
1096  + ( j23g(ke,i,j,i_xyz) * random_mx(ke,i,j) - j23g(ke-1,i,j,i_xyz) * random_mx(ke-1,i,j) ) * rfdz(ke-1) * 0.5_rp &
1097  ) / gsqrt(ke,i,j,i_xyz)
1098  end do
1099  end do
1100 
1101 #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) ) )
1102  !$omp parallel do
1103  do j = jjs, jje
1104  do i = iis, iie
1105  do k = ks, ke-1
1106  momz_t(k,i,j) = acc(k+1,i,j) * f2h(k,i,j,1) + acc(k,i,j) * f2h(k,i,j,2)
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  do j = jjs, jje
1115  do i = iis, iie+1
1116  do k = ks+1, ke-1
1117  acc(k,i,j) = ( ( gsqrt(k,i,j+1,i_xyz) * random_mz(k,i,j+1) - gsqrt(k,i,j-1,i_xyz) * random_mz(k,i,j-1) ) / ( fdy(j) + fdy(j-1) ) * mapf(i,j,2,i_xy) &
1118  + ( j23g(k+1,i,j,i_xyz) * random_mz(k+1,i,j) - j23g(k-1,i,j,i_xyz) * random_mz(k-1,i,j) ) / ( fdz(k) + fdz(k-1) ) &
1119  - j33g * ( random_my(k+1,i,j) - random_my(k-1,i,j) ) / ( fdz(k) + fdz(k-1) ) &
1120  ) / gsqrt(k,i,j,i_xyz)
1121  end do
1122  acc(ks,i,j) = ( ( gsqrt(ks,i,j+1,i_xyz) * random_mz(ks,i,j+1) - gsqrt(ks,i,j-1,i_xyz) * random_mz(ks,i,j-1) ) / ( fdy(j) + fdy(j-1) ) * mapf(i,j,2,i_xy) &
1123  + ( j23g(ks+1,i,j,i_xyz) * random_mz(ks+1,i,j) - j23g(ks,i,j,i_xyz) * random_mz(ks,i,j) ) * rfdz(ks) * 0.5_rp &
1124  - j33g * ( random_my(ks+1,i,j) - random_my(ks,i,j) ) * rfdz(ks) * 0.5_rp &
1125  ) / gsqrt(ks,i,j,i_xyz)
1126  acc(ke,i,j) = ( ( gsqrt(ke,i,j+1,i_xyz) * random_mz(ke,i,j+1) - gsqrt(ke,i,j-1,i_xyz) * random_mz(ke,i,j-1) ) / ( fdy(j) + fdy(j-1) ) * mapf(i,j,2,i_xy) &
1127  + ( j23g(ke,i,j,i_xyz) * random_mz(ke,i,j) - j23g(ke-1,i,j,i_xyz) * random_mz(ke-1,i,j) ) * rfdz(ke-1) * 0.5_rp &
1128  - j33g * ( random_my(ke,i,j) - random_my(ke-1,i,j) ) * rfdz(ke-1) * 0.5_rp &
1129  ) / gsqrt(ke,i,j,i_xyz)
1130  end do
1131  end do
1132 
1133  !$omp parallel do
1134  do j = jjs, jje
1135  do i = iis, iie
1136  do k = ks, ke
1137  momx_t(k,i,j) = ( acc(k,i+1,j) + acc(k,i,j) ) * 0.5_rp
1138  end do
1139  end do
1140  end do
1141 
1142  ! MOMY : dfx/dz - dfz/dx
1143  !$omp parallel do
1144  do j = jjs, jje+1
1145  do i = iis, iie
1146  do k = ks+1, ke-1
1147  acc(k,i,j) = ( j33g * ( random_mx(k+1,i,j) - random_mx(k-1,i,j) ) / ( fdz(k) + fdz(k-1) ) &
1148  - ( gsqrt(k,i+1,j,i_xyz) * random_mz(k,i+1,j) - gsqrt(k,i-1,j,i_xyz) * random_mz(k,i-1,j) ) / ( fdx(i) + fdx(i-1) ) * mapf(i,j,1,i_xy) &
1149  + ( j13g(k+1,i,j,i_xyz) * random_mz(k+1,i,j) - j13g(k-1,i,j,i_xyz) * random_mz(k-1,i,j) ) / ( fdz(k) + fdz(k-1) ) &
1150  ) / gsqrt(k,i,j,i_xyz)
1151  end do
1152  acc(ks,i,j) = ( j33g * ( random_mx(ks+1,i,j) - random_mx(ks,i,j) ) * rfdz(ks) * 0.5_rp &
1153  - ( gsqrt(ks,i+1,j,i_xyz) * random_mz(ks,i+1,j) - gsqrt(ks,i-1,j,i_xyz) * random_mz(ks,i-1,j) ) / ( fdx(i) + fdx(i-1) ) * mapf(i,j,1,i_xy) &
1154  + ( j13g(ks+1,i,j,i_xyz) * random_mz(ks+1,i,j) - j13g(ks,i,j,i_xyz) * random_mz(ks,i,j) ) * rfdz(ks) * 0.5_rp &
1155  ) / gsqrt(ks,i,j,i_xyz)
1156  acc(ke,i,j) = ( j33g * ( random_mx(ke,i,j) - random_mx(ke-1,i,j) ) * rfdz(ke-1) * 0.5_rp &
1157  - ( gsqrt(ke,i+1,j,i_xyz) * random_mz(ke,i+1,j) - gsqrt(ke,i-1,j,i_xyz) * random_mz(ke,i-1,j) ) / ( fdx(i) + fdx(i-1) ) * mapf(i,j,1,i_xy) &
1158  + ( j13g(ke,i,j,i_xyz) * random_mz(ke,i,j) - j13g(ke-1,i,j,i_xyz) * random_mz(ke-1,i,j) ) * rfdz(ke-1) * 0.5_rp &
1159  ) / gsqrt(ke,i,j,i_xyz)
1160  end do
1161  end do
1162 
1163  !$omp parallel do
1164  do j = jjs, jje
1165  do i = iis, iie
1166  do k = ks, ke
1167  momy_t(k,i,j) = ( acc(k,i,j+1) + acc(k,i,j) ) * 0.5_rp
1168  end do
1169  end do
1170  end do
1171 
1172  else
1173 
1174  !$omp parallel do
1175 !OCL XFILL
1176  do j = jjs, jje
1177  do i = iis, iie
1178  do k = ks, ke
1179  momz_t(k,i,j) = 0.0_rp
1180  end do
1181  end do
1182  end do
1183  !$omp parallel do
1184 !OCL XFILL
1185  do j = jjs, jje
1186  do i = iis, iie
1187  do k = ks, ke
1188  momx_t(k,i,j) = 0.0_rp
1189  end do
1190  end do
1191  end do
1192  !$omp parallel do
1193 !OCL XFILL
1194  do j = jjs, jje
1195  do i = iis, iie
1196  do k = ks, ke
1197  momy_t(k,i,j) = 0.0_rp
1198  end do
1199  end do
1200  end do
1201 
1202  end if
1203 
1204  !##### Thermodynamic Equation #####
1205 
1206  if ( atmos_phy_tb_smg_implicit ) then
1207 
1208  !$omp parallel do &
1209  !$omp private (ap,d)
1210  do j = jjs, jje
1211  do i = iis, iie
1212 
1213  ap = - dt * 0.25_rp * ( dens(ks,i,j)+dens(ks+1,i,j) ) &
1214  * ( kh(ks,i,j)+kh(ks+1,i,j) ) &
1215  * rfdz(ks) / gsqrt(ks,i,j,i_xyw)
1216  a(ks,i,j) = ap * rcdz(ks) / gsqrt(ks,i,j,i_xyz)
1217  c(ks,i,j) = 0.0_rp
1218  b(ks,i,j) = - a(ks,i,j) + dens(ks,i,j)
1219  do k = ks+1, ke-1
1220  c(k,i,j) = ap * rcdz(k) / gsqrt(k,i,j,i_xyz)
1221  ap = - dt * 0.25_rp * ( dens(k,i,j)+dens(k+1,i,j) ) &
1222  * ( kh(k,i,j)+kh(k+1,i,j) ) &
1223  * rfdz(k) / gsqrt(k,i,j,i_xyw)
1224  a(k,i,j) = ap * rcdz(k) / gsqrt(k,i,j,i_xyz)
1225  b(k,i,j) = - a(k,i,j) - c(k,i,j) + dens(k,i,j)
1226  end do
1227  a(ke,i,j) = 0.0_rp
1228  c(ke,i,j) = ap * rcdz(ke) / gsqrt(ke,i,j,i_xyz)
1229  b(ke,i,j) = - c(ke,i,j) + dens(ke,i,j)
1230 
1231  end do
1232  end do
1233 
1234  end if
1235 
1236  call calc_flux_phi( &
1237  qflx_sgs_rhot, &
1238  dens, pott, kh, 1.0_rp, &
1239  gsqrt, j13g, j23g, j33g, mapf, &
1240  atmos_phy_tb_smg_horizontal, &
1241  atmos_phy_tb_smg_implicit, &
1242  a, b, c, dt, &
1243  iis, iie, jjs, jje )
1244 
1245  if ( atmos_phy_tb_smg_backscatter ) then
1246 
1247  !$omp parallel do
1248  do j = jjs-1, jje+1
1249  do i = iis-1, iie+1
1250 
1251  do k = ks+1, ke-1
1252  dd(k,i,j) = sqrt( ( ( pott(k+1,i,j) - pott(k-1,i,j) ) * j33g / ( fdz(k) + fdz(k-1) ) )**2 &
1253  + ( ( ( 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) ) &
1254  + ( 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 &
1255  + ( ( ( 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) ) &
1256  + ( 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 &
1257  ) / gsqrt(k,i,j,i_xyz)
1258  end do
1259  dd(ks,i,j) = sqrt( ( ( pott(ks+1,i,j) - pott(ks,i,j) ) * j33g * rfdz(ks) )**2 &
1260  + ( ( ( 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) ) &
1261  + ( 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 &
1262  + ( ( ( 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) ) &
1263  + ( 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 &
1264  ) / gsqrt(ks,i,j,i_xyz)
1265  dd(ke,i,j) = sqrt( ( ( pott(ke,i,j) - pott(ke-1,i,j) ) * j33g * rfdz(ke-1) )**2 &
1266  + ( ( ( 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) ) &
1267  + ( 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 &
1268  + ( ( ( 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) ) &
1269  + ( 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 &
1270  ) / gsqrt(ke,i,j,i_xyz)
1271 
1272  end do
1273  end do
1274 
1275  !$omp parallel do
1276  do j = jjs, jje
1277  do i = iis, iie
1278  do k = ks+1, ke-1
1279  rhot_t(k,i,j) = ( j33g * ( random_qz(k+1,i,j) * dd(k+1,i,j) - random_qz(k-1,i,j) * dd(k-1,i,j) ) / ( fdz(k) + fdz(k-1) ) &
1280  + ( 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) &
1281  + ( j13g(k+1,i,j,i_xyz) * random_qx(k+1,i,j) * dd(k+1,i,j) - j13g(k-1,i,j,i_xyz) * random_qx(k-1,i,j) * dd(k-1,i,j) ) / ( fdz(k) + fdz(k-1) ) &
1282  + ( 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) &
1283  + ( j23g(k+1,i,j,i_xyz) * random_qy(k+1,i,j) * dd(k+1,i,j) - j23g(k-1,i,j,i_xyz) * random_qy(k-1,i,j) * dd(k-1,i,j) ) / ( fdz(k) + fdz(k-1) ) &
1284  ) / gsqrt(k,i,j,i_xyz)
1285  end do
1286  rhot_t(ks,i,j) = ( j33g * ( random_qz(ks+1,i,j) * dd(ks+1,i,j) - random_qz(ks,i,j) * dd(ks,i,j) ) * rfdz(ks) &
1287  + ( gsqrt(ks,i+1,j,i_xyz) * random_qx(ks,i+1,j) * dd(ks,i+1,j) - gsqrt(ks,i-1,j,i_xyz) * random_qx(ks,i-1,j) * dd(ks,i-1,j) ) / ( fdx(i) + fdx(i-1) ) * mapf(i,j,1,i_xy) &
1288  + ( j13g(ks+1,i,j,i_xyz) * random_qx(ks+1,i,j) * dd(ks+1,i,j) - j13g(ks,i,j,i_xyz) * random_qx(ks,i,j) * dd(ks,i,j) ) * rfdz(ks) * 0.5_rp &
1289  + ( gsqrt(ks,i,j+1,i_xyz) * random_qy(ks,i,j+1) * dd(ks,i,j+1) - gsqrt(ks,i,j-1,i_xyz) * random_qy(ks,i,j-1) * dd(ks,i,j-1) ) / ( fdy(j) + fdy(j-1) ) * mapf(i,j,2,i_xy) &
1290  + ( j23g(ks+1,i,j,i_xyz) * random_qy(ks+1,i,j) * dd(ks+1,i,j) - j23g(ks,i,j,i_xyz) * random_qy(ks,i,j) * dd(ks,i,j) ) * rfdz(ks) * 0.5_rp &
1291  ) / gsqrt(ks,i,j,i_xyz)
1292  rhot_t(ke,i,j) = ( j33g * ( random_qz(ke,i,j) * dd(ke,i,j) - random_qz(ke-1,i,j) * dd(ke-1,i,j) ) * rfdz(ke-1) * 0.5_rp &
1293  + ( gsqrt(ke,i+1,j,i_xyz) * random_qx(ke,i+1,j) * dd(ke,i+1,j) - gsqrt(ke,i-1,j,i_xyz) * random_qx(ke,i-1,j) * dd(ke,i-1,j) ) / ( fdx(i) + fdx(i-1) ) * mapf(i,j,1,i_xy) &
1294  + ( j13g(ke,i,j,i_xyz) * random_qx(ke,i,j) * dd(ke,i,j) - j13g(ke-1,i,j,i_xyz) * random_qx(ke-1,i,j) * dd(ke-1,i,j) ) * rfdz(ke-1) * 0.5_rp &
1295  + ( gsqrt(ke,i,j+1,i_xyz) * random_qy(ke,i,j+1) * dd(ke,i,j+1) - gsqrt(ke,i,j-1,i_xyz) * random_qy(ke,i,j-1) * dd(ke,i,j-1) ) / ( fdy(j) + fdy(j-1) ) * mapf(i,j,2,i_xy) &
1296  + ( j23g(ke,i,j,i_xyz) * random_qy(ke,i,j) * dd(ke,i,j) - j23g(ke-1,i,j,i_xyz) * random_qy(ke-1,i,j) * dd(ke-1,i,j) ) * rfdz(ke-1) * 0.5_rp &
1297  ) / gsqrt(ke,i,j,i_xyz)
1298  end do
1299  end do
1300 
1301  else
1302 
1303  !$omp parallel do
1304 !OCL XFILL
1305  do j = jjs, jje
1306  do i = iis, iie
1307  do k = ks, ke
1308  rhot_t(k,i,j) = 0.0_rp
1309  end do
1310  end do
1311  end do
1312 
1313  end if
1314 
1315  enddo
1316  enddo
1317 
1318  !##### Tracers #####
1319  do iq = 1, qa
1320 
1321  if ( .not. tracer_advc(iq) ) cycle
1322 
1323  do jjs = js, je, jblock
1324  jje = jjs+jblock-1
1325  do iis = is, ie, iblock
1326  iie = iis+iblock-1
1327 
1328  call calc_flux_phi( &
1329  qflx_sgs_rhoq(:,:,:,:,iq), &
1330  dens, qtrc(:,:,:,iq), kh, 1.0_rp, &
1331  gsqrt, j13g, j23g, j33g, mapf, &
1332  atmos_phy_tb_smg_horizontal, &
1333  atmos_phy_tb_smg_implicit, &
1334  a, b, c, dt, &
1335  iis, iie, jjs, jje )
1336 
1337 
1338  if ( atmos_phy_tb_smg_backscatter .and. iq == i_qv ) then
1339 
1340  !$omp parallel do
1341  do j = jjs-1, jje+1
1342  do i = iis-1, iie+1
1343 
1344  do k = ks+1, ke-1
1345  dd(k,i,j) = sqrt( ( ( qtrc(k+1,i,j,iq) - qtrc(k-1,i,j,iq) ) * j33g / ( fdz(k) + fdz(k-1) ) )**2 &
1346  + ( ( ( 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) ) &
1347  + ( 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 &
1348  + ( ( ( 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) ) &
1349  + ( 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 &
1350  ) / gsqrt(k,i,j,i_xyz)
1351  end do
1352  dd(ks,i,j) = sqrt( ( ( qtrc(ks+1,i,j,iq) - qtrc(ks,i,j,iq) ) * j33g * rfdz(ks) )**2 &
1353  + ( ( ( 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) ) &
1354  + ( 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 &
1355  + ( ( ( 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) ) &
1356  + ( 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 &
1357  ) / gsqrt(ks,i,j,i_xyz)
1358  dd(ke,i,j) = sqrt( ( ( qtrc(ke,i,j,iq) - qtrc(ke-1,i,j,iq) ) * j33g * rfdz(ke-1) )**2 &
1359  + ( ( ( 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) ) &
1360  + ( 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 &
1361  + ( ( ( 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) ) &
1362  + ( 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 &
1363  ) / gsqrt(ke,i,j,i_xyz)
1364 
1365  end do
1366  end do
1367 
1368  !$omp parallel do
1369  do j = jjs, jje
1370  do i = iis, iie
1371  do k = ks+1, ke-1
1372  rhoq_t(k,i,j,iq) = ( j33g * ( random_qz(k+1,i,j) * dd(k+1,i,j) - random_qz(k-1,i,j) * dd(k-1,i,j) ) / ( fdz(k) + fdz(k-1) ) &
1373  + ( 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) &
1374  + ( j13g(k+1,i,j,i_xyz) * random_qx(k+1,i,j) * dd(k+1,i,j) - j13g(k-1,i,j,i_xyz) * random_qx(k-1,i,j) * dd(k-1,i,j) ) / ( fdz(k) + fdz(k-1) ) &
1375  + ( 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) &
1376  + ( j23g(k+1,i,j,i_xyz) * random_qy(k+1,i,j) * dd(k+1,i,j) - j23g(k-1,i,j,i_xyz) * random_qy(k-1,i,j) * dd(k-1,i,j) ) / ( fdz(k) + fdz(k-1) ) &
1377  ) / gsqrt(k,i,j,i_xyz)
1378  end do
1379  rhoq_t(ks,i,j,iq) = ( j33g * ( random_qz(ks+1,i,j) * dd(ks+1,i,j) - random_qz(ks,i,j) * dd(ks,i,j) ) * rfdz(ks) * 0.5_rp &
1380  + ( gsqrt(ks,i+1,j,i_xyz) * random_qx(ks,i+1,j) * dd(ks,i+1,j) - gsqrt(ks,i-1,j,i_xyz) * random_qx(ks,i-1,j) * dd(ks,i-1,j) ) / ( fdx(i) + fdx(i-1) ) * mapf(i,j,1,i_xy) &
1381  + ( j13g(ks+1,i,j,i_xyz) * random_qx(ks+1,i,j) * dd(ks+1,i,j) - j13g(ks,i,j,i_xyz) * random_qx(ks,i,j) * dd(ks,i,j) ) * rfdz(ks) * 0.5_rp &
1382  + ( gsqrt(ks,i,j+1,i_xyz) * random_qy(ks,i,j+1) * dd(ks,i,j+1) - gsqrt(ks,i,j-1,i_xyz) * random_qy(ks,i,j-1) * dd(ks,i,j-1) ) / ( fdy(j) + fdy(j-1) ) * mapf(i,j,2,i_xy) &
1383  + ( j23g(ks+1,i,j,i_xyz) * random_qy(ks+1,i,j) * dd(ks+1,i,j) - j23g(ks,i,j,i_xyz) * random_qy(ks,i,j) * dd(ks,i,j) ) * rfdz(ks) * 0.5_rp &
1384  ) / gsqrt(ks,i,j,i_xyz)
1385  rhoq_t(ke,i,j,iq) = ( j33g * ( random_qz(ke,i,j) * dd(ke,i,j) - random_qz(ke-1,i,j) * dd(ke-1,i,j) ) * rfdz(ke-1) * 0.5_rp &
1386  + ( gsqrt(ke,i+1,j,i_xyz) * random_qx(ke,i+1,j) * dd(ke,i+1,j) - gsqrt(ke,i-1,j,i_xyz) * random_qx(ke,i-1,j) * dd(ke,i-1,j) ) / ( fdx(i) + fdx(i-1) ) * mapf(i,j,1,i_xy) &
1387  + ( j13g(ke,i,j,i_xyz) * random_qx(ke,i,j) * dd(ke,i,j) - j13g(ke-1,i,j,i_xyz) * random_qx(ke-1,i,j) * dd(ke-1,i,j) ) * rfdz(ke-1) * 0.5_rp &
1388  + ( gsqrt(ke,i,j+1,i_xyz) * random_qy(ke,i,j+1) * dd(ke,i,j+1) - gsqrt(ke,i,j-1,i_xyz) * random_qy(ke,i,j-1) * dd(ke,i,j-1) ) / ( fdy(j) + fdy(j-1) ) * mapf(i,j,2,i_xy) &
1389  + ( j23g(ke,i,j,i_xyz) * random_qy(ke,i,j) * dd(ke,i,j) - j23g(ke-1,i,j,i_xyz) * random_qy(ke-1,i,j) * dd(ke-1,i,j) ) * rfdz(ke-1) * 0.5_rp &
1390  ) / gsqrt(ke,i,j,i_xyz)
1391  end do
1392  end do
1393 
1394  else
1395 
1396  !$omp parallel do
1397 !OCL XFILL
1398  do j = jjs, jje
1399  do i = iis, iie
1400  do k = ks, ke
1401  rhoq_t(k,i,j,iq) = 0.0_rp
1402  end do
1403  end do
1404  end do
1405 
1406  end if
1407 
1408  enddo
1409  enddo
1410 
1411  enddo ! scalar quantities loop
1412 #ifdef DEBUG
1413  iq = iundef
1414 #endif
1415 
1416  call file_history_in( tke(:,:,:), 'TKE_SMG', 'turbulent kinetic energy (Smagorinsky)', 'm2/s2', fill_halo=.true. )
1417 
1418 
1419  return
scale_atmos_grid_cartesc_index::i_xyz
integer, public i_xyz
Definition: scale_atmos_grid_cartesC_index.F90:90
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_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_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_grid_cartesc_index::ia
integer, public ia
of whole cells: x, local, with HALO
Definition: scale_atmos_grid_cartesC_index.F90:48
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_xvw
integer, public i_xvw
Definition: scale_atmos_grid_cartesC_index.F90:93
scale_atmos_grid_cartesc_index::i_xvz
integer, public i_xvz
Definition: scale_atmos_grid_cartesC_index.F90:95
scale_atmos_grid_cartesc_index::iblock
integer, public iblock
block size for cache blocking: x
Definition: scale_atmos_grid_cartesC_index.F90:40
scale_tracer::qa
integer, public qa
Definition: scale_tracer.F90:37
scale_tracer::tracer_advc
logical, dimension(qa_max), public tracer_advc
Definition: scale_tracer.F90:45
scale_atmos_grid_cartesc_index::ja
integer, public ja
of whole cells: y, local, with HALO
Definition: scale_atmos_grid_cartesC_index.F90:49
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_atmos_grid_cartesc_index::ie
integer, public ie
end point of inner domain: x, local
Definition: scale_atmos_grid_cartesC_index.F90:54
scale_tracer
module TRACER
Definition: scale_tracer.F90:16
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::ydir
integer, parameter, public ydir
Definition: scale_atmos_grid_cartesC_index.F90:34
scale_atmos_hydrometeor
module atmosphere / hydrometeor
Definition: scale_atmos_hydrometeor.F90:12
scale_atmos_grid_cartesc_index
module atmosphere / grid / cartesC index
Definition: scale_atmos_grid_cartesC_index.F90:12
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::je
integer, public je
end point of inner domain: y, local
Definition: scale_atmos_grid_cartesC_index.F90:56
scale_const::const_grav
real(rp), public const_grav
standard acceleration of gravity [m/s2]
Definition: scale_const.F90:46
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_atmos_grid_cartesc_index::ks
integer, public ks
start point of inner domain: z, local
Definition: scale_atmos_grid_cartesC_index.F90:51
scale_atmos_grid_cartesc_index::jblock
integer, public jblock
block size for cache blocking: y
Definition: scale_atmos_grid_cartesC_index.F90:41
scale_const
module CONSTANT
Definition: scale_const.F90:11
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::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::i_uyw
integer, public i_uyw
Definition: scale_atmos_grid_cartesC_index.F90:92
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_uyz
integer, public i_uyz
Definition: scale_atmos_grid_cartesC_index.F90:94
scale_const::const_eps
real(rp), public const_eps
small number
Definition: scale_const.F90:33
scale_atmos_hydrometeor::i_qv
integer, public i_qv
Definition: scale_atmos_hydrometeor.F90:77
scale_random::random_normal
subroutine, public random_normal(var)
Get normal random number.
Definition: scale_random.F90:153
scale_precision
module PRECISION
Definition: scale_precision.F90:14
scale_atmos_grid_cartesc_index::ka
integer, public ka
of whole cells: z, local, with HALO
Definition: scale_atmos_grid_cartesC_index.F90:47
scale_random
module RANDOM
Definition: scale_random.F90:14
scale_atmos_grid_cartesc_index::zdir
integer, parameter, public zdir
Definition: scale_atmos_grid_cartesC_index.F90:32
scale_file_history
module file_history
Definition: scale_file_history.F90:15
scale_atmos_grid_cartesc_index::i_xyw
integer, public i_xyw
Definition: scale_atmos_grid_cartesC_index.F90:91
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 1424 of file scale_atmos_phy_tb_smg.F90.

Referenced by atmos_phy_tb_smg_setup().

1424  implicit none
1425  real(RP), intent(in) :: dz
1426  real(RP), intent(in) :: dx
1427  real(RP), intent(in) :: dy
1428  real(RP), intent(in) :: filter_fact
1429  real(RP) :: mixlen ! (out)
1430 
1431  mixlen = fact(dz, dx, dy) * filter_fact * ( dz * dx * dy )**oneoverthree ! Scotti et al. (1993)
1432 
1433  return
Here is the caller graph for this function:
Generated by   doxygen 1.8.13