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 (BULKFLUX_type, dx, TKE_MIN, PBL_MAX)
 ATMOS_PHY_BL_MYNN_setup Setup. More...
 
subroutine, public atmos_phy_bl_mynn_finalize
 
subroutine, public atmos_phy_bl_mynn_mkinit (KA, KS, KE, IA, IS, IE, JA, JS, JE, PROG, DENS, U, V, W, POTT, PRES, EXNER, N2, QDRY, QV, Qw, POTL, POTV, SFC_DENS, SFLX_MU, SFLX_MV, SFLX_SH, SFLX_QV, us, ts, qs, RLmo, frac_land, CZ, FZ, F2H, BULKFLUX_type)
 ATMOS_PHY_BL_MYNN_mkinit initialize TKE. More...
 
subroutine, public atmos_phy_bl_mynn_tendency (KA, KS, KE, IA, IS, IE, JA, JS, JE, DENS, U, V, W, POTT, PROG, PRES, EXNER, N2, QDRY, QV, Qw, POTL, POTV, SFC_DENS, SFLX_MU, SFLX_MV, SFLX_SH, SFLX_QV, us, ts, qs, RLmo, frac_land, CZ, FZ, F2H, dt_DP, BULKFLUX_type, RHOU_t, RHOV_t, RHOT_t, RHOQV_t, RPROG_t, Nu, Kh, Qlp, cldfrac, Zi, SFLX_BUOY)
 ATMOS_PHY_BL_MYNN_tendency calculate tendency by the vertical eddy viscosity. More...
 
subroutine calc_vertical_differece (KA, KS, KE, i, j, dudz2, dtldz, dqwdz, U, V, POTL, Qw, QDRY, CDZ, FDZ, F2H, us, ts, qs, phi_m, phi_h, z1, Uh, Vh, qw2, qh)
 
subroutine partial_condensation (KA, KS, KE, betat, betaq, Qlp, cldfrac, PRES, POTT, POTL, Qw, EXNER, tsq, qsq, cov, TEML, LHVL, psat)
 
subroutine get_phi (zeta, phi_m, phi_h, z1, RLmo, I_B_TYPE)
 

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
  • Kitamura, 2010: Modifications to the Mellor-Yamada-Nakanishi-Niino (MYNN) model for the stable stratification case. J. Meteorol. Soc. Japan, 88, 857-864
  • Olson et al., 2019: A description of the MYNN-EDMF scheme and the coupling to other components in WRF-ARW. NOAA Technical Memorandum OAR GSD-61
NAMELIST
  • PARAM_ATMOS_PHY_BL_MYNN
    nametypedefault valuecomment
    ATMOS_PHY_BL_MYNN_K2010 logical .false. > Kitamura (2010)
    ATMOS_PHY_BL_MYNN_O2019 logical .false. > Olson et al. (2019)
    ATMOS_PHY_BL_MYNN_PBL_MAX real(RP) 3000.0_RP > maximum height of the PBL
    ATMOS_PHY_BL_MYNN_N2_MAX real(RP) 1.E1_RP
    ATMOS_PHY_BL_MYNN_NU_MAX real(RP) 1.E4_RP
    ATMOS_PHY_BL_MYNN_KH_MAX real(RP) 1.E4_RP
    ATMOS_PHY_BL_MYNN_LT_MAX real(RP) 2000.0_RP
    ATMOS_PHY_BL_MYNN_USE_ZI logical .true.
    ATMOS_PHY_BL_MYNN_ZETA_LIM real(RP) 10.0_RP
    ATMOS_PHY_BL_MYNN_LEVEL character(len=H_SHORT) "2.5" ! "2.5" or "3"
    ATMOS_PHY_BL_MYNN_SQ_FACT real(RP) 3.0_RP
    ATMOS_PHY_BL_MYNN_CNS real(RP) -1.0_RP N01: 3.5, O19: 10.0
    ATMOS_PHY_BL_MYNN_ALPHA2 real(RP) -1.0_RP N01: 1.0, O19: 0.3
    ATMOS_PHY_BL_MYNN_ALPHA4 real(RP) -1.0_RP N01: 100.0, O19: 10.0
    ATMOS_PHY_BL_MYNN_DUMP_COEF real(RP) 0.5_RP
    ATMOS_PHY_BL_MYNN_DZ_SIM logical .true.
    ATMOS_PHY_BL_MYNN_SIMILARITY logical .true.

History Output
namedescriptionunitvariable
L_mix_MYNN minxing length m L_mix_MYNN
Pr_MYNN Prandtl number 1 Pr_MYNN
Ri_MYNN Richardson number 1 Ri_MYNN
TKE_diss_MYNN TKE dissipation m2/s3 TKE_diss_MYNN
TKE_prod_MYNN TKE production m2/s3 TKE_prod_MYNN
ZFLX_QV2_BL Z FLUX of RHOQV (MYNN) kg/m2/s ZFLX_QV2_BL
ZFLX_QV_BL Z FLUX of RHOQV (MYNN) kg/m2/s ZFLX_QV_BL
ZFLX_RHOT2_BL Z FLUX of RHOT (MYNN) K kg/m2/s ZFLX_RHOT2_BL
ZFLX_RHOT_BL Z FLUX of RHOT (MYNN) K kg/m2/s ZFLX_RHOT_BL
ZFLX_RHOU2_BL Z FLUX of RHOU (MYNN) kg/m/s2 ZFLX_RHOU2_BL
ZFLX_RHOU_BL Z FLUX of RHOU (MYNN) kg/m/s2 ZFLX_RHOU_BL
ZFLX_RHOV2_BL Z FLUX of RHOV (MYNN) kg/m/s2 ZFLX_RHOV2_BL
ZFLX_RHOV_BL Z FLUX of RHOV (MYNN) kg/m/s2 ZFLX_RHOV_BL
dUdZ2_MYNN dudz2 m2/s2 dUdZ2_MYNN

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 175 of file scale_atmos_phy_bl_mynn.F90.

175  use scale_prc, only: &
176  prc_abort
177  implicit none
178 
179  integer :: ierr
180 
181  log_newline
182  log_info("ATMOS_PHY_BL_MYNN_tracer_setup",*) 'Tracer Setup'
183  log_info("ATMOS_PHY_BL_MYNN_tracer_setup",*) 'Mellor-Yamada Nakanishi-Niino scheme'
184 
185  !--- read namelist
186  rewind(io_fid_conf)
187  read(io_fid_conf,nml=param_atmos_phy_bl_mynn,iostat=ierr)
188  if( ierr > 0 ) then !--- fatal error
189  log_error("ATMOS_PHY_BL_MYNN_tracer_setup",*) 'Not appropriate names in namelist PARAM_ATMOS_PHY_BL_MYNN. Check!'
190  call prc_abort
191  endif
192 
193  select case ( atmos_phy_bl_mynn_level )
194  case ( "2.5" )
195  atmos_phy_bl_mynn_ntracer = 1
196  allocate( atmos_phy_bl_mynn_name(1), atmos_phy_bl_mynn_desc(1), atmos_phy_bl_mynn_units(1) )
197  atmos_phy_bl_mynn_name(:) = (/ 'TKE_MYNN' /)
198  atmos_phy_bl_mynn_desc(:) = (/ 'turbulent kinetic energy (MYNN)' /)
199  atmos_phy_bl_mynn_units(:) = (/ 'm2/s2' /)
200  case ( "3" )
201  atmos_phy_bl_mynn_ntracer = 4
202  allocate( atmos_phy_bl_mynn_name(4), atmos_phy_bl_mynn_desc(4), atmos_phy_bl_mynn_units(4) )
203  atmos_phy_bl_mynn_name(:) = (/ 'TKE_MYNN', &
204  'TSQ_MYNN', &
205  'QSQ_MYNN', &
206  'COV_MYNN' /)
207  atmos_phy_bl_mynn_desc(:) = (/ 'turbulent kinetic energy (MYNN) ', &
208  'sub-grid variance of liquid water potential temperature (MYNN) ', &
209  'sub-grid variance of total water content (MYNN) ', &
210  'sub-grid covariance of liquid water potential temperature and total water content (MYNN)' /)
211  atmos_phy_bl_mynn_units(:) = (/ 'm2/s2 ', &
212  'K2 ', &
213  'kg2/kg2', &
214  'K kg ' /)
215  case default
216  log_error("ATMOS_PHY_BL_MYNN_tracer_setup",*) 'only level 2.5 and 3 are supported at this moment'
217  call prc_abort
218  end select
219 
220  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 ( character(len=*), intent(in)  BULKFLUX_type,
real(rp), intent(in), optional  dx,
real(rp), intent(in), optional  TKE_MIN,
real(rp), intent(in), optional  PBL_MAX 
)

ATMOS_PHY_BL_MYNN_setup Setup.

Definition at line 231 of file scale_atmos_phy_bl_mynn.F90.

231  use scale_prc, only: &
232  prc_abort
233  use scale_file_history, only: &
234  file_history_reg
235  implicit none
236  character(len=*), intent(in) :: BULKFLUX_type
237 
238  real(RP), intent(in), optional :: dx
239  real(RP), intent(in), optional :: TKE_MIN
240  real(RP), intent(in), optional :: PBL_MAX
241 
242  integer :: ierr
243  integer :: n
244  !---------------------------------------------------------------------------
245 
246  log_newline
247  log_info("ATMOS_PHY_BL_MYNN_setup",*) 'Setup'
248  log_info("ATMOS_PHY_BL_MYNN_setup",*) 'Mellor-Yamada Nakanishi-Niino scheme'
249 
250  if ( present(tke_min) ) atmos_phy_bl_mynn_tke_min = tke_min
251  if ( present(pbl_max) ) atmos_phy_bl_mynn_pbl_max = pbl_max
252 
253  !--- read namelist
254  rewind(io_fid_conf)
255  read(io_fid_conf,nml=param_atmos_phy_bl_mynn,iostat=ierr)
256  if( ierr < 0 ) then !--- missing
257  log_info("ATMOS_PHY_BL_MYNN_setup",*) 'Not found namelist. Default used.'
258  elseif( ierr > 0 ) then !--- fatal error
259  log_error("ATMOS_PHY_BL_MYNN_setup",*) 'Not appropriate names in namelist PARAM_ATMOS_PHY_BL_MYNN. Check!'
260  call prc_abort
261  endif
262  log_nml(param_atmos_phy_bl_mynn)
263 
264  atmos_phy_bl_mynn_mf = .false.
265  if ( atmos_phy_bl_mynn_o2019 ) then
266  if ( atmos_phy_bl_mynn_cns < 0.0_rp ) atmos_phy_bl_mynn_cns = 3.5_rp
267  if ( atmos_phy_bl_mynn_alpha2 < 0.0_rp ) atmos_phy_bl_mynn_alpha2 = 0.3_rp
268  if ( atmos_phy_bl_mynn_alpha4 < 0.0_rp ) atmos_phy_bl_mynn_alpha4 = 10.0_rp
269  if ( c2 < 0.0_rp ) c2 = 0.729_rp
270  if ( c3 < 0.0_rp ) c3 = 0.34_rp
271  atmos_phy_bl_mynn_k2010 = .true.
272  atmos_phy_bl_mynn_use_zi = .true.
273  if ( atmos_phy_bl_mynn_level .ne. "3" ) then
274  atmos_phy_bl_mynn_mf = .true.
275  end if
276  if ( .not. present(dx) ) then
277  log_error("ATMOS_PHY_BL_MYNN_setup",*) 'dx must be set with O2019'
278  call prc_abort
279  end if
280  if ( atmos_phy_bl_mynn_mf ) then
281  nplume = max_plume
282  do n = 1, max_plume
283  dplume(n) = 100.0_rp * n
284  pw(n) = 0.1_rp + ( 0.5_rp - 0.1_rp ) * ( n - 1 ) / ( max_plume - 1 ) ! not described in O2019
285  if ( dplume(n) > dx ) then
286  nplume = n - 1
287  exit
288  end if
289  end do
290  end if
291 
292  else
293  if ( atmos_phy_bl_mynn_cns < 0.0_rp ) atmos_phy_bl_mynn_cns = 2.7_rp
294  if ( atmos_phy_bl_mynn_alpha2 < 0.0_rp ) atmos_phy_bl_mynn_alpha2 = 1.0_rp
295  if ( atmos_phy_bl_mynn_alpha4 < 0.0_rp ) atmos_phy_bl_mynn_alpha4 = 100.0_rp
296  if ( c2 < 0.0_rp ) c2 = 0.75_rp
297  if ( c3 < 0.0_rp ) c3 = 0.352_rp
298  end if
299 
300  !$acc update device(nplume,dplume,pw,ATMOS_PHY_BL_MYNN_MF)
301 
302 
303  a1 = b1 * (1.0_rp - 3.0_rp * g1) / 6.0_rp
304  a2 = 1.0_rp / (3.0_rp * g1 * b1**(1.0_rp/3.0_rp) * prn )
305  c1 = g1 - 1.0_rp / ( 3.0_rp * a1 * b1**(1.0_rp/3.0_rp) )
306  g2 = ( 2.0_rp * a1 * (3.0_rp - 2.0_rp * c2) + b2 * (1.0_rp - c3) ) / b1
307  f2 = b1 * (g1 + g2) - 3.0_rp * a1 * (1.0_rp - c2)
308 
309  rf2 = b1 * g1 / f2
310  rfc = g1 / (g1 + g2)
311 
312  !$acc update device(A1, A2, C1, C2, C3, G2, F2, Rf2, Rfc)
313 
314  select case ( bulkflux_type )
315  case ( "B91", "B91W01" )
316  ! do nothing
317  case default
318  atmos_phy_bl_mynn_similarity = .false.
319  end select
320 
321  !$acc update device(ATMOS_PHY_BL_MYNN_K2010,ATMOS_PHY_BL_MYNN_O2019,ATMOS_PHY_BL_MYNN_PBL_MAX,ATMOS_PHY_BL_MYNN_TKE_MIN,ATMOS_PHY_BL_MYNN_N2_MAX,ATMOS_PHY_BL_MYNN_NU_MAX,ATMOS_PHY_BL_MYNN_KH_MAX,ATMOS_PHY_BL_MYNN_Lt_MAX,ATMOS_PHY_BL_MYNN_zeta_lim,ATMOS_PHY_BL_MYNN_Sq_fact,ATMOS_PHY_BL_MYNN_cns,ATMOS_PHY_BL_MYNN_alpha2,ATMOS_PHY_BL_MYNN_alpha4,ATMOS_PHY_BL_MYNN_DUMP_coef,ATMOS_PHY_BL_MYNN_dz_sim,ATMOS_PHY_BL_MYNN_similarity)
322 
323 
324  ! history
325  call file_history_reg('Ri_MYNN', 'Richardson number', '1', hist_ri, fill_halo=.true. )
326  call file_history_reg('Pr_MYNN', 'Prandtl number', '1', hist_pr, fill_halo=.true., dim_type="ZHXY" )
327  call file_history_reg('TKE_prod_MYNN', 'TKE production', 'm2/s3', hist_tke_pr, fill_halo=.true.)
328  call file_history_reg('TKE_diss_MYNN', 'TKE dissipation', 'm2/s3', hist_tke_di, fill_halo=.true.)
329  call file_history_reg('dUdZ2_MYNN', 'dudz2', 'm2/s2', hist_dudz2, fill_halo=.true.)
330  call file_history_reg('L_mix_MYNN', 'minxing length', 'm', hist_lmix, fill_halo=.true.)
331 
332  call file_history_reg('ZFLX_RHOU_BL', 'Z FLUX of RHOU (MYNN)', 'kg/m/s2', hist_flxu, fill_halo=.true., dim_type="ZHXY" )
333  call file_history_reg('ZFLX_RHOV_BL', 'Z FLUX of RHOV (MYNN)', 'kg/m/s2', hist_flxv, fill_halo=.true., dim_type="ZHXY" )
334  call file_history_reg('ZFLX_RHOT_BL', 'Z FLUX of RHOT (MYNN)', 'K kg/m2/s', hist_flxt, fill_halo=.true., dim_type="ZHXY" )
335  call file_history_reg('ZFLX_QV_BL', 'Z FLUX of RHOQV (MYNN)', 'kg/m2/s', hist_flxq, fill_halo=.true., dim_type="ZHXY" )
336 
337  call file_history_reg('ZFLX_RHOU2_BL', 'Z FLUX of RHOU (MYNN)', 'kg/m/s2', hist_flxu2, fill_halo=.true., dim_type="ZHXY" )
338  call file_history_reg('ZFLX_RHOV2_BL', 'Z FLUX of RHOV (MYNN)', 'kg/m/s2', hist_flxv2, fill_halo=.true., dim_type="ZHXY" )
339  call file_history_reg('ZFLX_RHOT2_BL', 'Z FLUX of RHOT (MYNN)', 'K kg/m2/s', hist_flxt2, fill_halo=.true., dim_type="ZHXY" )
340  call file_history_reg('ZFLX_QV2_BL', 'Z FLUX of RHOQV (MYNN)', 'kg/m2/s', hist_flxq2, fill_halo=.true., dim_type="ZHXY" )
341 
342  return

References scale_file_history::file_history_reg(), 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_finalize()

subroutine, public scale_atmos_phy_bl_mynn::atmos_phy_bl_mynn_finalize

Definition at line 349 of file scale_atmos_phy_bl_mynn.F90.

349  implicit none
350 
351  deallocate( atmos_phy_bl_mynn_name, atmos_phy_bl_mynn_desc, atmos_phy_bl_mynn_units )
352 
353  return

References atmos_phy_bl_mynn_desc, atmos_phy_bl_mynn_name, and atmos_phy_bl_mynn_units.

Referenced by mod_atmos_phy_bl_driver::atmos_phy_bl_driver_finalize().

Here is the caller graph for this function:

◆ atmos_phy_bl_mynn_mkinit()

subroutine, public scale_atmos_phy_bl_mynn::atmos_phy_bl_mynn_mkinit ( 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(ka,ia,ja,atmos_phy_bl_mynn_ntracer), intent(out)  PROG,
real(rp), dimension (ka,ia,ja), intent(in)  DENS,
real(rp), dimension (ka,ia,ja), intent(in)  U,
real(rp), dimension (ka,ia,ja), intent(in)  V,
real(rp), dimension (ka,ia,ja), intent(in)  W,
real(rp), dimension (ka,ia,ja), intent(in)  POTT,
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(in)  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(ia,ja), intent(in)  frac_land,
real(rp), dimension( ka,ia,ja), intent(in)  CZ,
real(rp), dimension(0:ka,ia,ja), intent(in)  FZ,
real(rp), dimension(ka,2,ia,ja), intent(in)  F2H,
character(len=*), intent(in)  BULKFLUX_type 
)

ATMOS_PHY_BL_MYNN_mkinit initialize TKE.

Definition at line 371 of file scale_atmos_phy_bl_mynn.F90.

