SCALE-RM
Functions/Subroutines
scale_atmos_phy_tb_mynn Module Reference

module ATMOSPHERE / Physics Turbulence More...

Functions/Subroutines

subroutine, public atmos_phy_tb_mynn_config (TYPE_TB, I_TKE_out)
 Config. More...
 
subroutine, public atmos_phy_tb_mynn_setup (CDZ, CDX, CDY, CZ)
 Setup. More...
 
subroutine, public atmos_phy_tb_mynn (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_in, SFLX_MW, SFLX_MU, SFLX_MV, SFLX_SH, SFLX_Q, GSQRT, J13G, J23G, J33G, MAPF, dt)
 
subroutine get_length (l, DENS, q, n2, SFLX_MU, SFLX_MV, SFLX_PT, PT0)
 
subroutine get_q2_level2 (q2_2, dudz2, Ri, l)
 

Detailed Description

module ATMOSPHERE / Physics Turbulence

Description
Boundary layer turbulence model Mellor-Yamada Nakanishi-Niino model
Author
Team SCALE
History
  • 2014-08-27 (S.Nishizawa) [new]
  • Reference
    • Nakanishi and Niino, 2009: Development of an improved turbulence closure model for the atmospheric boundary layer. J. Meteorol. Soc. Japan, 87, 895-912
NAMELIST
  • PARAM_ATMOS_PHY_TB_MYNN
    nametypedefault valuecomment
    ATMOS_PHY_TB_MYNN_TKE_INIT logical .false. set tke with that of level 2 at the first time if .true.
    ATMOS_PHY_TB_MYNN_PBL_MAX real(RP) 1.E+10_RP maximum height of the PBL
    ATMOS_PHY_TB_MYNN_N2_MAX real(RP) 1.E+3_RP
    ATMOS_PHY_TB_MYNN_NU_MIN real(RP) 1.E-6_RP
    ATMOS_PHY_TB_MYNN_NU_MAX real(RP) 10000.0_RP
    ATMOS_PHY_TB_MYNN_KH_MIN real(RP) 1.E-6_RP
    ATMOS_PHY_TB_MYNN_KH_MAX real(RP) 10000.0_RP
    ATMOS_PHY_TB_MYNN_LT_MAX real(RP) 700.0_RP ~ 0.23 * 3 km

History Output
namedescriptionunitvariable
L_mix minxing length m l
POTL liquid potential temperature K POTL
POTV virtual potential temperature K POTV
TKE_gen generation of TKE m2/s3 gen
dUdZ2 dudz2 m2/s2 dudz2

Function/Subroutine Documentation

◆ atmos_phy_tb_mynn_config()

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

Config.

Definition at line 109 of file scale_atmos_phy_tb_mynn.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().

109  use scale_process, only: &
110  prc_mpistop
111  use scale_tracer, only: &
112  tracer_regist
113  implicit none
114 
115  character(len=*), intent(in) :: TYPE_TB
116  integer, intent(out) :: I_TKE_out
117  !---------------------------------------------------------------------------
118 
119  if( io_l ) write(io_fid_log,*)
120  if( io_l ) write(io_fid_log,*) '++++++ Module[Turbulence Tracer] / Categ[ATMOS PHYSICS] / Origin[SCALElib]'
121  if( io_l ) write(io_fid_log,*) '*** Tracers for Mellor-Yamada Nakanishi-Niino scheme'
122 
123  if ( type_tb /= 'MYNN' ) then
124  write(*,*) 'xxx ATMOS_PHY_TB_TYPE is not MYNN. Check!'
125  call prc_mpistop
126  endif
127 
128  call tracer_regist( i_tke, & ! [OUT]
129  1, & ! [IN]
130  (/ 'TKE_MYNN' /), & ! [IN]
131  (/ 'turbulent kinetic energy (MYNN)' /), & ! [IN]
132  (/ 'm2/s2' /) ) ! [IN]
133 
134  i_tke_out = i_tke
135 
136  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_mynn_setup()

subroutine, public scale_atmos_phy_tb_mynn::atmos_phy_tb_mynn_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 143 of file scale_atmos_phy_tb_mynn.F90.

References scale_const::const_pi, 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::je, scale_grid_index::js, scale_grid_index::ke, scale_grid_index::ks, 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().

