SCALE-RM
scale_atmos_phy_tb_mynn.F90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
19 !-------------------------------------------------------------------------------
20 #include "inc_openmp.h"
21 module scale_atmos_phy_tb_mynn
22  !-----------------------------------------------------------------------------
23  !
24  !++ used modules
25  !
26  use scale_precision
27  use scale_stdio
28  use scale_prof
29  use scale_grid_index
30  use scale_tracer
31 
32 #include "macro_thermodyn.h"
33 
34 #if defined DEBUG || defined QUICKDEBUG
35  use scale_debug, only: &
36  check
37  use scale_const, only: &
38  undef => const_undef, &
39  iundef => const_undef2
40 #endif
41  !-----------------------------------------------------------------------------
42  implicit none
43  private
44  !-----------------------------------------------------------------------------
45  !
46  !++ Public procedure
47  !
48  public :: atmos_phy_tb_mynn_config
49  public :: atmos_phy_tb_mynn_setup
50  public :: atmos_phy_tb_mynn
51 
52  !-----------------------------------------------------------------------------
53  !
54  !++ Public parameters & variables
55  !
56  !-----------------------------------------------------------------------------
57  !
58  !++ Private procedure
59  !
60  !-----------------------------------------------------------------------------
61  !
62  !++ Private parameters & variables
63  !
64  real(RP), private, parameter :: OneOverThree = 1.0_rp / 3.0_rp
65  real(RP), private, parameter :: TKE_min = 1.e-10_rp
66  real(RP), private, parameter :: LT_min = 1.e-6_rp
67  real(RP), private, parameter :: Us_min = 1.e-6_rp
68 
69  real(RP), private :: A1
70  real(RP), private :: A2
71  real(RP), private, parameter :: B1 = 24.0_rp
72  real(RP), private, parameter :: B2 = 15.0_rp
73  real(RP), private :: C1
74  real(RP), private, parameter :: C2 = 0.75_rp
75  real(RP), private, parameter :: C3 = 0.352_rp
76  real(RP), private, parameter :: C5 = 0.2_rp
77  real(RP), private, parameter :: G1 = 0.235_rp
78  real(RP), private :: G2
79  real(RP), private :: F1
80  real(RP), private :: F2
81  real(RP), private :: Rf1
82  real(RP), private :: Rf2
83  real(RP), private :: Rfc
84  real(RP), private :: AF12
85  real(RP), private, parameter :: PrN = 0.74_rp
86 
87  real(RP), private :: SQRT_2PI
88  real(RP), private :: RSQRT_2PI
89  real(RP), private :: RSQRT_2
90 
91  integer, private :: I_TKE
92 
93  integer, private :: KE_PBL
94  logical, private :: ATMOS_PHY_TB_MYNN_TKE_INIT = .false.
95  real(RP), private :: ATMOS_PHY_TB_MYNN_N2_MAX = 1.e+3_rp
96  real(RP), private :: ATMOS_PHY_TB_MYNN_NU_MIN = 1.e-6_rp
97  real(RP), private :: ATMOS_PHY_TB_MYNN_NU_MAX = 10000.0_rp
98  real(RP), private :: ATMOS_PHY_TB_MYNN_KH_MIN = 1.e-6_rp
99  real(RP), private :: ATMOS_PHY_TB_MYNN_KH_MAX = 10000.0_rp
100  real(RP), private :: ATMOS_PHY_TB_MYNN_Lt_MAX = 700.0_rp ! ~ 0.23 * 3 km
101 
102  !-----------------------------------------------------------------------------
103 contains
104  !-----------------------------------------------------------------------------
106  subroutine atmos_phy_tb_mynn_config( &
107  TYPE_TB, &
108  I_TKE_out )
109  use scale_process, only: &
110  prc_mpistop
111  use scale_tracer, only: &
112  tracer_regist
113  implicit none
114 
115  character(len=*), intent(in) :: type_tb
116  integer, intent(out) :: i_tke_out
117  !---------------------------------------------------------------------------
118 
119  if( io_l ) write(io_fid_log,*)
120  if( io_l ) write(io_fid_log,*) '++++++ Module[Turbulence Tracer] / Categ[ATMOS PHYSICS] / Origin[SCALElib]'
121  if( io_l ) write(io_fid_log,*) '*** Tracers for Mellor-Yamada Nakanishi-Niino scheme'
122 
123  if ( type_tb /= 'MYNN' ) then
124  write(*,*) 'xxx ATMOS_PHY_TB_TYPE is not MYNN. Check!'
125  call prc_mpistop
126  endif
127 
128  call tracer_regist( i_tke, & ! [OUT]
129  1, & ! [IN]
130  (/ 'TKE_MYNN' /), & ! [IN]
131  (/ 'turbulent kinetic energy (MYNN)' /), & ! [IN]
132  (/ 'm2/s2' /) ) ! [IN]
133 
134  i_tke_out = i_tke
135 
136  return
137  end subroutine atmos_phy_tb_mynn_config
138 
139  !-----------------------------------------------------------------------------
141  subroutine atmos_phy_tb_mynn_setup( &
142  CDZ, CDX, CDY, CZ )
143  use scale_process, only: &
144  prc_mpistop
145  use scale_const, only: &
146  pi => const_pi
147  implicit none
148 
149  real(RP), intent(in) :: cdz(ka)
150  real(RP), intent(in) :: cdx(ia)
151  real(RP), intent(in) :: cdy(ja)
152  real(RP), intent(in) :: cz (ka,ia,ja)
153 
154  real(RP) :: atmos_phy_tb_mynn_pbl_max = 1.e+10_rp
155 
156  namelist / param_atmos_phy_tb_mynn / &
157  atmos_phy_tb_mynn_tke_init, &
158  atmos_phy_tb_mynn_pbl_max, &
159  atmos_phy_tb_mynn_n2_max, &
160  atmos_phy_tb_mynn_nu_min, &
161  atmos_phy_tb_mynn_nu_max, &
162  atmos_phy_tb_mynn_kh_min, &
163  atmos_phy_tb_mynn_kh_max, &
164  atmos_phy_tb_mynn_lt_max
165 
166  integer :: ierr
167  integer :: k, i, j
168  !---------------------------------------------------------------------------
169 
170  if( io_l ) write(io_fid_log,*)
171  if( io_l ) write(io_fid_log,*) '++++++ Module[Turbulence] / Categ[ATMOS PHYSICS] / Origin[SCALElib]'
172  if( io_l ) write(io_fid_log,*) '*** Mellor-Yamada Nakanishi-Niino scheme'
173 
174  !--- read namelist
175  rewind(io_fid_conf)
176  read(io_fid_conf,nml=param_atmos_phy_tb_mynn,iostat=ierr)
177  if( ierr < 0 ) then !--- missing
178  if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
179  elseif( ierr > 0 ) then !--- fatal error
180  write(*,*) 'xxx Not appropriate names in namelist PARAM_ATMOS_PHY_TB_MYNN. Check!'
181  call prc_mpistop
182  endif
183  if( io_nml ) write(io_fid_nml,nml=param_atmos_phy_tb_mynn)
184 
185  do k = ks, ke-1
186  do j = js, je
187  do i = is, ie
188  if ( atmos_phy_tb_mynn_pbl_max >= cz(k,i,j) ) then
189  ke_pbl = k
190  end if
191  end do
192  end do
193  end do
194 
195  a1 = b1 * (1.0_rp - 3.0_rp * g1) / 6.0_rp
196  a2 = 1.0_rp / (3.0_rp * g1 * b1**(1.0_rp/3.0_rp) * prn )
197  c1 = g1 - 1.0_rp / ( 3.0_rp * a1 * b1**(1.0_rp/3.0_rp) )
198  g2 = ( 2.0_rp * a1 * (3.0_rp - 2.0_rp * c2) + b2 * (1.0_rp - c3) ) / b1
199  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)
200  f2 = b1 * (g1 + g2) - 3.0_rp * a1 * (1.0_rp - c2)
201 
202  rf1 = b1 * (g1 - c1) / f1
203  rf2 = b1 * g1 / f2
204  rfc = g1 / (g1 + g2)
205 
206  af12 = a1 * f1 / ( a2 * f2 )
207 
208  sqrt_2pi = sqrt( 2.0_rp * pi )
209  rsqrt_2pi = 1.0_rp / sqrt_2pi
210  rsqrt_2 = 1.0_rp / sqrt( 2.0_rp )
211 
212  return
213  end subroutine atmos_phy_tb_mynn_setup
214 
215  !-----------------------------------------------------------------------------
216  subroutine atmos_phy_tb_mynn( &
217  qflx_sgs_momz, qflx_sgs_momx, qflx_sgs_momy, &
218  qflx_sgs_rhot, qflx_sgs_rhoq, &
219  RHOQ_t, &
220  Nu, Ri, Pr, &
221  MOMZ, MOMX, MOMY, RHOT, DENS, QTRC, N2_in, &
222  SFLX_MW, SFLX_MU, SFLX_MV, SFLX_SH, SFLX_Q, &
223  GSQRT, J13G, J23G, J33G, MAPF, dt )
224  use scale_precision
225  use scale_grid_index
226  use scale_tracer
227  use scale_const, only: &
228  grav => const_grav, &
229  r => const_rdry, &
230  rvap => const_rvap, &
231  cpdry => const_cpdry, &
232  cpvap => const_cpvap, &
233  lhv0 => const_lhv0, &
234  epstvap => const_epstvap
235  use scale_grid, only: &
236  rcdz => grid_rcdz, &
237  rfdz => grid_rfdz
238  use scale_grid_real, only: &
239  cz => real_cz
240  use scale_gridtrans, only: &
241  i_xyz, &
242  i_xyw, &
243  i_xvw, &
244  i_uyw
245  use scale_comm, only: &
246  comm_vars, &
247  comm_wait
248  use scale_atmos_phy_tb_common, only: &
249  diffusion_solver => atmos_phy_tb_diffusion_solver
250  use scale_atmos_thermodyn, only: &
251  atmos_thermodyn_temp_pres
252  use scale_atmos_hydrometeor, only: &
253  hydrometeor_lhv => atmos_hydrometeor_lhv, &
254  hydrometeor_lhs => atmos_hydrometeor_lhs, &
255  i_qv, &
256  i_qc, &
257  i_qi
258  use scale_atmos_saturation, only: &
259  atmos_saturation_dens2qsat => atmos_saturation_dens2qsat_all
260 #ifdef MORE_HIST
261  use scale_history, only: &
262  hist_in
263 #endif
264  implicit none
265 
266  real(RP), intent(out) :: qflx_sgs_momz(ka,ia,ja,3)
267  real(RP), intent(out) :: qflx_sgs_momx(ka,ia,ja,3)
268  real(RP), intent(out) :: qflx_sgs_momy(ka,ia,ja,3)
269  real(RP), intent(out) :: qflx_sgs_rhot(ka,ia,ja,3)
270  real(RP), intent(out) :: qflx_sgs_rhoq(ka,ia,ja,3,qa)
271 
272  real(RP), intent(inout) :: rhoq_t (ka,ia,ja,qa) ! tendency of rho * QTRC
273 
274  real(RP), intent(out) :: nu (ka,ia,ja) ! eddy viscosity (center)
275  real(RP), intent(out) :: ri (ka,ia,ja) ! Richardson number
276  real(RP), intent(out) :: pr (ka,ia,ja) ! Plandtle number
277 
278  real(RP), intent(in) :: momz (ka,ia,ja)
279  real(RP), intent(in) :: momx (ka,ia,ja)
280  real(RP), intent(in) :: momy (ka,ia,ja)
281  real(RP), intent(in) :: rhot (ka,ia,ja)
282  real(RP), intent(in) :: dens (ka,ia,ja)
283  real(RP), intent(in) :: qtrc (ka,ia,ja,qa)
284  real(RP), intent(in) :: n2_in (ka,ia,ja)
285 
286  real(RP), intent(in) :: sflx_mw (ia,ja)
287  real(RP), intent(in) :: sflx_mu (ia,ja)
288  real(RP), intent(in) :: sflx_mv (ia,ja)
289  real(RP), intent(in) :: sflx_sh (ia,ja)
290  real(RP), intent(in) :: sflx_q (ia,ja,qa)
291 
292  real(RP), intent(in) :: gsqrt (ka,ia,ja,7)
293  real(RP), intent(in) :: j13g (ka,ia,ja,7)
294  real(RP), intent(in) :: j23g (ka,ia,ja,7)
295  real(RP), intent(in) :: j33g
296  real(RP), intent(in) :: mapf (ia,ja,2,4)
297  real(DP), intent(in) :: dt
298 
299  real(RP) :: u (ka,ia,ja)
300  real(RP) :: v (ka,ia,ja)
301  real(RP) :: n2 (ka,ia,ja) ! squared Brunt-Vaisala frequency
302  real(RP) :: phin (ka,ia,ja)
303 
304  real(RP) :: sm (ka,ia,ja)
305  real(RP) :: sh (ka,ia,ja)
306  real(RP) :: l (ka,ia,ja)
307  real(RP) :: q (ka,ia,ja)
308  real(RP) :: dudz2(ka,ia,ja)
309  real(RP) :: q2_2 (ka,ia,ja)
310  real(RP) :: kh
311  real(RP) :: rhokh(ka,ia,ja)
312 
313  real(RP) :: sflx_pt(ia,ja) ! surface potential temperature flux
314 
315  real(RP) :: a(ka,ia,ja)
316  real(RP) :: b(ka,ia,ja)
317  real(RP) :: c(ka,ia,ja)
318  real(RP) :: d(ka,ia,ja)
319  real(RP) :: ap
320  real(RP) :: tke_n(ka,ia,ja)
321  real(RP) :: tke
322 
323  real(RP) :: pott(ka,ia,ja)
324  real(RP) :: potv(ka,ia,ja)
325  real(RP) :: potl(ka,ia,ja)
326  real(RP) :: teml(ka,ia,ja)
327 
328  real(RP) :: qw(ka,ia,ja)
329  real(RP) :: qdry(ka,ia,ja)
330  real(RP) :: qv
331  real(RP) :: ql
332  real(RP) :: qs
333 
334  real(RP) :: temp(ka,ia,ja)
335  real(RP) :: pres(ka,ia,ja)
336 
337  real(RP) :: lhv (ka,ia,ja)
338  real(RP) :: lhs (ka,ia,ja)
339 
340  real(RP) :: cptot
341 
342  real(RP) :: ac
343 
344  real(RP) :: q1
345  real(RP) :: qsl(ka)
346  real(RP) :: lhvl(ka)
347  real(RP) :: dqsl
348  real(RP) :: sigma_s
349  real(RP) :: rr
350  real(RP) :: rt
351  real(RP) :: betat
352  real(RP) :: betaq
353  real(RP) :: aa
354  real(RP) :: bb
355  real(RP) :: cc
356 
357 #ifdef MORE_HIST
358  real(RP) :: gen(ka,ia,ja)
359 #endif
360 
361  integer :: iis, iie
362  integer :: jjs, jje
363 
364  integer :: k, i, j, iq
365  !---------------------------------------------------------------------------
366 
367  if( io_l ) write(io_fid_log, *) "*** Atmos physics step: Turbulence (MYNN)"
368 
369 #ifdef DEBUG
370  pott(:,:,:) = undef
371  potl(:,:,:) = undef
372  temp(:,:,:) = undef
373 #endif
374 
375 #if defined DEBUG || defined QUICKDEBUG
376  qflx_sgs_momz(ks:ke, 1:is-1, : ,:) = undef
377  qflx_sgs_momz(ks:ke,ie+1:ia , : ,:) = undef
378  qflx_sgs_momz(ks:ke, : , 1:js-1,:) = undef
379  qflx_sgs_momz(ks:ke, : ,je+1:ja ,:) = undef
380  qflx_sgs_momx(ks:ke, 1:is-1, : ,:) = undef
381  qflx_sgs_momx(ks:ke,ie+1:ia , : ,:) = undef
382  qflx_sgs_momx(ks:ke, : , 1:js-1,:) = undef
383  qflx_sgs_momx(ks:ke, : ,je+1:ja ,:) = undef
384  qflx_sgs_momy(ks:ke, 1:is-1, : ,:) = undef
385  qflx_sgs_momy(ks:ke,ie+1:ia , : ,:) = undef
386  qflx_sgs_momy(ks:ke, : , 1:js-1,:) = undef
387  qflx_sgs_momy(ks:ke, : ,je+1:ja ,:) = undef
388 #endif
389 
390 !OCL XFILL
391  do j = js , je
392  do i = is-1, ie
393  do k = ks-1, ke
394  qflx_sgs_momz(k,i,j,xdir) = 0.0_rp
395  end do
396  end do
397  end do
398 !OCL XFILL
399  do j = js , je
400  do i = is , ie+1
401  do k = ks-1, ke+1
402  qflx_sgs_momx(k,i,j,xdir) = 0.0_rp
403  end do
404  end do
405  end do
406 
407 !OCL XFILL
408  do j = js , je
409  do i = is-1, ie
410  do k = ks-1, ke+1
411  qflx_sgs_momy(k,i,j,xdir) = 0.0_rp
412  end do
413  end do
414  end do
415 !OCL XFILL
416  do j = js , je
417  do i = is-1, ie
418  do k = ks-1, ke+1
419  qflx_sgs_rhot(k,i,j,xdir) = 0.0_rp
420  end do
421  end do
422  end do
423 !OCL XFILL
424  do iq = 1, qa
425  do j = js , je
426  do i = is-1, ie
427  do k = ks-1, ke+1
428  qflx_sgs_rhoq(k,i,j,xdir,iq) = 0.0_rp
429  end do
430  end do
431  end do
432  end do
433 
434 !OCL XFILL
435  do j = js-1, je
436  do i = is , ie
437  do k = ks-1, ke
438  qflx_sgs_momz(k,i,j,ydir) = 0.0_rp
439  end do
440  end do
441  end do
442 !OCL XFILL
443  do j = js-1, je
444  do i = is , ie
445  do k = ks-1, ke+1
446  qflx_sgs_momx(k,i,j,ydir) = 0.0_rp
447  end do
448  end do
449  end do
450 !OCL XFILL
451  do j = js , je+1
452  do i = is , ie
453  do k = ks-1, ke+1
454  qflx_sgs_momy(k,i,j,ydir) = 0.0_rp
455  end do
456  end do
457  end do
458 !OCL XFILL
459  do j = js-1, je
460  do i = is , ie
461  do k = ks-1, ke+1
462  qflx_sgs_rhot(k,i,j,ydir) = 0.0_rp
463  end do
464  end do
465  end do
466 !OCL XFILL
467  do iq = 1, qa
468  do j = js-1, je
469  do i = is , ie
470  do k = ks-1, ke+1
471  qflx_sgs_rhoq(k,i,j,ydir,iq) = 0.0_rp
472  end do
473  end do
474  end do
475  end do
476 !OCL XFILL
477  do j = js, je
478  do i = is, ie
479  do k = ks-1, ke+1
480  qflx_sgs_momz(k,i,j,zdir) = 0.0_rp
481  end do
482  end do
483  end do
484 
485 
486 
487  do jjs = js, je, jblock
488  jje = jjs+jblock-1
489  do iis = is, ie, iblock
490  iie = iis+iblock-1
491 
492 !OCL XFILL
493  do j = jjs , jje+1
494  do i = iis-1, iie+1
495  do k = ks, ke_pbl+1
496  u(k,i,j) = 0.5_rp * ( momx(k,i,j) + momx(k,i-1,j) ) / dens(k,i,j)
497  end do
498  end do
499  end do
500 
501 !OCL XFILL
502  do j = jjs-1, jje+1
503  do i = iis , iie+1
504  do k = ks, ke_pbl+1
505  v(k,i,j) = 0.5_rp * ( momy(k,i,j) + momy(k,i,j-1) ) / dens(k,i,j)
506  end do
507  end do
508  end do
509 
510  end do
511  end do
512 
513  call atmos_thermodyn_temp_pres( temp, pres, & ! (out)
514  dens, rhot, qtrc, & ! (in)
515  tracer_cv, tracer_r, tracer_mass ) ! (in)
516 
517 
518  call hydrometeor_lhv( lhv(:,:,:), temp(:,:,:) )
519  call hydrometeor_lhs( lhs(:,:,:), temp(:,:,:) )
520 
521 !OCL LOOP_NOFUSION,PREFETCH_SEQUENTIAL(SOFT),SWP
522  !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
523  !$omp shared(JS,JE,IS,IE,KS,KE_PBL,QA,I_QV,I_QC,I_QI) &
524  !$omp shared(GRAV,CPdry,EPSTvap,ATMOS_PHY_TB_MYNN_N2_MAX,CZ) &
525  !$omp shared(QTRC,TRACER_MASS,RHOT,DENS,N2_in,SFLX_PT,SFLX_SH) &
526  !$omp shared(Qdry,Qw,LHV,LHS,POTT,POTL,temp,TEML,POTV,n2,U,V,Ri,dudz2) &
527  !$omp private(i,j,k,iq) &
528  !$omp private(qv,ql,qs,CPtot)
529  do j = js, je
530  do i = is, ie
531 
532  do k = ks, ke_pbl+1
533 
534  qv = 0.0_rp
535  if ( i_qv > 0 ) qv = qtrc(k,i,j,i_qv)
536 
537  ql = 0.0_rp
538  if ( i_qc > 0 ) ql = qtrc(k,i,j,i_qc)
539 ! do iq = QWS, QWE
540 ! ql = ql + QTRC(k,i,j,iq)
541 ! end do
542  qs = 0.0_rp
543  if ( i_qi > 0 ) qs = qtrc(k,i,j,i_qi)
544 ! do iq = QIS, QIE
545 ! qs = qs + QTRC(k,i,j,iq)
546 ! end do
547  calc_qdry(qdry(k,i,j), qtrc, tracer_mass, k, i, j, iq)
548 
549  qw(k,i,j) = qv + ql + qs
550 
551  cptot = cpdry * qdry(k,i,j) + cpvap * qv
552 
553  pott(k,i,j) = rhot(k,i,j) / dens(k,i,j)
554  ! liquid water potential temperature
555  potl(k,i,j) = pott(k,i,j) * (1.0_rp - 1.0_rp * (lhv(k,i,j) * ql + lhs(k,i,j) * qs) / ( temp(k,i,j) * cptot ) )
556  teml(k,i,j) = potl(k,i,j) * temp(k,i,j) / pott(k,i,j)
557 
558  ! virtual potential temperature for derivertive
559 ! POTV(k,i,j) = ( 1.0_RP + EPSTvap * Qw(k,i,j) ) * POTL(k,i,j)
560  potv(k,i,j) = ( qdry(k,i,j) + (epstvap+1.0_rp) * qv ) * pott(k,i,j)
561 
562  end do
563 
564  sflx_pt(i,j) = sflx_sh(i,j) / ( cpdry * dens(ks,i,j) ) &
565  * pott(ks,i,j) / temp(ks,i,j)
566 
567  do k = ks, ke_pbl
568  n2(k,i,j) = min( atmos_phy_tb_mynn_n2_max, n2_in(k,i,j) )
569  end do
570 
571  dudz2(ks,i,j) = ( ( u(ks+1,i,j) - u(ks,i,j) )**2 + ( v(ks+1,i,j) - v(ks,i,j) )**2 ) &
572  / ( cz(ks+1,i,j) - cz(ks,i,j) )**2
573  do k = ks+1, ke_pbl
574  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 ) &
575  / ( cz(k+1,i,j) - cz(k-1,i,j) )**2
576  end do
577 
578  do k = ks, ke_pbl
579  ri(k,i,j) = n2(k,i,j) / max(dudz2(k,i,j), 1e-10_rp)
580  end do
581 
582  end do
583  end do
584 
585  if ( atmos_phy_tb_mynn_tke_init ) then
586 !OCL XFILL
587  do j = js, je
588  do i = is, ie
589  do k = ks, ke_pbl
590  q(k,i,j) = sqrt( tke_min*2.0_rp )
591  end do
592  end do
593  end do
594  else
595 !OCL XFILL
596  do j = js, je
597  do i = is, ie
598  do k = ks, ke_pbl
599  q(k,i,j) = sqrt( max(qtrc(k,i,j,i_tke), tke_min)*2.0_rp )
600  end do
601  end do
602  end do
603  end if
604 
605  ! length
606  call get_length( &
607  l, & ! (out)
608  dens, & ! (in)
609  q, n2, & ! (in)
610  sflx_mu, sflx_mv, sflx_pt, & ! (in)
611  pott ) ! (in)
612 
613  call get_q2_level2( &
614  q2_2, & ! (out)
615  dudz2, ri, l ) ! (in)
616 
617  if ( atmos_phy_tb_mynn_tke_init ) then
618 !OCL XFILL
619  do j = js, je
620  do i = is, ie
621  do k = ks, ke_pbl
622  tke = max(q2_2(k,i,j) * 0.5_rp, tke_min)
623  q(k,i,j) = sqrt( tke * 2.0_rp )
624  end do
625  end do
626  end do
627 
628  atmos_phy_tb_mynn_tke_init = .false.
629  end if
630 
631  call get_smsh( sm, sh, & ! (out)
632  q, q2_2, & ! (in)
633  l, n2, dudz2 ) ! (in)
634 
635 #if defined DEBUG || defined QUICKDEBUG
636  do j = 1, js-1
637  do i = 1, ia
638  do k = ks, ke
639  teml(k,i,j) = 300.0_rp
640  end do
641  end do
642  end do
643  do j = je+1, ja
644  do i = 1, ia
645  do k = ks, ke
646  teml(k,i,j) = 300.0_rp
647  end do
648  end do
649  end do
650  do j = 1, ja
651  do i = 1, is-1
652  do k = ks, ke
653  teml(k,i,j) = 300.0_rp
654  end do
655  end do
656  do i = ie+1, ia
657  do k = ks, ke
658  teml(k,i,j) = 300.0_rp
659  end do
660  end do
661  end do
662 #endif
663 
664 !OCL LOOP_NOFUSION,PREFETCH_SEQUENTIAL(SOFT),SWP
665  !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
666  !$omp private(i,j,k) &
667  !$omp private(CPtot,Qsl,dQsl,aa,bb,ac,sigma_s,Q1,RR,Ql,cc,Rt,betat,betaq,LHVL) &
668  !$omp shared(JS,JE,IS,IE,KS,KE,KE_PBL) &
669  !$omp shared(CZ,rsqrt_2,rsqrt_2pi,sqrt_2pi,GRAV,CPdry,EPSTvap,ATMOS_PHY_TB_MYNN_N2_MAX) &
670  !$omp shared(DENS,l,Qdry,TEMP,TEML,POTT,POTL,POTV,q2_2,Qw,q,dudz2,Ri,n2,sh)
671  do j = js, je
672  do i = is, ie
673 
674  call atmos_saturation_dens2qsat( qsl(:), & ! (out)
675  teml(:,i,j), dens(:,i,j) ) ! (in)
676 
677  call hydrometeor_lhv( lhvl(:), teml(:,i,j) )
678 
679  do k = ks+1, ke_pbl
680 
681  dqsl = ( qsl(k) * lhvl(k) / ( rvap * teml(k,i,j) ) - qsl(k) ) / teml(k,i,j)
682  cptot = qdry(k,i,j) * cpdry + qsl(k) * cpvap
683  aa = 1.0_rp / ( 1.0_rp + lhvl(k)/cptot * dqsl )
684  bb = temp(k,i,j) / pott(k,i,j) * dqsl
685  ac = min( q(k,i,j)/sqrt(q2_2(k,i,j)), 1.0_rp )
686  sigma_s = max( sqrt( 0.25_rp * aa**2 * l(k,i,j)**2 * ac * b2 * sh(k,i,j) ) &
687  * abs( qw(k+1,i,j) - qw(k-1,i,j) - bb * ( potl(k+1,i,j)-potl(k-1,i,j) ) ) &
688  / ( cz(k+1,i,j) - cz(k-1,i,j) ), &
689  1e-10_rp )
690  q1 = aa * ( qw(k,i,j) - qsl(k) ) * 0.5_rp / sigma_s
691  rr = min( max( 0.5_rp * ( 1.0_rp + erf(q1*rsqrt_2) ), 0.0_rp ), 1.0_rp )
692 #if defined(PGI) || defined(SX)
693  ql = min( max( 2.0_rp * sigma_s * ( rr * q1 + rsqrt_2pi * exp( -min( 0.5_rp*q1**2, 1.e+3_rp ) ) ), & ! apply exp limiter
694 #else
695  ql = min( max( 2.0_rp * sigma_s * ( rr * q1 + rsqrt_2pi * exp(-0.5_rp*q1**2) ), &
696 #endif
697  0.0_rp ), &
698  qw(k,i,j) * 0.5_rp )
699  cc = ( 1.0_rp + epstvap * qw(k,i,j) - (1.0_rp+epstvap) * ql ) * pott(k,i,j)/temp(k,i,j) * lhvl(k) / cptot &
700  - (1.0_rp+epstvap) * pott(k,i,j)
701 #if defined(PGI) || defined(SX)
702  rt = min( max( rr - ql / (2.0_rp*sigma_s*sqrt_2pi) * exp( -min( 0.5_rp*q1**2, 1.e+3_rp ) ), 0.0_rp ), 1.0_rp ) ! apply exp limiter
703 #else
704  rt = min( max( rr - ql / (2.0_rp*sigma_s*sqrt_2pi) * exp(-q1**2 * 0.5_rp), 0.0_rp ), 1.0_rp )
705 #endif
706  betat = 1.0_rp + epstvap * qw(k,i,j) - (1.0_rp+epstvap) * ql - rt * aa * bb * cc
707  betaq = epstvap * pott(k,i,j) + rt * aa * cc
708  n2(k,i,j) = min(atmos_phy_tb_mynn_n2_max, &
709  grav * ( ( potl(k+1,i,j) - potl(k-1,i,j) ) * betat &
710  + ( qw(k+1,i,j) - qw(k-1,i,j) ) * betaq ) &
711  / ( cz(k+1,i,j) - cz(k-1,i,j) ) / potv(k,i,j) )
712  end do
713  n2(ks,i,j) = n2(ks+1,i,j)
714 
715  do k = ks, ke_pbl
716  ri(k,i,j) = n2(k,i,j) / max(dudz2(k,i,j), 1e-10_rp)
717  end do
718  do k = ke_pbl+1, ke
719  ri(k,i,j) = 0.0_rp
720  n2(k,i,j) = 0.0_rp
721  end do
722 
723  end do
724  end do
725 
726  ! length
727  call get_length( &
728  l, & ! (out)
729  dens, & ! (in)
730  q, n2, & ! (in)
731  sflx_mu, sflx_mv, sflx_pt, & ! (in)
732  pott ) ! (in)
733 
734  call get_q2_level2( &
735  q2_2, & ! (out)
736  dudz2, ri, l ) ! (in)
737 
738  call get_smsh( &
739  sm, sh, & ! (out)
740  q, q2_2, & ! (in)
741  l, n2, dudz2 ) ! (in)
742 
743 
744  !$omp parallel do default(none) private(i,j,k,Kh) OMP_SCHEDULE_ collapse(2) &
745  !$omp shared(JS,JE,IS,IE,Nu,l,KE_PBL,q,sm,ATMOS_PHY_TB_MYNN_NU_MAX,ATMOS_PHY_TB_MYNN_NU_MIN) &
746  !$omp shared(KA,sh,ATMOS_PHY_TB_MYNN_KH_MAX,ATMOS_PHY_TB_MYNN_KH_MIN,RHOKh,DENS,Pr,KE,KS)
747  do j = js, je
748  do i = is, ie
749  do k = 1, ks-1
750  nu(k,i,j) = 0.0_rp
751  end do
752  do k = ks, ke_pbl
753  nu(k,i,j) = max( min( l(k,i,j) * q(k,i,j) * sm(k,i,j), &
754  atmos_phy_tb_mynn_nu_max ), &
755  atmos_phy_tb_mynn_nu_min )
756  end do
757  do k = ke_pbl+1, ka
758  nu(k,i,j) = 0.0_rp
759  end do
760 
761  do k = ks, ke_pbl
762  kh = max( min( l(k,i,j) * q(k,i,j) * sh(k,i,j), &
763  atmos_phy_tb_mynn_kh_max ), &
764  atmos_phy_tb_mynn_kh_min )
765  rhokh(k,i,j) = dens(k,i,j) * kh
766  pr(k,i,j) = nu(k,i,j) / kh
767  end do
768  do k = ke_pbl+1, ke
769  pr(k,i,j) = 1.0_rp
770  end do
771  end do
772  end do
773 
774  call comm_vars( nu, 1 )
775 
776  ! time integration
777 
778  ! for velocities
779 !OCL INDEPENDENT
780  !$omp parallel do default(none) &
781  !$omp shared(JS,JE,IS,IE,KS,KE_PBL,dt,DENS,Nu,RFDZ,GSQRT,I_XYW,a,c,b,d,RCDZ,I_XYZ,U,SFLX_MU) &
782  !$omp shared(phiN) &
783  !$omp private(i,j,k,ap) OMP_SCHEDULE_ collapse(2)
784  do j = js, je
785  do i = is, ie
786 
787  ap = - dt * 0.5_rp * ( dens(ks ,i,j)*nu(ks ,i,j) &
788  + dens(ks+1,i,j)*nu(ks+1,i,j) ) &
789  * rfdz(ks) / gsqrt(ks,i,j,i_xyw)
790  a(ks,i,j) = ap * rcdz(ks) / ( dens(ks,i,j) * gsqrt(ks,i,j,i_xyz) )
791  c(ks,i,j) = 0.0_rp
792  b(ks,i,j) = - a(ks,i,j) + 1.0_rp
793  do k = ks+1, ke_pbl-1
794  c(k,i,j) = ap * rcdz(k) / ( dens(k,i,j) * gsqrt(k,i,j,i_xyz) )
795  ap = - dt * 0.5_rp * ( dens(k ,i,j)*nu(k ,i,j) &
796  + dens(k+1,i,j)*nu(k+1,i,j) ) &
797  * rfdz(k) / gsqrt(k,i,j,i_xyw)
798  a(k,i,j) = ap * rcdz(k) / ( dens(k,i,j) * gsqrt(k,i,j,i_xyz) )
799  b(k,i,j) = - a(k,i,j) - c(k,i,j) + 1.0_rp
800  end do
801  a(ke_pbl,i,j) = 0.0_rp
802  c(ke_pbl,i,j) = ap * rcdz(ke_pbl) / ( dens(ke_pbl,i,j) * gsqrt(ke_pbl,i,j,i_xyz) )
803  b(ke_pbl,i,j) = - c(ke_pbl,i,j) + 1.0_rp
804 
805  d(ks,i,j) = u(ks,i,j) + dt * sflx_mu(i,j) * rcdz(ks) / ( dens(ks,i,j) * gsqrt(ks,i,j,i_xyz) )
806  do k = ks+1, ke_pbl
807  d(k,i,j) = u(k,i,j)
808  end do
809  call diffusion_solver( &
810  phin(:,i,j), & ! (out)
811  a(:,i,j), b(:,i,j), c(:,i,j), d(:,i,j), & ! (in)
812  ke_pbl ) ! (in)
813  end do
814  end do
815 
816  call comm_vars( phin, 2 )
817  call comm_wait( nu, 1 )
818  call comm_wait( phin, 2 )
819 
820 !OCL INDEPENDENT
821  do jjs = js, je, jblock
822  jje = jjs+jblock-1
823  do iis = is, ie, iblock
824  iie = iis+iblock-1
825 
826  do j = jjs, jje
827  do i = iis, iie
828  do k = ks, ke_pbl-1
829  qflx_sgs_momx(k,i,j,zdir) = - 0.03125_rp & ! 1/4/4/2
830  * ( dens(k,i,j) + dens(k+1,i,j) + dens(k,i+1,j) + dens(k+1,i+1,j) ) &
831  * ( nu(k,i,j) + nu(k+1,i,j) + nu(k,i+1,j) + nu(k+1,i+1,j) ) &
832  * ( (phin(k+1,i,j)+phin(k+1,i+1,j)) - (phin(k,i,j)+phin(k,i+1,j)) ) &
833  * j33g * rfdz(k) / gsqrt(k,i,j,i_uyw)
834  end do
835  end do
836  end do
837  do j = jjs, jje
838  do i = iis, iie
839  qflx_sgs_momx(ks-1,i,j,zdir) = 0.0_rp
840  end do
841  end do
842  do j = jjs, jje
843  do i = iis, iie
844  do k = ke_pbl, ke
845  qflx_sgs_momx(k,i,j,zdir) = 0.0_rp
846  end do
847  end do
848  end do
849 
850 
851  ! integration V
852 !OCL INDEPENDENT
853  !$omp parallel do default(none) &
854  !$omp shared(JJS,JJE,IIS,IIE,KS,KE_PBL,V,dt,SFLX_MV,RCDZ,DENS,GSQRT,I_XYZ,phiN,a,b,c,d) &
855  !$omp private(i,j,k) OMP_SCHEDULE_ collapse(2)
856  do j = jjs, jje
857  do i = iis, iie
858  d(ks,i,j) = v(ks,i,j) + dt * sflx_mv(i,j) * rcdz(ks) / ( dens(ks,i,j) * gsqrt(ks,i,j,i_xyz) )
859  do k = ks+1, ke_pbl
860  d(k,i,j) = v(k,i,j)
861  end do
862  call diffusion_solver( &
863  phin(:,i,j), & ! (out)
864  a(:,i,j), b(:,i,j), c(:,i,j), d(:,i,j), & ! (in)
865  ke_pbl ) ! (in)
866  end do
867  end do
868 
869  end do
870  end do
871 
872  call comm_vars( phin, 1 )
873  call comm_wait( phin, 1 )
874 
875  do jjs = js, je, jblock
876  jje = jjs+jblock-1
877  do iis = is, ie, iblock
878  iie = iis+iblock-1
879 
880  do j = jjs, jje
881  do i = iis, iie
882  qflx_sgs_momy(ks-1,i,j,zdir) = 0.0_rp
883  do k = ks, ke_pbl-1
884  qflx_sgs_momy(k,i,j,zdir) = - 0.03125_rp & ! 1/4/4/2
885  * ( dens(k,i,j) + dens(k+1,i,j) + dens(k,i,j+1) + dens(k+1,i,j+1) ) &
886  * ( nu(k,i,j) + nu(k+1,i,j) + nu(k,i,j+1) + nu(k+1,i,j+1) ) &
887  * ( (phin(k+1,i,j)+phin(k+1,i,j+1)) - (phin(k,i,j)+phin(k,i,j+1)) ) &
888  * j33g * rfdz(k) / gsqrt(k,i,j,i_xvw)
889  end do
890  do k = ke_pbl, ke
891  qflx_sgs_momy(k,i,j,zdir) = 0.0_rp
892  end do
893  end do
894  end do
895 
896 
897  ! for scalars
898  ! integration POTT
899 !OCL INDEPENDENT
900  do j = jjs, jje
901  do i = iis, iie
902  ap = - dt * 0.5_rp * ( rhokh(ks,i,j) + rhokh(ks+1,i,j) ) &
903  * rfdz(ks) / gsqrt(ks,i,j,i_xyw)
904  a(ks,i,j) = ap * rcdz(ks) / (dens(ks,i,j) * gsqrt(ks,i,j,i_xyz) )
905  c(ks,i,j) = 0.0_rp
906  b(ks,i,j) = - a(ks,i,j) + 1.0_rp
907  do k = ks+1, ke_pbl-1
908  c(k,i,j) = ap * rcdz(k) / (dens(k,i,j) * gsqrt(k,i,j,i_xyz) )
909  ap = - dt * 0.5_rp * ( rhokh(k,i,j) + rhokh(k+1,i,j) ) &
910  * rfdz(k) / gsqrt(k,i,j,i_xyw)
911  a(k,i,j) = ap * rcdz(k) / ( dens(k,i,j) * gsqrt(k,i,j,i_xyz) )
912  b(k,i,j) = - a(k,i,j) - c(k,i,j) + 1.0_rp
913  end do
914  a(ke_pbl,i,j) = 0.0_rp
915  c(ke_pbl,i,j) = ap * rcdz(ke_pbl) / (dens(ke_pbl,i,j) * gsqrt(ke_pbl,i,j,i_xyz) )
916  b(ke_pbl,i,j) = - c(ke_pbl,i,j) + 1.0_rp
917  d(ks,i,j) = pott(ks,i,j) + dt * sflx_pt(i,j) * rcdz(ks) / gsqrt(ks,i,j,i_xyz)
918  do k = ks+1, ke_pbl
919  d(k,i,j) = pott(k,i,j)
920  end do
921  call diffusion_solver( &
922  phin(:,i,j), & ! (out)
923  a(:,i,j), b(:,i,j), c(:,i,j), d(:,i,j), & ! (in)
924  ke_pbl ) ! (in)
925  end do
926  end do
927 
928  !$omp parallel do default(none) &
929  !$omp shared(JJS,JJE,IIS,IIE,KS,KE_PBL,RHOKh,J33G,phiN,RFDZ,GSQRT,I_XYW,KE,qflx_sgs_rhot) &
930  !$omp private(i,j,k) OMP_SCHEDULE_ collapse(2)
931  do j = jjs, jje
932  do i = iis, iie
933  qflx_sgs_rhot(ks-1,i,j,zdir) = 0.0_rp
934  do k = ks, ke_pbl-1
935  qflx_sgs_rhot(k,i,j,zdir) = - 0.5_rp & ! 1/2
936  * ( rhokh(k,i,j) + rhokh(k+1,i,j) ) &
937  * j33g * ( phin(k+1,i,j) - phin(k,i,j) ) * rfdz(k) / gsqrt(k,i,j,i_xyw)
938  end do
939  do k = ke_pbl, ke
940  qflx_sgs_rhot(k,i,j,zdir) = 0.0_rp
941  end do
942  end do
943  end do
944 
945 
946  ! TRACER
947  do iq = 1, qa
948 
949  if ( iq == i_tke ) then
950  qflx_sgs_rhoq(:,:,:,zdir,iq) = 0.0_rp
951  cycle
952  end if
953 
954  if( .NOT. tracer_advc(iq) ) cycle
955 
956 !OCL INDEPENDENT
957  !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
958  !$omp shared(JJS,JJE,IIS,IIE,iq,I_QV,d,QTRC,KS,dt,SFLX_Q,RCDZ,DENS,GSQRT,I_XYZ,KE_PBL) &
959  !$omp shared(a,b,c,phiN)
960  do j = jjs, jje
961  do i = iis, iie
962  d(ks,i,j) = qtrc(ks,i,j,iq) &
963  + dt * sflx_q(i,j,iq) * rcdz(ks) / ( dens(ks,i,j) * gsqrt(ks,i,j,i_xyz) )
964 
965  do k = ks+1, ke_pbl
966  d(k,i,j) = qtrc(k,i,j,iq)
967  end do
968 
969  call diffusion_solver( &
970  phin(:,i,j), & ! (out)
971  a(:,i,j), b(:,i,j), c(:,i,j), d(:,i,j), & ! (in)
972  ke_pbl ) ! (in)
973  end do
974  end do
975 
976  !$omp parallel do default(none) &
977  !$omp shared(JJS,JJE,IIS,IIE,KS,KE_PBL,iq,RHOKh,J33G,phiN,RFDZ,GSQRT,KE,I_XYW) &
978  !$omp shared(qflx_sgs_rhoq) &
979  !$omp private(i,j,k) OMP_SCHEDULE_ collapse(2)
980  do j = jjs, jje
981  do i = iis, iie
982  qflx_sgs_rhoq(ks-1,i,j,zdir,iq) = 0.0_rp
983  do k = ks, ke_pbl-1
984  qflx_sgs_rhoq(k,i,j,zdir,iq) = - 0.5_rp & ! 1/2
985  * ( rhokh(k,i,j) + rhokh(k+1,i,j) ) &
986  * j33g * ( phin(k+1,i,j) - phin(k,i,j) ) * rfdz(k) / gsqrt(k,i,j,i_xyw)
987  end do
988  do k = ke_pbl, ke
989  qflx_sgs_rhoq(k,i,j,zdir,iq) = 0.0_rp
990  end do
991  end do
992  end do
993 
994  end do
995 
996 
997  ! TKE
998  !$omp parallel do default(none) private(i,j,k,tke) OMP_SCHEDULE_ collapse(2) &
999  !$omp shared(JJS,JJE,IIS,IIE,KS,q,d,dt,Nu,dudz2,n2,Pr,KE_PBL)
1000  do j = jjs, jje
1001  do i = iis, iie
1002  tke = q(ks,i,j)**2 * 0.5
1003  d(ks,i,j) = tke &
1004  + dt * nu(ks,i,j) * ( dudz2(ks,i,j) - n2(ks,i,j)/pr(ks,i,j) )
1005 #ifdef MORE_HIST
1006  gen(ks,i,j) = nu(ks,i,j) * ( dudz2(ks,i,j) - n2(ks,i,j)/pr(ks,i,j) )
1007 #endif
1008  do k = ks+1, ke_pbl-1
1009  tke = q(k,i,j)**2 * 0.5
1010  d(k,i,j) = tke &
1011  + dt * nu(k,i,j) * ( dudz2(k,i,j) - n2(k,i,j)/pr(k,i,j) )
1012 #ifdef MORE_HIST
1013  gen(k,i,j) = nu(k,i,j) * ( dudz2(k,i,j) - n2(k,i,j)/pr(k,i,j) )
1014 #endif
1015  end do
1016  tke = q(ke_pbl,i,j)**2 * 0.5
1017  d(ke_pbl,i,j) = tke &
1018  + dt * nu(ke_pbl,i,j) * ( dudz2(ke_pbl,i,j) - n2(ke_pbl,i,j)/pr(ke_pbl,i,j) )
1019 #ifdef MORE_HIST
1020  gen(ke_pbl,i,j) = nu(ke_pbl,i,j) * ( dudz2(ke_pbl,i,j) - n2(ke_pbl,i,j)/pr(ke_pbl,i,j) )
1021 #endif
1022  end do
1023  end do
1024 
1025 !OCL INDEPENDENT
1026  !$omp parallel do default(none) &
1027  !$omp shared(JJS,JJE,IIS,IIE,KS,KE_PBL,dt,DENS,Nu,RFDZ,GSQRT,I_XYW,l,RCDZ,I_XYZ,I_TKE) &
1028  !$omp shared(tke_N,QTRC,KE,d,q,a,c,b,RHOQ_t) &
1029  !$omp private(i,j,k,ap) OMP_SCHEDULE_ collapse(2)
1030  do j = jjs, jje
1031  do i = iis, iie
1032  ap = - dt * 1.5_rp * ( dens(ks ,i,j)*nu(ks ,i,j) &
1033  + dens(ks+1,i,j)*nu(ks+1,i,j) ) &
1034  * rfdz(ks) / gsqrt(ks,i,j,i_xyw)
1035  a(ks,i,j) = ap * rcdz(ks) / ( dens(ks,i,j) * gsqrt(ks,i,j,i_xyz) )
1036  c(ks,i,j) = 0.0_rp
1037  b(ks,i,j) = - a(ks,i,j) + 1.0_rp + 2.0_rp * dt * q(ks,i,j) / ( b1 * l(ks,i,j) )
1038  do k = ks+1, ke_pbl-1
1039  c(k,i,j) = ap * rcdz(k) / ( dens(k,i,j) * gsqrt(k,i,j,i_xyz) )
1040  ap = - dt * 1.5_rp * ( dens(k ,i,j)*nu(k ,i,j) &
1041  + dens(k+1,i,j)*nu(k+1,i,j)) &
1042  * rfdz(k) / gsqrt(k,i,j,i_xyw)
1043  a(k,i,j) = ap * rcdz(k) / ( dens(k,i,j) * gsqrt(k,i,j,i_xyz) )
1044  b(k,i,j) = - a(k,i,j) - c(k,i,j) + 1.0_rp + 2.0_rp * dt * q(k,i,j) / ( b1 * l(k,i,j) )
1045  end do
1046 
1047  a(ke_pbl,i,j) = 0.0_rp
1048  c(ke_pbl,i,j) = ap * rcdz(ke_pbl) / ( dens(ke_pbl,i,j) * gsqrt(ke_pbl,i,j,i_xyz) )
1049  b(ke_pbl,i,j) = - c(ke_pbl,i,j) + 1.0_rp + 2.0_rp * dt * q(ke_pbl,i,j) / ( b1 * l(ke_pbl,i,j) )
1050  call diffusion_solver( &
1051  tke_n(:,i,j), & ! (out)
1052  a(:,i,j), b(:,i,j), c(:,i,j), d(:,i,j), & ! (in)
1053  ke_pbl ) ! (in)
1054  do k = ks, ke_pbl
1055  rhoq_t(k,i,j,i_tke) = ( max(tke_n(k,i,j), tke_min) - qtrc(k,i,j,i_tke) ) * dens(k,i,j) / dt
1056  end do
1057  do k = ke_pbl+1, ke
1058  rhoq_t(k,i,j,i_tke) = 0.0_rp
1059  end do
1060 
1061  end do
1062  end do
1063 
1064  end do
1065  end do
1066 
1067 #ifdef MORE_HIST
1068  gen(ke_pbl+1:ke,:,:) = 0.0_rp
1069  dudz2(ke_pbl+1:ke,:,:) = 0.0_rp
1070  l(ke_pbl+1:ke,:,:) = 0.0_rp
1071  potv(ke_pbl+1:ke,:,:) = 0.0_rp
1072  potl(ke_pbl+1:ke,:,:) = 0.0_rp
1073  call hist_in(gen, 'TKE_gen', 'generation of TKE', 'm2/s3', nohalo=.true.)
1074  call hist_in(dudz2, 'dUdZ2', 'dudz2', 'm2/s2', nohalo=.true.)
1075  call hist_in(l, 'L_mix', 'minxing length', 'm', nohalo=.true.)
1076  call hist_in(potv, 'POTV', 'virtual potential temperature', 'K', nohalo=.true.)
1077  call hist_in(potl, 'POTL', 'liquid potential temperature', 'K', nohalo=.true.)
1078 #endif
1079 
1080  return
1081  end subroutine atmos_phy_tb_mynn
1082 
1083  subroutine get_length( &
1084  l, &
1085  DENS, &
1086  q, n2, &
1087  SFLX_MU, SFLX_MV, SFLX_PT, &
1088  PT0 )
1089  use scale_grid_real, only: &
1090  fz => real_fz
1091  use scale_const, only: &
1092  grav => const_grav, &
1093  karman => const_karman, &
1094  cp => const_cpdry, &
1095  eps => const_eps
1096  implicit none
1097 
1098  real(RP), intent(out) :: l(KA,IA,JA)
1099  real(RP), intent(in) :: DENS(KA,IA,JA)
1100  real(RP), intent(in) :: q(KA,IA,JA)
1101  real(RP), intent(in) :: n2(KA,IA,JA)
1102  real(RP), intent(in) :: SFLX_MU(IA,JA)
1103  real(RP), intent(in) :: SFLX_MV(IA,JA)
1104  real(RP), intent(in) :: SFLX_PT(IA,JA)
1105  real(RP), intent(in) :: PT0(KA,IA,JA)
1106 
1107  real(RP) :: ls
1108  real(RP) :: lt
1109  real(RP) :: lb
1110  real(RP) :: rlm
1111  real(RP) :: rlt
1112 
1113  real(RP) :: qc
1114  real(RP) :: int_q
1115  real(RP) :: int_qz
1116  real(RP) :: rn2sr
1117  real(RP) :: us
1118  real(RP) :: zeta
1119 
1120  real(RP) :: z
1121  real(RP) :: qdz
1122 
1123  real(RP) :: sw
1124  integer :: k, i, j
1125 
1126 !OCL LOOP_NOFUSION,PREFETCH_SEQUENTIAL(SOFT),SWP
1127  !$omp parallel do default(none) &
1128  !$omp shared(JS,JE,IS,IE,KS,KE_PBL,q,FZ,EPS,ATMOS_PHY_TB_MYNN_Lt_MAX,SFLX_MU,SFLX_MV,DENS) &
1129  !$omp shared(SFLX_PT,PT0,n2,GRAV,l) &
1130  !$omp private(i,j,k,qdz,lt,rlt,us,rlm,qc,z,zeta,sw,ls,rn2sr,lb,int_qz,int_q) &
1131  !$omp OMP_SCHEDULE_ collapse(2)
1132  do j = js, je
1133  do i = is, ie
1134  int_qz = 0.0_rp
1135  int_q = 0.0_rp
1136  do k = ks, ke_pbl
1137  qdz = q(k,i,j) * ( fz(k,i,j) - fz(k-1,i,j) )
1138  int_qz = int_qz + ((fz(k,i,j)+fz(k-1,i,j))*0.5_rp-fz(ks-1,i,j)) * qdz
1139  int_q = int_q + qdz
1140  end do
1141  ! LT
1142  lt = min( max(0.23_rp * int_qz / (int_q + eps), &
1143  lt_min), &
1144  atmos_phy_tb_mynn_lt_max )
1145  rlt = 1.0_rp / lt
1146 
1147  us = sqrt( sqrt( sflx_mu(i,j)**2 + sflx_mv(i,j)**2 ) / dens(ks,i,j) ) ! friction velocity
1148  us = max(us, us_min)
1149  rlm = - karman * grav * sflx_pt(i,j) / (pt0(ks,i,j) * us**3 )
1150 
1151  qc = (grav/pt0(ks,i,j)*max(sflx_pt(i,j),0.0_rp)*lt)**oneoverthree
1152 
1153  do k = ks, ke_pbl
1154  z = ( fz(k,i,j)+fz(k-1,i,j) )*0.5_rp - fz(ks-1,i,j)
1155  zeta = z * rlm
1156 
1157  ! LS
1158  sw = sign(0.5_rp, zeta) + 0.5_rp ! 1 for zeta >= 0, 0 for zeta < 0
1159  ls = karman * z &
1160  * ( sw / (1.0_rp + 2.7_rp*min(zeta,1.0_rp)*sw ) &
1161  + ( (1.0_rp - 100.0_rp*zeta)*(1.0_rp-sw) )**0.2_rp )
1162 
1163  ! LB
1164  sw = sign(0.5_rp, n2(k,i,j)-eps) + 0.5_rp ! 1 for dptdz >0, 0 for dptdz <= 0
1165  rn2sr = 1.0_rp / ( sqrt(n2(k,i,j)*sw) + 1.0_rp-sw)
1166  lb = (1.0_rp + 5.0_rp * sqrt(qc*rn2sr/lt)) * q(k,i,j) * rn2sr * sw & ! qc=0 when SFLX_PT < 0
1167  + 999.e10_rp * (1.0_rp-sw)
1168 
1169  ! L
1170  l(k,i,j) = 1.0_rp / ( 1.0_rp/ls + rlt + 1.0_rp/lb )
1171  end do
1172  end do
1173  end do
1174 
1175  return
1176  end subroutine get_length
1177 
1178  subroutine get_q2_level2( &
1179  q2_2, &
1180  dudz2, Ri, l )
1181  implicit none
1182 
1183  real(RP), intent(out) :: q2_2(KA,IA,JA)
1184  real(RP), intent(in) :: dudz2(KA,IA,JA)
1185  real(RP), intent(in) :: Ri(KA,IA,JA)
1186  real(RP), intent(in) :: l(KA,IA,JA)
1187 
1188  real(RP) :: rf
1189  real(RP) :: sm_2
1190  real(RP) :: sh_2
1191 
1192  real(RP) :: q2
1193  integer :: k, i, j
1194 
1195 !OCL LOOP_NOFUSION,PREFETCH_SEQUENTIAL(SOFT),SWP
1196  !$omp parallel do default(none) &
1197  !$omp shared(JS,JE,IS,IE,KS,KE_PBL,AF12,Ri,Rf1,Rf2,Rfc,A2,G2,l,dudz2,q2_2) &
1198  !$omp private(i,j,k,rf,sh_2,sm_2,q2) OMP_SCHEDULE_ collapse(2)
1199  do j = js, je
1200  do i = is, ie
1201  do k = ks, ke_pbl
1202  rf = min(0.5_rp / af12 * ( ri(k,i,j) &
1203  + af12*rf1 &
1204  - sqrt(ri(k,i,j)**2 + 2.0_rp*af12*(rf1-2.0_rp*rf2)*ri(k,i,j) + (af12*rf1)**2) ), &
1205  rfc)
1206  sh_2 = 3.0_rp * a2 * (g1+g2) * (rfc-rf) / (1.0_rp-rf)
1207  sm_2 = sh_2 * af12 * (rf1-rf) / (rf2-rf)
1208  q2 = b1 * l(k,i,j)**2 * sm_2 * (1.0_rp-rf) * dudz2(k,i,j)
1209  q2_2(k,i,j) = max( q2, 1.e-10_rp )
1210  end do
1211  end do
1212  end do
1213 
1214  return
1215  end subroutine get_q2_level2
1216 
1217  subroutine get_smsh( &
1218  sm, sh, & ! (out)
1219  q, q2_2, & ! (in)
1220  l, n2, dudz2 ) ! (in)
1221  implicit none
1222  real(RP), intent(out) :: sm(KA,IA,JA)
1223  real(RP), intent(out) :: sh(KA,IA,JA)
1224  real(RP), intent(in) :: q(KA,IA,JA)
1225  real(RP), intent(in) :: q2_2(KA,IA,JA)
1226  real(RP), intent(in) :: l(KA,IA,JA)
1227  real(RP), intent(in) :: n2(KA,IA,JA)
1228  real(RP), intent(in) :: dudz2(KA,IA,JA)
1229 
1230  real(RP) :: l2q2
1231  real(RP) :: ac
1232  real(RP) :: ac2
1233  real(RP) :: p1
1234  real(RP) :: p2
1235  real(RP) :: p3
1236  real(RP) :: p4
1237  real(RP) :: p5
1238  real(RP) :: rd25
1239  real(RP) :: gh
1240 
1241  integer :: k, i, j
1242 
1243 !OCL LOOP_NOFUSION,PREFETCH_SEQUENTIAL(SOFT),SWP
1244  !$omp parallel do default(none) &
1245  !$omp shared(JS,JE,IS,IE,KS,KE_PBL,q,q2_2,l,n2,A2,A1,dudz2,C1,sm,sh) &
1246  !$omp private(i,j,k,ac,ac2,l2q2,gh,p1,p2,p3,p4,p5,rd25) OMP_SCHEDULE_ collapse(2)
1247  do j = js, je
1248  do i = is, ie
1249  do k = ks, ke_pbl
1250 
1251  ! level 2.5
1252  ac = min(q(k,i,j)/sqrt(q2_2(k,i,j)), 1.0_rp)
1253  ac2 = ac**2
1254  l2q2 = ( l(k,i,j) / q(k,i,j) )**2
1255  gh = - n2(k,i,j) * l2q2
1256 
1257  p1 = 1.0_rp - 3.0_rp * ac2 * a2 * b2 * (1.0_rp-c3) * gh
1258  p2 = 1.0_rp - 9.0_rp * ac2 * a1 * a2 * (1.0_rp-c2) * gh
1259  p3 = p1 + 9.0_rp * ac2 * a2**2 * (1.0_rp-c2) * (1.0_rp-c5) * gh
1260  p4 = p1 - 12.0_rp * ac2 * a1 * a2 * (1.0_rp-c2) * gh
1261  p5 = 6.0_rp * ac2 * a1**2 * dudz2(k,i,j) * l2q2
1262 
1263  rd25 = 1.0_rp / max(p2 * p4 + p5 * p3, 1.e-20_rp)
1264  sm(k,i,j) = max( ac * a1 * (p3 - 3.0_rp * c1 * p4) * rd25, 0.0_rp )
1265  sh(k,i,j) = max( ac * a2 * (p2 + 3.0_rp * c1 * p5) * rd25, 0.0_rp )
1266 
1267  end do
1268  end do
1269  end do
1270  return
1271  end subroutine get_smsh
1272 
1273  function erf(x)
1274  real(RP), parameter :: a = 0.1400122886866665_rp ! -8(pi-3)/(3pi(pi-4))
1275  real(RP), parameter :: fourpi = 1.2732395447351628_rp ! 4/pi
1276 
1277  real(RP), intent(in) :: x
1278  real(RP) :: erf
1279 
1280  real(RP) :: x2, tmp
1281 
1282  x2 = x*x
1283 
1284 #if defined(PGI) || defined(SX)
1285  tmp = min( x2 * ( fourpi + a*x2 ) / ( 1.0_rp + a*x2 ), 1.e+3_rp ) ! apply exp limiter
1286  erf = sign( sqrt( 1.0_rp - exp(-tmp) ), x )
1287 #else
1288  erf = sign( sqrt( 1.0_rp - exp(-x2 * (fourpi+a*x2)/(1.0_rp+a*x2) ) ), x )
1289 #endif
1290 
1291  return
1292  end function erf
1293 
1294 end module scale_atmos_phy_tb_mynn
scale_grid_index::is
integer, public is
start point of inner domain: x, local
Definition: scale_grid_index.F90:60
scale_debug
module DEBUG
Definition: scale_debug.F90:13
scale_grid_index::je
integer, public je
end point of inner domain: y, local
Definition: scale_grid_index.F90:63
scale_const::const_cpdry
real(rp), public const_cpdry
specific heat (dry air,constant pressure) [J/kg/K]
Definition: scale_const.F90:58
scale_atmos_saturation
module ATMOSPHERE / Saturation adjustment
Definition: scale_atmos_sub_saturation.F90:16
scale_process::prc_mpistop
subroutine, public prc_mpistop
Abort MPI.
Definition: scale_process.F90:267
scale_atmos_phy_tb_common::atmos_phy_tb_diffusion_solver
subroutine, public atmos_phy_tb_diffusion_solver(phi, a, b, c, d, KE_TB)
Definition: scale_atmos_phy_tb_common.F90:1459
scale_atmos_hydrometeor::i_qc
integer, public i_qc
Definition: scale_atmos_sub_hydrometeor.F90:79
scale_tracer::tracer_r
real(rp), dimension(qa_max), public tracer_r
Definition: scale_tracer.F90:41
scale_grid_index::iblock
integer, public iblock
block size for cache blocking: x
Definition: scale_grid_index.F90:47
scale_stdio::io_l
logical, public io_l
output log or not? (this process)
Definition: scale_stdio.F90:61
scale_grid_index::zdir
integer, parameter, public zdir
Definition: scale_grid_index.F90:35
scale_atmos_phy_tb_mynn::get_length
subroutine get_length(l, DENS, q, n2, SFLX_MU, SFLX_MV, SFLX_PT, PT0)
Definition: scale_atmos_phy_tb_mynn.F90:1089
scale_stdio
module STDIO
Definition: scale_stdio.F90:12
scale_grid_index::ydir
integer, parameter, public ydir
Definition: scale_grid_index.F90:37
scale_grid_index::ke
integer, public ke
end point of inner domain: z, local
Definition: scale_grid_index.F90:59
scale_grid_index::xdir
integer, parameter, public xdir
Definition: scale_grid_index.F90:36
scale_tracer::qa
integer, public qa
Definition: scale_tracer.F90:35
scale_grid_real::real_fz
real(rp), dimension(:,:,:), allocatable, public real_fz
geopotential height [m] (cell face )
Definition: scale_grid_real.F90:41
scale_atmos_phy_tb_mynn::get_q2_level2
subroutine get_q2_level2(q2_2, dudz2, Ri, l)
Definition: scale_atmos_phy_tb_mynn.F90:1181
scale_tracer::tracer_advc
logical, dimension(qa_max), public tracer_advc
Definition: scale_tracer.F90:43
scale_grid_real::real_cz
real(rp), dimension(:,:,:), allocatable, public real_cz
geopotential height [m] (cell center)
Definition: scale_grid_real.F90:40
scale_atmos_phy_tb_mynn::atmos_phy_tb_mynn_setup
subroutine, public atmos_phy_tb_mynn_setup(CDZ, CDX, CDY, CZ)
Setup.
Definition: scale_atmos_phy_tb_mynn.F90:143
scale_debug::check
subroutine, public check(current_line, v)
Undefined value checker.
Definition: scale_debug.F90:58
scale_grid::grid_rcdz
real(rp), dimension(:), allocatable, public grid_rcdz
reciprocal of center-dz
Definition: scale_grid_cartesian.F90:62
scale_const::const_karman
real(rp), parameter, public const_karman
von Karman constant
Definition: scale_const.F90:52
scale_tracer::tracer_cv
real(rp), dimension(qa_max), public tracer_cv
Definition: scale_tracer.F90:40
scale_gridtrans::i_xvw
integer, public i_xvw
Definition: scale_gridtrans.F90:52
scale_const::const_rdry
real(rp), public const_rdry
specific gas constant (dry air) [J/kg/K]
Definition: scale_const.F90:57
scale_const::const_undef
real(rp), public const_undef
Definition: scale_const.F90:43
scale_grid_index
module grid index
Definition: scale_grid_index.F90:14
scale_stdio::io_nml
logical, public io_nml
output log or not? (for namelist, this process)
Definition: scale_stdio.F90:62
scale_tracer
module TRACER
Definition: scale_tracer.F90:15
scale_grid_index::ia
integer, public ia
of whole cells: x, local, with HALO
Definition: scale_grid_index.F90:55
scale_atmos_hydrometeor
module Hydrometeor
Definition: scale_atmos_sub_hydrometeor.F90:14
scale_gridtrans
module GRIDTRANS
Definition: scale_gridtrans.F90:15
scale_grid_real
module GRID (real space)
Definition: scale_grid_real.F90:14
scale_grid_index::ka
integer, public ka
of whole cells: z, local, with HALO
Definition: scale_grid_index.F90:54
scale_const::const_lhv0
real(rp), parameter, public const_lhv0
latent heat of vaporizaion at 0C [J/kg]
Definition: scale_const.F90:77
scale_grid_index::jblock
integer, public jblock
block size for cache blocking: y
Definition: scale_grid_index.F90:48
scale_gridtrans::i_xyw
integer, public i_xyw
Definition: scale_gridtrans.F90:50
scale_atmos_phy_tb_mynn::atmos_phy_tb_mynn
subroutine, public atmos_phy_tb_mynn(qflx_sgs_momz, qflx_sgs_momx, qflx_sgs_momy, qflx_sgs_rhot, qflx_sgs_rhoq, RHOQ_t, Nu, Ri, Pr, MOMZ, MOMX, MOMY, RHOT, DENS, QTRC, N2_in, SFLX_MW, SFLX_MU, SFLX_MV, SFLX_SH, SFLX_Q, GSQRT, J13G, J23G, J33G, MAPF, dt)
Definition: scale_atmos_phy_tb_mynn.F90:224
scale_comm
module COMMUNICATION
Definition: scale_comm.F90:23
scale_const::const_grav
real(rp), public const_grav
standard acceleration of gravity [m/s2]
Definition: scale_const.F90:48
scale_grid_index::js
integer, public js
start point of inner domain: y, local
Definition: scale_grid_index.F90:62
scale_const::const_undef2
integer, parameter, public const_undef2
undefined value (INT2)
Definition: scale_const.F90:40
scale_atmos_phy_tb_mynn::atmos_phy_tb_mynn_config
subroutine, public atmos_phy_tb_mynn_config(TYPE_TB, I_TKE_out)
Config.
Definition: scale_atmos_phy_tb_mynn.F90:109
scale_process
module PROCESS
Definition: scale_process.F90:14
scale_gridtrans::i_uyw
integer, public i_uyw
Definition: scale_gridtrans.F90:51
scale_const::const_epstvap
real(rp), public const_epstvap
1 / epsilon - 1
Definition: scale_const.F90:72
scale_const::const_rvap
real(rp), parameter, public const_rvap
specific gas constant (water vapor) [J/kg/K]
Definition: scale_const.F90:65
scale_const
module CONSTANT
Definition: scale_const.F90:14
scale_grid_index::ks
integer, public ks
start point of inner domain: z, local
Definition: scale_grid_index.F90:58
scale_atmos_phy_tb_common
module ATMOSPHERE / Physics Turbulence
Definition: scale_atmos_phy_tb_common.F90:15
scale_grid
module GRID (cartesian)
Definition: scale_grid_cartesian.F90:17
scale_atmos_hydrometeor::i_qi
integer, public i_qi
Definition: scale_atmos_sub_hydrometeor.F90:81
scale_prof
module profiler
Definition: scale_prof.F90:10
scale_grid_index::ie
integer, public ie
end point of inner domain: x, local
Definition: scale_grid_index.F90:61
scale_const::const_eps
real(rp), public const_eps
small number
Definition: scale_const.F90:36
scale_atmos_thermodyn
module ATMOSPHERE / Thermodynamics
Definition: scale_atmos_sub_thermodyn.F90:19
scale_atmos_hydrometeor::i_qv
integer, public i_qv
Definition: scale_atmos_sub_hydrometeor.F90:68
scale_precision
module PRECISION
Definition: scale_precision.F90:13
scale_history
module HISTORY
Definition: scale_history.F90:16
scale_const::const_pi
real(rp), public const_pi
pi
Definition: scale_const.F90:34
scale_grid::grid_rfdz
real(rp), dimension(:), allocatable, public grid_rfdz
reciprocal of face-dz
Definition: scale_grid_cartesian.F90:72
scale_stdio::io_fid_conf
integer, public io_fid_conf
Config file ID.
Definition: scale_stdio.F90:55
scale_tracer::tracer_regist
subroutine, public tracer_regist(QS, NQ, NAME, DESC, UNIT, CV, CP, R, ADVC, MASS)
Regist tracer.
Definition: scale_tracer.F90:64
scale_gridtrans::i_xyz
integer, public i_xyz
Definition: scale_gridtrans.F90:49
scale_stdio::io_fid_log
integer, public io_fid_log
Log file ID.
Definition: scale_stdio.F90:56
scale_const::const_cpvap
real(rp), parameter, public const_cpvap
specific heat (water vapor, constant pressure) [J/kg/K]
Definition: scale_const.F90:66
scale_atmos_phy_tb_mynn
module ATMOSPHERE / Physics Turbulence
Definition: scale_atmos_phy_tb_mynn.F90:21
scale_stdio::io_fid_nml
integer, public io_fid_nml
Log file ID (only for output namelist)
Definition: scale_stdio.F90:57
scale_tracer::tracer_mass
real(rp), dimension(qa_max), public tracer_mass
Definition: scale_tracer.F90:42
scale_grid_index::ja
integer, public ja
of whole cells: y, local, with HALO
Definition: scale_grid_index.F90:56
Generated by   doxygen 1.8.14