371  use scale_const, only: &
372  cpdry => const_cpdry
373  use scale_prc, only: &
374  prc_abort
375  use scale_atmos_hydrometeor, only: &
376  hydrometeor_lhv => atmos_hydrometeor_lhv, &
377  cp_vapor
378  use scale_matrix, only: &
379  matrix_solver_tridiagonal
380  use scale_file_history, only: &
381  file_history_query, &
382  file_history_put
383  implicit none
384 
385  integer, intent(in) :: KA, KS, KE
386  integer, intent(in) :: IA, IS, IE
387  integer, intent(in) :: JA, JS, JE
388 
389  real(RP), intent(out) :: PROG(KA,IA,JA,ATMOS_PHY_BL_MYNN_ntracer)
390 
391  real(RP), intent(in) :: DENS (KA,IA,JA)
392  real(RP), intent(in) :: U (KA,IA,JA)
393  real(RP), intent(in) :: V (KA,IA,JA)
394  real(RP), intent(in) :: W (KA,IA,JA)
395  real(RP), intent(in) :: POTT (KA,IA,JA)
396  real(RP), intent(in) :: PRES (KA,IA,JA)
397  real(RP), intent(in) :: EXNER (KA,IA,JA)
398  real(RP), intent(in) :: N2 (KA,IA,JA)
399  real(RP), intent(in) :: QDRY (KA,IA,JA)
400  real(RP), intent(in) :: QV (KA,IA,JA)
401  real(RP), intent(in) :: Qw (KA,IA,JA)
402  real(RP), intent(in) :: POTL (KA,IA,JA)
403  real(RP), intent(in) :: POTV (KA,IA,JA)
404  real(RP), intent(in) :: SFC_DENS( IA,JA)
405  real(RP), intent(in) :: SFLX_MU ( IA,JA)
406  real(RP), intent(in) :: SFLX_MV ( IA,JA)
407  real(RP), intent(in) :: SFLX_SH ( IA,JA)
408  real(RP), intent(in) :: SFLX_QV ( IA,JA)
409  real(RP), intent(in) :: us ( IA,JA)
410  real(RP), intent(in) :: ts ( IA,JA)
411  real(RP), intent(in) :: qs ( IA,JA)
412  real(RP), intent(in) :: RLmo ( IA,JA)
413 
414  real(RP), intent(in) :: frac_land(IA,JA)
415 
416  real(RP), intent(in) :: CZ( KA,IA,JA)
417  real(RP), intent(in) :: FZ(0:KA,IA,JA)
418  real(RP), intent(in) :: F2H(KA,2,IA,JA)
419 
420  character(len=*), intent(in) :: BULKFLUX_type
421 
422  ! bulkflux
423  integer :: I_B_TYPE
424 
425  real(RP) :: Nu (KA)
426  real(RP) :: Kh (KA)
427  real(RP) :: Qlp (KA)
428  real(RP) :: cldfrac(KA)
429  real(RP) :: Zi
430  real(RP) :: SFLX_BUOY
431 
432  real(RP) :: Ri (KA)
433  real(RP) :: Pr (KA)
434  real(RP) :: prod (KA)
435  real(RP) :: diss (KA)
436  real(RP) :: dudz2(KA)
437  real(RP) :: l (KA)
438  real(RP) :: rho_h(KA)
439 
440  real(RP) :: RHO (KA)
441  real(RP) :: RHONu (KA)
442  real(RP) :: RHOKh (KA)
443  real(RP) :: N2_new(KA)
444  real(RP) :: SFLX_PT
445  real(RP) :: sm25 (KA)
446  real(RP) :: sh25 (KA)
447  real(RP) :: q (KA)
448  real(RP) :: lq (KA)
449 
450  real(RP) :: tke(KA)
451 
452  ! for level 3
453  real(RP) :: tsq (KA)
454  real(RP) :: qsq (KA)
455  real(RP) :: cov (KA)
456  real(RP) :: prod_t(KA)
457  real(RP) :: prod_q(KA)
458  real(RP) :: prod_c(KA)
459  real(RP) :: diss_p(KA)
460  real(RP) :: dtldz(KA)
461  real(RP) :: dqwdz(KA)
462  real(RP) :: smp (KA)
463  real(RP) :: shpgh(KA)
464  real(RP) :: gammat (KA)
465  real(RP) :: gammaq (KA)
466 
467  real(RP) :: CPtot
468 
469  real(RP) :: FDZ(KA)
470  real(RP) :: CDZ(KA)
471  real(RP) :: z (KA)
472 
473  logical :: mynn_level3
474 
475  real(RP), parameter :: dt = 1.0_rp
476 
477  ! for O2019
478  real(RP) :: tflux(0:KA)
479  real(RP) :: qflux(0:KA)
480  real(RP) :: uflux(0:KA)
481  real(RP) :: vflux(0:KA)
482  real(RP) :: eflux(0:KA)
483 
484  ! work
485  real(RP) :: PTLV (KA)
486  real(RP) :: Nu_f (KA)
487  real(RP) :: Kh_f (KA)
488  real(RP) :: q2_2 (KA)
489  real(RP) :: ac (KA)
490  real(RP) :: betat(KA)
491  real(RP) :: betaq(KA)
492 
493  real(RP) :: flx(0:KA)
494  real(RP) :: a(KA)
495  real(RP) :: b(KA)
496  real(RP) :: c(KA)
497  real(RP) :: d(KA)
498 
499  real(RP) :: f_smp (KA)
500  real(RP) :: f_shpgh(KA)
501  real(RP) :: f_gamma(KA)
502  real(RP) :: tltv25 (KA)
503  real(RP) :: qwtv25 (KA)
504  real(RP) :: tvsq25 (KA)
505  real(RP) :: tvsq_up(KA)
506  real(RP) :: tvsq_lo(KA)
507  real(RP) :: dtsq (KA)
508  real(RP) :: dqsq (KA)
509  real(RP) :: dcov (KA)
510 
511  real(RP) :: mflux(0:KA)
512 
513  real(RP) :: Uh(KA), Vh(KA), Wh(KA)
514  real(RP) :: qw2(KA), qh(KA)
515  real(RP) :: tlh(KA), tvh(KA)
516  real(RP) :: eh(KA), dh(KA)
517  real(RP) :: TEML(KA), LHVL(KA), psat(KA)
518 
519 #ifdef _OPENACC
520  real(RP) :: work(KA,4)
521 #endif
522 
523  integer :: KE_PBL
524  integer :: k, i, j
525 
526  !---------------------------------------------------------------------------
527 
528  mynn_level3 = ( atmos_phy_bl_mynn_level == "3" )
529 
530  select case ( bulkflux_type )
531  case ( 'B71' )
532  i_b_type = i_b71
533  case ( 'B91' )
534  i_b_type = i_b91
535  case ( 'B91W01' )
536  i_b_type = i_b91w01
537  case default
538  log_error("ATMOS_PHY_BL_MYNN_init_TKE",*) "BULKFLUX_type is invalid: ", trim(bulkflux_type)
539  call prc_abort
540  end select
541 
542  !$acc data copyout(PROG) copyin(DENS,U,V,POTT,PRES,EXNER,N2,QDRY,QV,Qw,POTL,POTV,SFC_DENS,SFLX_MU,SFLX_MV,SFLX_SH,SFLX_QV,us,ts,qs,RLmo,CZ,FZ,F2H)
543 
544 
545 !OCL INDEPENDENT
546  !$omp parallel do default(none) &
547  !$omp OMP_SCHEDULE_ collapse(2) &
548  !$omp shared(KA,KS,KE,IS,IE,JS,JE, &
549  !$omp CPdry,CP_VAPOR,UNDEF, &
550  !$omp ATMOS_PHY_BL_MYNN_N2_MAX,ATMOS_PHY_BL_MYNN_TKE_MIN, &
551  !$omp ATMOS_PHY_BL_MYNN_NU_MAX,ATMOS_PHY_BL_MYNN_KH_MAX, &
552  !$omp ATMOS_PHY_BL_MYNN_Sq_fact, ATMOS_PHY_BL_MYNN_zeta_lim, &
553  !$omp ATMOS_PHY_BL_MYNN_similarity,ATMOS_PHY_BL_MYNN_dz_sim, &
554  !$omp ATMOS_PHY_BL_MYNN_DUMP_coef, &
555  !$omp ATMOS_PHY_BL_MYNN_PBL_MAX, &
556  !$omp ATMOS_PHY_BL_MYNN_K2010,ATMOS_PHY_BL_MYNN_O2019,ATMOS_PHY_BL_MYNN_MF, &
557  !$omp DENS,PROG,U,V,W,POTT,PRES,QDRY,QV,Qw,POTV,POTL,EXNER,N2, &
558  !$omp SFC_DENS,SFLX_MU,SFLX_MV,SFLX_SH,SFLX_QV,us,ts,qs,RLmo, &
559  !$omp mynn_level3, &
560  !$omp frac_land, &
561  !$omp CZ,FZ,F2H, &
562  !$omp I_B_TYPE) &
563  !$omp private(N2_new,lq,sm25,sh25,q, &
564  !$omp SFLX_PT,CPtot,RHO,RHONu,RHOKh, &
565  !$omp Nu,Kh,Qlp,cldfrac,Zi,SFLX_BUOY, &
566  !$omp Ri,Pr,prod,diss,dudz2,l, &
567  !$omp smp,shpgh, &
568  !$omp dtldz,dqwdz,gammat,gammaq, &
569  !$omp rho_h,tke,CDZ,FDZ,z, &
570  !$omp tsq,qsq,cov, &
571  !$omp prod_t,prod_q,prod_c,diss_p, &
572  !$omp tflux,qflux,uflux,vflux,eflux, &
573  !$omp PTLV, Nu_f, Kh_f, q2_2, ac, betat, betaq, &
574  !$omp flx, a, b, c, d, &
575  !$omp f_smp, f_shpgh, f_gamma, tltv25, qwtv25, tvsq25, &
576  !$omp tvsq_up, tvsq_lo, dtsq, dqsq, dcov, &
577  !$omp mflux, &
578  !$omp Uh, Vh, Wh, qw2, qh, tlh, tvh, eh, dh, &
579  !$omp TEML, LHVL, psat, &
580  !$omp KE_PBL,k,i,j)
581  !$acc kernels
582  !$acc loop independent collapse(2) &
583  !$acc private(N2_new,lq,sm25,sh25,q, &
584  !$acc SFLX_PT,CPtot,RHO,RHONu,RHOKh, &
585  !$acc Nu,Kh,Qlp,cldfrac,Zi,SFLX_BUOY, &
586  !$acc Ri,Pr,prod,diss,dudz2,l, &
587  !$acc smp,shpgh, &
588  !$acc dtldz,dqwdz,gammat,gammaq, &
589  !$acc rho_h,tke,CDZ,FDZ,z, &
590  !$acc tsq,qsq,cov, &
591  !$acc prod_t,prod_q,prod_c,diss_p, &
592  !$acc tflux,qflux,uflux,vflux,eflux, &
593  !$acc work, &
594  !$acc PTLV, Nu_f, Kh_f, q2_2, ac, betat, betaq, &
595  !$acc flx, a, b, c, d, &
596  !$acc f_smp, f_shpgh, f_gamma, tltv25, qwtv25, tvsq25, &
597  !$acc tvsq_up, tvsq_lo, dtsq, dqsq, dcov, &
598  !$acc mflux, &
599  !$acc Uh, Vh, Wh, qw2, qh, tlh, tvh, eh, dh, &
600  !$acc TEML, LHVL, psat, &
601  !$acc KE_PBL,k,i,j)
602  do j = js, je
603  do i = is, ie
604 
605  do k = ks, ke-1
606  z(k) = cz(k,i,j) - fz(ks-1,i,j)
607  end do
608 
609  ke_pbl = ks+1
610  !$acc loop seq
611  do k = ks+2, ke-1
612  if ( atmos_phy_bl_mynn_pbl_max >= z(k) ) then
613  ke_pbl = k
614  else
615  exit
616  end if
617  end do
618 
619  do k = ks, ke_pbl+1
620  cdz(k) = fz(k,i,j) - fz(k-1,i,j)
621  end do
622  do k = ks, ke_pbl
623  fdz(k) = cz(k+1,i,j) - cz(k,i,j)
624  end do
625 
626  do k = ks, ke_pbl
627  tke(k) = 0.01_rp
628  end do
629 
630  cptot = cpdry + sflx_qv(i,j) * ( cp_vapor - cpdry )
631  sflx_pt = sflx_sh(i,j) / ( cptot * exner(ks,i,j) )
632 
633  call mynn_main( ka, ks, ke_pbl, &
634  i, j, &
635  tke(:), tsq(:), qsq(:), cov(:), & ! (inout)
636  q(:), l(:), lq(:), & ! (out)
637  nu(:), rhonu(:), kh(:), rhokh(:), pr(:), & ! (out)
638  prod(:), prod_t(:), prod_q(:), prod_c(:), & ! (out)
639  diss(:), diss_p(:), & ! (out)
640  sm25(:), smp(:), sh25(:), shpgh(:), & ! (out)
641  gammat(:), gammaq(:), & ! (out)
642  dudz2(:), n2_new(:), ri(:), & ! (out)
643  dtldz(:), dqwdz(:), & ! (out)
644  rho(:), rho_h(:), & ! (out)
645  uflux(:), vflux(:), tflux(:), qflux(:), eflux(:), & ! (out)
646  qlp(:), cldfrac(:), zi, sflx_buoy, & ! (out)
647  u(:,i,j), v(:,i,j), w(:,i,j), & ! (in)
648  dens(:,i,j), pres(:,i,j), & ! (in)
649  pott(:,i,j), potl(:,i,j), potv(:,i,j), & ! (in)
650  qw(:,i,j), n2(:,i,j), & ! (in)
651  exner(:,i,j), qdry(:,i,j), & ! (in)
652  sflx_pt, sflx_sh(i,j), sflx_qv(i,j), sfc_dens(i,j), & ! (in)
653  rlmo(i,j), us(i,j), ts(i,j), qs(i,j), & ! (in)
654  z(:), cdz(:), fdz(:), f2h(:,:,i,j), & ! (in)
655  frac_land(i,j), & ! (in)
656  dt, & ! (in)
657  i_b_type, & ! (in)
658  mynn_level3, .true., & ! (in)
659 #ifdef _OPENACC
660  work(:,:), & ! (work)
661 #endif
662  ptlv(:), nu_f(:), kh_f(:), q2_2(:), ac(:), & ! (work)
663  betat(:), betaq(:), & ! (work)
664  flx(:), a(:), b(:), c(:), d(:), & ! (work)
665  f_smp(:), f_shpgh(:), f_gamma(:), & ! (work)
666  tltv25(:), qwtv25(:), tvsq25(:), & ! (work)
667  tvsq_up(:), tvsq_lo(:), dtsq(:), dqsq(:), dcov(:), & ! (work)
668  mflux(:), & ! (work)
669  uh(:), vh(:), wh(:), qw2(:), qh(:), & ! (work)
670  tlh(:), tvh(:), eh(:), dh(:), & ! (work)
671  teml(:), lhvl(:), psat(:) ) ! (work)
672 
673  do k = ks, ke_pbl
674  prog(k,i,j,i_tke) = tke(k)
675  end do
676  do k = ke_pbl+1, ke
677  prog(k,i,j,i_tke) = 0.0_rp
678  end do
679 
680  if ( mynn_level3 ) then
681  do k = ks, ke_pbl
682  prog(k,i,j,i_tsq) = tsq(k)
683  prog(k,i,j,i_qsq) = qsq(k)
684  prog(k,i,j,i_cov) = cov(k)
685  end do
686  do k = ke_pbl+1, ke
687  prog(k,i,j,i_tsq) = 0.0_rp
688  prog(k,i,j,i_qsq) = 0.0_rp
689  prog(k,i,j,i_cov) = 0.0_rp
690  end do
691  end if
692 
693  end do
694  end do
695  !$acc end kernels
696 
697  !$acc end data
698 
699  return

References scale_const::const_cpdry, scale_atmos_hydrometeor::cp_vapor, 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, and scale_prc::prc_abort().

Referenced by mod_atmos_phy_bl_driver::atmos_phy_bl_driver_mkinit().

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 (ka,ia,ja), intent(out)  DENS,
real(rp), dimension (counter gradient or mass flux), intent(out)  U,
real(rp), dimension (counter gradient or mass flux), intent(out)  V,
real(rp), dimension (ka,ia,ja), intent(out)  W,
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 (counter gradient or mass flux), 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(ia,ja), intent(in)  frac_land,
real(rp), dimension( ka,ia,ja), intent(in)  CZ,
real(rp), dimension(0:ka,ia,ja), intent(in)  FZ,
real(rp), dimension(ka,2,ia,ja), intent(in)  F2H,
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,
real(rp), dimension (ia,ja), intent(out)  SFLX_BUOY 
)

ATMOS_PHY_BL_MYNN_tendency calculate tendency by the vertical eddy viscosity.

Definition at line 720 of file scale_atmos_phy_bl_mynn.F90.

