SCALE-RM
Functions/Subroutines
scale_atmos_phy_tb_mynn Module Reference

module ATMOSPHERE / Physics Turbulence More...

Functions/Subroutines

subroutine, public atmos_phy_tb_mynn_setup (TYPE_TB, CDZ, CDX, CDY, CZ)
 
subroutine, public atmos_phy_tb_mynn (qflx_sgs_momz, qflx_sgs_momx, qflx_sgs_momy, qflx_sgs_rhot, qflx_sgs_rhoq, tke, tke_t, Nu, Ri, Pr, N2, MOMZ, MOMX, MOMY, RHOT, DENS, QTRC, SFLX_MW, SFLX_MU, SFLX_MV, SFLX_SH, SFLX_QV, GSQRT, J13G, J23G, J33G, MAPF, dt)
 
subroutine get_length (l, DENS, q, n2, SFLX_MU, SFLX_MV, SFLX_PT, PT0)
 
subroutine get_q2_level2 (q2_2, dudz2, Ri, l)
 

Detailed Description

module ATMOSPHERE / Physics Turbulence

Description
Boundary layer turbulence model Mellor-Yamada Nakanishi-Niino model
Author
Team SCALE
History
  • 2014-08-27 (S.Nishizawa) [new]
  • Reference
    • Nakanishi and Niino, 2009: Development of an improved turbulence closure model for the atmospheric boundary layer. J. Meteorol. Soc. Japan, 87, 895-912
NAMELIST
  • PARAM_ATMOS_PHY_TB_MYNN
    nametypedefault valuecomment
    ATMOS_PHY_TB_MYNN_TKE_INIT logical .false. set tke with that of level 2 at the first time if .true.
    ATMOS_PHY_TB_MYNN_PBL_MAX real(RP) 1e10_RP maximum height of the PBL
    ATMOS_PHY_TB_MYNN_N2_MAX real(RP) 1.E3_RP
    ATMOS_PHY_TB_MYNN_NU_MIN real(RP) 1.E-6_RP
    ATMOS_PHY_TB_MYNN_NU_MAX real(RP) 10000._RP
    ATMOS_PHY_TB_MYNN_KH_MIN real(RP) 1.E-6_RP
    ATMOS_PHY_TB_MYNN_KH_MAX real(RP) 10000._RP
    ATMOS_PHY_TB_MYNN_LT_MAX real(RP) 700._RP ~ 0.23 * 3 km

History Output
namedescriptionunitvariable
L_mix minxing length m l
POTL liquid potential temperature K POTL
POTV virtual potential temperature K POTV
TKE_advc advection of TKE m2/s3 adv
TKE_gen generation of TKE m2/s3 gen
dUdZ2 dudz2 m2/s2 dudz2

Function/Subroutine Documentation

◆ atmos_phy_tb_mynn_setup()

subroutine, public scale_atmos_phy_tb_mynn::atmos_phy_tb_mynn_setup ( character(len=*), intent(in)  TYPE_TB,
real(rp), dimension(ka), intent(in)  CDZ,
real(rp), dimension(ia), intent(in)  CDX,
real(rp), dimension(ja), intent(in)  CDY,
real(rp), dimension (ka,ia,ja), intent(in)  CZ 
)

Definition at line 103 of file scale_atmos_phy_tb_mynn.F90.

References scale_const::const_pi, scale_grid_index::ie, scale_stdio::io_fid_conf, scale_stdio::io_fid_log, scale_stdio::io_l, scale_stdio::io_lnml, scale_grid_index::is, scale_grid_index::je, scale_grid_index::js, scale_grid_index::ke, scale_grid_index::ks, and scale_process::prc_mpistop().

Referenced by scale_atmos_phy_tb_hybrid::atmos_phy_tb_hybrid_setup(), and scale_atmos_phy_tb::atmos_phy_tb_setup().

103  use scale_process, only: &
104  prc_mpistop
105  use scale_const, only: &
106  pi => const_pi
107  implicit none
108 
109  character(len=*), intent(in) :: type_tb
110 
111  real(RP), intent(in) :: cdz(ka)
112  real(RP), intent(in) :: cdx(ia)
113  real(RP), intent(in) :: cdy(ja)
114  real(RP), intent(in) :: cz (ka,ia,ja)
115 
116  real(RP) :: atmos_phy_tb_mynn_pbl_max = 1e10_rp
117 
118  namelist / param_atmos_phy_tb_mynn / &
119  atmos_phy_tb_mynn_tke_init, &
120  atmos_phy_tb_mynn_pbl_max, &
121  atmos_phy_tb_mynn_n2_max, &
122  atmos_phy_tb_mynn_nu_min, &
123  atmos_phy_tb_mynn_nu_max, &
124  atmos_phy_tb_mynn_kh_min, &
125  atmos_phy_tb_mynn_kh_max, &
126  atmos_phy_tb_mynn_lt_max
127 
128  integer :: k, i, j
129  integer :: ierr
130  !---------------------------------------------------------------------------
131 
132  if( io_l ) write(io_fid_log,*)
133  if( io_l ) write(io_fid_log,*) '++++++ Module[TURBULENCE] / Categ[ATMOS PHYSICS] / Origin[SCALElib]'
134  if( io_l ) write(io_fid_log,*) '+++ Mellor-Yamada Nakanishi-Niino Model'
135 
136  if ( type_tb /= 'MYNN' ) then
137  write(*,*) 'xxx ATMOS_PHY_TB_TYPE is not MYNN. Check!'
138  call prc_mpistop
139  endif
140 
141 
142  !--- read namelist
143  rewind(io_fid_conf)
144  read(io_fid_conf,nml=param_atmos_phy_tb_mynn,iostat=ierr)
145  if( ierr < 0 ) then !--- missing
146  if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
147  elseif( ierr > 0 ) then !--- fatal error
148  write(*,*) 'xxx Not appropriate names in namelist PARAM_ATMOS_PHY_TB_MYNN. Check!'
149  call prc_mpistop
150  endif
151  if( io_lnml ) write(io_fid_log,nml=param_atmos_phy_tb_mynn)
152 
153  do k = ks, ke-1
154  do j = js, je
155  do i = is, ie
156  if ( atmos_phy_tb_mynn_pbl_max >= cz(k,i,j) ) then
157  ke_pbl = k
158  end if
159  end do
160  end do
161  end do
162 
163  a1 = b1 * (1.0_rp - 3.0_rp * g1) / 6.0_rp
164  a2 = 1.0_rp / (3.0_rp * g1 * b1**(1.0_rp/3.0_rp) * prn )
165  c1 = g1 - 1.0_rp / ( 3.0_rp * a1 * b1**(1.0_rp/3.0_rp) )
166  g2 = ( 2.0_rp * a1 * (3.0_rp - 2.0_rp * c2) + b2 * (1.0_rp - c3) ) / b1
167  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)
168  f2 = b1 * (g1 + g2) - 3.0_rp * a1 * (1.0_rp - c2)
169 
170  rf1 = b1 * (g1 - c1) / f1
171  rf2 = b1 * g1 / f2
172  rfc = g1 / (g1 + g2)
173 
174  af12 = a1 * f1 / ( a2 * f2 )
175 
176  sqrt_2pi = sqrt( 2.0_rp * pi )
177  rsqrt_2pi = 1.0_rp / sqrt_2pi
178  rsqrt_2 = 1.0_rp / sqrt( 2.0_rp )
179 
180  return
scale_grid_index::is
integer, public is
start point of inner domain: x, local
Definition: scale_grid_index.F90:60
scale_grid_index::je
integer, public je
end point of inner domain: y, local
Definition: scale_grid_index.F90:63
scale_process::prc_mpistop
subroutine, public prc_mpistop
Abort MPI.
Definition: scale_process.F90:234
scale_grid_index::ke
integer, public ke
end point of inner domain: z, local
Definition: scale_grid_index.F90:59
scale_grid_index::ia
integer, public ia
of x whole cells (local, with HALO)
Definition: scale_grid_index.F90:55
scale_grid_index::ka
integer, public ka
of z whole cells (local, with HALO)
Definition: scale_grid_index.F90:54
scale_grid_index::js
integer, public js
start point of inner domain: y, local
Definition: scale_grid_index.F90:62
scale_process
module PROCESS
Definition: scale_process.F90:14
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_grid_index::ie
integer, public ie
end point of inner domain: x, local
Definition: scale_grid_index.F90:61
scale_const::const_pi
real(rp), public const_pi
pi
Definition: scale_const.F90:34
scale_grid_index::ja
integer, public ja
of y whole cells (local, with HALO)
Definition: scale_grid_index.F90:56
Here is the call graph for this function:
Here is the caller graph for this function:

◆ atmos_phy_tb_mynn()

subroutine, public scale_atmos_phy_tb_mynn::atmos_phy_tb_mynn ( real(rp), dimension(ka,ia,ja,3), intent(out)  qflx_sgs_momz,
real(rp), dimension(ka,ia,ja,3), intent(out)  qflx_sgs_momx,
real(rp), dimension(ka,ia,ja,3), intent(out)  qflx_sgs_momy,
real(rp), dimension(ka,ia,ja,3), intent(out)  qflx_sgs_rhot,
real(rp), dimension(ka,ia,ja,3,qa), intent(out)  qflx_sgs_rhoq,
real(rp), dimension (ka,ia,ja), intent(inout)  tke,
real(rp), dimension(ka,ia,ja), intent(out)  tke_t,
real(rp), dimension(ka,ia,ja), intent(out)  Nu,
real(rp), dimension(ka,ia,ja), intent(out)  Ri,
real(rp), dimension(ka,ia,ja), intent(out)  Pr,
real(rp), dimension(ka,ia,ja), intent(out)  N2,
real(rp), dimension(ka,ia,ja), intent(in)  MOMZ,
real(rp), dimension(ka,ia,ja), intent(in)  MOMX,
real(rp), dimension(ka,ia,ja), intent(in)  MOMY,
real(rp), dimension(ka,ia,ja), intent(in)  RHOT,
real(rp), dimension(ka,ia,ja), intent(in)  DENS,
real(rp), dimension(ka,ia,ja,qa), intent(in)  QTRC,
real(rp), dimension(ia,ja), intent(in)  SFLX_MW,
real(rp), dimension(ia,ja), intent(in)  SFLX_MU,
real(rp), dimension(ia,ja), intent(in)  SFLX_MV,
real(rp), dimension(ia,ja), intent(in)  SFLX_SH,
real(rp), dimension(ia,ja), intent(in)  SFLX_QV,
real(rp), dimension (ka,ia,ja,7), intent(in)  GSQRT,
real(rp), dimension (ka,ia,ja,7), intent(in)  J13G,
real(rp), dimension (ka,ia,ja,7), intent(in)  J23G,
real(rp), intent(in)  J33G,
real(rp), dimension (ia,ja,2,4), intent(in)  MAPF,
real(dp), intent(in)  dt 
)
Parameters
[in]gsqrtvertical metrics {G}^1/2
[in]j13g(1,3) element of Jacobian matrix
[in]j23g(1,3) element of Jacobian matrix
[in]j33g(3,3) element of Jacobian matrix
[in]mapfmap factor

Definition at line 192 of file scale_atmos_phy_tb_mynn.F90.

References scale_atmos_phy_tb_common::atmos_phy_tb_diffusion_solver(), scale_const::const_cpdry, scale_const::const_epstvap, scale_const::const_grav, scale_const::const_rdry, scale_const::const_rvap, get_length(), get_q2_level2(), scale_grid::grid_rcdx, scale_grid::grid_rcdy, scale_grid::grid_rcdz, scale_grid::grid_rfdz, scale_tracer::i_qc, scale_tracer::i_qi, scale_tracer::i_qv, scale_gridtrans::i_uy, scale_gridtrans::i_uyw, scale_gridtrans::i_uyz, scale_gridtrans::i_xv, scale_gridtrans::i_xvw, scale_gridtrans::i_xvz, scale_gridtrans::i_xy, scale_gridtrans::i_xyw, scale_gridtrans::i_xyz, scale_grid_index::iblock, scale_grid_index::ie, scale_stdio::io_fid_log, scale_stdio::io_l, scale_grid_index::is, scale_grid_index::jblock, scale_grid_index::je, scale_grid_index::js, scale_grid_index::ka, scale_grid_index::ke, scale_grid_index::ks, scale_tracer::qa, scale_tracer::qqe, scale_tracer::qqs, scale_grid_real::real_cz, scale_grid_index::xdir, scale_grid_index::ydir, and scale_grid_index::zdir.

Referenced by scale_atmos_phy_tb_hybrid::atmos_phy_tb_hybrid_setup(), and scale_atmos_phy_tb::atmos_phy_tb_setup().

