SCALE-RM
Functions/Subroutines
scale_atmos_phy_tb_smg Module Reference

module ATMOSPHERE / Physics Turbulence More...

Functions/Subroutines

subroutine, public atmos_phy_tb_smg_config (TYPE_TB, I_TKE_out)
 Config. More...
 
subroutine, public atmos_phy_tb_smg_setup (CDZ, CDX, CDY, CZ)
 Setup. More...
 
subroutine, public atmos_phy_tb_smg (qflx_sgs_momz, qflx_sgs_momx, qflx_sgs_momy, qflx_sgs_rhot, qflx_sgs_rhoq, RHOQ_t, Nu, Ri, Pr, MOMZ, MOMX, MOMY, RHOT, DENS, QTRC, N2, SFLX_MW, SFLX_MU, SFLX_MV, SFLX_SH, SFLX_Q, GSQRT, J13G, J23G, J33G, MAPF, dt)
 
real(rp) function mixlen (dz, dx, dy, z, filter_fact)
 
real(rp) function fact (dz, dx, dy)
 

Detailed Description

module ATMOSPHERE / Physics Turbulence

Description
Sub-grid scale turbulence process Smagolinsky-type
Author
Team SCALE
History
  • 2011-11-29 (S.Iga) [new]
  • 2011-12-11 (H.Yashiro) [mod] integrate to SCALE-LES ver.3
  • 2012-03-23 (H.Yashiro) [mod] Explicit index parameter inclusion
  • 2012-03-27 (H.Yashiro) [mod] reconstruction
  • 2012-07-02 (S.Nishizawa) [mod] reconstruction with Brown et al. (1994)
  • 2012-10-26 (S.Nishizawa) [mod] remove surface flux
  • 2013-06-13 (S.Nishizawa) [mod] change mixing length by Brown et al. (1994) and Scotti et al. (1993)
  • 2014-04-02 (S.Nishizawa) [mod] modified for terrrain-following
  • 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
    • Sullivan et al., 1994: A subgrid-scale model for large-eddy simulation of planetary boundary-layer flows. Bound.-Layer Meteor., 71, 247-276
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_IMPLICIT logical .false.
    ATMOS_PHY_TB_SMG_CONSISTENT_TKE logical .true.
    ATMOS_PHY_TB_SMG_HORIZONTAL logical .false.
    ATMOS_PHY_TB_SMG_BOTTOM logical .false.

History Output
No history output

Function/Subroutine Documentation

◆ atmos_phy_tb_smg_config()

subroutine, public scale_atmos_phy_tb_smg::atmos_phy_tb_smg_config ( character(len=*), intent(in)  TYPE_TB,
integer, intent(out)  I_TKE_out 
)

Config.

Definition at line 107 of file scale_atmos_phy_tb_smg.F90.

References scale_stdio::io_fid_log, scale_stdio::io_l, scale_process::prc_mpistop(), and scale_tracer::tracer_regist().

Referenced by scale_atmos_phy_tb::atmos_phy_tb_config(), and scale_atmos_phy_tb_hybrid::atmos_phy_tb_hybrid_config().

107  use scale_process, only: &
108  prc_mpistop
109  use scale_tracer, only: &
110  tracer_regist
111  implicit none
112 
113  character(len=*), intent(in) :: TYPE_TB
114  integer, intent(out) :: I_TKE_out
115  !---------------------------------------------------------------------------
116 
117  if( io_l ) write(io_fid_log,*)
118  if( io_l ) write(io_fid_log,*) '++++++ Module[Turbulence Tracer] / Categ[ATMOS PHYSICS] / Origin[SCALElib]'
119  if( io_l ) write(io_fid_log,*) '*** Tracers for Smagorinsky-type Eddy Viscocity Model'
120 
121  if ( type_tb /= 'SMAGORINSKY' ) then
122  write(*,*) 'xxx ATMOS_PHY_TB_TYPE is not SMAGORINSKY. Check!'
123  call prc_mpistop
124  endif
125 
126  call tracer_regist( i_tke, & ! [OUT]
127  1, & ! [IN]
128  (/ 'TKE_SMG' /), & ! [IN]
129  (/ 'turbulent kinetic energy (Smagorinsky)' /), & ! [IN]
130  (/ 'm2/s2' /), & ! [IN]
131  advc = (/ .false. /) ) ! [IN], optional
132 
133  i_tke_out = i_tke
134 
135  return
scale_process::prc_mpistop
subroutine, public prc_mpistop
Abort MPI.
Definition: scale_process.F90:267
scale_tracer
module TRACER
Definition: scale_tracer.F90:15
scale_process
module PROCESS
Definition: scale_process.F90:14
scale_tracer::tracer_regist
subroutine, public tracer_regist(QS, NQ, NAME, DESC, UNIT, CV, CP, R, ADVC, MASS)
Regist tracer.
Definition: scale_tracer.F90:64
Here is the call graph for this function:
Here is the caller graph for this function:

◆ atmos_phy_tb_smg_setup()

subroutine, public scale_atmos_phy_tb_smg::atmos_phy_tb_smg_setup ( real(rp), dimension(ka), intent(in)  CDZ,
real(rp), dimension(ia), intent(in)  CDX,
real(rp), dimension(ja), intent(in)  CDY,
real(rp), dimension (ka,ia,ja), intent(in)  CZ 
)

Setup.

Definition at line 142 of file scale_atmos_phy_tb_smg.F90.

References scale_grid_index::ia, scale_grid_index::ie, scale_stdio::io_fid_conf, scale_stdio::io_fid_log, scale_stdio::io_fid_nml, scale_stdio::io_l, scale_stdio::io_nml, scale_grid_index::is, scale_grid_index::ja, scale_grid_index::je, scale_grid_index::js, scale_grid_index::ka, scale_grid_index::ke, scale_grid_index::ks, mixlen(), and scale_process::prc_mpistop().

Referenced by scale_atmos_phy_tb::atmos_phy_tb_config(), and scale_atmos_phy_tb_hybrid::atmos_phy_tb_hybrid_config().