720  use scale_const, only: &
721  cpdry => const_cpdry
722  use scale_prc, only: &
723  prc_abort
724  use scale_atmos_hydrometeor, only: &
725  hydrometeor_lhv => atmos_hydrometeor_lhv, &
726  cp_vapor
727  use scale_matrix, only: &
728  matrix_solver_tridiagonal
729  use scale_file_history, only: &
730  file_history_query, &
731  file_history_put
732  implicit none
733 
734  integer, intent(in) :: KA, KS, KE
735  integer, intent(in) :: IA, IS, IE
736  integer, intent(in) :: JA, JS, JE
737 
738  real(RP), intent(in) :: DENS (KA,IA,JA)
739  real(RP), intent(in) :: U (KA,IA,JA)
740  real(RP), intent(in) :: V (KA,IA,JA)
741  real(RP), intent(in) :: W (KA,IA,JA)
742  real(RP), intent(in) :: POTT (KA,IA,JA)
743  real(RP), intent(in) :: PROG (KA,IA,JA,ATMOS_PHY_BL_MYNN_ntracer)
744  real(RP), intent(in) :: PRES (KA,IA,JA)
745  real(RP), intent(in) :: EXNER (KA,IA,JA)
746  real(RP), intent(in) :: N2 (KA,IA,JA)
747  real(RP), intent(in) :: QDRY (KA,IA,JA)
748  real(RP), intent(in) :: QV (KA,IA,JA)
749  real(RP), intent(in) :: Qw (KA,IA,JA)
750  real(RP), intent(in) :: POTL (KA,IA,JA)
751  real(RP), intent(in) :: POTV (KA,IA,JA)
752  real(RP), intent(in) :: SFC_DENS( IA,JA)
753  real(RP), intent(in) :: SFLX_MU ( IA,JA)
754  real(RP), intent(in) :: SFLX_MV ( IA,JA)
755  real(RP), intent(in) :: SFLX_SH ( IA,JA)
756  real(RP), intent(in) :: SFLX_QV ( IA,JA)
757  real(RP), intent(in) :: us ( IA,JA)
758  real(RP), intent(in) :: ts ( IA,JA)
759  real(RP), intent(in) :: qs ( IA,JA)
760  real(RP), intent(in) :: RLmo ( IA,JA)
761 
762  real(RP), intent(in) :: frac_land(IA,JA)
763 
764  real(RP), intent(in) :: CZ( KA,IA,JA)
765  real(RP), intent(in) :: FZ(0:KA,IA,JA)
766  real(RP), intent(in) :: F2H(KA,2,IA,JA)
767  real(DP), intent(in) :: dt_DP
768 
769  character(len=*), intent(in) :: BULKFLUX_type
770 
771  real(RP), intent(out) :: RHOU_t (KA,IA,JA)
772  real(RP), intent(out) :: RHOV_t (KA,IA,JA)
773  real(RP), intent(out) :: RHOT_t (KA,IA,JA)
774  real(RP), intent(out) :: RHOQV_t(KA,IA,JA)
775  real(RP), intent(out) :: RPROG_t(KA,IA,JA,ATMOS_PHY_BL_MYNN_ntracer)
776  real(RP), intent(out) :: Nu (KA,IA,JA)
777  real(RP), intent(out) :: Kh (KA,IA,JA)
778  real(RP), intent(out) :: Qlp (KA,IA,JA)
779  real(RP), intent(out) :: cldfrac(KA,IA,JA)
780  real(RP), intent(out) :: Zi (IA,JA)
781  real(RP), intent(out) :: SFLX_BUOY (IA,JA)
782 
783  ! bulkflux
784  integer :: I_B_TYPE
785 
786  real(RP) :: Ri (KA,IA,JA)
787  real(RP) :: Pr (KA,IA,JA)
788  real(RP) :: prod (KA,IA,JA)
789  real(RP) :: diss (KA,IA,JA)
790  real(RP) :: dudz2(KA,IA,JA)
791  real(RP) :: l (KA,IA,JA)
792  real(RP) :: rho_h(KA)
793 
794  real(RP) :: flxU(0:KA,IA,JA)
795  real(RP) :: flxV(0:KA,IA,JA)
796  real(RP) :: flxT(0:KA,IA,JA)
797  real(RP) :: flxQ(0:KA,IA,JA)
798 
799  real(RP) :: flxU2(0:KA,IA,JA)
800  real(RP) :: flxV2(0:KA,IA,JA)
801  real(RP) :: flxT2(0:KA,IA,JA)
802  real(RP) :: flxQ2(0:KA,IA,JA)
803 
804  real(RP) :: RHO (KA)
805  real(RP) :: RHONu (KA)
806  real(RP) :: RHOKh (KA)
807  real(RP) :: N2_new(KA)
808  real(RP) :: SFLX_PT
809  real(RP) :: sm25 (KA)
810  real(RP) :: sh25 (KA)
811  real(RP) :: q (KA)
812  real(RP) :: lq (KA)
813 
814  real(RP) :: tke(KA)
815 
816  ! for level 3
817  real(RP) :: tsq (KA)
818  real(RP) :: qsq (KA)
819  real(RP) :: cov (KA)
820  real(RP) :: prod_t(KA)
821  real(RP) :: prod_q(KA)
822  real(RP) :: prod_c(KA)
823  real(RP) :: diss_p(KA)
824  real(RP) :: dtldz(KA)
825  real(RP) :: dqwdz(KA)
826  real(RP) :: smp (KA)
827  real(RP) :: shpgh(KA)
828  real(RP) :: gammat (KA)
829  real(RP) :: gammaq (KA)
830  real(RP) :: rlqsm_h(KA)
831 
832  real(RP) :: flx(0:KA)
833  real(RP) :: a(KA)
834  real(RP) :: b(KA)
835  real(RP) :: c(KA)
836  real(RP) :: d(KA)
837  real(RP) :: ap
838  real(RP) :: phi_N(KA)
839  real(RP) :: dummy(KA)
840 
841  real(RP) :: CPtot
842 
843  real(RP) :: sf_t
844 
845  real(RP) :: FDZ(KA)
846  real(RP) :: CDZ(KA)
847  real(RP) :: z (KA)
848 
849  logical :: mynn_level3
850 
851  real(RP) :: dt
852 
853  real(RP) :: fmin
854  integer :: kmin
855 
856  ! for O2019
857  real(RP) :: tflux(0:KA)
858  real(RP) :: qflux(0:KA)
859  real(RP) :: uflux(0:KA)
860  real(RP) :: vflux(0:KA)
861  real(RP) :: eflux(0:KA)
862 
863  real(RP) :: tmp
864 
865  logical :: do_put
866 
867  integer :: KE_PBL
868  integer :: k, i, j
869 
870 
871  ! work
872  real(RP) :: PTLV (KA)
873  real(RP) :: Nu_f (KA)
874  real(RP) :: Kh_f (KA)
875  real(RP) :: q2_2 (KA)
876  real(RP) :: ac (KA)
877  real(RP) :: betat(KA)
878  real(RP) :: betaq(KA)
879 
880  real(RP) :: f_smp (KA)
881  real(RP) :: f_shpgh(KA)
882  real(RP) :: f_gamma(KA)
883  real(RP) :: tltv25 (KA)
884  real(RP) :: qwtv25 (KA)
885  real(RP) :: tvsq25 (KA)
886  real(RP) :: tvsq_up(KA)
887  real(RP) :: tvsq_lo(KA)
888  real(RP) :: dtsq (KA)
889  real(RP) :: dqsq (KA)
890  real(RP) :: dcov (KA)
891 
892  real(RP) :: mflux(0:KA)
893 
894  real(RP) :: Uh(KA), Vh(KA), Wh(KA)
895  real(RP) :: qw2(KA), qh(KA)
896  real(RP) :: tlh(KA), tvh(KA)
897  real(RP) :: eh(KA), dh(KA)
898  real(RP) :: TEML(KA), LHVL(KA), psat(KA)
899 
900 #ifdef _OPENACC
901  real(RP) :: work(KA,4)
902 #endif
903 
904  !---------------------------------------------------------------------------
905 
906  dt = real(dt_dp, rp)
907 
908  log_progress(*) "atmosphere / physics / pbl / MYNN"
909 
910  mynn_level3 = ( atmos_phy_bl_mynn_level == "3" )
911 
912  select case ( bulkflux_type )
913  case ( 'B71' )
914  i_b_type = i_b71
915  case ( 'B91' )
916  i_b_type = i_b91
917  case ( 'B91W01' )
918  i_b_type = i_b91w01
919  case default
920  log_error("ATMOS_PHY_BL_MYNN_tendency",*) "BULKFLUX_type is invalid: ", trim(bulkflux_type)
921  call prc_abort
922  end select
923 
924  !$acc data copyin(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,F2H) &
925  !$acc copyout(RHOU_t,RHOV_t,RHOT_t,RHOQV_t,RPROG_t,Nu,Kh,Qlp,cldfrac,Zi,SFLX_BUOY), &
926  !$acc create(Ri,Pr,prod,diss,dudz2,l,flxU,flxV,flxT,flxQ)
927 
928 
929 
930 !OCL INDEPENDENT
931  !$omp parallel do default(none) &
932  !$omp OMP_SCHEDULE_ collapse(2) &
933  !$omp shared(KA,KS,KE,IS,IE,JS,JE, &
934  !$omp CPdry,CP_VAPOR,UNDEF, &
935  !$omp ATMOS_PHY_BL_MYNN_N2_MAX,ATMOS_PHY_BL_MYNN_TKE_MIN, &
936  !$omp ATMOS_PHY_BL_MYNN_NU_MAX,ATMOS_PHY_BL_MYNN_KH_MAX, &
937  !$omp ATMOS_PHY_BL_MYNN_Sq_fact, ATMOS_PHY_BL_MYNN_zeta_lim, &
938  !$omp ATMOS_PHY_BL_MYNN_similarity,ATMOS_PHY_BL_MYNN_dz_sim, &
939  !$omp ATMOS_PHY_BL_MYNN_DUMP_coef, &
940  !$omp ATMOS_PHY_BL_MYNN_PBL_MAX, &
941  !$omp ATMOS_PHY_BL_MYNN_K2010,ATMOS_PHY_BL_MYNN_O2019,ATMOS_PHY_BL_MYNN_MF, &
942  !$omp RHOU_t,RHOV_t,RHOT_t,RHOQV_t,RPROG_t,Nu,Kh,Qlp,cldfrac,Zi,SFLX_BUOY, &
943  !$omp DENS,PROG,U,V,W,POTT,PRES,QDRY,QV,Qw,POTV,POTL,EXNER,N2, &
944  !$omp SFC_DENS,SFLX_MU,SFLX_MV,SFLX_SH,SFLX_QV,us,ts,qs,RLmo, &
945  !$omp mynn_level3, &
946  !$omp frac_land, &
947  !$omp CZ,FZ,F2H,dt, &
948  !$omp I_B_TYPE, &
949  !$omp Ri,Pr,prod,diss,dudz2,l,flxU,flxV,flxT,flxQ,flxU2,flxV2,flxT2,flxQ2) &
950  !$omp private(N2_new,lq,sm25,sh25,rlqsm_h,q, &
951  !$omp SFLX_PT,CPtot,RHO,RHONu,RHOKh, &
952  !$omp smp,shpgh, &
953  !$omp dtldz,dqwdz,gammat,gammaq, &
954  !$omp flx,a,b,c,d,ap,rho_h,phi_n,tke,sf_t,CDZ,FDZ,z, &
955  !$omp dummy, &
956  !$omp tsq,qsq,cov, &
957  !$omp prod_t,prod_q,prod_c,diss_p, &
958  !$omp fmin,kmin, &
959  !$omp tflux,qflux,uflux,vflux,eflux, &
960  !$omp tmp, &
961  !$omp PTLV, Nu_f, Kh_f, q2_2, ac, betat, betaq, &
962  !$omp f_smp, f_shpgh, f_gamma, tltv25, qwtv25, tvsq25, &
963  !$omp tvsq_up, tvsq_lo, dtsq, dqsq, dcov, &
964  !$omp mflux, &
965  !$omp Uh, Vh, Wh, qw2, qh, tlh, tvh, eh, dh, &
966  !$omp TEML, LHVL, psat, &
967  !$omp KE_PBL,k,i,j)
968  !$acc kernels
969  !$acc loop independent collapse(2) &
970  !$acc private(N2_new,lq,sm25,sh25,rlqsm_h,q, &
971  !$acc SFLX_PT,CPtot,RHO,RHONu,RHOKh, &
972  !$acc smp,shpgh, &
973  !$acc dtldz,dqwdz,gammat,gammaq, &
974  !$acc flx,a,b,c,d,ap,rho_h,phi_n,tke,sf_t,CDZ,FDZ,z, &
975  !$acc dummy, &
976  !$acc tsq,qsq,cov, &
977  !$acc prod_t,prod_q,prod_c,diss_p, &
978  !$acc fmin,kmin, &
979  !$acc tflux,qflux,uflux,vflux,eflux, &
980  !$acc tmp, &
981  !$acc work, &
982  !$acc PTLV, Nu_f, Kh_f, q2_2, ac, betat, betaq, &
983  !$acc f_smp, f_shpgh, f_gamma, tltv25, qwtv25, tvsq25, &
984  !$acc tvsq_up, tvsq_lo, dtsq, dqsq, dcov, &
985  !$acc mflux, &
986  !$acc Uh, Vh, Wh, qw2, qh, tlh, tvh, eh, dh, &
987  !$acc TEML, LHVL, psat, &
988  !$acc KE_PBL,k,i,j)
989 
990  do j = js, je
991  do i = is, ie
992 
993  do k = ks, ke-1
994  z(k) = cz(k,i,j) - fz(ks-1,i,j)
995  end do
996 
997  ke_pbl = ks+1
998  !$acc loop seq
999  do k = ks+2, ke-1
1000  if ( atmos_phy_bl_mynn_pbl_max >= z(k) ) then
1001  ke_pbl = k
1002  else
1003  exit
1004  end if
1005  end do
1006 
1007  do k = ks, ke_pbl+1
1008  cdz(k) = fz(k,i,j) - fz(k-1,i,j)
1009  end do
1010  do k = ks, ke_pbl
1011  fdz(k) = cz(k+1,i,j) - cz(k,i,j)
1012  end do
1013 
1014  do k = ks, ke_pbl
1015  tke(k) = max( prog(k,i,j,i_tke), atmos_phy_bl_mynn_tke_min )
1016  end do
1017 
1018  if ( mynn_level3 ) then
1019  do k = ks, ke_pbl
1020  tsq(k) = max( prog(k,i,j,i_tsq), 0.0_rp )
1021  qsq(k) = max( prog(k,i,j,i_qsq), 0.0_rp )
1022  cov(k) = prog(k,i,j,i_cov)
1023  cov(k) = sign( min( abs(cov(k)), sqrt(tsq(k) * qsq(k))), cov(k) )
1024  end do
1025  end if
1026 
1027  cptot = cpdry + sflx_qv(i,j) * ( cp_vapor - cpdry )
1028  sflx_pt = sflx_sh(i,j) / ( cptot * exner(ks,i,j) )
1029 
1030  call mynn_main( ka, ks, ke_pbl, &
1031  i, j, &
1032  tke(:), tsq(:), qsq(:), cov(:), & ! (inout)
1033  q(:), l(:,i,j), lq(:), & ! (out)
1034  nu(:,i,j), rhonu(:), kh(:,i,j), rhokh(:), pr(:,i,j), & ! (out)
1035  prod(:,i,j), prod_t(:), prod_q(:), prod_c(:), & ! (out)
1036  diss(:,i,j), diss_p(:), & ! (out)
1037  sm25(:), smp(:), sh25(:), shpgh(:), & ! (out)
1038  gammat(:), gammaq(:), & ! (out)
1039  dudz2(:,i,j), n2_new(:), ri(:,i,j), & ! (out)
1040  dtldz(:), dqwdz(:), & ! (out)
1041  rho(:), rho_h(:), & ! (out)
1042  uflux(:), vflux(:), tflux(:), qflux(:), eflux(:), & ! (out)
1043  qlp(:,i,j), cldfrac(:,i,j), zi(i,j), sflx_buoy(i,j), & ! (out)
1044  u(:,i,j), v(:,i,j), w(:,i,j), & ! (in)
1045  dens(:,i,j), pres(:,i,j), & ! (in)
1046  pott(:,i,j), potl(:,i,j), potv(:,i,j), & ! (in)
1047  qw(:,i,j), n2(:,i,j), & ! (in)
1048  exner(:,i,j), qdry(:,i,j), & ! (in)
1049  sflx_pt, sflx_sh(i,j), sflx_qv(i,j), sfc_dens(i,j), & ! (in)
1050  rlmo(i,j), us(i,j), ts(i,j), qs(i,j), & ! (in)
1051  z(:), cdz(:), fdz(:), f2h(:,:,i,j), & ! (in)
1052  frac_land(i,j), & ! (in)
1053  dt, & ! (in)
1054  i_b_type, & ! (in)
1055  mynn_level3, .false., & ! (in)
1056 #ifdef _OPENACC
1057  work(:,:), & ! (work)
1058 #endif
1059  ptlv(:), nu_f(:), kh_f(:), q2_2(:), ac(:), & ! (work)
1060  betat(:), betaq(:), & ! (work)
1061  flx(:), a(:), b(:), c(:), d(:), & ! (work)
1062  f_smp(:), f_shpgh(:), f_gamma(:), & ! (work)
1063  tltv25(:), qwtv25(:), tvsq25(:), & ! (work)
1064  tvsq_up(:), tvsq_lo(:), dtsq(:), dqsq(:), dcov(:), & ! (work)
1065  mflux(:), & ! (work)
1066  uh(:), vh(:), wh(:), qw2(:), qh(:), & ! (work)
1067  tlh(:), tvh(:), eh(:), dh(:), & ! (work)
1068  teml(:), lhvl(:), psat(:) ) ! (work)
1069 
1070 
1071  ! time integration
1072 
1073  flx(ks-1 ) = 0.0_rp
1074  flx(ke_pbl) = 0.0_rp
1075 
1076  do k = ks, ke_pbl-1
1077  rlqsm_h(k) = rho_h(k) * ( f2h(k,1,i,j) * lq(k+1) * smp(k+1) + f2h(k,2,i,j) * lq(k) * smp(k) )
1078  end do
1079 
1080  ! dens * u
1081 
1082  if ( atmos_phy_bl_mynn_mf ) then
1083  do k = ks, ke_pbl-1
1084  flx(k) = uflux(k)
1085  end do
1086  else
1087  ! countergradient flux
1088  do k = ks, ke_pbl-1
1089  flx(k) = - rlqsm_h(k) * ( u(k+1,i,j) - u(k,i,j) ) / fdz(k)
1090  end do
1091  end if
1092 
1093  sf_t = sflx_mu(i,j) / cdz(ks)
1094  d(ks) = ( u(ks,i,j) * dens(ks,i,j) + dt * ( sf_t - flx(ks) / cdz(ks) ) ) / rho(ks)
1095  do k = ks+1, ke_pbl
1096  d(k) = u(k,i,j) - dt * ( flx(k) - flx(k-1) ) / ( cdz(k) * rho(k) )
1097  end do
1098  c(ks) = 0.0_rp
1099  do k = ks, ke_pbl-1
1100  ap = - dt * rhonu(k) / fdz(k)
1101  a(k) = ap / ( rho(k) * cdz(k) )
1102 #ifdef _OPENACC
1103  if ( k==ks ) then
1104  b(k) = - a(k) + 1.0_rp
1105  else
1106  b(k) = - a(k) + dt * rhonu(k-1) / ( fdz(k-1) * rho(k) * cdz(k) ) + 1.0_rp
1107  end if
1108 #else
1109  b(k) = - a(k) - c(k) + 1.0_rp
1110 #endif
1111  c(k+1) = ap / ( rho(k+1) * cdz(k+1) )
1112  end do
1113  a(ke_pbl) = 0.0_rp
1114  b(ke_pbl) = - c(ke_pbl) + 1.0_rp
1115 
1116  call matrix_solver_tridiagonal( &
1117  ka, ks, ke_pbl, &
1118 #ifdef _OPENACC
1119  work(:,:), & ! (work)
1120 #endif
1121  a(:), b(:), c(:), d(:), & ! (in)
1122  dummy(:) ) ! (out)
1123 ! phi_n(:) ) ! (out)
1124 
1125  phi_n(ks:ke_pbl) = dummy(ks:ke_pbl)
1126  rhou_t(ks,i,j) = ( phi_n(ks) * rho(ks) - u(ks,i,j) * dens(ks,i,j) ) / dt - sf_t
1127  do k = ks+1, ke_pbl
1128  rhou_t(k,i,j) = ( phi_n(k) - u(k,i,j) ) * rho(k) / dt
1129  end do
1130  do k = ke_pbl+1, ke
1131  rhou_t(k,i,j) = 0.0_rp
1132  end do
1133  flxu(ks-1,i,j) = 0.0_rp
1134  flxu2(ks-1,i,j) = 0.0_rp
1135  do k = ks, ke_pbl-1
1136  flxu(k,i,j) = flx(k) &
1137  - rhonu(k) * ( phi_n(k+1) - phi_n(k) ) / fdz(k)
1138  flxu2(k,i,j) = flx(k)
1139  end do
1140  do k = ke_pbl, ke
1141  flxu(k,i,j) = 0.0_rp
1142  flxu2(k,i,j) = 0.0_rp
1143  end do
1144 
1145 
1146  ! dens * v
1147 
1148  if ( atmos_phy_bl_mynn_mf ) then
1149  do k = ks, ke_pbl-1
1150  flx(k) = vflux(k)
1151  end do
1152  else
1153  ! countergradient flux
1154  do k = ks, ke_pbl-1
1155  flx(k) = - rlqsm_h(k) * ( v(k+1,i,j) - v(k,i,j) ) / fdz(k)
1156  end do
1157  end if
1158 
1159  sf_t = sflx_mv(i,j) / cdz(ks)
1160  d(ks) = ( v(ks,i,j) * dens(ks,i,j) + dt * ( sf_t - flx(ks) / cdz(ks) ) ) / rho(ks)
1161  do k = ks+1, ke_pbl
1162  d(k) = v(k,i,j) - dt * ( flx(k) - flx(k-1) ) / ( cdz(k) * rho(k) )
1163  end do
1164  ! a,b,c is the same as those for the u
1165 
1166  call matrix_solver_tridiagonal( &
1167  ka, ks, ke_pbl, &
1168 #ifdef _OPENACC
1169  work(:,:), & ! (work)
1170 #endif
1171  a(:), b(:), c(:), d(:), & ! (in)
1172  phi_n(:) ) ! (out)
1173 
1174  rhov_t(ks,i,j) = ( phi_n(ks) * rho(ks) - v(ks,i,j) * dens(ks,i,j) ) / dt - sf_t
1175  do k = ks+1, ke_pbl
1176  rhov_t(k,i,j) = ( phi_n(k) - v(k,i,j) ) * rho(k) / dt
1177  end do
1178  do k = ke_pbl+1, ke
1179  rhov_t(k,i,j) = 0.0_rp
1180  end do
1181  flxv(ks-1,i,j) = 0.0_rp
1182  flxv2(ks-1,i,j) = 0.0_rp
1183  do k = ks, ke_pbl-1
1184  flxv(k,i,j) = flx(k) &
1185  - rhonu(k) * ( phi_n(k+1) - phi_n(k) ) / fdz(k)
1186  flxv2(k,i,j) = flx(k)
1187  end do
1188  do k = ke_pbl, ke
1189  flxv(k,i,j) = 0.0_rp
1190  flxv2(k,i,j) = 0.0_rp
1191  end do
1192 
1193 
1194  ! dens * pott
1195 
1196  if ( atmos_phy_bl_mynn_mf ) then
1197  do k = ks, ke_pbl-1
1198  flx(k) = tflux(k)
1199  end do
1200  else
1201  ! countergradient flux
1202  do k = ks, ke_pbl-1
1203  flx(k) = - ( f2h(k,1,i,j) * lq(k+1) * gammat(k+1) + f2h(k,2,i,j) * lq(k) * gammat(k) ) &
1204  * rho_h(k)
1205  end do
1206  end if
1207 
1208  sf_t = sflx_pt / cdz(ks)
1209  ! assume that induced vapor has the same PT
1210  d(ks) = pott(ks,i,j) + dt * ( sf_t - flx(ks) / cdz(ks) ) / rho(ks)
1211  do k = ks+1, ke_pbl
1212  d(k) = pott(k,i,j) - dt * ( flx(k) - flx(k-1) ) / ( cdz(k) * rho(k) )
1213  end do
1214 
1215  c(ks) = 0.0_rp
1216  do k = ks, ke_pbl-1
1217  ap = - dt * rhokh(k) / fdz(k)
1218  a(k) = ap / ( rho(k) * cdz(k) )
1219 #ifdef _OPENACC
1220  if ( k==ks ) then
1221  b(k) = - a(k) + 1.0_rp
1222  else
1223  b(k) = - a(k) + dt * rhokh(k-1) / ( fdz(k-1) * rho(k) * cdz(k) ) + 1.0_rp
1224  end if
1225 #else
1226  b(k) = - a(k) - c(k) + 1.0_rp
1227 #endif
1228  c(k+1) = ap / ( rho(k+1) * cdz(k+1) )
1229  end do
1230  a(ke_pbl) = 0.0_rp
1231  b(ke_pbl) = - c(ke_pbl) + 1.0_rp
1232 
1233  call matrix_solver_tridiagonal( &
1234  ka, ks, ke_pbl, &
1235 #ifdef _OPENACC
1236  work(:,:), & ! (work)
1237 #endif
1238  a(:), b(:), c(:), d(:), & ! (in)
1239  phi_n(:) ) ! (out)
1240 
1241  rhot_t(ks,i,j) = ( phi_n(ks) - pott(ks,i,j) ) * rho(ks) / dt - sf_t
1242  do k = ks+1, ke_pbl
1243  rhot_t(k,i,j) = ( phi_n(k) - pott(k,i,j) ) * rho(k) / dt
1244  end do
1245  do k = ke_pbl+1, ke
1246  rhot_t(k,i,j) = 0.0_rp
1247  end do
1248  flxt(ks-1,i,j) = 0.0_rp
1249  flxt2(ks-1,i,j) = 0.0_rp
1250  do k = ks, ke_pbl-1
1251  flxt(k,i,j) = flx(k) &
1252  - rhokh(k) * ( phi_n(k+1) - phi_n(k) ) / fdz(k)
1253  flxt2(k,i,j) = flx(k)
1254  end do
1255  do k = ke_pbl, ke
1256  flxt(k,i,j) = 0.0_rp
1257  flxt2(k,i,j) = 0.0_rp
1258  end do
1259 
1260  kmin = ks-1
1261  if ( flxt(ks,i,j) > 1e-4_rp ) then
1262  fmin = flxt(ks,i,j) / rho_h(ks)
1263  !$acc loop seq
1264  do k = ks+1, ke_pbl-2
1265  tmp = ( flxt(k-1,i,j) + flxt(k,i,j) + flxt(k+1,i,j) ) &
1266  / ( rho_h(k-1) + rho_h(k) + rho_h(k+1) ) ! running mean
1267  if ( fmin < 0.0_rp .and. tmp > fmin ) exit
1268  if ( tmp < fmin ) then
1269  fmin = tmp
1270  kmin = k
1271  end if
1272  end do
1273  end if
1274  zi(i,j) = fz(kmin,i,j) - fz(ks-1,i,j)
1275 
1276  ! dens * qv
1277 
1278  if ( atmos_phy_bl_mynn_mf ) then
1279  do k = ks, ke_pbl-1
1280  flx(k) = qflux(k)
1281  end do
1282  else
1283  ! countergradient flux
1284  do k = ks, ke_pbl-1
1285  flx(k) = - ( f2h(k,1,i,j) * lq(k+1) * gammaq(k+1) + f2h(k,2,i,j) * lq(k) * gammaq(k) ) &
1286  * rho_h(k)
1287  end do
1288  end if
1289 
1290  sf_t = sflx_qv(i,j) / cdz(ks)
1291  d(ks) = ( qv(ks,i,j) * dens(ks,i,j) + dt * ( sf_t - flx(ks) / cdz(ks) ) ) / rho(ks)
1292  do k = ks+1, ke_pbl
1293  d(k) = qv(k,i,j) - dt * ( flx(k) - flx(k-1) ) / ( cdz(k) * rho(k) )
1294  end do
1295 
1296  call matrix_solver_tridiagonal( &
1297  ka, ks, ke_pbl, &
1298 #ifdef _OPENACC
1299  work(:,:), & ! (work)
1300 #endif
1301  a(:), b(:), c(:), d(:), & ! (in)
1302  phi_n(:) ) ! (out)
1303 
1304  rhoqv_t(ks,i,j) = ( phi_n(ks) * rho(ks) - qv(ks,i,j) * dens(ks,i,j) ) / dt - sf_t
1305  do k = ks+1, ke_pbl
1306  rhoqv_t(k,i,j) = ( phi_n(k) - qv(k,i,j) ) * rho(k) / dt
1307  end do
1308  do k = ke_pbl+1, ke
1309  rhoqv_t(k,i,j) = 0.0_rp
1310  end do
1311  flxq(ks-1,i,j) = 0.0_rp
1312  flxq2(ks-1,i,j) = 0.0_rp
1313  do k = ks, ke_pbl-1
1314  flxq(k,i,j) = flx(k) &
1315  - rhokh(k) * ( phi_n(k+1) - phi_n(k) ) / fdz(k)
1316  flxq2(k,i,j) = flx(k)
1317  end do
1318  do k = ke_pbl, ke
1319  flxq(k,i,j) = 0.0_rp
1320  flxq2(k,i,j) = 0.0_rp
1321  end do
1322 
1323  ! dens * TKE
1324 
1325 !!$ if ( ATMOS_PHY_BL_MYNN_MF ) then
1326 !!$ do k = KS, KE_PBL-1
1327 !!$ flx(k) = eflux(k)
1328 !!$ end do
1329 !!$ else
1330  do k = ks, ke_pbl-1
1331  flx(k) = 0.0_rp
1332  end do
1333 !!$ end if
1334 
1335  do k = ks, ke_pbl-1
1336  d(k) = tke(k) * dens(k,i,j) / rho(k) + dt * ( - ( flx(k) - flx(k-1) ) / ( cdz(k) * rho(k) ) + prod(k,i,j) )
1337  end do
1338  d(ke_pbl) = 0.0_rp
1339 
1340  c(ks) = 0.0_rp
1341  do k = ks, ke_pbl-1
1342  ap = - dt * atmos_phy_bl_mynn_sq_fact * rhonu(k) / fdz(k)
1343  a(k) = ap / ( rho(k) * cdz(k) )
1344 #ifdef _OPENACC
1345  if ( k==ks ) then
1346  b(k) = - a(k) + 1.0_rp - diss(k,i,j) * dt
1347  else
1348  b(k) = - a(k) + dt * atmos_phy_bl_mynn_sq_fact * rhonu(k-1) / ( fdz(k-1) * rho(k) * cdz(k) ) + 1.0_rp - diss(k,i,j) * dt
1349  end if
1350 #else
1351  b(k) = - a(k) - c(k) + 1.0_rp - diss(k,i,j) * dt
1352 #endif
1353  c(k+1) = ap / ( rho(k+1) * cdz(k+1) )
1354  end do
1355  a(ke_pbl) = 0.0_rp
1356 
1357  call matrix_solver_tridiagonal( &
1358  ka, ks, ke_pbl, &
1359 #ifdef _OPENACC
1360  work(:,:), & ! (work)
1361 #endif
1362  a(:), b(:), c(:), d(:), & ! (in)
1363  phi_n(:) ) ! (out)
1364 
1365  do k = ks, ke_pbl
1366  phi_n(k) = max( phi_n(k), atmos_phy_bl_mynn_tke_min )
1367  end do
1368 
1369  do k = ks, ke_pbl
1370  diss(k,i,j) = diss(k,i,j) * phi_n(k)
1371  rprog_t(k,i,j,i_tke) = ( phi_n(k) * rho(k) - prog(k,i,j,i_tke) * dens(k,i,j) ) / dt
1372  end do
1373  !$acc loop independent
1374  do k = ke_pbl+1, ke
1375  rprog_t(k,i,j,i_tke) = - atmos_phy_bl_mynn_dump_coef * prog(k,i,j,i_tke) * dens(k,i,j) / dt
1376  diss(k,i,j) = 0.0_rp
1377  prod(k,i,j) = 0.0_rp
1378  end do
1379 
1380 
1381  if ( mynn_level3 ) then
1382 
1383  ! dens * tsq
1384 
1385  do k = ks, ke_pbl
1386  diss_p(k) = dt * 2.0_rp * q(k) / ( b2 * l(k,i,j) )
1387  end do
1388  do k = ks, ke_pbl-1
1389  d(k) = tsq(k) * dens(k,i,j) / rho(k) + dt * prod_t(k)
1390  end do
1391  d(ke_pbl) = 0.0_rp
1392  c(ks) = 0.0_rp
1393  do k = ks, ke_pbl-1
1394  ap = - dt * rhonu(k) / fdz(k)
1395  a(k) = ap / ( rho(k) * cdz(k) )
1396 #ifdef _OPENACC
1397  if ( k==ks ) then
1398  b(k) = - a(k) + diss_p(k)
1399  else
1400  b(k) = - a(k) + dt * rhonu(k-1) / ( fdz(k-1) * rho(k) * cdz(k) ) + diss_p(k)
1401  end if
1402 #else
1403  b(k) = - a(k) - c(k) + 1.0_rp + diss_p(k)
1404 #endif
1405  c(k+1) = ap / ( rho(k+1) * cdz(k+1) )
1406  end do
1407  a(ke_pbl) = 0.0_rp
1408  b(ke_pbl) = - c(ke_pbl) + 1.0_rp + diss_p(ke_pbl)
1409 
1410  call matrix_solver_tridiagonal( &
1411  ka, ks, ke_pbl, &
1412 #ifdef _OPENACC
1413  work(:,:), & ! (work)
1414 #endif
1415  a(:), b(:), c(:), d(:), & ! (in)
1416  tsq(:) ) ! (out)
1417 
1418  do k = ks, ke_pbl
1419  tsq(k) = max( tsq(k), 0.0_rp )
1420  end do
1421 
1422  do k = ks, ke_pbl
1423  rprog_t(k,i,j,i_tsq) = ( tsq(k) * rho(k) - prog(k,i,j,i_tsq) * dens(k,i,j) ) / dt
1424  end do
1425  !$acc loop independent
1426  do k = ke_pbl+1, ke
1427  rprog_t(k,i,j,i_tsq) = - atmos_phy_bl_mynn_dump_coef * prog(k,i,j,i_tsq) * dens(k,i,j) / dt
1428  end do
1429 
1430 
1431  ! dens * qsq
1432 
1433  do k = ks, ke_pbl
1434  d(k) = qsq(k) * dens(k,i,j) / rho(k) + dt * prod_q(k)
1435  end do
1436  ! a, b, c are same as those for tsq
1437 
1438  call matrix_solver_tridiagonal( &
1439  ka, ks, ke_pbl, &
1440 #ifdef _OPENACC
1441  work(:,:), & ! (work)
1442 #endif
1443  a(:), b(:), c(:), d(:), & ! (in)
1444  qsq(:) ) ! (out)
1445 
1446  do k = ks, ke_pbl
1447  qsq(k) = max( qsq(k), 0.0_rp )
1448  end do
1449 
1450  do k = ks, ke_pbl
1451  rprog_t(k,i,j,i_qsq) = ( qsq(k) * rho(k) - prog(k,i,j,i_qsq) * dens(k,i,j) ) / dt
1452  end do
1453  !$acc loop independent
1454  do k = ke_pbl+1, ke
1455  rprog_t(k,i,j,i_qsq) = - atmos_phy_bl_mynn_dump_coef * prog(k,i,j,i_qsq) * dens(k,i,j) / dt
1456  end do
1457 
1458 
1459  ! dens * cov
1460 
1461  do k = ks, ke_pbl-1
1462  d(k) = cov(k) * dens(k,i,j) / rho(k) + dt * prod_c(k)
1463  end do
1464  d(ke_pbl) = 0.0_rp
1465  ! a, b, c are same as those for tsq
1466 
1467  call matrix_solver_tridiagonal( &
1468  ka, ks, ke_pbl, &
1469 #ifdef _OPENACC
1470  work(:,:), & ! (work)
1471 #endif
1472  a(:), b(:), c(:), d(:), & ! (in)
1473  cov(:) ) ! (out)
1474 
1475  do k = ks, ke_pbl
1476  cov(k) = sign( min( abs(cov(k)), sqrt(tsq(k)*qsq(k)) ), cov(k) )
1477  end do
1478 
1479  do k = ks, ke_pbl
1480  rprog_t(k,i,j,i_cov) = ( cov(k) * rho(k) - prog(k,i,j,i_cov) * dens(k,i,j) ) / dt
1481  end do
1482  !$acc loop independent
1483  do k = ke_pbl+1, ke
1484  rprog_t(k,i,j,i_cov) = - atmos_phy_bl_mynn_dump_coef * prog(k,i,j,i_cov) * dens(k,i,j) / dt
1485  end do
1486 
1487  end if
1488 
1489  nu(ks-1,i,j) = 0.0_rp
1490  kh(ks-1,i,j) = 0.0_rp
1491  do k = ke_pbl, ke
1492  nu(k,i,j) = 0.0_rp
1493  kh(k,i,j) = 0.0_rp
1494  pr(k,i,j) = 1.0_rp
1495  prod(k,i,j) = 0.0_rp
1496  diss(k,i,j) = 0.0_rp
1497  end do
1498  !$acc loop independent
1499  do k = ke_pbl+1, ke
1500  ri(k,i,j) = undef
1501  dudz2(k,i,j) = undef
1502  l(k,i,j) = undef
1503  qlp(k,i,j) = undef
1504  cldfrac(k,i,j) = undef
1505  end do
1506 
1507 
1508  end do
1509  end do
1510  !$acc end kernels
1511 
1512 
1513  call file_history_query(hist_ri, do_put)
1514  if ( do_put ) call file_history_put(hist_ri, ri(:,:,:))
1515  call file_history_query(hist_pr, do_put)
1516  if ( do_put ) call file_history_put(hist_pr, pr(:,:,:))
1517  call file_history_query(hist_tke_pr, do_put)
1518  if ( do_put ) call file_history_put(hist_tke_pr, prod(:,:,:))
1519  call file_history_query(hist_tke_di, do_put)
1520  if ( do_put ) call file_history_put(hist_tke_di, diss(:,:,:))
1521  call file_history_query(hist_dudz2, do_put)
1522  if ( do_put ) call file_history_put(hist_dudz2, dudz2(:,:,:))
1523  call file_history_query(hist_lmix, do_put)
1524  if ( do_put ) call file_history_put(hist_lmix, l(:,:,:))
1525 
1526  call file_history_query(hist_flxu, do_put)
1527  if ( do_put ) call file_history_put(hist_flxu, flxu(1:,:,:))
1528  call file_history_query(hist_flxv, do_put)
1529  if ( do_put ) call file_history_put(hist_flxv, flxv(1:,:,:))
1530  call file_history_query(hist_flxt, do_put)
1531  if ( do_put ) call file_history_put(hist_flxt, flxt(1:,:,:))
1532  call file_history_query(hist_flxq, do_put)
1533  if ( do_put ) call file_history_put(hist_flxq, flxq(1:,:,:))
1534 
1535  call file_history_query(hist_flxu2, do_put)
1536  if ( do_put ) call file_history_put(hist_flxu2, flxu2(1:,:,:))
1537  call file_history_query(hist_flxv2, do_put)
1538  if ( do_put ) call file_history_put(hist_flxv2, flxv2(1:,:,:))
1539  call file_history_query(hist_flxt2, do_put)
1540  if ( do_put ) call file_history_put(hist_flxt2, flxt2(1:,:,:))
1541  call file_history_query(hist_flxq2, do_put)
1542  if ( do_put ) call file_history_put(hist_flxq2, flxq2(1:,:,:))
1543 
1544  !$acc end data
1545 
1546  return
1547  end subroutine atmos_phy_bl_mynn_tendency
1548 
1549  !-----------------------------------------------------------------------------
1550  ! private routines
1551  !-----------------------------------------------------------------------------
1552 !OCL SERIAL
1553  subroutine mynn_main( &
1554  KA, KS, KE_PBL, &
1555  i, j, &
1556  tke, tsq, qsq, cov, &
1557  q, l, lq, &
1558  Nu, RHONu, Kh, RHOKh, Pr, &
1559  prod, prod_t, prod_q, prod_c, &
1560  diss, diss_p, &
1561  sm25, smp, sh25, shpgh, &
1562  gammat, gammaq, &
1563  dudz2, n2_new, Ri, &
1564  dtldz, dqwdz, &
1565  RHO, RHO_h, &
1566  uflux, vflux, tflux, qflux, eflux, &
1567  Qlp, cldfrac, Zi, SFLX_BUOY, &
1568  U, V, W, &
1569  DENS, PRES, &
1570  POTT, POTL, POTV, &
1571  Qw, N2, &
1572  EXNER, QDRY, &
1573  SFLX_PT, SFLX_SH, SFLX_QV, SFC_DENS, &
1574  RLmo, us, ts, qs, &
1575  z, CDZ, FDZ, F2H, &
1576  frac_land, &
1577  dt, &
1578  I_B_TYPE, &
1579  mynn_level3, initialize, &
1580 #ifdef _OPENACC
1581  work, &
1582 #endif
1583  PTLV, Nu_f, Kh_f, q2_2, ac, betat, betaq, &
1584  flx, a, b, c, d, &
1585  f_smp, f_shpgh, f_gamma, tltv25, qwtv25, tvsq25, &
1586  tvsq_up, tvsq_lo, dtsq, dqsq, dcov, &
1587  mflux, &
1588  Uh, Vh, Wh, qw2, qh, tlh, tvh, eh, dh, &
1589  TEML, LHVL, psat )
1590  !$acc routine vector
1591  use scale_const, only: &
1592  eps => const_eps, &
1593  karman => const_karman, &
1594  grav => const_grav, &
1595  epstvap => const_epstvap
1596  use scale_matrix, only: &
1597  matrix_solver_tridiagonal
1598  integer, intent(in) :: KA, KS, KE_PBL
1599  integer, intent(in) :: i, j
1600  real(RP), intent(inout) :: tke(KA)
1601  real(RP), intent(inout) :: tsq(KA)
1602  real(RP), intent(inout) :: qsq(KA)
1603  real(RP), intent(inout) :: cov(KA)
1604  real(RP), intent(out) :: q(KA)
1605  real(RP), intent(out) :: l(KA)
1606  real(RP), intent(out) :: lq(KA)
1607  real(RP), intent(out) :: Nu(KA)
1608  real(RP), intent(out) :: RHONu(KA)
1609  real(RP), intent(out) :: Kh(KA)
1610  real(RP), intent(out) :: RHOKh(KA)
1611  real(RP), intent(out) :: Pr(KA)
1612  real(RP), intent(out) :: prod(KA)
1613  real(RP), intent(out) :: prod_t(KA)
1614  real(RP), intent(out) :: prod_q(KA)
1615  real(RP), intent(out) :: prod_c(KA)
1616  real(RP), intent(out) :: diss(KA)
1617  real(RP), intent(out) :: diss_p(KA)
1618  real(RP), intent(out) :: sm25(KA)
1619  real(RP), intent(out) :: smp(KA)
1620  real(RP), intent(out) :: sh25(KA)
1621  real(RP), intent(out) :: shpgh(KA)
1622  real(RP), intent(out) :: gammat(KA)
1623  real(RP), intent(out) :: gammaq(KA)
1624  real(RP), intent(out) :: dudz2(KA)
1625  real(RP), intent(out) :: n2_new(KA)
1626  real(RP), intent(out) :: Ri(KA)
1627  real(RP), intent(out) :: dtldz(KA)
1628  real(RP), intent(out) :: dqwdz(KA)
1629  real(RP), intent(out) :: RHO(KA)
1630  real(RP), intent(out) :: RHO_h(KA)
1631  real(RP), intent(out) :: uflux(KA)
1632  real(RP), intent(out) :: vflux(KA)
1633  real(RP), intent(out) :: tflux(KA)
1634  real(RP), intent(out) :: qflux(KA)
1635  real(RP), intent(out) :: eflux(KA)
1636  real(RP), intent(out) :: Qlp(KA)
1637  real(RP), intent(out) :: cldfrac(KA)
1638  real(RP), intent(out) :: Zi
1639  real(RP), intent(out) :: SFLX_BUOY
1640  real(RP), intent(in) :: U(KA)
1641  real(RP), intent(in) :: V(KA)
1642  real(RP), intent(in) :: W(KA)
1643  real(RP), intent(in) :: DENS(KA)
1644  real(RP), intent(in) :: PRES(KA)
1645  real(RP), intent(in) :: POTT(KA)
1646  real(RP), intent(in) :: POTL(KA)
1647  real(RP), intent(in) :: POTV(KA)
1648  real(RP), intent(in) :: Qw(KA)
1649  real(RP), intent(in) :: N2(KA)
1650  real(RP), intent(in) :: EXNER(KA)
1651  real(RP), intent(in) :: QDRY(KA)
1652  real(RP), intent(in) :: SFLX_PT
1653  real(RP), intent(in) :: SFLX_SH
1654  real(RP), intent(in) :: SFLX_QV
1655  real(RP), intent(in) :: SFC_DENS
1656  real(RP), intent(in) :: RLmo
1657  real(RP), intent(in) :: us
1658  real(RP), intent(in) :: ts
1659  real(RP), intent(in) :: qs
1660  real(RP), intent(in) :: z(KA)
1661  real(RP), intent(in) :: CDZ(KA)
1662  real(RP), intent(in) :: FDZ(KA)
1663  real(RP), intent(in) :: F2H(KA,2)
1664  real(RP), intent(in) :: frac_land
1665  real(RP), intent(in) :: dt
1666  integer, intent(in) :: I_B_TYPE
1667  logical, intent(in) :: mynn_level3
1668  logical, intent(in) :: initialize
1669 
1670  ! work
1671  real(RP), intent(out) :: PTLV (KA)
1672  real(RP), intent(out) :: Nu_f (KA)
1673  real(RP), intent(out) :: Kh_f (KA)
1674  real(RP), intent(out) :: q2_2 (KA)
1675  real(RP), intent(out) :: ac (KA)
1676  real(RP), intent(out) :: betat(KA)
1677  real(RP), intent(out) :: betaq(KA)
1678 
1679  real(RP), intent(out) :: flx(0:KA)
1680  real(RP), intent(out) :: a(KA)
1681  real(RP), intent(out) :: b(KA)
1682  real(RP), intent(out) :: c(KA)
1683  real(RP), intent(out) :: d(KA)
1684 
1685  real(RP) :: ap
1686 
1687  real(RP) :: zeta
1688  real(RP) :: phi_m, phi_h
1689  real(RP) :: us3
1690 
1691  ! for level 3 (work)
1692  real(RP), intent(out) :: f_smp (KA)
1693  real(RP), intent(out) :: f_shpgh(KA)
1694  real(RP), intent(out) :: f_gamma(KA)
1695  real(RP), intent(out) :: tltv25 (KA)
1696  real(RP), intent(out) :: qwtv25 (KA)
1697  real(RP), intent(out) :: tvsq25 (KA)
1698  real(RP), intent(out) :: tvsq_up(KA)
1699  real(RP), intent(out) :: tvsq_lo(KA)
1700  real(RP), intent(out) :: dtsq (KA)
1701  real(RP), intent(out) :: dqsq (KA)
1702  real(RP), intent(out) :: dcov (KA)
1703  real(RP) :: tvsq
1704  real(RP) :: tltv
1705  real(RP) :: qwtv
1706  real(RP) :: wtl
1707  real(RP) :: wqw
1708 
1709  ! for O2019 (work)
1710  real(RP), intent(out) :: mflux(0:KA)
1711 
1712  ! work
1713  real(RP), intent(out) :: Uh(KA), Vh(KA), Wh(KA)
1714  real(RP), intent(out) :: qw2(KA), qh(KA)
1715  real(RP), intent(out) :: tlh(KA), tvh(KA)
1716  real(RP), intent(out) :: eh(KA), dh(KA)
1717  real(RP), intent(out) :: TEML(KA), LHVL(KA), psat(KA)
1718 
1719 #ifdef _OPENACC
1720  real(RP), intent(out) :: work(KA,4)
1721 #endif
1722 
1723  real(RP) :: sw, tmp
1724 
1725  integer :: k, it, nit
1726 
1727 
1728  if ( atmos_phy_bl_mynn_similarity .or. initialize ) then
1729  call get_phi( zeta, phi_m, phi_h, & ! (out)
1730  z(ks), rlmo, i_b_type ) ! (in)
1731  end if
1732 
1733  call calc_vertical_differece( ka, ks, ke_pbl, &
1734  i, j, &
1735  dudz2(:), dtldz(:), dqwdz(:), & ! (out)
1736  u(:), v(:), potl(:), & ! (in)
1737  qw(:), qdry(:), & ! (in)
1738  cdz(:), fdz(:), f2h(:,:), & ! (in)
1739  us, ts, qs, & ! (in)
1740  phi_m, phi_h, z(ks), & ! (in)
1741  uh(:), vh(:), qw2(:), qh(:) ) ! (work)
1742 
1743  us3 = us**3
1744 
1745  do k = ks, ke_pbl
1746  q(k) = sqrt( max( tke(k), atmos_phy_bl_mynn_tke_min ) * 2.0_rp )
1747  end do
1748 
1749  if ( atmos_phy_bl_mynn_use_zi ) then
1750  do k = ks, ke_pbl
1751  ptlv(k) = potl(k) * ( 1.0_rp + epstvap * qw(k) )
1752  end do
1753  call calc_zi_o2019( ka, ks, ke_pbl, &
1754  zi, & ! (out)
1755  ptlv(:), tke(:), & ! (in)
1756  z(:), frac_land, & ! (in)
1757  initialize )
1758  end if
1759 
1760  if ( initialize .or. (.not. mynn_level3) ) then
1761  ! estimate tsq, qsq, and cov
1762 
1763  do k = ks, ke_pbl
1764  n2_new(k) = min( max( n2(k), - atmos_phy_bl_mynn_n2_max ), atmos_phy_bl_mynn_n2_max )
1765  ri(k) = n2_new(k) / dudz2(k)
1766  end do
1767  if ( initialize ) then
1768  dudz2(ks) = max( dudz2(ks), 1e-4_rp )
1769  n2_new(ks) = min( n2_new(ks), 0.0_rp )
1770  ri(ks) = n2_new(ks) / dudz2(ks)
1771  end if
1772 
1773  sflx_buoy = - us3 * rlmo / karman
1774 
1775  if ( atmos_phy_bl_mynn_mf ) then
1776  do k = ks, ke_pbl
1777  cldfrac(k) = 0.0_rp
1778  end do
1779  call calc_mflux( &
1780  ka, ks, ke_pbl, &
1781  dens(:), potv(:), potl(:), qw(:), & ! (in)
1782  u(:), v(:), w(:), tke(:), & ! (in)
1783  cldfrac(:), & ! (in)
1784  exner(:), & ! (in)
1785  sflx_sh, sflx_buoy, & ! (in)
1786  sflx_pt, sflx_qv, & ! (in)
1787  sfc_dens, & ! (in)
1788  zi, z(:), cdz(:), f2h(:,:), & ! (in)
1789  mflux(:), & ! (out)
1790  tflux(:), qflux(:), & ! (out)
1791  uflux(:), vflux(:), eflux(:), & ! (out)
1792  tlh(:), tvh(:), qh(:), & ! (work)
1793  uh(:), vh(:), wh(:), eh(:), dh(:) ) ! (work)
1794  end if
1795 
1796  ! length
1797  call get_length( &
1798  ka, ks, ke_pbl, &
1799  q(:), n2_new(:), & ! (in)
1800  potv(:), dens(:), & ! (in)
1801  mflux(:), & ! (in)
1802  zi, & ! (in)
1803  sflx_buoy, rlmo, & ! (in)
1804  z(:), cdz(:), fdz(:), & ! (in)
1805  initialize, & ! (in)
1806  l(:) ) ! (out)
1807 
1808  call get_q2_level2( &
1809  ka, ks, ke_pbl, &
1810  dudz2(:), ri(:), l(:), & ! (in)
1811  q2_2(:) ) ! (out)
1812 
1813  if ( initialize ) then
1814  do k = ks, ke_pbl
1815  q(k) = sqrt( q2_2(k) )
1816  end do
1817  end if
1818 
1819  do k = ks, ke_pbl
1820  ac(k) = min( q(k) / sqrt( q2_2(k) + 1e-20_rp ), 1.0_rp )
1821  end do
1822 
1823  call get_smsh( &
1824  ka, ks, ke_pbl, &
1825  i, j, &
1826  q(:), ac(:), & ! (in)
1827  l(:), n2_new(:), & ! (in)
1828  ri(:), & ! (in)
1829  potv(:), dudz2(:), & ! (in)
1830  dtldz(:), dqwdz(:), & ! (in)
1831  betat(:), betaq(:), & ! (in) ! dummy
1832  .false., .false., & ! (in)
1833  tsq(:), qsq(:), cov(:), & ! (inout)
1834  sm25(:), f_smp(:), & ! (out) ! dummy
1835  sh25(:), f_shpgh(:), & ! (out) ! dymmy
1836  f_gamma(:), & ! (out) ! dummy
1837  tltv25(:), qwtv25(:), & ! (out) ! dummy
1838  tvsq25(:), & ! (out) ! dummy
1839  tvsq_up(:), tvsq_lo(:) ) ! (out) ! dummy
1840 
1841  end if
1842 
1843  flx(ks-1 ) = 0.0_rp
1844  flx(ke_pbl) = 0.0_rp
1845 
1846  if ( initialize ) then
1847  nit = ke_pbl - 1
1848  else
1849  nit = 1
1850  end if
1851 
1852  !$acc loop seq
1853  do it = 1, nit
1854 
1855  call partial_condensation( ka, ks, ke_pbl, &
1856  betat(:), betaq(:), & ! (out)
1857  qlp(:), cldfrac(:), & ! (out)
1858  pres(:), pott(:), & ! (in)
1859  potl(:), qw(:), & ! (in)
1860  exner(:), & ! (in)
1861  tsq(:), qsq(:), cov(:), & ! (in)
1862  teml(:), lhvl(:), psat(:) ) ! (work)
1863 
1864  ! update N2
1865  do k = ks, ke_pbl
1866  n2_new(k) = min(atmos_phy_bl_mynn_n2_max, &
1867  grav * ( dtldz(k) * betat(k) + dqwdz(k) * betaq(k) ) / potv(k) )
1868  ri(k) = n2_new(k) / dudz2(k)
1869  end do
1870 
1871  sflx_buoy = grav / potv(ks) * ( betat(ks) * sflx_pt + betaq(ks) * sflx_qv ) / sfc_dens
1872 
1873 
1874  if ( atmos_phy_bl_mynn_use_zi ) then
1875  do k = ks, ke_pbl
1876  ptlv(k) = potl(k) * betat(k) + qw(k) * betaq(k)
1877  end do
1878  call calc_zi_o2019( ka, ks, ke_pbl, &
1879  zi, & ! (out)
1880  ptlv(:), tke(:), & ! (in)
1881  z(:), frac_land, & ! (in)
1882  .false. ) ! (in)
1883  end if
1884 
1885  if ( atmos_phy_bl_mynn_mf ) then
1886  call calc_mflux( &
1887  ka, ks, ke_pbl, &
1888  dens(:), potv(:), potl(:), qw(:), & ! (in)
1889  u(:), v(:), w(:), tke(:), & ! (in)
1890  cldfrac(:), & ! (in)
1891  exner(:), & ! (in)
1892  sflx_sh, sflx_buoy, & ! (in)
1893  sflx_pt, sflx_qv, & ! (in)
1894  sfc_dens, & ! (in)
1895  zi, z(:), cdz(:), f2h(:,:), & ! (in)
1896  mflux(:), & ! (out)
1897  tflux(:), qflux(:), & ! (out)
1898  uflux(:), vflux(:), eflux(:), & ! (out)
1899  tlh(:), tvh(:), qh(:), & ! (work)
1900  uh(:), vh(:), wh(:), eh(:), dh(:) ) ! (work)
1901  end if
1902 
1903 
1904  ! length
1905  call get_length( &
1906  ka, ks, ke_pbl, &
1907  q(:), n2_new(:), & ! (in)
1908  potv(:), dens(:), & ! (in)
1909  mflux(:), & ! (in)
1910  zi, & ! (in)
1911  sflx_buoy, rlmo, & ! (in)
1912  z(:), cdz(:), fdz(:), & ! (in)
1913  .false., & ! (in)
1914  l(:) ) ! (out)
1915 
1916  call get_q2_level2( &
1917  ka, ks, ke_pbl, &
1918  dudz2(:), ri(:), l(:), & ! (in)
1919  q2_2(:) ) ! (out)
1920 
1921  do k = ks, ke_pbl
1922  ac(k) = min( q(k) / sqrt( q2_2(k) + 1e-20_rp ), 1.0_rp )
1923  end do
1924 
1925  call get_smsh( &
1926  ka, ks, ke_pbl, &
1927  i, j, & ! (in)
1928  q(:), ac(:), & ! (in)
1929  l(:), n2_new(:), & ! (in)
1930  ri(:), & ! (in)
1931  potv(:), dudz2(:), & ! (in)
1932  dtldz(:), dqwdz(:), & ! (in)
1933  betat(:), betaq(:), & ! (in)
1934  mynn_level3, & ! (in)
1935  initialize .and. it==1, & ! (in)
1936  tsq(:), qsq(:), cov(:), & ! (inout)
1937  sm25(:), f_smp(:), & ! (out)
1938  sh25(:), f_shpgh(:), & ! (out)
1939  f_gamma(:), & ! (out)
1940  tltv25(:), qwtv25(:), & ! (out)
1941  tvsq25(:), & ! (out)
1942  tvsq_up(:), tvsq_lo(:) ) ! (out)
1943 
1944  do k = ks, ke_pbl
1945  lq(k) = l(k) * q(k)
1946  nu_f(k) = lq(k) * sm25(k)
1947  kh_f(k) = lq(k) * sh25(k)
1948  end do
1949  if ( atmos_phy_bl_mynn_similarity ) then
1950  nu_f(ks) = karman * z(ks) * us / phi_m
1951  kh_f(ks) = karman * z(ks) * us / phi_h
1952  end if
1953 
1954  do k = ks, ke_pbl-1
1955  nu(k) = min( f2h(k,1) * nu_f(k+1) + f2h(k,2) * nu_f(k), &
1956  atmos_phy_bl_mynn_nu_max )
1957  kh(k) = min( f2h(k,1) * kh_f(k+1) + f2h(k,2) * kh_f(k), &
1958  atmos_phy_bl_mynn_kh_max )
1959  end do
1960 
1961  do k = ks, ke_pbl-1
1962  sw = 0.5_rp - sign(0.5_rp, abs(kh(k)) - eps)
1963  pr(k) = nu(k) / ( kh(k) + sw ) * ( 1.0_rp - sw ) &
1964  + 1.0_rp * sw
1965  end do
1966 
1967  rho(ks) = dens(ks) + dt * sflx_qv / cdz(ks)
1968  do k = ks+1, ke_pbl
1969  rho(k) = dens(k)
1970  end do
1971 
1972  do k = ks, ke_pbl-1
1973  rho_h(k) = f2h(k,1) * rho(k+1) + f2h(k,2) * rho(k)
1974  rhonu(k) = max( nu(k) * rho_h(k), 1e-20_rp )
1975  rhokh(k) = max( kh(k) * rho_h(k), 1e-20_rp )
1976 ! RHONu(k) = F2H(k,1) * RHO(k+1) * Nu_f(k+1) + F2H(k,2) * RHO(k) * Nu_f(k)
1977 ! RHOKh(k) = F2H(k,1) * RHO(k+1) * Kh_f(k+1) + F2H(k,2) * RHO(k) * Kh_f(k)
1978  end do
1979 
1980 
1981  if ( mynn_level3 ) then
1982 
1983  ! production term calculated by explicit scheme
1984  do k = ks, ke_pbl
1985  tltv = betat(k) * tsq(k) + betaq(k) * cov(k)
1986  qwtv = betat(k) * cov(k) + betaq(k) * qsq(k)
1987  gammat(k) = f_gamma(k) * ( tltv - tltv25(k) )
1988  gammaq(k) = f_gamma(k) * ( qwtv - qwtv25(k) )
1989 
1990  wtl = - lq(k) * ( sh25(k) * dtldz(k) + gammat(k) )
1991  wqw = - lq(k) * ( sh25(k) * dqwdz(k) + gammaq(k) )
1992  prod_t(k) = - 2.0_rp * wtl * dtldz(k)
1993  prod_q(k) = - 2.0_rp * wqw * dqwdz(k)
1994  prod_c(k) = - wtl * dqwdz(k) - wqw * dtldz(k)
1995  end do
1996 
1997  if ( atmos_phy_bl_mynn_similarity .or. initialize ) then
1998  tmp = 2.0_rp * us * phi_h / ( karman * z(ks) )
1999  tmp = tmp * ( zeta / ( z(ks) * rlmo ) )**2 ! correspoinding to the limitter for zeta
2000  ! TSQ
2001  prod_t(ks) = tmp * ts**2
2002  ! QSQ
2003  prod_q(ks) = tmp * ts * qs
2004  ! COV
2005  prod_c(ks) = tmp * qs**2
2006  end if
2007 
2008  ! diffusion (explicit)
2009  flx(ks-1) = 0.0_rp
2010  flx(ke_pbl) = 0.0_rp
2011  do k = ks, ke_pbl-1
2012  flx(k) = rhonu(k) * ( tsq(k+1) - tsq(k) ) / fdz(k)
2013  end do
2014  do k = ks, ke_pbl
2015  prod_t(k) = prod_t(k) + ( flx(k) - flx(k-1) ) / cdz(k)
2016  end do
2017  do k = ks, ke_pbl-1
2018  flx(k) = rhonu(k) * ( qsq(k+1) - qsq(k) ) / fdz(k)
2019  end do
2020  do k = ks, ke_pbl
2021  prod_q(k) = prod_q(k) + ( flx(k) - flx(k-1) ) / cdz(k)
2022  end do
2023  do k = ks, ke_pbl-1
2024  flx(k) = rhonu(k) * ( cov(k+1) - cov(k) ) / fdz(k)
2025  end do
2026  do k = ks, ke_pbl
2027  prod_c(k) = prod_c(k) + ( flx(k) - flx(k-1) ) / cdz(k)
2028  end do
2029 
2030  call get_gamma_implicit( &
2031  ka, ks, ke_pbl, &
2032  i, j, &
2033  tsq(:), qsq(:), cov(:), & ! (in)
2034  dtldz(:), dqwdz(:), potv(:), & ! (in)
2035  prod_t(:), prod_q(:), prod_c(:), & ! (in)
2036  betat(:), betaq(:), & ! (in)
2037  f_gamma(:), l(:), q(:), & ! (in)
2038  dt, & ! (in)
2039  dtsq(:), dqsq(:), dcov(:) ) ! (out)
2040 
2041  ! update
2042  do k = ks, ke_pbl
2043  tltv = betat(k) * ( tsq(k) + dtsq(k) ) + betaq(k) * ( cov(k) + dcov(k) )
2044  qwtv = betat(k) * ( cov(k) + dcov(k) ) + betaq(k) * ( qsq(k) + dqsq(k) )
2045  gammat(k) = f_gamma(k) * ( tltv - tltv25(k) )
2046  gammaq(k) = f_gamma(k) * ( qwtv - qwtv25(k) )
2047 
2048  tvsq = max( betat(k) * tltv + betaq(k) * qwtv, 0.0_rp )
2049  tvsq = tvsq - tvsq25(k)
2050  tvsq = min( max( tvsq, tvsq_lo(k) ), tvsq_up(k) )
2051  smp(k) = f_smp(k) * tvsq
2052  shpgh(k) = f_shpgh(k) * tvsq
2053 
2054  wtl = - lq(k) * ( sh25(k) * dtldz(k) + gammat(k) )
2055  wqw = - lq(k) * ( sh25(k) * dqwdz(k) + gammaq(k) )
2056  prod_t(k) = - 2.0_rp * wtl * dtldz(k)
2057  prod_q(k) = - 2.0_rp * wqw * dqwdz(k)
2058  prod_c(k) = - wtl * dqwdz(k) - wqw * dtldz(k)
2059  end do
2060 
2061  if ( atmos_phy_bl_mynn_similarity .or. initialize ) then
2062  smp(ks) = 0.0_rp
2063  shpgh(ks) = 0.0_rp
2064  gammat(ks) = 0.0_rp
2065  gammaq(ks) = 0.0_rp
2066 
2067  tmp = 2.0_rp * us * phi_h / ( karman * z(ks) )
2068  tmp = tmp * ( zeta / ( z(ks) * rlmo ) )**2 ! correspoinding to the limitter for zeta
2069  ! TSQ
2070  prod_t(ks) = tmp * ts**2
2071  ! QSQ
2072  prod_q(ks) = tmp * ts * qs
2073  ! COV
2074  prod_c(ks) = tmp * qs**2
2075  end if
2076 
2077  else
2078 
2079  do k = ks, ke_pbl
2080  smp(k) = 0.0_rp
2081  shpgh(k) = 0.0_rp
2082  gammat(k) = 0.0_rp
2083  gammaq(k) = 0.0_rp
2084  end do
2085 
2086  end if
2087 
2088 
2089  ! production of TKE
2090  do k = ks, ke_pbl
2091  prod(k) = lq(k) * ( ( sm25(k) + smp(k) ) * dudz2(k) &
2092  - ( sh25(k) * n2_new(k) - shpgh(k) ) )
2093  end do
2094  if ( atmos_phy_bl_mynn_similarity .or. initialize ) then
2095  prod(ks) = us3 / ( karman * z(ks) ) * ( phi_m - zeta )
2096  end if
2097 
2098  do k = ks, ke_pbl
2099  diss(k) = - 2.0_rp * q(k) / ( b1 * l(k) )
2100 ! prod(k) = max( prod(k), - tke(k) / dt - diss(k) * tke(k) )
2101  diss_p(k) = dt * 2.0_rp * q(k) / ( b2 * l(k) )
2102  end do
2103 
2104 
2105  if ( .not. initialize ) exit
2106 
2107  ! dens * TKE
2108 
2109 !!$ if ( ATMOS_PHY_BL_MYNN_MF ) then
2110 !!$ do k = KS, KE_PBL-1
2111 !!$ flx(k) = eflux(k)
2112 !!$ end do
2113 !!$ else
2114  do k = ks, ke_pbl-1
2115  flx(k) = 0.0_rp
2116  end do
2117 !!$ end if
2118 
2119  do k = ks, ke_pbl-1
2120  d(k) = dt * ( - ( flx(k) - flx(k-1) ) / ( cdz(k) * rho(k) ) + prod(k) )
2121  end do
2122  d(ke_pbl) = 0.0_rp
2123 
2124  c(ks) = 0.0_rp
2125  do k = ks, ke_pbl-1
2126  ap = - dt * atmos_phy_bl_mynn_sq_fact * rhonu(k) / fdz(k)
2127  a(k) = ap / ( rho(k) * cdz(k) )
2128 #ifdef _OPENACC
2129  if ( k==ks ) then
2130  b(k) = - a(k) - diss(k) * dt
2131  else
2132  b(k) = - a(k) + dt * atmos_phy_bl_mynn_sq_fact * rhonu(k-1) / ( fdz(k-1) * rho(k) * cdz(k) ) - diss(k) * dt
2133  end if
2134 #else
2135  b(k) = - a(k) - c(k) - diss(k) * dt
2136 #endif
2137  c(k+1) = ap / ( rho(k+1) * cdz(k+1) )
2138  end do
2139  a(ke_pbl) = 0.0_rp
2140  b(ke_pbl) = - c(ke_pbl) - diss(ke_pbl) * dt
2141 
2142  call matrix_solver_tridiagonal( &
2143  ka, ks, ke_pbl, &
2144 #ifdef _OPENACC
2145  work(:,:), & ! (work)
2146 #endif
2147  a(:), b(:), c(:), d(:), & ! (in)
2148  tke(:) ) ! (out)
2149 
2150  do k = ks, ke_pbl-1
2151  tke(k) = min( max( tke(k), atmos_phy_bl_mynn_tke_min ), 100.0_rp )
2152  end do
2153  tke(ke_pbl) = 0.0_rp
2154 
2155 
2156  do k = ks, ke_pbl
2157  q(k) = ( q(k) + sqrt( tke(k) * 2.0_rp ) ) * 0.5_rp ! to avoid oscillation
2158  end do
2159 
2160 
2161  if ( .not. mynn_level3 ) cycle
2162 
2163  ! dens * tsq
2164 
2165  do k = ks, ke_pbl-1
2166  d(k) = dt * prod_t(k)
2167  end do
2168  d(ke_pbl) = 0.0_rp
2169  c(ks) = 0.0_rp
2170  do k = ks, ke_pbl-1
2171  ap = - dt * rhonu(k) / fdz(k)
2172  a(k) = ap / ( rho(k) * cdz(k) )
2173 #ifdef _OPENACC
2174  if ( k==ks ) then
2175  b(k) = - a(k) + diss_p(k)
2176  else
2177  b(k) = - a(k) + dt * rhonu(k-1) / ( fdz(k-1) * rho(k) * cdz(k) ) + diss_p(k)
2178  end if
2179 #else
2180  b(k) = - a(k) - c(k) + 1.0_rp + diss_p(k)
2181 #endif
2182  c(k+1) = ap / ( rho(k+1) * cdz(k+1) )
2183  end do
2184  a(ke_pbl) = 0.0_rp
2185  b(ke_pbl) = - c(ke_pbl) + diss_p(ke_pbl)
2186 
2187  call matrix_solver_tridiagonal( &
2188  ka, ks, ke_pbl, &
2189 #ifdef _OPENACC
2190  work(:,:), & ! (work)
2191 #endif
2192  a(:), b(:), c(:), d(:), & ! (in)
2193  tsq(:) ) ! (out)
2194 
2195  do k = ks, ke_pbl-1
2196  tsq(k) = max( tsq(k), 0.0_rp )
2197  end do
2198  tsq(ke_pbl) = 0.0_rp
2199 
2200 
2201  ! dens * qsq
2202 
2203  do k = ks, ke_pbl-1
2204  d(k) = dt * prod_q(k)
2205  end do
2206  d(ke_pbl) = 0.0_rp
2207  ! a, b, c are same as those for tsq
2208 
2209  call matrix_solver_tridiagonal( &
2210  ka, ks, ke_pbl, &
2211 #ifdef _OPENACC
2212  work(:,:), & ! (work)
2213 #endif
2214  a(:), b(:), c(:), d(:), & ! (in)
2215  qsq(:) ) ! (out)
2216 
2217  do k = ks, ke_pbl-1
2218  qsq(k) = max( qsq(k), 0.0_rp )
2219  end do
2220  qsq(ke_pbl) = 0.0_rp
2221 
2222 
2223  ! dens * cov
2224 
2225  do k = ks, ke_pbl-1
2226  d(k) = dt * prod_c(k)
2227  end do
2228  d(ke_pbl) = 0.0_rp
2229  ! a, b, c are same as those for tsq
2230 
2231  call matrix_solver_tridiagonal( &
2232  ka, ks, ke_pbl, &
2233 #ifdef _OPENACC
2234  work(:,:), & ! (work)
2235 #endif
2236  a(:), b(:), c(:), d(:), & ! (in)
2237  cov(:) ) ! (out)
2238 
2239  do k = ks, ke_pbl-1
2240  cov(k) = sign( min( abs(cov(k)), sqrt(tsq(k)*qsq(k)) ), cov(k) )
2241  end do
2242  cov(ke_pbl) = 0.0_rp
2243 
2244  end do
2245 
2246  return
2247  end subroutine mynn_main
2248 
2249 !OCL SERIAL
2250  subroutine get_length( &
2251  KA, KS, KE_PBL, &
2252  q, n2, &
2253  POTV, DENS, &
2254  mflux_gl, &
2255  Zi, &
2256  SFLX_BUOY, RLmo, &
2257  z, CDZ, FDZ, &
2258  initialize, &
2259  l )
2260  !$acc routine vector
2261  use scale_const, only: &
2262 ! GRAV => CONST_GRAV, &
2263  karman => const_karman, &
2264  eps => const_eps
2265  implicit none
2266  integer, intent(in), value :: KA, KS, KE_PBL
2267 
2268  real(RP), intent(in) :: q(KA)
2269  real(RP), intent(in) :: n2(KA)
2270  real(RP), intent(in) :: POTV(KA)
2271  real(RP), intent(in) :: DENS(KA)
2272  real(RP), intent(in) :: mflux_gl(0:KA)
2273  real(RP), intent(in) :: Zi
2274  real(RP), intent(in) :: SFLX_BUOY
2275  real(RP), intent(in) :: RLmo
2276  real(RP), intent(in) :: z(KA)
2277  real(RP), intent(in) :: CDZ(KA)
2278  real(RP), intent(in) :: FDZ(KA)
2279  logical, intent(in) :: initialize
2280 
2281  real(RP), intent(out) :: l(KA)
2282 
2283  real(RP), parameter :: sqrt05 = sqrt( 0.5_rp )
2284 
2285  real(RP) :: ls
2286  real(RP) :: lb
2287  real(RP) :: lt
2288  real(RP) :: rlt
2289 
2290  real(RP) :: qc
2291  real(RP) :: int_q
2292  real(RP) :: int_qz
2293  real(RP) :: rn2sr
2294  real(RP) :: zeta
2295 
2296  real(RP) :: qdz
2297 
2298  ! for Olson et al. (2019)
2299 ! real(RP) :: lbl, lup, ldown
2300  real(RP) :: tau
2301 
2302  real(RP) :: sw
2303  integer :: kmax
2304  integer :: k
2305 
2306  if ( atmos_phy_bl_mynn_use_zi ) then
2307  kmax = ke_pbl
2308  !$acc loop seq
2309  do k = ks, ke_pbl
2310  if ( z(k) > zi * 1.3_rp ) then
2311  kmax = k - 1
2312  exit
2313  end if
2314  end do
2315  kmax = max(kmax, ks)
2316  else
2317  kmax = ke_pbl
2318  end if
2319 
2320  int_qz = 0.0_rp
2321  int_q = 0.0_rp
2322  !$acc loop private(qdz) reduction(+:int_qz,int_q)
2323  do k = ks, kmax
2324  qdz = q(k) * cdz(k)
2325  int_qz = int_qz + z(k) * qdz
2326  int_q = int_q + qdz
2327  end do
2328  ! LT
2329  lt = min( max(0.23_rp * int_qz / (int_q + 1e-20_rp), &
2330  lt_min), &
2331  atmos_phy_bl_mynn_lt_max )
2332  rlt = 1.0_rp / lt
2333 
2334  qc = ( lt * max(sflx_buoy,0.0_rp) )**oneoverthree ! qc=0 if SFLX_BUOY<0
2335  if ( atmos_phy_bl_mynn_o2019 ) then
2336  tau = 0.5_rp * zi / max(sflx_buoy,1.0e-10_rp)**oneoverthree
2337  end if
2338 
2339  !$acc loop private(zeta,sw,ls,rn2sr,tau,lb)
2340  do k = ks, ke_pbl
2341  zeta = z(k) * rlmo
2342 
2343  ! LS
2344  sw = sign(0.5_rp, zeta) + 0.5_rp ! 1 for zeta >= 0, 0 for zeta < 0
2345  ls = karman * z(k) &
2346  * ( sw / (1.0_rp + atmos_phy_bl_mynn_cns*zeta*sw ) &
2347  + ( (1.0_rp - atmos_phy_bl_mynn_alpha4*zeta)*(1.0_rp-sw) )**0.2_rp )
2348 
2349  ! LB
2350  sw = sign(0.5_rp, n2(k)-eps) + 0.5_rp ! 1 for dptdz >0, 0 for dptdz <= 0
2351  rn2sr = 1.0_rp / ( sqrt(n2(k)*sw) + 1.0_rp-sw)
2352  if ( atmos_phy_bl_mynn_o2019 ) then
2353  if ( z(k) > zi ) tau = 50.0_rp
2354 !!$ qtmp = q(k)**2 * 0.5_RP * POTV(k) / GRAV
2355 !!$ lup = z(KE_PBL) - z(k)
2356 !!$ do kk = k+1, KE_PBL
2357 !!$ qtmp = qtmp - ( POTV(kk) - POTV(k) ) * FDZ(kk-1)
2358 !!$ if ( qtmp < 0.0_RP ) then
2359 !!$ lup = z(kk) - z(k)
2360 !!$ exit
2361 !!$ end if
2362 !!$ end do
2363 !!$ qtmp = q(k)**2 * 0.5_RP * POTV(k) / GRAV
2364 !!$ ldown = z(k)
2365 !!$ do kk = k-1, KS
2366 !!$ qtmp = qtmp - ( POTV(k) - POTV(kk) ) * FDZ(kk)
2367 !!$ if ( qtmp < 0.0_RP ) then
2368 !!$ ldown = z(k) - z(kk)
2369 !!$ exit
2370 !!$ end if
2371 !!$ end do
2372 !!$ lbl = sqrt( lup**2 + ldown**2 )
2373 !!$ we = 0.5_RP * tanh( ( z(k) - zi * 1.3_RP ) / ( zi * 0.15_RP ) ) + 0.5_RP
2374 !!$ lb = lb * ( 1.0_RP - we ) + lbl * we
2375  lb = atmos_phy_bl_mynn_alpha2 * max( q(k), ( mflux_gl(k)+mflux_gl(k-1))/dens(k) ) * rn2sr * sw &
2376  + tau * q(k) * sqrt05 * (1.0_rp-sw)
2377  else
2378  lb = atmos_phy_bl_mynn_alpha2 * (1.0_rp + 5.0_rp * sqrt(qc*rn2sr*rlt)) * q(k) * rn2sr * sw & ! qc=0 when RLmo > 0
2379  + 1.e10_rp * (1.0_rp-sw)
2380  end if
2381 
2382  ! L
2383  if ( initialize ) then
2384  l(k) = min( ls, lb )
2385  else if ( atmos_phy_bl_mynn_o2019 ) then
2386  l(k) = min( 1.0_rp / ( 1.0_rp/ls + rlt ), lb )
2387  else
2388  l(k) = 1.0_rp / ( 1.0_rp/ls + rlt + 1.0_rp/(lb+1e-20_rp) )
2389  end if
2390  end do
2391 
2392  return
2393  end subroutine get_length
2394 
2395 !OCL SERIAL
2396  subroutine get_q2_level2( &
2397  KA, KS, KE_PBL, &
2398  dudz2, Ri, l, &
2399  q2_2 )
2400  !$acc routine vector
2401  implicit none
2402  integer, intent(in), value :: KA, KS, KE_PBL
2403 
2404  real(RP), intent(in) :: dudz2(KA)
2405  real(RP), intent(in) :: Ri(KA)
2406  real(RP), intent(in) :: l(KA)
2407 
2408  real(RP), intent(out) :: q2_2(KA)
2409 
2410  real(RP) :: rf
2411  real(RP) :: sm_2
2412  real(RP) :: sh_2
2413 
2414  ! for K2010
2415  real(RP) :: A2k, F1, Rf1, AF12
2416 
2417  integer :: k
2418 
2419  do k = ks, ke_pbl
2420  if ( atmos_phy_bl_mynn_k2010 .and. ri(k) > 0.0_rp ) then
2421  a2k = a2 / ( 1.0_rp + ri(k) )
2422  else
2423  a2k = a2
2424  end if
2425  f1 = b1 * ( g1 - c1 ) + 2.0 * a1 * ( 3.0_rp - 2.0_rp * c2 ) + 3.0 * a2k * ( 1.0_rp - c2 ) * ( 1.0_rp - c5 )
2426  rf1 = b1 * ( g1 - c1 ) / f1
2427  af12 = a1 * f1 / ( a2k * f2 )
2428  rf = min(0.5_rp / af12 * ( ri(k) &
2429  + af12*rf1 &
2430  - sqrt(ri(k)**2 + 2.0_rp*af12*(rf1-2.0_rp*rf2)*ri(k) + (af12*rf1)**2) ), &
2431  rfc)
2432  sh_2 = 3.0_rp * a2k * (g1+g2) * (rfc-rf) / (1.0_rp-rf)
2433  sm_2 = sh_2 * af12 * (rf1-rf) / (rf2-rf)
2434  q2_2(k) = b1 * l(k)**2 * sm_2 * (1.0_rp-rf) * dudz2(k)
2435  end do
2436 
2437  return
2438  end subroutine get_q2_level2
2439 
2440 !OCL SERIAL
2441  subroutine get_smsh( &
2442  KA, KS, KE_PBL, &
2443  i, j, &
2444  q, ac, &
2445  l, n2, Ri, &
2446  potv, dudz2, &
2447  dtldz, dqwdz, &
2448  betat, betaq, &
2449  mynn_level3, &
2450  initialize, &
2451  tsq, qsq, cov, &
2452  sm25, f_smp, &
2453  sh25, f_shpgh, &
2454  f_gamma, &
2455  tltv25, qwtv25, &
2456  tvsq25, tvsq_up, &
2457  tvsq_lo )
2458  !$acc routine vector
2459  use scale_const, only: &
2460  eps => const_eps, &
2461  huge => const_huge, &
2462  grav => const_grav
2463  implicit none
2464  integer, intent(in), value :: KA, KS, KE_PBL
2465  integer, intent(in), value :: i, j ! for debug
2466 
2467  real(RP), intent(in) :: q(KA)
2468  real(RP), intent(in) :: ac(KA)
2469  real(RP), intent(in) :: l(KA)
2470  real(RP), intent(in) :: n2(KA)
2471  real(RP), intent(in) :: Ri(KA)
2472  real(RP), intent(in) :: potv(KA)
2473  real(RP), intent(in) :: dudz2(KA)
2474  real(RP), intent(in) :: dtldz(KA)
2475  real(RP), intent(in) :: dqwdz(KA)
2476  real(RP), intent(in) :: betat(KA)
2477  real(RP), intent(in) :: betaq(KA)
2478  logical, intent(in), value :: mynn_level3
2479  logical, intent(in), value :: initialize
2480 
2481  real(RP), intent(inout) :: tsq(KA)
2482  real(RP), intent(inout) :: qsq(KA)
2483  real(RP), intent(inout) :: cov(KA)
2484 
2485  real(RP), intent(out) :: sm25 (KA) ! S_M2.5
2486  real(RP), intent(out) :: f_smp (KA) ! S_M' / ( <\theta_v^2> - <\theta_v^2>2.5 )
2487  real(RP), intent(out) :: sh25 (KA) ! S_H2.5
2488  real(RP), intent(out) :: f_shpgh(KA) ! S_H' G_H * (q/L)^2 / ( <\theta_v^2> - <\theta_v^2>2.5 )
2489  real(RP), intent(out) :: f_gamma(KA) ! - E_H / q^2 * GRAV / Theta_0
2490  real(RP), intent(out) :: tltv25 (KA)
2491  real(RP), intent(out) :: qwtv25 (KA)
2492  real(RP), intent(out) :: tvsq25 (KA)
2493  real(RP), intent(out) :: tvsq_up(KA)
2494  real(RP), intent(out) :: tvsq_lo(KA)
2495 
2496  real(RP) :: l2
2497  real(RP) :: q2
2498  real(RP) :: l2q2
2499  real(RP) :: ac2
2500  real(RP) :: p1q2
2501  real(RP) :: p2q2
2502  real(RP) :: p3q2
2503  real(RP) :: p4q2
2504  real(RP) :: p5q2
2505  real(RP) :: rd25q2
2506  real(RP) :: ghq2
2507  real(RP) :: gmq2
2508  real(RP) :: f1, f2, f3, f4
2509 
2510  ! for level 3
2511  real(RP) :: tvsq
2512  real(RP) :: tltv
2513  real(RP) :: qwtv
2514  real(RP) :: tsq25
2515  real(RP) :: qsq25
2516  real(RP) :: cov25
2517  real(RP) :: emq2
2518  real(RP) :: eh
2519  real(RP) :: ew
2520  real(RP) :: rdpq2
2521  real(RP) :: cw25
2522  real(RP) :: fact
2523  real(RP) :: tmp1, tmp2
2524 
2525  ! for K2010
2526  real(RP) :: A2k
2527 
2528  integer :: k
2529 
2530  ! level 2.5
2531  do k = ks, ke_pbl
2532 
2533  if ( atmos_phy_bl_mynn_k2010 .and. ri(k) > 0.0_rp ) then
2534  a2k = a2 / ( 1.0_rp + ri(k) )
2535  else
2536  a2k = a2
2537  end if
2538 
2539  ac2 = ac(k)**2
2540  l2 = l(k)**2
2541  q2 = q(k)**2
2542 
2543  f1 = - 3.0_rp * ac2 * a2k * b2 * ( 1.0_rp - c3 )
2544  f2 = - 9.0_rp * ac2 * a1 * a2k * ( 1.0_rp - c2 )
2545  f3 = 9.0_rp * ac2 * a2k**2 * ( 1.0_rp - c2 ) * ( 1.0_rp - c5 )
2546  f4 = - 12.0_rp * ac2 * a1 * a2k * ( 1.0_rp - c2 )
2547 
2548  if ( mynn_level3 ) then ! level 3
2549  ghq2 = max( - n2(k) * l2, -q2 ) ! L/q <= 1/N for N^2>0
2550  else
2551  ghq2 = - n2(k) * l2
2552  end if
2553  gmq2 = dudz2(k) * l2
2554 
2555  p1q2 = q2 + f1 * ghq2
2556  p2q2 = q2 + f2 * ghq2
2557  p3q2 = q2 + ( f1 + f3 ) * ghq2
2558  p4q2 = q2 + ( f1 + f4 ) * ghq2
2559  p5q2 = 6.0_rp * ac2 * a1**2 * gmq2
2560  rd25q2 = q2 / max( p2q2 * p4q2 + p5q2 * p3q2, 1e-20_rp )
2561 
2562  sm25(k) = ac(k) * a1 * ( p3q2 - 3.0_rp * c1 * p4q2 ) * rd25q2
2563  sh25(k) = ac(k) * a2k * ( p2q2 + 3.0_rp * c1 * p5q2 ) * rd25q2
2564 
2565  fact = ac(k) * b2 * l2 * sh25(k)
2566  tsq25 = fact * dtldz(k)**2
2567  qsq25 = fact * dqwdz(k)**2
2568  cov25 = fact * dtldz(k) * dqwdz(k)
2569 
2570  if ( initialize .or. (.not. mynn_level3 ) ) then
2571  tsq(k) = tsq25
2572  qsq(k) = qsq25
2573  cov(k) = cov25
2574  end if
2575 
2576  if ( mynn_level3 ) then ! level 3
2577 
2578  if ( q2 <= 1e-10_rp ) then
2579  tltv25(k) = 0.0_rp
2580  qwtv25(k) = 0.0_rp
2581  tvsq25(k) = 0.0_rp
2582  tvsq_up(k) = 0.0_rp
2583  tvsq_lo(k) = 0.0_rp
2584  f_smp(k) = 0.0_rp
2585  f_shpgh(k) = 0.0_rp
2586  f_gamma(k) = 0.0_rp
2587  else
2588  tltv25(k) = betat(k) * tsq25 + betaq(k) * cov25
2589  qwtv25(k) = betat(k) * cov25 + betaq(k) * qsq25
2590  tvsq25(k) = max(betat(k) * tltv25(k) + betaq(k) * qwtv25(k), 0.0_rp)
2591  cw25 = p1q2 * ( p2q2 + 3.0_rp * c1 * p5q2 ) * rd25q2 / ( 3.0_rp * q2 )
2592 
2593  l2q2 = l2 / q2
2594  l2q2 = min( l2q2, 1.0_rp / max(n2(k), eps) )
2595  q2 = l2 / l2q2
2596 
2597  rdpq2 = q2 / max( p2q2 * ( f4 * ghq2 + q2 ) + p5q2 * ( f3 * ghq2 + q2 ), 1e-20_rp )
2598 
2599  tltv = betat(k) * tsq(k) + betaq(k) * cov(k)
2600  qwtv = betat(k) * cov(k) + betaq(k) * qsq(k)
2601  tvsq = max(betat(k) * tltv + betaq(k) * qwtv, 0.0_rp)
2602 
2603  ew = ( 1.0_rp - c3 ) * ( - p2q2 * f4 - p5q2 * f3 ) * rdpq2 ! (p1-p4)/gh = -f4, (p1-p3)/gh = -f3
2604  if ( abs(ew) > eps ) then
2605  fact = q2 * potv(k)**2 / ( ew * l2q2 * grav**2 )
2606  if ( fact > 0.0_rp ) then
2607  tmp1 = 0.76_rp
2608  tmp2 = 0.12_rp
2609  else
2610  tmp1 = 0.12_rp
2611  tmp2 = 0.76_rp
2612  end if
2613  tvsq_up(k) = ( tmp1 - cw25 ) * fact
2614  tvsq_lo(k) = ( tmp2 - cw25 ) * fact
2615  else
2616  tvsq_up(k) = huge
2617  tvsq_lo(k) = -huge
2618  end if
2619 
2620  emq2 = 3.0_rp * ac(k) * a1 * ( 1.0_rp - c3 ) * ( f3 - f4 ) * rdpq2 ! (p3-p4)/gh = (f3-f4)
2621  eh = 3.0_rp * ac(k) * a2k * ( 1.0_rp - c3 ) * ( p2q2 + p5q2 ) * rdpq2
2622 
2623 ! q2 = l2 / l2q2
2624  fact = grav / potv(k)
2625  f_smp(k) = emq2 * fact**2 * l2q2
2626  f_shpgh(k) = eh * fact**2 / q2
2627  f_gamma(k) = - eh * fact / q2
2628  end if
2629  else ! level 2.5
2630  tvsq_up(k) = 0.0_rp
2631  tvsq_lo(k) = 0.0_rp
2632  f_smp(k) = 0.0_rp
2633  f_shpgh(k) = 0.0_rp
2634  f_gamma(k) = 0.0_rp
2635  end if
2636 
2637  end do
2638 
2639  return
2640  end subroutine get_smsh
2641 
2642 !OCL SERIAL
2643  subroutine get_gamma_implicit( &
2644  KA, KS, KE, &
2645  i, j, &
2646  tsq, qsq, cov, &
2647  dtldz, dqwdz, POTV, &
2648  prod_t, prod_q, prod_c, &
2649  betat, betaq, &
2650  f_gamma, l, q, &
2651  dt, &
2652  dtsq, dqsq, dcov )
2653  !$acc routine vector
2654  use scale_const, only: &
2655  grav => const_grav
2656  implicit none
2657  integer, intent(in), value :: KA, KS, KE
2658  integer, intent(in), value :: i, j ! for debug
2659 
2660  real(RP), intent(in) :: tsq (KA)
2661  real(RP), intent(in) :: qsq (KA)
2662  real(RP), intent(in) :: cov (KA)
2663  real(RP), intent(in) :: dtldz (KA)
2664  real(RP), intent(in) :: dqwdz (KA)
2665  real(RP), intent(in) :: POTV (KA)
2666  real(RP), intent(in) :: prod_t (KA)
2667  real(RP), intent(in) :: prod_q (KA)
2668  real(RP), intent(in) :: prod_c (KA)
2669  real(RP), intent(in) :: betat (KA)
2670  real(RP), intent(in) :: betaq (KA)
2671  real(RP), intent(in) :: f_gamma(KA)
2672  real(RP), intent(in) :: l (KA)
2673  real(RP), intent(in) :: q (KA)
2674  real(RP), intent(in), value :: dt
2675 
2676  real(RP), intent(out) :: dtsq(KA)
2677  real(RP), intent(out) :: dqsq(KA)
2678  real(RP), intent(out) :: dcov(KA)
2679 
2680  real(RP) :: a11, a12, a21, a22, a23, a32, a33
2681  real(RP) :: v1, v2, v3
2682  real(RP) :: f1, f2
2683  real(RP) :: a13, a31, det
2684 
2685  integer :: k
2686 
2687  ! calculate gamma by implicit scheme
2688 
2689  do k = ks, ke
2690 
2691  ! matrix coefficient
2692  f1 = f_gamma(k) * l(k) * q(k)
2693  f2 = - 2.0_rp * q(k) / ( b2 * l(k) )
2694 
2695  v1 = prod_t(k) + f2 * tsq(k)
2696  v2 = prod_c(k) + f2 * cov(k)
2697  v3 = prod_q(k) + f2 * qsq(k)
2698 
2699  a11 = - f1 * dtldz(k) * betat(k) * 2.0_rp + 1.0_rp / dt - f2
2700  a12 = - f1 * dtldz(k) * betaq(k) * 2.0_rp
2701  a21 = - f1 * dqwdz(k) * betat(k)
2702  a22 = - f1 * ( dtldz(k) * betat(k) + dqwdz(k) * betaq(k) ) + 1.0_rp / dt - f2
2703  a23 = - f1 * dtldz(k) * betaq(k)
2704  a32 = - f1 * dqwdz(k) * betat(k) * 2.0_rp
2705  a33 = - f1 * dqwdz(k) * betaq(k) * 2.0_rp + 1.0_rp / dt - f2
2706 
2707  if ( q(k) < 1e-20_rp .or. abs(f_gamma(k)) * dt >= q(k) ) then
2708  ! solve matrix
2709  f1 = a21 / a11
2710  f2 = a23 / a33
2711  dcov(k) = ( v2 - f1 * v1 - f2 * v3 ) / ( a22 - f1 * a12 - f2 * a32 )
2712  dtsq(k) = ( v1 - a12 * dcov(k) ) / a11
2713  dqsq(k) = ( v3 - a32 * dcov(k) ) / a33
2714  else
2715  ! consider change in q
2716  ! dq = - L * G / Theta0 * ( betat * dgammat + betaq * dgammaq )
2717  f1 = - l(k) * grav / potv(k) * f_gamma(k) / q(k)
2718  f2 = f1 * betat(k)**2
2719  a11 = a11 - f2 * v1
2720  a21 = a21 - f2 * v2
2721  a31 = - f2 * v3
2722  f2 = f1 * betat(k) * betaq(k) * 2.0_rp
2723  a12 = a12 - f2 * v1
2724  a22 = a22 - f2 * v2
2725  a32 = a32 - f2 * v3
2726  f2 = f1 * betaq(k)**2
2727  a13 = - f2 * v1
2728  a23 = a23 - f2 * v2
2729  a33 = a33 - f2 * v3
2730  ! solve matrix by Cramer's fomula
2731  det = a11 * ( a22 * a33 - a23 * a32 ) + a12 * ( a23 * a31 - a21 * a33 ) + a13 * ( a21 * a32 - a22 * a31 )
2732  dtsq(k) = ( v1 * ( a22 * a33 - a23 * a32 ) + a12 * ( a23 * v3 - v2 * a33 ) + a13 * ( v2 * a32 - a22 * v3 ) ) / det
2733  dcov(k) = ( a11 * ( v2 * a33 - a23 * v3 ) + v1 * ( a23 * a31 - a21 * a33 ) + a13 * ( a21 * v3 - v2 * a31 ) ) / det
2734  dqsq(k) = ( a11 * ( a22 * v3 - v2 * a32 ) + a12 * ( v2 * a31 - a21 * v3 ) + v3 * ( a21 * a32 - a22 * a31 ) ) / det
2735  end if
2736 
2737  end do
2738 
2739  return

