SCALE-RM
scale_atmos_phy_tb_smg.F90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
19 !-------------------------------------------------------------------------------
20 #include "scalelib.h"
21 module scale_atmos_phy_tb_smg
22  !-----------------------------------------------------------------------------
23  !
24  !++ used modules
25  !
26  use scale_precision
27  use scale_io
28  use scale_prof
29  use scale_atmos_grid_cartesc_index
30  use scale_tracer
31 
32 #if defined DEBUG || defined QUICKDEBUG
33  use scale_debug, only: &
34  check
35  use scale_const, only: &
36  undef => const_undef, &
37  iundef => const_undef2
38 #endif
39  !-----------------------------------------------------------------------------
40  implicit none
41  private
42  !-----------------------------------------------------------------------------
43  !
44  !++ Public procedure
45  !
46  public :: atmos_phy_tb_smg_setup
47  public :: atmos_phy_tb_smg
48 
49  !-----------------------------------------------------------------------------
50  !
51  !++ Public parameters & variables
52  !
53  !-----------------------------------------------------------------------------
54  !
55  !++ Private procedure
56  !
57  !-----------------------------------------------------------------------------
58  !
59  !++ Private parameters & variables
60  !
61  real(RP), private, parameter :: oneoverthree = 1.0_rp / 3.0_rp
62  real(RP), private, parameter :: twooverthree = 2.0_rp / 3.0_rp
63  real(RP), private, parameter :: fouroverthree = 4.0_rp / 3.0_rp
64 
65  real(RP), private :: cs = 0.13_rp ! Smagorinsky constant (Scotti et al. 1993)
66  real(RP), private, parameter :: prn = 0.7_rp ! Prandtl number in neutral conditions
67  real(RP), private, parameter :: ric = 0.25_rp ! critical Richardson number
68  real(RP), private, parameter :: fmc = 16.0_rp ! fum = sqrt(1 - c*Ri)
69  real(RP), private, parameter :: fhb = 40.0_rp ! fuh = sqrt(1 - b*Ri)/PrN
70  real(RP), private :: rprn ! 1 / PrN
71  real(RP), private :: rric ! 1 / RiC
72  real(RP), private :: prnovric ! PrN / RiC
73 
74  ! for backscatter
75  real(RP), private, parameter :: cb = 1.4_rp
76  real(RP), private, parameter :: cbt = 0.45_rp
77  real(RP), private, parameter :: an = 0.47958315233127197_rp ! a_N = sqrt(0.23)
78  real(RP), private, parameter :: atn4 = 0.09_rp ! a_{\theta N}^4 = 0.3**2
79  real(RP), private, parameter :: c1o = an**3
80  real(RP), private, parameter :: d1o = prn * atn4 / an
81 
82  real(RP), private, allocatable :: lambda0(:,:,:)
83  real(RP), private, allocatable :: lambda (:,:,:)
84 
85  real(RP), private :: atmos_phy_tb_smg_nu_max = 10000.0_rp
86  logical, private :: atmos_phy_tb_smg_backscatter = .false.
87  logical, private :: atmos_phy_tb_smg_bottom = .true.
88  logical, private :: atmos_phy_tb_smg_implicit = .false.
89  logical, private :: atmos_phy_tb_smg_horizontal = .false.
90 
91  real(RP), private :: tke_fact
92 
93  !-----------------------------------------------------------------------------
94 contains
95  !-----------------------------------------------------------------------------
97  subroutine atmos_phy_tb_smg_setup( &
98  FZ, CZ, CDX, CDY, MAPF )
99  use scale_prc, only: &
100  prc_abort
101  use scale_const, only: &
102  karman => const_karman
103  implicit none
104 
105  real(RP), intent(in) :: FZ (0:ka,ia,ja)
106  real(RP), intent(in) :: CZ ( ka,ia,ja)
107  real(RP), intent(in) :: CDX (ia)
108  real(RP), intent(in) :: CDY (ja)
109  real(RP), intent(in) :: MAPF(ia,ja,2)
110 
111  real(RP) :: ATMOS_PHY_TB_SMG_Cs
112  real(RP) :: ATMOS_PHY_TB_SMG_filter_fact = 2.0_rp
113  logical :: ATMOS_PHY_TB_SMG_consistent_tke = .true.
114 
115  namelist / param_atmos_phy_tb_smg / &
116  atmos_phy_tb_smg_cs, &
117  atmos_phy_tb_smg_nu_max, &
118  atmos_phy_tb_smg_filter_fact, &
119  atmos_phy_tb_smg_bottom, &
120  atmos_phy_tb_smg_backscatter, &
121  atmos_phy_tb_smg_implicit, &
122  atmos_phy_tb_smg_consistent_tke, &
123  atmos_phy_tb_smg_horizontal
124 
125  integer :: ierr
126  integer :: k, i, j
127  !---------------------------------------------------------------------------
128 
129  log_newline
130  log_info("ATMOS_PHY_TB_smg_setup",*) 'Setup'
131  log_info("ATMOS_PHY_TB_smg_setup",*) 'Smagorinsky-type Eddy Viscocity Model'
132 
133  atmos_phy_tb_smg_cs = cs
134 
135  !--- read namelist
136  rewind(io_fid_conf)
137  read(io_fid_conf,nml=param_atmos_phy_tb_smg,iostat=ierr)
138  if( ierr < 0 ) then !--- missing
139  log_info("ATMOS_PHY_TB_smg_setup",*) 'Not found namelist. Default used.'
140  elseif( ierr > 0 ) then !--- fatal error
141  log_error("ATMOS_PHY_TB_smg_setup",*) 'Not appropriate names in namelist PARAM_ATMOS_PHY_TB_SMG. Check!'
142  call prc_abort
143  endif
144  log_nml(param_atmos_phy_tb_smg)
145 
146  cs = atmos_phy_tb_smg_cs
147 
148  rprn = 1.0_rp / prn
149  rric = 1.0_rp / ric
150  prnovric = ( 1.0_rp - prn ) * rric
151 
152  allocate( lambda0(ka,ia,ja) )
153  allocate( lambda(ka,ia,ja) )
154 
155 #ifdef DEBUG
156  lambda0(:,:,:) = undef
157  lambda(:,:,:) = undef
158 #endif
159  if ( atmos_phy_tb_smg_horizontal ) then
160  !$omp parallel do
161  do j = js-1, je+1
162  do i = is-1, ie+1
163  do k = ks, ke
164  lambda0(k,i,j) = cs * sqrt( cdx(i) * cdy(j) / ( mapf(i,j,1) * mapf(i,j,2) ) )
165  lambda(k,i,j) = lambda0(k,i,j)
166  enddo
167  enddo
168  enddo
169 #ifdef DEBUG
170  i = iundef; j = iundef; k = iundef
171 #endif
172  atmos_phy_tb_smg_consistent_tke = .false.
173  atmos_phy_tb_smg_implicit = .false. ! flux in the z-direction is not necessary
174  atmos_phy_tb_smg_backscatter = .false.
175  else
176  !$omp parallel do
177  do j = js-1, je+1
178  do i = is-1, ie+1
179  do k = ks, ke
180  lambda0(k,i,j) = cs * mixlen( fz(k,i,j) - fz(k-1,i,j), &
181  cdx(i) / mapf(i,j,1), &
182  cdy(j) / mapf(i,j,2), &
183  atmos_phy_tb_smg_filter_fact )
184  enddo
185  enddo
186  enddo
187 #ifdef DEBUG
188  i = iundef; j = iundef; k = iundef
189 #endif
190 
191  if ( atmos_phy_tb_smg_bottom ) then
192  !$omp parallel do
193  do j = js-1, je+1
194  do i = is-1, ie+1
195  do k = ks, ke
196  lambda(k,i,j) = sqrt( 1.0_rp / ( 1.0_rp / lambda0(k,i,j)**2 + 1.0_rp / ( karman*(cz(k,i,j)-fz(ks-1,i,j)) )**2 ) )
197  enddo
198  enddo
199  enddo
200 #ifdef DEBUG
201  i = iundef; j = iundef; k = iundef
202 #endif
203  else
204  !$omp parallel do
205  do j = js-1, je+1
206  do i = is-1, ie+1
207  do k = ks, ke
208  lambda(k,i,j) = lambda0(k,i,j)
209  enddo
210  enddo
211  enddo
212 #ifdef DEBUG
213  i = iundef; j = iundef; k = iundef
214 #endif
215  end if
216  end if
217 
218  ! flag for isotropic stress tensor
219  if ( atmos_phy_tb_smg_consistent_tke ) then
220  tke_fact = 1.0_rp
221  else
222  tke_fact = 0.0_rp ! neglect
223  end if
224 
225  return
226  end subroutine atmos_phy_tb_smg_setup
227 
228  !-----------------------------------------------------------------------------
229  subroutine atmos_phy_tb_smg( &
230  qflx_sgs_momz, qflx_sgs_momx, qflx_sgs_momy, &
231  qflx_sgs_rhot, qflx_sgs_rhoq, &
232  MOMZ_t, MOMX_t, MOMY_t, RHOT_t, RHOQ_t, &
233  Nu, Ri, Pr, &
234  MOMZ, MOMX, MOMY, POTT, DENS, QTRC, N2, &
235  FZ, FDZ, RCDZ, RFDZ, CDX, FDX, CDY, FDY, &
236  GSQRT, J13G, J23G, J33G, MAPF, dt )
237  use scale_precision
238  use scale_atmos_grid_cartesc_index
239  use scale_tracer
240  use scale_const, only: &
241  eps => const_eps, &
242  grav => const_grav
243  use scale_file_history, only: &
244  file_history_in
245  use scale_atmos_phy_tb_common, only: &
246  calc_strain_tensor => atmos_phy_tb_calc_strain_tensor, &
247  diffusion_solver => atmos_phy_tb_diffusion_solver, &
248  calc_tend_momz => atmos_phy_tb_calc_tend_momz, &
249  calc_tend_momx => atmos_phy_tb_calc_tend_momx, &
250  calc_tend_momy => atmos_phy_tb_calc_tend_momy, &
251  calc_flux_phi => atmos_phy_tb_calc_flux_phi
252  use scale_atmos_hydrometeor, only: &
253  i_qv
254  use scale_random, only: &
255  random_normal
256  implicit none
257 
258  ! SGS flux
259  real(RP), intent(out) :: qflx_sgs_momz(ka,ia,ja,3)
260  real(RP), intent(out) :: qflx_sgs_momx(ka,ia,ja,3)
261  real(RP), intent(out) :: qflx_sgs_momy(ka,ia,ja,3)
262  real(RP), intent(out) :: qflx_sgs_rhot(ka,ia,ja,3)
263  real(RP), intent(out) :: qflx_sgs_rhoq(ka,ia,ja,3,qa)
264 
265  real(RP), intent(out) :: MOMZ_t (ka,ia,ja)
266  real(RP), intent(out) :: MOMX_t (ka,ia,ja)
267  real(RP), intent(out) :: MOMY_t (ka,ia,ja)
268  real(RP), intent(out) :: RHOT_t (ka,ia,ja)
269  real(RP), intent(out) :: RHOQ_t (ka,ia,ja,qa) ! tendency of rho * QTRC
270 
271  real(RP), intent(out) :: nu (ka,ia,ja) ! eddy viscosity (center)
272  real(RP), intent(out) :: Ri (ka,ia,ja) ! Richardson number
273  real(RP), intent(out) :: Pr (ka,ia,ja) ! Prantle number
274 
275  real(RP), intent(in) :: MOMZ (ka,ia,ja)
276  real(RP), intent(in) :: MOMX (ka,ia,ja)
277  real(RP), intent(in) :: MOMY (ka,ia,ja)
278  real(RP), intent(in) :: POTT (ka,ia,ja)
279  real(RP), intent(in) :: DENS (ka,ia,ja)
280  real(RP), intent(in) :: QTRC (ka,ia,ja,qa)
281  real(RP), intent(in) :: N2 (ka,ia,ja)
282 
283  real(RP), intent(in) :: FZ (0:ka,ia,ja)
284  real(RP), intent(in) :: FDZ (ka-1)
285  real(RP), intent(in) :: RCDZ (ka)
286  real(RP), intent(in) :: RFDZ (ka-1)
287  real(RP), intent(in) :: CDX (ia)
288  real(RP), intent(in) :: FDX (ia-1)
289  real(RP), intent(in) :: CDY (ja)
290  real(RP), intent(in) :: FDY (ja-1)
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  ! tke
300  real(RP) :: TKE(ka,ia,ja)
301 
302  ! deformation rate tensor
303  real(RP) :: S33_C(ka,ia,ja) ! (cell center)
304  real(RP) :: S11_C(ka,ia,ja)
305  real(RP) :: S22_C(ka,ia,ja)
306  real(RP) :: S31_C(ka,ia,ja)
307  real(RP) :: S12_C(ka,ia,ja)
308  real(RP) :: S23_C(ka,ia,ja)
309  real(RP) :: S12_Z(ka,ia,ja) ! (z edge or x-y plane)
310  real(RP) :: S23_X(ka,ia,ja) ! (x edge or y-z plane)
311  real(RP) :: S31_Y(ka,ia,ja) ! (y edge or z-x plane)
312  real(RP) :: S2 (ka,ia,ja) ! |S|^2
313 
314  real(RP) :: Kh(ka,ia,ja) ! eddy diffusion
315  real(RP) :: fm(ka)
316  real(RP) :: e(ka)
317  real(RP) :: et
318  real(RP) :: lambda_r(ka)
319  real(RP) :: Rf
320  real(RP) :: C1(ka)
321  real(RP) :: C2
322  real(RP) :: D2
323 
324  ! implicit scheme
325  real(RP) :: TEND(ka,ia,ja)
326  real(RP) :: a (ka,ia,ja)
327  real(RP) :: b (ka,ia,ja)
328  real(RP) :: c (ka,ia,ja)
329  real(RP) :: d (ka)
330  real(RP) :: ap
331 
332  ! backscatter
333  real(RP) :: random (ka,ia,ja)
334  real(RP) :: random_mz(ka,ia,ja)
335  real(RP) :: random_mx(ka,ia,ja)
336  real(RP) :: random_my(ka,ia,ja)
337  real(RP) :: random_qz(ka,ia,ja)
338  real(RP) :: random_qx(ka,ia,ja)
339  real(RP) :: random_qy(ka,ia,ja)
340  real(RP) :: acc (ka,ia,ja)
341  real(RP) :: dd (ka,ia,ja)
342  real(RP) :: leOvleo5
343  real(RP) :: dz, dx, dy
344  real(RP) :: fact
345 
346  integer :: IIS, IIE
347  integer :: JJS, JJE
348 
349  integer :: k, i, j, iq
350  !---------------------------------------------------------------------------
351 
352  log_progress(*) 'atmosphere / physics / turbulence / Smagorinsky'
353 
354 #ifdef DEBUG
355  qflx_sgs_momz(:,:,:,:) = undef
356  qflx_sgs_momx(:,:,:,:) = undef
357  qflx_sgs_momy(:,:,:,:) = undef
358  qflx_sgs_rhot(:,:,:,:) = undef
359  qflx_sgs_rhoq(:,:,:,:,:) = undef
360 
361  nu(:,:,:) = undef
362  tke(:,:,:) = undef
363  pr(:,:,:) = undef
364  ri(:,:,:) = undef
365  kh(:,:,:) = undef
366 #endif
367 
368 #ifdef QUICKDEBUG
369  qflx_sgs_momz(ks:ke, 1:is-1, : ,:) = undef
370  qflx_sgs_momz(ks:ke,ie+1:ia , : ,:) = undef
371  qflx_sgs_momz(ks:ke, : , 1:js-1,:) = undef
372  qflx_sgs_momz(ks:ke, : ,je+1:ja ,:) = undef
373  qflx_sgs_momx(ks:ke, 1:is-1, : ,:) = undef
374  qflx_sgs_momx(ks:ke,ie+1:ia , : ,:) = undef
375  qflx_sgs_momx(ks:ke, : , 1:js-1,:) = undef
376  qflx_sgs_momx(ks:ke, : ,je+1:ja ,:) = undef
377  qflx_sgs_momy(ks:ke, 1:is-1, : ,:) = undef
378  qflx_sgs_momy(ks:ke,ie+1:ia , : ,:) = undef
379  qflx_sgs_momy(ks:ke, : , 1:js-1,:) = undef
380  qflx_sgs_momy(ks:ke, : ,je+1:ja ,:) = undef
381 #endif
382 
383 
384  !##### Start Upadate #####
385 
386  call calc_strain_tensor( &
387  s33_c, s11_c, s22_c, & ! (out)
388  s31_c, s12_c, s23_c, & ! (out)
389  s12_z, s23_x, s31_y, & ! (out)
390  s2 , & ! (out)
391  dens, momz, momx, momy, & ! (in)
392  gsqrt, j13g, j23g, j33g, mapf ) ! (in)
393 
394  if ( atmos_phy_tb_smg_backscatter ) then
395 
396  call random_normal( random(:,:,:) )
397  ! 1:2:1 filter
398  !$omp parallel do
399  do j = js-1, je+1
400  do i = is-1, ie+1
401  do k = ks+1, ke-1
402  random_mz(k,i,j) = ( random(k,i,j) * 2.0_rp &
403  + random(k+1,i,j) + random(k-1,i,j) &
404  + random(k,i+1,j) + random(k,i-1,j) &
405  + random(k,i,j+1) + random(k,i,j-1) ) / 8.0_rp
406  end do
407  random_mz(ks,i,j) = ( random(ks,i,j) * 2.0_rp &
408  + random(ks+1,i,j) &
409  + random(ks,i+1,j) + random(ks,i-1,j) &
410  + random(ks,i,j+1) + random(ks,i,j-1) ) / 7.0_rp
411  random_mz(ke,i,j) = ( random(ke,i,j) * 2.0_rp &
412  + random(ke-1,i,j) &
413  + random(ke,i+1,j) + random(ke,i-1,j) &
414  + random(ke,i,j+1) + random(ke,i,j-1) ) / 7.0_rp
415  end do
416  end do
417 
418  call random_normal( random(:,:,:) )
419  ! 1:2:1 filter
420  !$omp parallel do
421  do j = js-1, je+1
422  do i = is-1, ie+1
423  do k = ks+1, ke-1
424  random_mx(k,i,j) = ( random(k,i,j) * 2.0_rp &
425  + random(k+1,i,j) + random(k-1,i,j) &
426  + random(k,i+1,j) + random(k,i-1,j) &
427  + random(k,i,j+1) + random(k,i,j-1) ) / 8.0_rp
428  end do
429  random_mx(ks,i,j) = ( random(ks,i,j) * 2.0_rp &
430  + random(ks+1,i,j) &
431  + random(ks,i+1,j) + random(ks,i-1,j) &
432  + random(ks,i,j+1) + random(ks,i,j-1) ) / 7.0_rp
433  random_mx(ke,i,j) = ( random(ke,i,j) * 2.0_rp &
434  + random(ke-1,i,j) &
435  + random(ke,i+1,j) + random(ke,i-1,j) &
436  + random(ke,i,j+1) + random(ke,i,j-1) ) / 7.0_rp
437  end do
438  end do
439 
440  call random_normal( random(:,:,:) )
441  ! 1:2:1 filter
442  !$omp parallel do
443  do j = js-1, je+1
444  do i = is-1, ie+1
445  do k = ks+1, ke-1
446  random_my(k,i,j) = ( random(k,i,j) * 2.0_rp &
447  + random(k+1,i,j) + random(k-1,i,j) &
448  + random(k,i+1,j) + random(k,i-1,j) &
449  + random(k,i,j+1) + random(k,i,j-1) ) / 8.0_rp
450  end do
451  random_my(ks,i,j) = ( random(ks,i,j) * 2.0_rp &
452  + random(ks+1,i,j) &
453  + random(ks,i+1,j) + random(ks,i-1,j) &
454  + random(ks,i,j+1) + random(ks,i,j-1) ) / 7.0_rp
455  random_my(ke,i,j) = ( random(ke,i,j) * 2.0_rp &
456  + random(ke-1,i,j) &
457  + random(ke,i+1,j) + random(ke,i-1,j) &
458  + random(ke,i,j+1) + random(ke,i,j-1) ) / 7.0_rp
459  end do
460  end do
461 
462  call random_normal( random(:,:,:) )
463  ! 1:2:1 filter
464  !$omp parallel do
465  do j = js-1, je+1
466  do i = is-1, ie+1
467  do k = ks+1, ke-1
468  random_qz(k,i,j) = ( random(k,i,j) * 2.0_rp &
469  + random(k+1,i,j) + random(k-1,i,j) &
470  + random(k,i+1,j) + random(k,i-1,j) &
471  + random(k,i,j+1) + random(k,i,j-1) ) / 8.0_rp
472  end do
473  random_qz(ks,i,j) = ( random(ks,i,j) * 2.0_rp &
474  + random(ks+1,i,j) &
475  + random(ks,i+1,j) + random(ks,i-1,j) &
476  + random(ks,i,j+1) + random(ks,i,j-1) ) / 7.0_rp
477  random_qz(ke,i,j) = ( random(ke,i,j) * 2.0_rp &
478  + random(ke-1,i,j) &
479  + random(ke,i+1,j) + random(ke,i-1,j) &
480  + random(ke,i,j+1) + random(ke,i,j-1) ) / 7.0_rp
481  end do
482  end do
483 
484  call random_normal( random(:,:,:) )
485  ! 1:2:1 filter
486  !$omp parallel do
487  do j = js-1, je+1
488  do i = is-1, ie+1
489  do k = ks+1, ke-1
490  random_qx(k,i,j) = ( random(k,i,j) * 2.0_rp &
491  + random(k+1,i,j) + random(k-1,i,j) &
492  + random(k,i+1,j) + random(k,i-1,j) &
493  + random(k,i,j+1) + random(k,i,j-1) ) / 8.0_rp
494  end do
495  random_qx(ks,i,j) = ( random(ks,i,j) * 2.0_rp &
496  + random(ks+1,i,j) &
497  + random(ks,i+1,j) + random(ks,i-1,j) &
498  + random(ks,i,j+1) + random(ks,i,j-1) ) / 7.0_rp
499  random_qx(ke,i,j) = ( random(ke,i,j) * 2.0_rp &
500  + random(ke-1,i,j) &
501  + random(ke,i+1,j) + random(ke,i-1,j) &
502  + random(ke,i,j+1) + random(ke,i,j-1) ) / 7.0_rp
503  end do
504  end do
505 
506  call random_normal( random(:,:,:) )
507  ! 1:2:1 filter
508  !$omp parallel do
509  do j = js-1, je+1
510  do i = is-1, ie+1
511  do k = ks+1, ke-1
512  random_qy(k,i,j) = ( random(k,i,j) * 2.0_rp &
513  + random(k+1,i,j) + random(k-1,i,j) &
514  + random(k,i+1,j) + random(k,i-1,j) &
515  + random(k,i,j+1) + random(k,i,j-1) ) / 8.0_rp
516  end do
517  random_qy(ks,i,j) = ( random(ks,i,j) * 2.0_rp &
518  + random(ks+1,i,j) &
519  + random(ks,i+1,j) + random(ks,i-1,j) &
520  + random(ks,i,j+1) + random(ks,i,j-1) ) / 7.0_rp
521  random_qy(ke,i,j) = ( random(ke,i,j) * 2.0_rp &
522  + random(ke-1,i,j) &
523  + random(ke,i+1,j) + random(ke,i-1,j) &
524  + random(ke,i,j+1) + random(ke,i,j-1) ) / 7.0_rp
525  end do
526  end do
527 
528 
529  end if
530 
531 
532  !$omp parallel do &
533  !$omp private(fm,Rf,lambda_r,leOvleo5,C1,C2,D2,e,et,fact,dz,dx,dy)
534  do j = js-1, je+1
535  do i = is-1, ie+1
536 
537  ! Ri = N^2 / |S|^2, N^2 = g / theta * dtheta/dz
538  do k = ks, ke
539  ri(k,i,j) = n2(k,i,j) / max(s2(k,i,j), eps)
540  enddo
541 
542  ! Nu
543  ! Pr = Nu / Kh = fm / fh
544  do k = ks, ke
545  if ( ri(k,i,j) < 0.0_rp ) then ! stable
546  fm(k) = sqrt( 1.0_rp - fmc * ri(k,i,j) )
547  nu(k,i,j) = lambda(k,i,j)**2 * sqrt( s2(k,i,j) ) * fm(k)
548  pr(k,i,j) = fm(k) / sqrt( 1.0_rp - fhb*ri(k,i,j) ) * prn
549  else if ( ri(k,i,j) < ric ) then ! weakly stable
550  fm(k) = ( 1.0_rp - ri(k,i,j)*rric )**4
551  nu(k,i,j) = lambda(k,i,j)**2 * sqrt( s2(k,i,j) ) * fm(k)
552  pr(k,i,j) = prn / ( 1.0_rp - prnovric * ri(k,i,j) )
553  else ! strongly stable
554  fm(k) = 0.0_rp
555  nu(k,i,j) = 0.0_rp
556  kh(k,i,j) = 0.0_rp
557  pr(k,i,j) = 1.0_rp
558  endif
559 
560  if ( ri(k,i,j) < ric ) then
561  kh(k,i,j) = max( min( nu(k,i,j) / pr(k,i,j), atmos_phy_tb_smg_nu_max ), eps )
562  nu(k,i,j) = max( min( nu(k,i,j), atmos_phy_tb_smg_nu_max ), eps )
563  pr(k,i,j) = nu(k,i,j) / kh(k,i,j)
564  rf = ri(k,i,j) / pr(k,i,j)
565  lambda_r(k) = lambda(k,i,j) * sqrt( fm(k) / sqrt( 1.0_rp - rf ) )
566  else
567  lambda_r(k) = 0.0_rp
568  end if
569 
570  enddo
571 
572  if ( atmos_phy_tb_smg_backscatter ) then
573 
574  do k = ks, ke
575  lambda_r(k) = min( 1.8_rp * lambda0(k,i,j), lambda_r(k) )
576  leovleo5 = ( lambda_r(k) / lambda0(k,i,j) )**5
577  c2 = cb * leovleo5 / ( 1.0_rp + cb * leovleo5 )
578  d2 = c2 * cbt / cb
579  c1(k) = c1o / sqrt( 1.0_rp - c2 )
580  ! D1 = D1o / sqrt( 1.0_RP - C2 )
581  e(k) = nu(k,i,j)**3 * ( 1.0_rp - c2 ) / ( lambda_r(k)**4 + eps )
582  et = kh(k,i,j) * ( 1.0_rp - d2 )
583 
584  dz = fz(k,i,j) - fz(k-1,i,j)
585  dx = cdx(i) / mapf(i,j,1,i_xy)
586  dy = cdy(j) / mapf(i,j,2,i_xy)
587 
588  fact = sqrt( cb * leovleo5 * e(k) / ( dt * ( dz**2 + dx**2 + dy**2 ) ) ) * dens(k,i,j)
589  random_mz(k,i,j) = random_mz(k,i,j) * fact * dx * dy
590  random_mx(k,i,j) = random_mx(k,i,j) * fact * dy * dz
591  random_my(k,i,j) = random_my(k,i,j) * fact * dz * dx
592 
593  fact = sqrt( twooverthree * cbt * leovleo5 * et / dt ) * dens(k,i,j)
594  random_qz(k,i,j) = random_qz(k,i,j) * fact * dz
595  random_qx(k,i,j) = random_qx(k,i,j) * fact * dx
596  random_qy(k,i,j) = random_qy(k,i,j) * fact * dy
597  end do
598 
599  else
600 
601  do k = ks, ke
602  e(k) = nu(k,i,j)**3 / ( lambda_r(k)**4 + eps )
603  c1(k) = c1o
604  end do
605 
606  end if
607 
608  ! TKE
609  do k = ks, ke
610  tke(k,i,j) = ( e(k) * lambda_r(k) / c1(k) )**twooverthree
611  enddo
612 
613  enddo
614  enddo
615 #ifdef DEBUG
616  i = iundef; j = iundef; k = iundef
617 #endif
618 
619 
620 
621  do jjs = js, je, jblock
622  jje = jjs+jblock-1
623  do iis = is, ie, iblock
624  iie = iis+iblock-1
625 
626  !##### momentum equation (z) #####
627  ! (cell center)
628  if ( atmos_phy_tb_smg_horizontal ) then
629  qflx_sgs_momz(:,:,:,zdir) = 0.0_rp
630  else
631  !$omp parallel do
632  do j = jjs, jje
633  do i = iis, iie
634  do k = ks+1, ke-1
635 #ifdef DEBUG
636  call check( __line__, dens(k,i,j) )
637  call check( __line__, nu(k,i,j) )
638  call check( __line__, s33_c(k,i,j) )
639  call check( __line__, s11_c(k,i,j) )
640  call check( __line__, s22_c(k,i,j) )
641  call check( __line__, tke(k,i,j) )
642 #endif
643  qflx_sgs_momz(k,i,j,zdir) = dens(k,i,j) * ( &
644  - 2.0_rp * nu(k,i,j) &
645  * ( s33_c(k,i,j) - ( s11_c(k,i,j) + s22_c(k,i,j) + s33_c(k,i,j) ) * oneoverthree ) &
646  + twooverthree * tke(k,i,j) * tke_fact )
647  enddo
648  enddo
649  enddo
650 #ifdef DEBUG
651  i = iundef; j = iundef; k = iundef
652 #endif
653  !$omp parallel do
654  do j = jjs, jje
655  do i = iis, iie
656  ! momentum will not be conserved
657  qflx_sgs_momz(ks,i,j,zdir) = dens(ks,i,j) * twooverthree * tke(ks,i,j) * tke_fact
658  qflx_sgs_momz(ke,i,j,zdir) = dens(ke,i,j) * twooverthree * tke(ke,i,j) * tke_fact
659  ! anti-isotropic stress is calculated by the surface scheme
660  enddo
661  enddo
662 #ifdef DEBUG
663  i = iundef; j = iundef; k = iundef
664 #endif
665  end if
666  ! (y edge)
667  !$omp parallel do
668  do j = jjs, jje
669  do i = iis-1, iie
670  do k = ks, ke-1
671 #ifdef DEBUG
672  call check( __line__, dens(k,i,j) )
673  call check( __line__, dens(k,i+1,j) )
674  call check( __line__, dens(k+1,i,j) )
675  call check( __line__, dens(k+1,i+1,j) )
676  call check( __line__, nu(k,i,j) )
677  call check( __line__, nu(k,i+1,j) )
678  call check( __line__, nu(k+1,i,j) )
679  call check( __line__, nu(k+1,i+1,j) )
680  call check( __line__, s31_y(k,i,j) )
681 #endif
682  qflx_sgs_momz(k,i,j,xdir) = - 0.125_rp & ! 2.0 / 4 / 4
683  * ( dens(k,i,j)+dens(k+1,i,j)+dens(k,i+1,j)+dens(k+1,i+1,j) ) &
684  * ( nu(k,i,j)+nu(k+1,i,j)+nu(k,i+1,j)+nu(k+1,i+1,j)) &
685  * s31_y(k,i,j)
686  enddo
687  enddo
688  enddo
689 #ifdef DEBUG
690  i = iundef; j = iundef; k = iundef
691 #endif
692  ! (x edge)
693  !$omp parallel do default(none) &
694  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,DENS,S23_X,nu,qflx_sgs_momz) &
695  !$omp private(i,j,k) OMP_SCHEDULE_ collapse(2)
696  do j = jjs-1, jje
697  do i = iis, iie
698  do k = ks, ke-1
699 #ifdef DEBUG
700  call check( __line__, dens(k,i,j) )
701  call check( __line__, dens(k,i,j+1) )
702  call check( __line__, dens(k+1,i,j) )
703  call check( __line__, dens(k+1,i,j+1) )
704  call check( __line__, nu(k,i,j) )
705  call check( __line__, nu(k,i,j+1) )
706  call check( __line__, nu(k+1,i,j) )
707  call check( __line__, nu(k+1,i,j+1) )
708  call check( __line__, s23_x(k,i,j) )
709 #endif
710  qflx_sgs_momz(k,i,j,ydir) = - 0.125_rp & ! 2/4/4
711  * ( dens(k,i,j)+dens(k+1,i,j)+dens(k,i,j+1)+dens(k+1,i,j+1) ) &
712  * ( nu(k,i,j)+nu(k+1,i,j)+nu(k,i,j+1)+nu(k+1,i,j+1) ) &
713  * s23_x(k,i,j)
714  enddo
715  enddo
716  enddo
717 #ifdef DEBUG
718  i = iundef; j = iundef; k = iundef
719 #endif
720 
721  if ( atmos_phy_tb_smg_implicit ) then
722 
723  call calc_tend_momz( tend, & ! (out)
724  qflx_sgs_momz, & ! (in)
725  gsqrt, j13g, j23g, j33g, mapf, & ! (in)
726  iis, iie, jjs, jje ) ! (in)
727 
728  !$omp parallel do &
729  !$omp private (ap,d)
730  do j = jjs, jje
731  do i = iis, iie
732 
733  ap = - fouroverthree * dt &
734  * dens(ks+1,i,j)*nu(ks+1,i,j) &
735  * rcdz(ks+1) / gsqrt(ks+1,i,j,i_xyz)
736  a(ks,i,j) = ap * rfdz(ks) / gsqrt(ks,i,j,i_xyw)
737  c(ks,i,j) = 0.0_rp
738  b(ks,i,j) = - a(ks,i,j) + 0.5_rp * ( dens(ks,i,j)+dens(ks+1,i,j) )
739  do k = ks+1, ke-2
740  c(k,i,j) = ap * rfdz(k+1) / gsqrt(k+1,i,j,i_xyw)
741  ap = - fouroverthree * dt &
742  * dens(k+1,i,j)*nu(k+1,i,j) &
743  * rcdz(k+1) / gsqrt(k+1,i,j,i_xyz)
744  a(k,i,j) = ap * rfdz(k) / gsqrt(k,i,j,i_xyw)
745  b(k,i,j) = - a(k,i,j) - c(k,i,j) + 0.5_rp * ( dens(k,i,j)+dens(k+1,i,j) )
746  end do
747  a(ke-1,i,j) = 0.0_rp
748  c(ke-1,i,j) = ap * rfdz(ke) / gsqrt(ke,i,j,i_xyw)
749  b(ke-1,i,j) = - c(ke-1,i,j) + 0.5_rp * ( dens(ke-1,i,j)+dens(ke,i,j) )
750 
751  do k = ks, ke-1
752  d(k) = tend(k,i,j)
753  end do
754 
755  call diffusion_solver( &
756  tend(:,i,j), & ! (out)
757  a(:,i,j), b(:,i,j), c(:,i,j), d, & ! (in)
758  ke-1 ) ! (in)
759 
760  do k = ks+1, ke-1
761  qflx_sgs_momz(k,i,j,zdir) = qflx_sgs_momz(k,i,j,zdir) &
762  - fouroverthree * dens(k,i,j) * nu(k,i,j) * dt &
763  * ( tend(k,i,j) - tend(k-1,i,j) ) * rcdz(k) / gsqrt(k,i,j,i_xyz)
764  end do
765 
766  end do
767  end do
768 
769  end if
770 
771  !##### momentum equation (x) #####
772  ! (y edge)
773  if ( atmos_phy_tb_smg_horizontal ) then
774  qflx_sgs_momx(:,:,:,zdir) = 0.0_rp
775  else
776  !$omp parallel do
777  do j = jjs, jje
778  do i = iis, iie
779  do k = ks, ke-1
780 #ifdef DEBUG
781  call check( __line__, dens(k,i,j) )
782  call check( __line__, dens(k,i+1,j) )
783  call check( __line__, dens(k+1,i,j) )
784  call check( __line__, dens(k+1,i+1,j) )
785  call check( __line__, nu(k,i,j) )
786  call check( __line__, nu(k,i+1,j) )
787  call check( __line__, nu(k+1,i,j) )
788  call check( __line__, nu(k+1,i+1,j) )
789  call check( __line__, s31_y(k,i,j) )
790 #endif
791  qflx_sgs_momx(k,i,j,zdir) = - 0.125_rp & ! 2/4/4
792  * ( dens(k,i,j)+dens(k+1,i,j)+dens(k,i+1,j)+dens(k+1,i+1,j) ) &
793  * ( nu(k,i,j)+nu(k+1,i,j)+nu(k,i+1,j)+nu(k+1,i+1,j) ) &
794  * s31_y(k,i,j)
795  enddo
796  enddo
797  enddo
798 #ifdef DEBUG
799  i = iundef; j = iundef; k = iundef
800 #endif
801  !$omp parallel do
802  do j = jjs, jje
803  do i = iis, iie
804  qflx_sgs_momx(ks-1,i,j,zdir) = 0.0_rp ! bottom boundary
805  qflx_sgs_momx(ke ,i,j,zdir) = 0.0_rp ! top boundary
806  enddo
807  enddo
808 #ifdef DEBUG
809  i = iundef; j = iundef; k = iundef
810 #endif
811  end if
812  ! (cell center)
813  !$omp parallel do default(none) &
814  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,DENS,nu,S11_C,S22_C,S33_C,TKE,tke_fact,qflx_sgs_momx) &
815  !$omp private(i,j,k) OMP_SCHEDULE_ collapse(2)
816  do j = jjs, jje
817  do i = iis, iie+1
818  do k = ks, ke
819 #ifdef DEBUG
820  call check( __line__, dens(k,i,j) )
821  call check( __line__, nu(k,i,j) )
822  call check( __line__, s11_c(k,i,j) )
823  call check( __line__, s22_c(k,i,j) )
824  call check( __line__, s33_c(k,i,j) )
825  call check( __line__, tke(k,i,j) )
826 #endif
827  qflx_sgs_momx(k,i,j,xdir) = dens(k,i,j) * ( &
828  - 2.0_rp * nu(k,i,j) &
829  * ( s11_c(k,i,j) - ( s11_c(k,i,j) + s22_c(k,i,j) + s33_c(k,i,j) ) * oneoverthree ) &
830  + twooverthree * tke(k,i,j) * tke_fact )
831  enddo
832  enddo
833  enddo
834 #ifdef DEBUG
835  i = iundef; j = iundef; k = iundef
836 #endif
837  ! (z edge)
838  !$omp parallel do default(none) &
839  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,DENS,nu,S12_Z,qflx_sgs_momx) &
840  !$omp private(i,j,k) OMP_SCHEDULE_ collapse(2)
841  do j = jjs-1, jje
842  do i = iis, iie
843  do k = ks, ke
844 #ifdef DEBUG
845  call check( __line__, dens(k,i,j) )
846  call check( __line__, dens(k,i+1,j) )
847  call check( __line__, dens(k,i,j+1) )
848  call check( __line__, dens(k,i+1,j+1) )
849  call check( __line__, nu(k,i,j) )
850  call check( __line__, nu(k,i+1,j) )
851  call check( __line__, nu(k,i,j+1) )
852  call check( __line__, nu(k,i+1,j+1) )
853  call check( __line__, s12_z(k,i,j) )
854 #endif
855  qflx_sgs_momx(k,i,j,ydir) = - 0.125_rp & ! 2/4/4
856  * ( dens(k,i,j)+dens(k,i+1,j)+dens(k,i,j+1)+dens(k,i+1,j+1) ) &
857  * ( nu(k,i,j)+nu(k,i+1,j)+nu(k,i,j+1)+nu(k,i+1,j+1) ) &
858  * s12_z(k,i,j)
859  enddo
860  enddo
861  enddo
862 #ifdef DEBUG
863  i = iundef; j = iundef; k = iundef
864 #endif
865 
866  if ( atmos_phy_tb_smg_implicit ) then
867  call calc_tend_momx( tend, & ! (out)
868  qflx_sgs_momx, & ! (in)
869  gsqrt, j13g, j23g, j33g, mapf, & ! (in)
870  iis, iie, jjs, jje ) ! (in)
871 
872  !$omp parallel do &
873  !$omp private(ap,d)
874  do j = jjs, jje
875  do i = iis, iie
876 
877  ap = - dt * 0.25_rp * ( dens(ks ,i ,j)*nu(ks ,i ,j) &
878  + dens(ks+1,i ,j)*nu(ks+1,i ,j) &
879  + dens(ks ,i+1,j)*nu(ks ,i+1,j) &
880  + dens(ks+1,i+1,j)*nu(ks+1,i+1,j) ) &
881  * rfdz(ks) / gsqrt(ks,i,j,i_uyw)
882  a(ks,i,j) = ap * rcdz(ks) / gsqrt(ks,i,j,i_uyz)
883  c(ks,i,j) = 0.0_rp
884  b(ks,i,j) = - a(ks,i,j) + 0.5_rp * ( dens(ks,i,j)+dens(ks,i+1,j) )
885  do k = ks+1, ke-1
886  c(k,i,j) = ap * rcdz(k) / gsqrt(k,i,j,i_uyz)
887  ap = - dt * 0.25_rp * ( dens(k ,i ,j)*nu(k ,i ,j) &
888  + dens(k+1,i ,j)*nu(k+1,i ,j) &
889  + dens(k ,i+1,j)*nu(k ,i+1,j) &
890  + dens(k+1,i+1,j)*nu(k+1,i+1,j) ) &
891  * rfdz(k) / gsqrt(k,i,j,i_uyw)
892  a(k,i,j) = ap * rcdz(k) / gsqrt(k,i,j,i_uyz)
893  b(k,i,j) = - a(k,i,j) - c(k,i,j) + 0.5_rp * ( dens(k,i,j)+dens(k,i+1,j) )
894  end do
895  a(ke,i,j) = 0.0_rp
896  c(ke,i,j) = ap * rcdz(ke) / gsqrt(ke,i,j,i_uyz)
897  b(ke,i,j) = - c(ke,i,j) + 0.5_rp * ( dens(ke,i,j)+dens(ke,i+1,j) )
898 
899  do k = ks, ke
900  d(k) = tend(k,i,j)
901  end do
902 
903  call diffusion_solver( &
904  tend(:,i,j), & ! (out)
905  a(:,i,j), b(:,i,j), c(:,i,j), d, & ! (in)
906  ke ) ! (in)
907 
908  do k = ks, ke-1
909  qflx_sgs_momx(k,i,j,zdir) = qflx_sgs_momx(k,i,j,zdir) &
910  - 0.25_rp * ( dens(k ,i ,j)*nu(k ,i ,j) &
911  + dens(k+1,i ,j)*nu(k+1,i ,j) &
912  + dens(k ,i+1,j)*nu(k ,i+1,j) &
913  + dens(k+1,i+1,j)*nu(k+1,i+1,j) ) &
914  * dt * ( tend(k+1,i,j) - tend(k,i,j) ) * rfdz(k) / gsqrt(k,i,j,i_uyw)
915  end do
916 
917  end do
918  end do
919 
920  end if
921 
922  !##### momentum equation (y) #####
923  ! (x edge)
924  if ( atmos_phy_tb_smg_horizontal ) then
925  qflx_sgs_momy(:,:,:,zdir) = 0.0_rp
926  else
927  !$omp parallel do
928  do j = jjs, jje
929  do i = iis, iie
930  do k = ks, ke-1
931 #ifdef DEBUG
932  call check( __line__, dens(k,i,j) )
933  call check( __line__, dens(k,i,j+1) )
934  call check( __line__, dens(k+1,i,j) )
935  call check( __line__, dens(k+1,i,j+1) )
936  call check( __line__, nu(k,i,j) )
937  call check( __line__, nu(k,i,j+1) )
938  call check( __line__, nu(k+1,i,j) )
939  call check( __line__, nu(k+1,i,j+1) )
940  call check( __line__, s23_x(k,i,j) )
941 #endif
942  qflx_sgs_momy(k,i,j,zdir) = - 0.125_rp & ! 2/4/4
943  * ( dens(k,i,j)+dens(k+1,i,j)+dens(k,i,j+1)+dens(k+1,i,j+1) ) &
944  * ( nu(k,i,j)+nu(k+1,i,j)+nu(k,i,j+1)+nu(k+1,i,j+1) ) &
945  * s23_x(k,i,j)
946  enddo
947  enddo
948  enddo
949 #ifdef DEBUG
950  i = iundef; j = iundef; k = iundef
951 #endif
952  !$omp parallel do
953  do j = jjs, jje
954  do i = iis, iie
955  qflx_sgs_momy(ks-1,i,j,zdir) = 0.0_rp ! bottom boundary
956  qflx_sgs_momy(ke ,i,j,zdir) = 0.0_rp ! top boundary
957  enddo
958  enddo
959 #ifdef DEBUG
960  i = iundef; j = iundef; k = iundef
961 #endif
962  end if
963 
964  ! (z edge)
965  !$omp parallel do default(none) &
966  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,DENS,nu,S12_Z,qflx_sgs_momy) &
967  !$omp private(i,j,k) OMP_SCHEDULE_ collapse(2)
968  do j = jjs, jje
969  do i = iis-1, iie
970  do k = ks, ke
971 #ifdef DEBUG
972  call check( __line__, dens(k,i,j) )
973  call check( __line__, dens(k,i+1,j) )
974  call check( __line__, dens(k,i,j+1) )
975  call check( __line__, dens(k,i+1,j+1) )
976  call check( __line__, nu(k,i,j) )
977  call check( __line__, nu(k,i+1,j) )
978  call check( __line__, nu(k,i,j+1) )
979  call check( __line__, nu(k,i+1,j+1) )
980  call check( __line__, s12_z(k,i,j) )
981 #endif
982  qflx_sgs_momy(k,i,j,xdir) = - 0.125_rp & !
983  * ( dens(k,i,j)+dens(k,i+1,j)+dens(k,i,j+1)+dens(k,i+1,j+1) ) &
984  * ( nu(k,i,j)+nu(k,i+1,j)+nu(k,i,j+1)+nu(k,i+1,j+1) ) &
985  * s12_z(k,i,j)
986  enddo
987  enddo
988  enddo
989 #ifdef DEBUG
990  i = iundef; j = iundef; k = iundef
991 #endif
992 
993  ! (z-x plane)
994  !$omp parallel do default(none) &
995  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,DENS,nu,S11_C,S22_C,S33_C,tke_fact,TKE,qflx_sgs_momy) &
996  !$omp private(i,j,k) OMP_SCHEDULE_ collapse(2)
997  do j = jjs, jje+1
998  do i = iis, iie
999  do k = ks, ke
1000 #ifdef DEBUG
1001  call check( __line__, dens(k,i,j) )
1002  call check( __line__, nu(k,i,j) )
1003  call check( __line__, s11_c(k,i,j) )
1004  call check( __line__, s22_c(k,i,j) )
1005  call check( __line__, s33_c(k,i,j) )
1006  call check( __line__, tke(k,i,j) )
1007 #endif
1008  qflx_sgs_momy(k,i,j,ydir) = dens(k,i,j) * ( &
1009  - 2.0_rp * nu(k,i,j) &
1010  * ( s22_c(k,i,j) - ( s11_c(k,i,j) + s22_c(k,i,j) + s33_c(k,i,j) ) * oneoverthree ) &
1011  + twooverthree * tke(k,i,j) * tke_fact)
1012  enddo
1013  enddo
1014  enddo
1015 #ifdef DEBUG
1016  i = iundef; j = iundef; k = iundef
1017 #endif
1018 
1019  if ( atmos_phy_tb_smg_implicit ) then
1020  call calc_tend_momy( tend, & ! (out)
1021  qflx_sgs_momy, & ! (in)
1022  gsqrt, j13g, j23g, j33g, mapf, & ! (in)
1023  iis, iie, jjs, jje ) ! (in)
1024 
1025  !$omp parallel do &
1026  !$omp private(ap,d)
1027  do j = jjs, jje
1028  do i = iis, iie
1029 
1030  ap = - dt * 0.25_rp * ( dens(ks ,i,j )*nu(ks ,i,j ) &
1031  + dens(ks+1,i,j )*nu(ks+1,i,j ) &
1032  + dens(ks ,i,j+1)*nu(ks ,i,j+1) &
1033  + dens(ks+1,i,j+1)*nu(ks+1,i,j+1) ) &
1034  * rfdz(ks) / gsqrt(ks,i,j,i_xvw)
1035  a(ks,i,j) = ap * rcdz(ks) / gsqrt(ks,i,j,i_xvz)
1036  c(ks,i,j) = 0.0_rp
1037  b(ks,i,j) = - a(ks,i,j) + 0.5_rp * ( dens(ks,i,j)+dens(ks,i,j+1) )
1038  do k = ks+1, ke-1
1039  c(k,i,j) = ap * rcdz(k) / gsqrt(k,i,j,i_xvz)
1040  ap = - dt * 0.25_rp * ( dens(k ,i,j )*nu(k ,i,j ) &
1041  + dens(k+1,i,j )*nu(k+1,i,j ) &
1042  + dens(k ,i,j+1)*nu(k ,i,j+1) &
1043  + dens(k+1,i,j+1)*nu(k+1,i,j+1) ) &
1044  * rfdz(k) / gsqrt(k,i,j,i_xvw)
1045  a(k,i,j) = ap * rcdz(k) / gsqrt(k,i,j,i_xvz)
1046  b(k,i,j) = - a(k,i,j) - c(k,i,j) + 0.5_rp * ( dens(k,i,j)+dens(k,i,j+1) )
1047  end do
1048  a(ke,i,j) = 0.0_rp
1049  c(ke,i,j) = ap * rcdz(ke) / gsqrt(ke,i,j,i_xvz)
1050  b(ke,i,j) = - c(ke,i,j) + 0.5_rp * ( dens(ke,i,j)+dens(ke,i,j+1) )
1051 
1052  do k = ks, ke
1053  d(k) = tend(k,i,j)
1054  end do
1055 
1056  call diffusion_solver( &
1057  tend(:,i,j), & ! (out)
1058  a(:,i,j), b(:,i,j), c(:,i,j), d, & ! (in)
1059  ke ) ! (in)
1060 
1061  do k = ks, ke-1
1062  qflx_sgs_momy(k,i,j,zdir) = qflx_sgs_momy(k,i,j,zdir) &
1063  - 0.25_rp * ( dens(k ,i,j )*nu(k ,i,j ) &
1064  + dens(k+1,i,j )*nu(k+1,i,j ) &
1065  + dens(k ,i,j+1)*nu(k ,i,j+1) &
1066  + dens(k+1,i,j+1)*nu(k+1,i,j+1) ) &
1067  * dt * ( tend(k+1,i,j) - tend(k,i,j) ) * rfdz(k) / gsqrt(k,i,j,i_xvw)
1068  end do
1069 
1070  end do
1071  end do
1072 
1073  end if
1074 
1075  if ( atmos_phy_tb_smg_backscatter ) then
1076 
1077  ! MOMZ : dfy/dx - dfx/dy
1078  !$omp parallel do
1079  do j = jjs, jje
1080  do i = iis, iie
1081  do k = ks+1, ke-1
1082  acc(k,i,j) = ( ( gsqrt(k,i+1,j,i_xyz) * random_my(k,i+1,j) - gsqrt(k,i-1,j,i_xyz) * random_my(k,i-1,j) ) / ( fdx(i) + fdx(i-1) ) * mapf(i,j,1,i_xy) &
1083  + ( j13g(k+1,i,j,i_xyz) * random_my(k+1,i,j) - j13g(k-1,i,j,i_xyz) * random_my(k-1,i,j) ) / ( fdz(k) + fdz(k-1) ) &
1084  - ( gsqrt(k,i,j+1,i_xyz) * random_mx(k,i,j+1) - gsqrt(k,i,j-1,i_xyz) * random_mx(k,i,j-1) ) / ( fdy(j) + fdy(j-1) ) * mapf(i,j,2,i_xy) &
1085  + ( j23g(k+1,i,j,i_xyz) * random_mx(k+1,i,j) - j23g(k-1,i,j,i_xyz) * random_mx(k-1,i,j) ) / ( fdz(k) + fdz(k-1) ) &
1086  ) / gsqrt(k,i,j,i_xyz)
1087  end do
1088  acc(ks,i,j) = ( ( gsqrt(ks,i+1,j,i_xyz) * random_my(ks,i+1,j) - gsqrt(ks,i-1,j,i_xyz) * random_my(ks,i-1,j) ) / ( fdx(i) + fdx(i-1) ) * mapf(i,j,1,i_xy) &
1089  + ( j13g(ks+1,i,j,i_xyz) * random_my(ks+1,i,j) - j13g(ks,i,j,i_xyz) * random_my(ks,i,j) ) * rfdz(ks) * 0.5_rp &
1090  - ( gsqrt(ks,i,j+1,i_xyz) * random_mx(ks,i,j+1) - gsqrt(ks,i,j-1,i_xyz) * random_mx(ks,i,j-1) ) / ( fdy(j) + fdy(j-1) ) * mapf(i,j,2,i_xy) &
1091  + ( j23g(ks+1,i,j,i_xyz) * random_mx(ks+1,i,j) - j23g(ks,i,j,i_xyz) * random_mx(ks,i,j) ) * rfdz(ks) * 0.5_rp &
1092  ) / gsqrt(ks,i,j,i_xyz)
1093  acc(ke,i,j) = ( ( gsqrt(ke,i+1,j,i_xyz) * random_my(ke,i+1,j) - gsqrt(ke,i-1,j,i_xyz) * random_my(ke,i-1,j) ) / ( fdx(i) + fdx(i-1) ) * mapf(i,j,1,i_xy) &
1094  + ( j13g(ke,i,j,i_xyz) * random_my(ke,i,j) - j13g(ke-1,i,j,i_xyz) * random_my(ke-1,i,j) ) * rfdz(ke-1) * 0.5_rp &
1095  - ( gsqrt(ke,i,j+1,i_xyz) * random_mx(ke,i,j+1) - gsqrt(ke,i,j-1,i_xyz) * random_mx(ke,i,j-1) ) / ( fdy(j) + fdy(j-1) ) * mapf(i,j,2,i_xy) &
1096  + ( j23g(ke,i,j,i_xyz) * random_mx(ke,i,j) - j23g(ke-1,i,j,i_xyz) * random_mx(ke-1,i,j) ) * rfdz(ke-1) * 0.5_rp &
1097  ) / gsqrt(ke,i,j,i_xyz)
1098  end do
1099  end do
1100 
1101 #define f2h(k,i,j,p) ( ( FZ(k+p-1,i,j) - FZ(k+p-2,i,j) ) / ( FZ(k+1,i,j) - FZ(k-1,i,j) ) )
1102  !$omp parallel do
1103  do j = jjs, jje
1104  do i = iis, iie
1105  do k = ks, ke-1
1106  momz_t(k,i,j) = acc(k+1,i,j) * f2h(k,i,j,1) + acc(k,i,j) * f2h(k,i,j,2)
1107  end do
1108  momz_t(ke,i,j) = 0.0_rp
1109  end do
1110  end do
1111 
1112  ! MOMX : dfz/dy - dfy/dz
1113  !$omp parallel do
1114  do j = jjs, jje
1115  do i = iis, iie+1
1116  do k = ks+1, ke-1
1117  acc(k,i,j) = ( ( gsqrt(k,i,j+1,i_xyz) * random_mz(k,i,j+1) - gsqrt(k,i,j-1,i_xyz) * random_mz(k,i,j-1) ) / ( fdy(j) + fdy(j-1) ) * mapf(i,j,2,i_xy) &
1118  + ( j23g(k+1,i,j,i_xyz) * random_mz(k+1,i,j) - j23g(k-1,i,j,i_xyz) * random_mz(k-1,i,j) ) / ( fdz(k) + fdz(k-1) ) &
1119  - j33g * ( random_my(k+1,i,j) - random_my(k-1,i,j) ) / ( fdz(k) + fdz(k-1) ) &
1120  ) / gsqrt(k,i,j,i_xyz)
1121  end do
1122  acc(ks,i,j) = ( ( gsqrt(ks,i,j+1,i_xyz) * random_mz(ks,i,j+1) - gsqrt(ks,i,j-1,i_xyz) * random_mz(ks,i,j-1) ) / ( fdy(j) + fdy(j-1) ) * mapf(i,j,2,i_xy) &
1123  + ( j23g(ks+1,i,j,i_xyz) * random_mz(ks+1,i,j) - j23g(ks,i,j,i_xyz) * random_mz(ks,i,j) ) * rfdz(ks) * 0.5_rp &
1124  - j33g * ( random_my(ks+1,i,j) - random_my(ks,i,j) ) * rfdz(ks) * 0.5_rp &
1125  ) / gsqrt(ks,i,j,i_xyz)
1126  acc(ke,i,j) = ( ( gsqrt(ke,i,j+1,i_xyz) * random_mz(ke,i,j+1) - gsqrt(ke,i,j-1,i_xyz) * random_mz(ke,i,j-1) ) / ( fdy(j) + fdy(j-1) ) * mapf(i,j,2,i_xy) &
1127  + ( j23g(ke,i,j,i_xyz) * random_mz(ke,i,j) - j23g(ke-1,i,j,i_xyz) * random_mz(ke-1,i,j) ) * rfdz(ke-1) * 0.5_rp &
1128  - j33g * ( random_my(ke,i,j) - random_my(ke-1,i,j) ) * rfdz(ke-1) * 0.5_rp &
1129  ) / gsqrt(ke,i,j,i_xyz)
1130  end do
1131  end do
1132 
1133  !$omp parallel do
1134  do j = jjs, jje
1135  do i = iis, iie
1136  do k = ks, ke
1137  momx_t(k,i,j) = ( acc(k,i+1,j) + acc(k,i,j) ) * 0.5_rp
1138  end do
1139  end do
1140  end do
1141 
1142  ! MOMY : dfx/dz - dfz/dx
1143  !$omp parallel do
1144  do j = jjs, jje+1
1145  do i = iis, iie
1146  do k = ks+1, ke-1
1147  acc(k,i,j) = ( j33g * ( random_mx(k+1,i,j) - random_mx(k-1,i,j) ) / ( fdz(k) + fdz(k-1) ) &
1148  - ( gsqrt(k,i+1,j,i_xyz) * random_mz(k,i+1,j) - gsqrt(k,i-1,j,i_xyz) * random_mz(k,i-1,j) ) / ( fdx(i) + fdx(i-1) ) * mapf(i,j,1,i_xy) &
1149  + ( j13g(k+1,i,j,i_xyz) * random_mz(k+1,i,j) - j13g(k-1,i,j,i_xyz) * random_mz(k-1,i,j) ) / ( fdz(k) + fdz(k-1) ) &
1150  ) / gsqrt(k,i,j,i_xyz)
1151  end do
1152  acc(ks,i,j) = ( j33g * ( random_mx(ks+1,i,j) - random_mx(ks,i,j) ) * rfdz(ks) * 0.5_rp &
1153  - ( gsqrt(ks,i+1,j,i_xyz) * random_mz(ks,i+1,j) - gsqrt(ks,i-1,j,i_xyz) * random_mz(ks,i-1,j) ) / ( fdx(i) + fdx(i-1) ) * mapf(i,j,1,i_xy) &
1154  + ( j13g(ks+1,i,j,i_xyz) * random_mz(ks+1,i,j) - j13g(ks,i,j,i_xyz) * random_mz(ks,i,j) ) * rfdz(ks) * 0.5_rp &
1155  ) / gsqrt(ks,i,j,i_xyz)
1156  acc(ke,i,j) = ( j33g * ( random_mx(ke,i,j) - random_mx(ke-1,i,j) ) * rfdz(ke-1) * 0.5_rp &
1157  - ( gsqrt(ke,i+1,j,i_xyz) * random_mz(ke,i+1,j) - gsqrt(ke,i-1,j,i_xyz) * random_mz(ke,i-1,j) ) / ( fdx(i) + fdx(i-1) ) * mapf(i,j,1,i_xy) &
1158  + ( j13g(ke,i,j,i_xyz) * random_mz(ke,i,j) - j13g(ke-1,i,j,i_xyz) * random_mz(ke-1,i,j) ) * rfdz(ke-1) * 0.5_rp &
1159  ) / gsqrt(ke,i,j,i_xyz)
1160  end do
1161  end do
1162 
1163  !$omp parallel do
1164  do j = jjs, jje
1165  do i = iis, iie
1166  do k = ks, ke
1167  momy_t(k,i,j) = ( acc(k,i,j+1) + acc(k,i,j) ) * 0.5_rp
1168  end do
1169  end do
1170  end do
1171 
1172  else
1173 
1174  !$omp parallel do
1175 !OCL XFILL
1176  do j = jjs, jje
1177  do i = iis, iie
1178  do k = ks, ke
1179  momz_t(k,i,j) = 0.0_rp
1180  end do
1181  end do
1182  end do
1183  !$omp parallel do
1184 !OCL XFILL
1185  do j = jjs, jje
1186  do i = iis, iie
1187  do k = ks, ke
1188  momx_t(k,i,j) = 0.0_rp
1189  end do
1190  end do
1191  end do
1192  !$omp parallel do
1193 !OCL XFILL
1194  do j = jjs, jje
1195  do i = iis, iie
1196  do k = ks, ke
1197  momy_t(k,i,j) = 0.0_rp
1198  end do
1199  end do
1200  end do
1201 
1202  end if
1203 
1204  !##### Thermodynamic Equation #####
1205 
1206  if ( atmos_phy_tb_smg_implicit ) then
1207 
1208  !$omp parallel do &
1209  !$omp private (ap,d)
1210  do j = jjs, jje
1211  do i = iis, iie
1212 
1213  ap = - dt * 0.25_rp * ( dens(ks,i,j)+dens(ks+1,i,j) ) &
1214  * ( kh(ks,i,j)+kh(ks+1,i,j) ) &
1215  * rfdz(ks) / gsqrt(ks,i,j,i_xyw)
1216  a(ks,i,j) = ap * rcdz(ks) / gsqrt(ks,i,j,i_xyz)
1217  c(ks,i,j) = 0.0_rp
1218  b(ks,i,j) = - a(ks,i,j) + dens(ks,i,j)
1219  do k = ks+1, ke-1
1220  c(k,i,j) = ap * rcdz(k) / gsqrt(k,i,j,i_xyz)
1221  ap = - dt * 0.25_rp * ( dens(k,i,j)+dens(k+1,i,j) ) &
1222  * ( kh(k,i,j)+kh(k+1,i,j) ) &
1223  * rfdz(k) / gsqrt(k,i,j,i_xyw)
1224  a(k,i,j) = ap * rcdz(k) / gsqrt(k,i,j,i_xyz)
1225  b(k,i,j) = - a(k,i,j) - c(k,i,j) + dens(k,i,j)
1226  end do
1227  a(ke,i,j) = 0.0_rp
1228  c(ke,i,j) = ap * rcdz(ke) / gsqrt(ke,i,j,i_xyz)
1229  b(ke,i,j) = - c(ke,i,j) + dens(ke,i,j)
1230 
1231  end do
1232  end do
1233 
1234  end if
1235 
1236  call calc_flux_phi( &
1237  qflx_sgs_rhot, &
1238  dens, pott, kh, 1.0_rp, &
1239  gsqrt, j13g, j23g, j33g, mapf, &
1240  atmos_phy_tb_smg_horizontal, &
1241  atmos_phy_tb_smg_implicit, &
1242  a, b, c, dt, &
1243  iis, iie, jjs, jje )
1244 
1245  if ( atmos_phy_tb_smg_backscatter ) then
1246 
1247  !$omp parallel do
1248  do j = jjs-1, jje+1
1249  do i = iis-1, iie+1
1250 
1251  do k = ks+1, ke-1
1252  dd(k,i,j) = sqrt( ( ( pott(k+1,i,j) - pott(k-1,i,j) ) * j33g / ( fdz(k) + fdz(k-1) ) )**2 &
1253  + ( ( ( gsqrt(k,i+1,j,i_xyz)*pott(k,i+1,j) - gsqrt(k,i-1,j,i_xyz)*pott(k,i-1,j) ) / ( fdx(i) + fdx(i-1) ) &
1254  + ( j13g(k+1,i,j,i_xyz)*pott(k+1,i,j) - j13g(k-1,i,j,i_xyz)*pott(k-1,i,j) ) / ( fdz(k) + fdz(k-1) ) ) * mapf(i,j,1,i_xy) )**2 &
1255  + ( ( ( gsqrt(k,i,j+1,i_xyz)*pott(k,i,j+1) - gsqrt(k,i,j-1,i_xyz)*pott(k,i,j-1) ) / ( fdy(j) + fdy(j-1) ) &
1256  + ( j23g(k+1,i,j,i_xyz)*pott(k+1,i,j) - j23g(k-1,i,j,i_xyz)*pott(k-1,i,j) ) / ( fdz(k) + fdz(k-1) ) ) * mapf(i,j,2,i_xy) )**2 &
1257  ) / gsqrt(k,i,j,i_xyz)
1258  end do
1259  dd(ks,i,j) = sqrt( ( ( pott(ks+1,i,j) - pott(ks,i,j) ) * j33g * rfdz(ks) )**2 &
1260  + ( ( ( gsqrt(ks,i+1,j,i_xyz)*pott(ks,i+1,j) - gsqrt(ks,i-1,j,i_xyz)*pott(ks,i-1,j) ) / ( fdx(i) + fdx(i-1) ) &
1261  + ( j13g(ks+1,i,j,i_xyz)*pott(ks+1,i,j) - j13g(ks,i,j,i_xyz)*pott(ks,i,j) ) * rfdz(ks) ) * mapf(i,j,1,i_xy) )**2 &
1262  + ( ( ( gsqrt(ks,i,j+1,i_xyz)*pott(ks,i,j+1) - gsqrt(ks,i,j-1,i_xyz)*pott(ks,i,j-1) ) / ( fdy(j) + fdy(j-1) ) &
1263  + ( j23g(ks+1,i,j,i_xyz)*pott(ks+1,i,j) - j23g(ks,i,j,i_xyz)*pott(ks,i,j) ) * rfdz(ks) ) * mapf(i,j,2,i_xy) )**2 &
1264  ) / gsqrt(ks,i,j,i_xyz)
1265  dd(ke,i,j) = sqrt( ( ( pott(ke,i,j) - pott(ke-1,i,j) ) * j33g * rfdz(ke-1) )**2 &
1266  + ( ( ( gsqrt(ke,i+1,j,i_xyz)*pott(ke,i+1,j) - gsqrt(ke,i-1,j,i_xyz)*pott(ke,i-1,j) ) / ( fdx(i) + fdx(i-1) ) &
1267  + ( j13g(ke,i,j,i_xyz)*pott(ke,i,j) - j13g(ke-1,i,j,i_xyz)*pott(ke-1,i,j) ) * rfdz(ke-1) ) * mapf(i,j,1,i_xy) )**2 &
1268  + ( ( ( gsqrt(ke,i,j+1,i_xyz)*pott(ke,i,j+1) - gsqrt(ke,i,j-1,i_xyz)*pott(ke,i,j-1) ) / ( fdy(j) + fdy(j-1) ) &
1269  + ( j23g(ke,i,j,i_xyz)*pott(ke,i,j) - j23g(ke-1,i,j,i_xyz)*pott(ke-1,i,j) ) * rfdz(ke-1) ) * mapf(i,j,2,i_xy) )**2 &
1270  ) / gsqrt(ke,i,j,i_xyz)
1271 
1272  end do
1273  end do
1274 
1275  !$omp parallel do
1276  do j = jjs, jje
1277  do i = iis, iie
1278  do k = ks+1, ke-1
1279  rhot_t(k,i,j) = ( j33g * ( random_qz(k+1,i,j) * dd(k+1,i,j) - random_qz(k-1,i,j) * dd(k-1,i,j) ) / ( fdz(k) + fdz(k-1) ) &
1280  + ( gsqrt(k,i+1,j,i_xyz) * random_qx(k,i+1,j) * dd(k,i+1,j) - gsqrt(k,i-1,j,i_xyz) * random_qx(k,i-1,j) * dd(k,i-1,j) ) / ( fdx(i) + fdx(i-1) ) * mapf(i,j,1,i_xy) &
1281  + ( j13g(k+1,i,j,i_xyz) * random_qx(k+1,i,j) * dd(k+1,i,j) - j13g(k-1,i,j,i_xyz) * random_qx(k-1,i,j) * dd(k-1,i,j) ) / ( fdz(k) + fdz(k-1) ) &
1282  + ( gsqrt(k,i,j+1,i_xyz) * random_qy(k,i,j+1) * dd(k,i,j+1) - gsqrt(k,i,j-1,i_xyz) * random_qy(k,i,j-1) * dd(k,i,j-1) ) / ( fdy(j) + fdy(j-1) ) * mapf(i,j,2,i_xy) &
1283  + ( j23g(k+1,i,j,i_xyz) * random_qy(k+1,i,j) * dd(k+1,i,j) - j23g(k-1,i,j,i_xyz) * random_qy(k-1,i,j) * dd(k-1,i,j) ) / ( fdz(k) + fdz(k-1) ) &
1284  ) / gsqrt(k,i,j,i_xyz)
1285  end do
1286  rhot_t(ks,i,j) = ( j33g * ( random_qz(ks+1,i,j) * dd(ks+1,i,j) - random_qz(ks,i,j) * dd(ks,i,j) ) * rfdz(ks) &
1287  + ( gsqrt(ks,i+1,j,i_xyz) * random_qx(ks,i+1,j) * dd(ks,i+1,j) - gsqrt(ks,i-1,j,i_xyz) * random_qx(ks,i-1,j) * dd(ks,i-1,j) ) / ( fdx(i) + fdx(i-1) ) * mapf(i,j,1,i_xy) &
1288  + ( j13g(ks+1,i,j,i_xyz) * random_qx(ks+1,i,j) * dd(ks+1,i,j) - j13g(ks,i,j,i_xyz) * random_qx(ks,i,j) * dd(ks,i,j) ) * rfdz(ks) * 0.5_rp &
1289  + ( gsqrt(ks,i,j+1,i_xyz) * random_qy(ks,i,j+1) * dd(ks,i,j+1) - gsqrt(ks,i,j-1,i_xyz) * random_qy(ks,i,j-1) * dd(ks,i,j-1) ) / ( fdy(j) + fdy(j-1) ) * mapf(i,j,2,i_xy) &
1290  + ( j23g(ks+1,i,j,i_xyz) * random_qy(ks+1,i,j) * dd(ks+1,i,j) - j23g(ks,i,j,i_xyz) * random_qy(ks,i,j) * dd(ks,i,j) ) * rfdz(ks) * 0.5_rp &
1291  ) / gsqrt(ks,i,j,i_xyz)
1292  rhot_t(ke,i,j) = ( j33g * ( random_qz(ke,i,j) * dd(ke,i,j) - random_qz(ke-1,i,j) * dd(ke-1,i,j) ) * rfdz(ke-1) * 0.5_rp &
1293  + ( gsqrt(ke,i+1,j,i_xyz) * random_qx(ke,i+1,j) * dd(ke,i+1,j) - gsqrt(ke,i-1,j,i_xyz) * random_qx(ke,i-1,j) * dd(ke,i-1,j) ) / ( fdx(i) + fdx(i-1) ) * mapf(i,j,1,i_xy) &
1294  + ( j13g(ke,i,j,i_xyz) * random_qx(ke,i,j) * dd(ke,i,j) - j13g(ke-1,i,j,i_xyz) * random_qx(ke-1,i,j) * dd(ke-1,i,j) ) * rfdz(ke-1) * 0.5_rp &
1295  + ( gsqrt(ke,i,j+1,i_xyz) * random_qy(ke,i,j+1) * dd(ke,i,j+1) - gsqrt(ke,i,j-1,i_xyz) * random_qy(ke,i,j-1) * dd(ke,i,j-1) ) / ( fdy(j) + fdy(j-1) ) * mapf(i,j,2,i_xy) &
1296  + ( j23g(ke,i,j,i_xyz) * random_qy(ke,i,j) * dd(ke,i,j) - j23g(ke-1,i,j,i_xyz) * random_qy(ke-1,i,j) * dd(ke-1,i,j) ) * rfdz(ke-1) * 0.5_rp &
1297  ) / gsqrt(ke,i,j,i_xyz)
1298  end do
1299  end do
1300 
1301  else
1302 
1303  !$omp parallel do
1304 !OCL XFILL
1305  do j = jjs, jje
1306  do i = iis, iie
1307  do k = ks, ke
1308  rhot_t(k,i,j) = 0.0_rp
1309  end do
1310  end do
1311  end do
1312 
1313  end if
1314 
1315  enddo
1316  enddo
1317 
1318  !##### Tracers #####
1319  do iq = 1, qa
1320 
1321  if ( .not. tracer_advc(iq) ) cycle
1322 
1323  do jjs = js, je, jblock
1324  jje = jjs+jblock-1
1325  do iis = is, ie, iblock
1326  iie = iis+iblock-1
1327 
1328  call calc_flux_phi( &
1329  qflx_sgs_rhoq(:,:,:,:,iq), &
1330  dens, qtrc(:,:,:,iq), kh, 1.0_rp, &
1331  gsqrt, j13g, j23g, j33g, mapf, &
1332  atmos_phy_tb_smg_horizontal, &
1333  atmos_phy_tb_smg_implicit, &
1334  a, b, c, dt, &
1335  iis, iie, jjs, jje )
1336 
1337 
1338  if ( atmos_phy_tb_smg_backscatter .and. iq == i_qv ) then
1339 
1340  !$omp parallel do
1341  do j = jjs-1, jje+1
1342  do i = iis-1, iie+1
1343 
1344  do k = ks+1, ke-1
1345  dd(k,i,j) = sqrt( ( ( qtrc(k+1,i,j,iq) - qtrc(k-1,i,j,iq) ) * j33g / ( fdz(k) + fdz(k-1) ) )**2 &
1346  + ( ( ( gsqrt(k,i+1,j,i_xyz)*qtrc(k,i+1,j,iq) - gsqrt(k,i-1,j,i_xyz)*qtrc(k,i-1,j,iq) ) / ( fdx(i) + fdx(i-1) ) &
1347  + ( j13g(k+1,i,j,i_xyz)*qtrc(k+1,i,j,iq) - j13g(k-1,i,j,i_xyz)*qtrc(k-1,i,j,iq) ) / ( fdz(k) + fdz(k-1) ) ) * mapf(i,j,1,i_xy) )**2 &
1348  + ( ( ( gsqrt(k,i,j+1,i_xyz)*qtrc(k,i,j+1,iq) - gsqrt(k,i,j-1,i_xyz)*qtrc(k,i,j-1,iq) ) / ( fdy(j) + fdy(j-1) ) &
1349  + ( j23g(k+1,i,j,i_xyz)*qtrc(k+1,i,j,iq) - j23g(k-1,i,j,i_xyz)*qtrc(k-1,i,j,iq) ) / ( fdz(k) + fdz(k-1) ) ) * mapf(i,j,2,i_xy) )**2 &
1350  ) / gsqrt(k,i,j,i_xyz)
1351  end do
1352  dd(ks,i,j) = sqrt( ( ( qtrc(ks+1,i,j,iq) - qtrc(ks,i,j,iq) ) * j33g * rfdz(ks) )**2 &
1353  + ( ( ( gsqrt(ks,i+1,j,i_xyz)*qtrc(ks,i+1,j,iq) - gsqrt(ks,i-1,j,i_xyz)*qtrc(ks,i-1,j,iq) ) / ( fdx(i) + fdx(i-1) ) &
1354  + ( j13g(ks+1,i,j,i_xyz)*qtrc(ks+1,i,j,iq) - j13g(ks,i,j,i_xyz)*qtrc(ks,i,j,iq) ) * rfdz(ks) ) * mapf(i,j,1,i_xy) )**2 &
1355  + ( ( ( gsqrt(ks,i,j+1,i_xyz)*qtrc(ks,i,j+1,iq) - gsqrt(ks,i,j-1,i_xyz)*qtrc(ks,i,j-1,iq) ) / ( fdy(j) + fdy(j-1) ) &
1356  + ( j23g(ks+1,i,j,i_xyz)*qtrc(ks+1,i,j,iq) - j23g(ks,i,j,i_xyz)*qtrc(ks,i,j,iq) ) * rfdz(ks) ) * mapf(i,j,2,i_xy) )**2 &
1357  ) / gsqrt(ks,i,j,i_xyz)
1358  dd(ke,i,j) = sqrt( ( ( qtrc(ke,i,j,iq) - qtrc(ke-1,i,j,iq) ) * j33g * rfdz(ke-1) )**2 &
1359  + ( ( ( gsqrt(ke,i+1,j,i_xyz)*qtrc(ke,i+1,j,iq) - gsqrt(ke,i-1,j,i_xyz)*qtrc(ke,i-1,j,iq) ) / ( fdx(i) + fdx(i-1) ) &
1360  + ( j13g(ke,i,j,i_xyz)*qtrc(ke,i,j,iq) - j13g(ke-1,i,j,i_xyz)*qtrc(ke-1,i,j,iq) ) * rfdz(ke-1) ) * mapf(i,j,1,i_xy) )**2 &
1361  + ( ( ( gsqrt(ke,i,j+1,i_xyz)*qtrc(ke,i,j+1,iq) - gsqrt(ke,i,j-1,i_xyz)*qtrc(ke,i,j-1,iq) ) / ( fdy(j) + fdy(j-1) ) &
1362  + ( j23g(ke,i,j,i_xyz)*qtrc(ke,i,j,iq) - j23g(ke-1,i,j,i_xyz)*qtrc(ke-1,i,j,iq) ) * rfdz(ke-1) ) * mapf(i,j,2,i_xy) )**2 &
1363  ) / gsqrt(ke,i,j,i_xyz)
1364 
1365  end do
1366  end do
1367 
1368  !$omp parallel do
1369  do j = jjs, jje
1370  do i = iis, iie
1371  do k = ks+1, ke-1
1372  rhoq_t(k,i,j,iq) = ( j33g * ( random_qz(k+1,i,j) * dd(k+1,i,j) - random_qz(k-1,i,j) * dd(k-1,i,j) ) / ( fdz(k) + fdz(k-1) ) &
1373  + ( gsqrt(k,i+1,j,i_xyz) * random_qx(k,i+1,j) * dd(k,i+1,j) - gsqrt(k,i-1,j,i_xyz) * random_qx(k,i-1,j) * dd(k,i-1,j) ) / ( fdx(i) + fdx(i-1) ) * mapf(i,j,1,i_xy) &
1374  + ( j13g(k+1,i,j,i_xyz) * random_qx(k+1,i,j) * dd(k+1,i,j) - j13g(k-1,i,j,i_xyz) * random_qx(k-1,i,j) * dd(k-1,i,j) ) / ( fdz(k) + fdz(k-1) ) &
1375  + ( gsqrt(k,i,j+1,i_xyz) * random_qy(k,i,j+1) * dd(k,i,j+1) - gsqrt(k,i,j-1,i_xyz) * random_qy(k,i,j-1) * dd(k,i,j-1) ) / ( fdy(j) + fdy(j-1) ) * mapf(i,j,2,i_xy) &
1376  + ( j23g(k+1,i,j,i_xyz) * random_qy(k+1,i,j) * dd(k+1,i,j) - j23g(k-1,i,j,i_xyz) * random_qy(k-1,i,j) * dd(k-1,i,j) ) / ( fdz(k) + fdz(k-1) ) &
1377  ) / gsqrt(k,i,j,i_xyz)
1378  end do
1379  rhoq_t(ks,i,j,iq) = ( j33g * ( random_qz(ks+1,i,j) * dd(ks+1,i,j) - random_qz(ks,i,j) * dd(ks,i,j) ) * rfdz(ks) * 0.5_rp &
1380  + ( gsqrt(ks,i+1,j,i_xyz) * random_qx(ks,i+1,j) * dd(ks,i+1,j) - gsqrt(ks,i-1,j,i_xyz) * random_qx(ks,i-1,j) * dd(ks,i-1,j) ) / ( fdx(i) + fdx(i-1) ) * mapf(i,j,1,i_xy) &
1381  + ( j13g(ks+1,i,j,i_xyz) * random_qx(ks+1,i,j) * dd(ks+1,i,j) - j13g(ks,i,j,i_xyz) * random_qx(ks,i,j) * dd(ks,i,j) ) * rfdz(ks) * 0.5_rp &
1382  + ( gsqrt(ks,i,j+1,i_xyz) * random_qy(ks,i,j+1) * dd(ks,i,j+1) - gsqrt(ks,i,j-1,i_xyz) * random_qy(ks,i,j-1) * dd(ks,i,j-1) ) / ( fdy(j) + fdy(j-1) ) * mapf(i,j,2,i_xy) &
1383  + ( j23g(ks+1,i,j,i_xyz) * random_qy(ks+1,i,j) * dd(ks+1,i,j) - j23g(ks,i,j,i_xyz) * random_qy(ks,i,j) * dd(ks,i,j) ) * rfdz(ks) * 0.5_rp &
1384  ) / gsqrt(ks,i,j,i_xyz)
1385  rhoq_t(ke,i,j,iq) = ( j33g * ( random_qz(ke,i,j) * dd(ke,i,j) - random_qz(ke-1,i,j) * dd(ke-1,i,j) ) * rfdz(ke-1) * 0.5_rp &
1386  + ( gsqrt(ke,i+1,j,i_xyz) * random_qx(ke,i+1,j) * dd(ke,i+1,j) - gsqrt(ke,i-1,j,i_xyz) * random_qx(ke,i-1,j) * dd(ke,i-1,j) ) / ( fdx(i) + fdx(i-1) ) * mapf(i,j,1,i_xy) &
1387  + ( j13g(ke,i,j,i_xyz) * random_qx(ke,i,j) * dd(ke,i,j) - j13g(ke-1,i,j,i_xyz) * random_qx(ke-1,i,j) * dd(ke-1,i,j) ) * rfdz(ke-1) * 0.5_rp &
1388  + ( gsqrt(ke,i,j+1,i_xyz) * random_qy(ke,i,j+1) * dd(ke,i,j+1) - gsqrt(ke,i,j-1,i_xyz) * random_qy(ke,i,j-1) * dd(ke,i,j-1) ) / ( fdy(j) + fdy(j-1) ) * mapf(i,j,2,i_xy) &
1389  + ( j23g(ke,i,j,i_xyz) * random_qy(ke,i,j) * dd(ke,i,j) - j23g(ke-1,i,j,i_xyz) * random_qy(ke-1,i,j) * dd(ke-1,i,j) ) * rfdz(ke-1) * 0.5_rp &
1390  ) / gsqrt(ke,i,j,i_xyz)
1391  end do
1392  end do
1393 
1394  else
1395 
1396  !$omp parallel do
1397 !OCL XFILL
1398  do j = jjs, jje
1399  do i = iis, iie
1400  do k = ks, ke
1401  rhoq_t(k,i,j,iq) = 0.0_rp
1402  end do
1403  end do
1404  end do
1405 
1406  end if
1407 
1408  enddo
1409  enddo
1410 
1411  enddo ! scalar quantities loop
1412 #ifdef DEBUG
1413  iq = iundef
1414 #endif
1415 
1416  call file_history_in( tke(:,:,:), 'TKE_SMG', 'turbulent kinetic energy (Smagorinsky)', 'm2/s2', fill_halo=.true. )
1417 
1418 
1419  return
1420  end subroutine atmos_phy_tb_smg
1421 
1422 
1423  function mixlen(dz, dx, dy, filter_fact)
1424  implicit none
1425  real(RP), intent(in) :: dz
1426  real(RP), intent(in) :: dx
1427  real(RP), intent(in) :: dy
1428  real(RP), intent(in) :: filter_fact
1429  real(RP) :: mixlen ! (out)
1430 
1431  mixlen = fact(dz, dx, dy) * filter_fact * ( dz * dx * dy )**oneoverthree ! Scotti et al. (1993)
1432 
1433  return
1434  end function mixlen
1435 
1436  function fact(dz, dx, dy)
1437  real(RP), intent(in) :: dz
1438  real(RP), intent(in) :: dx
1439  real(RP), intent(in) :: dy
1440  real(RP) :: fact ! (out)
1441 
1442  real(RP), parameter :: oot = -1.0_rp/3.0_rp
1443  real(RP), parameter :: fot = 5.0_rp/3.0_rp
1444  real(RP), parameter :: eot = 11.0_rp/3.0_rp
1445  real(RP), parameter :: tof = -3.0_rp/4.0_rp
1446  real(RP) :: a1, a2, b1, b2, dmax
1447 
1448 
1449  dmax = max(dz, dx, dy)
1450  if ( dz == dmax ) then
1451  a1 = dx / dmax
1452  a2 = dy / dmax
1453  else if ( dx == dmax ) then
1454  a1 = dz / dmax
1455  a2 = dy / dmax
1456  else ! dy == dmax
1457  a1 = dz / dmax
1458  a2 = dx / dmax
1459  end if
1460  b1 = atan( a1/a2 )
1461  b2 = atan( a2/a1 )
1462 
1463  fact = 1.736_rp * (a1*a2)**oot &
1464  * ( 4.0_rp*p1(b1)*a1**oot + 0.222_rp*p2(b1)*a1**fot + 0.077*p3(b1)*a1**eot - 3.0_rp*b1 &
1465  + 4.0_rp*p1(b2)*a2**oot + 0.222_rp*p2(b2)*a2**fot + 0.077*p3(b2)*a2**eot - 3.0_rp*b2 &
1466  )**tof
1467  return
1468  end function fact
1469  function p1(z)
1470  real(RP), intent(in) :: z
1471  real(RP) :: p1 ! (out)
1472 
1473  p1 = 2.5_rp * p2(z) - 1.5_rp * sin(z) * cos(z)**twooverthree
1474  return
1475  end function p1
1476  function p2(z)
1477  real(RP), intent(in) :: z
1478  real(RP) :: p2 ! (out)
1479 
1480  p2 = 0.986_rp * z + 0.073_rp * z**2 - 0.418_rp * z**3 + 0.120_rp * z**4
1481  return
1482  end function p2
1483  function p3(z)
1484  real(RP), intent(in) :: z
1485  real(RP) :: p3 ! (out)
1486 
1487  p3 = 0.976_rp * z + 0.188_rp * z**2 - 1.169_rp * z**3 + 0.755_rp * z**4 - 0.151_rp * z**5
1488  return
1489  end function p3
1490 
1491 end module scale_atmos_phy_tb_smg
scale_debug
module DEBUG
Definition: scale_debug.F90:11
scale_atmos_grid_cartesc_index::i_xyz
integer, public i_xyz
Definition: scale_atmos_grid_cartesC_index.F90:90
scale_atmos_phy_tb_smg::mixlen
real(rp) function mixlen(dz, dx, dy, filter_fact)
Definition: scale_atmos_phy_tb_smg.F90:1424
scale_atmos_phy_tb_common::atmos_phy_tb_calc_tend_momz
subroutine, public atmos_phy_tb_calc_tend_momz(MOMZ_t_TB, QFLX_MOMZ, GSQRT, J13G, J23G, J33G, MAPF, IIS, IIE, JJS, JJE)
Definition: scale_atmos_phy_tb_common.F90:1534
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:1492
scale_atmos_phy_tb_smg::atmos_phy_tb_smg_setup
subroutine, public atmos_phy_tb_smg_setup(FZ, CZ, CDX, CDY, MAPF)
Setup.
Definition: scale_atmos_phy_tb_smg.F90:99
scale_atmos_phy_tb_common::atmos_phy_tb_calc_tend_momy
subroutine, public atmos_phy_tb_calc_tend_momy(MOMY_t_TB, QFLX_MOMY, GSQRT, J13G, J23G, J33G, MAPF, IIS, IIE, JJS, JJE)
Definition: scale_atmos_phy_tb_common.F90:1653
scale_atmos_grid_cartesc_index::ia
integer, public ia
of whole cells: x, local, with HALO
Definition: scale_atmos_grid_cartesC_index.F90:48
scale_atmos_grid_cartesc_index::i_xy
integer, public i_xy
Definition: scale_atmos_grid_cartesC_index.F90:98
scale_atmos_grid_cartesc_index::i_xvw
integer, public i_xvw
Definition: scale_atmos_grid_cartesC_index.F90:93
scale_atmos_grid_cartesc_index::i_xvz
integer, public i_xvz
Definition: scale_atmos_grid_cartesC_index.F90:95
scale_atmos_grid_cartesc_index::iblock
integer, public iblock
block size for cache blocking: x
Definition: scale_atmos_grid_cartesC_index.F90:40
scale_tracer::qa
integer, public qa
Definition: scale_tracer.F90:37
scale_tracer::tracer_advc
logical, dimension(qa_max), public tracer_advc
Definition: scale_tracer.F90:45
scale_atmos_grid_cartesc_index::ja
integer, public ja
of whole cells: y, local, with HALO
Definition: scale_atmos_grid_cartesC_index.F90:49
scale_io::io_fid_conf
integer, public io_fid_conf
Config file ID.
Definition: scale_io.F90:55
scale_debug::check
subroutine, public check(current_line, v)
Undefined value checker.
Definition: scale_debug.F90:56
scale_const::const_karman
real(rp), parameter, public const_karman
von Karman constant
Definition: scale_const.F90:50
scale_const::const_undef
real(rp), public const_undef
Definition: scale_const.F90:41
scale_atmos_grid_cartesc_index::is
integer, public is
start point of inner domain: x, local
Definition: scale_atmos_grid_cartesC_index.F90:53
scale_atmos_grid_cartesc_index::ie
integer, public ie
end point of inner domain: x, local
Definition: scale_atmos_grid_cartesC_index.F90:54
scale_tracer
module TRACER
Definition: scale_tracer.F90:16
scale_atmos_phy_tb_smg
module ATMOSPHERE / Physics Turbulence
Definition: scale_atmos_phy_tb_smg.F90:21
scale_atmos_phy_tb_common::atmos_phy_tb_calc_tend_momx
subroutine, public atmos_phy_tb_calc_tend_momx(MOMX_t_TB, QFLX_MOMX, GSQRT, J13G, J23G, J33G, MAPF, IIS, IIE, JJS, JJE)
Definition: scale_atmos_phy_tb_common.F90:1594
scale_atmos_grid_cartesc_index::ydir
integer, parameter, public ydir
Definition: scale_atmos_grid_cartesC_index.F90:34
scale_atmos_hydrometeor
module atmosphere / hydrometeor
Definition: scale_atmos_hydrometeor.F90:12
scale_atmos_grid_cartesc_index
module atmosphere / grid / cartesC index
Definition: scale_atmos_grid_cartesC_index.F90:12
scale_atmos_grid_cartesc_index::ke
integer, public ke
end point of inner domain: z, local
Definition: scale_atmos_grid_cartesC_index.F90:52
scale_prc
module PROCESS
Definition: scale_prc.F90:11
scale_atmos_grid_cartesc_index::je
integer, public je
end point of inner domain: y, local
Definition: scale_atmos_grid_cartesC_index.F90:56
scale_atmos_phy_tb_smg::atmos_phy_tb_smg
subroutine, public atmos_phy_tb_smg(qflx_sgs_momz, qflx_sgs_momx, qflx_sgs_momy, qflx_sgs_rhot, qflx_sgs_rhoq, MOMZ_t, MOMX_t, MOMY_t, RHOT_t, RHOQ_t, Nu, Ri, Pr, MOMZ, MOMX, MOMY, POTT, DENS, QTRC, N2, FZ, FDZ, RCDZ, RFDZ, CDX, FDX, CDY, FDY, GSQRT, J13G, J23G, J33G, MAPF, dt)
Definition: scale_atmos_phy_tb_smg.F90:237
scale_const::const_grav
real(rp), public const_grav
standard acceleration of gravity [m/s2]
Definition: scale_const.F90:46
scale_atmos_phy_tb_common::atmos_phy_tb_calc_strain_tensor
subroutine, public atmos_phy_tb_calc_strain_tensor(S33_C, S11_C, S22_C, S31_C, S12_C, S23_C, S12_Z, S23_X, S31_Y, S2, DENS, MOMZ, MOMX, MOMY, GSQRT, J13G, J23G, J33G, MAPF)
Definition: scale_atmos_phy_tb_common.F90:67
scale_const::const_undef2
integer, parameter, public const_undef2
undefined value (INT2)
Definition: scale_const.F90:38
scale_atmos_grid_cartesc_index::ks
integer, public ks
start point of inner domain: z, local
Definition: scale_atmos_grid_cartesC_index.F90:51
scale_atmos_grid_cartesc_index::jblock
integer, public jblock
block size for cache blocking: y
Definition: scale_atmos_grid_cartesC_index.F90:41
scale_prc::prc_abort
subroutine, public prc_abort
Abort Process.
Definition: scale_prc.F90:338
scale_const
module CONSTANT
Definition: scale_const.F90:11
scale_atmos_grid_cartesc_index::js
integer, public js
start point of inner domain: y, local
Definition: scale_atmos_grid_cartesC_index.F90:55
scale_atmos_grid_cartesc_index::xdir
integer, parameter, public xdir
Definition: scale_atmos_grid_cartesC_index.F90:33
scale_atmos_phy_tb_common
module ATMOSPHERE / Physics Turbulence
Definition: scale_atmos_phy_tb_common.F90:12
scale_atmos_grid_cartesc_index::i_uyw
integer, public i_uyw
Definition: scale_atmos_grid_cartesC_index.F90:92
scale_atmos_phy_tb_common::atmos_phy_tb_calc_flux_phi
subroutine, public atmos_phy_tb_calc_flux_phi(qflx_phi, DENS, PHI, Kh, FACT, GSQRT, J13G, J23G, J33G, MAPF, horizontal, implicit, a, b, c, dt, IIS, IIE, JJS, JJE)
Definition: scale_atmos_phy_tb_common.F90:1234
scale_atmos_grid_cartesc_index::i_uyz
integer, public i_uyz
Definition: scale_atmos_grid_cartesC_index.F90:94
scale_prof
module profiler
Definition: scale_prof.F90:11
scale_const::const_eps
real(rp), public const_eps
small number
Definition: scale_const.F90:33
scale_atmos_hydrometeor::i_qv
integer, public i_qv
Definition: scale_atmos_hydrometeor.F90:77
scale_random::random_normal
subroutine, public random_normal(var)
Get normal random number.
Definition: scale_random.F90:153
scale_precision
module PRECISION
Definition: scale_precision.F90:14
scale_atmos_grid_cartesc_index::ka
integer, public ka
of whole cells: z, local, with HALO
Definition: scale_atmos_grid_cartesC_index.F90:47
scale_io
module STDIO
Definition: scale_io.F90:10
scale_random
module RANDOM
Definition: scale_random.F90:14
scale_atmos_grid_cartesc_index::zdir
integer, parameter, public zdir
Definition: scale_atmos_grid_cartesC_index.F90:32
scale_file_history
module file_history
Definition: scale_file_history.F90:15
scale_atmos_grid_cartesc_index::i_xyw
integer, public i_xyw
Definition: scale_atmos_grid_cartesC_index.F90:91
Generated by   doxygen 1.8.13