142  use scale_process, only: &
143  prc_mpistop
144  implicit none
145 
146  real(RP), intent(in) :: CDZ(KA)
147  real(RP), intent(in) :: CDX(IA)
148  real(RP), intent(in) :: CDY(JA)
149  real(RP), intent(in) :: CZ (KA,IA,JA)
150 
151  real(RP) :: ATMOS_PHY_TB_SMG_Cs
152  real(RP) :: ATMOS_PHY_TB_SMG_filter_fact = 2.0_rp
153  logical :: ATMOS_PHY_TB_SMG_consistent_tke = .true.
154 
155  namelist / param_atmos_phy_tb_smg / &
156  atmos_phy_tb_smg_cs, &
157  atmos_phy_tb_smg_nu_max, &
158  atmos_phy_tb_smg_filter_fact, &
159  atmos_phy_tb_smg_implicit, &
160  atmos_phy_tb_smg_consistent_tke, &
161  atmos_phy_tb_smg_horizontal, &
162  atmos_phy_tb_smg_bottom
163 
164  integer :: ierr
165  integer :: k, i, j
166  !---------------------------------------------------------------------------
167 
168  if( io_l ) write(io_fid_log,*)
169  if( io_l ) write(io_fid_log,*) '++++++ Module[Turbulence] / Categ[ATMOS PHYSICS] / Origin[SCALElib]'
170  if( io_l ) write(io_fid_log,*) '*** Smagorinsky-type Eddy Viscocity Model'
171 
172  atmos_phy_tb_smg_cs = cs
173 
174  !--- read namelist
175  rewind(io_fid_conf)
176  read(io_fid_conf,nml=param_atmos_phy_tb_smg,iostat=ierr)
177  if( ierr < 0 ) then !--- missing
178  if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
179  elseif( ierr > 0 ) then !--- fatal error
180  write(*,*) 'xxx Not appropriate names in namelist PARAM_ATMOS_PHY_TB_SMG. Check!'
181  call prc_mpistop
182  endif
183  if( io_nml ) write(io_fid_nml,nml=param_atmos_phy_tb_smg)
184 
185  cs = atmos_phy_tb_smg_cs
186 
187  rprn = 1.0_rp / prn
188  rric = 1.0_rp / ric
189  prnovric = ( 1.0_rp - prn ) * rric
190 
191  allocate( nu_fact(ka,ia,ja) )
192 
193 #ifdef DEBUG
194  nu_fact(:,:,:) = undef
195 #endif
196  if ( atmos_phy_tb_smg_horizontal ) then
197  do j = js-1, je+1
198  do i = is-1, ie+1
199  do k = ks, ke
200  nu_fact(k,i,j) = cs**2 * ( cdx(i) * cdy(j) )
201  enddo
202  enddo
203  enddo
204 #ifdef DEBUG
205  i = iundef; j = iundef; k = iundef
206 #endif
207  atmos_phy_tb_smg_consistent_tke = .false.
208  atmos_phy_tb_smg_implicit = .false. ! flux in the z-direction is not necessary
209  else
210  do j = js-1, je+1
211  do i = is-1, ie+1
212  do k = ks, ke
213  nu_fact(k,i,j) = ( cs * mixlen(cdz(k),cdx(i),cdy(j),cz(k,i,j),atmos_phy_tb_smg_filter_fact) )**2
214  enddo
215  enddo
216  enddo
217 #ifdef DEBUG
218  i = iundef; j = iundef; k = iundef
219 #endif
220  end if
221 
222  if ( atmos_phy_tb_smg_consistent_tke ) then
223  tke_fact = 1.0_rp
224  else
225  tke_fact = 0.0_rp
226  end if
227 
228  return
scale_grid_index::is
integer, public is
start point of inner domain: x, local
Definition: scale_grid_index.F90:60
scale_grid_index::je
integer, public je
end point of inner domain: y, local
Definition: scale_grid_index.F90:63
scale_process::prc_mpistop
subroutine, public prc_mpistop
Abort MPI.
Definition: scale_process.F90:267
scale_grid_index::js
integer, public js
start point of inner domain: y, local
Definition: scale_grid_index.F90:62
scale_process
module PROCESS
Definition: scale_process.F90:14
scale_grid_index::ie
integer, public ie
end point of inner domain: x, local
Definition: scale_grid_index.F90:61
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,qa), intent(inout)  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)  RHOT,
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 (ia,ja), intent(in)  SFLX_MW,
real(rp), dimension (ia,ja), intent(in)  SFLX_MU,
real(rp), dimension (ia,ja), intent(in)  SFLX_MV,
real(rp), dimension (ia,ja), intent(in)  SFLX_SH,
real(rp), dimension (ia,ja,qa), intent(in)  SFLX_Q,
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 240 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_calc_tend_phi(), scale_atmos_phy_tb_common::atmos_phy_tb_diffusion_solver(), scale_debug::check(), scale_const::const_eps, scale_const::const_grav, scale_grid::grid_rcdz, scale_grid::grid_rfdz, scale_gridtrans::i_uyw, scale_gridtrans::i_uyz, scale_gridtrans::i_xvw, scale_gridtrans::i_xvz, scale_gridtrans::i_xyw, scale_gridtrans::i_xyz, scale_grid_index::ia, scale_grid_index::iblock, scale_grid_index::ie, scale_stdio::io_fid_log, scale_stdio::io_l, scale_grid_index::is, scale_grid_index::ja, scale_grid_index::jblock, scale_grid_index::je, scale_grid_index::js, scale_grid_index::ke, scale_grid_index::ks, scale_tracer::qa, scale_tracer::tracer_advc, scale_grid_index::xdir, scale_grid_index::ydir, and scale_grid_index::zdir.

Referenced by scale_atmos_phy_tb::atmos_phy_tb_config(), and scale_atmos_phy_tb_hybrid::atmos_phy_tb_hybrid_config().

