SCALE-RM
Functions/Subroutines | Variables
scale_atmos_phy_bl_mynn Module Reference

module atmosphere / physics / pbl / mynn More...

Functions/Subroutines

subroutine, public atmos_phy_bl_mynn_tracer_setup
 ATMOS_PHY_BL_MYNN_tracer_setup Tracer Setup. More...
 
subroutine, public atmos_phy_bl_mynn_setup (KA, KS, KE, IA, IS, IE, JA, JS, JE, BULKFLUX_type, TKE_MIN, PBL_MAX)
 ATMOS_PHY_BL_MYNN_setup Setup. More...
 
subroutine, public atmos_phy_bl_mynn_tendency (KA, KS, KE, IA, IS, IE, JA, JS, JE, DENS, U, V, POTT, PROG, PRES, EXNER, N2, QDRY, QV, Qw, POTL, POTV, SFC_DENS, SFLX_MU, SFLX_MV, SFLX_SH, SFLX_QV, us, ts, qs, RLmo, CZ, FZ, dt_DP, BULKFLUX_type, RHOU_t, RHOV_t, RHOT_t, RHOQV_t, RPROG_t, Nu, Kh, Qlp, cldfrac, Zi)
 ATMOS_PHY_BL_MYNN_tendency calculate tendency by the virtical eddy viscosity. More...
 
subroutine get_gamma_implicit (KA, KS, KE, i, j, tsq, qsq, cov, dtldz, dqwdz, POTV, prod_t, prod_q, prod_c, betat, betaq, f_gamma, l, q, dt, gammat, gammaq)
 
subroutine calc_vertical_differece (KA, KS, KE, i, j, U, V, POTL, Qw, CDZ, FDZ, F2H, dudz2, dtldz, dqwdz)
 
subroutine partial_condensation (KA, KS, KE, PRES, POTT, POTL, Qw, Qdry, EXNER, dtldz, dqwdz, l, sh25, ac, tsq, qsq, cov, mynn_level3, betat, betaq, Qlp, cldfrac)
 
subroutine get_f2h (KA, KS, KE, CDZ, f2h)
 

Variables

integer, public atmos_phy_bl_mynn_ntracer
 
character(len=h_short), dimension(:), allocatable, public atmos_phy_bl_mynn_name
 
character(len=h_long), dimension(:), allocatable, public atmos_phy_bl_mynn_desc
 
character(len=h_short), dimension(:), allocatable, public atmos_phy_bl_mynn_units
 

Detailed Description

module atmosphere / physics / pbl / mynn

Description
Boundary layer turbulence model Mellor-Yamada Nakanishi-Niino model
Author
Team SCALE
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_BL_MYNN
    nametypedefault valuecomment
    ATMOS_PHY_BL_MYNN_PBL_MAX real(RP) 4.E+3_RP > maximum height of the PBL
    ATMOS_PHY_BL_MYNN_N2_MAX real(RP) 1.E1_RP
    ATMOS_PHY_BL_MYNN_NU_MIN real(RP) -1.E1_RP
    ATMOS_PHY_BL_MYNN_NU_MAX real(RP) 1.E4_RP
    ATMOS_PHY_BL_MYNN_KH_MIN real(RP) -1.E1_RP
    ATMOS_PHY_BL_MYNN_KH_MAX real(RP) 1.E4_RP
    ATMOS_PHY_BL_MYNN_LT_MAX real(RP) 700.0_RP ~ 0.23 * 3 km
    ATMOS_PHY_BL_MYNN_LEVEL character(len=H_SHORT) "2.5" ! "2.5" or "3", level 3 is under experimental yet.
    ATMOS_PHY_BL_MYNN_SQ_FACT real(RP) 3.0_RP
    ATMOS_PHY_BL_MYNN_INIT_TKE logical .false.
    ATMOS_PHY_BL_MYNN_SIMILARITY logical .true.

History Output
namedescriptionunitvariable
L_mix_MYNN minxing length m l
Pr_MYNN Prandtl number 1 Pr
Ri_MYNN Richardson number 1 Ri
TKE_diss_MYNN TKE dissipation m2/s3 diss
TKE_prod_MYNN TKE production m2/s3 prod
ZFLX_{TRACER_NAME}_MYNN Z FLUX of DENS * {TRACER_NAME} (MYNN);
{TRACER_NAME} depends on the physics schemes, e.g., QV, QC, QR. 		kg/m2/s flx
ZFLX_QV_MYNN Z FLUX of RHOQV (MYNN) kg/m2/s flxQ
ZFLX_RHOT_MYNN Z FLUX of RHOT (MYNN) K kg/m2/s flxT
ZFLX_RHOU_MYNN Z FLUX of RHOU (MYNN) kg/m/s2 flxU
ZFLX_RHOV_MYNN Z FLUX of RHOV (MYNN) kg/m/s2 flxV
dUdZ2_MYNN dudz2 m2/s2 dudz2

Function/Subroutine Documentation

◆ atmos_phy_bl_mynn_tracer_setup()

subroutine, public scale_atmos_phy_bl_mynn::atmos_phy_bl_mynn_tracer_setup

ATMOS_PHY_BL_MYNN_tracer_setup Tracer Setup.

Definition at line 134 of file scale_atmos_phy_bl_mynn.F90.

134  use scale_prc, only: &
135  prc_abort
136  implicit none
137 
138  integer :: ierr
139 
140  log_newline
141  log_info("ATMOS_PHY_BL_MYNN_tracer_setup",*) 'Tracer Setup'
142  log_info("ATMOS_PHY_BL_MYNN_tracer_setup",*) 'Mellor-Yamada Nakanishi-Niino scheme'
143 
144  !--- read namelist
145  rewind(io_fid_conf)
146  read(io_fid_conf,nml=param_atmos_phy_bl_mynn,iostat=ierr)
147  if( ierr > 0 ) then !--- fatal error
148  log_error("ATMOS_PHY_BL_MYNN_tracer_setup",*) 'Not appropriate names in namelist PARAM_ATMOS_PHY_BL_MYNN. Check!'
149  call prc_abort
150  endif
151 
152  select case ( atmos_phy_bl_mynn_level )
153  case ( "2.5" )
154  atmos_phy_bl_mynn_ntracer = 1
155  allocate( atmos_phy_bl_mynn_name(1), atmos_phy_bl_mynn_desc(1), atmos_phy_bl_mynn_units(1) )
156  atmos_phy_bl_mynn_name(:) = (/ 'TKE_MYNN' /)
157  atmos_phy_bl_mynn_desc(:) = (/ 'turbulent kinetic energy (MYNN)' /)
158  atmos_phy_bl_mynn_units(:) = (/ 'm2/s2' /)
159  case ( "3" )
160  atmos_phy_bl_mynn_ntracer = 4
161  allocate( atmos_phy_bl_mynn_name(4), atmos_phy_bl_mynn_desc(4), atmos_phy_bl_mynn_units(4) )
162  atmos_phy_bl_mynn_name(:) = (/ 'TKE_MYNN', &
163  'TSQ_MYNN', &
164  'QSQ_MYNN', &
165  'COV_MYNN' /)
166  atmos_phy_bl_mynn_desc(:) = (/ 'turbulent kinetic energy (MYNN) ', &
167  'sub-grid variance of liquid water potential temperature (MYNN) ', &
168  'sub-grid variance of total water content (MYNN) ', &
169  'sub-grid covariance of liquid water potential temperature and total water content (MYNN)' /)
170  atmos_phy_bl_mynn_units(:) = (/ 'm2/s2 ', &
171  'K2 ', &
172  'kg2/kg2', &
173  'K kg ' /)
174  case default
175  log_error("ATMOS_PHY_BL_MYNN_tracer_setup",*) 'only level 2.5 and 3 are supported at this moment'
176  call prc_abort
177  end select
178 
179  return

References atmos_phy_bl_mynn_desc, atmos_phy_bl_mynn_name, atmos_phy_bl_mynn_ntracer, atmos_phy_bl_mynn_units, scale_io::io_fid_conf, and scale_prc::prc_abort().

Referenced by mod_atmos_phy_bl_driver::atmos_phy_bl_driver_tracer_setup().

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

◆ atmos_phy_bl_mynn_setup()

subroutine, public scale_atmos_phy_bl_mynn::atmos_phy_bl_mynn_setup ( integer, intent(in)  KA,
integer, intent(in)  KS,
integer, intent(in)  KE,
integer, intent(in)  IA,
integer, intent(in)  IS,
integer, intent(in)  IE,
integer, intent(in)  JA,
integer, intent(in)  JS,
integer, intent(in)  JE,
character(len=*), intent(in)  BULKFLUX_type,
real(rp), intent(in), optional  TKE_MIN,
real(rp), intent(in), optional  PBL_MAX 
)

ATMOS_PHY_BL_MYNN_setup Setup.

Definition at line 190 of file scale_atmos_phy_bl_mynn.F90.

190  use scale_prc, only: &
191  prc_abort
192  use scale_const, only: &
193  pi => const_pi
194  implicit none
195 
196  integer, intent(in) :: KA, KS, KE
197  integer, intent(in) :: IA, IS, IE
198  integer, intent(in) :: JA, JS, JE
199 
200  character(len=*), intent(in) :: BULKFLUX_type
201 
202  real(RP), intent(in), optional :: TKE_MIN
203  real(RP), intent(in), optional :: PBL_MAX
204 
205  integer :: ierr
206  integer :: k, i, j
207  !---------------------------------------------------------------------------
208 
209  log_newline
210  log_info("ATMOS_PHY_BL_MYNN_setup",*) 'Setup'
211  log_info("ATMOS_PHY_BL_MYNN_setup",*) 'Mellor-Yamada Nakanishi-Niino scheme'
212 
213  if ( present(tke_min) ) atmos_phy_bl_mynn_tke_min = tke_min
214  if ( present(pbl_max) ) atmos_phy_bl_mynn_pbl_max = pbl_max
215 
216  !--- read namelist
217  rewind(io_fid_conf)
218  read(io_fid_conf,nml=param_atmos_phy_bl_mynn,iostat=ierr)
219  if( ierr < 0 ) then !--- missing
220  log_info("ATMOS_PHY_BL_MYNN_setup",*) 'Not found namelist. Default used.'
221  elseif( ierr > 0 ) then !--- fatal error
222  log_error("ATMOS_PHY_BL_MYNN_setup",*) 'Not appropriate names in namelist PARAM_ATMOS_PHY_BL_MYNN. Check!'
223  call prc_abort
224  endif
225  log_nml(param_atmos_phy_bl_mynn)
226 
227  a1 = b1 * (1.0_rp - 3.0_rp * g1) / 6.0_rp
228  a2 = 1.0_rp / (3.0_rp * g1 * b1**(1.0_rp/3.0_rp) * prn )
229  c1 = g1 - 1.0_rp / ( 3.0_rp * a1 * b1**(1.0_rp/3.0_rp) )
230  g2 = ( 2.0_rp * a1 * (3.0_rp - 2.0_rp * c2) + b2 * (1.0_rp - c3) ) / b1
231  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)
232  f2 = b1 * (g1 + g2) - 3.0_rp * a1 * (1.0_rp - c2)
233 
234  rf1 = b1 * (g1 - c1) / f1
235  rf2 = b1 * g1 / f2
236  rfc = g1 / (g1 + g2)
237 
238  af12 = a1 * f1 / ( a2 * f2 )
239 
240  sqrt_2pi = sqrt( 2.0_rp * pi )
241  rsqrt_2pi = 1.0_rp / sqrt_2pi
242  rsqrt_2 = 1.0_rp / sqrt( 2.0_rp )
243 
244  if ( atmos_phy_bl_mynn_level == "3" ) then
245  log_warn("ATMOS_PHY_BL_MYNN_setup", *) "At this moment, level 3 is still experimental"
246  end if
247 
248  initialize = atmos_phy_bl_mynn_init_tke
249 
250  select case ( bulkflux_type )
251  case ( "B91", "B91W01" )
252  ! do nothing
253  case default
254  atmos_phy_bl_mynn_similarity = .false.
255  end select
256 
257  return

References scale_const::const_pi, scale_io::io_fid_conf, and scale_prc::prc_abort().

Referenced by mod_atmos_phy_bl_driver::atmos_phy_bl_driver_setup().

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

◆ atmos_phy_bl_mynn_tendency()

