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 
71  real(RP), private :: a1
72  real(RP), private :: a2
73  real(RP), private, parameter :: b1 = 24.0_rp
74  real(RP), private, parameter :: b2 = 15.0_rp
75  real(RP), private :: c1
76  real(RP), private, parameter :: c2 = 0.75_rp
77  real(RP), private, parameter :: c3 = 0.352_rp
78  real(RP), private, parameter :: c5 = 0.2_rp
79  real(RP), private, parameter :: g1 = 0.235_rp
80  real(RP), private :: g2
81  real(RP), private :: f1
82  real(RP), private :: f2
83  real(RP), private :: rf1
84  real(RP), private :: rf2
85  real(RP), private :: rfc
86  real(RP), private :: af12
87  real(RP), private, parameter :: prn = 0.74_rp
88 
89  real(RP), private :: sqrt_2pi
90  real(RP), private :: rsqrt_2pi
91  real(RP), private :: rsqrt_2
92 
93  integer, private :: ke_pbl
94 
95  real(RP), private :: atmos_phy_bl_mynn_pbl_max = 1.e+10_rp
96  real(RP), private :: atmos_phy_bl_mynn_tke_min = 1.e-10_rp
97  real(RP), private :: atmos_phy_bl_mynn_n2_max = 1.e+3_rp
98  real(RP), private :: atmos_phy_bl_mynn_nu_min = 1.e-6_rp
99  real(RP), private :: atmos_phy_bl_mynn_nu_max = 10000.0_rp
100  real(RP), private :: atmos_phy_bl_mynn_kh_min = 1.e-6_rp
101  real(RP), private :: atmos_phy_bl_mynn_kh_max = 10000.0_rp
102  real(RP), private :: atmos_phy_bl_mynn_lt_max = 700.0_rp ! ~ 0.23 * 3 km
103 
104  character(len=H_SHORT), private :: atmos_phy_bl_mynn_level = "2.5" ! "2.5" or "3"
105 
106  namelist / param_atmos_phy_bl_mynn / &
107  atmos_phy_bl_mynn_pbl_max, &
108  atmos_phy_bl_mynn_n2_max, &
109  atmos_phy_bl_mynn_nu_min, &
110  atmos_phy_bl_mynn_nu_max, &
111  atmos_phy_bl_mynn_kh_min, &
112  atmos_phy_bl_mynn_kh_max, &
113  atmos_phy_bl_mynn_lt_max, &
114  atmos_phy_bl_mynn_level
115 
116 
117  !-----------------------------------------------------------------------------
118 contains
119  !-----------------------------------------------------------------------------
123  subroutine atmos_phy_bl_mynn_tracer_setup( )
124  use scale_prc, only: &
125  prc_abort
126 
127  integer :: ierr
128 
129  log_newline
130  log_info("ATMOS_PHY_BL_MYNN_tracer_setup",*) 'Tracer Setup'
131  log_info("ATMOS_PHY_BL_MYNN_tracer_setup",*) 'Mellor-Yamada Nakanishi-Niino scheme'
132 
133  !--- read namelist
134  rewind(io_fid_conf)
135  read(io_fid_conf,nml=param_atmos_phy_bl_mynn,iostat=ierr)
136  if( ierr > 0 ) then !--- fatal error
137  log_error("ATMOS_PHY_BL_MYNN_tracer_setup",*) 'Not appropriate names in namelist PARAM_ATMOS_PHY_BL_MYNN. Check!'
138  call prc_abort
139  endif
140 
141  select case ( atmos_phy_bl_mynn_level )
142  case ( "2.5" )
143  atmos_phy_bl_mynn_ntracer = 1
144  allocate( atmos_phy_bl_mynn_name(1), atmos_phy_bl_mynn_desc(1), atmos_phy_bl_mynn_units(1) )
145  atmos_phy_bl_mynn_name(:) = (/ 'TKE_MYNN' /)
146  atmos_phy_bl_mynn_desc(:) = (/ 'turbulent kinetic energy (MYNN)' /)
147  atmos_phy_bl_mynn_units(:) = (/ 'm2/s2' /)
148  case ( "3" )
149  atmos_phy_bl_mynn_ntracer = 4
150  allocate( atmos_phy_bl_mynn_name(4), atmos_phy_bl_mynn_desc(4), atmos_phy_bl_mynn_units(4) )
151  atmos_phy_bl_mynn_name(:) = (/ 'TKE_MYNN', &
152  'TSQ_MYNN', &
153  'QSQ_MYNN', &
154  'COV_MYNN' /)
155  atmos_phy_bl_mynn_desc(:) = (/ 'turbulent kinetic energy (MYNN) ', &
156  'sub-grid variance of liquid water potential temperature (MYNN) ', &
157  'sub-grid variance of total water content (MYNN) ', &
158  'sub-grid covariance of liquid water potential temperature and total water content (MYNN)' /)
159  atmos_phy_bl_mynn_units(:) = (/ 'm2/s2 ', &
160  'K2 ', &
161  'kg2/kg2', &
162  'K kg ' /)
163  case default
164  log_error("ATMOS_PHY_BL_MYNN_tracer_setup",*) 'only level 2.5 and 3 are supported at this moment'
165  call prc_abort
166  end select
167 
168  return
169  end subroutine atmos_phy_bl_mynn_tracer_setup
170 
171  !-----------------------------------------------------------------------------
175  subroutine atmos_phy_bl_mynn_setup( &
176  KA, KS, KE, IA, IS, IE, JA, JS, JE, &
177  CZ, &
178  TKE_MIN, PBL_MAX )
179  use scale_prc, only: &
180  prc_abort
181  use scale_const, only: &
182  pi => const_pi
183  implicit none
184 
185  integer, intent(in) :: ka, ks, ke
186  integer, intent(in) :: ia, is, ie
187  integer, intent(in) :: ja, js, je
188  real(RP), intent(in) :: cz (ka,ia,ja)
189 
190  real(RP), intent(in), optional :: tke_min
191  real(RP), intent(in), optional :: pbl_max
192 
193  integer :: ierr
194  integer :: k, i, j
195  !---------------------------------------------------------------------------
196 
197  log_newline
198  log_info("ATMOS_PHY_BL_MYNN_setup",*) 'Setup'
199  log_info("ATMOS_PHY_BL_MYNN_setup",*) 'Mellor-Yamada Nakanishi-Niino scheme'
200 
201  if ( present(tke_min) ) atmos_phy_bl_mynn_tke_min = tke_min
202  if ( present(pbl_max) ) atmos_phy_bl_mynn_pbl_max = pbl_max
203 
204  !--- read namelist
205  rewind(io_fid_conf)
206  read(io_fid_conf,nml=param_atmos_phy_bl_mynn,iostat=ierr)
207  if( ierr < 0 ) then !--- missing
208  log_info("ATMOS_PHY_BL_MYNN_setup",*) 'Not found namelist. Default used.'
209  elseif( ierr > 0 ) then !--- fatal error
210  log_error("ATMOS_PHY_BL_MYNN_setup",*) 'Not appropriate names in namelist PARAM_ATMOS_PHY_BL_MYNN. Check!'
211  call prc_abort
212  endif
213  log_nml(param_atmos_phy_bl_mynn)
214 
215  a1 = b1 * (1.0_rp - 3.0_rp * g1) / 6.0_rp
216  a2 = 1.0_rp / (3.0_rp * g1 * b1**(1.0_rp/3.0_rp) * prn )
217  c1 = g1 - 1.0_rp / ( 3.0_rp * a1 * b1**(1.0_rp/3.0_rp) )
218  g2 = ( 2.0_rp * a1 * (3.0_rp - 2.0_rp * c2) + b2 * (1.0_rp - c3) ) / b1
219  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)
220  f2 = b1 * (g1 + g2) - 3.0_rp * a1 * (1.0_rp - c2)
221 
222  rf1 = b1 * (g1 - c1) / f1
223  rf2 = b1 * g1 / f2
224  rfc = g1 / (g1 + g2)
225 
226  af12 = a1 * f1 / ( a2 * f2 )
227 
228  sqrt_2pi = sqrt( 2.0_rp * pi )
229  rsqrt_2pi = 1.0_rp / sqrt_2pi
230  rsqrt_2 = 1.0_rp / sqrt( 2.0_rp )
231 
232  do k = ks, ke-1
233  do j = js, je
234  do i = is, ie
235  if ( atmos_phy_bl_mynn_pbl_max >= cz(k,i,j) ) then
236  ke_pbl = k
237  end if
238  end do
239  end do
240  end do
241 
242  if ( atmos_phy_bl_mynn_level == "3" ) then
243  log_warn("ATMOS_PHY_BL_MYNN_setup", *) "At this moment, level 3 is still experimental"
244  end if
245 
246  return
247  end subroutine atmos_phy_bl_mynn_setup
248 
249  !-----------------------------------------------------------------------------
253  subroutine atmos_phy_bl_mynn_tendency( &
254  KA, KS, KE, IA, IS, IE, JA, JS, JE, &
255  DENS, U, V, POTT, PROG, &
256  PRES, EXNER, N2, &
257  QDRY, QV, Qw, POTL, POTV, &
258  SFLX_MU, SFLX_MV, SFLX_SH, SFLX_QV, &
259  l_mo, &
260  CZ, FZ, dt_DP, &
261  RHOU_t, RHOV_t, RHOT_t, RPROG_t, &
262  Nu, Kh )
263  use scale_const, only: &
264  grav => const_grav, &
265  karman => const_karman, &
266  rdry => const_rdry, &
267  rvap => const_rvap, &
268  cpdry => const_cpdry, &
269  cpvap => const_cpvap, &
270  epstvap => const_epstvap
271  use scale_prc, only: &
272  prc_abort
273  use scale_atmos_hydrometeor, only: & !! TODO
274  hydrometeor_lhv => atmos_hydrometeor_lhv, &
275  lhv
276  use scale_matrix, only: &
277  matrix_solver_tridiagonal
278  use scale_atmos_saturation, only: &
279  atmos_saturation_dens2qsat => atmos_saturation_dens2qsat_all
280  use scale_file_history, only: &
281  file_history_in
282  implicit none
283 
284  integer, intent(in) :: ka, ks, ke
285  integer, intent(in) :: ia, is, ie
286  integer, intent(in) :: ja, js, je
287 
288  real(RP), intent(in) :: dens (ka,ia,ja)
289  real(RP), intent(in) :: u (ka,ia,ja)
290  real(RP), intent(in) :: v (ka,ia,ja)
291  real(RP), intent(in) :: pott (ka,ia,ja)
292  real(RP), intent(in) :: prog (ka,ia,ja,atmos_phy_bl_mynn_ntracer)
293  real(RP), intent(in) :: pres (ka,ia,ja)
294  real(RP), intent(in) :: exner (ka,ia,ja)
295  real(RP), intent(in) :: n2 (ka,ia,ja)
296  real(RP), intent(in) :: qdry (ka,ia,ja)
297  real(RP), intent(in) :: qv (ka,ia,ja)
298  real(RP), intent(in) :: qw (ka,ia,ja)
299  real(RP), intent(in) :: potl (ka,ia,ja)
300  real(RP), intent(in) :: potv (ka,ia,ja)
301  real(RP), intent(in) :: sflx_mu( ia,ja)
302  real(RP), intent(in) :: sflx_mv( ia,ja)
303  real(RP), intent(in) :: sflx_sh( ia,ja)
304  real(RP), intent(in) :: sflx_qv( ia,ja)
305  real(RP), intent(in) :: l_mo ( ia,ja)
306 
307  real(RP), intent(in) :: cz( ka,ia,ja)
308  real(RP), intent(in) :: fz(0:ka,ia,ja)
309  real(DP), intent(in) :: dt_dp
310 
311  real(RP), intent(out) :: rhou_t(ka,ia,ja)
312  real(RP), intent(out) :: rhov_t(ka,ia,ja)
313  real(RP), intent(out) :: rhot_t(ka,ia,ja)
314  real(RP), intent(out) :: rprog_t(ka,ia,ja,atmos_phy_bl_mynn_ntracer)
315  real(RP), intent(out) :: nu (ka,ia,ja)
316  real(RP), intent(out) :: kh (ka,ia,ja)
317 
318  real(RP) :: ri (ka,ia,ja)
319  real(RP) :: pr (ka,ia,ja)
320  real(RP) :: prod (ka,ia,ja)
321  real(RP) :: diss (ka,ia,ja)
322  real(RP) :: dudz2(ka,ia,ja)
323  real(RP) :: l (ka,ia,ja)
324 
325  real(RP) :: flxu(ka,ia,ja)
326  real(RP) :: flxv(ka,ia,ja)
327  real(RP) :: flxt(ka,ia,ja)
328 
329  real(RP) :: teml (ka)
330  real(RP) :: rhonu (ka)
331  real(RP) :: rhokh (ka)
332  real(RP) :: lhvl (ka)
333  real(RP) :: cptot
334  real(RP) :: n2_new(ka)
335  real(RP) :: sflx_pt
336  real(RP) :: sm (ka)
337  real(RP) :: sh (ka)
338  real(RP) :: q (ka)
339  real(RP) :: q2_2 (ka)
340  real(RP) :: ac (ka)
341 
342  ! for level 3
343  real(RP) :: tvsq (ka)
344  real(RP) :: tvsq25(ka)
345  real(RP) :: tsq (ka)
346  real(RP) :: qsq (ka)
347  real(RP) :: cov (ka)
348  real(RP) :: tsq25
349  real(RP) :: qsq25
350  real(RP) :: cov25
351  real(RP) :: tltv
352  real(RP) :: qwtv
353  real(RP) :: prod_t1
354  real(RP) :: prod_q1
355  real(RP) :: prod_c1
356 
357  real(RP) :: qlp
358  real(RP) :: q1
359  real(RP) :: qsl(ka)
360  real(RP) :: dqsl
361  real(RP) :: dtldz(ka)
362  real(RP) :: dqwdz(ka)
363  real(RP) :: sigma_s
364  real(RP) :: rr
365  real(RP) :: rt
366  real(RP) :: betat
367  real(RP) :: betaq
368  real(RP) :: aa
369  real(RP) :: bb
370  real(RP) :: cc
371 
372  real(RP) :: a(ka)
373  real(RP) :: b(ka)
374  real(RP) :: c(ka)
375  real(RP) :: d(ka)
376  real(RP) :: ap
377  real(RP) :: phi_n(ka)
378  real(RP) :: dummy(ka)
379  real(RP) :: tke_p
380 
381  real(RP) :: sf_t
382  real(RP) :: us, us3, zeta, phi_m, phi_h
383 
384  real(RP) :: f2h(ka,2)
385  real(RP) :: z1
386 
387  logical :: mynn_level3
388 
389  real(RP) :: dt
390 
391  integer :: k, i, j
392  !---------------------------------------------------------------------------
393 
394  dt = real(dt_dp, rp)
395 
396  log_progress(*) "atmosphere / physics / pbl / MYNN"
397 
398  mynn_level3 = ( atmos_phy_bl_mynn_level == "3" )
399 
400 !OCL INDEPENDENT
401  !$omp parallel do default(none) &
402  !$omp OMP_SCHEDULE_ collapse(2) &
403  !$omp shared(KA,KS,KE_PBL,KE,IS,IE,JS,JE, &
404  !$omp GRAV,CPdry,LHV,EPSTvap,UNDEF,RSQRT_2,SQRT_2PI,RSQRT_2PI, &
405  !$omp ATMOS_PHY_BL_MYNN_N2_MAX,ATMOS_PHY_BL_MYNN_TKE_MIN, &
406  !$omp ATMOS_PHY_BL_MYNN_NU_MIN,ATMOS_PHY_BL_MYNN_NU_MAX, &
407  !$omp ATMOS_PHY_BL_MYNN_KH_MIN,ATMOS_PHY_BL_MYNN_KH_MAX, &
408  !$omp RHOU_t,RHOV_t,RHOT_t,RPROG_t,Nu,Kh, &
409  !$omp DENS,PROG,U,V,POTT,PRES,QDRY,QV,Qw,POTV,POTL,EXNER,N2,SFLX_MU,SFLX_MV,SFLX_SH,SFLX_QV,l_mo, &
410  !$omp mynn_level3, &
411  !$omp CZ,FZ,dt, &
412  !$omp Ri,Pr,prod,diss,dudz2,l,flxU,flxV,flxT) &
413  !$omp private(N2_new,sm,sh,q,q2_2,ac,SFLX_PT,TEML,RHONu,RHOKh,LHVL,CPtot,qlp, &
414  !$omp Q1,Qsl,dQsl,dtldz,dqwdz,sigma_s,RR,Rt,betat,betaq,aa,bb,cc, &
415  !$omp a,b,c,d,ap,phi_n,tke_P,sf_t,zeta,phi_m,phi_h,us,us3,f2h,z1,dummy, &
416  !$omp tvsq,tsq,qsq,cov,tvsq25,tsq25,qsq25,cov25,tltv,qwtv,prod_t1,prod_q1,prod_c1, &
417  !$omp k,i,j)
418  do j = js, je
419  do i = is, ie
420 
421  z1 = cz(ks,i,j) - fz(ks-1,i,j)
422 
423  sflx_pt = sflx_sh(i,j) / ( cpdry * dens(ks,i,j) * exner(ks,i,j) )
424 
425 ! dudz2(KS,i,j) = ( ( U(KS+1,i,j) )**2 + ( V(KS+1,i,j) )**2 ) &
426 ! / ( CZ(KS+1,i,j) - FZ(KS-1,i,j) )**2
427  dudz2(ks,i,j) = ( ( u(ks+1,i,j) - u(ks,i,j) )**2 + ( v(ks+1,i,j) - v(ks,i,j) )**2 ) &
428  / ( cz(ks+1,i,j) - cz(ks,i,j) )**2
429  do k = ks+1, ke_pbl
430  dudz2(k,i,j) = ( ( u(k+1,i,j) - u(k-1,i,j) )**2 + ( v(k+1,i,j) - v(k-1,i,j) )**2 ) &
431  / ( cz(k+1,i,j) - cz(k-1,i,j) )**2
432  end do
433  do k = ke_pbl+1, ke
434  dudz2(k,i,j) = undef
435  end do
436 
437  do k = ks, ke_pbl
438  n2_new(k) = min( atmos_phy_bl_mynn_n2_max, n2(k,i,j) )
439  ri(k,i,j) = n2_new(k) / max(dudz2(k,i,j), 1e-10_rp)
440  end do
441  do k = ke_pbl+1, ke
442  ri(k,i,j) = undef
443  end do
444 
445  do k = ks, ke_pbl
446  q(k) = sqrt( max(prog(k,i,j,i_tke), atmos_phy_bl_mynn_tke_min)*2.0_rp )
447  end do
448 
449  call get_f2h( &
450  ka, ks, ke, &
451  fz(:,i,j), & ! (in)
452  f2h(:,:) ) ! (out)
453 
454  dtldz(ks) = ( potl(ks+1,i,j) - potl(ks,i,j) ) / ( cz(ks+1,i,j) - cz(ks,i,j) )
455  dqwdz(ks) = ( qw(ks+1,i,j) - qw(ks,i,j) ) / ( cz(ks+1,i,j) - cz(ks,i,j) )
456  do k = ks+1, ke_pbl
457  dtldz(k) = ( potl(k+1,i,j) - potl(k-1,i,j) ) / ( cz(k+1,i,j) - cz(k-1,i,j) )
458  dqwdz(k) = ( qw(k+1,i,j) - qw(k-1,i,j) ) / ( cz(k+1,i,j) - cz(k-1,i,j) )
459  end do
460 
461 
462 
463  ! length
464  call get_length( &
465  ka, ks, ke_pbl, &
466  pott(ks,i,j), q(:), n2_new(:), & ! (in)
467  sflx_pt, l_mo(i,j), & ! (in)
468  fz(:,i,j), & ! (in)
469  l(:,i,j) ) ! (out)
470 
471  call get_q2_level2( &
472  ka, ks, ke_pbl, &
473  dudz2(:,i,j), ri(:,i,j), l(:,i,j), & ! (in)
474  q2_2(:) ) ! (out)
475 
476  do k = ks, ke_pbl
477  ac(k) = min( q(k) / sqrt(q2_2(k)), 1.0_rp )
478  end do
479 
480  call get_smsh( &
481  ka, ks, ke_pbl, & ! (in)
482  q(:), ac(:), & ! (in)
483  l(:,i,j), n2_new(:), & ! (in)
484  pott(:,i,j), dudz2(:,i,j), & ! (in)
485  tvsq(:), tvsq25(:), & ! (in) ! dummy
486  .false., & ! (in)
487  sm(:), sh(:) ) ! (out)
488 
489 
490  ! partial condensation
491 
492  do k = ks, ke_pbl+1
493  teml(k) = potl(k,i,j) * exner(k,i,j)
494  end do
495 
496  call atmos_saturation_dens2qsat( &
497  ka, ks, ke_pbl, &
498  teml(:), dens(:,i,j), & ! (in)
499  qsl(:) ) ! (out)
500 
501  call hydrometeor_lhv( &
502  ka, ks, ke_pbl, & ! (in)
503  teml(:), & ! (in)
504  lhvl(:) ) ! (out)
505 
506  do k = ks, ke_pbl
507 
508  dqsl = qsl(k) * lhvl(k) / ( rvap * teml(k)**2 )
509  cptot = qdry(k,i,j) * cpdry + qsl(k) * cpvap
510  aa = 1.0_rp / ( 1.0_rp + lhvl(k)/cptot * dqsl )
511  bb = exner(k,i,j) * dqsl
512 
513  if ( mynn_level3 ) then
514  ! level 3
515  tsq(k) = max( prog(k,i,j,i_tsq), 0.0_rp )
516  qsq(k) = max( prog(k,i,j,i_qsq), 0.0_rp )
517  cov(k) = prog(k,i,j,i_cov)
518  cov(k) = sign( min( abs(cov(k)), sqrt(tsq(k) * qsq(k))), cov(k) )
519  sigma_s = 0.5_rp * aa * sqrt( max( qsq(k) - 2.0_rp * bb * cov(k) + bb**2 * tsq(k), 1.0e-20_rp ) )
520  else
521  ! level 2.5
522  sigma_s = max( 0.5_rp * aa * l(k,i,j) * sqrt( ac(k) * b2 * sh(k) ) * abs( dqwdz(k) - bb * dtldz(k) ), &
523  1.0e-10_rp )
524  end if
525 
526  q1 = aa * ( qw(k,i,j) - qsl(k) ) * 0.5_rp / sigma_s
527  rr = min( max( 0.5_rp * ( 1.0_rp + erf(q1*rsqrt_2) ), 0.0_rp ), 1.0_rp )
528  qlp = min( max( 2.0_rp * sigma_s * ( rr * q1 + rsqrt_2pi &
529 #if defined(PGI) || defined(SX)
530  * exp( -min( 0.5_rp*q1**2, 1.e+3_rp ) ) & ! apply exp limiter
531 #else
532  * exp(-0.5_rp*q1**2) &
533 #endif
534  ), 0.0_rp ), qw(k,i,j) * 0.5_rp )
535  cc = ( 1.0_rp + epstvap * qw(k,i,j) - (1.0_rp+epstvap) * qlp ) / exner(k,i,j) * lhvl(k) / cptot &
536  - (1.0_rp+epstvap) * pott(k,i,j)
537  rt = min( max( rr - qlp / (2.0_rp*sigma_s*sqrt_2pi) &
538 #if defined(PGI) || defined(SX)
539  * exp( -min( 0.5_rp*q1**2, 1.e+3_rp ) ) & ! apply exp limiter
540 #else
541  * exp(-q1**2 * 0.5_rp) &
542 #endif
543  , 0.0_rp ), 1.0_rp )
544  betat = 1.0_rp + epstvap * qw(k,i,j) - (1.0_rp+epstvap) * qlp - rt * aa * bb * cc
545  betaq = epstvap * pott(k,i,j) + rt * aa * cc
546 
547  if ( mynn_level3 ) then
548  if ( k==ks ) then
549  ! production at KS
550  ! us3 = - l_mo(i,j) * KARMAN * GRAV * SFLX_PT / POTT(KS,i,j) ! u_*^3
551  us = max(- l_mo(i,j) * karman * grav * sflx_pt / pott(ks,i,j), 0.0_rp)**(1.0_rp/3.0_rp)
552  zeta = z1 / l_mo(i,j)
553  ! Businger et al. (1971)
554 !!$ if ( zeta > 0 ) then
555 !!$ phi_h = 4.7_RP * z1 / l_mo(i,j) / 0.74_RP + 1.0_RP
556 !!$ else
557 !!$ phi_h = 1.0_RP / sqrt( 1.0_RP - 9.0_RP * z1 / l_mo(i,j) )
558 !!$ end if
559  ! Beljaars and Holtslag (1991)
560  if ( zeta > 0 ) then
561  phi_h = - 2.0_rp / 3.0_rp * ( 0.35_rp * zeta - 6.0_rp ) * exp(-0.35_rp*zeta) + zeta * sqrt( 1.0_rp + 2.0_rp * zeta / 3.0_rp ) + 1.0_rp
562  else
563  phi_h = 1.0_rp / sqrt( 1.0_rp - 16.0_rp * zeta )
564  end if
565  ! TSQ
566  prod_t1 = 1.0_rp / us * phi_h / ( karman * z1 ) * sflx_pt**2
567  ! QSQ
568  prod_q1 = 1.0_rp / us * phi_h / ( karman * z1 ) * ( sflx_qv(i,j) / dens(ks,i,j) )**2
569  ! COV
570  prod_c1 = 1.0_rp / us * phi_h * ( karman * z1 ) * sflx_pt * sflx_qv(i,j) / dens(ks,i,j)
571  tsq25 = prod_t1 * b2 * l(k,i,j) / q(k) * 0.5_rp
572  qsq25 = prod_q1 * b2 * l(k,i,j) / q(k) * 0.5_rp
573  cov25 = prod_c1 * b2 * l(k,i,j) / q(k) * 0.5_rp
574  else
575  tsq25 = b2 * l(k,i,j)**2 * sh(k) * dtldz(k)**2
576  qsq25 = b2 * l(k,i,j)**2 * sh(k) * dqwdz(k)**2
577  cov25 = b2 * l(k,i,j)**2 * sh(k) * dtldz(k) * dqwdz(k)
578  end if
579  tltv = betat * tsq25 + betaq * cov25
580  qwtv = betat * cov25 + betaq * qsq25
581  tvsq25(k) = max(betat * tltv + betaq * qwtv, 0.0_rp)
582  ! level 3
583  if ( tsq(k) == 0.0_rp ) then
584  ! teq, qsq is not initialized
585  tsq(k) = tsq25
586  qsq(k) = qsq25
587  cov(k) = cov25
588  end if
589  tltv = betat * tsq(k) + betaq * cov(k)
590  qwtv = betat * cov(k) + betaq * qsq(k)
591  tvsq(k) = max(betat * tltv + betaq * qwtv, 0.0_rp)
592  end if
593 
594  n2_new(k) = min(atmos_phy_bl_mynn_n2_max, &
595  grav * ( dtldz(k) * betat + dqwdz(k) * betaq ) / potv(k,i,j) )
596  end do
597 
598  do k = ks, ke_pbl
599  ri(k,i,j) = n2_new(k) / max(dudz2(k,i,j), 1e-10_rp)
600  end do
601  do k = ke_pbl+1, ke
602  ri(k,i,j) = 0.0_rp
603  n2_new(k ) = 0.0_rp
604  end do
605 
606 
607  ! length
608  call get_length( &
609  ka, ks, ke_pbl, &
610  pott(ks,i,j), q(:), n2_new(:), & ! (in)
611  sflx_pt, l_mo(i,j), & ! (in)
612  fz(:,i,j), & ! (in)
613  l(:,i,j) ) ! (out)
614 
615  call get_q2_level2( &
616  ka, ks, ke_pbl, &
617  dudz2(:,i,j), ri(:,i,j), l(:,i,j), & ! (in)
618  q2_2(:) ) ! (out)
619 
620  do k = ks, ke_pbl
621  ac(k) = min( q(k) / sqrt(q2_2(k)), 1.0_rp )
622  end do
623 
624  call get_smsh( &
625  ka, ks, ke_pbl, & ! (in)
626  q(:), ac(:), & ! (in)
627  l(:,i,j), n2_new(:), & ! (in)
628  pott(:,i,j), dudz2(:,i,j), & ! (in)
629  tvsq(:), tvsq25(:), & ! (in)
630  mynn_level3, & ! (in)
631  sm(:), sh(:) ) ! (out)
632 
633 
634  do k = ks, ke_pbl
635  nu(k,i,j) = max( min( l(k,i,j) * q(k) * sm(k), &
636  atmos_phy_bl_mynn_nu_max ), &
637  atmos_phy_bl_mynn_nu_min )
638  kh(k,i,j) = max( min( l(k,i,j) * q(k) * sh(k), &
639  atmos_phy_bl_mynn_kh_max ), &
640  atmos_phy_bl_mynn_kh_min )
641  pr(k,i,j) = nu(k,i,j) / kh(k,i,j)
642  end do
643  do k = ke_pbl+1, ke
644  nu(k,i,j) = 0.0_rp
645  kh(k,i,j) = 0.0_rp
646  pr(k,i,j) = 1.0_rp
647  end do
648 
649  ! dens * coefficient at the half level
650  do k = ks, ke_pbl-1
651  rhonu(k) = f2h(k,1) * dens(k+1,i,j) * nu(k+1,i,j) &
652  + f2h(k,2) * dens(k ,i,j) * nu(k ,i,j)
653  rhokh(k) = f2h(k,1) * dens(k+1,i,j) * kh(k+1,i,j) &
654  + f2h(k,2) * dens(k ,i,j) * kh(k ,i,j)
655  end do
656 
657 
658 
659  ! time integration
660 
661  ! dens * u
662  sf_t = sflx_mu(i,j) / ( fz(ks,i,j) - fz(ks-1,i,j) )
663  d(ks) = u(ks,i,j) + dt * sf_t / dens(ks,i,j)
664  do k = ks+1, ke_pbl
665  d(k) = u(k,i,j)
666  end do
667  c(ks) = 0.0_rp
668  do k = ks, ke_pbl-1
669  ap = - dt * rhonu(k) / ( cz(k+1,i,j) - cz(k,i,j) )
670  a(k) = ap / ( dens(k,i,j) * ( fz(k,i,j) - fz(k-1,i,j) ) )
671  b(k) = - a(k) - c(k) + 1.0_rp
672  c(k+1) = ap / ( dens(k+1,i,j) * ( fz(k+1,i,j) - fz(k,i,j) ) )
673  end do
674  a(ke_pbl) = 0.0_rp
675  b(ke_pbl) = - c(ke_pbl) + 1.0_rp
676 
677  call matrix_solver_tridiagonal( &
678  ka, ks, ke_pbl, &
679  a(:), b(:), c(:), d(:), & ! (in)
680  dummy(:) ) ! (out)
681 ! phi_n(:) ) ! (out)
682 
683  phi_n(:) = dummy(:)
684  rhou_t(ks,i,j) = ( phi_n(ks) - u(ks,i,j) ) * dens(ks,i,j) / dt - sf_t
685  do k = ks+1, ke_pbl
686  rhou_t(k,i,j) = ( phi_n(k) - u(k,i,j) ) * dens(k,i,j) / dt
687  end do
688  do k = ke_pbl+1, ke
689  rhou_t(k,i,j) = 0.0_rp
690  end do
691  do k = ks, ke_pbl-1
692  flxu(k,i,j) = - rhonu(k) * ( phi_n(k+1) - phi_n(k) ) / ( cz(k+1,i,j) - cz(k,i,j) )
693  end do
694  do k = ke_pbl, ke
695  flxu(k,i,j) = 0.0_rp
696  end do
697 
698 
699  ! dens * v
700  sf_t = sflx_mv(i,j) / ( fz(ks,i,j) - fz(ks-1,i,j) )
701  d(ks) = v(ks,i,j) + dt * sf_t / dens(ks,i,j)
702  do k = ks+1, ke_pbl
703  d(k) = v(k,i,j)
704  end do
705  ! a,b,c is the same as those for the u
706 
707  call matrix_solver_tridiagonal( &
708  ka, ks, ke_pbl, &
709  a(:), b(:), c(:), d(:), & ! (in)
710  phi_n(:) ) ! (out)
711 
712  rhov_t(ks,i,j) = ( phi_n(ks) - v(ks,i,j) ) * dens(ks,i,j) / dt - sf_t
713  do k = ks+1, ke_pbl
714  rhov_t(k,i,j) = ( phi_n(k) - v(k,i,j) ) * dens(k,i,j) / dt
715  end do
716  do k = ke_pbl+1, ke
717  rhov_t(k,i,j) = 0.0_rp
718  end do
719  do k = ks, ke_pbl-1
720  flxv(k,i,j) = - rhonu(k) * ( phi_n(k+1) - phi_n(k) ) / ( cz(k+1,i,j) - cz(k,i,j) )
721  end do
722  do k = ke_pbl, ke
723  flxv(k,i,j) = 0.0_rp
724  end do
725 
726 
727  ! dens * pott
728  sf_t = sflx_pt / ( fz(ks,i,j) - fz(ks-1,i,j) )
729  d(ks) = pott(ks,i,j) + dt * sf_t
730  do k = ks+1, ke_pbl
731  d(k) = pott(k,i,j)
732  end do
733 
734  c(ks) = 0.0_rp
735  do k = ks, ke_pbl-1
736  ap = - dt * rhokh(k) / ( cz(k+1,i,j) - cz(k,i,j) )
737  a(k) = ap / ( dens(k,i,j) * ( fz(k,i,j) - fz(k-1,i,j) ) )
738  b(k) = - a(k) - c(k) + 1.0_rp
739  c(k+1) = ap / ( dens(k+1,i,j) * ( fz(k+1,i,j) - fz(k,i,j) ) )
740  end do
741  a(ke_pbl) = 0.0_rp
742  b(ke_pbl) = - c(ke_pbl) + 1.0_rp
743 
744  call matrix_solver_tridiagonal( &
745  ka, ks, ke_pbl, &
746  a(:), b(:), c(:), d(:), & ! (in)
747  phi_n(:) ) ! (out)
748 
749  rhot_t(ks,i,j) = ( ( phi_n(ks) - pott(ks,i,j) ) / dt - sf_t ) * dens(ks,i,j)
750  do k = ks+1, ke_pbl
751  rhot_t(k,i,j) = ( phi_n(k) - pott(k,i,j) ) * dens(k,i,j) / dt
752  end do
753  do k = ke_pbl+1, ke
754  rhot_t(k,i,j) = 0.0_rp
755  end do
756  do k = ks, ke_pbl-1
757  flxt(k,i,j) = - rhokh(k) * ( phi_n(k+1) - phi_n(k) ) / ( cz(k+1,i,j) - cz(k,i,j) )
758  end do
759  do k = ke_pbl, ke
760  flxt(k,i,j) = 0.0_rp
761  end do
762 
763 
764  ! dens * TKE
765  ! production at KS: us3 * phi_m(zeta) / ( KARMAN * z1 )
766  us3 = - l_mo(i,j) * karman * grav * sflx_pt / pott(ks,i,j) ! u_*^3
767  zeta = z1 / l_mo(i,j)
768 !!$ ! Businger et al. (1971)
769 !!$ if ( zeta > 0 ) then
770 !!$ phi_m = 4.7_RP * zeta + 1.0_RP
771 !!$ else
772 !!$ phi_m = 1.0_RP / sqrt(sqrt( 1.0_RP - 15.0_RP * zeta ))
773 !!$ end if
774  ! Beljaars and Holtslag (1991)
775  if ( zeta > 0 ) then
776  phi_m = - 2.0_rp / 3.0_rp * ( 0.35_rp * zeta - 6.0_rp ) * zeta * exp(-0.35_rp*zeta) + zeta + 1.0_rp
777  else
778  phi_m = 1.0_rp / sqrt(sqrt(1.0_rp - 16.0_rp * zeta))
779  end if
780  prod(ks,i,j) = us3 / ( karman * z1 ) * ( phi_m - zeta )
781  do k = ks+1, ke_pbl
782 ! do k = KS, KE_PBL
783  prod(k,i,j) = nu(k,i,j) * dudz2(k,i,j) - kh(k,i,j) * n2_new(k)
784  end do
785  do k = ks, ke_pbl
786  tke_p = q(k)**2 * 0.5_rp
787  prod(k,i,j) = max( prod(k,i,j), - tke_p / dt)
788  d(k) = tke_p + dt * prod(k,i,j)
789  end do
790  do k = ke_pbl+1, ke
791  prod(k,i,j) = 0.0_rp
792  end do
793  c(ks) = 0.0_rp
794  do k = ks, ke_pbl-1
795  ap = - dt * 3.0_rp * rhonu(k) / ( cz(k+1,i,j) - cz(k,i,j) )
796  a(k) = ap / ( dens(k,i,j) * ( fz(k,i,j) - fz(k-1,i,j) ) )
797  diss(k,i,j) = 2.0_rp * q(k) / ( b1 * l(k,i,j) )
798  b(k) = - a(k) - c(k) + 1.0_rp + diss(k,i,j) * dt
799  c(k+1) = ap / ( dens(k+1,i,j) * ( fz(k+1,i,j) - fz(k,i,j) ) )
800  end do
801  a(ke_pbl) = 0.0_rp
802  diss(ke_pbl,i,j) = 2.0_rp * q(ke_pbl) / ( b1 * l(ke_pbl,i,j) )
803  b(ke_pbl) = - c(ke_pbl) + 1.0_rp + diss(ke_pbl,i,j) * dt
804  do k = ke_pbl+1, ke
805  diss(k,i,j) = 0.0_rp
806  end do
807 
808  call matrix_solver_tridiagonal( &
809  ka, ks, ke_pbl, &
810  a(:), b(:), c(:), d(:), & ! (in)
811  phi_n(:) ) ! (out)
812 
813  do k = ks, ke_pbl
814  diss(k,i,j) = diss(k,i,j) * phi_n(k)
815  rprog_t(k,i,j,i_tke) = ( max(phi_n(k), atmos_phy_bl_mynn_tke_min) - prog(k,i,j,i_tke) ) * dens(k,i,j) / dt
816  end do
817  do k = ke_pbl+1, ke
818  rprog_t(k,i,j,i_tke) = 0.0_rp
819  end do
820 
821 
822 
823  if ( .not. mynn_level3 ) cycle
824 
825 
826  ! dens * tsq
827  d(ks) = max(tsq(ks) + dt * prod_t1, 0.0_rp)
828  do k = ks+1, ke_pbl
829 ! do k = KS, KE_PBL
830  d(k) = max(tsq(k) + dt * 2.0_rp * l(k,i,j) * q(k) * sh(k) * dtldz(k)**2, 0.0_rp)
831  end do
832  c(ks) = 0.0_rp
833  do k = ks, ke_pbl-1
834  ap = - dt * rhonu(k) / ( cz(k+1,i,j) - cz(k,i,j) )
835  a(k) = ap / ( dens(k,i,j) * ( fz(k,i,j) - fz(k-1,i,j) ) )
836  b(k) = - a(k) - c(k) + 1.0_rp + dt * 2.0_rp * q(k) / ( b2 * l(k,i,j) )
837  c(k+1) = ap / ( dens(k+1,i,j) * ( fz(k+1,i,j) - fz(k,i,j) ) )
838  end do
839  a(ke_pbl) = 0.0_rp
840  b(ke_pbl) = - c(ke_pbl) + 1.0_rp + dt * 2.0_rp * q(ke_pbl) / ( b2 * l(ke_pbl,i,j) )
841 
842  call matrix_solver_tridiagonal( &
843  ka, ks, ke_pbl, &
844  a(:), b(:), c(:), d(:), & ! (in)
845  tsq(:) ) ! (out)
846 
847  do k = ks, ke_pbl
848  rprog_t(k,i,j,i_tsq) = ( max(tsq(k), 0.0_rp) - prog(k,i,j,i_tsq) ) * dens(k,i,j) / dt
849  end do
850  do k = ke_pbl+1, ke
851  rprog_t(k,i,j,i_tsq) = 0.0_rp
852  end do
853 
854 
855  ! dens * qsq
856  d(ks) = max(qsq(ks) + dt * prod_q1, 0.0_rp)
857  do k = ks+1, ke_pbl
858 ! do k = KS, KE_PBL
859  d(k) = max(qsq(k) + dt * 2.0_rp * l(k,i,j) * q(k) * sh(k) * dqwdz(k)**2, 0.0_rp)
860  end do
861  ! a, b, c are same as those for tsq
862 
863  call matrix_solver_tridiagonal( &
864  ka, ks, ke_pbl, &
865  a(:), b(:), c(:), d(:), & ! (in)
866  qsq(:) ) ! (out)
867 
868  do k = ks, ke_pbl
869  rprog_t(k,i,j,i_qsq) = ( max(qsq(k), 0.0_rp) - prog(k,i,j,i_qsq) ) * dens(k,i,j) / dt
870  end do
871  do k = ke_pbl+1, ke
872  rprog_t(k,i,j,i_qsq) = 0.0_rp
873  end do
874 
875 
876  ! dens * cov
877  d(ks) = cov(ks) + dt * prod_c1
878  do k = ks+1, ke_pbl
879 ! do k = KS, KE_PBL
880  d(k) = cov(k) + dt * 2.0_rp * l(k,i,j) * q(k) * sh(k) * dtldz(k) * dqwdz(k)
881  end do
882  ! a, b, c are same as those for tsq
883 
884  call matrix_solver_tridiagonal( &
885  ka, ks, ke_pbl, &
886  a(:), b(:), c(:), d(:), & ! (in)
887  cov(:) ) ! (out)
888 
889  do k = ks, ke_pbl
890  rprog_t(k,i,j,i_cov) = ( cov(k) - prog(k,i,j,i_cov) ) * dens(k,i,j) / dt
891  end do
892  do k = ke_pbl+1, ke
893  rprog_t(k,i,j,i_cov) = 0.0_rp
894  end do
895 
896 
897  end do
898  end do
899 
900 
901  l(ke_pbl+1:ke,:,:) = 0.0_rp
902  call file_history_in(ri(:,:,:), 'Ri_MYNN', 'Richardson number', '1', fill_halo=.true. )
903  call file_history_in(pr(:,:,:), 'Pr_MYNN', 'Prandtl number', '1', fill_halo=.true. )
904  call file_history_in(prod(:,:,:), 'TKE_prod_MYNN', 'TKE production', 'm2/s3', fill_halo=.true.)
905  call file_history_in(diss(:,:,:), 'TKE_diss_MYNN', 'TKE dissipation', 'm2/s3', fill_halo=.true.)
906  call file_history_in(dudz2(:,:,:), 'dUdZ2_MYNN', 'dudz2', 'm2/s2', fill_halo=.true.)
907  call file_history_in(l(:,:,:), 'L_mix_MYNN', 'minxing length', 'm', fill_halo=.true.)
908 
909  call file_history_in(flxu(:,:,:), 'ZFLX_RHOU_MYNN', 'Z FLUX of RHOU (MYNN)', 'kg/m/s2', fill_halo=.true.)
910  call file_history_in(flxv(:,:,:), 'ZFLX_RHOV_MYNN', 'Z FLUX of RHOV (MYNN)', 'kg/m/s2', fill_halo=.true.)
911  call file_history_in(flxt(:,:,:), 'ZFLX_RHOT_MYNN', 'Z FLUX of RHOT (MYNN)', 'K kg/m2/s', fill_halo=.true.)
912 
913  return
914  end subroutine atmos_phy_bl_mynn_tendency
915 
916  !-----------------------------------------------------------------------------
920  subroutine atmos_phy_bl_mynn_tendency_tracer( &
921  KA, KS, KE, IA, IS, IE, JA, JS, JE, &
922  DENS, QTRC, SFLX_Q, Kh, &
923  CZ, FZ, DT, TRACER_NAME, &
924  RHOQ_t )
925  use scale_matrix, only: &
926  matrix_solver_tridiagonal
927  use scale_file_history, only: &
928  file_history_in
929  integer, intent(in) :: ka, ks, ke
930  integer, intent(in) :: ia, is, ie
931  integer, intent(in) :: ja, js, je
932 
933  real(RP), intent(in) :: dens (ka,ia,ja)
934  real(RP), intent(in) :: qtrc (ka,ia,ja)
935  real(RP), intent(in) :: sflx_q( ia,ja)
936  real(RP), intent(in) :: kh (ka,ia,ja)
937  real(RP), intent(in) :: cz ( ka,ia,ja)
938  real(RP), intent(in) :: fz (0:ka,ia,ja)
939  real(DP), intent(in) :: dt
940  character(len=*), intent(in) :: tracer_name
941 
942  real(RP), intent(out) :: rhoq_t(ka,ia,ja)
943 
944  real(RP) :: qtrc_n(ka)
945  real(RP) :: rhokh(ka)
946  real(RP) :: a(ka)
947  real(RP) :: b(ka)
948  real(RP) :: c(ka)
949  real(RP) :: d(ka)
950  real(RP) :: ap
951  real(RP) :: sf_t
952 
953  real(RP) :: flx(ka,ia,ja)
954 
955  real(RP) :: f2h(ka,2)
956 
957  integer :: k, i, j
958 
959 !OCL INDEPENDENT
960  !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
961  !$omp shared(KA,KS,KE_PBL,KE,IS,IE,JS,JE) &
962  !$omp shared(RHOQ_t,DENS,QTRC,SFLX_Q,Kh,CZ,FZ,DT,flx) &
963  !$omp private(QTRC_n,RHOKh,a,b,c,d,ap,sf_t,f2h) &
964  !$omp private(k,i,j)
965  do j = js, je
966  do i = is, ie
967 
968  call get_f2h( &
969  ka, ks, ke, &
970  fz(:,i,j), & ! (in)
971  f2h(:,:) ) ! (out)
972 
973  ! dens * coefficient at the half level
974  do k = ks, ke_pbl-1
975  rhokh(k) = f2h(k,1) * dens(k+1,i,j) * kh(k+1,i,j) &
976  + f2h(k,2) * dens(k ,i,j) * kh(k ,i,j)
977  end do
978 
979  sf_t = sflx_q(i,j) / ( fz(ks,i,j) - fz(ks-1,i,j) )
980  d(ks) = qtrc(ks,i,j) + dt * sf_t / dens(ks,i,j)
981  do k = ks+1, ke_pbl
982  d(k) = qtrc(k,i,j)
983  end do
984 
985  c(ks) = 0.0_rp
986  do k = ks, ke_pbl-1
987  ap = - dt * rhokh(k) / ( cz(k+1,i,j) - cz(k,i,j) )
988  a(k) = ap / ( dens(k,i,j) * ( fz(k,i,j) - fz(k-1,i,j) ) )
989  b(k) = - a(k) - c(k) + 1.0_rp
990  c(k+1) = ap / ( dens(k+1,i,j) * ( fz(k+1,i,j) - fz(k,i,j) ) )
991  end do
992  a(ke_pbl) = 0.0_rp
993  b(ke_pbl) = - c(ke_pbl) + 1.0_rp
994 
995  call matrix_solver_tridiagonal( &
996  ka, ks, ke_pbl, &
997  a(:), b(:), c(:), d(:), & ! (in)
998  qtrc_n(:) ) ! (out)
999 
1000  rhoq_t(ks,i,j) = ( qtrc_n(ks) - qtrc(ks,i,j) ) * dens(ks,i,j) / dt - sf_t
1001  do k = ks+1, ke_pbl
1002  rhoq_t(k,i,j) = ( qtrc_n(k) - qtrc(k,i,j) ) * dens(k,i,j) / dt
1003  end do
1004  do k = ke_pbl+1, ke
1005  rhoq_t(k,i,j) = 0.0_rp
1006  end do
1007 
1008  do k = ks, ke_pbl-1
1009  flx(k,i,j) = - rhokh(k) * ( qtrc_n(k+1) - qtrc_n(k) ) / ( cz(k+1,i,j) - cz(k,i,j) )
1010  end do
1011  do k = ke_pbl, ke
1012  flx(k,i,j) = 0.0_rp
1013  end do
1014 
1015  end do
1016  end do
1017 
1018  call file_history_in(flx(:,:,:), 'ZFLX_'//trim(tracer_name)//'_MYNN', 'Z FLUX of DENS * '//trim(tracer_name)//' (MYNN)', 'kg/m2/s', fill_halo=.true.)
1019 
1020  return
1021  end subroutine atmos_phy_bl_mynn_tendency_tracer
1022 
1023  !-----------------------------------------------------------------------------
1024  ! private routines
1025  !-----------------------------------------------------------------------------
1026 
1027 !OCL SERIAL
1028  subroutine get_length( &
1029  KA, KS, KE_PBL, &
1030  PT0, q, n2, &
1031  SFLX_PT, l_mo, &
1032  FZ, &
1033  l )
1034  use scale_const, only: &
1035  grav => const_grav, &
1036  karman => const_karman, &
1037  cp => const_cpdry, &
1038  eps => const_eps
1039  implicit none
1040  integer, intent(in) :: ka, ks, ke_pbl
1041 
1042  real(RP), intent(in) :: pt0
1043  real(RP), intent(in) :: q(ka)
1044  real(RP), intent(in) :: n2(ka)
1045  real(RP), intent(in) :: sflx_pt
1046  real(RP), intent(in) :: l_mo
1047  real(RP), intent(in) :: fz(0:ka)
1048 
1049  real(RP), intent(out) :: l(ka)
1050 
1051  real(RP) :: ls
1052  real(RP) :: lt
1053  real(RP) :: lb
1054  real(RP) :: rlm
1055  real(RP) :: rlt
1056 
1057  real(RP) :: qc
1058  real(RP) :: int_q
1059  real(RP) :: int_qz
1060  real(RP) :: rn2sr
1061  real(RP) :: zeta
1062 
1063  real(RP) :: z
1064  real(RP) :: qdz
1065 
1066  real(RP) :: sw
1067  integer :: k
1068 
1069  int_qz = 0.0_rp
1070  int_q = 0.0_rp
1071  do k = ks, ke_pbl
1072  qdz = q(k) * ( fz(k) - fz(k-1) )
1073  int_qz = int_qz + ((fz(k)+fz(k-1))*0.5_rp-fz(ks-1)) * qdz
1074  int_q = int_q + qdz
1075  end do
1076  ! LT
1077  lt = min( max(0.23_rp * int_qz / (int_q + eps), &
1078  lt_min), &
1079  atmos_phy_bl_mynn_lt_max )
1080  rlt = 1.0_rp / lt
1081 
1082  rlm = 1.0_rp / l_mo
1083 
1084  qc = ( grav / pt0 * max(sflx_pt,0.0_rp) * lt )**oneoverthree
1085 
1086  do k = ks, ke_pbl
1087  z = ( fz(k)+fz(k-1) )*0.5_rp - fz(ks-1)
1088  zeta = z * rlm
1089 
1090  ! LS
1091  sw = sign(0.5_rp, zeta) + 0.5_rp ! 1 for zeta >= 0, 0 for zeta < 0
1092  ls = karman * z &
1093  * ( sw / (1.0_rp + 2.7_rp*min(zeta,1.0_rp)*sw ) &
1094  + ( (1.0_rp - 100.0_rp*zeta)*(1.0_rp-sw) )**0.2_rp )
1095 
1096  ! LB
1097  sw = sign(0.5_rp, n2(k)-eps) + 0.5_rp ! 1 for dptdz >0, 0 for dptdz <= 0
1098  rn2sr = 1.0_rp / ( sqrt(n2(k)*sw) + 1.0_rp-sw)
1099  lb = (1.0_rp + 5.0_rp * sqrt(qc*rn2sr/lt)) * q(k) * rn2sr * sw & ! qc=0 when l_mo > 0
1100  + 999.e10_rp * (1.0_rp-sw)
1101 
1102  ! L
1103  l(k) = 1.0_rp / ( 1.0_rp/ls + rlt + 1.0_rp/lb )
1104  end do
1105 
1106  return
1107  end subroutine get_length
1108 
1109 !OCL SERIAL
1110  subroutine get_q2_level2( &
1111  KA, KS, KE_PBL, &
1112  dudz2, Ri, l, &
1113  q2_2 )
1114  implicit none
1115  integer, intent(in) :: ka, ks, ke_pbl
1116 
1117  real(RP), intent(in) :: dudz2(ka)
1118  real(RP), intent(in) :: ri(ka)
1119  real(RP), intent(in) :: l(ka)
1120 
1121  real(RP), intent(out) :: q2_2(ka)
1122 
1123  real(RP) :: rf
1124  real(RP) :: sm_2
1125  real(RP) :: sh_2
1126  real(RP) :: q2
1127 
1128  integer :: k
1129 
1130  do k = ks, ke_pbl
1131  rf = min(0.5_rp / af12 * ( ri(k) &
1132  + af12*rf1 &
1133  - sqrt(ri(k)**2 + 2.0_rp*af12*(rf1-2.0_rp*rf2)*ri(k) + (af12*rf1)**2) ), &
1134  rfc)
1135  sh_2 = 3.0_rp * a2 * (g1+g2) * (rfc-rf) / (1.0_rp-rf)
1136  sm_2 = sh_2 * af12 * (rf1-rf) / (rf2-rf)
1137  q2 = b1 * l(k)**2 * sm_2 * (1.0_rp-rf) * dudz2(k)
1138  q2_2(k) = max( q2, 1.e-10_rp )
1139  end do
1140 
1141  return
1142  end subroutine get_q2_level2
1143 
1144 !OCL SERIAL
1145  subroutine get_smsh( &
1146  KA, KS, KE_PBL, &
1147  q, ac, &
1148  l, n2, &
1149  pott, dudz2, &
1150  tvsq, tvsq25, &
1151  mynn_level3, &
1152  sm, sh )
1153  use scale_const, only: &
1154  eps => const_eps, &
1155  grav => const_grav
1156  implicit none
1157  integer, intent(in) :: ka, ks, ke_pbl
1158 
1159  real(RP), intent(in) :: q(ka)
1160  real(RP), intent(in) :: ac(ka)
1161  real(RP), intent(in) :: l(ka)
1162  real(RP), intent(in) :: n2(ka)
1163  real(RP), intent(in) :: pott(ka)
1164  real(RP), intent(in) :: dudz2(ka)
1165  real(RP), intent(in) :: tvsq(ka)
1166  real(RP), intent(in) :: tvsq25(ka)
1167  logical, intent(in) :: mynn_level3
1168 
1169  real(RP), intent(out) :: sm (ka) ! S_M2.5 + S_M'
1170  real(RP), intent(out) :: sh (ka) ! S_H2.5 + S_H'
1171 
1172 
1173  real(RP) :: l2q2
1174  real(RP) :: ac2
1175  real(RP) :: p1
1176  real(RP) :: p2
1177  real(RP) :: p3
1178  real(RP) :: p4
1179  real(RP) :: p5
1180  real(RP) :: rd25
1181  real(RP) :: gh
1182 
1183  ! for level 3
1184  real(RP) :: em
1185  real(RP) :: eh
1186  real(RP) :: ew
1187  real(RP) :: rdp
1188  real(RP) :: cw25
1189  real(RP) :: fact
1190 
1191  integer :: k
1192 
1193  ! level 2.5
1194  do k = ks, ke_pbl
1195 
1196  ac2 = ac(k)**2
1197  l2q2 = ( l(k) / max(q(k),eps) )**2
1198 
1199  gh = - n2(k) * l2q2
1200 
1201  p1 = 1.0_rp - 3.0_rp * ac2 * a2 * b2 * ( 1.0_rp - c3 ) * gh
1202  p2 = 1.0_rp - 9.0_rp * ac2 * a1 * a2 * ( 1.0_rp - c2 ) * gh
1203  p3 = p1 + 9.0_rp * ac2 * a2**2 * ( 1.0_rp - c2 ) * ( 1.0_rp - c5 ) * gh
1204  p4 = p1 - 12.0_rp * ac2 * a1 * a2 * ( 1.0_rp - c2 ) * gh
1205  p5 = 6.0_rp * ac2 * a1**2 * dudz2(k) * l2q2
1206 
1207  rd25 = 1.0_rp / max( p2 * p4 + p5 * p3, eps )
1208 
1209  sm(k) = max( ac(k) * a1 * ( p3 - 3.0_rp * c1 * p4 ) * rd25, 0.0_rp )
1210  sh(k) = max( ac(k) * a2 * ( p2 + 3.0_rp * c1 * p5 ) * rd25, 0.0_rp )
1211 
1212  end do
1213 
1214  ! level 3
1215  if ( mynn_level3 ) then
1216 
1217  do k = ks, ke_pbl
1218 
1219  ac2 = ac(k)**2
1220  l2q2 = ( l(k) / max(q(k),eps) )**2
1221  ! restriction: L/q <= 1/N
1222  if ( n2(k) > 0 ) l2q2 = min( l2q2, 1.0_rp/n2(k) )
1223  !l2q2 = min( l2q2, 1.0_RP/abs(n2(k)) )
1224 
1225  gh = - n2(k) * l2q2
1226  if ( abs(gh) < eps ) gh = eps
1227 
1228  p1 = 1.0_rp - 3.0_rp * ac2 * a2 * b2 * ( 1.0_rp - c3 ) * gh
1229  p2 = 1.0_rp - 9.0_rp * ac2 * a1 * a2 * ( 1.0_rp - c2 ) * gh
1230  p3 = p1 + 9.0_rp * ac2 * a2**2 * ( 1.0_rp - c2 ) * ( 1.0_rp - c5 ) * gh
1231  p4 = p1 - 12.0_rp * ac2 * a1 * a2 * ( 1.0_rp - c2 ) * gh
1232  p5 = 6.0_rp * ac2 * a1**2 * dudz2(k) * l2q2
1233 
1234  rd25 = 1.0_rp / max( p2 * p4 + p5 * p3, eps )
1235  rdp = 1.0_rp / max( p2 * ( p4 - p1 + 1.0_rp ) + p5 * ( p3 - p1 + 1.0_rp ), eps )
1236 
1237  cw25 = p1 * ( p2 + 3.0_rp * c1 * p5 ) * rd25 / 3.0_rp
1238 
1239  ew = ( 1.0_rp - c3 ) * ( p2 * ( p1 - p4 ) + p5 * ( p1 - p3 ) ) * rdp
1240  ew = sign( max(abs(ew),eps), ew )
1241  fact = ( l2q2 * grav / ( l(k) * pott(k) ) )**2 * ( tvsq(k) - tvsq25(k) ) / gh
1242  fact = fact * ew
1243  fact = min( max( fact, 0.12_rp - cw25 ), 0.76_rp - cw25 )
1244  fact = fact / ew
1245 
1246  em = 3.0_rp * ac(k) * a1 * ( 1.0_rp - c3 ) * ( p3 - p4 ) * rdp
1247  eh = 3.0_rp * ac(k) * a2 * ( 1.0_rp - c3 ) * ( p2 + p5 ) * rdp
1248 
1249 
1250  sm(k) = sm(k) + em * fact
1251  sh(k) = sh(k) + eh * fact
1252 
1253  end do
1254 
1255  end if
1256 
1257  return
1258  end subroutine get_smsh
1259 
1260 !OCL SERIAL
1261  subroutine get_f2h( &
1262  KA, KS, KE, &
1263  FZ, &
1264  f2h )
1265  integer, intent(in) :: ka, ks, ke
1266  real(RP), intent(in) :: fz(0:ka)
1267  real(RP), intent(out) :: f2h(ka,2)
1268 
1269  real(RP) :: dz1, dz2
1270  integer :: k
1271 
1272  do k = ks, ke-1
1273  dz1 = fz(k+1) - fz(k )
1274  dz2 = fz(k) - fz(k-1)
1275  f2h(k,1) = dz2 / ( dz1 + dz2 )
1276  f2h(k,2) = dz1 / ( dz1 + dz2 )
1277  end do
1278 
1279  return
1280  end subroutine get_f2h
1281 
1282 end module scale_atmos_phy_bl_mynn
scale_debug
module DEBUG
Definition: scale_debug.F90:11
scale_const::const_cpdry
real(rp), public const_cpdry
specific heat (dry air,constant pressure) [J/kg/K]
Definition: scale_const.F90:56
scale_atmos_saturation
module atmosphere / saturation
Definition: scale_atmos_saturation.F90:12
scale_atmos_grid_cartesc_index::ia
integer, public ia
of whole cells: x, local, with HALO
Definition: scale_atmos_grid_cartesC_index.F90:48
scale_atmos_phy_bl_mynn
Definition: scale_atmos_phy_bl_mynn.dox:123
scale_atmos_grid_cartesc_index::ja
integer, public ja
of whole cells: y, local, with HALO
Definition: scale_atmos_grid_cartesC_index.F90:49
scale_io::io_fid_conf
integer, public io_fid_conf
Config file ID.
Definition: scale_io.F90:55
scale_debug::check
subroutine, public check(current_line, v)
Undefined value checker.
Definition: scale_debug.F90:56
scale_const::const_karman
real(rp), parameter, public const_karman
von Karman constant
Definition: scale_const.F90:50
scale_tracer::tracer_name
character(len=h_short), dimension(qa_max), public tracer_name
Definition: scale_tracer.F90:46
scale_const::const_rdry
real(rp), public const_rdry
specific gas constant (dry air) [J/kg/K]
Definition: scale_const.F90:55
scale_const::const_undef
real(rp), public const_undef
Definition: scale_const.F90:41
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_atmos_grid_cartesc_index::ie
integer, public ie
end point of inner domain: x, local
Definition: scale_atmos_grid_cartesC_index.F90:54
scale_atmos_hydrometeor
module atmosphere / hydrometeor
Definition: scale_atmos_hydrometeor.F90:12
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_prc
module PROCESS
Definition: scale_prc.F90:11
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_const::const_grav
real(rp), public const_grav
standard acceleration of gravity [m/s2]
Definition: scale_const.F90:46
scale_matrix
module MATRIX
Definition: scale_matrix.F90:11
scale_const::const_undef2
integer, parameter, public const_undef2
undefined value (INT2)
Definition: scale_const.F90:38
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_const::const_epstvap
real(rp), public const_epstvap
1 / epsilon - 1
Definition: scale_const.F90:70
scale_const::const_rvap
real(rp), parameter, public const_rvap
specific gas constant (water vapor) [J/kg/K]
Definition: scale_const.F90:63
scale_prc::prc_abort
subroutine, public prc_abort
Abort Process.
Definition: scale_prc.F90:338
scale_atmos_hydrometeor::lhv
real(rp), public lhv
latent heat of vaporization for use [J/kg]
Definition: scale_atmos_hydrometeor.F90:126
scale_const
module CONSTANT
Definition: scale_const.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_prof
module profiler
Definition: scale_prof.F90:11
scale_const::const_eps
real(rp), public const_eps
small number
Definition: scale_const.F90:33
scale_precision
module PRECISION
Definition: scale_precision.F90:14
scale_atmos_grid_cartesc_index::ka
integer, public ka
of whole cells: z, local, with HALO
Definition: scale_atmos_grid_cartesC_index.F90:47
scale_const::const_pi
real(rp), public const_pi
pi
Definition: scale_const.F90:31
scale_io
module STDIO
Definition: scale_io.F90:10
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_precision::rp
integer, parameter, public rp
Definition: scale_precision.F90:41
scale_file_history
module file_history
Definition: scale_file_history.F90:15
Generated by   doxygen 1.8.13