SCALE-RM
scale_atmos_phy_bl_mynn.F90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
21 !-------------------------------------------------------------------------------
22 #include "scalelib.h"
23 module scale_atmos_phy_bl_mynn
24  !-----------------------------------------------------------------------------
25  !
26  !++ used modules
27  !
28  use scale_precision
29  use scale_io
30  use scale_prof
31 
32 #if defined DEBUG || defined QUICKDEBUG
33  use scale_debug, only: &
34  check
35 #endif
36  use scale_const, only: &
37  pi => const_pi, &
38  undef => const_undef, &
39  iundef => const_undef2
40  !-----------------------------------------------------------------------------
41  implicit none
42  private
43  !-----------------------------------------------------------------------------
44  !
45  !++ Public procedure
46  !
47  public :: atmos_phy_bl_mynn_tracer_setup
48  public :: atmos_phy_bl_mynn_setup
49  public :: atmos_phy_bl_mynn_finalize
50  public :: atmos_phy_bl_mynn_mkinit
51  public :: atmos_phy_bl_mynn_tendency
52 
53  !-----------------------------------------------------------------------------
54  !
55  !++ Public parameters & variables
56  !
57  integer, public :: atmos_phy_bl_mynn_ntracer
58  character(len=H_SHORT), public, allocatable :: atmos_phy_bl_mynn_name(:)
59  character(len=H_LONG), public, allocatable :: atmos_phy_bl_mynn_desc(:)
60  character(len=H_SHORT), public, allocatable :: atmos_phy_bl_mynn_units(:)
61 
62  !-----------------------------------------------------------------------------
63  !
64  !++ Private procedure
65  !
66  !-----------------------------------------------------------------------------
67  !
68  !++ Private parameters & variables
69  !
70  integer, private, parameter :: i_tke = 1
71  integer, private, parameter :: i_tsq = 2
72  integer, private, parameter :: i_qsq = 3
73  integer, private, parameter :: i_cov = 4
74 
75  integer, private, parameter :: i_b71 = 1
76  integer, private, parameter :: i_b91 = 2
77  integer, private, parameter :: i_b91w01 = 3
78 
79  real(rp), private, parameter :: oneoverthree = 1.0_rp / 3.0_rp
80  real(rp), private, parameter :: lt_min = 1.e-6_rp
81  real(rp), private, parameter :: flx_lim_fact = 0.5_rp
82 
83  real(rp), private :: a1
84  real(rp), private :: a2
85  real(rp), private, parameter :: b1 = 24.0_rp
86  real(rp), private, parameter :: b2 = 15.0_rp
87  real(rp), private :: c1
88  real(rp), private :: c2 = -1.0_rp ! N01: 0.65, N04: 0.70, N09: 0.75, O19: 0.729
89  real(rp), private :: c3 = -1.0_rp ! N01: 0.294, N04: 0.323, N09: 0.352, O19: 0.34
90  real(rp), private, parameter :: c5 = 0.2_rp
91  real(rp), private, parameter :: g1 = 0.235_rp
92  real(rp), private :: g2
93  real(rp), private :: f2
94  real(rp), private :: rf2
95  real(rp), private :: rfc
96  real(rp), private, parameter :: prn = 0.74_rp
97  !$acc declare create(A1, A2, C1, C2, C3, G2, F2, Rf2, Rfc)
98 
99  real(rp), private, parameter :: sqrt_2pi = sqrt( 2.0_rp * pi )
100  real(rp), private, parameter :: rsqrt_2pi = 1.0_rp / sqrt_2pi
101  real(rp), private, parameter :: rsqrt_2 = 1.0_rp / sqrt( 2.0_rp )
102 
103  ! for O2019
104  integer, parameter :: max_plume = 10
105  integer :: nplume
106  real(rp) :: dplume(max_plume)
107  real(rp) :: pw(max_plume)
108  logical :: atmos_phy_bl_mynn_mf
109  !$acc declare create(nplume,dplume,pw,ATMOS_PHY_BL_MYNN_MF)
110 
111  ! history
112  integer, private :: hist_ri
113  integer, private :: hist_pr
114  integer, private :: hist_tke_pr
115  integer, private :: hist_tke_di
116  integer, private :: hist_dudz2
117  integer, private :: hist_lmix
118  integer, private :: hist_flxu
119  integer, private :: hist_flxv
120  integer, private :: hist_flxt
121  integer, private :: hist_flxq
122  integer, private :: hist_flxu2
123  integer, private :: hist_flxv2
124  integer, private :: hist_flxt2
125  integer, private :: hist_flxq2
126 
127  ! namelist
128  logical, private :: atmos_phy_bl_mynn_k2010 = .false.
129  logical, private :: atmos_phy_bl_mynn_o2019 = .false.
130  real(rp), private :: atmos_phy_bl_mynn_pbl_max = 3000.0_rp
131  real(rp), private :: atmos_phy_bl_mynn_tke_min = 1.e-20_rp
132  real(rp), private :: atmos_phy_bl_mynn_n2_max = 1.e1_rp
133  real(rp), private :: atmos_phy_bl_mynn_nu_max = 1.e4_rp
134  real(rp), private :: atmos_phy_bl_mynn_kh_max = 1.e4_rp
135  real(rp), private :: atmos_phy_bl_mynn_lt_max = 2000.0_rp
136  logical, private :: atmos_phy_bl_mynn_use_zi = .true.
137  real(rp), private :: atmos_phy_bl_mynn_zeta_lim = 10.0_rp
138  real(rp), private :: atmos_phy_bl_mynn_sq_fact = 3.0_rp
139  real(rp), private :: atmos_phy_bl_mynn_cns = -1.0_rp ! N01: 3.5, O19: 10.0
140  real(rp), private :: atmos_phy_bl_mynn_alpha2 = -1.0_rp ! N01: 1.0, O19: 0.3
141  real(rp), private :: atmos_phy_bl_mynn_alpha4 = -1.0_rp ! N01: 100.0, O19: 10.0
142  real(rp), private :: atmos_phy_bl_mynn_dump_coef = 0.5_rp
143  logical, private :: atmos_phy_bl_mynn_dz_sim = .true.
144  logical, private :: atmos_phy_bl_mynn_similarity = .true.
145  !$acc declare create(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_use_Zi,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)
146 
147  character(len=H_SHORT), private :: atmos_phy_bl_mynn_level = "2.5" ! "2.5" or "3"
148 
149  namelist / param_atmos_phy_bl_mynn / &
150  atmos_phy_bl_mynn_k2010, &
151  atmos_phy_bl_mynn_o2019, &
152  atmos_phy_bl_mynn_pbl_max, &
153  atmos_phy_bl_mynn_n2_max, &
154  atmos_phy_bl_mynn_nu_max, &
155  atmos_phy_bl_mynn_kh_max, &
156  atmos_phy_bl_mynn_lt_max, &
157  atmos_phy_bl_mynn_use_zi, &
158  atmos_phy_bl_mynn_zeta_lim, &
159  atmos_phy_bl_mynn_level, &
160  atmos_phy_bl_mynn_sq_fact, &
161  atmos_phy_bl_mynn_cns, &
162  atmos_phy_bl_mynn_alpha2, &
163  atmos_phy_bl_mynn_alpha4, &
164  atmos_phy_bl_mynn_dump_coef,&
165  atmos_phy_bl_mynn_dz_sim, &
166  atmos_phy_bl_mynn_similarity
167 
168  !-----------------------------------------------------------------------------
169 contains
170  !-----------------------------------------------------------------------------
174  subroutine atmos_phy_bl_mynn_tracer_setup
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
221  end subroutine atmos_phy_bl_mynn_tracer_setup
222 
223  !-----------------------------------------------------------------------------
227  subroutine atmos_phy_bl_mynn_setup( &
228  BULKFLUX_type, &
229  dx, &
230  TKE_MIN, PBL_MAX )
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
343  end subroutine atmos_phy_bl_mynn_setup
344 
345  !-----------------------------------------------------------------------------
346  !! Finalize
347 
348  subroutine atmos_phy_bl_mynn_finalize
349  implicit none
350 
351  deallocate( atmos_phy_bl_mynn_name, atmos_phy_bl_mynn_desc, atmos_phy_bl_mynn_units )
352 
353  return
354  end subroutine atmos_phy_bl_mynn_finalize
355 
356  !-----------------------------------------------------------------------------
360  subroutine atmos_phy_bl_mynn_mkinit( &
361  KA, KS, KE, IA, IS, IE, JA, JS, JE, &
362  PROG, &
363  DENS, U, V, W, POTT, &
364  PRES, EXNER, N2, &
365  QDRY, QV, Qw, POTL, POTV, SFC_DENS, &
366  SFLX_MU, SFLX_MV, SFLX_SH, SFLX_QV, &
367  us, ts, qs, RLmo, &
368  frac_land, &
369  CZ, FZ, F2H, &
370  BULKFLUX_type )
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
700  end subroutine atmos_phy_bl_mynn_mkinit
701 
702  !-----------------------------------------------------------------------------
706  subroutine atmos_phy_bl_mynn_tendency( &
707  KA, KS, KE, IA, IS, IE, JA, JS, JE, &
708  DENS, U, V, W, POTT, PROG, &
709  PRES, EXNER, N2, &
710  QDRY, QV, Qw, POTL, POTV, SFC_DENS, &
711  SFLX_MU, SFLX_MV, SFLX_SH, SFLX_QV, &
712  us, ts, qs, RLmo, &
713  frac_land, &
714  CZ, FZ, F2H, dt_DP, &
715  BULKFLUX_type, &
716  RHOU_t, RHOV_t, RHOT_t, RHOQV_t, &
717  RPROG_t, &
718  Nu, Kh, Qlp, cldfrac, &
719  Zi, SFLX_BUOY )
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
2740  end subroutine get_gamma_implicit
2741 
2742 !OCL SERIAL
2743  subroutine calc_vertical_differece( &
2744  KA, KS, KE, &
2745  i, j, &
2746  dudz2, dtldz, dqwdz, &
2747  U, V, POTL, Qw, QDRY, &
2748  CDZ, FDZ, F2H, &
2749  us, ts, qs, &
2750  phi_m, phi_h, z1, &
2751  Uh, Vh, qw2, qh )
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
2842  end subroutine calc_vertical_differece
2843 
2844 !OCL SERIAL
2845  subroutine partial_condensation( &
2846  KA, KS, KE, &
2847  betat, betaq, Qlp, cldfrac, &
2848  PRES, POTT, POTL, Qw, EXNER, &
2849  tsq, qsq, cov, &
2850  TEML, LHVL, psat )
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
2941  end subroutine partial_condensation
2942 
2943  subroutine get_phi( &
2944  zeta, phi_m, phi_h, &
2945  z1, RLmo, I_B_TYPE )
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
2991  end subroutine get_phi
2992 
2993 !OCL SERIAL
2994  subroutine calc_zi_o2019( &
2995  KA, KS, KE, &
2996  Zi, &
2997  PTLV, &
2998  tke, &
2999  z, frac_land, &
3000  initialize )
3001  !$acc routine seq
3002  integer, intent(in) :: KA, KS, KE
3003  real(RP), intent(out) :: Zi
3004  real(RP), intent(in) :: PTLV(KA)
3005  real(RP), intent(in) :: tke(KA)
3006  real(RP), intent(in) :: z(KA)
3007  real(RP), intent(in) :: frac_land
3008  logical, intent(in) :: initialize
3009 
3010  real(RP) :: zit, zie
3011  real(RP) :: tmin, dtv
3012  real(RP) :: tke_min
3013  real(RP) :: we
3014  integer :: k, k2
3015 
3016  tmin = 999e10_rp
3017  k2 = ke
3018  do k = ks, ke
3019  tmin = min( tmin, ptlv(k) )
3020  if ( z(k) > 200.0_rp ) then
3021  k2 = k
3022  exit
3023  end if
3024  end do
3025  if ( frac_land >= 0.5_rp ) then ! over land
3026  dtv = 1.25_rp
3027  else ! over water
3028  dtv = 0.75_rp
3029  end if
3030  zit = z(ke)
3031  do k = k2+1, ke
3032  if ( tmin + dtv <= ptlv(k) ) then
3033  zit = z(k)
3034  exit
3035  end if
3036  end do
3037  zie = 1000.0_rp
3038  if ( .not. initialize ) then
3039  tke_min = max( tke(ks), 0.02_rp )
3040  do k = ks+1, ke
3041  if ( tke(k) <= tke_min * 0.05_rp ) then
3042  zie = z(k)
3043  exit
3044  end if
3045  end do
3046  end if
3047  we = 0.5_rp * tanh( ( zit - 200.0_rp ) / 400.0_rp ) + 0.5_rp
3048  zi = zit * ( 1.0_rp - we ) + zie * we
3049 
3050  return
3051  end subroutine calc_zi_o2019
3052 
3053 !OCL SERIAL
3054  subroutine calc_mflux( &
3055  KA, KS, KE, &
3056  DENS, &
3057  POTV, POTL, Qw, &
3058  U, V, W, &
3059  tke, &
3060  cldfrac, &
3061  EXNER, &
3062  SHFLX, SFLX_BUOY, SFLX_PT, SFLX_QV, &
3063  SFC_DENS, &
3064  Zi, &
3065  Z, CDZ, F2H, &
3066  mflux, tflux, qflux, uflux, vflux, eflux, &
3067  tlh, tvh, qh, uh, vh, wh, eh, dh )
3068  !$acc routine vector
3069  use scale_const, only: &
3070  grav => const_grav, &
3071  cpdry => const_cpdry, &
3072  cvdry => const_cvdry, &
3073  epstvap => const_epstvap
3074  use scale_atmos_hydrometeor, only: &
3075  cp_vapor, &
3076  cv_vapor, &
3077  lhv
3078  use scale_atmos_saturation, only: &
3079  atmos_saturation_moist_conversion_dens_liq
3080  integer, intent(in) :: KA, KS, KE
3081 
3082  real(RP), intent(in) :: DENS(KA)
3083  real(RP), intent(in) :: POTV(KA)
3084  real(RP), intent(in) :: POTL(KA)
3085  real(RP), intent(in) :: Qw(KA)
3086  real(RP), intent(in) :: U(KA)
3087  real(RP), intent(in) :: V(KA)
3088  real(RP), intent(in) :: W(KA)
3089  real(RP), intent(in) :: tke(KA)
3090  real(RP), intent(in) :: cldfrac(KA)
3091  real(RP), intent(in) :: EXNER(KA)
3092  real(RP), intent(in) :: SHFLX
3093  real(RP), intent(in) :: SFLX_BUOY
3094  real(RP), intent(in) :: SFLX_PT
3095  real(RP), intent(in) :: SFLX_QV
3096  real(RP), intent(in) :: SFC_DENS
3097  real(RP), intent(in) :: Zi
3098  real(RP), intent(in) :: Z(KA)
3099  real(RP), intent(in) :: CDZ(KA)
3100  real(RP), intent(in) :: F2H(KA,2)
3101 
3102  real(RP), intent(out) :: mflux(0:KA)
3103  real(RP), intent(out) :: tflux(0:KA)
3104  real(RP), intent(out) :: qflux(0:KA)
3105  real(RP), intent(out) :: uflux(0:KA)
3106  real(RP), intent(out) :: vflux(0:KA)
3107  real(RP), intent(out) :: eflux(0:KA)
3108 
3109  ! work
3110  real(RP), intent(out) :: tlh(KA)
3111  real(RP), intent(out) :: tvh(KA)
3112  real(RP), intent(out) :: qh(KA)
3113  real(RP), intent(out) :: uh(KA)
3114  real(RP), intent(out) :: vh(KA)
3115  real(RP), intent(out) :: wh(KA)
3116  real(RP), intent(out) :: eh(KA)
3117  real(RP), intent(out) :: dh(KA)
3118 
3119  real(RP), parameter :: amax = 0.1_rp
3120  real(RP), parameter :: Ce = 0.35_rp
3121  real(RP), parameter :: x = -1.9_rp
3122  real(RP), parameter :: Cwt = 0.58_rp
3123  real(RP), parameter :: Cs = 1.34_rp
3124  real(RP), parameter :: zs = 50.0_rp
3125  real(RP), parameter :: a = 2.0_rp
3126 
3127  real(RP) :: au ! total area
3128  real(RP) :: ap(nplume) ! area
3129  real(RP) :: wp, tp, qp, up, vp, ep
3130  real(RP) :: zmax
3131  real(RP) :: er ! entrainment rate
3132  real(RP) :: sw, sq, st ! sigma_w,q,t
3133  real(RP) :: ws, ts, qs ! scales
3134  real(RP) :: buoy ! buoyancy
3135  real(RP) :: mf
3136  real(RP) :: b
3137  real(RP) :: dd
3138  real(RP) :: temp, qc
3139  real(RP) :: tps, qps
3140  real(RP) :: Emoist
3141  real(RP) :: cp, cv
3142  real(RP) :: tmp, fact
3143  logical :: converged
3144  integer :: np
3145  integer :: k, n
3146 
3147  do k = ks-1, ke
3148  mflux(k) = 0.0_rp
3149  tflux(k) = 0.0_rp
3150  qflux(k) = 0.0_rp
3151  uflux(k) = 0.0_rp
3152  vflux(k) = 0.0_rp
3153  eflux(k) = 0.0_rp
3154  end do
3155 
3156  ! must be a positive surface buoyancy flux
3157  if ( sflx_buoy <= 0.0_rp ) then
3158  return
3159  end if
3160 
3161  ! must be superadiabatic in the lowest 50 m
3162  !$acc loop seq
3163  do k = ks, ke-1
3164  if ( z(k) > 50.0_rp ) exit
3165  if ( potv(k+1) >= potv(k) ) then
3166  return
3167  end if
3168  end do
3169 
3170  zmax = zi
3171  !$acc loop seq
3172  do k = ks, ke
3173  if ( cldfrac(k) > 0.5_rp .and. z(k) <= 2000.0_rp ) then
3174  zmax = min( zmax, z(k) * 0.5_rp ) ! zc
3175  exit
3176  end if
3177  end do
3178  np = nplume
3179  !$acc loop seq
3180  do n = 1, nplume
3181  if ( dplume(n) > zmax ) then
3182  np = n - 1
3183  exit
3184  end if
3185  end do
3186  if ( np == 0 ) return
3187 
3188  ! updraft area
3189  au = amax * ( tanh( ( shflx - 20.0_rp ) / 90.0_rp ) * 0.5_rp + 0.5_rp )
3190  fact = 0.0_rp
3191  !$acc loop reduction(+:fact)
3192  do n = 1, np
3193  if ( n == 1 ) then
3194  dd = dplume(2) - dplume(1)
3195  else if ( n == np ) then
3196  dd = dplume(np) - dplume(np-1)
3197  else
3198  dd = ( dplume(n+1) - dplume(n-1) ) * 0.5_rp
3199  end if
3200  ap(n) = dplume(n)**(x + 2.0_rp) * dd
3201  fact = fact + ap(n)
3202  end do
3203  fact = au / fact
3204  do n = 1, np
3205  ap(n) = ap(n) * fact
3206  end do
3207 
3208  ! @half levels
3209  do k = ks, ke-1
3210  tlh(k) = f2h(k,1) * potl(k) + f2h(k,2) * potl(k+1)
3211  tvh(k) = f2h(k,1) * potv(k) + f2h(k,2) * potv(k+1)
3212  qh(k) = f2h(k,1) * qw(k) + f2h(k,2) * qw(k+1)
3213  uh(k) = f2h(k,1) * u(k) + f2h(k,2) * u(k+1)
3214  vh(k) = f2h(k,1) * v(k) + f2h(k,2) * v(k+1)
3215  wh(k) = f2h(k,1) * w(k) + f2h(k,2) * w(k+1)
3216  eh(k) = f2h(k,1) * tke(k) + f2h(k,2) * tke(k+1)
3217  dh(k) = f2h(k,1) * dens(k) + f2h(k,2) * dens(k+1)
3218  end do
3219 
3220  ws = max( ( zi * sflx_buoy )**oneoverthree, 1.0e-10_rp )
3221  ts = sflx_pt / ( sfc_dens * ws )
3222  qs = max( sflx_qv, 0.0_rp ) / ( sfc_dens * ws )
3223  tmp = ( zs / zi )**oneoverthree
3224  sw = cs * ws * tmp * ( 1.0_rp - 0.8_rp*zs/zi )
3225  sq = cs * qs * tmp
3226  st = cs * ts * tmp
3227  !$acc loop seq
3228  do n = 1, np
3229  ! initial properties of the plume
3230  wp = min( pw(n) * sw, 0.5_rp )
3231  tp = tlh(ks) + wp * cwt * st / sw
3232  qp = qh(ks) + wp * cwt * sq / sw
3233  up = uh(ks)
3234  vp = vh(ks)
3235  ep = eh(ks)
3236  !$acc loop seq
3237  do k = ks, ke-1
3238  mf = ap(n) * ( wp - wh(k) ) * dh(k)
3239  mflux(k) = mflux(k) + mf
3240  tflux(k) = tflux(k) + mf * ( tp - tlh(k) )
3241  qflux(k) = qflux(k) + mf * ( qp - qh(k) )
3242  uflux(k) = uflux(k) + mf * ( up - uh(k) )
3243  vflux(k) = vflux(k) + mf * ( vp - vh(k) )
3244  eflux(k) = eflux(k) + mf * ( ep - eh(k) )
3245 
3246  er = ce / ( wp * dplume(n) )
3247 
3248  ! d/dz phi = - er ( phi - phi_env )
3249  fact = exp( - er * cdz(k+1) )
3250  tp = potl(k+1) + ( tp - potl(k+1) ) * fact
3251  qp = qw(k+1) + ( qp - qw(k+1) ) * fact
3252  up = u(k+1) + ( up - u(k+1) ) * fact
3253  vp = v(k+1) + ( vp - v(k+1) ) * fact
3254  ep = tke(k+1) + ( ep - tke(k+1) ) * fact
3255 
3256  qps = qp
3257  qc = 0.0_rp
3258  temp = tp * exner(k+1)
3259  cp = cpdry * ( 1.0_rp - qp ) + cp_vapor * qp
3260  cv = cvdry * ( 1.0_rp - qp ) + cv_vapor * qp
3261  emoist = temp * cv + qp * lhv
3262  call atmos_saturation_moist_conversion_dens_liq( dens(k+1), emoist, & ! (in)
3263  temp, qps, qc, & ! (inout)
3264  cp, cv, & ! (inout)
3265  converged ) ! (out)
3266  if ( converged ) then
3267  tps = temp / exner(k+1)
3268  else
3269  log_warn("calc_mflux",*) "moist_conversion did not converged"
3270  ! tentative
3271  qps = qp
3272  tps = tp
3273  end if
3274  buoy = grav * ( tps * ( 1.0_rp + epstvap * qps ) / potv(k+1) - 1.0_rp )
3275  if ( buoy > 0.0_rp ) then
3276  b = 0.15_rp
3277  else
3278  b = 0.2_rp
3279  end if
3280  ! d/dz w = - er a w + b buoy / w
3281  wp = sqrt( max( 0.0_rp, ( er * a * wp**2 - b * buoy ) * exp( - 2.0_rp * er * a * cdz(k+1) ) + b * buoy ) / ( er * a ) )
3282  if ( wp <= 0.0_rp ) then
3283  exit
3284  end if
3285  end do
3286  end do
3287 
3288 
3289  return
3290  end subroutine calc_mflux
3291 
3292 end module scale_atmos_phy_bl_mynn
scale_const::const_grav
real(rp), public const_grav
standard acceleration of gravity [m/s2]
Definition: scale_const.F90:49
scale_atmos_phy_bl_mynn::calc_vertical_differece
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)
Definition: scale_atmos_phy_bl_mynn.F90:2752
scale_atmos_grid_cartesc_index::ke
integer, public ke
end point of inner domain: z, local
Definition: scale_atmos_grid_cartesC_index.F90:52
scale_atmos_phy_bl_mynn::atmos_phy_bl_mynn_desc
character(len=h_long), dimension(:), allocatable, public atmos_phy_bl_mynn_desc
Definition: scale_atmos_phy_bl_mynn.F90:59
scale_prc::prc_abort
subroutine, public prc_abort
Abort Process.
Definition: scale_prc.F90:350
scale_atmos_phy_bl_mynn::get_phi
subroutine get_phi(zeta, phi_m, phi_h, z1, RLmo, I_B_TYPE)
Definition: scale_atmos_phy_bl_mynn.F90:2946
scale_const::const_epstvap
real(rp), public const_epstvap
1 / epsilon - 1
Definition: scale_const.F90:76
scale_const::const_undef2
integer, parameter, public const_undef2
undefined value (INT2)
Definition: scale_const.F90:40
scale_const::const_epsvap
real(rp), public const_epsvap
Rdry / Rvap.
Definition: scale_const.F90:75
scale_atmos_phy_bl_mynn::atmos_phy_bl_mynn_units
character(len=h_short), dimension(:), allocatable, public atmos_phy_bl_mynn_units
Definition: scale_atmos_phy_bl_mynn.F90:60
scale_precision
module PRECISION
Definition: scale_precision.F90:14
scale_atmos_grid_cartesc_index::ka
integer, public ka
Definition: scale_atmos_grid_cartesC_index.F90:47
scale_atmos_phy_bl_mynn
module atmosphere / physics / pbl / mynn
Definition: scale_atmos_phy_bl_mynn.F90:23
scale_const::const_rvap
real(rp), parameter, public const_rvap
specific gas constant (water vapor) [J/kg/K]
Definition: scale_const.F90:68
scale_atmos_phy_bl_mynn::atmos_phy_bl_mynn_mkinit
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.
Definition: scale_atmos_phy_bl_mynn.F90:371
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_precision::rp
integer, parameter, public rp
Definition: scale_precision.F90:41
scale_atmos_hydrometeor::cv_vapor
real(rp), public cv_vapor
CV for vapor [J/kg/K].
Definition: scale_atmos_hydrometeor.F90:149
scale_atmos_grid_cartesc_index::ie
integer, public ie
end point of inner domain: x, local
Definition: scale_atmos_grid_cartesC_index.F90:54
scale_io
module STDIO
Definition: scale_io.F90:10
scale_tracer::k
real(rp), public k
Definition: scale_tracer.F90:45
scale_const
module CONSTANT
Definition: scale_const.F90:11
scale_atmos_grid_cartesc_index::ia
integer, public ia
Definition: scale_atmos_grid_cartesC_index.F90:48
scale_debug::check
subroutine, public check(current_line, v)
Undefined value checker.
Definition: scale_debug.F90:59
scale_atmos_phy_bl_mynn::atmos_phy_bl_mynn_finalize
subroutine, public atmos_phy_bl_mynn_finalize
Definition: scale_atmos_phy_bl_mynn.F90:349
scale_atmos_hydrometeor::lhv
real(rp), public lhv
latent heat of vaporization for use [J/kg]
Definition: scale_atmos_hydrometeor.F90:144
scale_atmos_grid_cartesc_index::kmax
integer, public kmax
Definition: scale_atmos_grid_cartesC_index.F90:36
scale_const::const_cvdry
real(rp), public const_cvdry
specific heat (dry air,constant volume) [J/kg/K]
Definition: scale_const.F90:61
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
scale_prof
module profiler
Definition: scale_prof.F90:11
scale_atmos_grid_cartesc_index::is
integer, public is
start point of inner domain: x, local
Definition: scale_atmos_grid_cartesC_index.F90:53
scale_precision::dp
integer, parameter, public dp
Definition: scale_precision.F90:32
scale_atmos_grid_cartesc_index::ja
integer, public ja
Definition: scale_atmos_grid_cartesC_index.F90:49
scale_const::const_pi
real(rp), parameter, public const_pi
pi
Definition: scale_const.F90:32
scale_atmos_grid_cartesc_index::ks
integer, public ks
start point of inner domain: z, local
Definition: scale_atmos_grid_cartesC_index.F90:51
scale_debug
module DEBUG
Definition: scale_debug.F90:11
scale_atmos_grid_cartesc_index::js
integer, public js
start point of inner domain: y, local
Definition: scale_atmos_grid_cartesC_index.F90:55
scale_const::const_karman
real(rp), parameter, public const_karman
von Karman constant
Definition: scale_const.F90:54
scale_matrix
module MATRIX
Definition: scale_matrix.F90:17
scale_atmos_phy_bl_mynn::atmos_phy_bl_mynn_ntracer
integer, public atmos_phy_bl_mynn_ntracer
Definition: scale_atmos_phy_bl_mynn.F90:57
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_phy_bl_mynn::partial_condensation
subroutine partial_condensation(KA, KS, KE, betat, betaq, Qlp, cldfrac, PRES, POTT, POTL, Qw, EXNER, tsq, qsq, cov, TEML, LHVL, psat)
Definition: scale_atmos_phy_bl_mynn.F90:2851
scale_const::const_undef
real(rp), public const_undef
Definition: scale_const.F90:43
scale_atmos_saturation
module atmosphere / saturation
Definition: scale_atmos_saturation.F90:12
scale_atmos_phy_bl_mynn::atmos_phy_bl_mynn_name
character(len=h_short), dimension(:), allocatable, public atmos_phy_bl_mynn_name
Definition: scale_atmos_phy_bl_mynn.F90:58
scale_io::io_fid_conf
integer, public io_fid_conf
Config file ID.
Definition: scale_io.F90:57
scale_atmos_grid_cartesc_index::je
integer, public je
end point of inner domain: y, local
Definition: scale_atmos_grid_cartesC_index.F90:56
scale_atmos_hydrometeor::cp_vapor
real(rp), public cp_vapor
CP for vapor [J/kg/K].
Definition: scale_atmos_hydrometeor.F90:150
scale_atmos_phy_bl_mynn::atmos_phy_bl_mynn_setup
subroutine, public atmos_phy_bl_mynn_setup(BULKFLUX_type, dx, TKE_MIN, PBL_MAX)
ATMOS_PHY_BL_MYNN_setup Setup.
Definition: scale_atmos_phy_bl_mynn.F90:231
scale_atmos_phy_bl_mynn::atmos_phy_bl_mynn_tracer_setup
subroutine, public atmos_phy_bl_mynn_tracer_setup
ATMOS_PHY_BL_MYNN_tracer_setup Tracer Setup.
Definition: scale_atmos_phy_bl_mynn.F90:175
Generated by   doxygen 1.8.17