subroutine, public scale_atmos_phy_bl_mynn::atmos_phy_bl_mynn_tendency ( integer, intent(in)  KA,
integer, intent(in)  KS,
integer, intent(in)  KE,
integer, intent(in)  IA,
integer, intent(in)  IS,
integer, intent(in)  IE,
integer, intent(in)  JA,
integer, intent(in)  JS,
integer, intent(in)  JE,
real(rp), dimension (mynn), intent(out)  DENS,
real(rp), dimension (ka,ia,ja), intent(out)  U,
real(rp), dimension (ka,ia,ja), intent(out)  V,
real(rp), dimension (ka,ia,ja), intent(in)  POTT,
real(rp), dimension (ka,ia,ja,atmos_phy_bl_mynn_ntracer), intent(in)  PROG,
real(rp), dimension (ka,ia,ja), intent(in)  PRES,
real(rp), dimension (ka,ia,ja), intent(in)  EXNER,
real(rp), dimension (ka,ia,ja), intent(in)  N2,
real(rp), dimension (ka,ia,ja), intent(in)  QDRY,
real(rp), dimension (ka,ia,ja), intent(out)  QV,
real(rp), dimension (ka,ia,ja), intent(in)  Qw,
real(rp), dimension (ka,ia,ja), intent(in)  POTL,
real(rp), dimension (ka,ia,ja), intent(in)  POTV,
real(rp), dimension( ia,ja), intent(in)  SFC_DENS,
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), intent(in)  SFLX_QV,
real(rp), dimension ( ia,ja), intent(in)  us,
real(rp), dimension ( ia,ja), intent(in)  ts,
real(rp), dimension ( ia,ja), intent(in)  qs,
real(rp), dimension ( ia,ja), intent(in)  RLmo,
real(rp), dimension( ka,ia,ja), intent(in)  CZ,
real(rp), dimension(0:ka,ia,ja), intent(in)  FZ,
real(dp), intent(in)  dt_DP,
character(len=*), intent(in)  BULKFLUX_type,
real(rp), dimension (ka,ia,ja), intent(out)  RHOU_t,
real(rp), dimension (ka,ia,ja), intent(out)  RHOV_t,
real(rp), dimension (ka,ia,ja), intent(out)  RHOT_t,
real(rp), dimension(ka,ia,ja), intent(out)  RHOQV_t,
real(rp), dimension(ka,ia,ja,atmos_phy_bl_mynn_ntracer), intent(out)  RPROG_t,
real(rp), dimension (ka,ia,ja), intent(out)  Nu,
real(rp), dimension (ka,ia,ja), intent(out)  Kh,
real(rp), dimension (ka,ia,ja), intent(out)  Qlp,
real(rp), dimension(ka,ia,ja), intent(out)  cldfrac,
real(rp), dimension (ia,ja), intent(out)  Zi 
)

ATMOS_PHY_BL_MYNN_tendency calculate tendency by the virtical eddy viscosity.

Definition at line 276 of file scale_atmos_phy_bl_mynn.F90.