References calc_vertical_differece(), scale_const::const_cpdry, scale_const::const_eps, scale_const::const_epstvap, scale_const::const_grav, scale_const::const_huge, scale_const::const_karman, scale_atmos_hydrometeor::cp_vapor, get_phi(), 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::kmax, scale_atmos_grid_cartesc_index::ks, partial_condensation(), scale_prc::prc_abort(), and scale_precision::rp.

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:

◆ calc_vertical_differece()

subroutine scale_atmos_phy_bl_mynn::calc_vertical_differece ( integer, intent(in), value  KA,
integer, intent(in), value  KS,
integer, intent(in), value  KE,
integer, intent(in), value  i,
integer, intent(in), value  j,
real(rp), dimension(ka), intent(out)  dudz2,
real(rp), dimension(ka), intent(out)  dtldz,
real(rp), dimension(ka), intent(out)  dqwdz,
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)  QDRY,
real(rp), dimension (ka), intent(in)  CDZ,
real(rp), dimension (ka), intent(in)  FDZ,
real(rp), dimension (ka,2), intent(in)  F2H,
real(rp), intent(in)  us,
real(rp), intent(in)  ts,
real(rp), intent(in)  qs,
real(rp), intent(in)  phi_m,
real(rp), intent(in)  phi_h,
real(rp), intent(in)  z1,
real(rp), dimension (ka), intent(out)  Uh,
real(rp), dimension (ka), intent(out)  Vh,
real(rp), dimension(ka), intent(out)  qw2,
real(rp), dimension (ka), intent(out)  qh 
)

