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_finalize
 finalize 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)
    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
    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 102 of file scale_atmos_phy_tb_smg.F90.

102  use scale_prc, only: &
103  prc_abort
104  use scale_const, only: &
105  karman => const_karman
106  implicit none
107 
108  real(RP), intent(in) :: FZ (0:KA,IA,JA)
109  real(RP), intent(in) :: CZ ( KA,IA,JA)
110  real(RP), intent(in) :: CDX (IA)
111  real(RP), intent(in) :: CDY (JA)
112  real(RP), intent(in) :: MAPF(IA,JA,2)
113 
114  logical, intent(in), optional :: horizontal
115 
116  real(RP) :: ATMOS_PHY_TB_SMG_Cs
117  real(RP) :: ATMOS_PHY_TB_SMG_filter_fact
118  logical :: ATMOS_PHY_TB_SMG_consistent_tke
119 
120  namelist / param_atmos_phy_tb_smg / &
121  atmos_phy_tb_smg_cs, &
122  atmos_phy_tb_smg_nu_max, &
123  atmos_phy_tb_smg_filter_fact, &
124  atmos_phy_tb_smg_bottom, &
125  atmos_phy_tb_smg_backscatter, &
126  atmos_phy_tb_smg_implicit, &
127  atmos_phy_tb_smg_consistent_tke, &
128  atmos_phy_tb_smg_horizontal
129 
130  integer :: ierr
131  integer :: k, i, j
132  !---------------------------------------------------------------------------
133 
134  !$acc data copyin(FZ, CZ, CDX, CDY, MAPF)
135 
136  log_newline
137  log_info("ATMOS_PHY_TB_smg_setup",*) 'Setup'
138  log_info("ATMOS_PHY_TB_smg_setup",*) 'Smagorinsky-type Eddy Viscocity Model'
139 
140  atmos_phy_tb_smg_cs = cs
141  atmos_phy_tb_smg_filter_fact = 2.0_rp
142  atmos_phy_tb_smg_consistent_tke = .true.
143 
144  if ( present(horizontal) ) atmos_phy_tb_smg_horizontal = horizontal
145  !--- read namelist
146  rewind(io_fid_conf)
147  read(io_fid_conf,nml=param_atmos_phy_tb_smg,iostat=ierr)
148  if( ierr < 0 ) then !--- missing
149  log_info("ATMOS_PHY_TB_smg_setup",*) 'Not found namelist. Default used.'
150  elseif( ierr > 0 ) then !--- fatal error
151  log_error("ATMOS_PHY_TB_smg_setup",*) 'Not appropriate names in namelist PARAM_ATMOS_PHY_TB_SMG. Check!'
152  call prc_abort
153  endif
154  log_nml(param_atmos_phy_tb_smg)
155 
156  cs = atmos_phy_tb_smg_cs
157 
158  rprn = 1.0_rp / prn
159  rric = 1.0_rp / ric
160  prnovric = ( 1.0_rp - prn ) * rric
161 
162  allocate( lambda0(ka,ia,ja) )
163  allocate( lambda(ka,ia,ja) )
164  !$acc enter data create(lambda0,lambda)
165 
166 #ifdef DEBUG
167  lambda0(:,:,:) = undef
168  lambda(:,:,:) = undef
169 #endif
170  if ( atmos_phy_tb_smg_horizontal ) then
171  !$omp parallel do collapse(2)
172  !$acc kernels
173  do j = js-1, je+1
174  do i = is-1, ie+1
175  do k = ks, ke
176  lambda0(k,i,j) = cs * sqrt( cdx(i) * cdy(j) / ( mapf(i,j,1) * mapf(i,j,2) ) )
177  lambda(k,i,j) = lambda0(k,i,j)
178  enddo
179  enddo
180  enddo
181  !$acc end kernels
182 
183 #ifdef DEBUG
184  i = iundef; j = iundef; k = iundef
185 #endif
186  atmos_phy_tb_smg_consistent_tke = .false.
187  atmos_phy_tb_smg_implicit = .false. ! flux in the z-direction is not necessary
188  atmos_phy_tb_smg_backscatter = .false.
189  else
190  !$omp parallel do collapse(2)
191  !$acc kernels
192  do j = js-1, je+1
193  do i = is-1, ie+1
194  do k = ks, ke
195  lambda0(k,i,j) = cs * mixlen( fz(k,i,j) - fz(k-1,i,j), &
196  cdx(i) / mapf(i,j,1), &
197  cdy(j) / mapf(i,j,2), &
198  atmos_phy_tb_smg_filter_fact )
199  enddo
200  enddo
201  enddo
202  !$acc end kernels
203 #ifdef DEBUG
204  i = iundef; j = iundef; k = iundef
205 #endif
206  if ( atmos_phy_tb_smg_bottom ) then
207  !$omp parallel do collapse(2)
208  !$acc kernels
209  do j = js-1, je+1
210  do i = is-1, ie+1
211  do k = ks, ke
212  lambda(k,i,j) = sqrt( 1.0_rp / ( 1.0_rp / lambda0(k,i,j)**2 + 1.0_rp / ( karman*(cz(k,i,j)-fz(ks-1,i,j)) )**2 ) )
213  enddo
214  enddo
215  enddo
216  !$acc end kernels
217 #ifdef DEBUG
218  i = iundef; j = iundef; k = iundef
219 #endif
220  else
221  !$omp parallel do collapse(2)
222  !$acc kernels
223  do j = js-1, je+1
224  do i = is-1, ie+1
225  do k = ks, ke
226  lambda(k,i,j) = lambda0(k,i,j)
227  enddo
228  enddo
229  enddo
230  !$acc end kernels
231 #ifdef DEBUG
232  i = iundef; j = iundef; k = iundef
233 #endif
234  end if
235  end if
236 
237  ! flag for isotropic stress tensor
238  if ( atmos_phy_tb_smg_consistent_tke ) then
239  tke_fact = 1.0_rp
240  else
241  tke_fact = 0.0_rp ! neglect
242  end if
243 
244  !$acc end data
245 
246  return

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