276  use scale_const, only: &
277  eps => const_eps, &
278  grav => const_grav, &
279  karman => const_karman, &
280  cpdry => const_cpdry, &
281  epstvap => const_epstvap
282  use scale_prc, only: &
283  prc_abort
284  use scale_atmos_hydrometeor, only: &
285  hydrometeor_lhv => atmos_hydrometeor_lhv
286  use scale_matrix, only: &
287  matrix_solver_tridiagonal
288  use scale_file_history, only: &
289  file_history_in
290  implicit none
291 
292  integer, intent(in) :: KA, KS, KE
293  integer, intent(in) :: IA, IS, IE
294  integer, intent(in) :: JA, JS, JE
295 
296  real(RP), intent(in) :: DENS (KA,IA,JA)
297  real(RP), intent(in) :: U (KA,IA,JA)
298  real(RP), intent(in) :: V (KA,IA,JA)
299  real(RP), intent(in) :: POTT (KA,IA,JA)
300  real(RP), intent(in) :: PROG (KA,IA,JA,ATMOS_PHY_BL_MYNN_ntracer)
301  real(RP), intent(in) :: PRES (KA,IA,JA)
302  real(RP), intent(in) :: EXNER (KA,IA,JA)
303  real(RP), intent(in) :: N2 (KA,IA,JA)
304  real(RP), intent(in) :: QDRY (KA,IA,JA)
305  real(RP), intent(in) :: QV (KA,IA,JA)
306  real(RP), intent(in) :: Qw (KA,IA,JA)
307  real(RP), intent(in) :: POTL (KA,IA,JA)
308  real(RP), intent(in) :: POTV (KA,IA,JA)
309  real(RP), intent(in) :: SFC_DENS( IA,JA)
310  real(RP), intent(in) :: SFLX_MU ( IA,JA)
311  real(RP), intent(in) :: SFLX_MV ( IA,JA)
312  real(RP), intent(in) :: SFLX_SH ( IA,JA)
313  real(RP), intent(in) :: SFLX_QV ( IA,JA)
314  real(RP), intent(in) :: us ( IA,JA)
315  real(RP), intent(in) :: ts ( IA,JA)
316  real(RP), intent(in) :: qs ( IA,JA)
317  real(RP), intent(in) :: RLmo ( IA,JA)
318 
319  real(RP), intent(in) :: CZ( KA,IA,JA)
320  real(RP), intent(in) :: FZ(0:KA,IA,JA)
321  real(DP), intent(in) :: dt_DP
322 
323  character(len=*), intent(in) :: BULKFLUX_type
324 
325  real(RP), intent(out) :: RHOU_t (KA,IA,JA)
326  real(RP), intent(out) :: RHOV_t (KA,IA,JA)
327  real(RP), intent(out) :: RHOT_t (KA,IA,JA)
328  real(RP), intent(out) :: RHOQV_t(KA,IA,JA)
329  real(RP), intent(out) :: RPROG_t(KA,IA,JA,ATMOS_PHY_BL_MYNN_ntracer)
330  real(RP), intent(out) :: Nu (KA,IA,JA)
331  real(RP), intent(out) :: Kh (KA,IA,JA)
332  real(RP), intent(out) :: Qlp (KA,IA,JA)
333  real(RP), intent(out) :: cldfrac(KA,IA,JA)
334  real(RP), intent(out) :: Zi (IA,JA)
335 
336  real(RP) :: Ri (KA,IA,JA)
337  real(RP) :: Pr (KA,IA,JA)
338  real(RP) :: prod (KA,IA,JA)
339  real(RP) :: diss (KA,IA,JA)
340  real(RP) :: dudz2(KA,IA,JA)
341  real(RP) :: l (KA,IA,JA)
342  real(RP) :: rho_h(KA)
343 
344  real(RP) :: flxU(KA,IA,JA)
345  real(RP) :: flxV(KA,IA,JA)
346  real(RP) :: flxT(KA,IA,JA)
347  real(RP) :: flxQ(KA,IA,JA)
348 
349  real(RP) :: RHO (KA)
350  real(RP) :: RHONu (KA)
351  real(RP) :: RHOKh (KA)
352  real(RP) :: N2_new(KA)
353  real(RP) :: SFLX_PT
354  real(RP) :: SFLX_PTV
355  real(RP) :: sm25 (KA)
356  real(RP) :: smp (KA)
357  real(RP) :: sh25 (KA)
358  real(RP) :: shpgh (KA)
359  real(RP) :: Nu_f (KA)
360  real(RP) :: Kh_f (KA)
361  real(RP) :: q (KA)
362  real(RP) :: q2_2 (KA)
363  real(RP) :: ac (KA)
364  real(RP) :: lq (KA)
365 
366  ! for level 3
367  real(RP) :: tsq (KA)
368  real(RP) :: qsq (KA)
369  real(RP) :: cov (KA)
370  real(RP) :: wtl
371  real(RP) :: wqw
372  real(RP) :: prod_t(KA)
373  real(RP) :: prod_q(KA)
374  real(RP) :: prod_c(KA)
375  real(RP) :: diss_p(KA)
376 
377  real(RP) :: dtldz(KA)
378  real(RP) :: dqwdz(KA)
379  real(RP) :: betat(KA)
380  real(RP) :: betaq(KA)
381 
382  real(RP) :: gammat (KA)
383  real(RP) :: gammaq (KA)
384  real(RP) :: f_gamma(KA)
385 
386  real(RP) :: rlqsm_h(KA)
387 
388  real(RP) :: flx(KA)
389  real(RP) :: a(KA)
390  real(RP) :: b(KA)
391  real(RP) :: c(KA)
392  real(RP) :: d(KA)
393  real(RP) :: ap
394  real(RP) :: phi_N(KA)
395  real(RP) :: tke_P(KA)
396  real(RP) :: dummy(KA)
397 
398  real(RP) :: sf_t
399  real(RP) :: us3
400  real(RP) :: zeta
401  real(RP) :: phi_m, phi_h
402 
403  real(RP) :: FDZ(KA)
404  real(RP) :: CDZ(KA)
405  real(RP) :: f2h(KA,2)
406  real(RP) :: z1
407 
408  logical :: mynn_level3
409 
410  real(RP) :: dt
411 
412  real(RP) :: fmin
413  integer :: kmin
414 
415  real(RP) :: tmp, sw
416 
417  integer :: KE_PBL
418  integer :: k, i, j
419  integer :: nit, it
420  !---------------------------------------------------------------------------
421 
422  dt = real(dt_dp, rp)
423 
424  log_progress(*) "atmosphere / physics / pbl / MYNN"
425 
426  mynn_level3 = ( atmos_phy_bl_mynn_level == "3" )
427 
428 
429 !OCL INDEPENDENT
430  !$omp parallel do default(none) &
431  !$omp OMP_SCHEDULE_ collapse(2) &
432  !$omp shared(KA,KS,KE,IS,IE,JS,JE, &
433  !$omp EPS,GRAV,CPdry,EPSTvap,UNDEF,RSQRT_2,SQRT_2PI,RSQRT_2PI, &
434  !$omp ATMOS_PHY_BL_MYNN_N2_MAX,ATMOS_PHY_BL_MYNN_TKE_MIN, &
435  !$omp ATMOS_PHY_BL_MYNN_NU_MIN,ATMOS_PHY_BL_MYNN_NU_MAX, &
436  !$omp ATMOS_PHY_BL_MYNN_KH_MIN,ATMOS_PHY_BL_MYNN_KH_MAX, &
437  !$omp ATMOS_PHY_BL_MYNN_Sq_fact,ATMOS_PHY_BL_MYNN_similarity, &
438  !$omp ATMOS_PHY_BL_MYNN_PBL_MAX, &
439  !$omp RHOU_t,RHOV_t,RHOT_t,RHOQV_t,RPROG_t,Nu,Kh,Qlp,cldfrac,Zi, &
440  !$omp DENS,PROG,U,V,POTT,PRES,QDRY,QV,Qw,POTV,POTL,EXNER,N2, &
441  !$omp SFC_DENS,SFLX_MU,SFLX_MV,SFLX_SH,SFLX_QV,us,ts,qs,RLmo, &
442  !$omp mynn_level3,initialize, &
443  !$omp CZ,FZ,dt, &
444  !$omp BULKFLUX_type, &
445  !$omp Ri,Pr,prod,diss,dudz2,l,flxU,flxV,flxT,flxQ) &
446  !$omp private(N2_new,lq,sm25,smp,sh25,shpgh,rlqsm_h,Nu_f,Kh_f,q,q2_2,ac, &
447  !$omp SFLX_PT,SFLX_PTV,RHO,RHONu,RHOKh, &
448  !$omp dtldz,dqwdz,betat,betaq,gammat,gammaq,f_gamma,wtl,wqw, &
449  !$omp flx,a,b,c,d,ap,rho_h,phi_n,tke_P,sf_t,zeta,phi_m,phi_h,us3,CDZ,FDZ,f2h,z1, &
450  !$omp dummy, &
451  !$omp tsq,qsq,cov, &
452  !$omp prod_t,prod_q,prod_c,diss_p, &
453  !$omp fmin,kmin, &
454  !$omp sw,tmp, &
455  !$omp KE_PBL,k,i,j,it,nit)
456  do j = js, je
457  do i = is, ie
458 
459  ke_pbl = ks+1
460  do k = ks+2, ke-1
461  if ( atmos_phy_bl_mynn_pbl_max >= cz(k,i,j) - fz(ks-1,i,j) ) then
462  ke_pbl = k
463  else
464  exit
465  end if
466  end do
467 
468  if ( initialize ) then
469  nit = ke_pbl - 1
470  else
471  nit = 1
472  end if
473 
474 
475  z1 = cz(ks,i,j) - fz(ks-1,i,j)
476 
477  do k = ks, ke_pbl
478  fdz(k) = cz(k+1,i,j) - cz(k ,i,j)
479  end do
480  do k = ks, ke_pbl+1
481  cdz(k) = fz(k ,i,j) - fz(k-1,i,j)
482  end do
483 
484  call get_f2h( &
485  ka, ks, ke_pbl, &
486  cdz(:), & ! (in)
487  f2h(:,:) ) ! (out)
488 
489  call calc_vertical_differece( ka, ks, ke_pbl, &
490  i, j, &
491  u(:,i,j), v(:,i,j), potl(:,i,j), qw(:,i,j), & ! (in)
492  cdz(:), fdz(:), f2h(:,:), & ! (in)
493  dudz2(:,i,j), dtldz(:), dqwdz(:) ) ! (out)
494 
495  do k = ks, ke_pbl
496  n2_new(k) = min( max( n2(k,i,j), - atmos_phy_bl_mynn_n2_max ), atmos_phy_bl_mynn_n2_max )
497 ! n2_new(k) = GRAV * POTV(k,i,j) * dtldz(k)
498  ri(k,i,j) = n2_new(k) / dudz2(k,i,j)
499  end do
500 
501  do k = ks, ke_pbl
502  q(k) = sqrt( max( prog(k,i,j,i_tke), atmos_phy_bl_mynn_tke_min ) * 2.0_rp )
503  end do
504 
505  if ( mynn_level3 ) then
506  do k = ks, ke_pbl
507  tsq(k) = max( prog(k,i,j,i_tsq), 0.0_rp )
508  qsq(k) = max( prog(k,i,j,i_qsq), 0.0_rp )
509  cov(k) = prog(k,i,j,i_cov)
510  cov(k) = sign( min( abs(cov(k)), sqrt(tsq(k) * qsq(k))), cov(k) )
511  end do
512  end if
513 
514 
515  flx(ks-1 ) = 0.0_rp
516  flx(ke_pbl) = 0.0_rp
517 
518  do it = 1, nit
519 
520  if ( initialize .and. it==1 ) then
521  do k = ks, ke_pbl
522  q(k) = 1e10_rp
523  end do
524  end if
525 
526  ! length
527  us3 = us(i,j)**3
528  sflx_ptv = - us3 * rlmo(i,j) / karman
529  call get_length( &
530  ka, ks, ke_pbl, &
531  i, j, &
532  q(:), n2_new(:), & ! (in)
533  sflx_ptv, rlmo(i,j), & ! (in)
534  cz(:,i,j), & ! (in)
535  fz(:,i,j), & ! (in)
536  l(:,i,j) ) ! (out)
537 
538  call get_q2_level2( &
539  ka, ks, ke_pbl, &
540  dudz2(:,i,j), ri(:,i,j), l(:,i,j), & ! (in)
541  q2_2(:) ) ! (out)
542 
543  if ( initialize .and. it==1 ) then
544  do k = ks, ke_pbl
545  q(k) = sqrt( q2_2(k) )
546  end do
547  end if
548 
549  do k = ks, ke_pbl
550  ac(k) = min( q(k) / sqrt( q2_2(k) + 1e-20_rp ), 1.0_rp )
551  end do
552 
553  call get_smsh( &
554  ka, ks, ke_pbl, & ! (in)
555  i, j, & ! (in)
556  q(:), ac(:), & ! (in)
557  l(:,i,j), n2_new(:), & ! (in)
558  potv(:,i,j), dudz2(:,i,j), & ! (in)
559  dtldz(:), dqwdz(:), & ! (in)
560  betat(:), betaq(:), & ! (in)
561  .false., .false., & ! (in)
562  tsq(:), qsq(:), cov(:), & ! (inout)
563  sm25(:), smp(:), & ! (out)
564  sh25(:), shpgh(:), & ! (out) ! dummy
565  gammat(:), gammaq(:), & ! (out) ! dummy
566  f_gamma(:) ) ! (out) ! dummy
567 
568 
569  call partial_condensation( ka, ks, ke_pbl, &
570  pres(:,i,j), pott(:,i,j), & ! (in)
571  potl(:,i,j), qw(:,i,j), & ! (in)
572  qdry(:,i,j), exner(:,i,j), & ! (in)
573  dtldz(:), dqwdz(:), & ! (in)
574  l(:,i,j), sh25(:), ac(:), & ! (in)
575  tsq(:), qsq(:), cov(:), & ! (in)
576  mynn_level3, & ! (in)
577  betat(:), betaq(:), & ! (out)
578  qlp(:,i,j), cldfrac(:,i,j) ) ! (out)
579 
580  if ( atmos_phy_bl_mynn_similarity ) then
581 
582  zeta = min( max( z1 * rlmo(i,j), zeta_min ), zeta_max )
583 
584  select case ( bulkflux_type )
585 !!$ case ( 'B71' )
586 !!$ ! Businger et al. (1971)
587 !!$ if ( zeta >= 0 ) then
588 !!$ phi_m = 4.7_RP * zeta + 1.0_RP
589 !!$ phi_h = 4.7_RP * zeta + 0.74_RP
590 !!$ else
591 !!$ phi_m = 1.0_RP / sqrt(sqrt( 1.0_RP - 15.0_RP * zeta ))
592 !!$ phi_h = 0.47_RP / sqrt( 1.0_RP - 9.0_RP * zeta )
593 !!$ end if
594  case ( 'B91', 'B91W01' )
595  ! Beljaars and Holtslag (1991)
596  if ( zeta >= 0 ) then
597  tmp = - 2.0_rp / 3.0_rp * ( 0.35_rp * zeta - 6.0_rp ) * exp(-0.35_rp*zeta)
598  phi_m = tmp * zeta + zeta + 1.0_rp
599  phi_h = tmp + zeta * sqrt( 1.0_rp + 2.0_rp * zeta / 3.0_rp ) + 1.0_rp
600  else
601  if ( bulkflux_type == 'B91W01' ) then
602  ! Wilson (2001)
603 ! tmp = (-zeta)**(2.0_RP/3.0_RP)
604  tmp = abs(zeta)**(2.0_rp/3.0_rp)
605  phi_m = 1.0_rp / sqrt( 1.0_rp + 3.6_rp * tmp )
606  phi_h = 0.95_rp / sqrt( 1.0_rp + 7.9_rp * tmp )
607  else
608 ! tmp = sqrt( 1.0_RP - 16.0_RP * zeta )
609  tmp = sqrt( 1.0_rp + 16.0_rp * abs(zeta) )
610  phi_m = 1.0_rp / sqrt(tmp)
611  phi_h = 1.0_rp / tmp
612  end if
613  end if
614  end select
615 
616  end if
617 
618  ! update N2
619  do k = ks, ke_pbl
620  n2_new(k) = min(atmos_phy_bl_mynn_n2_max, &
621  grav * ( dtldz(k) * betat(k) + dqwdz(k) * betaq(k) ) / potv(k,i,j) )
622  end do
623 
624  do k = ks, ke_pbl
625  ri(k,i,j) = n2_new(k) / dudz2(k,i,j)
626  end do
627 
628  ! length
629  sflx_pt = sflx_sh(i,j) / ( cpdry * exner(ks,i,j) )
630  sflx_ptv = grav / potv(ks,i,j) * ( betat(ks) * sflx_pt + betaq(ks) * sflx_qv(i,j) ) / sfc_dens(i,j)
631  call get_length( &
632  ka, ks, ke_pbl, &
633  i, j, &
634  q(:), n2_new(:), & ! (in)
635  sflx_ptv, rlmo(i,j), & ! (in)
636  cz(:,i,j), & ! (in)
637  fz(:,i,j), & ! (in)
638  l(:,i,j) ) ! (out)
639 
640  call get_q2_level2( &
641  ka, ks, ke_pbl, &
642  dudz2(:,i,j), ri(:,i,j), l(:,i,j), & ! (in)
643  q2_2(:) ) ! (out)
644 
645  do k = ks, ke_pbl
646  ac(k) = min( q(k) / sqrt( q2_2(k) + 1e-20_rp ), 1.0_rp )
647  end do
648 
649  call get_smsh( &
650  ka, ks, ke_pbl, & ! (in)
651  i, j, & ! (in)
652  q(:), ac(:), & ! (in)
653  l(:,i,j), n2_new(:), & ! (in)
654  potv(:,i,j), dudz2(:,i,j), & ! (in)
655  dtldz(:), dqwdz(:), & ! (in)
656  betat(:), betaq(:), & ! (in)
657  mynn_level3, & ! (in)
658  initialize .and. it==1, & ! (in)
659  tsq(:), qsq(:), cov(:), & ! (inout)
660  sm25(:), smp(:), & ! (out)
661  sh25(:), shpgh(:), & ! (out)
662  gammat(:), gammaq(:), & ! (out)
663  f_gamma(:) ) ! (out)
664 
665  do k = ks, ke_pbl
666  lq(k) = l(k,i,j) * q(k)
667  nu_f(k) = lq(k) * sm25(k)
668  kh_f(k) = lq(k) * sh25(k)
669  end do
670  if ( atmos_phy_bl_mynn_similarity ) then
671 ! Nu_f(KS) = KARMAN * z1 * us(i,j) / phi_m
672 ! Kh_f(KS) = KARMAN * z1 * us(i,j) / phi_h
673  end if
674 
675  do k = ks, ke_pbl-1
676  nu(k,i,j) = min( f2h(k,1) * nu_f(k+1) + f2h(k,2) * nu_f(k), &
677  atmos_phy_bl_mynn_nu_max )
678  kh(k,i,j) = min( f2h(k,1) * kh_f(k+1) + f2h(k,2) * kh_f(k), &
679  atmos_phy_bl_mynn_kh_max )
680  end do
681 
682  do k = ks, ke_pbl-1
683  sw = 0.5_rp - sign(0.5_rp, abs(kh(k,i,j)) - eps)
684  pr(k,i,j) = nu(k,i,j) / ( kh(k,i,j) + sw ) * ( 1.0_rp - sw ) &
685  + 1.0_rp * sw
686  end do
687 
688  rho(ks) = dens(ks,i,j) + dt * sflx_qv(i,j) / cdz(ks)
689  do k = ks+1, ke_pbl
690  rho(k) = dens(k,i,j)
691  end do
692 
693  do k = ks, ke_pbl-1
694  rho_h(k) = f2h(k,1) * rho(k+1) + f2h(k,2) * rho(k)
695  rhonu(k) = rho_h(k) * nu(k,i,j)
696  rhokh(k) = rho_h(k) * kh(k,i,j)
697  end do
698 
699 
700  if ( mynn_level3 ) then
701 
702  ! production term calculated by explicit scheme
703  do k = ks, ke_pbl
704  wtl = - lq(k) * ( sh25(k) * dtldz(k) + gammat(k) )
705  wqw = - lq(k) * ( sh25(k) * dqwdz(k) + gammaq(k) )
706  prod_t(k) = - 2.0_rp * wtl * dtldz(k)
707  prod_q(k) = - 2.0_rp * wqw * dqwdz(k)
708  prod_c(k) = - wtl * dqwdz(k) - wqw * dtldz(k)
709  end do
710 
711  call get_gamma_implicit( &
712  ka, ks, ke_pbl, &
713  i, j, &
714  tsq(:), qsq(:), cov(:), & ! (in)
715  dtldz(:), dqwdz(:), potv(:,i,j), & ! (in)
716  prod_t(:), prod_q(:), prod_c(:), & ! (in)
717  betat(:), betaq(:), & ! (in)
718  f_gamma(:), l(:,i,j), q(:), & ! (in)
719  dt, & ! (in)
720  gammat(:), gammaq(:) ) ! (inout)
721 
722  ! update production terms
723  do k = ks, ke_pbl
724  wtl = - lq(k) * ( sh25(k) * dtldz(k) + gammat(k) )
725  wqw = - lq(k) * ( sh25(k) * dqwdz(k) + gammaq(k) )
726  prod_t(k) = - 2.0_rp * wtl * dtldz(k)
727  prod_q(k) = - 2.0_rp * wqw * dqwdz(k)
728  prod_c(k) = - wtl * dqwdz(k) - wqw * dtldz(k)
729  end do
730 
731  end if
732 
733 
734 
735  ! time integration
736 
737  if ( it == nit ) then
738 
739  do k = ks, ke_pbl-1
740  rlqsm_h(k) = rho_h(k) * ( f2h(k,1) * lq(k+1) * smp(k+1) + f2h(k,2) * lq(k) * smp(k) )
741  end do
742 
743  ! dens * u
744 
745  ! countergradient flux
746  do k = ks, ke_pbl-1
747  flx(k) = - rlqsm_h(k) * ( u(k+1,i,j) - u(k,i,j) ) / fdz(k)
748  end do
749 
750  sf_t = sflx_mu(i,j) / cdz(ks)
751  d(ks) = ( u(ks,i,j) * dens(ks,i,j) + dt * ( sf_t - flx(ks) / cdz(ks) ) ) / rho(ks)
752  do k = ks+1, ke_pbl
753  d(k) = u(k,i,j) - dt * ( flx(k) - flx(k-1) ) / ( cdz(k) * rho(k) )
754  end do
755  c(ks) = 0.0_rp
756  do k = ks, ke_pbl-1
757  ap = - dt * rhonu(k) / fdz(k)
758  a(k) = ap / ( rho(k) * cdz(k) )
759  b(k) = - a(k) - c(k) + 1.0_rp
760  c(k+1) = ap / ( rho(k+1) * cdz(k+1) )
761  end do
762  a(ke_pbl) = 0.0_rp
763  b(ke_pbl) = - c(ke_pbl) + 1.0_rp
764 
765  call matrix_solver_tridiagonal( &
766  ka, ks, ke_pbl, &
767  a(:), b(:), c(:), d(:), & ! (in)
768  dummy(:) ) ! (out)
769 ! phi_n(:) ) ! (out)
770 
771  phi_n(ks:ke_pbl) = dummy(ks:ke_pbl)
772  rhou_t(ks,i,j) = ( phi_n(ks) * rho(ks) - u(ks,i,j) * dens(ks,i,j) ) / dt - sf_t
773  do k = ks+1, ke_pbl
774  rhou_t(k,i,j) = ( phi_n(k) - u(k,i,j) ) * rho(k) / dt
775  end do
776  do k = ke_pbl+1, ke
777  rhou_t(k,i,j) = 0.0_rp
778  end do
779  flxu(ks-1,i,j) = 0.0_rp
780  do k = ks, ke_pbl-1
781  flxu(k,i,j) = flx(k) &
782  - rhonu(k) * ( phi_n(k+1) - phi_n(k) ) / fdz(k)
783  end do
784  do k = ke_pbl, ke
785  flxu(k,i,j) = 0.0_rp
786  end do
787 
788 
789  ! dens * v
790 
791  ! countergradient flux
792  do k = ks, ke_pbl-1
793  flx(k) = - rlqsm_h(k) * ( v(k+1,i,j) - v(k,i,j) ) / fdz(k)
794  end do
795 
796  sf_t = sflx_mv(i,j) / cdz(ks)
797  d(ks) = ( v(ks,i,j) * dens(ks,i,j) + dt * ( sf_t - flx(ks) / cdz(ks) ) ) / rho(ks)
798  do k = ks+1, ke_pbl
799  d(k) = v(k,i,j) - dt * ( flx(k) - flx(k-1) ) / ( cdz(k) * rho(k) )
800  end do
801  ! a,b,c is the same as those for the u
802 
803  call matrix_solver_tridiagonal( &
804  ka, ks, ke_pbl, &
805  a(:), b(:), c(:), d(:), & ! (in)
806  phi_n(:) ) ! (out)
807 
808  rhov_t(ks,i,j) = ( phi_n(ks) * rho(ks) - v(ks,i,j) * dens(ks,i,j) ) / dt - sf_t
809  do k = ks+1, ke_pbl
810  rhov_t(k,i,j) = ( phi_n(k) - v(k,i,j) ) * rho(k) / dt
811  end do
812  do k = ke_pbl+1, ke
813  rhov_t(k,i,j) = 0.0_rp
814  end do
815  flxv(ks-1,i,j) = 0.0_rp
816  do k = ks, ke_pbl-1
817  flxv(k,i,j) = flx(k) &
818  - rhonu(k) * ( phi_n(k+1) - phi_n(k) ) / fdz(k)
819  end do
820  do k = ke_pbl, ke
821  flxv(k,i,j) = 0.0_rp
822  end do
823 
824 
825  ! dens * pott
826 
827  ! countergradient flux
828  do k = ks, ke_pbl-1
829  flx(k) = - ( f2h(k,1) * lq(k+1) * gammat(k+1) + f2h(k,2) * lq(k) * gammat(k) ) &
830  * rho_h(k)
831  end do
832 
833  sf_t = sflx_pt / cdz(ks)
834  ! assume that induced vapor has the same PT
835  d(ks) = pott(ks,i,j) + dt * ( sf_t - flx(ks) / cdz(ks) ) / rho(ks)
836  do k = ks+1, ke_pbl
837  d(k) = pott(k,i,j) - dt * ( flx(k) - flx(k-1) ) / ( cdz(k) * rho(k) )
838  end do
839 
840  c(ks) = 0.0_rp
841  do k = ks, ke_pbl-1
842  ap = - dt * rhokh(k) / fdz(k)
843  a(k) = ap / ( rho(k) * cdz(k) )
844  b(k) = - a(k) - c(k) + 1.0_rp
845  c(k+1) = ap / ( rho(k+1) * cdz(k+1) )
846  end do
847  a(ke_pbl) = 0.0_rp
848  b(ke_pbl) = - c(ke_pbl) + 1.0_rp
849 
850  call matrix_solver_tridiagonal( &
851  ka, ks, ke_pbl, &
852  a(:), b(:), c(:), d(:), & ! (in)
853  phi_n(:) ) ! (out)
854 
855  rhot_t(ks,i,j) = ( phi_n(ks) - pott(ks,i,j) ) * rho(ks) / dt - sf_t
856  do k = ks+1, ke_pbl
857  rhot_t(k,i,j) = ( phi_n(k) - pott(k,i,j) ) * rho(k) / dt
858  end do
859  do k = ke_pbl+1, ke
860  rhot_t(k,i,j) = 0.0_rp
861  end do
862  flxt(ks-1,i,j) = 0.0_rp
863  do k = ks, ke_pbl-1
864  flxt(k,i,j) = flx(k) &
865  - rhokh(k) * ( phi_n(k+1) - phi_n(k) ) / cdz(k)
866  end do
867  do k = ke_pbl, ke
868  flxt(k,i,j) = 0.0_rp
869  end do
870 
871  kmin = ks-1
872  if ( flxt(ks,i,j) > 1e-4_rp ) then
873  fmin = flxt(ks,i,j) / rho_h(ks)
874  do k = ks+1, ke_pbl-2
875  tmp = ( flxt(k-1,i,j) + flxt(k,i,j) + flxt(k+1,i,j) ) &
876  / ( rho_h(k-1) + rho_h(k) + rho_h(k+1) ) ! running mean
877  if ( fmin < 0.0_rp .and. tmp > fmin ) exit
878  if ( tmp < fmin ) then
879  fmin = tmp
880  kmin = k
881  end if
882  end do
883  end if
884  zi(i,j) = fz(kmin,i,j) - fz(ks-1,i,j)
885 
886  ! dens * qv
887 
888  ! countergradient flux
889  do k = ks, ke_pbl-1
890  flx(k) = - ( f2h(k,1) * lq(k+1) * gammaq(k+1) + f2h(k,2) * lq(k) * gammaq(k) ) &
891  * rho_h(k)
892  end do
893 
894  sf_t = sflx_qv(i,j) / cdz(ks)
895  d(ks) = ( qv(ks,i,j) * dens(ks,i,j) + dt * ( sf_t - flx(ks) / cdz(ks) ) ) / rho(ks)
896  do k = ks+1, ke_pbl
897  d(k) = qv(k,i,j) - dt * ( flx(k) - flx(k-1) ) / ( cdz(k) * rho(k) )
898  end do
899 
900  call matrix_solver_tridiagonal( &
901  ka, ks, ke_pbl, &
902  a(:), b(:), c(:), d(:), & ! (in)
903  phi_n(:) ) ! (out)
904 
905  rhoqv_t(ks,i,j) = ( phi_n(ks) * rho(ks) - qv(ks,i,j) * dens(ks,i,j) ) / dt - sf_t
906  do k = ks+1, ke_pbl
907  rhoqv_t(k,i,j) = ( phi_n(k) - qv(k,i,j) ) * rho(k) / dt
908  end do
909  do k = ke_pbl+1, ke
910  rhoqv_t(k,i,j) = 0.0_rp
911  end do
912  flxq(ks-1,i,j) = 0.0_rp
913  do k = ks, ke_pbl-1
914  flxq(k,i,j) = flx(k) &
915  - rhokh(k) * ( phi_n(k+1) - phi_n(k) ) / cdz(k)
916  end do
917  do k = ke_pbl, ke
918  flxq(k,i,j) = 0.0_rp
919  end do
920 
921  end if
922 
923 
924  ! dens * TKE
925 
926  ! production
927  nu(ks-1,i,j) = 0.0_rp
928  kh(ks-1,i,j) = 0.0_rp
929  do k = ks, ke_pbl
930  prod(k,i,j) = lq(k) * ( ( sm25(k) + smp(k) ) * dudz2(k,i,j) &
931  - ( sh25(k) * n2_new(k) - shpgh(k) ) )
932  end do
933  if ( atmos_phy_bl_mynn_similarity ) then
934  prod(ks,i,j) = us3 / ( karman * z1 ) * ( phi_m - zeta )
935  end if
936 
937  do k = ks, ke_pbl
938  tke_p(k) = q(k)**2 * 0.5_rp
939  diss(k,i,j) = - 2.0_rp * q(k) / ( b1 * l(k,i,j) )
940 ! prod(k,i,j) = max( prod(k,i,j), - tke_p(k) / dt - diss(k,i,j) * tke_p(k) )
941  end do
942  do k = ke_pbl+1, ke
943  diss(k,i,j) = 0.0_rp
944  prod(k,i,j) = 0.0_rp
945  end do
946 
947  do k = ks, ke_pbl-1
948  d(k) = ( tke_p(k) * dens(k,i,j) + dt * prod(k,i,j) * rho(k) ) / rho(k)
949  end do
950  d(ke_pbl) = 0.0_rp
951 
952  c(ks) = 0.0_rp
953  do k = ks, ke_pbl-1
954  ap = - dt * atmos_phy_bl_mynn_sq_fact * rhonu(k) / fdz(k)
955  a(k) = ap / ( rho(k) * cdz(k) )
956  b(k) = - a(k) - c(k) + 1.0_rp - diss(k,i,j) * dt
957  c(k+1) = ap / ( rho(k+1) * cdz(k+1) )
958  end do
959  a(ke_pbl) = 0.0_rp
960  b(ke_pbl) = - c(ke_pbl) + 1.0_rp - diss(ke_pbl,i,j) * dt
961 
962  call matrix_solver_tridiagonal( &
963  ka, ks, ke_pbl, &
964  a(:), b(:), c(:), d(:), & ! (in)
965  phi_n(:) ) ! (out)
966 
967  do k = ks, ke_pbl-1
968  phi_n(k) = max( phi_n(k), atmos_phy_bl_mynn_tke_min )
969  end do
970  phi_n(ke_pbl) = 0.0_rp
971 
972 
973  if ( it == nit ) then
974  do k = ks, ke_pbl
975  diss(k,i,j) = diss(k,i,j) * phi_n(k)
976  rprog_t(k,i,j,i_tke) = ( phi_n(k) * rho(k) - prog(k,i,j,i_tke) * dens(k,i,j) ) / dt
977  end do
978  do k = ke_pbl+1, ke
979  rprog_t(k,i,j,i_tke) = 0.0_rp
980  end do
981  else
982  do k = ks, ke_pbl
983  q(k) = sqrt( phi_n(k) * 2.0_rp )
984  end do
985  end if
986 
987 
988  if ( .not. mynn_level3 ) cycle
989 
990 !!$ if ( ATMOS_PHY_BL_MYNN_similarity ) then
991 !!$ tmp = 2.0_RP * us(i,j) * phi_h / ( KARMAN * z1 )
992 !!$ tmp = tmp * ( zeta / ( z1 * RLmo(i,j) ) )**2 ! correspoindint to the limitter for zeta
993 !!$ ! TSQ
994 !!$ prod_t(KS) = tmp * ts(i,j)**2
995 !!$ ! QSQ
996 !!$ prod_q(KS) = tmp * ts(i,j) * qs(i,j)
997 !!$ ! COV
998 !!$ prod_c(KS) = tmp * qs(i,j)**2
999 !!$ end if
1000 
1001  ! dens * tsq
1002 
1003  do k = ks, ke_pbl
1004  diss_p(k) = dt * 2.0_rp * q(k) / ( b2 * l(k,i,j) )
1005  end do
1006  do k = ks, ke_pbl-1
1007  d(k) = ( tsq(k) * dens(k,i,j) + dt * prod_t(k) * rho(k) ) / rho(k)
1008  end do
1009  d(ke_pbl) = 0.0_rp
1010  c(ks) = 0.0_rp
1011  do k = ks, ke_pbl-1
1012  ap = - dt * rhonu(k) / fdz(k)
1013  a(k) = ap / ( rho(k) * cdz(k) )
1014  b(k) = - a(k) - c(k) + 1.0_rp + diss_p(k)
1015  c(k+1) = ap / ( rho(k+1) * cdz(k+1) )
1016  end do
1017  a(ke_pbl) = 0.0_rp
1018  b(ke_pbl) = - c(ke_pbl) + 1.0_rp + diss_p(ke_pbl)
1019 
1020  call matrix_solver_tridiagonal( &
1021  ka, ks, ke_pbl, &
1022  a(:), b(:), c(:), d(:), & ! (in)
1023  tsq(:) ) ! (out)
1024 
1025  do k = ks, ke_pbl-1
1026  tsq(k) = max( tsq(k), 0.0_rp )
1027  end do
1028  tsq(ke_pbl) = 0.0_rp
1029 
1030  if ( it == nit ) then
1031  do k = ks, ke_pbl
1032  rprog_t(k,i,j,i_tsq) = ( tsq(k) * rho(k) - prog(k,i,j,i_tsq) * dens(k,i,j) ) / dt
1033  end do
1034  do k = ke_pbl+1, ke
1035  rprog_t(k,i,j,i_tsq) = 0.0_rp
1036  end do
1037  end if
1038 
1039 
1040  ! dens * qsq
1041 
1042  do k = ks, ke_pbl-1
1043  d(k) = ( qsq(k) * dens(k,i,j) + dt * prod_q(k) * rho(k) ) / rho(k)
1044  end do
1045  d(ke_pbl) = 0.0_rp
1046  ! a, b, c are same as those for tsq
1047 
1048  call matrix_solver_tridiagonal( &
1049  ka, ks, ke_pbl, &
1050  a(:), b(:), c(:), d(:), & ! (in)
1051  qsq(:) ) ! (out)
1052 
1053  do k = ks, ke_pbl-1
1054  qsq(k) = max( qsq(k), 0.0_rp )
1055  end do
1056  qsq(ke_pbl) = 0.0_rp
1057 
1058  if ( it == nit ) then
1059  do k = ks, ke_pbl
1060  rprog_t(k,i,j,i_qsq) = ( qsq(k) * rho(k) - prog(k,i,j,i_qsq) * dens(k,i,j) ) / dt
1061  end do
1062  do k = ke_pbl+1, ke
1063  rprog_t(k,i,j,i_qsq) = 0.0_rp
1064  end do
1065  end if
1066 
1067 
1068  ! dens * cov
1069 
1070  do k = ks, ke_pbl-1
1071  d(k) = ( cov(k) * dens(k,i,j) + dt * prod_c(k) * rho(k) ) / rho(k)
1072  end do
1073  d(ke_pbl) = 0.0_rp
1074  ! a, b, c are same as those for tsq
1075 
1076  call matrix_solver_tridiagonal( &
1077  ka, ks, ke_pbl, &
1078  a(:), b(:), c(:), d(:), & ! (in)
1079  cov(:) ) ! (out)
1080 
1081  do k = ks, ke_pbl-1
1082  cov(k) = sign( min( abs(cov(k)), sqrt(tsq(k)*qsq(k)) ), cov(k) )
1083  end do
1084  cov(ke_pbl) = 0.0_rp
1085 
1086  if ( it == nit ) then
1087  do k = ks, ke_pbl
1088  rprog_t(k,i,j,i_cov) = ( cov(k) * rho(k) - prog(k,i,j,i_cov) * dens(k,i,j) ) / dt
1089  end do
1090  do k = ke_pbl+1, ke
1091  rprog_t(k,i,j,i_cov) = 0.0_rp
1092  end do
1093  end if
1094 
1095  end do
1096 
1097  do k = ke_pbl, ke
1098  nu(k,i,j) = 0.0_rp
1099  kh(k,i,j) = 0.0_rp
1100  pr(k,i,j) = 1.0_rp
1101  end do
1102  do k = ke_pbl+1, ke
1103  ri(k,i,j) = undef
1104  prod(k,i,j) = undef
1105  diss(k,i,j) = undef
1106  dudz2(k,i,j) = undef
1107  l(k,i,j) = undef
1108  qlp(k,i,j) = undef
1109  cldfrac(k,i,j) = undef
1110  end do
1111 
1112  end do
1113  end do
1114 
1115 
1116  call file_history_in(ri(:,:,:), 'Ri_MYNN', 'Richardson number', '1', fill_halo=.true. )
1117  call file_history_in(pr(:,:,:), 'Pr_MYNN', 'Prandtl number', '1', fill_halo=.true., dim_type="ZHXY" )
1118  call file_history_in(prod(:,:,:), 'TKE_prod_MYNN', 'TKE production', 'm2/s3', fill_halo=.true.)
1119  call file_history_in(diss(:,:,:), 'TKE_diss_MYNN', 'TKE dissipation', 'm2/s3', fill_halo=.true.)
1120  call file_history_in(dudz2(:,:,:), 'dUdZ2_MYNN', 'dudz2', 'm2/s2', fill_halo=.true.)
1121  call file_history_in(l(:,:,:), 'L_mix_MYNN', 'minxing length', 'm', fill_halo=.true.)
1122 
1123  call file_history_in(flxu(:,:,:), 'ZFLX_RHOU_MYNN', 'Z FLUX of RHOU (MYNN)', 'kg/m/s2', fill_halo=.true., dim_type="ZHXY" )
1124  call file_history_in(flxv(:,:,:), 'ZFLX_RHOV_MYNN', 'Z FLUX of RHOV (MYNN)', 'kg/m/s2', fill_halo=.true., dim_type="ZHXY" )
1125  call file_history_in(flxt(:,:,:), 'ZFLX_RHOT_MYNN', 'Z FLUX of RHOT (MYNN)', 'K kg/m2/s', fill_halo=.true., dim_type="ZHXY" )
1126  call file_history_in(flxq(:,:,:), 'ZFLX_QV_MYNN', 'Z FLUX of RHOQV (MYNN)', 'kg/m2/s', fill_halo=.true., dim_type="ZHXY" )
1127 
1128 
1129  initialize = .false.
1130 
1131  return
1132  end subroutine atmos_phy_bl_mynn_tendency
1133 
1134  !-----------------------------------------------------------------------------
1138  subroutine atmos_phy_bl_mynn_tendency_tracer( &
1139  KA, KS, KE, IA, IS, IE, JA, JS, JE, &
1140  DENS, QTRC, SFLX_Q, &
1141  Kh, MASS, &
1142  CZ, FZ, DDT, &
1143  TRACER_NAME, &
1144  RHOQ_t )
1145  use scale_matrix, only: &
1146  matrix_solver_tridiagonal
1147  use scale_file_history, only: &
1148  file_history_in
1149  integer, intent(in) :: KA, KS, KE
1150  integer, intent(in) :: IA, IS, IE
1151  integer, intent(in) :: JA, JS, JE
1152 
1153  real(RP), intent(in) :: DENS (KA,IA,JA)
1154  real(RP), intent(in) :: QTRC (KA,IA,JA)
1155  real(RP), intent(in) :: SFLX_Q( IA,JA)
1156  real(RP), intent(in) :: Kh (KA,IA,JA)
1157  real(RP), intent(in) :: MASS
1158  real(RP), intent(in) :: CZ ( KA,IA,JA)
1159  real(RP), intent(in) :: FZ (0:KA,IA,JA)
1160  real(DP), intent(in) :: DDT
1161  character(len=*), intent(in) :: TRACER_NAME
1162 
1163  real(RP), intent(out) :: RHOQ_t(KA,IA,JA)
1164 
1165  real(RP) :: QTRC_n(KA)
1166  real(RP) :: RHO (KA)
1167  real(RP) :: RHOKh (KA)
1168  real(RP) :: a(KA)
1169  real(RP) :: b(KA)
1170  real(RP) :: c(KA)
1171  real(RP) :: d(KA)
1172  real(RP) :: rho_h
1173  real(RP) :: ap
1174  real(RP) :: sf_t
1175 
1176  real(RP) :: flx(KA,IA,JA)
1177 
1178  real(RP) :: CDZ(KA)
1179  real(RP) :: FDZ(KA)
1180  real(RP) :: f2h(KA,2)
1181 
1182  real(RP) :: dt
1183 
1184  integer :: KE_PBL
1185  integer :: k, i, j
1186 
1187  dt = real( ddt, kind=rp )
1188 
1189 !OCL INDEPENDENT
1190  !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
1191  !$omp shared(KA,KS,KE,IS,IE,JS,JE) &
1192  !$omp shared(ATMOS_PHY_BL_MYNN_PBL_MAX) &
1193  !$omp shared(RHOQ_t,DENS,QTRC,SFLX_Q,Kh,MASS,CZ,FZ,DT,flx) &
1194  !$omp private(QTRC_n,RHO,RHOKh,rho_h,a,b,c,d,ap,sf_t,CDZ,FDZ,f2h) &
1195  !$omp private(KE_PBL,k,i,j)
1196  do j = js, je
1197  do i = is, ie
1198 
1199  ke_pbl = ks+1
1200  do k = ks+2, ke-1
1201  if ( atmos_phy_bl_mynn_pbl_max >= cz(k,i,j) - fz(ks-1,i,j) ) then
1202  ke_pbl = k
1203  else
1204  exit
1205  end if
1206  end do
1207 
1208  do k = ks, ke_pbl
1209  cdz(k) = fz(k ,i,j) - fz(k-1,i,j)
1210  fdz(k) = cz(k+1,i,j) - cz(k ,i,j)
1211  end do
1212 
1213  call get_f2h( &
1214  ka, ks, ke_pbl-1, &
1215  cdz(:), & ! (in)
1216  f2h(:,:) ) ! (out)
1217 
1218  sf_t = sflx_q(i,j) / cdz(ks)
1219 
1220  rho(ks) = dens(ks,i,j) + dt * sf_t * mass
1221  do k = ks+1, ke_pbl
1222  rho(k) = dens(k,i,j)
1223  end do
1224 
1225  ! dens * coefficient at the half level
1226  do k = ks, ke_pbl-1
1227  rho_h = f2h(k,1) * dens(k+1,i,j) + f2h(k,2) * dens(k,i,j)
1228  rhokh(k) = rho_h * kh(k,i,j)
1229  end do
1230 
1231  d(ks) = ( qtrc(ks,i,j) * dens(ks,i,j) + dt * sf_t ) / rho(ks)
1232  do k = ks+1, ke_pbl
1233  d(k) = qtrc(k,i,j)
1234  end do
1235 
1236  c(ks) = 0.0_rp
1237  do k = ks, ke_pbl-1
1238  ap = - dt * rhokh(k) / fdz(k)
1239  a(k) = ap / ( rho(k) * cdz(k) )
1240  b(k) = - a(k) - c(k) + 1.0_rp
1241  c(k+1) = ap / ( rho(k+1) * cdz(k+1) )
1242  end do
1243  a(ke_pbl) = 0.0_rp
1244  b(ke_pbl) = - c(ke_pbl) + 1.0_rp
1245 
1246  call matrix_solver_tridiagonal( &
1247  ka, ks, ke_pbl, &
1248  a(:), b(:), c(:), d(:), & ! (in)
1249  qtrc_n(:) ) ! (out)
1250 
1251  rhoq_t(ks,i,j) = ( qtrc_n(ks) * rho(ks) - qtrc(ks,i,j) * dens(ks,i,j) ) / dt - sf_t
1252  do k = ks+1, ke_pbl
1253  rhoq_t(k,i,j) = ( qtrc_n(k) - qtrc(k,i,j) ) * rho(k) / dt
1254  end do
1255  do k = ke_pbl+1, ke
1256  rhoq_t(k,i,j) = 0.0_rp
1257  end do
1258 
1259  do k = ks, ke_pbl-1
1260  flx(k,i,j) = - rhokh(k) * ( qtrc_n(k+1) - qtrc_n(k) ) / fdz(k)
1261  end do
1262 
1263  end do
1264  end do
1265 
1266  call file_history_in(flx(:,:,:), 'ZFLX_'//trim(tracer_name)//'_MYNN', 'Z FLUX of DENS * '//trim(tracer_name)//' (MYNN)', 'kg/m2/s', fill_halo=.true.)
1267 
1268  return
1269  end subroutine atmos_phy_bl_mynn_tendency_tracer
1270 
1271  !-----------------------------------------------------------------------------
1272  ! private routines
1273  !-----------------------------------------------------------------------------
1274 
1275 !OCL SERIAL
1276  subroutine get_length( &
1277  KA, KS, KE_PBL, &
1278  i, j, &
1279  q, n2, &
1280  SFLX_PTV, RLmo, &
1281  CZ, &
1282  FZ, &
1283  l )
1284  use scale_const, only: &
1285  grav => const_grav, &
1286  karman => const_karman, &
1287  eps => const_eps
1288  implicit none
1289  integer, intent(in) :: KA, KS, KE_PBL
1290  integer, intent(in) :: i, j ! for debug
1291 
1292  real(RP), intent(in) :: q(KA)
1293  real(RP), intent(in) :: n2(KA)
1294  real(RP), intent(in) :: SFLX_PTV
1295  real(RP), intent(in) :: RLmo
1296  real(RP), intent(in) :: CZ(KA)
1297  real(RP), intent(in) :: FZ(0:KA)
1298 
1299  real(RP), intent(out) :: l(KA)
1300 
1301  real(RP), parameter :: ls_fact_max = 2.0_rp
1302 
1303  real(RP) :: ls
1304  real(RP) :: lb
1305  real(RP) :: lt
1306  real(RP) :: rlt
1307 
1308  real(RP) :: qc
1309  real(RP) :: int_q
1310  real(RP) :: int_qz
1311  real(RP) :: rn2sr
1312  real(RP) :: zeta
1313 
1314  real(RP) :: z
1315  real(RP) :: qdz
1316 
1317  real(RP) :: sw
1318  integer :: k
1319 
1320  int_qz = 0.0_rp
1321  int_q = 0.0_rp
1322  do k = ks, ke_pbl
1323  qdz = q(k) * ( fz(k) - fz(k-1) )
1324  z = cz(k) - fz(ks-1)
1325  int_qz = int_qz + z * qdz
1326  int_q = int_q + qdz
1327  end do
1328  ! LT
1329  lt = min( max(0.23_rp * int_qz / (int_q + 1e-20_rp), &
1330  lt_min), &
1331  atmos_phy_bl_mynn_lt_max )
1332  rlt = 1.0_rp / lt
1333 
1334  qc = ( lt * max(sflx_ptv,0.0_rp) )**oneoverthree ! qc=0 if SFLX_PTV<0
1335 
1336  do k = ks, ke_pbl
1337  z = cz(k) - fz(ks-1)
1338  zeta = min( max( z * rlmo, zeta_min ), zeta_max )
1339 
1340  ! LS
1341  sw = sign(0.5_rp, zeta) + 0.5_rp ! 1 for zeta >= 0, 0 for zeta < 0
1342  ls = karman * z &
1343  * ( sw / (1.0_rp + 2.7_rp*zeta*sw ) &
1344  + min( ( (1.0_rp - 100.0_rp*zeta)*(1.0_rp-sw) )**0.2_rp, ls_fact_max) )
1345 
1346  ! LB
1347  sw = sign(0.5_rp, n2(k)-eps) + 0.5_rp ! 1 for dptdz >0, 0 for dptdz <= 0
1348  rn2sr = 1.0_rp / ( sqrt(n2(k)*sw) + 1.0_rp-sw)
1349  lb = (1.0_rp + 5.0_rp * sqrt(qc*rn2sr*rlt)) * q(k) * rn2sr * sw & ! qc=0 when RLmo > 0
1350  + 1.e10_rp * (1.0_rp-sw)
1351 
1352  ! L
1353  l(k) = 1.0_rp / ( 1.0_rp/ls + rlt + 1.0_rp/(lb+1e-20_rp) )
1354  end do
1355 
1356  return
1357  end subroutine get_length
1358 
1359 !OCL SERIAL
1360  subroutine get_q2_level2( &
1361  KA, KS, KE_PBL, &
1362  dudz2, Ri, l, &
1363  q2_2 )
1364  implicit none
1365  integer, intent(in) :: KA, KS, KE_PBL
1366 
1367  real(RP), intent(in) :: dudz2(KA)
1368  real(RP), intent(in) :: Ri(KA)
1369  real(RP), intent(in) :: l(KA)
1370 
1371  real(RP), intent(out) :: q2_2(KA)
1372 
1373  real(RP) :: rf
1374  real(RP) :: sm_2
1375  real(RP) :: sh_2
1376 
1377  integer :: k
1378 
1379  do k = ks, ke_pbl
1380  rf = min(0.5_rp / af12 * ( ri(k) &
1381  + af12*rf1 &
1382  - sqrt(ri(k)**2 + 2.0_rp*af12*(rf1-2.0_rp*rf2)*ri(k) + (af12*rf1)**2) ), &
1383  rfc)
1384  sh_2 = 3.0_rp * a2 * (g1+g2) * (rfc-rf) / (1.0_rp-rf)
1385  sm_2 = sh_2 * af12 * (rf1-rf) / (rf2-rf)
1386  q2_2(k) = b1 * l(k)**2 * sm_2 * (1.0_rp-rf) * dudz2(k)
1387  end do
1388 
1389  return
1390  end subroutine get_q2_level2
1391 
1392 !OCL SERIAL
1393  subroutine get_smsh( &
1394  KA, KS, KE_PBL, &
1395  i, j, &
1396  q, ac, &
1397  l, n2, &
1398  potv, dudz2, &
1399  dtldz, dqwdz, &
1400  betat, betaq, &
1401  mynn_level3, &
1402  initialize, &
1403  tsq, qsq, cov, &
1404  sm25, smp, &
1405  sh25, shpgh, &
1406  gammat, gammaq, &
1407  f_gamma )
1408  use scale_const, only: &
1409  eps => const_eps, &
1410  grav => const_grav
1411  implicit none
1412  integer, intent(in) :: KA, KS, KE_PBL
1413  integer, intent(in) :: i, j ! for debug
1414 
1415  real(RP), intent(in) :: q(KA)
1416  real(RP), intent(in) :: ac(KA)
1417  real(RP), intent(in) :: l(KA)
1418  real(RP), intent(in) :: n2(KA)
1419  real(RP), intent(in) :: potv(KA)
1420  real(RP), intent(in) :: dudz2(KA)
1421  real(RP), intent(in) :: dtldz(KA)
1422  real(RP), intent(in) :: dqwdz(KA)
1423  real(RP), intent(in) :: betat(KA)
1424  real(RP), intent(in) :: betaq(KA)
1425  logical, intent(in) :: mynn_level3
1426  logical, intent(in) :: initialize
1427 
1428  real(RP), intent(inout) :: tsq(KA)
1429  real(RP), intent(inout) :: qsq(KA)
1430  real(RP), intent(inout) :: cov(KA)
1431 
1432  real(RP), intent(out) :: sm25 (KA) ! S_M2.5
1433  real(RP), intent(out) :: smp (KA) ! S_M'
1434  real(RP), intent(out) :: sh25 (KA) ! S_H2.5
1435  real(RP), intent(out) :: shpgh (KA) ! S_H' G_H * (q/L)^2
1436  real(RP), intent(out) :: gammat (KA) ! \Gamma_\theta
1437  real(RP), intent(out) :: gammaq (KA) ! \Gamma_q
1438  real(RP), intent(out) :: f_gamma(KA) ! - E_H / q^2 * GRAV / Theta_0
1439 
1440  real(RP) :: l2
1441  real(RP) :: q2
1442  real(RP) :: ac2
1443  real(RP) :: p1q2
1444  real(RP) :: p2q2
1445  real(RP) :: p3q2
1446  real(RP) :: p4q2
1447  real(RP) :: p5q2
1448  real(RP) :: rd25
1449  real(RP) :: ghq2
1450  real(RP) :: gmq2
1451  real(RP) :: f1, f2, f3, f4
1452 
1453  ! for level 3
1454  real(RP) :: tvsq
1455  real(RP) :: tvsq25
1456  real(RP) :: tltv
1457  real(RP) :: tltv25
1458  real(RP) :: qwtv
1459  real(RP) :: qwtv25
1460  real(RP) :: tsq25
1461  real(RP) :: qsq25
1462  real(RP) :: cov25
1463  real(RP) :: emq2
1464  real(RP) :: eh
1465  real(RP) :: ew
1466  real(RP) :: rdp
1467  real(RP) :: cw25q2
1468  real(RP) :: fact
1469 
1470  integer :: k
1471 
1472  ! level 2.5
1473  do k = ks, ke_pbl
1474 
1475  ac2 = ac(k)**2
1476  l2 = l(k)**2
1477  q2 = q(k)**2
1478 
1479  f1 = - 3.0_rp * ac2 * a2 * b2 * ( 1.0_rp - c3 )
1480  f2 = - 9.0_rp * ac2 * a1 * a2 * ( 1.0_rp - c2 )
1481  f3 = 9.0_rp * ac2 * a2**2 * ( 1.0_rp - c2 ) * ( 1.0_rp - c5 )
1482  f4 = - 12.0_rp * ac2 * a1 * a2 * ( 1.0_rp - c2 )
1483 
1484  ghq2 = max( - n2(k) * l2, -q2 ) ! L/q <= 1/N for N^2>0
1485 
1486  gmq2 = dudz2(k) * l2
1487 
1488  p1q2 = q2 + f1 * ghq2
1489  p2q2 = q2 + f2 * ghq2
1490  p3q2 = q2 + ( f1 + f3 ) * ghq2
1491  p4q2 = q2 + ( f1 + f4 ) * ghq2
1492  p5q2 = 6.0_rp * ac2 * a1**2 * gmq2
1493  rd25 = q2 / max( p2q2 * p4q2 + p5q2 * p3q2, 1e-20_rp )
1494 
1495  sm25(k) = ac(k) * a1 * ( p3q2 - 3.0_rp * c1 * p4q2 ) * rd25
1496  sh25(k) = ac(k) * a2 * ( p2q2 + 3.0_rp * c1 * p5q2 ) * rd25
1497 
1498 
1499  if ( mynn_level3 ) then ! level 3
1500 
1501  fact = ac(k) * b2 * l2 * sh25(k)
1502  tsq25 = fact * dtldz(k)**2
1503  qsq25 = fact * dqwdz(k)**2
1504  cov25 = fact * dtldz(k) * dqwdz(k)
1505 
1506  ! initialize tsq, qsq, and cov with those of level 2.5
1507  if ( initialize ) then
1508  tsq(k) = tsq25
1509  qsq(k) = qsq25
1510  cov(k) = cov25
1511  end if
1512 
1513  if ( q2 <= 1e-10_rp ) then
1514  smp(k) = 0.0_rp
1515  shpgh(k) = 0.0_rp
1516  f_gamma(k) = 0.0_rp
1517  gammat(k) = 0.0_rp
1518  gammaq(k) = 0.0_rp
1519  else
1520  tltv25 = betat(k) * tsq25 + betaq(k) * cov25
1521  qwtv25 = betat(k) * cov25 + betaq(k) * qsq25
1522  tvsq25 = betat(k) * tltv25 + betaq(k) * qwtv25
1523 
1524  tltv = betat(k) * tsq(k) + betaq(k) * cov(k)
1525  qwtv = betat(k) * cov(k) + betaq(k) * qsq(k)
1526  tvsq = betat(k) * tltv + betaq(k) * qwtv
1527 
1528  rdp = q2 / max( p2q2 * ( f4 * ghq2 + q2 ) + p5q2 * ( f3 * ghq2 + q2 ), 1e-20_rp )
1529 
1530  fact = l2 / q2 * ( grav / potv(k) )**2 * ( tvsq - tvsq25 )
1531 
1532  ew = ( 1.0_rp - c3 ) * ( - p2q2 * f4 - p5q2 * f3 ) * rdp ! (p1-p4)/gh = -f4, (p1-p3)/gh = -f3
1533  if ( abs(ew) > 1e-20_rp ) then
1534  cw25q2 = p1q2 * ( p2q2 + 3.0_rp * c1 * p5q2 ) * rd25 / 3.0_rp
1535  fact = fact * ew
1536  fact = min( max( fact, 0.12_rp * q2 - cw25q2 ), 0.76_rp * q2 - cw25q2 )
1537  fact = fact / ew
1538  end if
1539 
1540  emq2 = 3.0_rp * ac(k) * a1 * ( 1.0_rp - c3 ) * ( f3 - f4 ) * rdp ! (p3-p4)/gh = (f3-f4)
1541  eh = 3.0_rp * ac(k) * a2 * ( 1.0_rp - c3 ) * ( p2q2 + p5q2 ) * rdp
1542 
1543  smp(k) = emq2 * fact
1544  shpgh(k) = eh * fact / l2
1545  f_gamma(k) = - eh * grav / ( q2 * potv(k) )
1546  gammat(k) = ( tltv - tltv25 ) * f_gamma(k)
1547  gammaq(k) = ( qwtv - qwtv25 ) * f_gamma(k)
1548  end if
1549  else ! level 2.5
1550  smp(k) = 0.0_rp
1551  shpgh(k) = 0.0_rp
1552  f_gamma(k) = 0.0_rp
1553  gammat(k) = 0.0_rp
1554  gammaq(k) = 0.0_rp
1555  end if
1556 
1557  end do
1558 
1559  return

References calc_vertical_differece(), scale_const::const_cpdry, scale_const::const_eps, scale_const::const_epstvap, scale_const::const_grav, scale_const::const_karman, get_f2h(), get_gamma_implicit(), scale_atmos_grid_cartesc_index::ie, scale_atmos_grid_cartesc_index::is, 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::mass, partial_condensation(), scale_prc::prc_abort(), scale_precision::rp, and scale_tracer::tracer_name.

Referenced by mod_atmos_phy_bl_driver::atmos_phy_bl_driver_calc_tendency().

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

◆ get_gamma_implicit()

subroutine scale_atmos_phy_bl_mynn::get_gamma_implicit ( integer, intent(in)  KA,
integer, intent(in)  KS,
integer, intent(in)  KE,
integer, intent(in)  i,
integer, intent(in)  j,
real(rp), dimension (ka), intent(in)  tsq,
real(rp), dimension (ka), intent(in)  qsq,
real(rp), dimension (ka), intent(in)  cov,
real(rp), dimension (ka), intent(in)  dtldz,
real(rp), dimension (ka), intent(in)  dqwdz,
real(rp), dimension (ka), intent(in)  POTV,
real(rp), dimension (ka), intent(in)  prod_t,
real(rp), dimension (ka), intent(in)  prod_q,
real(rp), dimension (ka), intent(in)  prod_c,
real(rp), dimension (ka), intent(in)  betat,
real(rp), dimension (ka), intent(in)  betaq,
real(rp), dimension(ka), intent(in)  f_gamma,
real(rp), dimension (ka), intent(in)  l,
real(rp), dimension (ka), intent(in)  q,
real(rp), intent(in)  dt,
real(rp), dimension(ka), intent(inout)  gammat,
real(rp), dimension(ka), intent(inout)  gammaq 
)

Definition at line 1573 of file scale_atmos_phy_bl_mynn.F90.

1573  use scale_const, only: &
1574  grav => const_grav
1575  implicit none
1576  integer, intent(in) :: KA, KS, KE
1577  integer, intent(in) :: i, j ! for debug
1578 
1579  real(RP), intent(in) :: tsq (KA)
1580  real(RP), intent(in) :: qsq (KA)
1581  real(RP), intent(in) :: cov (KA)
1582  real(RP), intent(in) :: dtldz (KA)
1583  real(RP), intent(in) :: dqwdz (KA)
1584  real(RP), intent(in) :: POTV (KA)
1585  real(RP), intent(in) :: prod_t (KA)
1586  real(RP), intent(in) :: prod_q (KA)
1587  real(RP), intent(in) :: prod_c (KA)
1588  real(RP), intent(in) :: betat (KA)
1589  real(RP), intent(in) :: betaq (KA)
1590  real(RP), intent(in) :: f_gamma(KA)
1591  real(RP), intent(in) :: l (KA)
1592  real(RP), intent(in) :: q (KA)
1593  real(RP), intent(in) :: dt
1594 
1595  real(RP), intent(inout) :: gammat(KA)
1596  real(RP), intent(inout) :: gammaq(KA)
1597 
1598  real(RP) :: dtsq
1599  real(RP) :: dqsq
1600  real(RP) :: dcov
1601  real(RP) :: a11, a12, a21, a22, a23, a32, a33
1602  real(RP) :: v1, v2, v3
1603  real(RP) :: f1, f2
1604  real(RP) :: a13, a31, det
1605 
1606  integer :: k
1607 
1608  ! calculate gamma by implicit scheme
1609 
1610  do k = ks, ke
1611 
1612  ! matrix coefficient
1613  f1 = f_gamma(k) * l(k) * q(k)
1614  f2 = - 2.0_rp * q(k) / ( b2 * l(k) )
1615 
1616  v1 = prod_t(k) + f2 * tsq(k)
1617  v2 = prod_c(k) + f2 * cov(k)
1618  v3 = prod_q(k) + f2 * qsq(k)
1619 
1620  a11 = - f1 * dtldz(k) * betat(k) * 2.0_rp + 1.0_rp / dt - f2
1621  a12 = - f1 * dtldz(k) * betaq(k) * 2.0_rp
1622  a21 = - f1 * dqwdz(k) * betat(k)
1623  a22 = - f1 * ( dtldz(k) * betat(k) + dqwdz(k) * betaq(k) ) + 1.0_rp / dt - f2
1624  a23 = - f1 * dtldz(k) * betaq(k)
1625  a32 = - f1 * dqwdz(k) * betat(k) * 2.0_rp
1626  a33 = - f1 * dqwdz(k) * betaq(k) * 2.0_rp + 1.0_rp / dt - f2
1627 
1628  if ( q(k) < 1e-20_rp .or. abs(f_gamma(k)) * dt >= q(k) ) then
1629  ! solve matrix
1630  f1 = a21 / a11
1631  f2 = a23 / a33
1632  dcov = ( v2 - f1 * v1 - f2 * v3 ) / ( a22 - f1 * a12 - f2 * a32 )
1633  dtsq = ( v1 - a12 * dcov ) / a11
1634  dqsq = ( v3 - a32 * dcov ) / a33
1635  else
1636  ! consider change in q
1637  ! dq = - L * G / Theta0 * ( betat * dgammat + betaq * dgammaq )
1638  f1 = - l(k) * grav / potv(k) * f_gamma(k) / q(k)
1639  f2 = f1 * betat(k)**2
1640  a11 = a11 - f2 * v1
1641  a21 = a21 - f2 * v2
1642  a31 = - f2 * v3
1643  f2 = f1 * betat(k) * betaq(k) * 2.0_rp
1644  a12 = a12 - f2 * v1
1645  a22 = a22 - f2 * v2
1646  a32 = a32 - f2 * v3
1647  f2 = f1 * betaq(k)**2
1648  a13 = - f2 * v1
1649  a23 = a23 - f2 * v2
1650  a33 = a33 - f2 * v3
1651  ! solve matrix by Cramer's fomula
1652  det = a11 * ( a22 * a33 - a23 * a32 ) + a12 * ( a23 * a31 - a21 * a33 ) + a13 * ( a21 * a32 - a22 * a31 )
1653  dtsq = ( v1 * ( a22 * a33 - a23 * a32 ) + a12 * ( a23 * v3 - v2 * a33 ) + a13 * ( v2 * a32 - a22 * v3 ) ) / det
1654  dcov = ( a11 * ( v2 * a33 - a23 * v3 ) + v1 * ( a23 * a31 - a21 * a33 ) + a13 * ( a21 * v3 - v2 * a31 ) ) / det
1655  dqsq = ( a11 * ( a22 * v3 - v2 * a32 ) + a12 * ( v2 * a31 - a21 * v3 ) + v3 * ( a21 * a32 - a22 * a31 ) ) / det
1656  end if
1657 
1658  ! modifiy the explicit solution
1659  gammat(k) = gammat(k) + f_gamma(k) * ( betat(k) * dtsq + betaq(k) * dcov )
1660  gammaq(k) = gammaq(k) + f_gamma(k) * ( betat(k) * dcov + betaq(k) * dqsq )
1661  end do
1662 
1663  return

References scale_const::const_grav.

Referenced by atmos_phy_bl_mynn_tendency().

Here is the caller graph for this function:

◆ calc_vertical_differece()

subroutine scale_atmos_phy_bl_mynn::calc_vertical_differece ( integer, intent(in)  KA,
integer, intent(in)  KS,
integer, intent(in)  KE,
integer, intent(in)  i,
integer, intent(in)  j,
real(rp), dimension (ka), intent(in)  U,
real(rp), dimension (ka), intent(in)  V,
real(rp), dimension(ka), intent(in)  POTL,
real(rp), dimension (ka), intent(in)  Qw,
real(rp), dimension (ka), intent(in)  CDZ,
real(rp), dimension (ka), intent(in)  FDZ,
real(rp), dimension (ka,2), intent(in)  F2H,
real(rp), dimension(ka), intent(out)  dudz2,
real(rp), dimension(ka), intent(out)  dtldz,
real(rp), dimension(ka), intent(out)  dqwdz 
)

Definition at line 1673 of file scale_atmos_phy_bl_mynn.F90.

1673  use scale_const, only: &
1674  eps => const_eps
1675  integer, intent(in) :: KA, KS, KE
1676  integer, intent(in) :: i, j ! for debug
1677 
1678  real(RP), intent(in) :: U (KA)
1679  real(RP), intent(in) :: V (KA)
1680  real(RP), intent(in) :: POTL(KA)
1681  real(RP), intent(in) :: Qw (KA)
1682  real(RP), intent(in) :: CDZ (KA)
1683  real(RP), intent(in) :: FDZ (KA)
1684  real(RP), intent(in) :: F2H (KA,2)
1685 
1686  real(RP), intent(out) :: dudz2(KA)
1687  real(RP), intent(out) :: dtldz(KA)
1688  real(RP), intent(out) :: dqwdz(KA)
1689 
1690  real(RP) :: Uh(KA), Vh(KA)
1691  real(RP) :: qh(KA)
1692 
1693  integer :: k
1694 
1695  do k = ks, ke
1696  uh(k) = f2h(k,1) * u(k+1) + f2h(k,2) * u(k)
1697  vh(k) = f2h(k,1) * v(k+1) + f2h(k,2) * v(k)
1698  end do
1699 
1700  dudz2(ks) = ( ( uh(ks) - u(ks) )**2 + ( vh(ks) - v(ks) )**2 ) / ( cdz(ks) * 0.5_rp )**2
1701 ! dudz2(KS) = ( ( Uh(KS) )**2 + ( Vh(KS) )**2 ) / CDZ(KS)**2
1702  dudz2(ks) = max( dudz2(ks), 1e-20_rp )
1703  do k = ks+1, ke
1704  dudz2(k) = ( ( uh(k) - uh(k-1) )**2 + ( vh(k) - vh(k-1) )**2 ) / cdz(k)**2
1705  dudz2(k) = max( dudz2(k), 1e-20_rp )
1706  end do
1707 
1708 
1709  do k = ks, ke
1710  qh(k) = f2h(k,1) * potl(k+1) + f2h(k,2) * potl(k)
1711  end do
1712  dtldz(ks) = ( qh(ks) - potl(ks) ) / ( cdz(ks) * 0.5_rp )
1713 ! dtldz(KS) = ( POTL(KS+1) - POTL(KS) ) / FDZ(KS)
1714  do k = ks+1, ke
1715  dtldz(k) = ( qh(k) - qh(k-1) ) / cdz(k)
1716  end do
1717 
1718  do k = ks, ke
1719  qh(k) = f2h(k,1) * qw(k+1) + f2h(k,2) * qw(k)
1720  end do
1721  dqwdz(ks) = ( qh(ks) - qw(ks) ) / ( cdz(ks) * 0.5_rp )
1722 ! dqwdz(KS) = ( Qw(KS+1) - Qw(KS) ) / FDZ(KS)
1723  do k = ks+1, ke
1724  dqwdz(k) = ( qh(k) - qh(k-1) ) / cdz(k)
1725  end do
1726 
1727  return

References scale_const::const_eps.

Referenced by atmos_phy_bl_mynn_tendency().

Here is the caller graph for this function:

◆ partial_condensation()

subroutine scale_atmos_phy_bl_mynn::partial_condensation ( integer, intent(in)  KA,
integer, intent(in)  KS,
integer, intent(in)  KE,
real(rp), dimension (ka), intent(in)  PRES,
real(rp), dimension (ka), intent(in)  POTT,
real(rp), dimension (ka), intent(in)  POTL,
real(rp), dimension (ka), intent(in)  Qw,
real(rp), dimension (ka), intent(in)  Qdry,
real(rp), dimension(ka), intent(in)  EXNER,
real(rp), dimension(ka), intent(in)  dtldz,
real(rp), dimension(ka), intent(in)  dqwdz,
real(rp), dimension (ka), intent(in)  l,
real(rp), dimension (ka), intent(in)  sh25,
real(rp), dimension (ka), intent(in)  ac,
real(rp), dimension (ka), intent(in)  tsq,
real(rp), dimension (ka), intent(in)  qsq,
real(rp), dimension (ka), intent(in)  cov,
logical, intent(in)  mynn_level3,
real(rp), dimension (ka), intent(out)  betat,
real(rp), dimension (ka), intent(out)  betaq,
real(rp), dimension (ka), intent(out)  Qlp,
real(rp), dimension(ka), intent(out)  cldfrac 
)

Definition at line 1738 of file scale_atmos_phy_bl_mynn.F90.

1738  use scale_const, only: &
1739  cpdry => const_cpdry, &
1740  rvap => const_rvap, &
1741  cpvap => const_cpvap, &
1742  epsvap => const_epsvap, &
1743  epstvap => const_epstvap
1744  use scale_atmos_saturation, only: &
1745  atmos_saturation_pres2qsat => atmos_saturation_pres2qsat_liq
1746  use scale_atmos_hydrometeor, only: &
1747  hydrometeor_lhv => atmos_hydrometeor_lhv
1748  integer, intent(in) :: KA, KS, KE
1749  real(RP), intent(in) :: PRES (KA)
1750  real(RP), intent(in) :: POTT (KA)
1751  real(RP), intent(in) :: POTL (KA)
1752  real(RP), intent(in) :: Qw (KA)
1753  real(RP), intent(in) :: Qdry (KA)
1754  real(RP), intent(in) :: EXNER(KA)
1755  real(RP), intent(in) :: dtldz(KA)
1756  real(RP), intent(in) :: dqwdz(KA)
1757  real(RP), intent(in) :: l (KA)
1758  real(RP), intent(in) :: sh25 (KA)
1759  real(RP), intent(in) :: ac (KA)
1760  real(RP), intent(in) :: tsq (KA)
1761  real(RP), intent(in) :: qsq (KA)
1762  real(RP), intent(in) :: cov (KA)
1763  logical, intent(in) :: mynn_level3
1764 
1765  real(RP), intent(out) :: betat (KA)
1766  real(RP), intent(out) :: betaq (KA)
1767  real(RP), intent(out) :: Qlp (KA)
1768  real(RP), intent(out) :: cldfrac(KA)
1769 
1770  real(RP) :: TEML(KA)
1771  real(RP) :: Qsl (KA)
1772  real(RP) :: LHVL(KA)
1773  real(RP) :: dQsl
1774  real(RP) :: CPtot
1775  real(RP) :: aa, bb, cc
1776  real(RP) :: sigma_s
1777  real(RP) :: Q1
1778  real(RP) :: Rt
1779 
1780  integer :: k
1781 
1782  do k = ks, ke
1783  teml(k) = potl(k) * exner(k)
1784  end do
1785 
1786  call atmos_saturation_pres2qsat( &
1787  ka, ks, ke, &
1788  teml(:), pres(:), qdry(:), & ! (in)
1789  qsl(:) ) ! (out)
1790 
1791  call hydrometeor_lhv( &
1792  ka, ks, ke, & ! (in)
1793  teml(:), & ! (in)
1794  lhvl(:) ) ! (out)
1795 
1796  do k = ks, ke
1797 
1798  dqsl = ( 1.0_rp + qsl(k) / ( qdry(k) * epsvap ) ) * qsl(k) * lhvl(k) / ( rvap * teml(k)**2 )
1799  cptot = qdry(k) * cpdry + qsl(k) * cpvap
1800  aa = 1.0_rp / ( 1.0_rp + lhvl(k)/cptot * dqsl )
1801  bb = exner(k) * dqsl
1802 
1803  if ( mynn_level3 ) then
1804  sigma_s = 0.5_rp * aa * sqrt( max( qsq(k) - 2.0_rp * bb * cov(k) + bb**2 * tsq(k), 1.0e-20_rp ) )
1805  else
1806  ! level 2.5
1807  sigma_s = max( 0.5_rp * aa * l(k) * sqrt( ac(k) * b2 * sh25(k) ) * abs( dqwdz(k) - bb * dtldz(k) ), &
1808  1.0e-10_rp )
1809  end if
1810 
1811  q1 = aa * ( qw(k) - qsl(k) ) * 0.5_rp / sigma_s
1812  cldfrac(k) = min( max( 0.5_rp * ( 1.0_rp + erf(q1*rsqrt_2) ), 0.0_rp ), 1.0_rp )
1813  qlp(k) = min( max( 2.0_rp * sigma_s * ( cldfrac(k) * q1 + rsqrt_2pi &
1814 #if defined(PGI) || defined(SX)
1815  * exp( -min( 0.5_rp*q1**2, 1.e+3_rp ) ) & ! apply exp limiter
1816 #else
1817  * exp(-0.5_rp*q1**2) &
1818 #endif
1819  ), 0.0_rp ), qw(k) * 0.5_rp )
1820  cc = ( 1.0_rp + epstvap * qw(k) - qlp(k) / epsvap ) / exner(k) * lhvl(k) / cptot &
1821  - pott(k) / epsvap
1822  rt = cldfrac(k) - qlp(k) / (2.0_rp*sigma_s*sqrt_2pi) &
1823 #if defined(PGI) || defined(SX)
1824  * exp( -min( q1**2 * 0.5_rp, 1.e+3_rp ) ) ! apply exp limiter
1825 #else
1826  * exp(-q1**2 * 0.5_rp)
1827 #endif
1828  betat(k) = 1.0_rp + epstvap * qw(k) - qlp(k) / epsvap - rt * aa * bb * cc
1829  betaq(k) = epstvap * pott(k) + rt * aa * cc
1830 
1831  end do
1832 
1833  return

References scale_const::const_cpdry, scale_const::const_cpvap, scale_const::const_epstvap, scale_const::const_epsvap, and scale_const::const_rvap.

Referenced by atmos_phy_bl_mynn_tendency().

Here is the caller graph for this function:

◆ get_f2h()

subroutine scale_atmos_phy_bl_mynn::get_f2h ( integer, intent(in)  KA,
integer, intent(in)  KS,
integer, intent(in)  KE,
real(rp), dimension(ka), intent(in)  CDZ,
real(rp), dimension(ka,2), intent(out)  f2h 
)

Definition at line 1841 of file scale_atmos_phy_bl_mynn.F90.

1841  integer, intent(in) :: KA, KS, KE
1842  real(RP), intent(in) :: CDZ(KA)
1843  real(RP), intent(out) :: f2h(KA,2)
1844 
1845  real(RP) :: dz1, dz2
1846  integer :: k
1847 
1848  do k = ks, ke
1849  dz1 = cdz(k+1)
1850  dz2 = cdz(k)
1851  f2h(k,1) = dz2 / ( dz1 + dz2 )
1852  f2h(k,2) = dz1 / ( dz1 + dz2 )
1853  end do
1854 
1855  return

Referenced by atmos_phy_bl_mynn_tendency().

Here is the caller graph for this function:

Variable Documentation

◆ atmos_phy_bl_mynn_ntracer

integer, public scale_atmos_phy_bl_mynn::atmos_phy_bl_mynn_ntracer

Definition at line 50 of file scale_atmos_phy_bl_mynn.F90.

50  integer, public :: ATMOS_PHY_BL_MYNN_NTRACER

Referenced by mod_atmos_phy_bl_driver::atmos_phy_bl_driver_tracer_setup(), and atmos_phy_bl_mynn_tracer_setup().

◆ atmos_phy_bl_mynn_name

character(len=h_short), dimension(:), allocatable, public scale_atmos_phy_bl_mynn::atmos_phy_bl_mynn_name

Definition at line 51 of file scale_atmos_phy_bl_mynn.F90.

51  character(len=H_SHORT), public, allocatable :: ATMOS_PHY_BL_MYNN_NAME(:)

Referenced by mod_atmos_phy_bl_driver::atmos_phy_bl_driver_tracer_setup(), and atmos_phy_bl_mynn_tracer_setup().

◆ atmos_phy_bl_mynn_desc

character(len=h_long), dimension(:), allocatable, public scale_atmos_phy_bl_mynn::atmos_phy_bl_mynn_desc

Definition at line 52 of file scale_atmos_phy_bl_mynn.F90.

52  character(len=H_LONG), public, allocatable :: ATMOS_PHY_BL_MYNN_DESC(:)

Referenced by mod_atmos_phy_bl_driver::atmos_phy_bl_driver_tracer_setup(), and atmos_phy_bl_mynn_tracer_setup().

◆ atmos_phy_bl_mynn_units

character(len=h_short), dimension(:), allocatable, public scale_atmos_phy_bl_mynn::atmos_phy_bl_mynn_units

Definition at line 53 of file scale_atmos_phy_bl_mynn.F90.

53  character(len=H_SHORT), public, allocatable :: ATMOS_PHY_BL_MYNN_UNITS(:)

Referenced by mod_atmos_phy_bl_driver::atmos_phy_bl_driver_tracer_setup(), and atmos_phy_bl_mynn_tracer_setup().

scale_const::const_grav
real(rp), public const_grav
standard acceleration of gravity [m/s2]
Definition: scale_const.F90:46
scale_prc::prc_abort
subroutine, public prc_abort
Abort Process.
Definition: scale_prc.F90:342
scale_const::const_epstvap
real(rp), public const_epstvap
1 / epsilon - 1
Definition: scale_const.F90:70
scale_atmos_phy_bl_mynn::atmos_phy_bl_mynn_tendency
subroutine, public atmos_phy_bl_mynn_tendency(KA, KS, KE, IA, IS, IE, JA, JS, JE, DENS, U, V, POTT, PROG, PRES, EXNER, N2, QDRY, QV, Qw, POTL, POTV, SFC_DENS, SFLX_MU, SFLX_MV, SFLX_SH, SFLX_QV, us, ts, qs, RLmo, CZ, FZ, dt_DP, BULKFLUX_type, RHOU_t, RHOV_t, RHOT_t, RHOQV_t, RPROG_t, Nu, Kh, Qlp, cldfrac, Zi)
ATMOS_PHY_BL_MYNN_tendency calculate tendency by the virtical eddy viscosity.
Definition: scale_atmos_phy_bl_mynn.F90:276
scale_const::const_epsvap
real(rp), public const_epsvap
Rdry / Rvap.
Definition: scale_const.F90:69
scale_const::const_rvap
real(rp), parameter, public const_rvap
specific gas constant (water vapor) [J/kg/K]
Definition: scale_const.F90:63
scale_const::const_eps
real(rp), public const_eps
small number
Definition: scale_const.F90:33
scale_atmos_hydrometeor
module atmosphere / hydrometeor
Definition: scale_atmos_hydrometeor.F90:12
scale_file_history
module file_history
Definition: scale_file_history.F90:15
scale_const::const_cpvap
real(rp), parameter, public const_cpvap
specific heat (water vapor, constant pressure) [J/kg/K]
Definition: scale_const.F90:64
scale_prc
module PROCESS
Definition: scale_prc.F90:11
scale_const::const_pi
real(rp), public const_pi
pi
Definition: scale_const.F90:31
scale_const
module CONSTANT
Definition: scale_const.F90:11
scale_const::const_cpdry
real(rp), public const_cpdry
specific heat (dry air,constant pressure) [J/kg/K]
Definition: scale_const.F90:56
mod_atmos_phy_tb_vars::i_tke
integer, public i_tke
Definition: mod_atmos_phy_tb_vars.F90:61
scale_const::const_karman
real(rp), parameter, public const_karman
von Karman constant
Definition: scale_const.F90:50
scale_matrix
module MATRIX
Definition: scale_matrix.F90:11
scale_atmos_saturation
module atmosphere / saturation
Definition: scale_atmos_saturation.F90:12
Generated by   doxygen 1.8.17