Definition at line 2752 of file scale_atmos_phy_bl_mynn.F90.

2752  !$acc routine vector
2753  use scale_const, only: &
2754  karman => const_karman
2755  integer, intent(in), value :: KA, KS, KE
2756  integer, intent(in), value :: i, j ! for debug
2757 
2758  real(RP), intent(out) :: dudz2(KA)
2759  real(RP), intent(out) :: dtldz(KA)
2760  real(RP), intent(out) :: dqwdz(KA)
2761 
2762  real(RP), intent(in) :: U (KA)
2763  real(RP), intent(in) :: V (KA)
2764  real(RP), intent(in) :: POTL(KA)
2765  real(RP), intent(in) :: Qw (KA)
2766  real(RP), intent(in) :: QDRY(KA)
2767  real(RP), intent(in) :: CDZ (KA)
2768  real(RP), intent(in) :: FDZ (KA)
2769  real(RP), intent(in) :: F2H (KA,2)
2770  real(RP), intent(in) :: us
2771  real(RP), intent(in) :: ts
2772  real(RP), intent(in) :: qs
2773  real(RP), intent(in) :: phi_m
2774  real(RP), intent(in) :: phi_h
2775  real(RP), intent(in) :: z1
2776 
2777  real(RP), intent(out) :: Uh (KA) ! work
2778  real(RP), intent(out) :: Vh (KA) ! work
2779  real(RP), intent(out) :: qw2(KA) ! work
2780  real(RP), intent(out) :: qh (KA) ! work
2781 
2782  integer :: k
2783 
2784  do k = ks, ke
2785  uh(k) = f2h(k,1) * u(k+1) + f2h(k,2) * u(k)
2786  vh(k) = f2h(k,1) * v(k+1) + f2h(k,2) * v(k)
2787  end do
2788 
2789  if ( atmos_phy_bl_mynn_dz_sim ) then
2790  dudz2(ks) = ( us * phi_m / ( karman * z1 ) )**2
2791  else
2792  dudz2(ks) = ( ( uh(ks) - u(ks) )**2 + ( vh(ks) - v(ks) )**2 ) / ( cdz(ks) * 0.5_rp )**2
2793 ! dudz2(KS) = ( ( Uh(KS) )**2 + ( Vh(KS) )**2 ) / CDZ(KS)**2
2794  end if
2795  dudz2(ks) = max( dudz2(ks), 1e-20_rp )
2796  do k = ks+1, ke
2797  dudz2(k) = ( ( uh(k) - uh(k-1) )**2 + ( vh(k) - vh(k-1) )**2 ) / cdz(k)**2
2798 ! dudz2(k) = ( &
2799 ! ( U(k+1) * FDZ(k-1)**2 - U(k-1) * FDZ(k)**2 + U(k) * ( FDZ(k)**2 - FDZ(k-1)**2 ) )**2 &
2800 ! + ( V(k+1) * FDZ(k-1)**2 - V(k-1) * FDZ(k)**2 + V(k) * ( FDZ(k)**2 - FDZ(k-1)**2 ) )**2 &
2801 ! ) / ( FDZ(k-1) * FDZ(k) * ( FDZ(k-1) + FDZ(k) ) )**2
2802  dudz2(k) = max( dudz2(k), 1e-20_rp )
2803  end do
2804 
2805 
2806  do k = ks, ke
2807  qh(k) = f2h(k,1) * potl(k+1) + f2h(k,2) * potl(k)
2808  end do
2809  if ( atmos_phy_bl_mynn_dz_sim ) then
2810  dtldz(ks) = ts * phi_h / ( karman * z1 )
2811  else
2812  dtldz(ks) = ( qh(ks) - potl(ks) ) / ( cdz(ks) * 0.5_rp )
2813 ! dtldz(KS) = ( POTL(KS+1) - POTL(KS) ) / FDZ(KS)
2814  end if
2815  do k = ks+1, ke
2816  dtldz(k) = ( qh(k) - qh(k-1) ) / cdz(k)
2817 ! dtldz(k) = ( POTL(k+1) * FDZ(k-1)**2 - POTL(k-1) * FDZ(k)**2 + POTL(k) * ( FDZ(k)**2 - FDZ(k-1)**2 ) ) &
2818 ! / ( FDZ(k-1) * FDZ(k) * ( FDZ(k-1) + FDZ(k) ) )
2819  end do
2820 
2821  do k = ks, ke+1
2822  qw2(k) = qw(k) / ( qdry(k) + qw(k) )
2823  end do
2824  do k = ks, ke
2825 ! qh(k) = f2h(k,1) * Qw(k+1) + f2h(k,2) * Qw(k)
2826  qh(k) = f2h(k,1) * qw2(k+1) + f2h(k,2) * qw2(k)
2827  end do
2828  if ( atmos_phy_bl_mynn_dz_sim ) then
2829  dqwdz(ks) = qs * phi_h / ( karman * z1 )
2830  else
2831  dqwdz(ks) = ( qh(ks) - qw2(ks) ) / ( cdz(ks) * 0.5_rp )
2832 ! dqwdz(KS) = ( qh(KS) - Qw(KS) ) / ( CDZ(KS) * 0.5_RP )
2833 ! dqwdz(KS) = ( Qw(KS+1) - Qw(KS) ) / FDZ(KS)
2834  end if
2835  do k = ks+1, ke
2836  dqwdz(k) = ( qh(k) - qh(k-1) ) / cdz(k)
2837 ! dqwdz(k) = ( Qw(k+1) * FDZ(k-1)**2 - Qw(k-1) * FDZ(k)**2 + Qw(k) * ( FDZ(k)**2 - FDZ(k-1)**2 ) ) &
2838 ! / ( FDZ(k-1) * FDZ(k) * ( FDZ(k-1) + FDZ(k) ) )
2839  end do
2840 
2841  return

