SCALE-RM
scale_atmos_phy_bl_mynn.F90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
16 !-------------------------------------------------------------------------------
17 #include "scalelib.h"
18 module scale_atmos_phy_bl_mynn
19  !-----------------------------------------------------------------------------
20  !
21  !++ used modules
22  !
23  use scale_precision
24  use scale_io
25  use scale_prof
26 
27 #if defined DEBUG || defined QUICKDEBUG
28  use scale_debug, only: &
29  check
30 #endif
31  use scale_const, only: &
32  undef => const_undef, &
33  iundef => const_undef2
34  !-----------------------------------------------------------------------------
35  implicit none
36  private
37  !-----------------------------------------------------------------------------
38  !
39  !++ Public procedure
40  !
41  public :: atmos_phy_bl_mynn_tracer_setup
42  public :: atmos_phy_bl_mynn_setup
43  public :: atmos_phy_bl_mynn_tendency
44  public :: atmos_phy_bl_mynn_tendency_tracer
45 
46  !-----------------------------------------------------------------------------
47  !
48  !++ Public parameters & variables
49  !
50  integer, public :: atmos_phy_bl_mynn_ntracer
51  character(len=H_SHORT), public, allocatable :: atmos_phy_bl_mynn_name(:)
52  character(len=H_LONG), public, allocatable :: atmos_phy_bl_mynn_desc(:)
53  character(len=H_SHORT), public, allocatable :: atmos_phy_bl_mynn_units(:)
54 
55  !-----------------------------------------------------------------------------
56  !
57  !++ Private procedure
58  !
59  !-----------------------------------------------------------------------------
60  !
61  !++ Private parameters & variables
62  !
63  integer, private, parameter :: i_tke = 1
64  integer, private, parameter :: i_tsq = 2
65  integer, private, parameter :: i_qsq = 3
66  integer, private, parameter :: i_cov = 4
67 
68  real(rp), private, parameter :: oneoverthree = 1.0_rp / 3.0_rp
69  real(rp), private, parameter :: lt_min = 1.e-6_rp
70  real(rp), private, parameter :: flx_lim_fact = 0.5_rp
71 
72  real(rp), private :: a1
73  real(rp), private :: a2
74  real(rp), private, parameter :: b1 = 24.0_rp
75  real(rp), private, parameter :: b2 = 15.0_rp
76  real(rp), private :: c1
77  real(rp), private, parameter :: c2 = 0.75_rp
78  real(rp), private, parameter :: c3 = 0.352_rp
79  real(rp), private, parameter :: c5 = 0.2_rp
80  real(rp), private, parameter :: g1 = 0.235_rp
81  real(rp), private :: g2
82  real(rp), private :: f1
83  real(rp), private :: f2
84  real(rp), private :: rf1
85  real(rp), private :: rf2
86  real(rp), private :: rfc
87  real(rp), private :: af12
88  real(rp), private, parameter :: prn = 0.74_rp
89 
90  real(rp), private, parameter :: zeta_min = -3.0_rp
91  real(rp), private, parameter :: zeta_max = 1.0_rp
92 
93  real(rp), private :: sqrt_2pi
94  real(rp), private :: rsqrt_2pi
95  real(rp), private :: rsqrt_2
96 
97  logical, private :: initialize
98 
99  real(rp), private :: atmos_phy_bl_mynn_pbl_max = 4.e+3_rp
100  real(rp), private :: atmos_phy_bl_mynn_tke_min = 1.e-20_rp
101  real(rp), private :: atmos_phy_bl_mynn_n2_max = 1.e1_rp
102  real(rp), private :: atmos_phy_bl_mynn_nu_min = -1.e1_rp
103  real(rp), private :: atmos_phy_bl_mynn_nu_max = 1.e4_rp
104  real(rp), private :: atmos_phy_bl_mynn_kh_min = -1.e1_rp
105  real(rp), private :: atmos_phy_bl_mynn_kh_max = 1.e4_rp
106  real(rp), private :: atmos_phy_bl_mynn_lt_max = 700.0_rp ! ~ 0.23 * 3 km
107  real(rp), private :: atmos_phy_bl_mynn_sq_fact = 3.0_rp
108  logical, private :: atmos_phy_bl_mynn_init_tke = .false.
109  logical, private :: atmos_phy_bl_mynn_similarity = .true.
110 
111  character(len=H_SHORT), private :: atmos_phy_bl_mynn_level = "2.5" ! "2.5" or "3", level 3 is under experimental yet.
112 
113  namelist / param_atmos_phy_bl_mynn / &
114  atmos_phy_bl_mynn_pbl_max, &
115  atmos_phy_bl_mynn_n2_max, &
116  atmos_phy_bl_mynn_nu_min, &
117  atmos_phy_bl_mynn_nu_max, &
118  atmos_phy_bl_mynn_kh_min, &
119  atmos_phy_bl_mynn_kh_max, &
120  atmos_phy_bl_mynn_lt_max, &
121  atmos_phy_bl_mynn_level, &
122  atmos_phy_bl_mynn_sq_fact, &
123  atmos_phy_bl_mynn_init_tke, &
124  atmos_phy_bl_mynn_similarity
125 
126 
127  !-----------------------------------------------------------------------------
128 contains
129  !-----------------------------------------------------------------------------
133  subroutine atmos_phy_bl_mynn_tracer_setup
134  use scale_prc, only: &
135  prc_abort
136  implicit none
137 
138  integer :: ierr
139 
140  log_newline
141  log_info("ATMOS_PHY_BL_MYNN_tracer_setup",*) 'Tracer Setup'
142  log_info("ATMOS_PHY_BL_MYNN_tracer_setup",*) 'Mellor-Yamada Nakanishi-Niino scheme'
143 
144  !--- read namelist
145  rewind(io_fid_conf)
146  read(io_fid_conf,nml=param_atmos_phy_bl_mynn,iostat=ierr)
147  if( ierr > 0 ) then !--- fatal error
148  log_error("ATMOS_PHY_BL_MYNN_tracer_setup",*) 'Not appropriate names in namelist PARAM_ATMOS_PHY_BL_MYNN. Check!'
149  call prc_abort
150  endif
151 
152  select case ( atmos_phy_bl_mynn_level )
153  case ( "2.5" )
154  atmos_phy_bl_mynn_ntracer = 1
155  allocate( atmos_phy_bl_mynn_name(1), atmos_phy_bl_mynn_desc(1), atmos_phy_bl_mynn_units(1) )
156  atmos_phy_bl_mynn_name(:) = (/ 'TKE_MYNN' /)
157  atmos_phy_bl_mynn_desc(:) = (/ 'turbulent kinetic energy (MYNN)' /)
158  atmos_phy_bl_mynn_units(:) = (/ 'm2/s2' /)
159  case ( "3" )
160  atmos_phy_bl_mynn_ntracer = 4
161  allocate( atmos_phy_bl_mynn_name(4), atmos_phy_bl_mynn_desc(4), atmos_phy_bl_mynn_units(4) )
162  atmos_phy_bl_mynn_name(:) = (/ 'TKE_MYNN', &
163  'TSQ_MYNN', &
164  'QSQ_MYNN', &
165  'COV_MYNN' /)
166  atmos_phy_bl_mynn_desc(:) = (/ 'turbulent kinetic energy (MYNN) ', &
167  'sub-grid variance of liquid water potential temperature (MYNN) ', &
168  'sub-grid variance of total water content (MYNN) ', &
169  'sub-grid covariance of liquid water potential temperature and total water content (MYNN)' /)
170  atmos_phy_bl_mynn_units(:) = (/ 'm2/s2 ', &
171  'K2 ', &
172  'kg2/kg2', &
173  'K kg ' /)
174  case default
175  log_error("ATMOS_PHY_BL_MYNN_tracer_setup",*) 'only level 2.5 and 3 are supported at this moment'
176  call prc_abort
177  end select
178 
179  return
180  end subroutine atmos_phy_bl_mynn_tracer_setup
181 
182  !-----------------------------------------------------------------------------
186  subroutine atmos_phy_bl_mynn_setup( &
187  KA, KS, KE, IA, IS, IE, JA, JS, JE, &
188  BULKFLUX_type, &
189  TKE_MIN, PBL_MAX )
190  use scale_prc, only: &
191  prc_abort
192  use scale_const, only: &
193  pi => const_pi
194  implicit none
195 
196  integer, intent(in) :: ka, ks, ke
197  integer, intent(in) :: ia, is, ie
198  integer, intent(in) :: ja, js, je
199 
200  character(len=*), intent(in) :: bulkflux_type
201 
202  real(rp), intent(in), optional :: tke_min
203  real(rp), intent(in), optional :: pbl_max
204 
205  integer :: ierr
206  integer :: k, i, j
207  !---------------------------------------------------------------------------
208 
209  log_newline
210  log_info("ATMOS_PHY_BL_MYNN_setup",*) 'Setup'
211  log_info("ATMOS_PHY_BL_MYNN_setup",*) 'Mellor-Yamada Nakanishi-Niino scheme'
212 
213  if ( present(tke_min) ) atmos_phy_bl_mynn_tke_min = tke_min
214  if ( present(pbl_max) ) atmos_phy_bl_mynn_pbl_max = pbl_max
215 
216  !--- read namelist
217  rewind(io_fid_conf)
218  read(io_fid_conf,nml=param_atmos_phy_bl_mynn,iostat=ierr)
219  if( ierr < 0 ) then !--- missing
220  log_info("ATMOS_PHY_BL_MYNN_setup",*) 'Not found namelist. Default used.'
221  elseif( ierr > 0 ) then !--- fatal error
222  log_error("ATMOS_PHY_BL_MYNN_setup",*) 'Not appropriate names in namelist PARAM_ATMOS_PHY_BL_MYNN. Check!'
223  call prc_abort
224  endif
225  log_nml(param_atmos_phy_bl_mynn)
226 
227  a1 = b1 * (1.0_rp - 3.0_rp * g1) / 6.0_rp
228  a2 = 1.0_rp / (3.0_rp * g1 * b1**(1.0_rp/3.0_rp) * prn )
229  c1 = g1 - 1.0_rp / ( 3.0_rp * a1 * b1**(1.0_rp/3.0_rp) )
230  g2 = ( 2.0_rp * a1 * (3.0_rp - 2.0_rp * c2) + b2 * (1.0_rp - c3) ) / b1
231  f1 = b1 * (g1 - c1) + 2.0_rp * a1 * (3.0_rp - 2.0_rp * c2) + 3.0_rp * a2 * (1.0_rp - c2) * (1.0_rp - c5)
232  f2 = b1 * (g1 + g2) - 3.0_rp * a1 * (1.0_rp - c2)
233 
234  rf1 = b1 * (g1 - c1) / f1
235  rf2 = b1 * g1 / f2
236  rfc = g1 / (g1 + g2)
237 
238  af12 = a1 * f1 / ( a2 * f2 )
239 
240  sqrt_2pi = sqrt( 2.0_rp * pi )
241  rsqrt_2pi = 1.0_rp / sqrt_2pi
242  rsqrt_2 = 1.0_rp / sqrt( 2.0_rp )
243 
244  if ( atmos_phy_bl_mynn_level == "3" ) then
245  log_warn("ATMOS_PHY_BL_MYNN_setup", *) "At this moment, level 3 is still experimental"
246  end if
247 
248  initialize = atmos_phy_bl_mynn_init_tke
249 
250  select case ( bulkflux_type )
251  case ( "B91", "B91W01" )
252  ! do nothing
253  case default
254  atmos_phy_bl_mynn_similarity = .false.
255  end select
256 
257  return
258  end subroutine atmos_phy_bl_mynn_setup
259 
260  !-----------------------------------------------------------------------------
264  subroutine atmos_phy_bl_mynn_tendency( &
265  KA, KS, KE, IA, IS, IE, JA, JS, JE, &
266  DENS, U, V, POTT, PROG, &
267  PRES, EXNER, N2, &
268  QDRY, QV, Qw, POTL, POTV, SFC_DENS, &
269  SFLX_MU, SFLX_MV, SFLX_SH, SFLX_QV, &
270  us, ts, qs, RLmo, &
271  CZ, FZ, dt_DP, &
272  BULKFLUX_type, &
273  RHOU_t, RHOV_t, RHOT_t, RHOQV_t, &
274  RPROG_t, &
275  Nu, Kh, Qlp, cldfrac, Zi )
276  use scale_const, only: &
277  eps => const_eps, &
278  grav => const_grav, &
279  karman => const_karman, &
280  cpdry => const_cpdry, &
281  epstvap => const_epstvap
282  use scale_prc, only: &
283  prc_abort
284  use scale_atmos_hydrometeor, only: &
285  hydrometeor_lhv => atmos_hydrometeor_lhv
286  use scale_matrix, only: &
287  matrix_solver_tridiagonal
288  use scale_file_history, only: &
289  file_history_in
290  implicit none
291 
292  integer, intent(in) :: ka, ks, ke
293  integer, intent(in) :: ia, is, ie
294  integer, intent(in) :: ja, js, je
295 
296  real(rp), intent(in) :: dens (ka,ia,ja)
297  real(rp), intent(in) :: u (ka,ia,ja)
298  real(rp), intent(in) :: v (ka,ia,ja)
299  real(rp), intent(in) :: pott (ka,ia,ja)
300  real(rp), intent(in) :: prog (ka,ia,ja,atmos_phy_bl_mynn_ntracer)
301  real(rp), intent(in) :: pres (ka,ia,ja)
302  real(rp), intent(in) :: exner (ka,ia,ja)
303  real(rp), intent(in) :: n2 (ka,ia,ja)
304  real(rp), intent(in) :: qdry (ka,ia,ja)
305  real(rp), intent(in) :: qv (ka,ia,ja)
306  real(rp), intent(in) :: qw (ka,ia,ja)
307  real(rp), intent(in) :: potl (ka,ia,ja)
308  real(rp), intent(in) :: potv (ka,ia,ja)
309  real(rp), intent(in) :: sfc_dens( ia,ja)
310  real(rp), intent(in) :: sflx_mu ( ia,ja)
311  real(rp), intent(in) :: sflx_mv ( ia,ja)
312  real(rp), intent(in) :: sflx_sh ( ia,ja)
313  real(rp), intent(in) :: sflx_qv ( ia,ja)
314  real(rp), intent(in) :: us ( ia,ja)
315  real(rp), intent(in) :: ts ( ia,ja)
316  real(rp), intent(in) :: qs ( ia,ja)
317  real(rp), intent(in) :: rlmo ( ia,ja)
318 
319  real(rp), intent(in) :: cz( ka,ia,ja)
320  real(rp), intent(in) :: fz(0:ka,ia,ja)
321  real(dp), intent(in) :: dt_dp
322 
323  character(len=*), intent(in) :: bulkflux_type
324 
325  real(rp), intent(out) :: rhou_t (ka,ia,ja)
326  real(rp), intent(out) :: rhov_t (ka,ia,ja)
327  real(rp), intent(out) :: rhot_t (ka,ia,ja)
328  real(rp), intent(out) :: rhoqv_t(ka,ia,ja)
329  real(rp), intent(out) :: rprog_t(ka,ia,ja,atmos_phy_bl_mynn_ntracer)
330  real(rp), intent(out) :: nu (ka,ia,ja)
331  real(rp), intent(out) :: kh (ka,ia,ja)
332  real(rp), intent(out) :: qlp (ka,ia,ja)
333  real(rp), intent(out) :: cldfrac(ka,ia,ja)
334  real(rp), intent(out) :: zi (ia,ja)
335 
336  real(rp) :: ri (ka,ia,ja)
337  real(rp) :: pr (ka,ia,ja)
338  real(rp) :: prod (ka,ia,ja)
339  real(rp) :: diss (ka,ia,ja)
340  real(rp) :: dudz2(ka,ia,ja)
341  real(rp) :: l (ka,ia,ja)
342  real(rp) :: rho_h(ka)
343 
344  real(rp) :: flxu(ka,ia,ja)
345  real(rp) :: flxv(ka,ia,ja)
346  real(rp) :: flxt(ka,ia,ja)
347  real(rp) :: flxq(ka,ia,ja)
348 
349  real(rp) :: rho (ka)
350  real(rp) :: rhonu (ka)
351  real(rp) :: rhokh (ka)
352  real(rp) :: n2_new(ka)
353  real(rp) :: sflx_pt
354  real(rp) :: sflx_ptv
355  real(rp) :: sm25 (ka)
356  real(rp) :: smp (ka)
357  real(rp) :: sh25 (ka)
358  real(rp) :: shpgh (ka)
359  real(rp) :: nu_f (ka)
360  real(rp) :: kh_f (ka)
361  real(rp) :: q (ka)
362  real(rp) :: q2_2 (ka)
363  real(rp) :: ac (ka)
364  real(rp) :: lq (ka)
365 
366  ! for level 3
367  real(rp) :: tsq (ka)
368  real(rp) :: qsq (ka)
369  real(rp) :: cov (ka)
370  real(rp) :: wtl
371  real(rp) :: wqw
372  real(rp) :: prod_t(ka)
373  real(rp) :: prod_q(ka)
374  real(rp) :: prod_c(ka)
375  real(rp) :: diss_p(ka)
376 
377  real(rp) :: dtldz(ka)
378  real(rp) :: dqwdz(ka)
379  real(rp) :: betat(ka)
380  real(rp) :: betaq(ka)
381 
382  real(rp) :: gammat (ka)
383  real(rp) :: gammaq (ka)
384  real(rp) :: f_gamma(ka)
385 
386  real(rp) :: rlqsm_h(ka)
387 
388  real(rp) :: flx(ka)
389  real(rp) :: a(ka)
390  real(rp) :: b(ka)
391  real(rp) :: c(ka)
392  real(rp) :: d(ka)
393  real(rp) :: ap
394  real(rp) :: phi_n(ka)
395  real(rp) :: tke_p(ka)
396  real(rp) :: dummy(ka)
397 
398  real(rp) :: sf_t
399  real(rp) :: us3
400  real(rp) :: zeta
401  real(rp) :: phi_m, phi_h
402 
403  real(rp) :: fdz(ka)
404  real(rp) :: cdz(ka)
405  real(rp) :: f2h(ka,2)
406  real(rp) :: z1
407 
408  logical :: mynn_level3
409 
410  real(rp) :: dt
411 
412  real(rp) :: fmin
413  integer :: kmin
414 
415  real(rp) :: tmp, sw
416 
417  integer :: ke_pbl
418  integer :: k, i, j
419  integer :: nit, it
420  !---------------------------------------------------------------------------
421 
422  dt = real(dt_dp, rp)
423 
424  log_progress(*) "atmosphere / physics / pbl / MYNN"
425 
426  mynn_level3 = ( atmos_phy_bl_mynn_level == "3" )
427 
428 
429 !OCL INDEPENDENT
430  !$omp parallel do default(none) &
431  !$omp OMP_SCHEDULE_ collapse(2) &
432  !$omp shared(KA,KS,KE,IS,IE,JS,JE, &
433  !$omp EPS,GRAV,CPdry,EPSTvap,UNDEF,RSQRT_2,SQRT_2PI,RSQRT_2PI, &
434  !$omp ATMOS_PHY_BL_MYNN_N2_MAX,ATMOS_PHY_BL_MYNN_TKE_MIN, &
435  !$omp ATMOS_PHY_BL_MYNN_NU_MIN,ATMOS_PHY_BL_MYNN_NU_MAX, &
436  !$omp ATMOS_PHY_BL_MYNN_KH_MIN,ATMOS_PHY_BL_MYNN_KH_MAX, &
437  !$omp ATMOS_PHY_BL_MYNN_Sq_fact,ATMOS_PHY_BL_MYNN_similarity, &
438  !$omp ATMOS_PHY_BL_MYNN_PBL_MAX, &
439  !$omp RHOU_t,RHOV_t,RHOT_t,RHOQV_t,RPROG_t,Nu,Kh,Qlp,cldfrac,Zi, &
440  !$omp DENS,PROG,U,V,POTT,PRES,QDRY,QV,Qw,POTV,POTL,EXNER,N2, &
441  !$omp SFC_DENS,SFLX_MU,SFLX_MV,SFLX_SH,SFLX_QV,us,ts,qs,RLmo, &
442  !$omp mynn_level3,initialize, &
443  !$omp CZ,FZ,dt, &
444  !$omp BULKFLUX_type, &
445  !$omp Ri,Pr,prod,diss,dudz2,l,flxU,flxV,flxT,flxQ) &
446  !$omp private(N2_new,lq,sm25,smp,sh25,shpgh,rlqsm_h,Nu_f,Kh_f,q,q2_2,ac, &
447  !$omp SFLX_PT,SFLX_PTV,RHO,RHONu,RHOKh, &
448  !$omp dtldz,dqwdz,betat,betaq,gammat,gammaq,f_gamma,wtl,wqw, &
449  !$omp flx,a,b,c,d,ap,rho_h,phi_n,tke_P,sf_t,zeta,phi_m,phi_h,us3,CDZ,FDZ,f2h,z1, &
450  !$omp dummy, &
451  !$omp tsq,qsq,cov, &
452  !$omp prod_t,prod_q,prod_c,diss_p, &
453  !$omp fmin,kmin, &
454  !$omp sw,tmp, &
455  !$omp KE_PBL,k,i,j,it,nit)
456  do j = js, je
457  do i = is, ie
458 
459  ke_pbl = ks+1
460  do k = ks+2, ke-1
461  if ( atmos_phy_bl_mynn_pbl_max >= cz(k,i,j) - fz(ks-1,i,j) ) then
462  ke_pbl = k
463  else
464  exit
465  end if
466  end do
467 
468  if ( initialize ) then
469  nit = ke_pbl - 1
470  else
471  nit = 1
472  end if
473 
474 
475  z1 = cz(ks,i,j) - fz(ks-1,i,j)
476 
477  do k = ks, ke_pbl
478  fdz(k) = cz(k+1,i,j) - cz(k ,i,j)
479  end do
480  do k = ks, ke_pbl+1
481  cdz(k) = fz(k ,i,j) - fz(k-1,i,j)
482  end do
483 
484  call get_f2h( &
485  ka, ks, ke_pbl, &
486  cdz(:), & ! (in)
487  f2h(:,:) ) ! (out)
488 
489  call calc_vertical_differece( ka, ks, ke_pbl, &
490  i, j, &
491  u(:,i,j), v(:,i,j), potl(:,i,j), qw(:,i,j), & ! (in)
492  cdz(:), fdz(:), f2h(:,:), & ! (in)
493  dudz2(:,i,j), dtldz(:), dqwdz(:) ) ! (out)
494 
495  do k = ks, ke_pbl
496  n2_new(k) = min( max( n2(k,i,j), - atmos_phy_bl_mynn_n2_max ), atmos_phy_bl_mynn_n2_max )
497 ! n2_new(k) = GRAV * POTV(k,i,j) * dtldz(k)
498  ri(k,i,j) = n2_new(k) / dudz2(k,i,j)
499  end do
500 
501  do k = ks, ke_pbl
502  q(k) = sqrt( max( prog(k,i,j,i_tke), atmos_phy_bl_mynn_tke_min ) * 2.0_rp )
503  end do
504 
505  if ( mynn_level3 ) then
506  do k = ks, ke_pbl
507  tsq(k) = max( prog(k,i,j,i_tsq), 0.0_rp )
508  qsq(k) = max( prog(k,i,j,i_qsq), 0.0_rp )
509  cov(k) = prog(k,i,j,i_cov)
510  cov(k) = sign( min( abs(cov(k)), sqrt(tsq(k) * qsq(k))), cov(k) )
511  end do
512  end if
513 
514 
515  flx(ks-1 ) = 0.0_rp
516  flx(ke_pbl) = 0.0_rp
517 
518  do it = 1, nit
519 
520  if ( initialize .and. it==1 ) then
521  do k = ks, ke_pbl
522  q(k) = 1e10_rp
523  end do
524  end if
525 
526  ! length
527  us3 = us(i,j)**3
528  sflx_ptv = - us3 * rlmo(i,j) / karman
529  call get_length( &
530  ka, ks, ke_pbl, &
531  i, j, &
532  q(:), n2_new(:), & ! (in)
533  sflx_ptv, rlmo(i,j), & ! (in)
534  cz(:,i,j), & ! (in)
535  fz(:,i,j), & ! (in)
536  l(:,i,j) ) ! (out)
537 
538  call get_q2_level2( &
539  ka, ks, ke_pbl, &
540  dudz2(:,i,j), ri(:,i,j), l(:,i,j), & ! (in)
541  q2_2(:) ) ! (out)
542 
543  if ( initialize .and. it==1 ) then
544  do k = ks, ke_pbl
545  q(k) = sqrt( q2_2(k) )
546  end do
547  end if
548 
549  do k = ks, ke_pbl
550  ac(k) = min( q(k) / sqrt( q2_2(k) + 1e-20_rp ), 1.0_rp )
551  end do
552 
553  call get_smsh( &
554  ka, ks, ke_pbl, & ! (in)
555  i, j, & ! (in)
556  q(:), ac(:), & ! (in)
557  l(:,i,j), n2_new(:), & ! (in)
558  potv(:,i,j), dudz2(:,i,j), & ! (in)
559  dtldz(:), dqwdz(:), & ! (in)
560  betat(:), betaq(:), & ! (in)
561  .false., .false., & ! (in)
562  tsq(:), qsq(:), cov(:), & ! (inout)
563  sm25(:), smp(:), & ! (out)
564  sh25(:), shpgh(:), & ! (out) ! dummy
565  gammat(:), gammaq(:), & ! (out) ! dummy
566  f_gamma(:) ) ! (out) ! dummy
567 
568 
569  call partial_condensation( ka, ks, ke_pbl, &
570  pres(:,i,j), pott(:,i,j), & ! (in)
571  potl(:,i,j), qw(:,i,j), & ! (in)
572  qdry(:,i,j), exner(:,i,j), & ! (in)
573  dtldz(:), dqwdz(:), & ! (in)
574  l(:,i,j), sh25(:), ac(:), & ! (in)
575  tsq(:), qsq(:), cov(:), & ! (in)
576  mynn_level3, & ! (in)
577  betat(:), betaq(:), & ! (out)
578  qlp(:,i,j), cldfrac(:,i,j) ) ! (out)
579 
580  if ( atmos_phy_bl_mynn_similarity ) then
581 
582  zeta = min( max( z1 * rlmo(i,j), zeta_min ), zeta_max )
583 
584  select case ( bulkflux_type )
585 !!$ case ( 'B71' )
586 !!$ ! Businger et al. (1971)
587 !!$ if ( zeta >= 0 ) then
588 !!$ phi_m = 4.7_RP * zeta + 1.0_RP
589 !!$ phi_h = 4.7_RP * zeta + 0.74_RP
590 !!$ else
591 !!$ phi_m = 1.0_RP / sqrt(sqrt( 1.0_RP - 15.0_RP * zeta ))
592 !!$ phi_h = 0.47_RP / sqrt( 1.0_RP - 9.0_RP * zeta )
593 !!$ end if
594  case ( 'B91', 'B91W01' )
595  ! Beljaars and Holtslag (1991)
596  if ( zeta >= 0 ) then
597  tmp = - 2.0_rp / 3.0_rp * ( 0.35_rp * zeta - 6.0_rp ) * exp(-0.35_rp*zeta)
598  phi_m = tmp * zeta + zeta + 1.0_rp
599  phi_h = tmp + zeta * sqrt( 1.0_rp + 2.0_rp * zeta / 3.0_rp ) + 1.0_rp
600  else
601  if ( bulkflux_type == 'B91W01' ) then
602  ! Wilson (2001)
603 ! tmp = (-zeta)**(2.0_RP/3.0_RP)
604  tmp = abs(zeta)**(2.0_rp/3.0_rp)
605  phi_m = 1.0_rp / sqrt( 1.0_rp + 3.6_rp * tmp )
606  phi_h = 0.95_rp / sqrt( 1.0_rp + 7.9_rp * tmp )
607  else
608 ! tmp = sqrt( 1.0_RP - 16.0_RP * zeta )
609  tmp = sqrt( 1.0_rp + 16.0_rp * abs(zeta) )
610  phi_m = 1.0_rp / sqrt(tmp)
611  phi_h = 1.0_rp / tmp
612  end if
613  end if
614  end select
615 
616  end if
617 
618  ! update N2
619  do k = ks, ke_pbl
620  n2_new(k) = min(atmos_phy_bl_mynn_n2_max, &
621  grav * ( dtldz(k) * betat(k) + dqwdz(k) * betaq(k) ) / potv(k,i,j) )
622  end do
623 
624  do k = ks, ke_pbl
625  ri(k,i,j) = n2_new(k) / dudz2(k,i,j)
626  end do
627 
628  ! length
629  sflx_pt = sflx_sh(i,j) / ( cpdry * exner(ks,i,j) )
630  sflx_ptv = grav / potv(ks,i,j) * ( betat(ks) * sflx_pt + betaq(ks) * sflx_qv(i,j) ) / sfc_dens(i,j)
631  call get_length( &
632  ka, ks, ke_pbl, &
633  i, j, &
634  q(:), n2_new(:), & ! (in)
635  sflx_ptv, rlmo(i,j), & ! (in)
636  cz(:,i,j), & ! (in)
637  fz(:,i,j), & ! (in)
638  l(:,i,j) ) ! (out)
639 
640  call get_q2_level2( &
641  ka, ks, ke_pbl, &
642  dudz2(:,i,j), ri(:,i,j), l(:,i,j), & ! (in)
643  q2_2(:) ) ! (out)
644 
645  do k = ks, ke_pbl
646  ac(k) = min( q(k) / sqrt( q2_2(k) + 1e-20_rp ), 1.0_rp )
647  end do
648 
649  call get_smsh( &
650  ka, ks, ke_pbl, & ! (in)
651  i, j, & ! (in)
652  q(:), ac(:), & ! (in)
653  l(:,i,j), n2_new(:), & ! (in)
654  potv(:,i,j), dudz2(:,i,j), & ! (in)
655  dtldz(:), dqwdz(:), & ! (in)
656  betat(:), betaq(:), & ! (in)
657  mynn_level3, & ! (in)
658  initialize .and. it==1, & ! (in)
659  tsq(:), qsq(:), cov(:), & ! (inout)
660  sm25(:), smp(:), & ! (out)
661  sh25(:), shpgh(:), & ! (out)
662  gammat(:), gammaq(:), & ! (out)
663  f_gamma(:) ) ! (out)
664 
665  do k = ks, ke_pbl
666  lq(k) = l(k,i,j) * q(k)
667  nu_f(k) = lq(k) * sm25(k)
668  kh_f(k) = lq(k) * sh25(k)
669  end do
670  if ( atmos_phy_bl_mynn_similarity ) then
671 ! Nu_f(KS) = KARMAN * z1 * us(i,j) / phi_m
672 ! Kh_f(KS) = KARMAN * z1 * us(i,j) / phi_h
673  end if
674 
675  do k = ks, ke_pbl-1
676  nu(k,i,j) = min( f2h(k,1) * nu_f(k+1) + f2h(k,2) * nu_f(k), &
677  atmos_phy_bl_mynn_nu_max )
678  kh(k,i,j) = min( f2h(k,1) * kh_f(k+1) + f2h(k,2) * kh_f(k), &
679  atmos_phy_bl_mynn_kh_max )
680  end do
681 
682  do k = ks, ke_pbl-1
683  sw = 0.5_rp - sign(0.5_rp, abs(kh(k,i,j)) - eps)
684  pr(k,i,j) = nu(k,i,j) / ( kh(k,i,j) + sw ) * ( 1.0_rp - sw ) &
685  + 1.0_rp * sw
686  end do
687 
688  rho(ks) = dens(ks,i,j) + dt * sflx_qv(i,j) / cdz(ks)
689  do k = ks+1, ke_pbl
690  rho(k) = dens(k,i,j)
691  end do
692 
693  do k = ks, ke_pbl-1
694  rho_h(k) = f2h(k,1) * rho(k+1) + f2h(k,2) * rho(k)
695  rhonu(k) = rho_h(k) * nu(k,i,j)
696  rhokh(k) = rho_h(k) * kh(k,i,j)
697  end do
698 
699 
700  if ( mynn_level3 ) then
701 
702  ! production term calculated by explicit scheme
703  do k = ks, ke_pbl
704  wtl = - lq(k) * ( sh25(k) * dtldz(k) + gammat(k) )
705  wqw = - lq(k) * ( sh25(k) * dqwdz(k) + gammaq(k) )
706  prod_t(k) = - 2.0_rp * wtl * dtldz(k)
707  prod_q(k) = - 2.0_rp * wqw * dqwdz(k)
708  prod_c(k) = - wtl * dqwdz(k) - wqw * dtldz(k)
709  end do
710 
711  call get_gamma_implicit( &
712  ka, ks, ke_pbl, &
713  i, j, &
714  tsq(:), qsq(:), cov(:), & ! (in)
715  dtldz(:), dqwdz(:), potv(:,i,j), & ! (in)
716  prod_t(:), prod_q(:), prod_c(:), & ! (in)
717  betat(:), betaq(:), & ! (in)
718  f_gamma(:), l(:,i,j), q(:), & ! (in)
719  dt, & ! (in)
720  gammat(:), gammaq(:) ) ! (inout)
721 
722  ! update production terms
723  do k = ks, ke_pbl
724  wtl = - lq(k) * ( sh25(k) * dtldz(k) + gammat(k) )
725  wqw = - lq(k) * ( sh25(k) * dqwdz(k) + gammaq(k) )
726  prod_t(k) = - 2.0_rp * wtl * dtldz(k)
727  prod_q(k) = - 2.0_rp * wqw * dqwdz(k)
728  prod_c(k) = - wtl * dqwdz(k) - wqw * dtldz(k)
729  end do
730 
731  end if
732 
733 
734 
735  ! time integration
736 
737  if ( it == nit ) then
738 
739  do k = ks, ke_pbl-1
740  rlqsm_h(k) = rho_h(k) * ( f2h(k,1) * lq(k+1) * smp(k+1) + f2h(k,2) * lq(k) * smp(k) )
741  end do
742 
743  ! dens * u
744 
745  ! countergradient flux
746  do k = ks, ke_pbl-1
747  flx(k) = - rlqsm_h(k) * ( u(k+1,i,j) - u(k,i,j) ) / fdz(k)
748  end do
749 
750  sf_t = sflx_mu(i,j) / cdz(ks)
751  d(ks) = ( u(ks,i,j) * dens(ks,i,j) + dt * ( sf_t - flx(ks) / cdz(ks) ) ) / rho(ks)
752  do k = ks+1, ke_pbl
753  d(k) = u(k,i,j) - dt * ( flx(k) - flx(k-1) ) / ( cdz(k) * rho(k) )
754  end do
755  c(ks) = 0.0_rp
756  do k = ks, ke_pbl-1
757  ap = - dt * rhonu(k) / fdz(k)
758  a(k) = ap / ( rho(k) * cdz(k) )
759  b(k) = - a(k) - c(k) + 1.0_rp
760  c(k+1) = ap / ( rho(k+1) * cdz(k+1) )
761  end do
762  a(ke_pbl) = 0.0_rp
763  b(ke_pbl) = - c(ke_pbl) + 1.0_rp
764 
765  call matrix_solver_tridiagonal( &
766  ka, ks, ke_pbl, &
767  a(:), b(:), c(:), d(:), & ! (in)
768  dummy(:) ) ! (out)
769 ! phi_n(:) ) ! (out)
770 
771  phi_n(ks:ke_pbl) = dummy(ks:ke_pbl)
772  rhou_t(ks,i,j) = ( phi_n(ks) * rho(ks) - u(ks,i,j) * dens(ks,i,j) ) / dt - sf_t
773  do k = ks+1, ke_pbl
774  rhou_t(k,i,j) = ( phi_n(k) - u(k,i,j) ) * rho(k) / dt
775  end do
776  do k = ke_pbl+1, ke
777  rhou_t(k,i,j) = 0.0_rp
778  end do
779  flxu(ks-1,i,j) = 0.0_rp
780  do k = ks, ke_pbl-1
781  flxu(k,i,j) = flx(k) &
782  - rhonu(k) * ( phi_n(k+1) - phi_n(k) ) / fdz(k)
783  end do
784  do k = ke_pbl, ke
785  flxu(k,i,j) = 0.0_rp
786  end do
787 
788 
789  ! dens * v
790 
791  ! countergradient flux
792  do k = ks, ke_pbl-1
793  flx(k) = - rlqsm_h(k) * ( v(k+1,i,j) - v(k,i,j) ) / fdz(k)
794  end do
795 
796  sf_t = sflx_mv(i,j) / cdz(ks)
797  d(ks) = ( v(ks,i,j) * dens(ks,i,j) + dt * ( sf_t - flx(ks) / cdz(ks) ) ) / rho(ks)
798  do k = ks+1, ke_pbl
799  d(k) = v(k,i,j) - dt * ( flx(k) - flx(k-1) ) / ( cdz(k) * rho(k) )
800  end do
801  ! a,b,c is the same as those for the u
802 
803  call matrix_solver_tridiagonal( &
804  ka, ks, ke_pbl, &
805  a(:), b(:), c(:), d(:), & ! (in)
806  phi_n(:) ) ! (out)
807 
808  rhov_t(ks,i,j) = ( phi_n(ks) * rho(ks) - v(ks,i,j) * dens(ks,i,j) ) / dt - sf_t
809  do k = ks+1, ke_pbl
810  rhov_t(k,i,j) = ( phi_n(k) - v(k,i,j) ) * rho(k) / dt
811  end do
812  do k = ke_pbl+1, ke
813  rhov_t(k,i,j) = 0.0_rp
814  end do
815  flxv(ks-1,i,j) = 0.0_rp
816  do k = ks, ke_pbl-1
817  flxv(k,i,j) = flx(k) &
818  - rhonu(k) * ( phi_n(k+1) - phi_n(k) ) / fdz(k)
819  end do
820  do k = ke_pbl, ke
821  flxv(k,i,j) = 0.0_rp
822  end do
823 
824 
825  ! dens * pott
826 
827  ! countergradient flux
828  do k = ks, ke_pbl-1
829  flx(k) = - ( f2h(k,1) * lq(k+1) * gammat(k+1) + f2h(k,2) * lq(k) * gammat(k) ) &
830  * rho_h(k)
831  end do
832 
833  sf_t = sflx_pt / cdz(ks)
834  ! assume that induced vapor has the same PT
835  d(ks) = pott(ks,i,j) + dt * ( sf_t - flx(ks) / cdz(ks) ) / rho(ks)
836  do k = ks+1, ke_pbl
837  d(k) = pott(k,i,j) - dt * ( flx(k) - flx(k-1) ) / ( cdz(k) * rho(k) )
838  end do
839 
840  c(ks) = 0.0_rp
841  do k = ks, ke_pbl-1
842  ap = - dt * rhokh(k) / fdz(k)
843  a(k) = ap / ( rho(k) * cdz(k) )
844  b(k) = - a(k) - c(k) + 1.0_rp
845  c(k+1) = ap / ( rho(k+1) * cdz(k+1) )
846  end do
847  a(ke_pbl) = 0.0_rp
848  b(ke_pbl) = - c(ke_pbl) + 1.0_rp
849 
850  call matrix_solver_tridiagonal( &
851  ka, ks, ke_pbl, &
852  a(:), b(:), c(:), d(:), & ! (in)
853  phi_n(:) ) ! (out)
854 
855  rhot_t(ks,i,j) = ( phi_n(ks) - pott(ks,i,j) ) * rho(ks) / dt - sf_t
856  do k = ks+1, ke_pbl
857  rhot_t(k,i,j) = ( phi_n(k) - pott(k,i,j) ) * rho(k) / dt
858  end do
859  do k = ke_pbl+1, ke
860  rhot_t(k,i,j) = 0.0_rp
861  end do
862  flxt(ks-1,i,j) = 0.0_rp
863  do k = ks, ke_pbl-1
864  flxt(k,i,j) = flx(k) &
865  - rhokh(k) * ( phi_n(k+1) - phi_n(k) ) / cdz(k)
866  end do
867  do k = ke_pbl, ke
868  flxt(k,i,j) = 0.0_rp
869  end do
870 
871  kmin = ks-1
872  if ( flxt(ks,i,j) > 1e-4_rp ) then
873  fmin = flxt(ks,i,j) / rho_h(ks)
874  do k = ks+1, ke_pbl-2
875  tmp = ( flxt(k-1,i,j) + flxt(k,i,j) + flxt(k+1,i,j) ) &
876  / ( rho_h(k-1) + rho_h(k) + rho_h(k+1) ) ! running mean
877  if ( fmin < 0.0_rp .and. tmp > fmin ) exit
878  if ( tmp < fmin ) then
879  fmin = tmp
880  kmin = k
881  end if
882  end do
883  end if
884  zi(i,j) = fz(kmin,i,j) - fz(ks-1,i,j)
885 
886  ! dens * qv
887 
888  ! countergradient flux
889  do k = ks, ke_pbl-1
890  flx(k) = - ( f2h(k,1) * lq(k+1) * gammaq(k+1) + f2h(k,2) * lq(k) * gammaq(k) ) &
891  * rho_h(k)
892  end do
893 
894  sf_t = sflx_qv(i,j) / cdz(ks)
895  d(ks) = ( qv(ks,i,j) * dens(ks,i,j) + dt * ( sf_t - flx(ks) / cdz(ks) ) ) / rho(ks)
896  do k = ks+1, ke_pbl
897  d(k) = qv(k,i,j) - dt * ( flx(k) - flx(k-1) ) / ( cdz(k) * rho(k) )
898  end do
899 
900  call matrix_solver_tridiagonal( &
901  ka, ks, ke_pbl, &
902  a(:), b(:), c(:), d(:), & ! (in)
903  phi_n(:) ) ! (out)
904 
905  rhoqv_t(ks,i,j) = ( phi_n(ks) * rho(ks) - qv(ks,i,j) * dens(ks,i,j) ) / dt - sf_t
906  do k = ks+1, ke_pbl
907  rhoqv_t(k,i,j) = ( phi_n(k) - qv(k,i,j) ) * rho(k) / dt
908  end do
909  do k = ke_pbl+1, ke
910  rhoqv_t(k,i,j) = 0.0_rp
911  end do
912  flxq(ks-1,i,j) = 0.0_rp
913  do k = ks, ke_pbl-1
914  flxq(k,i,j) = flx(k) &
915  - rhokh(k) * ( phi_n(k+1) - phi_n(k) ) / cdz(k)
916  end do
917  do k = ke_pbl, ke
918  flxq(k,i,j) = 0.0_rp
919  end do
920 
921  end if
922 
923 
924  ! dens * TKE
925 
926  ! production
927  nu(ks-1,i,j) = 0.0_rp
928  kh(ks-1,i,j) = 0.0_rp
929  do k = ks, ke_pbl
930  prod(k,i,j) = lq(k) * ( ( sm25(k) + smp(k) ) * dudz2(k,i,j) &
931  - ( sh25(k) * n2_new(k) - shpgh(k) ) )
932  end do
933  if ( atmos_phy_bl_mynn_similarity ) then
934  prod(ks,i,j) = us3 / ( karman * z1 ) * ( phi_m - zeta )
935  end if
936 
937  do k = ks, ke_pbl
938  tke_p(k) = q(k)**2 * 0.5_rp
939  diss(k,i,j) = - 2.0_rp * q(k) / ( b1 * l(k,i,j) )
940 ! prod(k,i,j) = max( prod(k,i,j), - tke_p(k) / dt - diss(k,i,j) * tke_p(k) )
941  end do
942  do k = ke_pbl+1, ke
943  diss(k,i,j) = 0.0_rp
944  prod(k,i,j) = 0.0_rp
945  end do
946 
947  do k = ks, ke_pbl-1
948  d(k) = ( tke_p(k) * dens(k,i,j) + dt * prod(k,i,j) * rho(k) ) / rho(k)
949  end do
950  d(ke_pbl) = 0.0_rp
951 
952  c(ks) = 0.0_rp
953  do k = ks, ke_pbl-1
954  ap = - dt * atmos_phy_bl_mynn_sq_fact * rhonu(k) / fdz(k)
955  a(k) = ap / ( rho(k) * cdz(k) )
956  b(k) = - a(k) - c(k) + 1.0_rp - diss(k,i,j) * dt
957  c(k+1) = ap / ( rho(k+1) * cdz(k+1) )
958  end do
959  a(ke_pbl) = 0.0_rp
960  b(ke_pbl) = - c(ke_pbl) + 1.0_rp - diss(ke_pbl,i,j) * dt
961 
962  call matrix_solver_tridiagonal( &
963  ka, ks, ke_pbl, &
964  a(:), b(:), c(:), d(:), & ! (in)
965  phi_n(:) ) ! (out)
966 
967  do k = ks, ke_pbl-1
968  phi_n(k) = max( phi_n(k), atmos_phy_bl_mynn_tke_min )
969  end do
970  phi_n(ke_pbl) = 0.0_rp
971 
972 
973  if ( it == nit ) then
974  do k = ks, ke_pbl
975  diss(k,i,j) = diss(k,i,j) * phi_n(k)
976  rprog_t(k,i,j,i_tke) = ( phi_n(k) * rho(k) - prog(k,i,j,i_tke) * dens(k,i,j) ) / dt
977  end do
978  do k = ke_pbl+1, ke
979  rprog_t(k,i,j,i_tke) = 0.0_rp
980  end do
981  else
982  do k = ks, ke_pbl
983  q(k) = sqrt( phi_n(k) * 2.0_rp )
984  end do
985  end if
986 
987 
988  if ( .not. mynn_level3 ) cycle
989 
990 !!$ if ( ATMOS_PHY_BL_MYNN_similarity ) then
991 !!$ tmp = 2.0_RP * us(i,j) * phi_h / ( KARMAN * z1 )
992 !!$ tmp = tmp * ( zeta / ( z1 * RLmo(i,j) ) )**2 ! correspoindint to the limitter for zeta
993 !!$ ! TSQ
994 !!$ prod_t(KS) = tmp * ts(i,j)**2
995 !!$ ! QSQ
996 !!$ prod_q(KS) = tmp * ts(i,j) * qs(i,j)
997 !!$ ! COV
998 !!$ prod_c(KS) = tmp * qs(i,j)**2
999 !!$ end if
1000 
1001  ! dens * tsq
1002 
1003  do k = ks, ke_pbl
1004  diss_p(k) = dt * 2.0_rp * q(k) / ( b2 * l(k,i,j) )
1005  end do
1006  do k = ks, ke_pbl-1
1007  d(k) = ( tsq(k) * dens(k,i,j) + dt * prod_t(k) * rho(k) ) / rho(k)
1008  end do
1009  d(ke_pbl) = 0.0_rp
1010  c(ks) = 0.0_rp
1011  do k = ks, ke_pbl-1
1012  ap = - dt * rhonu(k) / fdz(k)
1013  a(k) = ap / ( rho(k) * cdz(k) )
1014  b(k) = - a(k) - c(k) + 1.0_rp + diss_p(k)
1015  c(k+1) = ap / ( rho(k+1) * cdz(k+1) )
1016  end do
1017  a(ke_pbl) = 0.0_rp
1018  b(ke_pbl) = - c(ke_pbl) + 1.0_rp + diss_p(ke_pbl)
1019 
1020  call matrix_solver_tridiagonal( &
1021  ka, ks, ke_pbl, &
1022  a(:), b(:), c(:), d(:), & ! (in)
1023  tsq(:) ) ! (out)
1024 
1025  do k = ks, ke_pbl-1
1026  tsq(k) = max( tsq(k), 0.0_rp )
1027  end do
1028  tsq(ke_pbl) = 0.0_rp
1029 
1030  if ( it == nit ) then
1031  do k = ks, ke_pbl
1032  rprog_t(k,i,j,i_tsq) = ( tsq(k) * rho(k) - prog(k,i,j,i_tsq) * dens(k,i,j) ) / dt
1033  end do
1034  do k = ke_pbl+1, ke
1035  rprog_t(k,i,j,i_tsq) = 0.0_rp
1036  end do
1037  end if
1038 
1039 
1040  ! dens * qsq
1041 
1042  do k = ks, ke_pbl-1
1043  d(k) = ( qsq(k) * dens(k,i,j) + dt * prod_q(k) * rho(k) ) / rho(k)
1044  end do
1045  d(ke_pbl) = 0.0_rp
1046  ! a, b, c are same as those for tsq
1047 
1048  call matrix_solver_tridiagonal( &
1049  ka, ks, ke_pbl, &
1050  a(:), b(:), c(:), d(:), & ! (in)
1051  qsq(:) ) ! (out)
1052 
1053  do k = ks, ke_pbl-1
1054  qsq(k) = max( qsq(k), 0.0_rp )
1055  end do
1056  qsq(ke_pbl) = 0.0_rp
1057 
1058  if ( it == nit ) then
1059  do k = ks, ke_pbl
1060  rprog_t(k,i,j,i_qsq) = ( qsq(k) * rho(k) - prog(k,i,j,i_qsq) * dens(k,i,j) ) / dt
1061  end do
1062  do k = ke_pbl+1, ke
1063  rprog_t(k,i,j,i_qsq) = 0.0_rp
1064  end do
1065  end if
1066 
1067 
1068  ! dens * cov
1069 
1070  do k = ks, ke_pbl-1
1071  d(k) = ( cov(k) * dens(k,i,j) + dt * prod_c(k) * rho(k) ) / rho(k)
1072  end do
1073  d(ke_pbl) = 0.0_rp
1074  ! a, b, c are same as those for tsq
1075 
1076  call matrix_solver_tridiagonal( &
1077  ka, ks, ke_pbl, &
1078  a(:), b(:), c(:), d(:), & ! (in)
1079  cov(:) ) ! (out)
1080 
1081  do k = ks, ke_pbl-1
1082  cov(k) = sign( min( abs(cov(k)), sqrt(tsq(k)*qsq(k)) ), cov(k) )
1083  end do
1084  cov(ke_pbl) = 0.0_rp
1085 
1086  if ( it == nit ) then
1087  do k = ks, ke_pbl
1088  rprog_t(k,i,j,i_cov) = ( cov(k) * rho(k) - prog(k,i,j,i_cov) * dens(k,i,j) ) / dt
1089  end do
1090  do k = ke_pbl+1, ke
1091  rprog_t(k,i,j,i_cov) = 0.0_rp
1092  end do
1093  end if
1094 
1095  end do
1096 
1097  do k = ke_pbl, ke
1098  nu(k,i,j) = 0.0_rp
1099  kh(k,i,j) = 0.0_rp
1100  pr(k,i,j) = 1.0_rp
1101  end do
1102  do k = ke_pbl+1, ke
1103  ri(k,i,j) = undef
1104  prod(k,i,j) = undef
1105  diss(k,i,j) = undef
1106  dudz2(k,i,j) = undef
1107  l(k,i,j) = undef
1108  qlp(k,i,j) = undef
1109  cldfrac(k,i,j) = undef
1110  end do
1111 
1112  end do
1113  end do
1114 
1115 
1116  call file_history_in(ri(:,:,:), 'Ri_MYNN', 'Richardson number', '1', fill_halo=.true. )
1117  call file_history_in(pr(:,:,:), 'Pr_MYNN', 'Prandtl number', '1', fill_halo=.true., dim_type="ZHXY" )
1118  call file_history_in(prod(:,:,:), 'TKE_prod_MYNN', 'TKE production', 'm2/s3', fill_halo=.true.)
1119  call file_history_in(diss(:,:,:), 'TKE_diss_MYNN', 'TKE dissipation', 'm2/s3', fill_halo=.true.)
1120  call file_history_in(dudz2(:,:,:), 'dUdZ2_MYNN', 'dudz2', 'm2/s2', fill_halo=.true.)
1121  call file_history_in(l(:,:,:), 'L_mix_MYNN', 'minxing length', 'm', fill_halo=.true.)
1122 
1123  call file_history_in(flxu(:,:,:), 'ZFLX_RHOU_MYNN', 'Z FLUX of RHOU (MYNN)', 'kg/m/s2', fill_halo=.true., dim_type="ZHXY" )
1124  call file_history_in(flxv(:,:,:), 'ZFLX_RHOV_MYNN', 'Z FLUX of RHOV (MYNN)', 'kg/m/s2', fill_halo=.true., dim_type="ZHXY" )
1125  call file_history_in(flxt(:,:,:), 'ZFLX_RHOT_MYNN', 'Z FLUX of RHOT (MYNN)', 'K kg/m2/s', fill_halo=.true., dim_type="ZHXY" )
1126  call file_history_in(flxq(:,:,:), 'ZFLX_QV_MYNN', 'Z FLUX of RHOQV (MYNN)', 'kg/m2/s', fill_halo=.true., dim_type="ZHXY" )
1127 
1128 
1129  initialize = .false.
1130 
1131  return
1132  end subroutine atmos_phy_bl_mynn_tendency
1133 
1134  !-----------------------------------------------------------------------------
1138  subroutine atmos_phy_bl_mynn_tendency_tracer( &
1139  KA, KS, KE, IA, IS, IE, JA, JS, JE, &
1140  DENS, QTRC, SFLX_Q, &
1141  Kh, MASS, &
1142  CZ, FZ, DDT, &
1143  TRACER_NAME, &
1144  RHOQ_t )
1145  use scale_matrix, only: &
1146  matrix_solver_tridiagonal
1147  use scale_file_history, only: &
1148  file_history_in
1149  integer, intent(in) :: ka, ks, ke
1150  integer, intent(in) :: ia, is, ie
1151  integer, intent(in) :: ja, js, je
1152 
1153  real(rp), intent(in) :: dens (ka,ia,ja)
1154  real(rp), intent(in) :: qtrc (ka,ia,ja)
1155  real(rp), intent(in) :: sflx_q( ia,ja)
1156  real(rp), intent(in) :: kh (ka,ia,ja)
1157  real(rp), intent(in) :: mass
1158  real(rp), intent(in) :: cz ( ka,ia,ja)
1159  real(rp), intent(in) :: fz (0:ka,ia,ja)
1160  real(dp), intent(in) :: ddt
1161  character(len=*), intent(in) :: tracer_name
1162 
1163  real(rp), intent(out) :: rhoq_t(ka,ia,ja)
1164 
1165  real(rp) :: qtrc_n(ka)
1166  real(rp) :: rho (ka)
1167  real(rp) :: rhokh (ka)
1168  real(rp) :: a(ka)
1169  real(rp) :: b(ka)
1170  real(rp) :: c(ka)
1171  real(rp) :: d(ka)
1172  real(rp) :: rho_h
1173  real(rp) :: ap
1174  real(rp) :: sf_t
1175 
1176  real(rp) :: flx(ka,ia,ja)
1177 
1178  real(rp) :: cdz(ka)
1179  real(rp) :: fdz(ka)
1180  real(rp) :: f2h(ka,2)
1181 
1182  real(rp) :: dt
1183 
1184  integer :: ke_pbl
1185  integer :: k, i, j
1186 
1187  dt = real( ddt, kind=rp )
1188 
1189 !OCL INDEPENDENT
1190  !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
1191  !$omp shared(KA,KS,KE,IS,IE,JS,JE) &
1192  !$omp shared(ATMOS_PHY_BL_MYNN_PBL_MAX) &
1193  !$omp shared(RHOQ_t,DENS,QTRC,SFLX_Q,Kh,MASS,CZ,FZ,DT,flx) &
1194  !$omp private(QTRC_n,RHO,RHOKh,rho_h,a,b,c,d,ap,sf_t,CDZ,FDZ,f2h) &
1195  !$omp private(KE_PBL,k,i,j)
1196  do j = js, je
1197  do i = is, ie
1198 
1199  ke_pbl = ks+1
1200  do k = ks+2, ke-1
1201  if ( atmos_phy_bl_mynn_pbl_max >= cz(k,i,j) - fz(ks-1,i,j) ) then
1202  ke_pbl = k
1203  else
1204  exit
1205  end if
1206  end do
1207 
1208  do k = ks, ke_pbl
1209  cdz(k) = fz(k ,i,j) - fz(k-1,i,j)
1210  fdz(k) = cz(k+1,i,j) - cz(k ,i,j)
1211  end do
1212 
1213  call get_f2h( &
1214  ka, ks, ke_pbl-1, &
1215  cdz(:), & ! (in)
1216  f2h(:,:) ) ! (out)
1217 
1218  sf_t = sflx_q(i,j) / cdz(ks)
1219 
1220  rho(ks) = dens(ks,i,j) + dt * sf_t * mass
1221  do k = ks+1, ke_pbl
1222  rho(k) = dens(k,i,j)
1223  end do
1224 
1225  ! dens * coefficient at the half level
1226  do k = ks, ke_pbl-1
1227  rho_h = f2h(k,1) * dens(k+1,i,j) + f2h(k,2) * dens(k,i,j)
1228  rhokh(k) = rho_h * kh(k,i,j)
1229  end do
1230 
1231  d(ks) = ( qtrc(ks,i,j) * dens(ks,i,j) + dt * sf_t ) / rho(ks)
1232  do k = ks+1, ke_pbl
1233  d(k) = qtrc(k,i,j)
1234  end do
1235 
1236  c(ks) = 0.0_rp
1237  do k = ks, ke_pbl-1
1238  ap = - dt * rhokh(k) / fdz(k)
1239  a(k) = ap / ( rho(k) * cdz(k) )
1240  b(k) = - a(k) - c(k) + 1.0_rp
1241  c(k+1) = ap / ( rho(k+1) * cdz(k+1) )
1242  end do
1243  a(ke_pbl) = 0.0_rp
1244  b(ke_pbl) = - c(ke_pbl) + 1.0_rp
1245 
1246  call matrix_solver_tridiagonal( &
1247  ka, ks, ke_pbl, &
1248  a(:), b(:), c(:), d(:), & ! (in)
1249  qtrc_n(:) ) ! (out)
1250 
1251  rhoq_t(ks,i,j) = ( qtrc_n(ks) * rho(ks) - qtrc(ks,i,j) * dens(ks,i,j) ) / dt - sf_t
1252  do k = ks+1, ke_pbl
1253  rhoq_t(k,i,j) = ( qtrc_n(k) - qtrc(k,i,j) ) * rho(k) / dt
1254  end do
1255  do k = ke_pbl+1, ke
1256  rhoq_t(k,i,j) = 0.0_rp
1257  end do
1258 
1259  do k = ks, ke_pbl-1
1260  flx(k,i,j) = - rhokh(k) * ( qtrc_n(k+1) - qtrc_n(k) ) / fdz(k)
1261  end do
1262 
1263  end do
1264  end do
1265 
1266  call file_history_in(flx(:,:,:), 'ZFLX_'//trim(tracer_name)//'_MYNN', 'Z FLUX of DENS * '//trim(tracer_name)//' (MYNN)', 'kg/m2/s', fill_halo=.true.)
1267 
1268  return
1269  end subroutine atmos_phy_bl_mynn_tendency_tracer
1270 
1271  !-----------------------------------------------------------------------------
1272  ! private routines
1273  !-----------------------------------------------------------------------------
1274 
1275 !OCL SERIAL
1276  subroutine get_length( &
1277  KA, KS, KE_PBL, &
1278  i, j, &
1279  q, n2, &
1280  SFLX_PTV, RLmo, &
1281  CZ, &
1282  FZ, &
1283  l )
1284  use scale_const, only: &
1285  grav => const_grav, &
1286  karman => const_karman, &
1287  eps => const_eps
1288  implicit none
1289  integer, intent(in) :: ka, ks, ke_pbl
1290  integer, intent(in) :: i, j ! for debug
1291 
1292  real(rp), intent(in) :: q(ka)
1293  real(rp), intent(in) :: n2(ka)
1294  real(rp), intent(in) :: sflx_ptv
1295  real(rp), intent(in) :: rlmo
1296  real(rp), intent(in) :: cz(ka)
1297  real(rp), intent(in) :: fz(0:ka)
1298 
1299  real(rp), intent(out) :: l(ka)
1300 
1301  real(rp), parameter :: ls_fact_max = 2.0_rp
1302 
1303  real(rp) :: ls
1304  real(rp) :: lb
1305  real(rp) :: lt
1306  real(rp) :: rlt
1307 
1308  real(rp) :: qc
1309  real(rp) :: int_q
1310  real(rp) :: int_qz
1311  real(rp) :: rn2sr
1312  real(rp) :: zeta
1313 
1314  real(rp) :: z
1315  real(rp) :: qdz
1316 
1317  real(rp) :: sw
1318  integer :: k
1319 
1320  int_qz = 0.0_rp
1321  int_q = 0.0_rp
1322  do k = ks, ke_pbl
1323  qdz = q(k) * ( fz(k) - fz(k-1) )
1324  z = cz(k) - fz(ks-1)
1325  int_qz = int_qz + z * qdz
1326  int_q = int_q + qdz
1327  end do
1328  ! LT
1329  lt = min( max(0.23_rp * int_qz / (int_q + 1e-20_rp), &
1330  lt_min), &
1331  atmos_phy_bl_mynn_lt_max )
1332  rlt = 1.0_rp / lt
1333 
1334  qc = ( lt * max(sflx_ptv,0.0_rp) )**oneoverthree ! qc=0 if SFLX_PTV<0
1335 
1336  do k = ks, ke_pbl
1337  z = cz(k) - fz(ks-1)
1338  zeta = min( max( z * rlmo, zeta_min ), zeta_max )
1339 
1340  ! LS
1341  sw = sign(0.5_rp, zeta) + 0.5_rp ! 1 for zeta >= 0, 0 for zeta < 0
1342  ls = karman * z &
1343  * ( sw / (1.0_rp + 2.7_rp*zeta*sw ) &
1344  + min( ( (1.0_rp - 100.0_rp*zeta)*(1.0_rp-sw) )**0.2_rp, ls_fact_max) )
1345 
1346  ! LB
1347  sw = sign(0.5_rp, n2(k)-eps) + 0.5_rp ! 1 for dptdz >0, 0 for dptdz <= 0
1348  rn2sr = 1.0_rp / ( sqrt(n2(k)*sw) + 1.0_rp-sw)
1349  lb = (1.0_rp + 5.0_rp * sqrt(qc*rn2sr*rlt)) * q(k) * rn2sr * sw & ! qc=0 when RLmo > 0
1350  + 1.e10_rp * (1.0_rp-sw)
1351 
1352  ! L
1353  l(k) = 1.0_rp / ( 1.0_rp/ls + rlt + 1.0_rp/(lb+1e-20_rp) )
1354  end do
1355 
1356  return
1357  end subroutine get_length
1358 
1359 !OCL SERIAL
1360  subroutine get_q2_level2( &
1361  KA, KS, KE_PBL, &
1362  dudz2, Ri, l, &
1363  q2_2 )
1364  implicit none
1365  integer, intent(in) :: ka, ks, ke_pbl
1366 
1367  real(rp), intent(in) :: dudz2(ka)
1368  real(rp), intent(in) :: ri(ka)
1369  real(rp), intent(in) :: l(ka)
1370 
1371  real(rp), intent(out) :: q2_2(ka)
1372 
1373  real(rp) :: rf
1374  real(rp) :: sm_2
1375  real(rp) :: sh_2
1376 
1377  integer :: k
1378 
1379  do k = ks, ke_pbl
1380  rf = min(0.5_rp / af12 * ( ri(k) &
1381  + af12*rf1 &
1382  - sqrt(ri(k)**2 + 2.0_rp*af12*(rf1-2.0_rp*rf2)*ri(k) + (af12*rf1)**2) ), &
1383  rfc)
1384  sh_2 = 3.0_rp * a2 * (g1+g2) * (rfc-rf) / (1.0_rp-rf)
1385  sm_2 = sh_2 * af12 * (rf1-rf) / (rf2-rf)
1386  q2_2(k) = b1 * l(k)**2 * sm_2 * (1.0_rp-rf) * dudz2(k)
1387  end do
1388 
1389  return
1390  end subroutine get_q2_level2
1391 
1392 !OCL SERIAL
1393  subroutine get_smsh( &
1394  KA, KS, KE_PBL, &
1395  i, j, &
1396  q, ac, &
1397  l, n2, &
1398  potv, dudz2, &
1399  dtldz, dqwdz, &
1400  betat, betaq, &
1401  mynn_level3, &
1402  initialize, &
1403  tsq, qsq, cov, &
1404  sm25, smp, &
1405  sh25, shpgh, &
1406  gammat, gammaq, &
1407  f_gamma )
1408  use scale_const, only: &
1409  eps => const_eps, &
1410  grav => const_grav
1411  implicit none
1412  integer, intent(in) :: ka, ks, ke_pbl
1413  integer, intent(in) :: i, j ! for debug
1414 
1415  real(rp), intent(in) :: q(ka)
1416  real(rp), intent(in) :: ac(ka)
1417  real(rp), intent(in) :: l(ka)
1418  real(rp), intent(in) :: n2(ka)
1419  real(rp), intent(in) :: potv(ka)
1420  real(rp), intent(in) :: dudz2(ka)
1421  real(rp), intent(in) :: dtldz(ka)
1422  real(rp), intent(in) :: dqwdz(ka)
1423  real(rp), intent(in) :: betat(ka)
1424  real(rp), intent(in) :: betaq(ka)
1425  logical, intent(in) :: mynn_level3
1426  logical, intent(in) :: initialize
1427 
1428  real(rp), intent(inout) :: tsq(ka)
1429  real(rp), intent(inout) :: qsq(ka)
1430  real(rp), intent(inout) :: cov(ka)
1431 
1432  real(rp), intent(out) :: sm25 (ka) ! S_M2.5
1433  real(rp), intent(out) :: smp (ka) ! S_M'
1434  real(rp), intent(out) :: sh25 (ka) ! S_H2.5
1435  real(rp), intent(out) :: shpgh (ka) ! S_H' G_H * (q/L)^2
1436  real(rp), intent(out) :: gammat (ka) ! \Gamma_\theta
1437  real(rp), intent(out) :: gammaq (ka) ! \Gamma_q
1438  real(rp), intent(out) :: f_gamma(ka) ! - E_H / q^2 * GRAV / Theta_0
1439 
1440  real(rp) :: l2
1441  real(rp) :: q2
1442  real(rp) :: ac2
1443  real(rp) :: p1q2
1444  real(rp) :: p2q2
1445  real(rp) :: p3q2
1446  real(rp) :: p4q2
1447  real(rp) :: p5q2
1448  real(rp) :: rd25
1449  real(rp) :: ghq2
1450  real(rp) :: gmq2
1451  real(rp) :: f1, f2, f3, f4
1452 
1453  ! for level 3
1454  real(rp) :: tvsq
1455  real(rp) :: tvsq25
1456  real(rp) :: tltv
1457  real(rp) :: tltv25
1458  real(rp) :: qwtv
1459  real(rp) :: qwtv25
1460  real(rp) :: tsq25
1461  real(rp) :: qsq25
1462  real(rp) :: cov25
1463  real(rp) :: emq2
1464  real(rp) :: eh
1465  real(rp) :: ew
1466  real(rp) :: rdp
1467  real(rp) :: cw25q2
1468  real(rp) :: fact
1469 
1470  integer :: k
1471 
1472  ! level 2.5
1473  do k = ks, ke_pbl
1474 
1475  ac2 = ac(k)**2
1476  l2 = l(k)**2
1477  q2 = q(k)**2
1478 
1479  f1 = - 3.0_rp * ac2 * a2 * b2 * ( 1.0_rp - c3 )
1480  f2 = - 9.0_rp * ac2 * a1 * a2 * ( 1.0_rp - c2 )
1481  f3 = 9.0_rp * ac2 * a2**2 * ( 1.0_rp - c2 ) * ( 1.0_rp - c5 )
1482  f4 = - 12.0_rp * ac2 * a1 * a2 * ( 1.0_rp - c2 )
1483 
1484  ghq2 = max( - n2(k) * l2, -q2 ) ! L/q <= 1/N for N^2>0
1485 
1486  gmq2 = dudz2(k) * l2
1487 
1488  p1q2 = q2 + f1 * ghq2
1489  p2q2 = q2 + f2 * ghq2
1490  p3q2 = q2 + ( f1 + f3 ) * ghq2
1491  p4q2 = q2 + ( f1 + f4 ) * ghq2
1492  p5q2 = 6.0_rp * ac2 * a1**2 * gmq2
1493  rd25 = q2 / max( p2q2 * p4q2 + p5q2 * p3q2, 1e-20_rp )
1494 
1495  sm25(k) = ac(k) * a1 * ( p3q2 - 3.0_rp * c1 * p4q2 ) * rd25
1496  sh25(k) = ac(k) * a2 * ( p2q2 + 3.0_rp * c1 * p5q2 ) * rd25
1497 
1498 
1499  if ( mynn_level3 ) then ! level 3
1500 
1501  fact = ac(k) * b2 * l2 * sh25(k)
1502  tsq25 = fact * dtldz(k)**2
1503  qsq25 = fact * dqwdz(k)**2
1504  cov25 = fact * dtldz(k) * dqwdz(k)
1505 
1506  ! initialize tsq, qsq, and cov with those of level 2.5
1507  if ( initialize ) then
1508  tsq(k) = tsq25
1509  qsq(k) = qsq25
1510  cov(k) = cov25
1511  end if
1512 
1513  if ( q2 <= 1e-10_rp ) then
1514  smp(k) = 0.0_rp
1515  shpgh(k) = 0.0_rp
1516  f_gamma(k) = 0.0_rp
1517  gammat(k) = 0.0_rp
1518  gammaq(k) = 0.0_rp
1519  else
1520  tltv25 = betat(k) * tsq25 + betaq(k) * cov25
1521  qwtv25 = betat(k) * cov25 + betaq(k) * qsq25
1522  tvsq25 = betat(k) * tltv25 + betaq(k) * qwtv25
1523 
1524  tltv = betat(k) * tsq(k) + betaq(k) * cov(k)
1525  qwtv = betat(k) * cov(k) + betaq(k) * qsq(k)
1526  tvsq = betat(k) * tltv + betaq(k) * qwtv
1527 
1528  rdp = q2 / max( p2q2 * ( f4 * ghq2 + q2 ) + p5q2 * ( f3 * ghq2 + q2 ), 1e-20_rp )
1529 
1530  fact = l2 / q2 * ( grav / potv(k) )**2 * ( tvsq - tvsq25 )
1531 
1532  ew = ( 1.0_rp - c3 ) * ( - p2q2 * f4 - p5q2 * f3 ) * rdp ! (p1-p4)/gh = -f4, (p1-p3)/gh = -f3
1533  if ( abs(ew) > 1e-20_rp ) then
1534  cw25q2 = p1q2 * ( p2q2 + 3.0_rp * c1 * p5q2 ) * rd25 / 3.0_rp
1535  fact = fact * ew
1536  fact = min( max( fact, 0.12_rp * q2 - cw25q2 ), 0.76_rp * q2 - cw25q2 )
1537  fact = fact / ew
1538  end if
1539 
1540  emq2 = 3.0_rp * ac(k) * a1 * ( 1.0_rp - c3 ) * ( f3 - f4 ) * rdp ! (p3-p4)/gh = (f3-f4)
1541  eh = 3.0_rp * ac(k) * a2 * ( 1.0_rp - c3 ) * ( p2q2 + p5q2 ) * rdp
1542 
1543  smp(k) = emq2 * fact
1544  shpgh(k) = eh * fact / l2
1545  f_gamma(k) = - eh * grav / ( q2 * potv(k) )
1546  gammat(k) = ( tltv - tltv25 ) * f_gamma(k)
1547  gammaq(k) = ( qwtv - qwtv25 ) * f_gamma(k)
1548  end if
1549  else ! level 2.5
1550  smp(k) = 0.0_rp
1551  shpgh(k) = 0.0_rp
1552  f_gamma(k) = 0.0_rp
1553  gammat(k) = 0.0_rp
1554  gammaq(k) = 0.0_rp
1555  end if
1556 
1557  end do
1558 
1559  return
1560  end subroutine get_smsh
1561 
1562 !OCL SERIAL
1563  subroutine get_gamma_implicit( &
1564  KA, KS, KE, &
1565  i, j, &
1566  tsq, qsq, cov, &
1567  dtldz, dqwdz, POTV, &
1568  prod_t, prod_q, prod_c, &
1569  betat, betaq, &
1570  f_gamma, l, q, &
1571  dt, &
1572  gammat, gammaq )
1573  use scale_const, only: &
1574  grav => const_grav
1575  implicit none
1576  integer, intent(in) :: KA, KS, KE
1577  integer, intent(in) :: i, j ! for debug
1578 
1579  real(RP), intent(in) :: tsq (KA)
1580  real(RP), intent(in) :: qsq (KA)
1581  real(RP), intent(in) :: cov (KA)
1582  real(RP), intent(in) :: dtldz (KA)
1583  real(RP), intent(in) :: dqwdz (KA)
1584  real(RP), intent(in) :: POTV (KA)
1585  real(RP), intent(in) :: prod_t (KA)
1586  real(RP), intent(in) :: prod_q (KA)
1587  real(RP), intent(in) :: prod_c (KA)
1588  real(RP), intent(in) :: betat (KA)
1589  real(RP), intent(in) :: betaq (KA)
1590  real(RP), intent(in) :: f_gamma(KA)
1591  real(RP), intent(in) :: l (KA)
1592  real(RP), intent(in) :: q (KA)
1593  real(RP), intent(in) :: dt
1594 
1595  real(RP), intent(inout) :: gammat(KA)
1596  real(RP), intent(inout) :: gammaq(KA)
1597 
1598  real(RP) :: dtsq
1599  real(RP) :: dqsq
1600  real(RP) :: dcov
1601  real(RP) :: a11, a12, a21, a22, a23, a32, a33
1602  real(RP) :: v1, v2, v3
1603  real(RP) :: f1, f2
1604  real(RP) :: a13, a31, det
1605 
1606  integer :: k
1607 
1608  ! calculate gamma by implicit scheme
1609 
1610  do k = ks, ke
1611 
1612  ! matrix coefficient
1613  f1 = f_gamma(k) * l(k) * q(k)
1614  f2 = - 2.0_rp * q(k) / ( b2 * l(k) )
1615 
1616  v1 = prod_t(k) + f2 * tsq(k)
1617  v2 = prod_c(k) + f2 * cov(k)
1618  v3 = prod_q(k) + f2 * qsq(k)
1619 
1620  a11 = - f1 * dtldz(k) * betat(k) * 2.0_rp + 1.0_rp / dt - f2
1621  a12 = - f1 * dtldz(k) * betaq(k) * 2.0_rp
1622  a21 = - f1 * dqwdz(k) * betat(k)
1623  a22 = - f1 * ( dtldz(k) * betat(k) + dqwdz(k) * betaq(k) ) + 1.0_rp / dt - f2
1624  a23 = - f1 * dtldz(k) * betaq(k)
1625  a32 = - f1 * dqwdz(k) * betat(k) * 2.0_rp
1626  a33 = - f1 * dqwdz(k) * betaq(k) * 2.0_rp + 1.0_rp / dt - f2
1627 
1628  if ( q(k) < 1e-20_rp .or. abs(f_gamma(k)) * dt >= q(k) ) then
1629  ! solve matrix
1630  f1 = a21 / a11
1631  f2 = a23 / a33
1632  dcov = ( v2 - f1 * v1 - f2 * v3 ) / ( a22 - f1 * a12 - f2 * a32 )
1633  dtsq = ( v1 - a12 * dcov ) / a11
1634  dqsq = ( v3 - a32 * dcov ) / a33
1635  else
1636  ! consider change in q
1637  ! dq = - L * G / Theta0 * ( betat * dgammat + betaq * dgammaq )
1638  f1 = - l(k) * grav / potv(k) * f_gamma(k) / q(k)
1639  f2 = f1 * betat(k)**2
1640  a11 = a11 - f2 * v1
1641  a21 = a21 - f2 * v2
1642  a31 = - f2 * v3
1643  f2 = f1 * betat(k) * betaq(k) * 2.0_rp
1644  a12 = a12 - f2 * v1
1645  a22 = a22 - f2 * v2
1646  a32 = a32 - f2 * v3
1647  f2 = f1 * betaq(k)**2
1648  a13 = - f2 * v1
1649  a23 = a23 - f2 * v2
1650  a33 = a33 - f2 * v3
1651  ! solve matrix by Cramer's fomula
1652  det = a11 * ( a22 * a33 - a23 * a32 ) + a12 * ( a23 * a31 - a21 * a33 ) + a13 * ( a21 * a32 - a22 * a31 )
1653  dtsq = ( v1 * ( a22 * a33 - a23 * a32 ) + a12 * ( a23 * v3 - v2 * a33 ) + a13 * ( v2 * a32 - a22 * v3 ) ) / det
1654  dcov = ( a11 * ( v2 * a33 - a23 * v3 ) + v1 * ( a23 * a31 - a21 * a33 ) + a13 * ( a21 * v3 - v2 * a31 ) ) / det
1655  dqsq = ( a11 * ( a22 * v3 - v2 * a32 ) + a12 * ( v2 * a31 - a21 * v3 ) + v3 * ( a21 * a32 - a22 * a31 ) ) / det
1656  end if
1657 
1658  ! modifiy the explicit solution
1659  gammat(k) = gammat(k) + f_gamma(k) * ( betat(k) * dtsq + betaq(k) * dcov )
1660  gammaq(k) = gammaq(k) + f_gamma(k) * ( betat(k) * dcov + betaq(k) * dqsq )
1661  end do
1662 
1663  return
1664  end subroutine get_gamma_implicit
1665 
1666 !OCL SERIAL
1667  subroutine calc_vertical_differece( &
1668  KA, KS, KE, &
1669  i, j, &
1670  U, V, POTL, Qw, &
1671  CDZ, FDZ, F2H, &
1672  dudz2, dtldz, dqwdz )
1673  use scale_const, only: &
1674  eps => const_eps
1675  integer, intent(in) :: KA, KS, KE
1676  integer, intent(in) :: i, j ! for debug
1677 
1678  real(RP), intent(in) :: U (KA)
1679  real(RP), intent(in) :: V (KA)
1680  real(RP), intent(in) :: POTL(KA)
1681  real(RP), intent(in) :: Qw (KA)
1682  real(RP), intent(in) :: CDZ (KA)
1683  real(RP), intent(in) :: FDZ (KA)
1684  real(RP), intent(in) :: F2H (KA,2)
1685 
1686  real(RP), intent(out) :: dudz2(KA)
1687  real(RP), intent(out) :: dtldz(KA)
1688  real(RP), intent(out) :: dqwdz(KA)
1689 
1690  real(RP) :: Uh(KA), Vh(KA)
1691  real(RP) :: qh(KA)
1692 
1693  integer :: k
1694 
1695  do k = ks, ke
1696  uh(k) = f2h(k,1) * u(k+1) + f2h(k,2) * u(k)
1697  vh(k) = f2h(k,1) * v(k+1) + f2h(k,2) * v(k)
1698  end do
1699 
1700  dudz2(ks) = ( ( uh(ks) - u(ks) )**2 + ( vh(ks) - v(ks) )**2 ) / ( cdz(ks) * 0.5_rp )**2
1701 ! dudz2(KS) = ( ( Uh(KS) )**2 + ( Vh(KS) )**2 ) / CDZ(KS)**2
1702  dudz2(ks) = max( dudz2(ks), 1e-20_rp )
1703  do k = ks+1, ke
1704  dudz2(k) = ( ( uh(k) - uh(k-1) )**2 + ( vh(k) - vh(k-1) )**2 ) / cdz(k)**2
1705  dudz2(k) = max( dudz2(k), 1e-20_rp )
1706  end do
1707 
1708 
1709  do k = ks, ke
1710  qh(k) = f2h(k,1) * potl(k+1) + f2h(k,2) * potl(k)
1711  end do
1712  dtldz(ks) = ( qh(ks) - potl(ks) ) / ( cdz(ks) * 0.5_rp )
1713 ! dtldz(KS) = ( POTL(KS+1) - POTL(KS) ) / FDZ(KS)
1714  do k = ks+1, ke
1715  dtldz(k) = ( qh(k) - qh(k-1) ) / cdz(k)
1716  end do
1717 
1718  do k = ks, ke
1719  qh(k) = f2h(k,1) * qw(k+1) + f2h(k,2) * qw(k)
1720  end do
1721  dqwdz(ks) = ( qh(ks) - qw(ks) ) / ( cdz(ks) * 0.5_rp )
1722 ! dqwdz(KS) = ( Qw(KS+1) - Qw(KS) ) / FDZ(KS)
1723  do k = ks+1, ke
1724  dqwdz(k) = ( qh(k) - qh(k-1) ) / cdz(k)
1725  end do
1726 
1727  return
1728  end subroutine calc_vertical_differece
1729 
1730 !OCL SERIAL
1731  subroutine partial_condensation( &
1732  KA, KS, KE, &
1733  PRES, POTT, POTL, Qw, Qdry, EXNER, &
1734  dtldz, dqwdz, l, sh25, ac, &
1735  tsq, qsq, cov, &
1736  mynn_level3, &
1737  betat, betaq, Qlp, cldfrac )
1738  use scale_const, only: &
1739  cpdry => const_cpdry, &
1740  rvap => const_rvap, &
1741  cpvap => const_cpvap, &
1742  epsvap => const_epsvap, &
1743  epstvap => const_epstvap
1744  use scale_atmos_saturation, only: &
1745  atmos_saturation_pres2qsat => atmos_saturation_pres2qsat_liq
1746  use scale_atmos_hydrometeor, only: &
1747  hydrometeor_lhv => atmos_hydrometeor_lhv
1748  integer, intent(in) :: KA, KS, KE
1749  real(RP), intent(in) :: PRES (KA)
1750  real(RP), intent(in) :: POTT (KA)
1751  real(RP), intent(in) :: POTL (KA)
1752  real(RP), intent(in) :: Qw (KA)
1753  real(RP), intent(in) :: Qdry (KA)
1754  real(RP), intent(in) :: EXNER(KA)
1755  real(RP), intent(in) :: dtldz(KA)
1756  real(RP), intent(in) :: dqwdz(KA)
1757  real(RP), intent(in) :: l (KA)
1758  real(RP), intent(in) :: sh25 (KA)
1759  real(RP), intent(in) :: ac (KA)
1760  real(RP), intent(in) :: tsq (KA)
1761  real(RP), intent(in) :: qsq (KA)
1762  real(RP), intent(in) :: cov (KA)
1763  logical, intent(in) :: mynn_level3
1764 
1765  real(RP), intent(out) :: betat (KA)
1766  real(RP), intent(out) :: betaq (KA)
1767  real(RP), intent(out) :: Qlp (KA)
1768  real(RP), intent(out) :: cldfrac(KA)
1769 
1770  real(RP) :: TEML(KA)
1771  real(RP) :: Qsl (KA)
1772  real(RP) :: LHVL(KA)
1773  real(RP) :: dQsl
1774  real(RP) :: CPtot
1775  real(RP) :: aa, bb, cc
1776  real(RP) :: sigma_s
1777  real(RP) :: Q1
1778  real(RP) :: Rt
1779 
1780  integer :: k
1781 
1782  do k = ks, ke
1783  teml(k) = potl(k) * exner(k)
1784  end do
1785 
1786  call atmos_saturation_pres2qsat( &
1787  ka, ks, ke, &
1788  teml(:), pres(:), qdry(:), & ! (in)
1789  qsl(:) ) ! (out)
1790 
1791  call hydrometeor_lhv( &
1792  ka, ks, ke, & ! (in)
1793  teml(:), & ! (in)
1794  lhvl(:) ) ! (out)
1795 
1796  do k = ks, ke
1797 
1798  dqsl = ( 1.0_rp + qsl(k) / ( qdry(k) * epsvap ) ) * qsl(k) * lhvl(k) / ( rvap * teml(k)**2 )
1799  cptot = qdry(k) * cpdry + qsl(k) * cpvap
1800  aa = 1.0_rp / ( 1.0_rp + lhvl(k)/cptot * dqsl )
1801  bb = exner(k) * dqsl
1802 
1803  if ( mynn_level3 ) then
1804  sigma_s = 0.5_rp * aa * sqrt( max( qsq(k) - 2.0_rp * bb * cov(k) + bb**2 * tsq(k), 1.0e-20_rp ) )
1805  else
1806  ! level 2.5
1807  sigma_s = max( 0.5_rp * aa * l(k) * sqrt( ac(k) * b2 * sh25(k) ) * abs( dqwdz(k) - bb * dtldz(k) ), &
1808  1.0e-10_rp )
1809  end if
1810 
1811  q1 = aa * ( qw(k) - qsl(k) ) * 0.5_rp / sigma_s
1812  cldfrac(k) = min( max( 0.5_rp * ( 1.0_rp + erf(q1*rsqrt_2) ), 0.0_rp ), 1.0_rp )
1813  qlp(k) = min( max( 2.0_rp * sigma_s * ( cldfrac(k) * q1 + rsqrt_2pi &
1814 #if defined(PGI) || defined(SX)
1815  * exp( -min( 0.5_rp*q1**2, 1.e+3_rp ) ) & ! apply exp limiter
1816 #else
1817  * exp(-0.5_rp*q1**2) &
1818 #endif
1819  ), 0.0_rp ), qw(k) * 0.5_rp )
1820  cc = ( 1.0_rp + epstvap * qw(k) - qlp(k) / epsvap ) / exner(k) * lhvl(k) / cptot &
1821  - pott(k) / epsvap
1822  rt = cldfrac(k) - qlp(k) / (2.0_rp*sigma_s*sqrt_2pi) &
1823 #if defined(PGI) || defined(SX)
1824  * exp( -min( q1**2 * 0.5_rp, 1.e+3_rp ) ) ! apply exp limiter
1825 #else
1826  * exp(-q1**2 * 0.5_rp)
1827 #endif
1828  betat(k) = 1.0_rp + epstvap * qw(k) - qlp(k) / epsvap - rt * aa * bb * cc
1829  betaq(k) = epstvap * pott(k) + rt * aa * cc
1830 
1831  end do
1832 
1833  return
1834  end subroutine partial_condensation
1835 
1836 !OCL SERIAL
1837  subroutine get_f2h( &
1838  KA, KS, KE, &
1839  CDZ, &
1840  f2h )
1841  integer, intent(in) :: KA, KS, KE
1842  real(RP), intent(in) :: CDZ(KA)
1843  real(RP), intent(out) :: f2h(KA,2)
1844 
1845  real(RP) :: dz1, dz2
1846  integer :: k
1847 
1848  do k = ks, ke
1849  dz1 = cdz(k+1)
1850  dz2 = cdz(k)
1851  f2h(k,1) = dz2 / ( dz1 + dz2 )
1852  f2h(k,2) = dz1 / ( dz1 + dz2 )
1853  end do
1854 
1855  return
1856  end subroutine get_f2h
1857 
1858 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:46
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:52
scale_prc::prc_abort
subroutine, public prc_abort
Abort Process.
Definition: scale_prc.F90:342
scale_const::const_epstvap
real(rp), public const_epstvap
1 / epsilon - 1
Definition: scale_const.F90:70
scale_atmos_phy_bl_mynn::atmos_phy_bl_mynn_tendency
subroutine, public atmos_phy_bl_mynn_tendency(KA, KS, KE, IA, IS, IE, JA, JS, JE, DENS, U, V, POTT, PROG, PRES, EXNER, N2, QDRY, QV, Qw, POTL, POTV, SFC_DENS, SFLX_MU, SFLX_MV, SFLX_SH, SFLX_QV, us, ts, qs, RLmo, CZ, FZ, dt_DP, BULKFLUX_type, RHOU_t, RHOV_t, RHOT_t, RHOQV_t, RPROG_t, Nu, Kh, Qlp, cldfrac, Zi)
ATMOS_PHY_BL_MYNN_tendency calculate tendency by the virtical eddy viscosity.
Definition: scale_atmos_phy_bl_mynn.F90:276
scale_const::const_undef2
integer, parameter, public const_undef2
undefined value (INT2)
Definition: scale_const.F90:38
scale_const::const_epsvap
real(rp), public const_epsvap
Rdry / Rvap.
Definition: scale_const.F90:69
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:53
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:18
scale_atmos_phy_bl_mynn::get_f2h
subroutine get_f2h(KA, KS, KE, CDZ, f2h)
Definition: scale_atmos_phy_bl_mynn.F90:1841
scale_const::const_rvap
real(rp), parameter, public const_rvap
specific gas constant (water vapor) [J/kg/K]
Definition: scale_const.F90:63
scale_tracer::mass
real(rp), public mass
Definition: scale_tracer.F90:46
scale_const::const_eps
real(rp), public const_eps
small number
Definition: scale_const.F90:33
scale_atmos_hydrometeor
module atmosphere / hydrometeor
Definition: scale_atmos_hydrometeor.F90:12
scale_file_history
module file_history
Definition: scale_file_history.F90:15
scale_const::const_cpvap
real(rp), parameter, public const_cpvap
specific heat (water vapor, constant pressure) [J/kg/K]
Definition: scale_const.F90:64
scale_prc
module PROCESS
Definition: scale_prc.F90:11
scale_const::const_pi
real(rp), public const_pi
pi
Definition: scale_const.F90:31
scale_precision::rp
integer, parameter, public rp
Definition: scale_precision.F90:41
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:44
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:56
scale_const::const_cpdry
real(rp), public const_cpdry
specific heat (dry air,constant pressure) [J/kg/K]
Definition: scale_const.F90:56
scale_tracer::tracer_name
character(len=h_short), dimension(qa_max), public tracer_name
Definition: scale_tracer.F90:38
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_phy_bl_mynn::calc_vertical_differece
subroutine calc_vertical_differece(KA, KS, KE, i, j, U, V, POTL, Qw, CDZ, FDZ, F2H, dudz2, dtldz, dqwdz)
Definition: scale_atmos_phy_bl_mynn.F90:1673
scale_atmos_grid_cartesc_index::ja
integer, public ja
Definition: scale_atmos_grid_cartesC_index.F90:49
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:50
scale_matrix
module MATRIX
Definition: scale_matrix.F90:11
scale_atmos_phy_bl_mynn::atmos_phy_bl_mynn_ntracer
integer, public atmos_phy_bl_mynn_ntracer
Definition: scale_atmos_phy_bl_mynn.F90:50
scale_atmos_phy_bl_mynn::atmos_phy_bl_mynn_setup
subroutine, public atmos_phy_bl_mynn_setup(KA, KS, KE, IA, IS, IE, JA, JS, JE, BULKFLUX_type, TKE_MIN, PBL_MAX)
ATMOS_PHY_BL_MYNN_setup Setup.
Definition: scale_atmos_phy_bl_mynn.F90:190
scale_const::const_undef
real(rp), public const_undef
Definition: scale_const.F90:41
scale_atmos_phy_bl_mynn::get_gamma_implicit
subroutine get_gamma_implicit(KA, KS, KE, i, j, tsq, qsq, cov, dtldz, dqwdz, POTV, prod_t, prod_q, prod_c, betat, betaq, f_gamma, l, q, dt, gammat, gammaq)
Definition: scale_atmos_phy_bl_mynn.F90:1573
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:51
scale_io::io_fid_conf
integer, public io_fid_conf
Config file ID.
Definition: scale_io.F90:56
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_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:134
scale_atmos_phy_bl_mynn::partial_condensation
subroutine partial_condensation(KA, KS, KE, PRES, POTT, POTL, Qw, Qdry, EXNER, dtldz, dqwdz, l, sh25, ac, tsq, qsq, cov, mynn_level3, betat, betaq, Qlp, cldfrac)
Definition: scale_atmos_phy_bl_mynn.F90:1738
Generated by   doxygen 1.8.17