240  use scale_precision
241  use scale_grid_index
242  use scale_tracer
243  use scale_const, only: &
244  eps => const_eps, &
245  grav => const_grav
246  use scale_grid, only: &
247  rcdz => grid_rcdz, &
248  rfdz => grid_rfdz
249  use scale_gridtrans, only: &
250  i_xyz, &
251  i_xyw, &
252  i_uyw, &
253  i_xvw, &
254  i_uyz, &
255  i_xvz
256  use scale_atmos_phy_tb_common, only: &
257  calc_strain_tensor => atmos_phy_tb_calc_strain_tensor, &
258  diffusion_solver => atmos_phy_tb_diffusion_solver, &
259  calc_tend_momz => atmos_phy_tb_calc_tend_momz, &
260  calc_tend_momx => atmos_phy_tb_calc_tend_momx, &
261  calc_tend_momy => atmos_phy_tb_calc_tend_momy, &
262  calc_tend_phi => atmos_phy_tb_calc_tend_phi, &
263  calc_flux_phi => atmos_phy_tb_calc_flux_phi
264  implicit none
265 
266  ! SGS flux
267  real(RP), intent(out) :: qflx_sgs_momz(KA,IA,JA,3)
268  real(RP), intent(out) :: qflx_sgs_momx(KA,IA,JA,3)
269  real(RP), intent(out) :: qflx_sgs_momy(KA,IA,JA,3)
270  real(RP), intent(out) :: qflx_sgs_rhot(KA,IA,JA,3)
271  real(RP), intent(out) :: qflx_sgs_rhoq(KA,IA,JA,3,QA)
272 
273  real(RP), intent(inout) :: RHOQ_t (KA,IA,JA,QA) ! tendency of rho * QTRC
274 
275  real(RP), intent(out) :: nu (KA,IA,JA) ! eddy viscosity (center)
276  real(RP), intent(out) :: Ri (KA,IA,JA) ! Richardson number
277  real(RP), intent(out) :: Pr (KA,IA,JA) ! Prantle number
278 
279  real(RP), intent(in) :: MOMZ (KA,IA,JA)
280  real(RP), intent(in) :: MOMX (KA,IA,JA)
281  real(RP), intent(in) :: MOMY (KA,IA,JA)
282  real(RP), intent(in) :: RHOT (KA,IA,JA)
283  real(RP), intent(in) :: DENS (KA,IA,JA)
284  real(RP), intent(in) :: QTRC (KA,IA,JA,QA)
285  real(RP), intent(in) :: N2 (KA,IA,JA)
286 
287  real(RP), intent(in) :: SFLX_MW (IA,JA)
288  real(RP), intent(in) :: SFLX_MU (IA,JA)
289  real(RP), intent(in) :: SFLX_MV (IA,JA)
290  real(RP), intent(in) :: SFLX_SH (IA,JA)
291  real(RP), intent(in) :: SFLX_Q (IA,JA,QA)
292 
293  real(RP), intent(in) :: GSQRT (KA,IA,JA,7)
294  real(RP), intent(in) :: J13G (KA,IA,JA,7)
295  real(RP), intent(in) :: J23G (KA,IA,JA,7)
296  real(RP), intent(in) :: J33G
297  real(RP), intent(in) :: MAPF(IA,JA,2,4)
298  real(DP), intent(in) :: dt
299 
300  ! diagnostic variables
301  real(RP) :: TKE (KA,IA,JA)
302  real(RP) :: POTT (KA,IA,JA)
303 
304  ! deformation rate tensor
305  real(RP) :: S33_C(KA,IA,JA) ! (cell center)
306  real(RP) :: S11_C(KA,IA,JA)
307  real(RP) :: S22_C(KA,IA,JA)
308  real(RP) :: S31_C(KA,IA,JA)
309  real(RP) :: S12_C(KA,IA,JA)
310  real(RP) :: S23_C(KA,IA,JA)
311  real(RP) :: S12_Z(KA,IA,JA) ! (z edge or x-y plane)
312  real(RP) :: S23_X(KA,IA,JA) ! (x edge or y-z plane)
313  real(RP) :: S31_Y(KA,IA,JA) ! (y edge or z-x plane)
314  real(RP) :: S2 (KA,IA,JA) ! |S|^2
315 
316  real(RP) :: Kh (KA,IA,JA) ! eddy diffusion
317 
318  real(RP) :: TEND (KA,IA,JA)
319  real(RP) :: a (KA,IA,JA)
320  real(RP) :: b (KA,IA,JA)
321  real(RP) :: c (KA,IA,JA)
322  real(RP) :: d (KA)
323  real(RP) :: ap
324 
325  integer :: IIS, IIE
326  integer :: JJS, JJE
327 
328  integer :: k, i, j, iq
329  !---------------------------------------------------------------------------
330 
331  if( io_l ) write(io_fid_log,*) '*** Atmos physics step: Turbulence(smagorinsky)'
332 
333 #ifdef DEBUG
334  qflx_sgs_momz(:,:,:,:) = undef
335  qflx_sgs_momx(:,:,:,:) = undef
336  qflx_sgs_momy(:,:,:,:) = undef
337  qflx_sgs_rhot(:,:,:,:) = undef
338  qflx_sgs_rhoq(:,:,:,:,:) = undef
339 
340  nu(:,:,:) = undef
341  tke(:,:,:) = undef
342  pr(:,:,:) = undef
343  ri(:,:,:) = undef
344  kh(:,:,:) = undef
345 
346  pott(:,:,:) = undef
347 #endif
348 
349 #ifdef QUICKDEBUG
350  qflx_sgs_momz(ks:ke, 1:is-1, : ,:) = undef
351  qflx_sgs_momz(ks:ke,ie+1:ia , : ,:) = undef
352  qflx_sgs_momz(ks:ke, : , 1:js-1,:) = undef
353  qflx_sgs_momz(ks:ke, : ,je+1:ja ,:) = undef
354  qflx_sgs_momx(ks:ke, 1:is-1, : ,:) = undef
355  qflx_sgs_momx(ks:ke,ie+1:ia , : ,:) = undef
356  qflx_sgs_momx(ks:ke, : , 1:js-1,:) = undef
357  qflx_sgs_momx(ks:ke, : ,je+1:ja ,:) = undef
358  qflx_sgs_momy(ks:ke, 1:is-1, : ,:) = undef
359  qflx_sgs_momy(ks:ke,ie+1:ia , : ,:) = undef
360  qflx_sgs_momy(ks:ke, : , 1:js-1,:) = undef
361  qflx_sgs_momy(ks:ke, : ,je+1:ja ,:) = undef
362 #endif
363 
364  ! potential temperature
365  !$omp parallel do default(none) &
366  !$omp shared(JS,JE,IS,IE,KS,KE,RHOT,DENS,POTT) &
367  !$omp private(i,j,k) OMP_SCHEDULE_ collapse(2)
368  do j = js-1, je+1
369  do i = is-1, ie+1
370  do k = ks, ke
371 #ifdef DEBUG
372  call check( __line__, rhot(k,i,j) )
373  call check( __line__, dens(k,i,j) )
374 #endif
375  pott(k,i,j) = rhot(k,i,j) / dens(k,i,j)
376  enddo
377  enddo
378  enddo
379 #ifdef DEBUG
380  i = iundef; j = iundef; k = iundef
381 #endif
382 
383  !##### Start Upadate #####
384 
385  call calc_strain_tensor( &
386  s33_c, s11_c, s22_c, & ! (out)
387  s31_c, s12_c, s23_c, & ! (out)
388  s12_z, s23_x, s31_y, & ! (out)
389  s2 , & ! (out)
390  dens, momz, momx, momy, & ! (in)
391  gsqrt, j13g, j23g, j33g, mapf ) ! (in)
392 
393  do jjs = js, je, jblock
394  jje = jjs+jblock-1
395  do iis = is, ie, iblock
396  iie = iis+iblock-1
397 
398  ! Ri = N^2 / |S|^2, N^2 = g / theta * dtheta/dz
399  do j = jjs-1, jje+1
400  do i = iis-1, iie+1
401  do k = ks, ke
402 #ifdef DEBUG
403  call check( __line__, s2(k,i,j) )
404  call check( __line__, n2(k,i,j) )
405 #endif
406  ri(k,i,j) = n2(k,i,j) / s2(k,i,j)
407  enddo
408  enddo
409  enddo
410 #ifdef DEBUG
411  i = iundef; j = iundef; k = iundef
412 #endif
413 
414  do j = jjs-1, jje+1
415  do i = iis-1, iie+1
416  do k = ks, ke
417 #ifdef DEBUG
418  call check( __line__, ri(k,i,j) )
419  call check( __line__, nu_fact(k,i,j) )
420  call check( __line__, s2(k,i,j) )
421 #endif
422  if ( ri(k,i,j) < 0.0_rp ) then ! stable
423  nu(k,i,j) = nu_fact(k,i,j) &
424  * sqrt( s2(k,i,j) * (1.0_rp - fmc*ri(k,i,j)) )
425  else if ( ri(k,i,j) < ric ) then ! weakly stable
426  nu(k,i,j) = nu_fact(k,i,j) &
427  * sqrt( s2(k,i,j) ) * ( 1.0_rp - ri(k,i,j)*rric )**4
428  else ! strongly stable
429  nu(k,i,j) = 0.0_rp
430  endif
431  enddo
432  enddo
433  enddo
434 #ifdef DEBUG
435  i = iundef; j = iundef; k = iundef
436 #endif
437  ! Pr = nu_m / nu_h = fm / fh
438  do j = jjs-1, jje+1
439  do i = iis-1, iie+1
440  do k = ks, ke
441 #ifdef DEBUG
442  call check( __line__, ri(k,i,j) )
443 #endif
444  if ( ri(k,i,j) < 0.0_rp ) then ! stable
445  pr(k,i,j) = sqrt( ( 1.0_rp - fmc*ri(k,i,j) ) &
446  / ( 1.0_rp - fhb*ri(k,i,j) ) ) * prn
447  else if ( ri(k,i,j) < ric ) then ! weakly stable
448  pr(k,i,j) = prn / ( 1.0_rp - prnovric * ri(k,i,j) )
449  else ! strongly stable
450  pr(k,i,j) = 1.0_rp
451  endif
452  kh(k,i,j) = min( nu(k,i,j)/pr(k,i,j), atmos_phy_tb_smg_nu_max )
453  nu(k,i,j) = min( nu(k,i,j), atmos_phy_tb_smg_nu_max )
454  pr(k,i,j) = nu(k,i,j) / ( kh(k,i,j) + ( 0.5_rp - sign(0.5_rp, kh(k,i,j)-eps) ) )
455  enddo
456  enddo
457  enddo
458 #ifdef DEBUG
459  i = iundef; j = iundef; k = iundef
460 #endif
461 
462  ! tke = (nu/(Ck * Delta))^2 = ( nu * Cs / Ck )^2 / ( Cs * Delta )^2
463  ! Sullivan et al. (1994)
464  !$omp parallel do default(none) &
465  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,nu,nu_fact,Cs,TKE) &
466  !$omp private(i,j,k) OMP_SCHEDULE_ collapse(2)
467  do j = jjs, jje+1
468  do i = iis, iie+1
469  do k = ks, ke
470 #ifdef DEBUG
471  call check( __line__, nu(k,i,j) )
472  call check( __line__, nu_fact(k,i,j) )
473 #endif
474  tke(k,i,j) = ( nu(k,i,j) * cs / ck )**2 / nu_fact(k,i,j)
475  enddo
476  enddo
477  enddo
478 #ifdef DEBUG
479  i = iundef; j = iundef; k = iundef
480 #endif
481 
482  !##### momentum equation (z) #####
483  ! (cell center)
484  if ( atmos_phy_tb_smg_horizontal ) then
485  qflx_sgs_momz(:,:,:,zdir) = 0.0_rp
486  else
487  do j = jjs, jje
488  do i = iis, iie
489  do k = ks+1, ke-1
490 #ifdef DEBUG
491  call check( __line__, dens(k,i,j) )
492  call check( __line__, nu(k,i,j) )
493  call check( __line__, s33_c(k,i,j) )
494  call check( __line__, s11_c(k,i,j) )
495  call check( __line__, s22_c(k,i,j) )
496  call check( __line__, tke(k,i,j) )
497 #endif
498  qflx_sgs_momz(k,i,j,zdir) = dens(k,i,j) * ( &
499  - 2.0_rp * nu(k,i,j) &
500  * ( s33_c(k,i,j) - ( s11_c(k,i,j) + s22_c(k,i,j) + s33_c(k,i,j) ) * oneoverthree * tke_fact ) &
501  + twooverthree * tke(k,i,j) * tke_fact )
502  enddo
503  enddo
504  enddo
505 #ifdef DEBUG
506  i = iundef; j = iundef; k = iundef
507 #endif
508  do j = jjs, jje
509  do i = iis, iie
510  qflx_sgs_momz(ks,i,j,zdir) = 0.0_rp ! just above bottom boundary
511  qflx_sgs_momz(ke,i,j,zdir) = 0.0_rp ! just below top boundary
512  enddo
513  enddo
514 #ifdef DEBUG
515  i = iundef; j = iundef; k = iundef
516 #endif
517  end if
518  ! (y edge)
519  do j = jjs, jje
520  do i = iis-1, iie
521  do k = ks, ke-1
522 #ifdef DEBUG
523  call check( __line__, dens(k,i,j) )
524  call check( __line__, dens(k,i+1,j) )
525  call check( __line__, dens(k+1,i,j) )
526  call check( __line__, dens(k+1,i+1,j) )
527  call check( __line__, nu(k,i,j) )
528  call check( __line__, nu(k,i+1,j) )
529  call check( __line__, nu(k+1,i,j) )
530  call check( __line__, nu(k+1,i+1,j) )
531  call check( __line__, s31_y(k,i,j) )
532 #endif
533  qflx_sgs_momz(k,i,j,xdir) = - 0.125_rp & ! 2.0 / 4 / 4
534  * ( dens(k,i,j)+dens(k+1,i,j)+dens(k,i+1,j)+dens(k+1,i+1,j) ) &
535  * ( nu(k,i,j)+nu(k+1,i,j)+nu(k,i+1,j)+nu(k+1,i+1,j)) &
536  * s31_y(k,i,j)
537  enddo
538  enddo
539  enddo
540 #ifdef DEBUG
541  i = iundef; j = iundef; k = iundef
542 #endif
543  ! (x edge)
544  !$omp parallel do default(none) &
545  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,DENS,S23_X,nu,qflx_sgs_momz) &
546  !$omp private(i,j,k) OMP_SCHEDULE_ collapse(2)
547  do j = jjs-1, jje
548  do i = iis, iie
549  do k = ks, ke-1
550 #ifdef DEBUG
551  call check( __line__, dens(k,i,j) )
552  call check( __line__, dens(k,i,j+1) )
553  call check( __line__, dens(k+1,i,j) )
554  call check( __line__, dens(k+1,i,j+1) )
555  call check( __line__, nu(k,i,j) )
556  call check( __line__, nu(k,i,j+1) )
557  call check( __line__, nu(k+1,i,j) )
558  call check( __line__, nu(k+1,i,j+1) )
559  call check( __line__, s23_x(k,i,j) )
560 #endif
561  qflx_sgs_momz(k,i,j,ydir) = - 0.125_rp & ! 2/4/4
562  * ( dens(k,i,j)+dens(k+1,i,j)+dens(k,i,j+1)+dens(k+1,i,j+1) ) &
563  * ( nu(k,i,j)+nu(k+1,i,j)+nu(k,i,j+1)+nu(k+1,i,j+1) ) &
564  * s23_x(k,i,j)
565  enddo
566  enddo
567  enddo
568 #ifdef DEBUG
569  i = iundef; j = iundef; k = iundef
570 #endif
571 
572  if ( atmos_phy_tb_smg_implicit ) then
573 
574  call calc_tend_momz( tend, & ! (out)
575  qflx_sgs_momz, & ! (in)
576  gsqrt, j13g, j23g, j33g, mapf, & ! (in)
577  iis, iie, jjs, jje ) ! (in)
578 
579  do j = jjs, jje
580  do i = iis, iie
581 
582  ap = - fouroverthree * dt &
583  * dens(ks+1,i,j)*nu(ks+1,i,j) &
584  * rcdz(ks+1) / gsqrt(ks+1,i,j,i_xyz)
585  a(ks,i,j) = ap * rfdz(ks) / gsqrt(ks,i,j,i_xyw)
586  c(ks,i,j) = 0.0_rp
587  b(ks,i,j) = - a(ks,i,j) + 0.5_rp * ( dens(ks,i,j)+dens(ks+1,i,j) )
588  do k = ks+1, ke-2
589  c(k,i,j) = ap * rfdz(k+1) / gsqrt(k+1,i,j,i_xyw)
590  ap = - fouroverthree * dt &
591  * dens(k+1,i,j)*nu(k+1,i,j) &
592  * rcdz(k+1) / gsqrt(k+1,i,j,i_xyz)
593  a(k,i,j) = ap * rfdz(k) / gsqrt(k,i,j,i_xyw)
594  b(k,i,j) = - a(k,i,j) - c(k,i,j) + 0.5_rp * ( dens(k,i,j)+dens(k+1,i,j) )
595  end do
596  a(ke-1,i,j) = 0.0_rp
597  c(ke-1,i,j) = ap * rfdz(ke) / gsqrt(ke,i,j,i_xyw)
598  b(ke-1,i,j) = - c(ke-1,i,j) + 0.5_rp * ( dens(ke-1,i,j)+dens(ke,i,j) )
599 
600  do k = ks, ke-1
601  d(k) = tend(k,i,j)
602  end do
603 
604  call diffusion_solver( &
605  tend(:,i,j), & ! (out)
606  a(:,i,j), b(:,i,j), c(:,i,j), d, & ! (in)
607  ke-1 ) ! (in)
608 
609  do k = ks+1, ke-1
610  qflx_sgs_momz(k,i,j,zdir) = qflx_sgs_momz(k,i,j,zdir) &
611  - fouroverthree * dens(k,i,j) * nu(k,i,j) * dt &
612  * ( tend(k,i,j) - tend(k-1,i,j) ) * rcdz(k) / gsqrt(k,i,j,i_xyz)
613  end do
614 
615  end do
616  end do
617 
618  end if
619 
620  !##### momentum equation (x) #####
621  ! (y edge)
622  if ( atmos_phy_tb_smg_horizontal ) then
623  qflx_sgs_momx(:,:,:,zdir) = 0.0_rp
624  else
625  do j = jjs, jje
626  do i = iis, iie
627  do k = ks, ke-1
628 #ifdef DEBUG
629  call check( __line__, dens(k,i,j) )
630  call check( __line__, dens(k,i+1,j) )
631  call check( __line__, dens(k+1,i,j) )
632  call check( __line__, dens(k+1,i+1,j) )
633  call check( __line__, nu(k,i,j) )
634  call check( __line__, nu(k,i+1,j) )
635  call check( __line__, nu(k+1,i,j) )
636  call check( __line__, nu(k+1,i+1,j) )
637  call check( __line__, s31_y(k,i,j) )
638 #endif
639  qflx_sgs_momx(k,i,j,zdir) = - 0.125_rp & ! 2/4/4
640  * ( dens(k,i,j)+dens(k+1,i,j)+dens(k,i+1,j)+dens(k+1,i+1,j) ) &
641  * ( nu(k,i,j)+nu(k+1,i,j)+nu(k,i+1,j)+nu(k+1,i+1,j) ) &
642  * s31_y(k,i,j)
643  enddo
644  enddo
645  enddo
646 #ifdef DEBUG
647  i = iundef; j = iundef; k = iundef
648 #endif
649  do j = jjs, jje
650  do i = iis, iie
651  qflx_sgs_momx(ks-1,i,j,zdir) = 0.0_rp ! bottom boundary
652  qflx_sgs_momx(ke ,i,j,zdir) = 0.0_rp ! top boundary
653  enddo
654  enddo
655 #ifdef DEBUG
656  i = iundef; j = iundef; k = iundef
657 #endif
658  end if
659  ! (cell center)
660  !$omp parallel do default(none) &
661  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,DENS,nu,S11_C,S22_C,S33_C,TKE,tke_fact,qflx_sgs_momx) &
662  !$omp private(i,j,k) OMP_SCHEDULE_ collapse(2)
663  do j = jjs, jje
664  do i = iis, iie+1
665  do k = ks, ke
666 #ifdef DEBUG
667  call check( __line__, dens(k,i,j) )
668  call check( __line__, nu(k,i,j) )
669  call check( __line__, s11_c(k,i,j) )
670  call check( __line__, s22_c(k,i,j) )
671  call check( __line__, s33_c(k,i,j) )
672  call check( __line__, tke(k,i,j) )
673 #endif
674  qflx_sgs_momx(k,i,j,xdir) = dens(k,i,j) * ( &
675  - 2.0_rp * nu(k,i,j) &
676  * ( s11_c(k,i,j) - ( s11_c(k,i,j) + s22_c(k,i,j) + s33_c(k,i,j) ) * oneoverthree * tke_fact ) &
677  + twooverthree * tke(k,i,j) * tke_fact )
678  enddo
679  enddo
680  enddo
681 #ifdef DEBUG
682  i = iundef; j = iundef; k = iundef
683 #endif
684  ! (z edge)
685  !$omp parallel do default(none) &
686  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,DENS,nu,S12_Z,qflx_sgs_momx) &
687  !$omp private(i,j,k) OMP_SCHEDULE_ collapse(2)
688  do j = jjs-1, jje
689  do i = iis, iie
690  do k = ks, ke
691 #ifdef DEBUG
692  call check( __line__, dens(k,i,j) )
693  call check( __line__, dens(k,i+1,j) )
694  call check( __line__, dens(k,i,j+1) )
695  call check( __line__, dens(k,i+1,j+1) )
696  call check( __line__, nu(k,i,j) )
697  call check( __line__, nu(k,i+1,j) )
698  call check( __line__, nu(k,i,j+1) )
699  call check( __line__, nu(k,i+1,j+1) )
700  call check( __line__, s12_z(k,i,j) )
701 #endif
702  qflx_sgs_momx(k,i,j,ydir) = - 0.125_rp & ! 2/4/4
703  * ( dens(k,i,j)+dens(k,i+1,j)+dens(k,i,j+1)+dens(k,i+1,j+1) ) &
704  * ( nu(k,i,j)+nu(k,i+1,j)+nu(k,i,j+1)+nu(k,i+1,j+1) ) &
705  * s12_z(k,i,j)
706  enddo
707  enddo
708  enddo
709 #ifdef DEBUG
710  i = iundef; j = iundef; k = iundef
711 #endif
712 
713  if ( atmos_phy_tb_smg_implicit ) then
714  call calc_tend_momx( tend, & ! (out)
715  qflx_sgs_momx, & ! (in)
716  gsqrt, j13g, j23g, j33g, mapf, & ! (in)
717  iis, iie, jjs, jje ) ! (in)
718 
719  do j = jjs, jje
720  do i = iis, iie
721 
722  ap = - dt * 0.25_rp * ( dens(ks ,i ,j)*nu(ks ,i ,j) &
723  + dens(ks+1,i ,j)*nu(ks+1,i ,j) &
724  + dens(ks ,i+1,j)*nu(ks ,i+1,j) &
725  + dens(ks+1,i+1,j)*nu(ks+1,i+1,j) ) &
726  * rfdz(ks) / gsqrt(ks,i,j,i_uyw)
727  a(ks,i,j) = ap * rcdz(ks) / gsqrt(ks,i,j,i_uyz)
728  c(ks,i,j) = 0.0_rp
729  b(ks,i,j) = - a(ks,i,j) + 0.5_rp * ( dens(ks,i,j)+dens(ks,i+1,j) )
730  do k = ks+1, ke-1
731  c(k,i,j) = ap * rcdz(k) / gsqrt(k,i,j,i_uyz)
732  ap = - dt * 0.25_rp * ( dens(k ,i ,j)*nu(k ,i ,j) &
733  + dens(k+1,i ,j)*nu(k+1,i ,j) &
734  + dens(k ,i+1,j)*nu(k ,i+1,j) &
735  + dens(k+1,i+1,j)*nu(k+1,i+1,j) ) &
736  * rfdz(k) / gsqrt(k,i,j,i_uyw)
737  a(k,i,j) = ap * rcdz(k) / gsqrt(k,i,j,i_uyz)
738  b(k,i,j) = - a(k,i,j) - c(k,i,j) + 0.5_rp * ( dens(k,i,j)+dens(k,i+1,j) )
739  end do
740  a(ke,i,j) = 0.0_rp
741  c(ke,i,j) = ap * rcdz(ke) / gsqrt(ke,i,j,i_uyz)
742  b(ke,i,j) = - c(ke,i,j) + 0.5_rp * ( dens(ke,i,j)+dens(ke,i+1,j) )
743 
744  do k = ks, ke
745  d(k) = tend(k,i,j)
746  end do
747 
748  call diffusion_solver( &
749  tend(:,i,j), & ! (out)
750  a(:,i,j), b(:,i,j), c(:,i,j), d, & ! (in)
751  ke ) ! (in)
752 
753  do k = ks, ke-1
754  qflx_sgs_momx(k,i,j,zdir) = qflx_sgs_momx(k,i,j,zdir) &
755  - 0.25_rp * ( dens(k ,i ,j)*nu(k ,i ,j) &
756  + dens(k+1,i ,j)*nu(k+1,i ,j) &
757  + dens(k ,i+1,j)*nu(k ,i+1,j) &
758  + dens(k+1,i+1,j)*nu(k+1,i+1,j) ) &
759  * dt * ( tend(k+1,i,j) - tend(k,i,j) ) * rfdz(k) / gsqrt(k,i,j,i_uyw)
760  end do
761 
762  end do
763  end do
764 
765  end if
766 
767  !##### momentum equation (y) #####
768  ! (x edge)
769  if ( atmos_phy_tb_smg_horizontal ) then
770  qflx_sgs_momy(:,:,:,zdir) = 0.0_rp
771  else
772  do j = jjs, jje
773  do i = iis, iie
774  do k = ks, ke-1
775 #ifdef DEBUG
776  call check( __line__, dens(k,i,j) )
777  call check( __line__, dens(k,i,j+1) )
778  call check( __line__, dens(k+1,i,j) )
779  call check( __line__, dens(k+1,i,j+1) )
780  call check( __line__, nu(k,i,j) )
781  call check( __line__, nu(k,i,j+1) )
782  call check( __line__, nu(k+1,i,j) )
783  call check( __line__, nu(k+1,i,j+1) )
784  call check( __line__, s23_x(k,i,j) )
785 #endif
786  qflx_sgs_momy(k,i,j,zdir) = - 0.125_rp & ! 2/4/4
787  * ( dens(k,i,j)+dens(k+1,i,j)+dens(k,i,j+1)+dens(k+1,i,j+1) ) &
788  * ( nu(k,i,j)+nu(k+1,i,j)+nu(k,i,j+1)+nu(k+1,i,j+1) ) &
789  * s23_x(k,i,j)
790  enddo
791  enddo
792  enddo
793 #ifdef DEBUG
794  i = iundef; j = iundef; k = iundef
795 #endif
796  do j = jjs, jje
797  do i = iis, iie
798  qflx_sgs_momy(ks-1,i,j,zdir) = 0.0_rp ! bottom boundary
799  qflx_sgs_momy(ke ,i,j,zdir) = 0.0_rp ! top boundary
800  enddo
801  enddo
802 #ifdef DEBUG
803  i = iundef; j = iundef; k = iundef
804 #endif
805  end if
806 
807  ! (z edge)
808  !$omp parallel do default(none) &
809  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,DENS,nu,S12_Z,qflx_sgs_momy) &
810  !$omp private(i,j,k) OMP_SCHEDULE_ collapse(2)
811  do j = jjs, jje
812  do i = iis-1, iie
813  do k = ks, ke
814 #ifdef DEBUG
815  call check( __line__, dens(k,i,j) )
816  call check( __line__, dens(k,i+1,j) )
817  call check( __line__, dens(k,i,j+1) )
818  call check( __line__, dens(k,i+1,j+1) )
819  call check( __line__, nu(k,i,j) )
820  call check( __line__, nu(k,i+1,j) )
821  call check( __line__, nu(k,i,j+1) )
822  call check( __line__, nu(k,i+1,j+1) )
823  call check( __line__, s12_z(k,i,j) )
824 #endif
825  qflx_sgs_momy(k,i,j,xdir) = - 0.125_rp & !
826  * ( dens(k,i,j)+dens(k,i+1,j)+dens(k,i,j+1)+dens(k,i+1,j+1) ) &
827  * ( nu(k,i,j)+nu(k,i+1,j)+nu(k,i,j+1)+nu(k,i+1,j+1) ) &
828  * s12_z(k,i,j)
829  enddo
830  enddo
831  enddo
832 #ifdef DEBUG
833  i = iundef; j = iundef; k = iundef
834 #endif
835 
836  ! (z-x plane)
837  !$omp parallel do default(none) &
838  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,DENS,nu,S11_C,S22_C,S33_C,tke_fact,TKE,qflx_sgs_momy) &
839  !$omp private(i,j,k) OMP_SCHEDULE_ collapse(2)
840  do j = jjs, jje+1
841  do i = iis, iie
842  do k = ks, ke
843 #ifdef DEBUG
844  call check( __line__, dens(k,i,j) )
845  call check( __line__, nu(k,i,j) )
846  call check( __line__, s11_c(k,i,j) )
847  call check( __line__, s22_c(k,i,j) )
848  call check( __line__, s33_c(k,i,j) )
849  call check( __line__, tke(k,i,j) )
850 #endif
851  qflx_sgs_momy(k,i,j,ydir) = dens(k,i,j) * ( &
852  - 2.0_rp * nu(k,i,j) &
853  * ( s22_c(k,i,j) - ( s11_c(k,i,j) + s22_c(k,i,j) + s33_c(k,i,j) ) * oneoverthree * tke_fact ) &
854  + twooverthree * tke(k,i,j) * tke_fact)
855  enddo
856  enddo
857  enddo
858 #ifdef DEBUG
859  i = iundef; j = iundef; k = iundef
860 #endif
861 
862  if ( atmos_phy_tb_smg_implicit ) then
863  call calc_tend_momy( tend, & ! (out)
864  qflx_sgs_momy, & ! (in)
865  gsqrt, j13g, j23g, j33g, mapf, & ! (in)
866  iis, iie, jjs, jje ) ! (in)
867 
868  do j = jjs, jje
869  do i = iis, iie
870 
871  ap = - dt * 0.25_rp * ( dens(ks ,i,j )*nu(ks ,i,j ) &
872  + dens(ks+1,i,j )*nu(ks+1,i,j ) &
873  + dens(ks ,i,j+1)*nu(ks ,i,j+1) &
874  + dens(ks+1,i,j+1)*nu(ks+1,i,j+1) ) &
875  * rfdz(ks) / gsqrt(ks,i,j,i_xvw)
876  a(ks,i,j) = ap * rcdz(ks) / gsqrt(ks,i,j,i_xvz)
877  c(ks,i,j) = 0.0_rp
878  b(ks,i,j) = - a(ks,i,j) + 0.5_rp * ( dens(ks,i,j)+dens(ks,i,j+1) )
879  do k = ks+1, ke-1
880  c(k,i,j) = ap * rcdz(k) / gsqrt(k,i,j,i_xvz)
881  ap = - dt * 0.25_rp * ( dens(k ,i,j )*nu(k ,i,j ) &
882  + dens(k+1,i,j )*nu(k+1,i,j ) &
883  + dens(k ,i,j+1)*nu(k ,i,j+1) &
884  + dens(k+1,i,j+1)*nu(k+1,i,j+1) ) &
885  * rfdz(k) / gsqrt(k,i,j,i_xvw)
886  a(k,i,j) = ap * rcdz(k) / gsqrt(k,i,j,i_xvz)
887  b(k,i,j) = - a(k,i,j) - c(k,i,j) + 0.5_rp * ( dens(k,i,j)+dens(k,i,j+1) )
888  end do
889  a(ke,i,j) = 0.0_rp
890  c(ke,i,j) = ap * rcdz(ke) / gsqrt(ke,i,j,i_xvz)
891  b(ke,i,j) = - c(ke,i,j) + 0.5_rp * ( dens(ke,i,j)+dens(ke,i,j+1) )
892 
893  do k = ks, ke
894  d(k) = tend(k,i,j)
895  end do
896 
897  call diffusion_solver( &
898  tend(:,i,j), & ! (out)
899  a(:,i,j), b(:,i,j), c(:,i,j), d, & ! (in)
900  ke ) ! (in)
901 
902  do k = ks, ke-1
903  qflx_sgs_momy(k,i,j,zdir) = qflx_sgs_momy(k,i,j,zdir) &
904  - 0.25_rp * ( dens(k ,i,j )*nu(k ,i,j ) &
905  + dens(k+1,i,j )*nu(k+1,i,j ) &
906  + dens(k ,i,j+1)*nu(k ,i,j+1) &
907  + dens(k+1,i,j+1)*nu(k+1,i,j+1) ) &
908  * dt * ( tend(k+1,i,j) - tend(k,i,j) ) * rfdz(k) / gsqrt(k,i,j,i_xvw)
909  end do
910 
911  end do
912  end do
913 
914  end if
915 
916  !##### Thermodynamic Equation #####
917 
918  if ( atmos_phy_tb_smg_implicit ) then
919 
920  do j = jjs, jje
921  do i = iis, iie
922 
923  ap = - dt * 0.25_rp * ( dens(ks,i,j)+dens(ks+1,i,j) ) &
924  * ( kh(ks,i,j)+kh(ks+1,i,j) ) &
925  * rfdz(ks) / gsqrt(ks,i,j,i_xyw)
926  a(ks,i,j) = ap * rcdz(ks) / gsqrt(ks,i,j,i_xyz)
927  c(ks,i,j) = 0.0_rp
928  b(ks,i,j) = - a(ks,i,j) + dens(ks,i,j)
929  do k = ks+1, ke-1
930  c(k,i,j) = ap * rcdz(k) / gsqrt(k,i,j,i_xyz)
931  ap = - dt * 0.25_rp * ( dens(k,i,j)+dens(k+1,i,j) ) &
932  * ( kh(k,i,j)+kh(k+1,i,j) ) &
933  * rfdz(k) / gsqrt(k,i,j,i_xyw)
934  a(k,i,j) = ap * rcdz(k) / gsqrt(k,i,j,i_xyz)
935  b(k,i,j) = - a(k,i,j) - c(k,i,j) + dens(k,i,j)
936  end do
937  a(ke,i,j) = 0.0_rp
938  c(ke,i,j) = ap * rcdz(ke) / gsqrt(ke,i,j,i_xyz)
939  b(ke,i,j) = - c(ke,i,j) + dens(ke,i,j)
940 
941  end do
942  end do
943 
944  end if
945 
946  call calc_flux_phi( &
947  qflx_sgs_rhot, &
948  dens, pott, kh, 1.0_rp, &
949  gsqrt, j13g, j23g, j33g, mapf, &
950  a, b, c, dt, &
951  atmos_phy_tb_smg_implicit, &
952  iis, iie, jjs, jje )
953 
954  enddo
955  enddo
956 
957 
958  !##### Tracers #####
959  do iq = 1, qa
960 
961  if ( iq == i_tke .or. .not. tracer_advc(iq) ) cycle
962 
963  do jjs = js, je, jblock
964  jje = jjs+jblock-1
965  do iis = is, ie, iblock
966  iie = iis+iblock-1
967 
968  call calc_flux_phi( &
969  qflx_sgs_rhoq(:,:,:,:,iq), &
970  dens, qtrc(:,:,:,iq), kh, 1.0_rp, &
971  gsqrt, j13g, j23g, j33g, mapf, &
972  a, b, c, dt, &
973  atmos_phy_tb_smg_implicit, &
974  iis, iie, jjs, jje )
975 
976 
977  enddo
978  enddo
979 #ifdef DEBUG
980  iis = iundef; iie = iundef; jjs = iundef; jje = iundef
981 #endif
982 
983  enddo ! scalar quantities loop
984 #ifdef DEBUG
985  iq = iundef
986 #endif
987 
988  do j = js, je
989  do i = is, ie
990  do k = ks, ke
991  rhoq_t(k,i,j,i_tke) = ( tke(k,i,j) - qtrc(k,i,j,i_tke) ) * dens(k,i,j) / dt
992  end do
993  end do
994  end do
995 
996  return
scale_grid_index::is
integer, public is
start point of inner domain: x, local
Definition: scale_grid_index.F90:60
scale_gridtrans::i_xvz
integer, public i_xvz
Definition: scale_gridtrans.F90:54
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:1501
scale_grid_index::je
integer, public je
end point of inner domain: y, local
Definition: scale_grid_index.F90:63
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:1459
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, a, b, c, dt, implicit, IIS, IIE, JJS, JJE)
Definition: scale_atmos_phy_tb_common.F90:1204
scale_grid_index::iblock
integer, public iblock
block size for cache blocking: x
Definition: scale_grid_index.F90:47
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:1688
scale_grid_index::zdir
integer, parameter, public zdir
Definition: scale_grid_index.F90:35
scale_grid_index::ydir
integer, parameter, public ydir
Definition: scale_grid_index.F90:37
scale_grid_index::xdir
integer, parameter, public xdir
Definition: scale_grid_index.F90:36
scale_tracer::tracer_advc
logical, dimension(qa_max), public tracer_advc
Definition: scale_tracer.F90:43
scale_grid::grid_rcdz
real(rp), dimension(:), allocatable, public grid_rcdz
reciprocal of center-dz
Definition: scale_grid_cartesian.F90:62
scale_gridtrans::i_xvw
integer, public i_xvw
Definition: scale_gridtrans.F90:52
scale_grid_index
module grid index
Definition: scale_grid_index.F90:14
scale_atmos_phy_tb_common::atmos_phy_tb_calc_tend_phi
subroutine, public atmos_phy_tb_calc_tend_phi(phi_t_TB, QFLX_phi, GSQRT, J13G, J23G, J33G, MAPF, IIS, IIE, JJS, JJE)
Definition: scale_atmos_phy_tb_common.F90:1780
scale_tracer
module TRACER
Definition: scale_tracer.F90:15
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:1595
scale_gridtrans
module GRIDTRANS
Definition: scale_gridtrans.F90:15
scale_grid_index::jblock
integer, public jblock
block size for cache blocking: y
Definition: scale_grid_index.F90:48
scale_gridtrans::i_xyw
integer, public i_xyw
Definition: scale_gridtrans.F90:50
scale_const::const_grav
real(rp), public const_grav
standard acceleration of gravity [m/s2]
Definition: scale_const.F90:48
scale_grid_index::js
integer, public js
start point of inner domain: y, local
Definition: scale_grid_index.F90:62
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:70
scale_gridtrans::i_uyw
integer, public i_uyw
Definition: scale_gridtrans.F90:51
scale_const
module CONSTANT
Definition: scale_const.F90:14
scale_atmos_phy_tb_common
module ATMOSPHERE / Physics Turbulence
Definition: scale_atmos_phy_tb_common.F90:15
scale_grid
module GRID (cartesian)
Definition: scale_grid_cartesian.F90:17
scale_grid_index::ie
integer, public ie
end point of inner domain: x, local
Definition: scale_grid_index.F90:61
scale_const::const_eps
real(rp), public const_eps
small number
Definition: scale_const.F90:36
scale_gridtrans::i_uyz
integer, public i_uyz
Definition: scale_gridtrans.F90:53
scale_precision
module PRECISION
Definition: scale_precision.F90:13
scale_grid::grid_rfdz
real(rp), dimension(:), allocatable, public grid_rfdz
reciprocal of face-dz
Definition: scale_grid_cartesian.F90:72
scale_gridtrans::i_xyz
integer, public i_xyz
Definition: scale_gridtrans.F90:49
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)  z,
real(rp), intent(in)  filter_fact 
)