References scale_const::const_karman.

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), value  KA,
integer, intent(in), value  KS,
integer, intent(in), value  KE,
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,
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)  EXNER,
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(out)  TEML,
real(rp), dimension(ka), intent(out)  LHVL,
real(rp), dimension(ka), intent(out)  psat 
)

Definition at line 2851 of file scale_atmos_phy_bl_mynn.F90.

2851  use scale_const, only: &
2852  cpdry => const_cpdry, &
2853  rvap => const_rvap, &
2854  epsvap => const_epsvap, &
2855  epstvap => const_epstvap
2856  !$acc routine vector
2857  use scale_atmos_saturation, only: &
2858  atmos_saturation_psat => atmos_saturation_psat_liq
2859 ! ATMOS_SATURATION_pres2qsat => ATMOS_SATURATION_pres2qsat_liq
2860  use scale_atmos_hydrometeor, only: &
2861  hydrometeor_lhv => atmos_hydrometeor_lhv, &
2862  cp_vapor
2863  integer, intent(in), value :: KA, KS, KE
2864  real(RP), intent(in) :: PRES (KA)
2865  real(RP), intent(in) :: POTT (KA)
2866  real(RP), intent(in) :: POTL (KA)
2867  real(RP), intent(in) :: Qw (KA)
2868  real(RP), intent(in) :: EXNER(KA)
2869  real(RP), intent(in) :: tsq (KA)
2870  real(RP), intent(in) :: qsq (KA)
2871  real(RP), intent(in) :: cov (KA)
2872 
2873  real(RP), intent(out) :: betat (KA)
2874  real(RP), intent(out) :: betaq (KA)
2875  real(RP), intent(out) :: Qlp (KA)
2876  real(RP), intent(out) :: cldfrac(KA)
2877 
2878  real(RP), intent(out) :: TEML(KA)
2879  real(RP), intent(out) :: LHVL(KA)
2880  real(RP), intent(out) :: psat(KA)
2881 
2882  real(RP) :: Qsl
2883  real(RP) :: dQsl
2884  real(RP) :: CPtot
2885  real(RP) :: aa, bb, cc
2886  real(RP) :: sigma_s
2887  real(RP) :: Q1
2888  real(RP) :: Rt
2889 
2890  integer :: k
2891 
2892  do k = ks, ke
2893  teml(k) = potl(k) * exner(k)
2894  end do
2895 
2896  call atmos_saturation_psat( &
2897  ka, ks, ke, &
2898  teml(:), & ! (in)
2899  psat(:) ) ! (out)
2900 
2901  call hydrometeor_lhv( &
2902  ka, ks, ke, & ! (in)
2903  teml(:), & ! (in)
2904  lhvl(:) ) ! (out)
2905 
2906  do k = ks, ke
2907 
2908  qsl = epsvap * psat(k) / ( pres(k) - ( 1.0_rp - epsvap ) * psat(k) )
2909  dqsl = pres(k) * qsl**2 * lhvl(k) / ( epsvap * psat(k) * rvap * teml(k)**2 )
2910 
2911  cptot = ( 1.0_rp - qsl ) * cpdry + qsl * cp_vapor
2912  aa = 1.0_rp / ( 1.0_rp + lhvl(k)/cptot * dqsl )
2913  bb = exner(k) * dqsl
2914 
2915  sigma_s = min( max( &
2916  0.5_rp * aa * sqrt( max( qsq(k) - 2.0_rp * bb * cov(k) + bb**2 * tsq(k), 1.0e-20_rp ) ), &
2917  aa * qsl * 0.09_rp), aa * qsl )
2918  q1 = aa * ( qw(k) - qsl ) * 0.5_rp / sigma_s
2919  cldfrac(k) = min( max( 0.5_rp * ( 1.0_rp + erf(q1*rsqrt_2) ), 0.0_rp ), 1.0_rp )
2920  qlp(k) = min( max( 2.0_rp * sigma_s * ( cldfrac(k) * q1 + rsqrt_2pi &
2921 #if defined(NVIDIA) || defined(SX)
2922  * exp( -min( 0.5_rp*q1**2, 1.e+3_rp ) ) & ! apply exp limiter
2923 #else
2924  * exp(-0.5_rp*q1**2) &
2925 #endif
2926  ), 0.0_rp ), qw(k) * 0.5_rp )
2927  cc = ( 1.0_rp + epstvap * qw(k) - qlp(k) / epsvap ) / exner(k) * lhvl(k) / cptot &
2928  - pott(k) / epsvap
2929  rt = cldfrac(k) - qlp(k) / (2.0_rp*sigma_s*sqrt_2pi) &
2930 #if defined(NVIDIA) || defined(SX)
2931  * exp( -min( q1**2 * 0.5_rp, 1.e+3_rp ) ) ! apply exp limiter
2932 #else
2933  * exp(-q1**2 * 0.5_rp)
2934 #endif
2935  betat(k) = 1.0_rp + epstvap * qw(k) - qlp(k) / epsvap - rt * aa * bb * cc
2936  betaq(k) = epstvap * pott(k) + rt * aa * cc
2937 
2938  end do
2939 
2940  return