143  use scale_process, only: &
144  prc_mpistop
145  use scale_const, only: &
146  pi => const_pi
147  implicit none
148 
149  real(RP), intent(in) :: CDZ(KA)
150  real(RP), intent(in) :: CDX(IA)
151  real(RP), intent(in) :: CDY(JA)
152  real(RP), intent(in) :: CZ (KA,IA,JA)
153 
154  real(RP) :: ATMOS_PHY_TB_MYNN_PBL_MAX = 1.e+10_rp
155 
156  namelist / param_atmos_phy_tb_mynn / &
157  atmos_phy_tb_mynn_tke_init, &
158  atmos_phy_tb_mynn_pbl_max, &
159  atmos_phy_tb_mynn_n2_max, &
160  atmos_phy_tb_mynn_nu_min, &
161  atmos_phy_tb_mynn_nu_max, &
162  atmos_phy_tb_mynn_kh_min, &
163  atmos_phy_tb_mynn_kh_max, &
164  atmos_phy_tb_mynn_lt_max
165 
166  integer :: ierr
167  integer :: k, i, j
168  !---------------------------------------------------------------------------
169 
170  if( io_l ) write(io_fid_log,*)
171  if( io_l ) write(io_fid_log,*) '++++++ Module[Turbulence] / Categ[ATMOS PHYSICS] / Origin[SCALElib]'
172  if( io_l ) write(io_fid_log,*) '*** Mellor-Yamada Nakanishi-Niino scheme'
173 
174  !--- read namelist
175  rewind(io_fid_conf)
176  read(io_fid_conf,nml=param_atmos_phy_tb_mynn,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_MYNN. Check!'
181  call prc_mpistop
182  endif
183  if( io_nml ) write(io_fid_nml,nml=param_atmos_phy_tb_mynn)
184 
185  do k = ks, ke-1
186  do j = js, je
187  do i = is, ie
188  if ( atmos_phy_tb_mynn_pbl_max >= cz(k,i,j) ) then
189  ke_pbl = k
190  end if
191  end do
192  end do
193  end do
194 
195  a1 = b1 * (1.0_rp - 3.0_rp * g1) / 6.0_rp
196  a2 = 1.0_rp / (3.0_rp * g1 * b1**(1.0_rp/3.0_rp) * prn )
197  c1 = g1 - 1.0_rp / ( 3.0_rp * a1 * b1**(1.0_rp/3.0_rp) )
198  g2 = ( 2.0_rp * a1 * (3.0_rp - 2.0_rp * c2) + b2 * (1.0_rp - c3) ) / b1
199  f1 = b1 * (g1 - c1) + 2.0_rp * a1 * (3.0_rp - 2.0_rp * c2) + 3.0_rp * a2 * (1.0_rp - c2) * (1.0_rp - c5)
200  f2 = b1 * (g1 + g2) - 3.0_rp * a1 * (1.0_rp - c2)
201 
202  rf1 = b1 * (g1 - c1) / f1
203  rf2 = b1 * g1 / f2
204  rfc = g1 / (g1 + g2)
205 
206  af12 = a1 * f1 / ( a2 * f2 )
207 
208  sqrt_2pi = sqrt( 2.0_rp * pi )
209  rsqrt_2pi = 1.0_rp / sqrt_2pi
210  rsqrt_2 = 1.0_rp / sqrt( 2.0_rp )
211 
212  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_const
module CONSTANT
Definition: scale_const.F90:14
scale_grid_index::ie
integer, public ie
end point of inner domain: x, local
Definition: scale_grid_index.F90:61
scale_const::const_pi
real(rp), public const_pi
pi
Definition: scale_const.F90:34
Here is the call graph for this function:
Here is the caller graph for this function:

◆ atmos_phy_tb_mynn()

subroutine, public scale_atmos_phy_tb_mynn::atmos_phy_tb_mynn ( 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_in,
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 224 of file scale_atmos_phy_tb_mynn.F90.

References scale_atmos_phy_tb_common::atmos_phy_tb_diffusion_solver(), scale_const::const_cpdry, scale_const::const_cpvap, scale_const::const_epstvap, scale_const::const_grav, scale_const::const_lhv0, scale_const::const_rdry, scale_const::const_rvap, get_length(), get_q2_level2(), scale_grid::grid_rcdz, scale_grid::grid_rfdz, scale_atmos_hydrometeor::i_qc, scale_atmos_hydrometeor::i_qi, scale_atmos_hydrometeor::i_qv, scale_gridtrans::i_uyw, scale_gridtrans::i_xvw, 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::ka, scale_grid_index::ke, scale_grid_index::ks, scale_tracer::qa, scale_grid_real::real_cz, scale_tracer::tracer_advc, scale_tracer::tracer_cv, scale_tracer::tracer_mass, scale_tracer::tracer_r, 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().

224  use scale_precision
225  use scale_grid_index
226  use scale_tracer
227  use scale_const, only: &
228  grav => const_grav, &
229  r => const_rdry, &
230  rvap => const_rvap, &
231  cpdry => const_cpdry, &
232  cpvap => const_cpvap, &
233  lhv0 => const_lhv0, &
234  epstvap => const_epstvap
235  use scale_grid, only: &
236  rcdz => grid_rcdz, &
237  rfdz => grid_rfdz
238  use scale_grid_real, only: &
239  cz => real_cz
240  use scale_gridtrans, only: &
241  i_xyz, &
242  i_xyw, &
243  i_xvw, &
244  i_uyw
245  use scale_comm, only: &
246  comm_vars, &
247  comm_wait
248  use scale_atmos_phy_tb_common, only: &
249  diffusion_solver => atmos_phy_tb_diffusion_solver
250  use scale_atmos_thermodyn, only: &
251  atmos_thermodyn_temp_pres
252  use scale_atmos_hydrometeor, only: &
253  hydrometeor_lhv => atmos_hydrometeor_lhv, &
254  hydrometeor_lhs => atmos_hydrometeor_lhs, &
255  i_qv, &
256  i_qc, &
257  i_qi
258  use scale_atmos_saturation, only: &
259  atmos_saturation_dens2qsat => atmos_saturation_dens2qsat_all
260 #ifdef MORE_HIST
261  use scale_history, only: &
262  hist_in
263 #endif
264  implicit none
265 
266  real(RP), intent(out) :: qflx_sgs_momz(KA,IA,JA,3)
267  real(RP), intent(out) :: qflx_sgs_momx(KA,IA,JA,3)
268  real(RP), intent(out) :: qflx_sgs_momy(KA,IA,JA,3)
269  real(RP), intent(out) :: qflx_sgs_rhot(KA,IA,JA,3)
270  real(RP), intent(out) :: qflx_sgs_rhoq(KA,IA,JA,3,QA)
271 
272  real(RP), intent(inout) :: RHOQ_t (KA,IA,JA,QA) ! tendency of rho * QTRC
273 
274  real(RP), intent(out) :: Nu (KA,IA,JA) ! eddy viscosity (center)
275  real(RP), intent(out) :: Ri (KA,IA,JA) ! Richardson number
276  real(RP), intent(out) :: Pr (KA,IA,JA) ! Plandtle number
277 
278  real(RP), intent(in) :: MOMZ (KA,IA,JA)
279  real(RP), intent(in) :: MOMX (KA,IA,JA)
280  real(RP), intent(in) :: MOMY (KA,IA,JA)
281  real(RP), intent(in) :: RHOT (KA,IA,JA)
282  real(RP), intent(in) :: DENS (KA,IA,JA)
283  real(RP), intent(in) :: QTRC (KA,IA,JA,QA)
284  real(RP), intent(in) :: N2_in (KA,IA,JA)
285 
286  real(RP), intent(in) :: SFLX_MW (IA,JA)
287  real(RP), intent(in) :: SFLX_MU (IA,JA)
288  real(RP), intent(in) :: SFLX_MV (IA,JA)
289  real(RP), intent(in) :: SFLX_SH (IA,JA)
290  real(RP), intent(in) :: SFLX_Q (IA,JA,QA)
291 
292  real(RP), intent(in) :: GSQRT (KA,IA,JA,7)
293  real(RP), intent(in) :: J13G (KA,IA,JA,7)
294  real(RP), intent(in) :: J23G (KA,IA,JA,7)
295  real(RP), intent(in) :: J33G
296  real(RP), intent(in) :: MAPF (IA,JA,2,4)
297  real(DP), intent(in) :: dt
298 
299  real(RP) :: U (KA,IA,JA)
300  real(RP) :: V (KA,IA,JA)
301  real(RP) :: N2 (KA,IA,JA) ! squared Brunt-Vaisala frequency
302  real(RP) :: phiN (KA,IA,JA)
303 
304  real(RP) :: sm (KA,IA,JA)
305  real(RP) :: sh (KA,IA,JA)
306  real(RP) :: l (KA,IA,JA)
307  real(RP) :: q (KA,IA,JA)
308  real(RP) :: dudz2(KA,IA,JA)
309  real(RP) :: q2_2 (KA,IA,JA)
310  real(RP) :: Kh
311  real(RP) :: RHOKh(KA,IA,JA)
312 
313  real(RP) :: SFLX_PT(IA,JA) ! surface potential temperature flux
314 
315  real(RP) :: a(KA,IA,JA)
316  real(RP) :: b(KA,IA,JA)
317  real(RP) :: c(KA,IA,JA)
318  real(RP) :: d(KA,IA,JA)
319  real(RP) :: ap
320  real(RP) :: tke_N(KA,IA,JA)
321  real(RP) :: tke
322 
323  real(RP) :: POTT(KA,IA,JA)
324  real(RP) :: POTV(KA,IA,JA)
325  real(RP) :: POTL(KA,IA,JA)
326  real(RP) :: TEML(KA,IA,JA)
327 
328  real(RP) :: Qw(KA,IA,JA)
329  real(RP) :: qdry(KA,IA,JA)
330  real(RP) :: qv
331  real(RP) :: ql
332  real(RP) :: qs
333 
334  real(RP) :: temp(KA,IA,JA)
335  real(RP) :: pres(KA,IA,JA)
336 
337  real(RP) :: LHV (KA,IA,JA)
338  real(RP) :: LHS (KA,IA,JA)
339 
340  real(RP) :: CPtot
341 
342  real(RP) :: ac
343 
344  real(RP) :: Q1
345  real(RP) :: Qsl(KA)
346  real(RP) :: LHVL(KA)
347  real(RP) :: dQsl
348  real(RP) :: sigma_s
349  real(RP) :: RR
350  real(RP) :: Rt
351  real(RP) :: betat
352  real(RP) :: betaq
353  real(RP) :: aa
354  real(RP) :: bb
355  real(RP) :: cc
356 
357 #ifdef MORE_HIST
358  real(RP) :: gen(KA,IA,JA)
359 #endif
360 
361  integer :: IIS, IIE
362  integer :: JJS, JJE
363 
364  integer :: k, i, j, iq
365  !---------------------------------------------------------------------------
366 
367  if( io_l ) write(io_fid_log, *) "*** Atmos physics step: Turbulence (MYNN)"
368 
369 #ifdef DEBUG
370  pott(:,:,:) = undef
371  potl(:,:,:) = undef
372  temp(:,:,:) = undef
373 #endif
374 
375 #if defined DEBUG || defined QUICKDEBUG
376  qflx_sgs_momz(ks:ke, 1:is-1, : ,:) = undef
377  qflx_sgs_momz(ks:ke,ie+1:ia , : ,:) = undef
378  qflx_sgs_momz(ks:ke, : , 1:js-1,:) = undef
379  qflx_sgs_momz(ks:ke, : ,je+1:ja ,:) = undef
380  qflx_sgs_momx(ks:ke, 1:is-1, : ,:) = undef
381  qflx_sgs_momx(ks:ke,ie+1:ia , : ,:) = undef
382  qflx_sgs_momx(ks:ke, : , 1:js-1,:) = undef
383  qflx_sgs_momx(ks:ke, : ,je+1:ja ,:) = undef
384  qflx_sgs_momy(ks:ke, 1:is-1, : ,:) = undef
385  qflx_sgs_momy(ks:ke,ie+1:ia , : ,:) = undef
386  qflx_sgs_momy(ks:ke, : , 1:js-1,:) = undef
387  qflx_sgs_momy(ks:ke, : ,je+1:ja ,:) = undef
388 #endif
389 
390 !OCL XFILL
391  do j = js , je
392  do i = is-1, ie
393  do k = ks-1, ke
394  qflx_sgs_momz(k,i,j,xdir) = 0.0_rp
395  end do
396  end do
397  end do
398 !OCL XFILL
399  do j = js , je
400  do i = is , ie+1
401  do k = ks-1, ke+1
402  qflx_sgs_momx(k,i,j,xdir) = 0.0_rp
403  end do
404  end do
405  end do
406 
407 !OCL XFILL
408  do j = js , je
409  do i = is-1, ie
410  do k = ks-1, ke+1
411  qflx_sgs_momy(k,i,j,xdir) = 0.0_rp
412  end do
413  end do
414  end do
415 !OCL XFILL
416  do j = js , je
417  do i = is-1, ie
418  do k = ks-1, ke+1
419  qflx_sgs_rhot(k,i,j,xdir) = 0.0_rp
420  end do
421  end do
422  end do
423 !OCL XFILL
424  do iq = 1, qa
425  do j = js , je
426  do i = is-1, ie
427  do k = ks-1, ke+1
428  qflx_sgs_rhoq(k,i,j,xdir,iq) = 0.0_rp
429  end do
430  end do
431  end do
432  end do
433 
434 !OCL XFILL
435  do j = js-1, je
436  do i = is , ie
437  do k = ks-1, ke
438  qflx_sgs_momz(k,i,j,ydir) = 0.0_rp
439  end do
440  end do
441  end do
442 !OCL XFILL
443  do j = js-1, je
444  do i = is , ie
445  do k = ks-1, ke+1
446  qflx_sgs_momx(k,i,j,ydir) = 0.0_rp
447  end do
448  end do
449  end do
450 !OCL XFILL
451  do j = js , je+1
452  do i = is , ie
453  do k = ks-1, ke+1
454  qflx_sgs_momy(k,i,j,ydir) = 0.0_rp
455  end do
456  end do
457  end do
458 !OCL XFILL
459  do j = js-1, je
460  do i = is , ie
461  do k = ks-1, ke+1
462  qflx_sgs_rhot(k,i,j,ydir) = 0.0_rp
463  end do
464  end do
465  end do
466 !OCL XFILL
467  do iq = 1, qa
468  do j = js-1, je
469  do i = is , ie
470  do k = ks-1, ke+1
471  qflx_sgs_rhoq(k,i,j,ydir,iq) = 0.0_rp
472  end do
473  end do
474  end do
475  end do
476 !OCL XFILL
477  do j = js, je
478  do i = is, ie
479  do k = ks-1, ke+1
480  qflx_sgs_momz(k,i,j,zdir) = 0.0_rp
481  end do
482  end do
483  end do
484 
485 
486 
487  do jjs = js, je, jblock
488  jje = jjs+jblock-1
489  do iis = is, ie, iblock
490  iie = iis+iblock-1
491 
492 !OCL XFILL
493  do j = jjs , jje+1
494  do i = iis-1, iie+1
495  do k = ks, ke_pbl+1
496  u(k,i,j) = 0.5_rp * ( momx(k,i,j) + momx(k,i-1,j) ) / dens(k,i,j)
497  end do
498  end do
499  end do
500 
501 !OCL XFILL
502  do j = jjs-1, jje+1
503  do i = iis , iie+1
504  do k = ks, ke_pbl+1
505  v(k,i,j) = 0.5_rp * ( momy(k,i,j) + momy(k,i,j-1) ) / dens(k,i,j)
506  end do
507  end do
508  end do
509 
510  end do
511  end do
512 
513  call atmos_thermodyn_temp_pres( temp, pres, & ! (out)
514  dens, rhot, qtrc, & ! (in)
515  tracer_cv, tracer_r, tracer_mass ) ! (in)
516 
517 
518  call hydrometeor_lhv( lhv(:,:,:), temp(:,:,:) )
519  call hydrometeor_lhs( lhs(:,:,:), temp(:,:,:) )
520 
521 !OCL LOOP_NOFUSION,PREFETCH_SEQUENTIAL(SOFT),SWP
522  !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
523  !$omp shared(JS,JE,IS,IE,KS,KE_PBL,QA,I_QV,I_QC,I_QI) &
524  !$omp shared(GRAV,CPdry,EPSTvap,ATMOS_PHY_TB_MYNN_N2_MAX,CZ) &
525  !$omp shared(QTRC,TRACER_MASS,RHOT,DENS,N2_in,SFLX_PT,SFLX_SH) &
526  !$omp shared(Qdry,Qw,LHV,LHS,POTT,POTL,temp,TEML,POTV,n2,U,V,Ri,dudz2) &
527  !$omp private(i,j,k,iq) &
528  !$omp private(qv,ql,qs,CPtot)
529  do j = js, je
530  do i = is, ie
531 
532  do k = ks, ke_pbl+1
533 
534  qv = 0.0_rp
535  if ( i_qv > 0 ) qv = qtrc(k,i,j,i_qv)
536 
537  ql = 0.0_rp
538  if ( i_qc > 0 ) ql = qtrc(k,i,j,i_qc)
539 ! do iq = QWS, QWE
540 ! ql = ql + QTRC(k,i,j,iq)
541 ! end do
542  qs = 0.0_rp
543  if ( i_qi > 0 ) qs = qtrc(k,i,j,i_qi)
544 ! do iq = QIS, QIE
545 ! qs = qs + QTRC(k,i,j,iq)
546 ! end do
547  calc_qdry(qdry(k,i,j), qtrc, tracer_mass, k, i, j, iq)
548 
549  qw(k,i,j) = qv + ql + qs
550 
551  cptot = cpdry * qdry(k,i,j) + cpvap * qv
552 
553  pott(k,i,j) = rhot(k,i,j) / dens(k,i,j)
554  ! liquid water potential temperature
555  potl(k,i,j) = pott(k,i,j) * (1.0_rp - 1.0_rp * (lhv(k,i,j) * ql + lhs(k,i,j) * qs) / ( temp(k,i,j) * cptot ) )
556  teml(k,i,j) = potl(k,i,j) * temp(k,i,j) / pott(k,i,j)
557 
558  ! virtual potential temperature for derivertive
559 ! POTV(k,i,j) = ( 1.0_RP + EPSTvap * Qw(k,i,j) ) * POTL(k,i,j)
560  potv(k,i,j) = ( qdry(k,i,j) + (epstvap+1.0_rp) * qv ) * pott(k,i,j)
561 
562  end do
563 
564  sflx_pt(i,j) = sflx_sh(i,j) / ( cpdry * dens(ks,i,j) ) &
565  * pott(ks,i,j) / temp(ks,i,j)
566 
567  do k = ks, ke_pbl
568  n2(k,i,j) = min( atmos_phy_tb_mynn_n2_max, n2_in(k,i,j) )
569  end do
570 
571  dudz2(ks,i,j) = ( ( u(ks+1,i,j) - u(ks,i,j) )**2 + ( v(ks+1,i,j) - v(ks,i,j) )**2 ) &
572  / ( cz(ks+1,i,j) - cz(ks,i,j) )**2
573  do k = ks+1, ke_pbl
574  dudz2(k,i,j) = ( ( u(k+1,i,j) - u(k-1,i,j) )**2 + ( v(k+1,i,j) - v(k-1,i,j) )**2 ) &
575  / ( cz(k+1,i,j) - cz(k-1,i,j) )**2
576  end do
577 
578  do k = ks, ke_pbl
579  ri(k,i,j) = n2(k,i,j) / max(dudz2(k,i,j), 1e-10_rp)
580  end do
581 
582  end do
583  end do
584 
585  if ( atmos_phy_tb_mynn_tke_init ) then
586 !OCL XFILL
587  do j = js, je
588  do i = is, ie
589  do k = ks, ke_pbl
590  q(k,i,j) = sqrt( tke_min*2.0_rp )
591  end do
592  end do
593  end do
594  else
595 !OCL XFILL
596  do j = js, je
597  do i = is, ie
598  do k = ks, ke_pbl
599  q(k,i,j) = sqrt( max(qtrc(k,i,j,i_tke), tke_min)*2.0_rp )
600  end do
601  end do
602  end do
603  end if
604 
605  ! length
606  call get_length( &
607  l, & ! (out)
608  dens, & ! (in)
609  q, n2, & ! (in)
610  sflx_mu, sflx_mv, sflx_pt, & ! (in)
611  pott ) ! (in)
612 
613  call get_q2_level2( &
614  q2_2, & ! (out)
615  dudz2, ri, l ) ! (in)
616 
617  if ( atmos_phy_tb_mynn_tke_init ) then
618 !OCL XFILL
619  do j = js, je
620  do i = is, ie
621  do k = ks, ke_pbl
622  tke = max(q2_2(k,i,j) * 0.5_rp, tke_min)
623  q(k,i,j) = sqrt( tke * 2.0_rp )
624  end do
625  end do
626  end do
627 
628  atmos_phy_tb_mynn_tke_init = .false.
629  end if
630 
631  call get_smsh( sm, sh, & ! (out)
632  q, q2_2, & ! (in)
633  l, n2, dudz2 ) ! (in)
634 
635 #if defined DEBUG || defined QUICKDEBUG
636  do j = 1, js-1
637  do i = 1, ia
638  do k = ks, ke
639  teml(k,i,j) = 300.0_rp
640  end do
641  end do
642  end do
643  do j = je+1, ja
644  do i = 1, ia
645  do k = ks, ke
646  teml(k,i,j) = 300.0_rp
647  end do
648  end do
649  end do
650  do j = 1, ja
651  do i = 1, is-1
652  do k = ks, ke
653  teml(k,i,j) = 300.0_rp
654  end do
655  end do
656  do i = ie+1, ia
657  do k = ks, ke
658  teml(k,i,j) = 300.0_rp
659  end do
660  end do
661  end do
662 #endif
663 
664 !OCL LOOP_NOFUSION,PREFETCH_SEQUENTIAL(SOFT),SWP
665  !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
666  !$omp private(i,j,k) &
667  !$omp private(CPtot,Qsl,dQsl,aa,bb,ac,sigma_s,Q1,RR,Ql,cc,Rt,betat,betaq,LHVL) &
668  !$omp shared(JS,JE,IS,IE,KS,KE,KE_PBL) &
669  !$omp shared(CZ,rsqrt_2,rsqrt_2pi,sqrt_2pi,GRAV,CPdry,EPSTvap,ATMOS_PHY_TB_MYNN_N2_MAX) &
670  !$omp shared(DENS,l,Qdry,TEMP,TEML,POTT,POTL,POTV,q2_2,Qw,q,dudz2,Ri,n2,sh)
671  do j = js, je
672  do i = is, ie
673 
674  call atmos_saturation_dens2qsat( qsl(:), & ! (out)
675  teml(:,i,j), dens(:,i,j) ) ! (in)
676 
677  call hydrometeor_lhv( lhvl(:), teml(:,i,j) )
678 
679  do k = ks+1, ke_pbl
680 
681  dqsl = ( qsl(k) * lhvl(k) / ( rvap * teml(k,i,j) ) - qsl(k) ) / teml(k,i,j)
682  cptot = qdry(k,i,j) * cpdry + qsl(k) * cpvap
683  aa = 1.0_rp / ( 1.0_rp + lhvl(k)/cptot * dqsl )
684  bb = temp(k,i,j) / pott(k,i,j) * dqsl
685  ac = min( q(k,i,j)/sqrt(q2_2(k,i,j)), 1.0_rp )
686  sigma_s = max( sqrt( 0.25_rp * aa**2 * l(k,i,j)**2 * ac * b2 * sh(k,i,j) ) &
687  * abs( qw(k+1,i,j) - qw(k-1,i,j) - bb * ( potl(k+1,i,j)-potl(k-1,i,j) ) ) &
688  / ( cz(k+1,i,j) - cz(k-1,i,j) ), &
689  1e-10_rp )
690  q1 = aa * ( qw(k,i,j) - qsl(k) ) * 0.5_rp / sigma_s
691  rr = min( max( 0.5_rp * ( 1.0_rp + erf(q1*rsqrt_2) ), 0.0_rp ), 1.0_rp )
692 #if defined(PGI) || defined(SX)
693  ql = min( max( 2.0_rp * sigma_s * ( rr * q1 + rsqrt_2pi * exp( -min( 0.5_rp*q1**2, 1.e+3_rp ) ) ), & ! apply exp limiter
694 #else
695  ql = min( max( 2.0_rp * sigma_s * ( rr * q1 + rsqrt_2pi * exp(-0.5_rp*q1**2) ), &
696 #endif
697  0.0_rp ), &
698  qw(k,i,j) * 0.5_rp )
699  cc = ( 1.0_rp + epstvap * qw(k,i,j) - (1.0_rp+epstvap) * ql ) * pott(k,i,j)/temp(k,i,j) * lhvl(k) / cptot &
700  - (1.0_rp+epstvap) * pott(k,i,j)
701 #if defined(PGI) || defined(SX)
702  rt = min( max( rr - ql / (2.0_rp*sigma_s*sqrt_2pi) * exp( -min( 0.5_rp*q1**2, 1.e+3_rp ) ), 0.0_rp ), 1.0_rp ) ! apply exp limiter
703 #else
704  rt = min( max( rr - ql / (2.0_rp*sigma_s*sqrt_2pi) * exp(-q1**2 * 0.5_rp), 0.0_rp ), 1.0_rp )
705 #endif
706  betat = 1.0_rp + epstvap * qw(k,i,j) - (1.0_rp+epstvap) * ql - rt * aa * bb * cc
707  betaq = epstvap * pott(k,i,j) + rt * aa * cc
708  n2(k,i,j) = min(atmos_phy_tb_mynn_n2_max, &
709  grav * ( ( potl(k+1,i,j) - potl(k-1,i,j) ) * betat &
710  + ( qw(k+1,i,j) - qw(k-1,i,j) ) * betaq ) &
711  / ( cz(k+1,i,j) - cz(k-1,i,j) ) / potv(k,i,j) )
712  end do
713  n2(ks,i,j) = n2(ks+1,i,j)
714 
715  do k = ks, ke_pbl
716  ri(k,i,j) = n2(k,i,j) / max(dudz2(k,i,j), 1e-10_rp)
717  end do
718  do k = ke_pbl+1, ke
719  ri(k,i,j) = 0.0_rp
720  n2(k,i,j) = 0.0_rp
721  end do
722 
723  end do
724  end do
725 
726  ! length
727  call get_length( &
728  l, & ! (out)
729  dens, & ! (in)
730  q, n2, & ! (in)
731  sflx_mu, sflx_mv, sflx_pt, & ! (in)
732  pott ) ! (in)
733 
734  call get_q2_level2( &
735  q2_2, & ! (out)
736  dudz2, ri, l ) ! (in)
737 
738  call get_smsh( &
739  sm, sh, & ! (out)
740  q, q2_2, & ! (in)
741  l, n2, dudz2 ) ! (in)
742 
743 
744  !$omp parallel do default(none) private(i,j,k,Kh) OMP_SCHEDULE_ collapse(2) &
745  !$omp shared(JS,JE,IS,IE,Nu,l,KE_PBL,q,sm,ATMOS_PHY_TB_MYNN_NU_MAX,ATMOS_PHY_TB_MYNN_NU_MIN) &
746  !$omp shared(KA,sh,ATMOS_PHY_TB_MYNN_KH_MAX,ATMOS_PHY_TB_MYNN_KH_MIN,RHOKh,DENS,Pr,KE,KS)
747  do j = js, je
748  do i = is, ie
749  do k = 1, ks-1
750  nu(k,i,j) = 0.0_rp
751  end do
752  do k = ks, ke_pbl
753  nu(k,i,j) = max( min( l(k,i,j) * q(k,i,j) * sm(k,i,j), &
754  atmos_phy_tb_mynn_nu_max ), &
755  atmos_phy_tb_mynn_nu_min )
756  end do
757  do k = ke_pbl+1, ka
758  nu(k,i,j) = 0.0_rp
759  end do
760 
761  do k = ks, ke_pbl
762  kh = max( min( l(k,i,j) * q(k,i,j) * sh(k,i,j), &
763  atmos_phy_tb_mynn_kh_max ), &
764  atmos_phy_tb_mynn_kh_min )
765  rhokh(k,i,j) = dens(k,i,j) * kh
766  pr(k,i,j) = nu(k,i,j) / kh
767  end do
768  do k = ke_pbl+1, ke
769  pr(k,i,j) = 1.0_rp
770  end do
771  end do
772  end do
773 
774  call comm_vars( nu, 1 )
775 
776  ! time integration
777 
778  ! for velocities
779 !OCL INDEPENDENT
780  !$omp parallel do default(none) &
781  !$omp shared(JS,JE,IS,IE,KS,KE_PBL,dt,DENS,Nu,RFDZ,GSQRT,I_XYW,a,c,b,d,RCDZ,I_XYZ,U,SFLX_MU) &
782  !$omp shared(phiN) &
783  !$omp private(i,j,k,ap) OMP_SCHEDULE_ collapse(2)
784  do j = js, je
785  do i = is, ie
786 
787  ap = - dt * 0.5_rp * ( dens(ks ,i,j)*nu(ks ,i,j) &
788  + dens(ks+1,i,j)*nu(ks+1,i,j) ) &
789  * rfdz(ks) / gsqrt(ks,i,j,i_xyw)
790  a(ks,i,j) = ap * rcdz(ks) / ( dens(ks,i,j) * gsqrt(ks,i,j,i_xyz) )
791  c(ks,i,j) = 0.0_rp
792  b(ks,i,j) = - a(ks,i,j) + 1.0_rp
793  do k = ks+1, ke_pbl-1
794  c(k,i,j) = ap * rcdz(k) / ( dens(k,i,j) * gsqrt(k,i,j,i_xyz) )
795  ap = - dt * 0.5_rp * ( dens(k ,i,j)*nu(k ,i,j) &
796  + dens(k+1,i,j)*nu(k+1,i,j) ) &
797  * rfdz(k) / gsqrt(k,i,j,i_xyw)
798  a(k,i,j) = ap * rcdz(k) / ( dens(k,i,j) * gsqrt(k,i,j,i_xyz) )
799  b(k,i,j) = - a(k,i,j) - c(k,i,j) + 1.0_rp
800  end do
801  a(ke_pbl,i,j) = 0.0_rp
802  c(ke_pbl,i,j) = ap * rcdz(ke_pbl) / ( dens(ke_pbl,i,j) * gsqrt(ke_pbl,i,j,i_xyz) )
803  b(ke_pbl,i,j) = - c(ke_pbl,i,j) + 1.0_rp
804 
805  d(ks,i,j) = u(ks,i,j) + dt * sflx_mu(i,j) * rcdz(ks) / ( dens(ks,i,j) * gsqrt(ks,i,j,i_xyz) )
806  do k = ks+1, ke_pbl
807  d(k,i,j) = u(k,i,j)
808  end do
809  call diffusion_solver( &
810  phin(:,i,j), & ! (out)
811  a(:,i,j), b(:,i,j), c(:,i,j), d(:,i,j), & ! (in)
812  ke_pbl ) ! (in)
813  end do
814  end do
815 
816  call comm_vars( phin, 2 )
817  call comm_wait( nu, 1 )
818  call comm_wait( phin, 2 )
819 
820 !OCL INDEPENDENT
821  do jjs = js, je, jblock
822  jje = jjs+jblock-1
823  do iis = is, ie, iblock
824  iie = iis+iblock-1
825 
826  do j = jjs, jje
827  do i = iis, iie
828  do k = ks, ke_pbl-1
829  qflx_sgs_momx(k,i,j,zdir) = - 0.03125_rp & ! 1/4/4/2
830  * ( dens(k,i,j) + dens(k+1,i,j) + dens(k,i+1,j) + dens(k+1,i+1,j) ) &
831  * ( nu(k,i,j) + nu(k+1,i,j) + nu(k,i+1,j) + nu(k+1,i+1,j) ) &
832  * ( (phin(k+1,i,j)+phin(k+1,i+1,j)) - (phin(k,i,j)+phin(k,i+1,j)) ) &
833  * j33g * rfdz(k) / gsqrt(k,i,j,i_uyw)
834  end do
835  end do
836  end do
837  do j = jjs, jje
838  do i = iis, iie
839  qflx_sgs_momx(ks-1,i,j,zdir) = 0.0_rp
840  end do
841  end do
842  do j = jjs, jje
843  do i = iis, iie
844  do k = ke_pbl, ke
845  qflx_sgs_momx(k,i,j,zdir) = 0.0_rp
846  end do
847  end do
848  end do
849 
850 
851  ! integration V
852 !OCL INDEPENDENT
853  !$omp parallel do default(none) &
854  !$omp shared(JJS,JJE,IIS,IIE,KS,KE_PBL,V,dt,SFLX_MV,RCDZ,DENS,GSQRT,I_XYZ,phiN,a,b,c,d) &
855  !$omp private(i,j,k) OMP_SCHEDULE_ collapse(2)
856  do j = jjs, jje
857  do i = iis, iie
858  d(ks,i,j) = v(ks,i,j) + dt * sflx_mv(i,j) * rcdz(ks) / ( dens(ks,i,j) * gsqrt(ks,i,j,i_xyz) )
859  do k = ks+1, ke_pbl
860  d(k,i,j) = v(k,i,j)
861  end do
862  call diffusion_solver( &
863  phin(:,i,j), & ! (out)
864  a(:,i,j), b(:,i,j), c(:,i,j), d(:,i,j), & ! (in)
865  ke_pbl ) ! (in)
866  end do
867  end do
868 
869  end do
870  end do
871 
872  call comm_vars( phin, 1 )
873  call comm_wait( phin, 1 )
874 
875  do jjs = js, je, jblock
876  jje = jjs+jblock-1
877  do iis = is, ie, iblock
878  iie = iis+iblock-1
879 
880  do j = jjs, jje
881  do i = iis, iie
882  qflx_sgs_momy(ks-1,i,j,zdir) = 0.0_rp
883  do k = ks, ke_pbl-1
884  qflx_sgs_momy(k,i,j,zdir) = - 0.03125_rp & ! 1/4/4/2
885  * ( dens(k,i,j) + dens(k+1,i,j) + dens(k,i,j+1) + dens(k+1,i,j+1) ) &
886  * ( nu(k,i,j) + nu(k+1,i,j) + nu(k,i,j+1) + nu(k+1,i,j+1) ) &
887  * ( (phin(k+1,i,j)+phin(k+1,i,j+1)) - (phin(k,i,j)+phin(k,i,j+1)) ) &
888  * j33g * rfdz(k) / gsqrt(k,i,j,i_xvw)
889  end do
890  do k = ke_pbl, ke
891  qflx_sgs_momy(k,i,j,zdir) = 0.0_rp
892  end do
893  end do
894  end do
895 
896 
897  ! for scalars
898  ! integration POTT
899 !OCL INDEPENDENT
900  do j = jjs, jje
901  do i = iis, iie
902  ap = - dt * 0.5_rp * ( rhokh(ks,i,j) + rhokh(ks+1,i,j) ) &
903  * rfdz(ks) / gsqrt(ks,i,j,i_xyw)
904  a(ks,i,j) = ap * rcdz(ks) / (dens(ks,i,j) * gsqrt(ks,i,j,i_xyz) )
905  c(ks,i,j) = 0.0_rp
906  b(ks,i,j) = - a(ks,i,j) + 1.0_rp
907  do k = ks+1, ke_pbl-1
908  c(k,i,j) = ap * rcdz(k) / (dens(k,i,j) * gsqrt(k,i,j,i_xyz) )
909  ap = - dt * 0.5_rp * ( rhokh(k,i,j) + rhokh(k+1,i,j) ) &
910  * rfdz(k) / gsqrt(k,i,j,i_xyw)
911  a(k,i,j) = ap * rcdz(k) / ( dens(k,i,j) * gsqrt(k,i,j,i_xyz) )
912  b(k,i,j) = - a(k,i,j) - c(k,i,j) + 1.0_rp
913  end do
914  a(ke_pbl,i,j) = 0.0_rp
915  c(ke_pbl,i,j) = ap * rcdz(ke_pbl) / (dens(ke_pbl,i,j) * gsqrt(ke_pbl,i,j,i_xyz) )
916  b(ke_pbl,i,j) = - c(ke_pbl,i,j) + 1.0_rp
917  d(ks,i,j) = pott(ks,i,j) + dt * sflx_pt(i,j) * rcdz(ks) / gsqrt(ks,i,j,i_xyz)
918  do k = ks+1, ke_pbl
919  d(k,i,j) = pott(k,i,j)
920  end do
921  call diffusion_solver( &
922  phin(:,i,j), & ! (out)
923  a(:,i,j), b(:,i,j), c(:,i,j), d(:,i,j), & ! (in)
924  ke_pbl ) ! (in)
925  end do
926  end do
927 
928  !$omp parallel do default(none) &
929  !$omp shared(JJS,JJE,IIS,IIE,KS,KE_PBL,RHOKh,J33G,phiN,RFDZ,GSQRT,I_XYW,KE,qflx_sgs_rhot) &
930  !$omp private(i,j,k) OMP_SCHEDULE_ collapse(2)
931  do j = jjs, jje
932  do i = iis, iie
933  qflx_sgs_rhot(ks-1,i,j,zdir) = 0.0_rp
934  do k = ks, ke_pbl-1
935  qflx_sgs_rhot(k,i,j,zdir) = - 0.5_rp & ! 1/2
936  * ( rhokh(k,i,j) + rhokh(k+1,i,j) ) &
937  * j33g * ( phin(k+1,i,j) - phin(k,i,j) ) * rfdz(k) / gsqrt(k,i,j,i_xyw)
938  end do
939  do k = ke_pbl, ke
940  qflx_sgs_rhot(k,i,j,zdir) = 0.0_rp
941  end do
942  end do
943  end do
944 
945 
946  ! TRACER
947  do iq = 1, qa
948 
949  if ( iq == i_tke ) then
950  qflx_sgs_rhoq(:,:,:,zdir,iq) = 0.0_rp
951  cycle
952  end if
953 
954  if( .NOT. tracer_advc(iq) ) cycle
955 
956 !OCL INDEPENDENT
957  !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
958  !$omp shared(JJS,JJE,IIS,IIE,iq,I_QV,d,QTRC,KS,dt,SFLX_Q,RCDZ,DENS,GSQRT,I_XYZ,KE_PBL) &
959  !$omp shared(a,b,c,phiN)
960  do j = jjs, jje
961  do i = iis, iie
962  d(ks,i,j) = qtrc(ks,i,j,iq) &
963  + dt * sflx_q(i,j,iq) * rcdz(ks) / ( dens(ks,i,j) * gsqrt(ks,i,j,i_xyz) )
964 
965  do k = ks+1, ke_pbl
966  d(k,i,j) = qtrc(k,i,j,iq)
967  end do
968 
969  call diffusion_solver( &
970  phin(:,i,j), & ! (out)
971  a(:,i,j), b(:,i,j), c(:,i,j), d(:,i,j), & ! (in)
972  ke_pbl ) ! (in)
973  end do
974  end do
975 
976  !$omp parallel do default(none) &
977  !$omp shared(JJS,JJE,IIS,IIE,KS,KE_PBL,iq,RHOKh,J33G,phiN,RFDZ,GSQRT,KE,I_XYW) &
978  !$omp shared(qflx_sgs_rhoq) &
979  !$omp private(i,j,k) OMP_SCHEDULE_ collapse(2)
980  do j = jjs, jje
981  do i = iis, iie
982  qflx_sgs_rhoq(ks-1,i,j,zdir,iq) = 0.0_rp
983  do k = ks, ke_pbl-1
984  qflx_sgs_rhoq(k,i,j,zdir,iq) = - 0.5_rp & ! 1/2
985  * ( rhokh(k,i,j) + rhokh(k+1,i,j) ) &
986  * j33g * ( phin(k+1,i,j) - phin(k,i,j) ) * rfdz(k) / gsqrt(k,i,j,i_xyw)
987  end do
988  do k = ke_pbl, ke
989  qflx_sgs_rhoq(k,i,j,zdir,iq) = 0.0_rp
990  end do
991  end do
992  end do
993 
994  end do
995 
996 
997  ! TKE
998  !$omp parallel do default(none) private(i,j,k,tke) OMP_SCHEDULE_ collapse(2) &
999  !$omp shared(JJS,JJE,IIS,IIE,KS,q,d,dt,Nu,dudz2,n2,Pr,KE_PBL)
1000  do j = jjs, jje
1001  do i = iis, iie
1002  tke = q(ks,i,j)**2 * 0.5
1003  d(ks,i,j) = tke &
1004  + dt * nu(ks,i,j) * ( dudz2(ks,i,j) - n2(ks,i,j)/pr(ks,i,j) )
1005 #ifdef MORE_HIST
1006  gen(ks,i,j) = nu(ks,i,j) * ( dudz2(ks,i,j) - n2(ks,i,j)/pr(ks,i,j) )
1007 #endif
1008  do k = ks+1, ke_pbl-1
1009  tke = q(k,i,j)**2 * 0.5
1010  d(k,i,j) = tke &
1011  + dt * nu(k,i,j) * ( dudz2(k,i,j) - n2(k,i,j)/pr(k,i,j) )
1012 #ifdef MORE_HIST
1013  gen(k,i,j) = nu(k,i,j) * ( dudz2(k,i,j) - n2(k,i,j)/pr(k,i,j) )
1014 #endif
1015  end do
1016  tke = q(ke_pbl,i,j)**2 * 0.5
1017  d(ke_pbl,i,j) = tke &
1018  + dt * nu(ke_pbl,i,j) * ( dudz2(ke_pbl,i,j) - n2(ke_pbl,i,j)/pr(ke_pbl,i,j) )
1019 #ifdef MORE_HIST
1020  gen(ke_pbl,i,j) = nu(ke_pbl,i,j) * ( dudz2(ke_pbl,i,j) - n2(ke_pbl,i,j)/pr(ke_pbl,i,j) )
1021 #endif
1022  end do
1023  end do
1024 
1025 !OCL INDEPENDENT
1026  !$omp parallel do default(none) &
1027  !$omp shared(JJS,JJE,IIS,IIE,KS,KE_PBL,dt,DENS,Nu,RFDZ,GSQRT,I_XYW,l,RCDZ,I_XYZ,I_TKE) &
1028  !$omp shared(tke_N,QTRC,KE,d,q,a,c,b,RHOQ_t) &
1029  !$omp private(i,j,k,ap) OMP_SCHEDULE_ collapse(2)
1030  do j = jjs, jje
1031  do i = iis, iie
1032  ap = - dt * 1.5_rp * ( dens(ks ,i,j)*nu(ks ,i,j) &
1033  + dens(ks+1,i,j)*nu(ks+1,i,j) ) &
1034  * rfdz(ks) / gsqrt(ks,i,j,i_xyw)
1035  a(ks,i,j) = ap * rcdz(ks) / ( dens(ks,i,j) * gsqrt(ks,i,j,i_xyz) )
1036  c(ks,i,j) = 0.0_rp
1037  b(ks,i,j) = - a(ks,i,j) + 1.0_rp + 2.0_rp * dt * q(ks,i,j) / ( b1 * l(ks,i,j) )
1038  do k = ks+1, ke_pbl-1
1039  c(k,i,j) = ap * rcdz(k) / ( dens(k,i,j) * gsqrt(k,i,j,i_xyz) )
1040  ap = - dt * 1.5_rp * ( dens(k ,i,j)*nu(k ,i,j) &
1041  + dens(k+1,i,j)*nu(k+1,i,j)) &
1042  * rfdz(k) / gsqrt(k,i,j,i_xyw)
1043  a(k,i,j) = ap * rcdz(k) / ( dens(k,i,j) * gsqrt(k,i,j,i_xyz) )
1044  b(k,i,j) = - a(k,i,j) - c(k,i,j) + 1.0_rp + 2.0_rp * dt * q(k,i,j) / ( b1 * l(k,i,j) )
1045  end do
1046 
1047  a(ke_pbl,i,j) = 0.0_rp
1048  c(ke_pbl,i,j) = ap * rcdz(ke_pbl) / ( dens(ke_pbl,i,j) * gsqrt(ke_pbl,i,j,i_xyz) )
1049  b(ke_pbl,i,j) = - c(ke_pbl,i,j) + 1.0_rp + 2.0_rp * dt * q(ke_pbl,i,j) / ( b1 * l(ke_pbl,i,j) )
1050  call diffusion_solver( &
1051  tke_n(:,i,j), & ! (out)
1052  a(:,i,j), b(:,i,j), c(:,i,j), d(:,i,j), & ! (in)
1053  ke_pbl ) ! (in)
1054  do k = ks, ke_pbl
1055  rhoq_t(k,i,j,i_tke) = ( max(tke_n(k,i,j), tke_min) - qtrc(k,i,j,i_tke) ) * dens(k,i,j) / dt
1056  end do
1057  do k = ke_pbl+1, ke
1058  rhoq_t(k,i,j,i_tke) = 0.0_rp
1059  end do
1060 
1061  end do
1062  end do
1063 
1064  end do
1065  end do
1066 
1067 #ifdef MORE_HIST
1068  gen(ke_pbl+1:ke,:,:) = 0.0_rp
1069  dudz2(ke_pbl+1:ke,:,:) = 0.0_rp
1070  l(ke_pbl+1:ke,:,:) = 0.0_rp
1071  potv(ke_pbl+1:ke,:,:) = 0.0_rp
1072  potl(ke_pbl+1:ke,:,:) = 0.0_rp
1073  call hist_in(gen, 'TKE_gen', 'generation of TKE', 'm2/s3', nohalo=.true.)
1074  call hist_in(dudz2, 'dUdZ2', 'dudz2', 'm2/s2', nohalo=.true.)
1075  call hist_in(l, 'L_mix', 'minxing length', 'm', nohalo=.true.)
1076  call hist_in(potv, 'POTV', 'virtual potential temperature', 'K', nohalo=.true.)
1077  call hist_in(potl, 'POTL', 'liquid potential temperature', 'K', nohalo=.true.)
1078 #endif
1079 
1080  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_const::const_cpdry
real(rp), public const_cpdry
specific heat (dry air,constant pressure) [J/kg/K]
Definition: scale_const.F90:58
scale_atmos_saturation
module ATMOSPHERE / Saturation adjustment
Definition: scale_atmos_sub_saturation.F90:16
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_tracer::tracer_r
real(rp), dimension(qa_max), public tracer_r
Definition: scale_tracer.F90:41
scale_grid_index::iblock
integer, public iblock
block size for cache blocking: x
Definition: scale_grid_index.F90:47
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_real::real_cz
real(rp), dimension(:,:,:), allocatable, public real_cz
geopotential height [m] (cell center)
Definition: scale_grid_real.F90:40
scale_grid::grid_rcdz
real(rp), dimension(:), allocatable, public grid_rcdz
reciprocal of center-dz
Definition: scale_grid_cartesian.F90:62
scale_tracer::tracer_cv
real(rp), dimension(qa_max), public tracer_cv
Definition: scale_tracer.F90:40
scale_gridtrans::i_xvw
integer, public i_xvw
Definition: scale_gridtrans.F90:52
scale_const::const_rdry
real(rp), public const_rdry
specific gas constant (dry air) [J/kg/K]
Definition: scale_const.F90:57
scale_grid_index
module grid index
Definition: scale_grid_index.F90:14
scale_tracer
module TRACER
Definition: scale_tracer.F90:15
scale_atmos_hydrometeor
module Hydrometeor
Definition: scale_atmos_sub_hydrometeor.F90:14
scale_gridtrans
module GRIDTRANS
Definition: scale_gridtrans.F90:15
scale_grid_real
module GRID (real space)
Definition: scale_grid_real.F90:14
scale_const::const_lhv0
real(rp), parameter, public const_lhv0
latent heat of vaporizaion at 0C [J/kg]
Definition: scale_const.F90:77
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_comm
module COMMUNICATION
Definition: scale_comm.F90:23
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_gridtrans::i_uyw
integer, public i_uyw
Definition: scale_gridtrans.F90:51
scale_const::const_epstvap
real(rp), public const_epstvap
1 / epsilon - 1
Definition: scale_const.F90:72
scale_const::const_rvap
real(rp), parameter, public const_rvap
specific gas constant (water vapor) [J/kg/K]
Definition: scale_const.F90:65
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_atmos_thermodyn
module ATMOSPHERE / Thermodynamics
Definition: scale_atmos_sub_thermodyn.F90:19
scale_precision
module PRECISION
Definition: scale_precision.F90:13
scale_history
module HISTORY
Definition: scale_history.F90:16
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
scale_const::const_cpvap
real(rp), parameter, public const_cpvap
specific heat (water vapor, constant pressure) [J/kg/K]
Definition: scale_const.F90:66
scale_tracer::tracer_mass
real(rp), dimension(qa_max), public tracer_mass
Definition: scale_tracer.F90:42
Here is the call graph for this function:
Here is the caller graph for this function:

◆ get_length()

subroutine scale_atmos_phy_tb_mynn::get_length ( real(rp), dimension(ka,ia,ja), intent(out)  l,
real(rp), dimension(ka,ia,ja), intent(in)  DENS,
real(rp), dimension(ka,ia,ja), intent(in)  q,
real(rp), dimension(ka,ia,ja), intent(in)  n2,
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_PT,
real(rp), dimension(ka,ia,ja), intent(in)  PT0 
)

Definition at line 1089 of file scale_atmos_phy_tb_mynn.F90.

References scale_const::const_cpdry, scale_const::const_eps, scale_const::const_grav, scale_const::const_karman, scale_grid_index::ie, scale_grid_index::is, scale_grid_index::je, scale_grid_index::js, scale_grid_index::ks, and scale_grid_real::real_fz.

Referenced by atmos_phy_tb_mynn().

1089  use scale_grid_real, only: &
1090  fz => real_fz
1091  use scale_const, only: &
1092  grav => const_grav, &
1093  karman => const_karman, &
1094  cp => const_cpdry, &
1095  eps => const_eps
1096  implicit none
1097 
1098  real(RP), intent(out) :: l(KA,IA,JA)
1099  real(RP), intent(in) :: DENS(KA,IA,JA)
1100  real(RP), intent(in) :: q(KA,IA,JA)
1101  real(RP), intent(in) :: n2(KA,IA,JA)
1102  real(RP), intent(in) :: SFLX_MU(IA,JA)
1103  real(RP), intent(in) :: SFLX_MV(IA,JA)
1104  real(RP), intent(in) :: SFLX_PT(IA,JA)
1105  real(RP), intent(in) :: PT0(KA,IA,JA)
1106 
1107  real(RP) :: ls
1108  real(RP) :: lt
1109  real(RP) :: lb
1110  real(RP) :: rlm
1111  real(RP) :: rlt
1112 
1113  real(RP) :: qc
1114  real(RP) :: int_q
1115  real(RP) :: int_qz
1116  real(RP) :: rn2sr
1117  real(RP) :: us
1118  real(RP) :: zeta
1119 
1120  real(RP) :: z
1121  real(RP) :: qdz
1122 
1123  real(RP) :: sw
1124  integer :: k, i, j
1125 
1126 !OCL LOOP_NOFUSION,PREFETCH_SEQUENTIAL(SOFT),SWP
1127  !$omp parallel do default(none) &
1128  !$omp shared(JS,JE,IS,IE,KS,KE_PBL,q,FZ,EPS,ATMOS_PHY_TB_MYNN_Lt_MAX,SFLX_MU,SFLX_MV,DENS) &
1129  !$omp shared(SFLX_PT,PT0,n2,GRAV,l) &
1130  !$omp private(i,j,k,qdz,lt,rlt,us,rlm,qc,z,zeta,sw,ls,rn2sr,lb,int_qz,int_q) &
1131  !$omp OMP_SCHEDULE_ collapse(2)
1132  do j = js, je
1133  do i = is, ie
1134  int_qz = 0.0_rp
1135  int_q = 0.0_rp
1136  do k = ks, ke_pbl
1137  qdz = q(k,i,j) * ( fz(k,i,j) - fz(k-1,i,j) )
1138  int_qz = int_qz + ((fz(k,i,j)+fz(k-1,i,j))*0.5_rp-fz(ks-1,i,j)) * qdz
1139  int_q = int_q + qdz
1140  end do
1141  ! LT
1142  lt = min( max(0.23_rp * int_qz / (int_q + eps), &
1143  lt_min), &
1144  atmos_phy_tb_mynn_lt_max )
1145  rlt = 1.0_rp / lt
1146 
1147  us = sqrt( sqrt( sflx_mu(i,j)**2 + sflx_mv(i,j)**2 ) / dens(ks,i,j) ) ! friction velocity
1148  us = max(us, us_min)
1149  rlm = - karman * grav * sflx_pt(i,j) / (pt0(ks,i,j) * us**3 )
1150 
1151  qc = (grav/pt0(ks,i,j)*max(sflx_pt(i,j),0.0_rp)*lt)**oneoverthree
1152 
1153  do k = ks, ke_pbl
1154  z = ( fz(k,i,j)+fz(k-1,i,j) )*0.5_rp - fz(ks-1,i,j)
1155  zeta = z * rlm
1156 
1157  ! LS
1158  sw = sign(0.5_rp, zeta) + 0.5_rp ! 1 for zeta >= 0, 0 for zeta < 0
1159  ls = karman * z &
1160  * ( sw / (1.0_rp + 2.7_rp*min(zeta,1.0_rp)*sw ) &
1161  + ( (1.0_rp - 100.0_rp*zeta)*(1.0_rp-sw) )**0.2_rp )
1162 
1163  ! LB
1164  sw = sign(0.5_rp, n2(k,i,j)-eps) + 0.5_rp ! 1 for dptdz >0, 0 for dptdz <= 0
1165  rn2sr = 1.0_rp / ( sqrt(n2(k,i,j)*sw) + 1.0_rp-sw)
1166  lb = (1.0_rp + 5.0_rp * sqrt(qc*rn2sr/lt)) * q(k,i,j) * rn2sr * sw & ! qc=0 when SFLX_PT < 0
1167  + 999.e10_rp * (1.0_rp-sw)
1168 
1169  ! L
1170  l(k,i,j) = 1.0_rp / ( 1.0_rp/ls + rlt + 1.0_rp/lb )
1171  end do
1172  end do
1173  end do
1174 
1175  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_const::const_cpdry
real(rp), public const_cpdry
specific heat (dry air,constant pressure) [J/kg/K]
Definition: scale_const.F90:58
scale_grid_real::real_fz
real(rp), dimension(:,:,:), allocatable, public real_fz
geopotential height [m] (cell face )
Definition: scale_grid_real.F90:41
scale_const::const_karman
real(rp), parameter, public const_karman
von Karman constant
Definition: scale_const.F90:52
scale_grid_real
module GRID (real space)
Definition: scale_grid_real.F90:14
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_const
module CONSTANT
Definition: scale_const.F90:14
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
Here is the caller graph for this function:

◆ get_q2_level2()

subroutine scale_atmos_phy_tb_mynn::get_q2_level2 ( real(rp), dimension(ka,ia,ja), intent(out)  q2_2,
real(rp), dimension(ka,ia,ja), intent(in)  dudz2,
real(rp), dimension(ka,ia,ja), intent(in)  Ri,
real(rp), dimension(ka,ia,ja), intent(in)  l 
)

Definition at line 1181 of file scale_atmos_phy_tb_mynn.F90.

References scale_grid_index::ie, scale_grid_index::is, scale_grid_index::je, scale_grid_index::js, and scale_grid_index::ks.

Referenced by atmos_phy_tb_mynn().

1181  implicit none
1182 
1183  real(RP), intent(out) :: q2_2(KA,IA,JA)
1184  real(RP), intent(in) :: dudz2(KA,IA,JA)
1185  real(RP), intent(in) :: Ri(KA,IA,JA)
1186  real(RP), intent(in) :: l(KA,IA,JA)
1187 
1188  real(RP) :: rf
1189  real(RP) :: sm_2
1190  real(RP) :: sh_2
1191 
1192  real(RP) :: q2
1193  integer :: k, i, j
1194 
1195 !OCL LOOP_NOFUSION,PREFETCH_SEQUENTIAL(SOFT),SWP
1196  !$omp parallel do default(none) &
1197  !$omp shared(JS,JE,IS,IE,KS,KE_PBL,AF12,Ri,Rf1,Rf2,Rfc,A2,G2,l,dudz2,q2_2) &
1198  !$omp private(i,j,k,rf,sh_2,sm_2,q2) OMP_SCHEDULE_ collapse(2)
1199  do j = js, je
1200  do i = is, ie
1201  do k = ks, ke_pbl
1202  rf = min(0.5_rp / af12 * ( ri(k,i,j) &
1203  + af12*rf1 &
1204  - sqrt(ri(k,i,j)**2 + 2.0_rp*af12*(rf1-2.0_rp*rf2)*ri(k,i,j) + (af12*rf1)**2) ), &
1205  rfc)
1206  sh_2 = 3.0_rp * a2 * (g1+g2) * (rfc-rf) / (1.0_rp-rf)
1207  sm_2 = sh_2 * af12 * (rf1-rf) / (rf2-rf)
1208  q2 = b1 * l(k,i,j)**2 * sm_2 * (1.0_rp-rf) * dudz2(k,i,j)
1209  q2_2(k,i,j) = max( q2, 1.e-10_rp )
1210  end do
1211  end do
1212  end do
1213 
1214  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_grid_index::js
integer, public js
start point of inner domain: y, local
Definition: scale_grid_index.F90:62
scale_grid_index::ie
integer, public ie
end point of inner domain: x, local
Definition: scale_grid_index.F90:61
Here is the caller graph for this function:
Generated by   doxygen 1.8.14