Definition at line 1001 of file scale_atmos_phy_tb_smg.F90.

References scale_const::const_karman, and fact().

Referenced by atmos_phy_tb_smg_setup().

1001  use scale_const, only: &
1002  karman => const_karman
1003  implicit none
1004  real(RP), intent(in) :: dz
1005  real(RP), intent(in) :: dx
1006  real(RP), intent(in) :: dy
1007  real(RP), intent(in) :: z
1008  real(RP), intent(in) :: filter_fact
1009  real(RP) :: mixlen ! (out)
1010 
1011  real(RP) :: d0
1012 
1013  d0 = fact(dz, dx, dy) * filter_fact * ( dz * dx * dy )**oneoverthree ! Scotti et al. (1993)
1014  if ( atmos_phy_tb_smg_bottom ) then
1015  mixlen = sqrt( 1.0_rp / ( 1.0_rp/d0**2 + 1.0_rp/(karman*z)**2 ) ) ! Brown et al. (1994)
1016  else
1017  mixlen = d0
1018  end if
1019 
1020  return
scale_const::const_karman
real(rp), parameter, public const_karman
von Karman constant
Definition: scale_const.F90:52
scale_const
module CONSTANT
Definition: scale_const.F90:14
Here is the call graph for this function:
Here is the caller graph for this function:

◆ fact()

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

Definition at line 1024 of file scale_atmos_phy_tb_smg.F90.

Referenced by mixlen().

1024  real(RP), intent(in) :: dz
1025  real(RP), intent(in) :: dx
1026  real(RP), intent(in) :: dy
1027  real(RP) :: fact ! (out)
1028 
1029  real(RP), parameter :: oot = -1.0_rp/3.0_rp
1030  real(RP), parameter :: fot = 5.0_rp/3.0_rp
1031  real(RP), parameter :: eot = 11.0_rp/3.0_rp
1032  real(RP), parameter :: tof = -3.0_rp/4.0_rp
1033  real(RP) :: a1, a2, b1, b2, dmax
1034 
1035 
1036  dmax = max(dz, dx, dy)
1037  if ( dz == dmax ) then
1038  a1 = dx / dmax
1039  a2 = dy / dmax
1040  else if ( dx == dmax ) then
1041  a1 = dz / dmax
1042  a2 = dy / dmax
1043  else ! dy == dmax
1044  a1 = dz / dmax
1045  a2 = dx / dmax
1046  end if
1047  b1 = atan( a1/a2 )
1048  b2 = atan( a2/a1 )
1049 
1050  fact = 1.736_rp * (a1*a2)**oot &
1051  * ( 4.0_rp*p1(b1)*a1**oot + 0.222_rp*p2(b1)*a1**fot + 0.077*p3(b1)*a1**eot - 3.0_rp*b1 &
1052  + 4.0_rp*p1(b2)*a2**oot + 0.222_rp*p2(b2)*a2**fot + 0.077*p3(b2)*a2**eot - 3.0_rp*b2 &
1053  )**tof
1054  return
Here is the caller graph for this function:
Generated by   doxygen 1.8.14