References scale_const::const_cpdry, scale_const::const_epstvap, scale_const::const_epsvap, scale_const::const_rvap, and scale_atmos_hydrometeor::cp_vapor.

Referenced by atmos_phy_bl_mynn_tendency().

Here is the caller graph for this function:

◆ get_phi()

subroutine scale_atmos_phy_bl_mynn::get_phi ( real(rp), intent(out)  zeta,
real(rp), intent(out)  phi_m,
real(rp), intent(out)  phi_h,
real(rp), intent(in)  z1,
real(rp), intent(in)  RLmo,
integer, intent(in)  I_B_TYPE 
)

Definition at line 2946 of file scale_atmos_phy_bl_mynn.F90.

2946  !$acc routine seq
2947  real(RP), intent(out) :: zeta
2948  real(RP), intent(out) :: phi_m
2949  real(RP), intent(out) :: phi_h
2950  real(RP), intent(in) :: z1
2951  real(RP), intent(in) :: RLmo
2952  integer, intent(in) :: I_B_TYPE
2953 
2954  real(RP) :: tmp
2955 
2956  zeta = min( max( z1 * rlmo, -atmos_phy_bl_mynn_zeta_lim ), atmos_phy_bl_mynn_zeta_lim )
2957 
2958  select case ( i_b_type )
2959  case ( i_b71 )
2960  ! Businger et al. (1971)
2961  if ( zeta >= 0 ) then
2962  phi_m = 4.7_rp * zeta + 1.0_rp
2963  phi_h = 4.7_rp * zeta + 0.74_rp
2964  else
2965  phi_m = 1.0_rp / sqrt(sqrt( 1.0_rp - 15.0_rp * zeta ))
2966  phi_h = 0.47_rp / sqrt( 1.0_rp - 9.0_rp * zeta )
2967  end if
2968  case ( i_b91, i_b91w01 )
2969  ! Beljaars and Holtslag (1991)
2970  if ( zeta >= 0 ) then
2971  tmp = - 2.0_rp / 3.0_rp * ( 0.35_rp * zeta - 6.0_rp ) * exp(-0.35_rp*zeta) * zeta
2972  phi_m = tmp + zeta + 1.0_rp
2973  phi_h = tmp + zeta * sqrt( 1.0_rp + 2.0_rp * zeta / 3.0_rp ) + 1.0_rp
2974  else
2975  if ( i_b_type == i_b91w01 ) then
2976  ! Wilson (2001)
2977  !tmp = (-zeta)**(2.0_RP/3.0_RP)
2978  tmp = abs(zeta)**(2.0_rp/3.0_rp)
2979  phi_m = 1.0_rp / sqrt( 1.0_rp + 3.6_rp * tmp )
2980  phi_h = 0.95_rp / sqrt( 1.0_rp + 7.9_rp * tmp )
2981  else
2982  ! tmp = sqrt( 1.0_RP - 16.0_RP * zeta )
2983  tmp = sqrt( 1.0_rp + 16.0_rp * abs(zeta) )
2984  phi_m = 1.0_rp / sqrt(tmp)
2985  phi_h = 1.0_rp / tmp
2986  end if
2987  end if
2988  end select
2989 
2990  return