192  use scale_precision
193  use scale_grid_index
194  use scale_tracer
195  use scale_const, only: &
196  grav => const_grav, &
197  r => const_rdry, &
198  rvap => const_rvap, &
199  cp => const_cpdry, &
200  epstvap => const_epstvap
201  use scale_grid, only: &
202  rcdz => grid_rcdz, &
203  rfdz => grid_rfdz, &
204  rcdx => grid_rcdx, &
205  rcdy => grid_rcdy
206  use scale_grid_real, only: &
207  cz => real_cz
208  use scale_gridtrans, only: &
209  i_xyz, &
210  i_xyw, &
211  i_xvw, &
212  i_uyw, &
213  i_uyz, &
214  i_xvz, &
215  i_xy, &
216  i_uy, &
217  i_xv
218  use scale_comm, only: &
219  comm_vars8, &
220  comm_vars, &
221  comm_wait
222  use scale_atmos_phy_tb_common, only: &
223  diffusion_solver => atmos_phy_tb_diffusion_solver
224  use scale_atmos_thermodyn, only: &
225  atmos_thermodyn_temp_pres, &
226  atmos_thermodyn_templhv, &
227  atmos_thermodyn_templhs
228  use scale_atmos_saturation, only: &
229  atmos_saturation_alpha, &
230  atmos_saturation_pres2qsat => atmos_saturation_pres2qsat_all
231 #ifdef MORE_HIST
232  use scale_history, only: &
233  hist_in
234 #endif
235  implicit none
236 
237  real(RP), intent(out) :: qflx_sgs_momz(ka,ia,ja,3)
238  real(RP), intent(out) :: qflx_sgs_momx(ka,ia,ja,3)
239  real(RP), intent(out) :: qflx_sgs_momy(ka,ia,ja,3)
240  real(RP), intent(out) :: qflx_sgs_rhot(ka,ia,ja,3)
241  real(RP), intent(out) :: qflx_sgs_rhoq(ka,ia,ja,3,qa)
242 
243  real(RP), intent(inout) :: tke (ka,ia,ja) ! TKE
244  real(RP), intent(out) :: tke_t(ka,ia,ja) ! tendency of TKE
245  real(RP), intent(out) :: nu(ka,ia,ja) ! eddy viscosity (center)
246  real(RP), intent(out) :: ri(ka,ia,ja) ! Richardson number
247  real(RP), intent(out) :: pr(ka,ia,ja) ! Plandtle number
248  real(RP), intent(out) :: n2(ka,ia,ja) ! squared Brunt-Vaisala frequency
249 
250  real(RP), intent(in) :: momz(ka,ia,ja)
251  real(RP), intent(in) :: momx(ka,ia,ja)
252  real(RP), intent(in) :: momy(ka,ia,ja)
253  real(RP), intent(in) :: rhot(ka,ia,ja)
254  real(RP), intent(in) :: dens(ka,ia,ja)
255  real(RP), intent(in) :: qtrc(ka,ia,ja,qa)
256 
257  real(RP), intent(in) :: sflx_mw(ia,ja)
258  real(RP), intent(in) :: sflx_mu(ia,ja)
259  real(RP), intent(in) :: sflx_mv(ia,ja)
260  real(RP), intent(in) :: sflx_sh(ia,ja)
261  real(RP), intent(in) :: sflx_qv(ia,ja)
262 
263  real(RP), intent(in) :: gsqrt (ka,ia,ja,7)
264  real(RP), intent(in) :: j13g (ka,ia,ja,7)
265  real(RP), intent(in) :: j23g (ka,ia,ja,7)
266  real(RP), intent(in) :: j33g
267  real(RP), intent(in) :: mapf (ia,ja,2,4)
268  real(DP), intent(in) :: dt
269 
270  real(RP) :: u(ka,ia,ja)
271  real(RP) :: v(ka,ia,ja)
272  real(RP) :: phin(ka,ia,ja)
273 
274 
275  real(RP) :: sm(ka,ia,ja)
276  real(RP) :: sh(ka,ia,ja)
277  real(RP) :: l(ka,ia,ja)
278  real(RP) :: q(ka,ia,ja)
279  real(RP) :: dudz2(ka,ia,ja)
280  real(RP) :: q2_2(ka,ia,ja)
281  real(RP) :: kh
282  real(RP) :: rhokh (ka,ia,ja)
283 
284  real(RP) :: sflx_pt(ia,ja) ! surface potential temperature flux
285 
286  real(RP) :: a(ka,ia,ja)
287  real(RP) :: b(ka,ia,ja)
288  real(RP) :: c(ka,ia,ja)
289  real(RP) :: d(ka,ia,ja)
290  real(RP) :: ap
291  real(RP) :: tke_n(ka,ia,ja)
292 
293  real(RP) :: pott(ka,ia,ja)
294  real(RP) :: potv(ka,ia,ja)
295  real(RP) :: potl(ka,ia,ja)
296  real(RP) :: teml(ka,ia,ja)
297 
298  real(RP) :: qw(ka,ia,ja)
299  real(RP) :: ql
300  real(RP) :: qs
301  real(RP) :: qdry
302 
303  real(RP) :: temp(ka,ia,ja)
304  real(RP) :: pres(ka,ia,ja)
305 
306  real(RP) :: lh(ka,ia,ja)
307  real(RP) :: lhv(ka,ia,ja)
308  real(RP) :: lhs(ka,ia,ja)
309  real(RP) :: alpha(ka,ia,ja)
310 
311  real(RP) :: ac
312 
313  real(RP) :: q1
314  real(RP) :: qsl(ka,ia,ja)
315  real(RP) :: dqsl
316  real(RP) :: sigma_s
317  real(RP) :: rr
318  real(RP) :: rt
319  real(RP) :: betat
320  real(RP) :: betaq
321  real(RP) :: aa
322  real(RP) :: bb
323  real(RP) :: cc
324 
325  real(RP) :: flux_z(ka,ia,ja)
326  real(RP) :: flux_x(ka,ia,ja)
327  real(RP) :: flux_y(ka,ia,ja)
328  real(RP) :: mflx
329  real(RP) :: advc
330 
331 #ifdef MORE_HIST
332  real(RP) :: adv(ka,ia,ja)
333  real(RP) :: gen(ka,ia,ja)
334 #endif
335 
336  integer :: k, i, j, iq
337  integer :: iis, iie, jjs, jje
338 
339  if ( io_l ) write(io_fid_log, *) "*** Physics step: Turbulence (MYNN)"
340 
341 #ifdef DEBUG
342  pott(:,:,:) = undef
343  potl(:,:,:) = undef
344  temp(:,:,:) = undef
345 #endif
346 
347 !OCL XFILL
348  do j = js , je
349  do i = is-1, ie
350  do k = ks-1, ke
351  qflx_sgs_momz(k,i,j,xdir) = 0.0_rp
352  end do
353  end do
354  end do
355 !OCL XFILL
356  do j = js , je
357  do i = is , ie+1
358  do k = ks-1, ke+1
359  qflx_sgs_momx(k,i,j,xdir) = 0.0_rp
360  end do
361  end do
362  end do
363 
364 !OCL XFILL
365  do j = js , je
366  do i = is-1, ie
367  do k = ks-1, ke+1
368  qflx_sgs_momy(k,i,j,xdir) = 0.0_rp
369  end do
370  end do
371  end do
372 !OCL XFILL
373  do j = js , je
374  do i = is-1, ie
375  do k = ks-1, ke+1
376  qflx_sgs_rhot(k,i,j,xdir) = 0.0_rp
377  end do
378  end do
379  end do
380 !OCL XFILL
381  do iq = 1, qa
382  do j = js , je
383  do i = is-1, ie
384  do k = ks-1, ke+1
385  qflx_sgs_rhoq(k,i,j,xdir,iq) = 0.0_rp
386  end do
387  end do
388  end do
389  end do
390 
391 !OCL XFILL
392  do j = js-1, je
393  do i = is , ie
394  do k = ks-1, ke
395  qflx_sgs_momz(k,i,j,ydir) = 0.0_rp
396  end do
397  end do
398  end do
399 !OCL XFILL
400  do j = js-1, je
401  do i = is , ie
402  do k = ks-1, ke+1
403  qflx_sgs_momx(k,i,j,ydir) = 0.0_rp
404  end do
405  end do
406  end do
407 !OCL XFILL
408  do j = js , je+1
409  do i = is , ie
410  do k = ks-1, ke+1
411  qflx_sgs_momy(k,i,j,ydir) = 0.0_rp
412  end do
413  end do
414  end do
415 !OCL XFILL
416  do j = js-1, je
417  do i = is , ie
418  do k = ks-1, ke+1
419  qflx_sgs_rhot(k,i,j,ydir) = 0.0_rp
420  end do
421  end do
422  end do
423 !OCL XFILL
424  do iq = 1, qa
425  do j = js-1, je
426  do i = is , ie
427  do k = ks-1, ke+1
428  qflx_sgs_rhoq(k,i,j,ydir,iq) = 0.0_rp
429  end do
430  end do
431  end do
432  end do
433 !OCL XFILL
434  do j = js, je
435  do i = is, ie
436  do k = ks-1, ke+1
437  qflx_sgs_momz(k,i,j,zdir) = 0.0_rp
438  end do
439  end do
440  end do
441 !OCL XFILL
442  do iq = i_qv+1, qa
443  do j = js-1, je
444  do i = is , ie
445  do k = ks-1, ke+1
446  qflx_sgs_rhoq(k,i,j,zdir,iq) = 0.0_rp
447  end do
448  end do
449  end do
450  end do
451 
452 
453 
454  do jjs = js, je, jblock
455  jje = jjs+jblock-1
456  do iis = is, ie, iblock
457  iie = iis+iblock-1
458 
459 !OCL XFILL
460  do j = jjs , jje+1
461  do i = iis-1, iie+1
462  do k = ks, ke_pbl+1
463  u(k,i,j) = 0.5_rp * ( momx(k,i,j) + momx(k,i-1,j) ) / dens(k,i,j)
464  end do
465  end do
466  end do
467 
468 !OCL XFILL
469  do j = jjs-1, jje+1
470  do i = iis , iie+1
471  do k = ks, ke_pbl+1
472  v(k,i,j) = 0.5_rp * ( momy(k,i,j) + momy(k,i,j-1) ) / dens(k,i,j)
473  end do
474  end do
475  end do
476 
477  end do
478  end do
479 
480  call atmos_thermodyn_temp_pres( temp, pres, & ! (out)
481  dens, rhot, qtrc ) ! (in)
482 
483 
484  call atmos_thermodyn_templhv( lhv, temp )
485  call atmos_thermodyn_templhs( lhs, temp )
486  call atmos_saturation_alpha( alpha, temp )
487 
488 !OCL LOOP_NOFUSION,PREFETCH_SEQUENTIAL(SOFT),SWP
489  do j = js, je
490  do i = is, ie
491  do k = ks, ke_pbl+1
492  ql = 0.0_rp
493  if ( i_qc > 0 ) ql = qtrc(k,i,j,i_qc)
494 ! do iq = QWS, QWE
495 ! ql = ql + QTRC(k,i,j,iq)
496 ! end do
497  qs = 0.0_rp
498  if ( i_qi > 0 ) qs = qtrc(k,i,j,i_qi)
499 ! do iq = QIS, QIE
500 ! qs = qs + QTRC(k,i,j,iq)
501 ! end do
502  qdry = 1.0_rp
503  do iq = qqs, qqe
504  qdry = qdry - qtrc(k,i,j,iq)
505  end do
506 
507  qw(k,i,j) = qtrc(k,i,j,i_qv) + ql + qs
508 
509  lh(k,i,j) = lhv(k,i,j) * alpha(k,i,j) + lhs(k,i,j) * ( 1.0_rp-alpha(k,i,j) )
510 
511  pott(k,i,j) = rhot(k,i,j) / dens(k,i,j)
512  ! liquid water potential temperature
513  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) * cp ) )
514  teml(k,i,j) = potl(k,i,j) * temp(k,i,j) / pott(k,i,j)
515 
516  ! virtual potential temperature for derivertive
517 ! POTV(k,i,j) = ( 1.0_RP + EPSTvap * Qw(k,i,j) ) * POTL(k,i,j)
518  potv(k,i,j) = ( qdry + (epstvap+1.0_rp) * qtrc(k,i,j,i_qv) ) * potl(k,i,j)
519 
520  end do
521 
522  sflx_pt(i,j) = sflx_sh(i,j) / ( cp * dens(ks,i,j) ) &
523  * pott(ks,i,j) / temp(ks,i,j)
524 
525  n2(ks,i,j) = min(atmos_phy_tb_mynn_n2_max, &
526  grav * ( potv(ks+1,i,j) - potv(ks,i,j) ) &
527  / ( ( cz(ks+1,i,j)-cz(ks,i,j) ) * potv(ks,i,j) ) )
528  do k = ks+1, ke_pbl
529  n2(k,i,j) = min(atmos_phy_tb_mynn_n2_max, &
530  grav * ( potv(k+1,i,j) - potv(k-1,i,j) ) &
531  / ( ( cz(k+1,i,j)-cz(k-1,i,j) ) * potv(k,i,j) ) )
532  end do
533 
534  dudz2(ks,i,j) = ( ( u(ks+1,i,j) - u(ks,i,j) )**2 + ( v(ks+1,i,j) - v(ks,i,j) )**2 ) &
535  / ( cz(ks+1,i,j) - cz(ks,i,j) )**2
536  do k = ks+1, ke_pbl
537  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 ) &
538  / ( cz(k+1,i,j) - cz(k-1,i,j) )**2
539  end do
540 
541  do k = ks, ke_pbl
542  ri(k,i,j) = n2(k,i,j) / max(dudz2(k,i,j), 1e-10_rp)
543  end do
544 
545  end do
546  end do
547 
548  if ( atmos_phy_tb_mynn_tke_init ) then
549 !OCL XFILL
550  do j = js, je
551  do i = is, ie
552  do k = ks, ke_pbl
553  q(k,i,j) = sqrt( tke_min*2.0_rp )
554  end do
555  end do
556  end do
557  else
558 !OCL XFILL
559  do j = js, je
560  do i = is, ie
561  do k = ks, ke_pbl
562  q(k,i,j) = sqrt( max(tke(k,i,j), tke_min)*2.0_rp )
563  end do
564  end do
565  end do
566  end if
567 
568  ! length
569  call get_length( &
570  l, & ! (out)
571  dens, & ! (in)
572  q, n2, & ! (in)
573  sflx_mu, sflx_mv, sflx_pt, & ! (in)
574  pott ) ! (in)
575 
576  call get_q2_level2( &
577  q2_2, & ! (out)
578  dudz2, ri, l ) ! (in)
579 
580  if ( atmos_phy_tb_mynn_tke_init ) then
581 !OCL XFILL
582  do j = js, je
583  do i = is, ie
584  do k = ks, ke_pbl
585  tke(k,i,j) = max(q2_2(k,i,j) * 0.5_rp, tke_min)
586  q(k,i,j) = sqrt( tke(k,i,j) * 2.0_rp )
587  end do
588  end do
589  end do
590 !OCL XFILL
591  do j = js, je
592  do i = is, ie
593  do k = ke_pbl+1, ke
594  tke(k,i,j) = 0.0_rp
595  end do
596  end do
597  end do
598  call comm_vars8(tke, 1)
599  call comm_wait (tke, 1)
600 
601  atmos_phy_tb_mynn_tke_init = .false.
602  end if
603 
604  call get_smsh( sm, sh, & ! (out)
605  q, q2_2, & ! (in)
606  l, n2, dudz2 ) ! (in)
607 
608  call atmos_saturation_pres2qsat( qsl(ks:ke_pbl,:,:), & ! (out)
609  teml(ks:ke_pbl,:,:), pres(ks:ke_pbl,:,:), & ! (in)
610  ke_pbl-ks+1 ) ! (in)
611 
612 !OCL LOOP_NOFUSION,PREFETCH_SEQUENTIAL(SOFT),SWP
613  do j = js, je
614  do i = is, ie
615  do k = ks+1, ke_pbl
616 
617  dqsl = qsl(k,i,j) * lh(k,i,j) / ( rvap * potl(k,i,j)**2 )
618  aa = 1.0_rp / ( 1.0_rp + lh(k,i,j)/cp * dqsl )
619  bb = temp(k,i,j) / pott(k,i,j) * dqsl
620  ac = min( q(k,i,j)/sqrt(q2_2(k,i,j)), 1.0_rp )
621  sigma_s = max( sqrt( 0.25_rp * aa**2 * l(k,i,j)**2 * ac * b2 * sh(k,i,j) ) &
622  * abs( qw(k+1,i,j) - qw(k-1,i,j) - bb * ( potl(k+1,i,j)-potl(k-1,i,j) ) ) &
623  / ( cz(k+1,i,j) - cz(k-1,i,j) ), &
624  1e-10_rp )
625  q1 = aa * ( qw(k,i,j) - qsl(k,i,j) ) * 0.5_rp / sigma_s
626  rr = min( max( 0.5_rp * ( 1.0_rp + erf(q1*rsqrt_2) ), 0.0_rp ), 1.0_rp )
627  ql = min( max( 2.0_rp * sigma_s * ( rr * q1 + rsqrt_2pi * exp(-0.5_rp*q1**2) ), &
628  0.0_rp ), &
629  qw(k,i,j) * 0.5_rp )
630  cc = ( 1.0_rp + epstvap * qw(k,i,j) - (1.0_rp+epstvap) * ql ) * pott(k,i,j)/temp(k,i,j) * lh(k,i,j) / cp &
631  - (1.0_rp+epstvap) * pott(k,i,j)
632  rt = min( max( rr - ql / (2.0_rp*sigma_s*sqrt_2pi) * exp(-q1**2 * 0.5_rp), 0.0_rp ), 1.0_rp )
633  betat = 1.0_rp + epstvap * qw(k,i,j) - (1.0_rp+epstvap) * ql - rt * aa * bb * cc
634  betaq = epstvap * pott(k,i,j) + rt * aa * cc
635  n2(k,i,j) = min(atmos_phy_tb_mynn_n2_max, &
636  grav * ( ( potl(k+1,i,j) - potl(k-1,i,j) ) * betat &
637  + ( qw(k+1,i,j) - qw(k-1,i,j) ) * betaq ) &
638  / ( cz(k+1,i,j) - cz(k-1,i,j) ) / potv(k,i,j) )
639  end do
640  n2(ks,i,j) = n2(ks+1,i,j)
641 
642  do k = ks, ke_pbl
643  ri(k,i,j) = n2(k,i,j) / max(dudz2(k,i,j), 1e-10_rp)
644  end do
645  do k = ke_pbl+1, ke
646  ri(k,i,j) = 0.0_rp
647  n2(k,i,j) = 0.0_rp
648  end do
649 
650  end do
651  end do
652 
653  ! length
654  call get_length( &
655  l, & ! (out)
656  dens, & ! (in)
657  q, n2, & ! (in)
658  sflx_mu, sflx_mv, sflx_pt, & ! (in)
659  pott ) ! (in)
660 
661  call get_q2_level2( &
662  q2_2, & ! (out)
663  dudz2, ri, l ) ! (in)
664 
665  call get_smsh( &
666  sm, sh, & ! (out)
667  q, q2_2, & ! (in)
668  l, n2, dudz2 ) ! (in)
669 
670 
671  do j = js, je
672  do i = is, ie
673  do k = 1, ks-1
674  nu(k,i,j) = 0.0_rp
675  end do
676  do k = ks, ke_pbl
677  nu(k,i,j) = max( min( l(k,i,j) * q(k,i,j) * sm(k,i,j), &
678  atmos_phy_tb_mynn_nu_max ), &
679  atmos_phy_tb_mynn_nu_min )
680  end do
681  do k = ke_pbl+1, ka
682  nu(k,i,j) = 0.0_rp
683  end do
684 
685  do k = ks, ke_pbl
686  kh = max( min( l(k,i,j) * q(k,i,j) * sh(k,i,j), &
687  atmos_phy_tb_mynn_kh_max ), &
688  atmos_phy_tb_mynn_kh_min )
689  rhokh(k,i,j) = dens(k,i,j) * kh
690  pr(k,i,j) = nu(k,i,j) / kh
691  end do
692  do k = ke_pbl+1, ke
693  pr(k,i,j) = 1.0_rp
694  end do
695  end do
696  end do
697 
698  call comm_vars( nu, 1 )
699 
700  ! time integration
701 
702  ! for velocities
703 !OCL INDEPENDENT
704  do j = js, je
705  do i = is, ie
706 
707  ap = - dt * 0.5_rp * ( dens(ks ,i,j)*nu(ks ,i,j) &
708  + dens(ks+1,i,j)*nu(ks+1,i,j) ) &
709  * rfdz(ks) / gsqrt(ks,i,j,i_xyw)
710  a(ks,i,j) = ap * rcdz(ks) / ( dens(ks,i,j) * gsqrt(ks,i,j,i_xyz) )
711  c(ks,i,j) = 0.0_rp
712  b(ks,i,j) = - a(ks,i,j) + 1.0_rp
713  do k = ks+1, ke_pbl-1
714  c(k,i,j) = ap * rcdz(k) / ( dens(k,i,j) * gsqrt(k,i,j,i_xyz) )
715  ap = - dt * 0.5_rp * ( dens(k ,i,j)*nu(k ,i,j) &
716  + dens(k+1,i,j)*nu(k+1,i,j) ) &
717  * rfdz(k) / gsqrt(k,i,j,i_xyw)
718  a(k,i,j) = ap * rcdz(k) / ( dens(k,i,j) * gsqrt(k,i,j,i_xyz) )
719  b(k,i,j) = - a(k,i,j) - c(k,i,j) + 1.0_rp
720  end do
721  a(ke_pbl,i,j) = 0.0_rp
722  c(ke_pbl,i,j) = ap * rcdz(ke_pbl) / ( dens(ke_pbl,i,j) * gsqrt(ke_pbl,i,j,i_xyz) )
723  b(ke_pbl,i,j) = - c(ke_pbl,i,j) + 1.0_rp
724 
725  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) )
726  do k = ks+1, ke_pbl
727  d(k,i,j) = u(k,i,j)
728  end do
729  call diffusion_solver( &
730  phin(:,i,j), & ! (out)
731  a(:,i,j), b(:,i,j), c(:,i,j), d(:,i,j), & ! (in)
732  ke_pbl ) ! (in)
733  end do
734  end do
735 
736  call comm_vars( phin, 2 )
737  call comm_wait( nu, 1 )
738  call comm_wait( phin, 2 )
739 
740 !OCL INDEPENDENT
741  do jjs = js, je, jblock
742  jje = jjs+jblock-1
743  do iis = is, ie, iblock
744  iie = iis+iblock-1
745 
746  do j = jjs, jje
747  do i = iis, iie
748  do k = ks, ke_pbl-1
749  qflx_sgs_momx(k,i,j,zdir) = - 0.03125_rp & ! 1/4/4/2
750  * ( dens(k,i,j) + dens(k+1,i,j) + dens(k,i+1,j) + dens(k+1,i+1,j) ) &
751  * ( nu(k,i,j) + nu(k+1,i,j) + nu(k,i+1,j) + nu(k+1,i+1,j) ) &
752  * ( (phin(k+1,i,j)+phin(k+1,i+1,j)) - (phin(k,i,j)+phin(k,i+1,j)) ) &
753  * j33g * rfdz(k) / gsqrt(k,i,j,i_uyw)
754  end do
755  end do
756  end do
757  do j = jjs, jje
758  do i = iis, iie
759  qflx_sgs_momx(ks-1,i,j,zdir) = 0.0_rp
760  end do
761  end do
762  do j = jjs, jje
763  do i = iis, iie
764  do k = ke_pbl, ke
765  qflx_sgs_momx(k,i,j,zdir) = 0.0_rp
766  end do
767  end do
768  end do
769 
770 
771  ! integration V
772 !OCL INDEPENDENT
773  do j = jjs, jje
774  do i = iis, iie
775  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) )
776  do k = ks+1, ke_pbl
777  d(k,i,j) = v(k,i,j)
778  end do
779  call diffusion_solver( &
780  phin(:,i,j), & ! (out)
781  a(:,i,j), b(:,i,j), c(:,i,j), d(:,i,j), & ! (in)
782  ke_pbl ) ! (in)
783  end do
784  end do
785 
786  end do
787  end do
788 
789  call comm_vars( phin, 1 )
790  call comm_wait( phin, 1 )
791 
792  do jjs = js, je, jblock
793  jje = jjs+jblock-1
794  do iis = is, ie, iblock
795  iie = iis+iblock-1
796 
797  do j = jjs, jje
798  do i = iis, iie
799  qflx_sgs_momy(ks-1,i,j,zdir) = 0.0_rp
800  do k = ks, ke_pbl-1
801  qflx_sgs_momy(k,i,j,zdir) = - 0.03125_rp & ! 1/4/4/2
802  * ( dens(k,i,j) + dens(k+1,i,j) + dens(k,i,j+1) + dens(k+1,i,j+1) ) &
803  * ( nu(k,i,j) + nu(k+1,i,j) + nu(k,i,j+1) + nu(k+1,i,j+1) ) &
804  * ( (phin(k+1,i,j)+phin(k+1,i,j+1)) - (phin(k,i,j)+phin(k,i,j+1)) ) &
805  * j33g * rfdz(k) / gsqrt(k,i,j,i_xvw)
806  end do
807  do k = ke_pbl, ke
808  qflx_sgs_momy(k,i,j,zdir) = 0.0_rp
809  end do
810  end do
811  end do
812 
813 
814  ! for scalars
815  ! integration POTT
816 !OCL INDEPENDENT
817  do j = jjs, jje
818  do i = iis, iie
819  ap = - dt * 0.5_rp * ( rhokh(ks,i,j) + rhokh(ks+1,i,j) ) &
820  * rfdz(ks) / gsqrt(ks,i,j,i_xyw)
821  a(ks,i,j) = ap * rcdz(ks) / (dens(ks,i,j) * gsqrt(ks,i,j,i_xyz) )
822  c(ks,i,j) = 0.0_rp
823  b(ks,i,j) = - a(ks,i,j) + 1.0_rp
824  do k = ks+1, ke_pbl-1
825  c(k,i,j) = ap * rcdz(k) / (dens(k,i,j) * gsqrt(k,i,j,i_xyz) )
826  ap = - dt * 0.5_rp * ( rhokh(k,i,j) + rhokh(k+1,i,j) ) &
827  * rfdz(k) / gsqrt(k,i,j,i_xyw)
828  a(k,i,j) = ap * rcdz(k) / ( dens(k,i,j) * gsqrt(k,i,j,i_xyz) )
829  b(k,i,j) = - a(k,i,j) - c(k,i,j) + 1.0_rp
830  end do
831  a(ke_pbl,i,j) = 0.0_rp
832  c(ke_pbl,i,j) = ap * rcdz(ke_pbl) / (dens(ke_pbl,i,j) * gsqrt(ke_pbl,i,j,i_xyz) )
833  b(ke_pbl,i,j) = - c(ke_pbl,i,j) + 1.0_rp
834  d(ks,i,j) = potl(ks,i,j) + dt * sflx_pt(i,j) * rcdz(ks) / gsqrt(ks,i,j,i_xyz)
835  do k = ks+1, ke_pbl
836  d(k,i,j) = potl(k,i,j)
837  end do
838  call diffusion_solver( &
839  phin(:,i,j), & ! (out)
840  a(:,i,j), b(:,i,j), c(:,i,j), d(:,i,j), & ! (in)
841  ke_pbl ) ! (in)
842  end do
843  end do
844 
845  do j = jjs, jje
846  do i = iis, iie
847  qflx_sgs_rhot(ks-1,i,j,zdir) = 0.0_rp
848  do k = ks, ke_pbl-1
849  qflx_sgs_rhot(k,i,j,zdir) = - 0.5_rp & ! 1/2
850  * ( rhokh(k,i,j) + rhokh(k+1,i,j) ) &
851  * j33g * ( phin(k+1,i,j) - phin(k,i,j) ) * rfdz(k) / gsqrt(k,i,j,i_xyw)
852  end do
853  do k = ke_pbl, ke
854  qflx_sgs_rhot(k,i,j,zdir) = 0.0_rp
855  end do
856  end do
857  end do
858 
859 
860  ! integration QV
861  iq = i_qv
862 !OCL INDEPENDENT
863  do j = jjs, jje
864  do i = iis, iie
865  d(ks,i,j) = qw(ks,i,j) + dt * sflx_qv(i,j) * rcdz(ks) / ( dens(ks,i,j) * gsqrt(ks,i,j,i_xyz) )
866  do k = ks+1, ke_pbl
867  d(k,i,j) = qw(k,i,j)
868  end do
869  call diffusion_solver( &
870  phin(:,i,j), & ! (out)
871  a(:,i,j), b(:,i,j), c(:,i,j), d(:,i,j), & ! (in)
872  ke_pbl ) ! (in)
873  end do
874  end do
875 
876  do j = jjs, jje
877  do i = iis, iie
878  qflx_sgs_rhoq(ks-1,i,j,zdir,iq) = 0.0_rp
879  do k = ks, ke_pbl-1
880  qflx_sgs_rhoq(k,i,j,zdir,iq) = - 0.5_rp & ! 1/2
881  * ( rhokh(k,i,j) + rhokh(k+1,i,j) ) &
882  * j33g * ( phin(k+1,i,j) - phin(k,i,j) ) * rfdz(k) / gsqrt(k,i,j,i_xyw)
883  end do
884  do k = ke_pbl, ke
885  qflx_sgs_rhoq(k,i,j,zdir,iq) = 0.0_rp
886  end do
887  end do
888  end do
889 
890 
891  ! time integration tke
892 
893  do j = jjs, jje
894  do i = iis, iie
895  do k = ks , ke_pbl-1
896  mflx = j33g * momz(k,i,j) / ( mapf(i,j,1,i_xy)*mapf(i,j,2,i_xy) ) &
897  + j13g(k,i,j,i_xyw) * 0.25_rp * ( momx(k,i-1,j)+momx(k,i,j)+momx(k+1,i-1,j)+momx(k+1,i,j) ) / mapf(i,j,2,i_xy) &
898  + j23g(k,i,j,i_xyw) * 0.25_rp * ( momy(k,i,j-1)+momy(k,i,j)+momy(k+1,i,j-1)+momy(k+1,i,j) ) / mapf(i,j,1,i_xy)
899  flux_z(k,i,j) = 0.5_rp * ( mflx * ( tke(k+1,i,j)+tke(k,i,j) ) &
900  -abs(mflx) * ( tke(k+1,i,j)-tke(k,i,j) ) )
901  end do
902  end do
903  end do
904  do j = jjs , jje
905  do i = iis-1, iie
906  do k = ks , ke_pbl
907  mflx = momx(k,i,j) * gsqrt(k,i,j,i_uyz) / mapf(i,j,2,i_uy)
908  flux_x(k,i,j) = 0.5_rp * ( mflx * ( tke(k,i+1,j)+tke(k,i,j) ) &
909  -abs(mflx) * ( tke(k,i+1,j)-tke(k,i,j) ) )
910  end do
911  end do
912  end do
913  do j = jjs-1, jje
914  do i = iis , iie
915  do k = ks , ke_pbl
916  mflx = momy(k,i,j) * gsqrt(k,i,j,i_xvz) / mapf(i,j,1,i_xv)
917  flux_y(k,i,j) = 0.5_rp * ( mflx * ( tke(k,i,j+1)+tke(k,i,j) ) &
918  -abs(mflx) * ( tke(k,i,j+1)-tke(k,i,j) ) )
919  end do
920  end do
921  end do
922 
923  do j = jjs, jje
924  do i = iis, iie
925  advc = ( ( flux_z(ks,i,j) ) * rcdz(ks) &
926  + ( flux_x(ks,i,j) - flux_x(ks,i-1,j) ) * rcdx(i) &
927  + ( flux_y(ks,i,j) - flux_y(ks,i,j-1) ) * rcdy(j) ) &
928  * mapf(i,j,1,i_xy) * mapf(i,j,2,i_xy) / ( gsqrt(ks,i,j,i_xyz) * dens(ks,i,j) )
929  d(ks,i,j) = tke(ks,i,j) + dt * ( nu(ks,i,j) * (dudz2(ks,i,j) - n2(ks,i,j)/pr(ks,i,j)) &
930  - advc )
931 #ifdef MORE_HIST
932  adv(ks,i,j) = advc
933  gen(ks,i,j) = nu(ks,i,j) * (dudz2(ks,i,j) - n2(ks,i,j)/pr(ks,i,j))
934 #endif
935  do k = ks+1, ke_pbl-1
936  advc = ( ( flux_z(k,i,j) - flux_z(k-1,i,j) ) * rcdz(k) &
937  + ( flux_x(k,i,j) - flux_x(k,i-1,j) ) * rcdx(i) &
938  + ( flux_y(k,i,j) - flux_y(k,i,j-1) ) * rcdy(j) ) &
939  * mapf(i,j,1,i_xy) * mapf(i,j,2,i_xy) / ( gsqrt(k,i,j,i_xyz) * dens(k,i,j) )
940  d(k,i,j) = tke(k,i,j) &
941  + dt * ( nu(k,i,j) * (dudz2(k,i,j) - n2(k,i,j)/pr(k,i,j)) &
942  - advc )
943 #ifdef MORE_HIST
944  adv(k,i,j) = advc
945  gen(k,i,j) = nu(k,i,j) * (dudz2(k,i,j) - n2(k,i,j)/pr(k,i,j))
946 #endif
947  end do
948  advc = ( ( - flux_z(ke_pbl-1,i,j) ) * rcdz(ke_pbl) &
949  + ( flux_x(ke_pbl,i,j) - flux_x(ke_pbl,i-1,j) ) * rcdx(i) &
950  + ( flux_y(ke_pbl,i,j) - flux_y(ke_pbl,i,j-1) ) * rcdy(j) ) &
951  * mapf(i,j,1,i_xy) * mapf(i,j,2,i_xy) / ( gsqrt(ke_pbl,i,j,i_xyz) * dens(ke_pbl,i,j) )
952  d(ke_pbl,i,j) = tke(ke_pbl,i,j) + dt * ( nu(ke_pbl,i,j) * (dudz2(ke_pbl,i,j) - n2(ke_pbl,i,j)/pr(ke_pbl,i,j)) &
953  - advc )
954 #ifdef MORE_HIST
955  adv(ke_pbl,i,j) = advc
956  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))
957 #endif
958  end do
959  end do
960 
961 !OCL INDEPENDENT
962  do j = jjs, jje
963  do i = iis, iie
964  ap = - dt * 1.5_rp * ( dens(ks ,i,j)*nu(ks ,i,j) &
965  + dens(ks+1,i,j)*nu(ks+1,i,j) ) &
966  * rfdz(ks) / gsqrt(ks,i,j,i_xyw)
967  a(ks,i,j) = ap * rcdz(ks) / ( dens(ks,i,j) * gsqrt(ks,i,j,i_xyz) )
968  c(ks,i,j) = 0.0_rp
969  b(ks,i,j) = - a(ks,i,j) + 1.0_rp + 2.0_rp * dt * q(ks,i,j) / ( b1 * l(ks,i,j) )
970  do k = ks+1, ke_pbl-1
971  c(k,i,j) = ap * rcdz(k) / ( dens(k,i,j) * gsqrt(k,i,j,i_xyz) )
972  ap = - dt * 1.5_rp * ( dens(k ,i,j)*nu(k ,i,j) &
973  + dens(k+1,i,j)*nu(k+1,i,j)) &
974  * rfdz(k) / gsqrt(k,i,j,i_xyw)
975  a(k,i,j) = ap * rcdz(k) / ( dens(k,i,j) * gsqrt(k,i,j,i_xyz) )
976  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) )
977  end do
978 
979  a(ke_pbl,i,j) = 0.0_rp
980  c(ke_pbl,i,j) = ap * rcdz(ke_pbl) / ( dens(ke_pbl,i,j) * gsqrt(ke_pbl,i,j,i_xyz) )
981  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) )
982  call diffusion_solver( &
983  tke_n(:,i,j), & ! (out)
984  a(:,i,j), b(:,i,j), c(:,i,j), d(:,i,j), & ! (in)
985  ke_pbl ) ! (in)
986 #ifdef DEBUG
987  do k = ks+1, ke_pbl-1
988  tke(1,i,j) = d(k,i,j) &
989  + ( ( (tke_n(k+1,i,j)-tke_n(k,i,j)) * rfdz(k) / gsqrt(k,i,j,i_xyw) &
990  * (nu(k+1,i,j)*dens(k+1,i,j)+nu(k,i,j)*dens(k,i,j))*1.5_rp &
991  - (tke_n(k,i,j)-tke_n(k-1,i,j)) * rfdz(k-1) / gsqrt(k-1,i,j,i_xyw) &
992  * (nu(k,i,j)*dens(k,i,j)+nu(k-1,i,j)*dens(k-1,i,j))*1.5_rp ) &
993  * rcdz(k) / gsqrt(k,i,j,i_xyz) / dens(k,i,j) &
994  - 2.0_rp*tke_n(k,i,j) * q(k,i,j) / (b1 * l(k,i,j)) ) &
995  * dt
996  if ( tke_n(k,i,j) > 0.1_rp .AND. abs(tke(1,i,j) - tke_n(k,i,j))/tke_n(k,i,j) > 1e-10_rp ) then
997  advc = ( tke(k,i,j) - d(k,i,j) ) / dt + nu(k,i,j) * (dudz2(k,i,j) - n2(k,i,j)/pr(k,i,j))
998  write(*,*)k,i,j,tke(1,i,j),tke_n(k,i,j), tke(k,i,j),nu(k,i,j),dudz2(k,i,j),n2(k,i,j),pr(k,i,j),advc
999  open(90, file="mynn.dat")
1000  write(90,*)ke_pbl-ks+1
1001  write(90,*)dt, b1
1002  write(90,*)tke(ks:ke_pbl,i,j)
1003  write(90,*)q(ks:ke_pbl,i,j)
1004  write(90,*)l(ks:ke_pbl,i,j)
1005  write(90,*)rcdz(ks:ke_pbl)
1006  write(90,*)rfdz(ks:ke_pbl)
1007  write(90,*)nu(ks:ke_pbl,i,j)
1008  write(90,*)dens(ks:ke_pbl,i,j)
1009  write(90,*)dudz2(ks:ke_pbl,i,j)
1010  write(90,*)n2(ks:ke_pbl,i,j)
1011  write(90,*)pr(ks:ke_pbl,i,j)
1012  close(90)
1013  call abort
1014  end if
1015  end do
1016 #endif
1017  do k = ks, ke_pbl
1018  tke_t(k,i,j) = ( max(tke_n(k,i,j), tke_min) - tke(k,i,j) ) / dt
1019  end do
1020  do k = ke_pbl+1, ke
1021  tke_t(k,i,j) = 0.0_rp
1022  end do
1023 
1024  end do
1025  end do
1026  end do
1027  end do
1028 
1029 #ifdef MORE_HIST
1030  adv(ke_pbl+1:ke,:,:) = 0.0_rp
1031  gen(ke_pbl+1:ke,:,:) = 0.0_rp
1032  dudz2(ke_pbl+1:ke,:,:) = 0.0_rp
1033  l(ke_pbl+1:ke,:,:) = 0.0_rp
1034  potv(ke_pbl+1:ke,:,:) = 0.0_rp
1035  potl(ke_pbl+1:ke,:,:) = 0.0_rp
1036  call hist_in(adv, 'TKE_advc', 'advection of TKE', 'm2/s3', nohalo=.true.)
1037  call hist_in(gen, 'TKE_gen', 'generation of TKE', 'm2/s3', nohalo=.true.)
1038  call hist_in(dudz2, 'dUdZ2', 'dudz2', 'm2/s2', nohalo=.true.)
1039  call hist_in(l, 'L_mix', 'minxing length', 'm', nohalo=.true.)
1040  call hist_in(potv, 'POTV', 'virtual potential temperature', 'K', nohalo=.true.)
1041  call hist_in(potl, 'POTL', 'liquid potential temperature', 'K', nohalo=.true.)
1042 #endif
1043 
1044  return
scale_grid_index::is
integer, public is
start point of inner domain: x, local
Definition: scale_grid_index.F90:60
scale_atmos_dyn_fvm_flux::flux_z
Definition: scale_atmos_dyn_fvm_flux.F90:94
scale_gridtrans::i_xvz
integer, public i_xvz
Definition: scale_gridtrans.F90:54
scale_grid_index::je
integer, public je
end point of inner domain: y, local
Definition: scale_grid_index.F90:63
scale_grid::grid_rcdy
real(rp), dimension(:), allocatable, public grid_rcdy
reciprocal of center-dy
Definition: scale_grid_cartesian.F90:64
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_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:1453
scale_grid_index::iblock
integer, public iblock
block size for cache blocking: x
Definition: scale_grid_index.F90:47
scale_tracer::qqe
integer, public qqe
Definition: scale_tracer.F90:61
scale_grid::grid_rcdx
real(rp), dimension(:), allocatable, public grid_rcdx
reciprocal of center-dx
Definition: scale_grid_cartesian.F90:63
scale_grid_index::zdir
integer, parameter, public zdir
Definition: scale_grid_index.F90:35
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:44
scale_grid_real::real_cz
real(rp), dimension(:,:,:), allocatable, public real_cz
geopotential height [m] (cell center)
Definition: scale_grid_real.F90:40
scale_gridtrans::i_xy
integer, public i_xy
Definition: scale_gridtrans.F90:57
scale_grid::grid_rcdz
real(rp), dimension(:), allocatable, public grid_rcdz
reciprocal of center-dz
Definition: scale_grid_cartesian.F90:62
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_grid_index
module grid index
Definition: scale_grid_index.F90:14
scale_tracer
module TRACER
Definition: scale_tracer.F90:18
scale_grid_index::ia
integer, public ia
of x whole cells (local, with HALO)
Definition: scale_grid_index.F90:55
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 z whole cells (local, with HALO)
Definition: scale_grid_index.F90:54
scale_gridtrans::i_uy
integer, public i_uy
Definition: scale_gridtrans.F90:58
scale_tracer::i_qv
integer, public i_qv
Definition: scale_tracer.F90:47
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_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_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:14
scale_grid
module GRID (cartesian)
Definition: scale_grid_cartesian.F90:17
scale_gridtrans::i_xv
integer, public i_xv
Definition: scale_gridtrans.F90:59
scale_tracer::i_qi
integer, public i_qi
Definition: scale_tracer.F90:50
scale_grid_index::ie
integer, public ie
end point of inner domain: x, local
Definition: scale_grid_index.F90:61
scale_atmos_thermodyn
module ATMOSPHERE / Thermodynamics
Definition: scale_atmos_sub_thermodyn.F90:19
scale_tracer::qqs
integer, public qqs
Definition: scale_tracer.F90:60
scale_gridtrans::i_uyz
integer, public i_uyz
Definition: scale_gridtrans.F90:53
scale_precision
module PRECISION
Definition: scale_precision.F90:13
scale_history
module HISTORY
Definition: scale_history.F90:16
scale_grid::grid_rfdz
real(rp), dimension(:), allocatable, public grid_rfdz
reciprocal of face-dz
Definition: scale_grid_cartesian.F90:72
scale_gridtrans::i_xyz
integer, public i_xyz
Definition: scale_gridtrans.F90:49
scale_tracer::i_qc
integer, public i_qc
Definition: scale_tracer.F90:48
scale_grid_index::ja
integer, public ja
of y whole cells (local, with HALO)
Definition: scale_grid_index.F90:56
Here is the call graph for this function:
Here is the caller graph for this function:

◆ get_length()

subroutine scale_atmos_phy_tb_mynn::get_length ( real(rp), dimension(ka,ia,ja), intent(out)  l,
real(rp), dimension(ka,ia,ja), intent(in)  DENS,
real(rp), dimension(ka,ia,ja), intent(in)  q,
real(rp), dimension(ka,ia,ja), intent(in)  n2,
real(rp), dimension(ia,ja), intent(in)  SFLX_MU,
real(rp), dimension(ia,ja), intent(in)  SFLX_MV,
real(rp), dimension(ia,ja), intent(in)  SFLX_PT,
real(rp), dimension(ka,ia,ja), intent(in)  PT0 
)

Definition at line 1053 of file scale_atmos_phy_tb_mynn.F90.

References scale_const::const_cpdry, scale_const::const_eps, scale_const::const_grav, scale_const::const_karman, scale_grid_index::ie, scale_grid_index::is, scale_grid_index::je, scale_grid_index::js, scale_grid_index::ks, and scale_grid_real::real_fz.

Referenced by atmos_phy_tb_mynn().

1053  use scale_grid_real, only: &
1054  fz => real_fz
1055  use scale_const, only: &
1056  grav => const_grav, &
1057  karman => const_karman, &
1058  cp => const_cpdry, &
1059  eps => const_eps
1060  implicit none
1061 
1062  real(RP), intent(out) :: l(ka,ia,ja)
1063  real(RP), intent(in) :: dens(ka,ia,ja)
1064  real(RP), intent(in) :: q(ka,ia,ja)
1065  real(RP), intent(in) :: n2(ka,ia,ja)
1066  real(RP), intent(in) :: sflx_mu(ia,ja)
1067  real(RP), intent(in) :: sflx_mv(ia,ja)
1068  real(RP), intent(in) :: sflx_pt(ia,ja)
1069  real(RP), intent(in) :: pt0(ka,ia,ja)
1070 
1071  real(RP) :: ls
1072  real(RP) :: lt
1073  real(RP) :: lb
1074  real(RP) :: rlm
1075  real(RP) :: rlt
1076 
1077  real(RP) :: qc
1078  real(RP) :: int_q
1079  real(RP) :: int_qz
1080  real(RP) :: rn2sr
1081  real(RP) :: us
1082  real(RP) :: zeta
1083 
1084  real(RP) :: z
1085  real(RP) :: qdz
1086 
1087  real(RP) :: sw
1088  integer :: k, i, j
1089 
1090 !OCL LOOP_NOFUSION,PREFETCH_SEQUENTIAL(SOFT),SWP
1091  do j = js, je
1092  do i = is, ie
1093  int_qz = 0.0_rp
1094  int_q = 0.0_rp
1095  do k = ks, ke_pbl
1096  qdz = q(k,i,j) * ( fz(k,i,j) - fz(k-1,i,j) )
1097  int_qz = int_qz + ((fz(k,i,j)+fz(k-1,i,j))*0.5_rp-fz(ks-1,i,j)) * qdz
1098  int_q = int_q + qdz
1099  end do
1100  ! LT
1101  lt = min( max(0.23_rp * int_qz / (int_q + eps), &
1102  lt_min), &
1103  atmos_phy_tb_mynn_lt_max )
1104  rlt = 1.0_rp / lt
1105 
1106  us = ( sflx_mu(i,j)**2 + sflx_mv(i,j)**2 )**0.25_rp / dens(ks,i,j) ! friction velocity
1107  us = max(us, us_min)
1108  rlm = - karman * grav * sflx_pt(i,j) / (pt0(ks,i,j) * us**3 )
1109 
1110  qc = (grav/pt0(ks,i,j)*max(sflx_pt(i,j),0.0_rp)*lt)**oneoverthree
1111 
1112  do k = ks, ke_pbl
1113  z = ( fz(k,i,j)+fz(k-1,i,j) )*0.5_rp - fz(ks-1,i,j)
1114  zeta = z * rlm
1115 
1116  ! LS
1117  sw = sign(0.5_rp, zeta) + 0.5_rp ! 1 for zeta >= 0, 0 for zeta < 0
1118  ls = karman * z &
1119  * ( sw / (1.0_rp + 2.7_rp*min(zeta,1.0_rp)*sw ) &
1120  + ( (1.0_rp - 100.0_rp*zeta)*(1.0_rp-sw) )**0.2_rp )
1121 
1122  ! LB
1123  sw = sign(0.5_rp, n2(k,i,j)) + 0.5_rp ! 1 for dptdz >0, 0 for dptdz < 0
1124  rn2sr = 1.0_rp / ( sqrt(n2(k,i,j)*sw) + 1.0_rp-sw)
1125  lb = (1.0_rp + 5.0_rp * sqrt(qc*rn2sr/lt)) * q(k,i,j) * rn2sr * sw & ! qc=0 when SFLX_PT < 0
1126  + 999.e10_rp * (1.0_rp-sw)
1127 
1128  ! L
1129  l(k,i,j) = 1.0_rp / ( 1.0_rp/ls + rlt + 1.0_rp/lb )
1130  end do
1131  end do
1132  end do
1133 
1134  return
scale_grid_index::is
integer, public is
start point of inner domain: x, local
Definition: scale_grid_index.F90:60
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_grid_real::real_fz
real(rp), dimension(:,:,:), allocatable, public real_fz
geopotential height [m] (cell face )
Definition: scale_grid_real.F90:41
scale_const::const_karman
real(rp), parameter, public const_karman
von Karman constant
Definition: scale_const.F90:52
scale_grid_index::ia
integer, public ia
of x whole cells (local, with HALO)
Definition: scale_grid_index.F90:55
scale_grid_real
module GRID (real space)
Definition: scale_grid_real.F90:14
scale_grid_index::ka
integer, public ka
of z whole cells (local, with HALO)
Definition: scale_grid_index.F90:54
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
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_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_grid_index::ja
integer, public ja
of y whole cells (local, with HALO)
Definition: scale_grid_index.F90:56
Here is the caller graph for this function:

◆ get_q2_level2()

subroutine scale_atmos_phy_tb_mynn::get_q2_level2 ( real(rp), dimension(ka,ia,ja), intent(out)  q2_2,
real(rp), dimension(ka,ia,ja), intent(in)  dudz2,
real(rp), dimension(ka,ia,ja), intent(in)  Ri,
real(rp), dimension(ka,ia,ja), intent(in)  l 
)

Definition at line 1140 of file scale_atmos_phy_tb_mynn.F90.

References scale_grid_index::ie, scale_grid_index::is, scale_grid_index::je, scale_grid_index::js, and scale_grid_index::ks.

Referenced by atmos_phy_tb_mynn().

1140  implicit none
1141 
1142  real(RP), intent(out) :: q2_2(ka,ia,ja)
1143  real(RP), intent(in) :: dudz2(ka,ia,ja)
1144  real(RP), intent(in) :: ri(ka,ia,ja)
1145  real(RP), intent(in) :: l(ka,ia,ja)
1146 
1147  real(RP) :: rf
1148  real(RP) :: sm_2
1149  real(RP) :: sh_2
1150 
1151  real(RP) :: q2
1152  integer :: k, i, j
1153 
1154 !OCL LOOP_NOFUSION,PREFETCH_SEQUENTIAL(SOFT),SWP
1155  do j = js, je
1156  do i = is, ie
1157  do k = ks, ke_pbl
1158  rf = min(0.5_rp / af12 * ( ri(k,i,j) &
1159  + af12*rf1 &
1160  - sqrt(ri(k,i,j)**2 + 2.0_rp*af12*(rf1-2.0_rp*rf2)*ri(k,i,j) + (af12*rf1)**2) ), &
1161  rfc)
1162  sh_2 = 3.0_rp * a2 * (g1+g2) * (rfc-rf) / (1.0_rp-rf)
1163  sm_2 = sh_2 * af12 * (rf1-rf) / (rf2-rf)
1164  q2 = b1 * l(k,i,j)**2 * sm_2 * (1.0_rp-rf) * dudz2(k,i,j)
1165  q2_2(k,i,j) = max( q2, 1.e-10_rp )
1166  end do
1167  end do
1168  end do
1169 
1170  return
scale_grid_index::is
integer, public is
start point of inner domain: x, local
Definition: scale_grid_index.F90:60
scale_grid_index::je
integer, public je
end point of inner domain: y, local
Definition: scale_grid_index.F90:63
scale_grid_index::ia
integer, public ia
of x whole cells (local, with HALO)
Definition: scale_grid_index.F90:55
scale_grid_index::ka
integer, public ka
of z whole cells (local, with HALO)
Definition: scale_grid_index.F90:54
scale_grid_index::js
integer, public js
start point of inner domain: y, local
Definition: scale_grid_index.F90:62
scale_grid_index::ks
integer, public ks
start point of inner domain: z, local
Definition: scale_grid_index.F90:58
scale_grid_index::ie
integer, public ie
end point of inner domain: x, local
Definition: scale_grid_index.F90:61
scale_grid_index::ja
integer, public ja
of y whole cells (local, with HALO)
Definition: scale_grid_index.F90:56
Here is the caller graph for this function:
Generated by   doxygen 1.8.13