subroutine, public scale_atmos_phy_tb_smg::atmos_phy_tb_smg_finalize

finalize

Definition at line 252 of file scale_atmos_phy_tb_smg.F90.

252 
253  !$acc exit data delete(lambda0, lambda)
254  deallocate( lambda0 )
255  deallocate( lambda )
256 
257  return

Referenced by mod_atmos_phy_tb_driver::atmos_phy_tb_driver_finalize().

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 269 of file scale_atmos_phy_tb_smg.F90.

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

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_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::ka, 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), value  dz,
real(rp), intent(in), value  dx,
real(rp), intent(in), value  dy,
real(rp), intent(in), value  filter_fact 
)

Definition at line 1604 of file scale_atmos_phy_tb_smg.F90.

1604  !$acc routine seq
1605  implicit none
1606  real(RP), intent(in), value :: dz
1607  real(RP), intent(in), value :: dx
1608  real(RP), intent(in), value :: dy
1609  real(RP), intent(in), value :: filter_fact
1610  real(RP) :: mixlen ! (out)
1611 
1612  mixlen = fact(dz, dx, dy) * filter_fact * ( dz * dx * dy )**oneoverthree ! Scotti et al. (1993)
1613 
1614  return

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:49
scale_atmos_grid_cartesc_index::i_uy
integer, public i_uy
Definition: scale_atmos_grid_cartesC_index.F90:100
scale_atmos_grid_cartesc_index::i_xv
integer, public i_xv
Definition: scale_atmos_grid_cartesC_index.F90:101
scale_prc::prc_abort
subroutine, public prc_abort
Abort Process.
Definition: scale_prc.F90:350
scale_atmos_phy_tb_common::atmos_phy_tb_calc_strain_tensor
subroutine, public atmos_phy_tb_calc_strain_tensor(S33_C, S11_C, S22_C, S31_C, S12_C, S23_C, S12_Z, S23_X, S31_Y, S2, DENS, MOMZ, MOMX, MOMY, GSQRT, J13G, J23G, J33G, MAPF)
Definition: scale_atmos_phy_tb_common.F90:66
scale_precision
module PRECISION
Definition: scale_precision.F90:14
scale_atmos_grid_cartesc_index::jblock
integer, public jblock
block size for cache blocking: y
Definition: scale_atmos_grid_cartesC_index.F90:41
scale_tracer::tracer_advc
logical, dimension(qa_max), public tracer_advc
Definition: scale_tracer.F90:46
scale_const::const_eps
real(rp), public const_eps
small number
Definition: scale_const.F90:35
scale_atmos_phy_tb_common::atmos_phy_tb_calc_tend_momx
subroutine, public atmos_phy_tb_calc_tend_momx(MOMX_t_TB, QFLX_MOMX, GSQRT, J13G, J23G, J33G, MAPF, IIS, IIE, JJS, JJE)
Definition: scale_atmos_phy_tb_common.F90:1679
scale_atmos_grid_cartesc_index::iblock
integer, public iblock
block size for cache blocking: x
Definition: scale_atmos_grid_cartesC_index.F90:40
scale_atmos_hydrometeor
module atmosphere / hydrometeor
Definition: scale_atmos_hydrometeor.F90:12
scale_atmos_phy_tb_common::atmos_phy_tb_calc_tend_momy
subroutine, public atmos_phy_tb_calc_tend_momy(MOMY_t_TB, QFLX_MOMY, GSQRT, J13G, J23G, J33G, MAPF, IIS, IIE, JJS, JJE)
Definition: scale_atmos_phy_tb_common.F90:1741
scale_random
module RANDOM
Definition: scale_random.F90:11
scale_file_history
module file_history
Definition: scale_file_history.F90:15
scale_atmos_grid_cartesc_index::i_xyz
integer, public i_xyz
Definition: scale_atmos_grid_cartesC_index.F90:91
scale_prc
module PROCESS
Definition: scale_prc.F90:11
scale_atmos_grid_cartesc_index::i_xy
integer, public i_xy
Definition: scale_atmos_grid_cartesC_index.F90:99
scale_atmos_grid_cartesc_index::i_uyz
integer, public i_uyz
Definition: scale_atmos_grid_cartesC_index.F90:95
scale_atmos_grid_cartesc_index::i_xvw
integer, public i_xvw
Definition: scale_atmos_grid_cartesC_index.F90:94
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:96
scale_atmos_phy_tb_common::atmos_phy_tb_calc_flux_phi
subroutine, public atmos_phy_tb_calc_flux_phi(qflx_phi, DENS, PHI, Kh, FACT, GSQRT, J13G, J23G, J33G, MAPF, horizontal, implicit, a, b, c, dt, IIS, IIE, JJS, JJE)
Definition: scale_atmos_phy_tb_common.F90:1337
scale_atmos_grid_cartesc_index::i_uyw
integer, public i_uyw
Definition: scale_atmos_grid_cartesC_index.F90:93
scale_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:54
scale_matrix
module MATRIX
Definition: scale_matrix.F90:17
scale_atmos_phy_tb_common::atmos_phy_tb_calc_tend_momz
subroutine, public atmos_phy_tb_calc_tend_momz(MOMZ_t_TB, QFLX_MOMZ, GSQRT, J13G, J23G, J33G, MAPF, IIS, IIE, JJS, JJE)
Definition: scale_atmos_phy_tb_common.F90:1616
scale_atmos_grid_cartesc_index::i_xyw
integer, public i_xyw
Definition: scale_atmos_grid_cartesC_index.F90:92
Generated by   doxygen 1.8.17