References scale_const::const_cpdry, scale_const::const_cvdry, scale_const::const_epstvap, scale_const::const_grav, scale_atmos_hydrometeor::cp_vapor, scale_atmos_hydrometeor::cv_vapor, and scale_atmos_hydrometeor::lhv.

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 57 of file scale_atmos_phy_bl_mynn.F90.

57  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 58 of file scale_atmos_phy_bl_mynn.F90.

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

Referenced by mod_atmos_phy_bl_driver::atmos_phy_bl_driver_tracer_setup(), atmos_phy_bl_mynn_finalize(), 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 59 of file scale_atmos_phy_bl_mynn.F90.

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

Referenced by mod_atmos_phy_bl_driver::atmos_phy_bl_driver_tracer_setup(), atmos_phy_bl_mynn_finalize(), 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 60 of file scale_atmos_phy_bl_mynn.F90.

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

Referenced by mod_atmos_phy_bl_driver::atmos_phy_bl_driver_tracer_setup(), atmos_phy_bl_mynn_finalize(), 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:49
scale_prc::prc_abort
subroutine, public prc_abort
Abort Process.
Definition: scale_prc.F90:349
scale_const::const_epstvap
real(rp), public const_epstvap
1 / epsilon - 1
Definition: scale_const.F90:76
scale_const::const_epsvap
real(rp), public const_epsvap
Rdry / Rvap.
Definition: scale_const.F90:75
scale_const::const_rvap
real(rp), parameter, public const_rvap
specific gas constant (water vapor) [J/kg/K]
Definition: scale_const.F90:68
scale_const::const_eps
real(rp), public const_eps
small number
Definition: scale_const.F90:35
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_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, W, POTT, PROG, PRES, EXNER, N2, QDRY, QV, Qw, POTL, POTV, SFC_DENS, SFLX_MU, SFLX_MV, SFLX_SH, SFLX_QV, us, ts, qs, RLmo, frac_land, CZ, FZ, F2H, dt_DP, BULKFLUX_type, RHOU_t, RHOV_t, RHOT_t, RHOQV_t, RPROG_t, Nu, Kh, Qlp, cldfrac, Zi, SFLX_BUOY)
ATMOS_PHY_BL_MYNN_tendency calculate tendency by the vertical eddy viscosity.
Definition: scale_atmos_phy_bl_mynn.F90:720
scale_prc
module PROCESS
Definition: scale_prc.F90:11
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:60
scale_const::const_huge
real(rp), public const_huge
huge number
Definition: scale_const.F90:37
mod_atmos_phy_tb_vars::i_tke
integer, public i_tke
Definition: mod_atmos_phy_tb_vars.F90:62
scale_const::const_karman
real(rp), parameter, public const_karman
von Karman constant
Definition: scale_const.F90:54
scale_matrix
module MATRIX
Definition: scale_matrix.F90:17
scale_file_history::file_history_reg
subroutine, public file_history_reg(name, desc, unit, itemid, standard_name, ndims, dim_type, cell_measures, fill_halo)
Register/Append variable to history file.
Definition: scale_file_history.F90:685
scale_atmos_saturation
module atmosphere / saturation
Definition: scale_atmos_saturation.F90:12
scale_atmos_hydrometeor::cp_vapor
real(rp), public cp_vapor
CP for vapor [J/kg/K].
Definition: scale_atmos_hydrometeor.F90:150
Generated by   doxygen 1.8.17