SCALE-RM
scale_atmos_phy_tb_d1980.F90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
19 !-------------------------------------------------------------------------------
21  !-----------------------------------------------------------------------------
22  !
23  !++ used modules
24  !
25  use scale_precision
26  use scale_stdio
27  use scale_prof
29  use scale_tracer
30 
31 #ifdef DEBUG
32  use scale_debug, only: &
33  check
34  use scale_const, only: &
35  undef => const_undef, &
36  iundef => const_undef2
37 #endif
38  !-----------------------------------------------------------------------------
39  implicit none
40  private
41  !-----------------------------------------------------------------------------
42  !
43  !++ Public procedure
44  !
45  public :: atmos_phy_tb_d1980_setup
46  public :: atmos_phy_tb_d1980
47 
48  !-----------------------------------------------------------------------------
49  !
50  !++ Public parameters & variables
51  !
52  !-----------------------------------------------------------------------------
53  !
54  !++ Private procedure
55  !
56  !-----------------------------------------------------------------------------
57  !
58  !++ Private parameters & variables
59  !
60  real(RP), private, parameter :: c1 = 0.19_rp
61  real(RP), private, parameter :: c2 = 0.51_rp
62 
63  real(RP), private, parameter :: oneoverthree = 1.0_rp / 3.0_rp
64  real(RP), private, parameter :: twooverthree = 2.0_rp / 3.0_rp
65  real(RP), private, parameter :: fouroverthree = 4.0_rp / 3.0_rp
66 
67  real(RP), private, allocatable :: delta(:,:,:) ! (dx*dy*dz)^(1/3)
68 
69  real(RP), private :: atmos_phy_tb_d1980_km_max = 1000.0_rp
70  logical, private :: atmos_phy_tb_d1980_implicit = .false.
71 
72  !-----------------------------------------------------------------------------
73 contains
74  !-----------------------------------------------------------------------------
75  subroutine atmos_phy_tb_d1980_setup( &
76  TYPE_TB, &
77  CDZ, CDX, CDY, &
78  CZ )
79  use scale_process, only: &
81  implicit none
82 
83  character(len=*), intent(in) :: TYPE_TB
84 
85  real(RP), intent(in) :: CDZ(ka)
86  real(RP), intent(in) :: CDX(ia)
87  real(RP), intent(in) :: CDY(ja)
88  real(RP), intent(in) :: CZ (ka,ia,ja)
89 
90  namelist / param_atmos_phy_tb_d1980 / &
91  atmos_phy_tb_d1980_km_max, &
92  atmos_phy_tb_d1980_implicit
93  integer :: k, i, j
94 
95  integer :: ierr
96  !---------------------------------------------------------------------------
97 
98  if( io_l ) write(io_fid_log,*)
99  if( io_l ) write(io_fid_log,*) '++++++ Module[TURBULENCE] / Categ[ATMOS PHYSICS] / Origin[SCALElib]'
100  if( io_l ) write(io_fid_log,*) '+++ 1.5th TKE Model'
101 
102  if ( type_tb /= 'D1980' ) then
103  write(*,*) 'xxx ATMOS_PHY_TB_TYPE is not D1980. Check!'
104  call prc_mpistop
105  endif
106 
107  !--- read namelist
108  rewind(io_fid_conf)
109  read(io_fid_conf,nml=param_atmos_phy_tb_d1980,iostat=ierr)
110  if( ierr < 0 ) then !--- missing
111  if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
112  elseif( ierr > 0 ) then !--- fatal error
113  write(*,*) 'xxx Not appropriate names in namelist PARAM_ATMOS_PHY_TB_D1980. Check!'
114  call prc_mpistop
115  endif
116  if( io_lnml ) write(io_fid_log,nml=param_atmos_phy_tb_d1980)
117 
118  allocate( delta(ka,ia,ja) )
119 
120 #ifdef DEBUG
121  delta(:,:,:) = undef
122 #endif
123  do j = js-1, je+1
124  do i = is-1, ie+1
125  do k = ks, ke
126  delta(k,i,j) = ( cdz(k) * cdx(i) * cdy(j) )**(1.0_rp/3.0_rp)
127  enddo
128  enddo
129  enddo
130 
131  return
132  end subroutine atmos_phy_tb_d1980_setup
133 
134  !-----------------------------------------------------------------------------
135  subroutine atmos_phy_tb_d1980( &
136  qflx_sgs_momz, qflx_sgs_momx, qflx_sgs_momy, &
137  qflx_sgs_rhot, qflx_sgs_rhoq, &
138  tke, tke_t, Km, Ri, Pr, N2, &
139  MOMZ, MOMX, MOMY, RHOT, DENS, QTRC, &
140  SFLX_MW, SFLX_MU, SFLX_MV, SFLX_SH, SFLX_QV, &
141  GSQRT, J13G, J23G, J33G, MAPF, dt )
143  use scale_tracer
144  use scale_const, only: &
145  eps => const_eps, &
146  grav => const_grav
147  use scale_grid, only: &
148  fdz => grid_fdz, &
149  fdx => grid_fdx, &
150  fdy => grid_fdy, &
151  rcdz => grid_rcdz, &
152  rcdx => grid_rcdx, &
153  rcdy => grid_rcdy, &
154  rfdz => grid_rfdz, &
155  rfdx => grid_rfdx, &
156  rfdy => grid_rfdy
157  use scale_gridtrans, only: &
158  i_xyz, &
159  i_xyw, &
160  i_uyw, &
161  i_xvw, &
162  i_uyz, &
163  i_xvz, &
164  i_uvz, &
165  i_xy, &
166  i_uy, &
167  i_xv, &
168  i_uv
169  use scale_atmos_phy_tb_common, only: &
170  calc_strain_tensor => atmos_phy_tb_calc_strain_tensor, &
171  diffusion_solver => atmos_phy_tb_diffusion_solver, &
172  calc_tend_momz => atmos_phy_tb_calc_tend_momz, &
173  calc_tend_momx => atmos_phy_tb_calc_tend_momx, &
174  calc_tend_momy => atmos_phy_tb_calc_tend_momy, &
175  calc_tend_phi => atmos_phy_tb_calc_tend_phi, &
176  calc_flux_phi => atmos_phy_tb_calc_flux_phi
177  implicit none
178 
179  ! SGS flux
180  real(RP), intent(out) :: qflx_sgs_momz(ka,ia,ja,3)
181  real(RP), intent(out) :: qflx_sgs_momx(ka,ia,ja,3)
182  real(RP), intent(out) :: qflx_sgs_momy(ka,ia,ja,3)
183  real(RP), intent(out) :: qflx_sgs_rhot(ka,ia,ja,3)
184  real(RP), intent(out) :: qflx_sgs_rhoq(ka,ia,ja,3,qa)
185 
186  real(RP), intent(inout) :: tke(ka,ia,ja) ! TKE
187  real(RP), intent(out) :: tke_t(ka,ia,ja) ! tendency TKE
188  real(RP), intent(out) :: Km (ka,ia,ja) ! eddy viscosity (center)
189  real(RP), intent(out) :: Ri (ka,ia,ja) ! Richardson number
190  real(RP), intent(out) :: Pr (ka,ia,ja) ! Prantle number
191  real(RP), intent(out) :: N2 (ka,ia,ja) ! squared Brunt-Vaisala frequency
192 
193  real(RP), intent(in) :: MOMZ(ka,ia,ja)
194  real(RP), intent(in) :: MOMX(ka,ia,ja)
195  real(RP), intent(in) :: MOMY(ka,ia,ja)
196  real(RP), intent(in) :: RHOT(ka,ia,ja)
197  real(RP), intent(in) :: DENS(ka,ia,ja)
198  real(RP), intent(in) :: QTRC(ka,ia,ja,qa)
199 
200  real(RP), intent(in) :: SFLX_MW(ia,ja)
201  real(RP), intent(in) :: SFLX_MU(ia,ja)
202  real(RP), intent(in) :: SFLX_MV(ia,ja)
203  real(RP), intent(in) :: SFLX_SH(ia,ja)
204  real(RP), intent(in) :: SFLX_QV(ia,ja)
205 
206  real(RP), intent(in) :: GSQRT (ka,ia,ja,7)
207  real(RP), intent(in) :: J13G (ka,ia,ja,7)
208  real(RP), intent(in) :: J23G (ka,ia,ja,7)
209  real(RP), intent(in) :: J33G
210  real(RP), intent(in) :: MAPF(ia,ja,2,4)
211  real(DP), intent(in) :: dt
212 
213  ! diagnostic variables
214  real(RP) :: POTT (ka,ia,ja)
215 
216  ! deformation rate tensor
217  real(RP) :: S33_C (ka,ia,ja) ! (cell center)
218  real(RP) :: S11_C (ka,ia,ja)
219  real(RP) :: S22_C (ka,ia,ja)
220  real(RP) :: S31_C (ka,ia,ja)
221  real(RP) :: S12_C (ka,ia,ja)
222  real(RP) :: S23_C (ka,ia,ja)
223  real(RP) :: S12_Z (ka,ia,ja) ! (z edge or x-y plane)
224  real(RP) :: S23_X (ka,ia,ja) ! (x edge or y-z plane)
225  real(RP) :: S31_Y (ka,ia,ja) ! (y edge or z-x plane)
226  real(RP) :: S2 (ka,ia,ja) ! |S|^2
227 
228  real(RP) :: Kh (ka,ia,ja) ! eddy diffusion
229  real(RP) :: l (ka,ia,ja) ! mixing length
230 
231  real(RP) :: qflx_tke(ka,ia,ja,3)
232 
233  real(RP) :: TEND(ka,ia,ja)
234  real(RP) :: a(ka,ia,ja)
235  real(RP) :: b(ka,ia,ja)
236  real(RP) :: c(ka,ia,ja)
237  real(RP) :: d(ka)
238  real(RP) :: ap
239 
240  integer :: IIS, IIE
241  integer :: JJS, JJE
242 
243  integer :: k, i, j, iq
244  !---------------------------------------------------------------------------
245 
246  if( io_l ) write(io_fid_log,*) '*** Physics step: Turbulence(D1980)'
247 
248 #ifdef DEBUG
249  qflx_sgs_momz(:,:,:,:) = undef
250  qflx_sgs_momx(:,:,:,:) = undef
251  qflx_sgs_momy(:,:,:,:) = undef
252  qflx_sgs_rhot(:,:,:,:) = undef
253  qflx_sgs_rhoq(:,:,:,:,:) = undef
254 
255  tke_t(:,:,:) = undef
256 
257  pr(:,:,:) = undef
258  ri(:,:,:) = undef
259  kh(:,:,:) = undef
260  km(:,:,:) = undef
261  n2(:,:,:) = undef
262 
263  pott(:,:,:) = undef
264 #endif
265 
266  ! potential temperature
267  do j = js-1, je+1
268  do i = is-1, ie+1
269  do k = ks, ke
270 #ifdef DEBUG
271  call check( __line__, rhot(k,i,j) )
272  call check( __line__, dens(k,i,j) )
273 #endif
274  pott(k,i,j) = rhot(k,i,j) / dens(k,i,j)
275  enddo
276  enddo
277  enddo
278 #ifdef DEBUG
279  i = iundef; j = iundef; k = iundef
280 #endif
281 
282  !##### Start Upadate #####
283 
284  call calc_strain_tensor( &
285  s33_c, s11_c, s22_c, & ! (out)
286  s31_c, s12_c, s23_c, & ! (out)
287  s12_z, s23_x, s31_y, & ! (out)
288  s2 , & ! (out)
289  dens, momz, momx, momy, & ! (in)
290  gsqrt, j13g, j23g, j33g, mapf ) ! (in)
291 
292  do jjs = js, je, jblock
293  jje = jjs+jblock-1
294  do iis = is, ie, iblock
295  iie = iis+iblock-1
296 
297  ! Ri = N^2 / |S|^2, N^2 = g / theta * dtheta/dz
298  do j = jjs-1, jje+1
299  do i = iis-1, iie+1
300  do k = ks+1, ke-1
301 #ifdef DEBUG
302  call check( __line__, pott(k+1,i,j) )
303  call check( __line__, pott(k,i,j) )
304  call check( __line__, pott(k-1,i,j) )
305  call check( __line__, fdz(k) )
306  call check( __line__, fdz(k-1) )
307  call check( __line__, s2(k,i,j) )
308 #endif
309  n2(k,i,j) = grav * ( pott(k+1,i,j) - pott(k-1,i,j) ) * j33g &
310  / ( ( fdz(k) + fdz(k-1) ) * gsqrt(k,i,j,i_xyz) * pott(k,i,j) )
311  ri(k,i,j) = n2(k,i,j) / s2(k,i,j)
312  enddo
313  enddo
314  enddo
315 #ifdef DEBUG
316  i = iundef; j = iundef; k = iundef
317 #endif
318  do j = jjs-1, jje+1
319  do i = iis-1, iie+1
320 #ifdef DEBUG
321  call check( __line__, pott(ks+1,i,j) )
322  call check( __line__, pott(ks,i,j) )
323  call check( __line__, rfdz(ks) )
324  call check( __line__, s2(ks,i,j) )
325 #endif
326  n2(ks,i,j) = grav * ( pott(ks+1,i,j) - pott(ks,i,j) ) * j33g &
327  / ( fdz(ks) * gsqrt(ks,i,j,i_xyz) * pott(ks,i,j) )
328  ri(ks,i,j) = grav * ( pott(ks+1,i,j) - pott(ks,i,j) ) * j33g * rfdz(ks) &
329  / ( gsqrt(ks,i,j,i_xyz) * pott(ks,i,j) * s2(ks,i,j) )
330  enddo
331  enddo
332 #ifdef DEBUG
333  i = iundef; j = iundef; k = iundef
334 #endif
335  do j = jjs-1, jje+1
336  do i = iis-1, iie+1
337 #ifdef DEBUG
338  call check( __line__, pott(ke,i,j) )
339  call check( __line__, pott(ke-1,i,j) )
340  call check( __line__, rfdz(ke-1) )
341  call check( __line__, s2(ke,i,j) )
342 #endif
343  n2(ke,i,j) = grav * ( pott(ke,i,j) - pott(ke-1,i,j) ) * j33g &
344  / ( fdz(ke-1) * gsqrt(ke,i,j,i_xyz) * pott(ke,i,j) )
345  ri(ke,i,j) = grav * ( pott(ke,i,j) - pott(ke-1,i,j) ) * j33g * rfdz(ke-1) &
346  / ( gsqrt(ke,i,j,i_xyz) * pott(ke,i,j) * s2(ke,i,j) )
347  enddo
348  enddo
349 #ifdef DEBUG
350  i = iundef; j = iundef; k = iundef
351 #endif
352 
353  ! mixing length
354  do j = jjs-1, jje+1
355  do i = iis-1, iie+1
356  do k = ks, ke
357  if ( n2(k,i,j) > 1e-10_rp ) then
358  l(k,i,j) = max( min( 0.76_rp * sqrt(tke(k,i,j)/n2(k,i,j)), delta(k,i,j) ), 1e-10_rp )
359  else
360  l(k,i,j) = delta(k,i,j)
361  end if
362  enddo
363  enddo
364  enddo
365 #ifdef DEBUG
366  i = iundef; j = iundef; k = iundef
367 #endif
368 
369  do j = jjs-1, jje+1
370  do i = iis-1, iie+1
371  do k = ks, ke
372  km(k,i,j) = min( 0.10_rp * l(k,i,j) * sqrt(tke(k,i,j)), atmos_phy_tb_d1980_km_max )
373  pr(k,i,j) = 1.0_rp / ( 1.0_rp + 2.0_rp * l(k,i,j) / delta(k,i,j) )
374  kh(k,i,j) = km(k,i,j) / pr(k,i,j)
375  enddo
376  enddo
377  enddo
378 #ifdef DEBUG
379  i = iundef; j = iundef; k = iundef
380 #endif
381 
382  !##### momentum equation (z) #####
383  ! (cell center)
384  do j = jjs, jje
385  do i = iis, iie
386  do k = ks+1, ke-1
387 #ifdef DEBUG
388  call check( __line__, dens(k,i,j) )
389  call check( __line__, km(k,i,j) )
390  call check( __line__, s33_c(k,i,j) )
391 #endif
392  qflx_sgs_momz(k,i,j,zdir) = dens(k,i,j) * ( - 2.0_rp * km(k,i,j) * s33_c(k,i,j) )
393  enddo
394  enddo
395  enddo
396 #ifdef DEBUG
397  i = iundef; j = iundef; k = iundef
398 #endif
399  do j = jjs, jje
400  do i = iis, iie
401  qflx_sgs_momz(ks,i,j,zdir) = 0.0_rp ! just above bottom boundary
402  qflx_sgs_momz(ke,i,j,zdir) = 0.0_rp ! just below top boundary
403  enddo
404  enddo
405 #ifdef DEBUG
406  i = iundef; j = iundef; k = iundef
407 #endif
408  ! (y edge)
409  do j = jjs, jje
410  do i = iis-1, iie
411  do k = ks, ke-1
412 #ifdef DEBUG
413  call check( __line__, dens(k,i,j) )
414  call check( __line__, dens(k,i+1,j) )
415  call check( __line__, dens(k+1,i,j) )
416  call check( __line__, dens(k+1,i+1,j) )
417  call check( __line__, km(k,i,j) )
418  call check( __line__, km(k,i+1,j) )
419  call check( __line__, km(k+1,i,j) )
420  call check( __line__, km(k+1,i+1,j) )
421  call check( __line__, s31_y(k,i,j) )
422 #endif
423  qflx_sgs_momz(k,i,j,xdir) = - 0.125_rp & ! 2.0 / 4 / 4
424  * ( dens(k,i,j)+dens(k+1,i,j)+dens(k,i+1,j)+dens(k+1,i+1,j) ) &
425  * ( km(k,i,j)+km(k+1,i,j)+km(k,i+1,j)+km(k+1,i+1,j)) &
426  * s31_y(k,i,j)
427  enddo
428  enddo
429  enddo
430 #ifdef DEBUG
431  i = iundef; j = iundef; k = iundef
432 #endif
433  ! (x edge)
434  do j = jjs-1, jje
435  do i = iis, iie
436  do k = ks, ke-1
437 #ifdef DEBUG
438  call check( __line__, dens(k,i,j) )
439  call check( __line__, dens(k,i,j+1) )
440  call check( __line__, dens(k+1,i,j) )
441  call check( __line__, dens(k+1,i,j+1) )
442  call check( __line__, km(k,i,j) )
443  call check( __line__, km(k,i,j+1) )
444  call check( __line__, km(k+1,i,j) )
445  call check( __line__, km(k+1,i,j+1) )
446  call check( __line__, s23_x(k,i,j) )
447 #endif
448  qflx_sgs_momz(k,i,j,ydir) = - 0.125_rp & ! 2/4/4
449  * ( dens(k,i,j)+dens(k+1,i,j)+dens(k,i,j+1)+dens(k+1,i,j+1) ) &
450  * ( km(k,i,j)+km(k+1,i,j)+km(k,i,j+1)+km(k+1,i,j+1) ) &
451  * s23_x(k,i,j)
452  enddo
453  enddo
454  enddo
455 #ifdef DEBUG
456  i = iundef; j = iundef; k = iundef
457 #endif
458 
459  if ( atmos_phy_tb_d1980_implicit ) then
460  call calc_tend_momz( tend, & ! (out)
461  qflx_sgs_momz, & ! (in)
462  gsqrt, j13g, j23g, j33g, mapf, & ! (in)
463  iis, iie, jjs, jje ) ! (in)
464 
465  do j = jjs, jje
466  do i = iis, iie
467 
468  ap = - fouroverthree * dt &
469  * dens(ks+1,i,j)*km(ks+1,i,j) &
470  * rcdz(ks+1) / gsqrt(ks+1,i,j,i_xyz)
471  a(ks,i,j) = ap * rfdz(ks) / gsqrt(ks,i,j,i_xyw)
472  c(ks,i,j) = 0.0_rp
473  b(ks,i,j) = - a(ks,i,j) + 0.5_rp * ( dens(ks,i,j)+dens(ks+1,i,j) )
474  do k = ks+1, ke-2
475  c(k,i,j) = ap * rfdz(k+1) / gsqrt(k+1,i,j,i_xyw)
476  ap = - fouroverthree * dt &
477  * dens(k+1,i,j)*km(k+1,i,j) &
478  * rcdz(k+1) / gsqrt(k+1,i,j,i_xyz)
479  a(k,i,j) = ap * rfdz(k) / gsqrt(k,i,j,i_xyw)
480  b(k,i,j) = - a(k,i,j) - c(k,i,j) + 0.5_rp * ( dens(k,i,j)+dens(k+1,i,j) )
481  end do
482  a(ke-1,i,j) = 0.0_rp
483  c(ke-1,i,j) = ap * rfdz(ke) / gsqrt(ke,i,j,i_xyw)
484  b(ke-1,i,j) = - c(ke-1,i,j) + 0.5_rp * ( dens(ke-1,i,j)+dens(ke,i,j) )
485 
486  do k = ks, ke-1
487  d(k) = tend(k,i,j)
488  end do
489 
490  call diffusion_solver( &
491  tend(:,i,j), & ! (out)
492  a(:,i,j), b(:,i,j), c(:,i,j), d, & ! (in)
493  ke-1 ) ! (in)
494 
495  do k = ks+1, ke-1
496  qflx_sgs_momz(k,i,j,zdir) = qflx_sgs_momz(k,i,j,zdir) &
497  - fouroverthree * dens(k,i,j) * km(k,i,j) * dt &
498  * ( tend(k,i,j) - tend(k-1,i,j) ) * rcdz(k) / gsqrt(k,i,j,i_xyz)
499  end do
500 
501  end do
502  end do
503 
504  end if
505 
506  !##### momentum equation (x) #####
507  ! (y edge)
508  do j = jjs, jje
509  do i = iis, iie
510  do k = ks, ke-1
511 #ifdef DEBUG
512  call check( __line__, dens(k,i,j) )
513  call check( __line__, dens(k,i+1,j) )
514  call check( __line__, dens(k+1,i,j) )
515  call check( __line__, dens(k+1,i+1,j) )
516  call check( __line__, km(k,i,j) )
517  call check( __line__, km(k,i+1,j) )
518  call check( __line__, km(k+1,i,j) )
519  call check( __line__, km(k+1,i+1,j) )
520  call check( __line__, s31_y(k,i,j) )
521 #endif
522  qflx_sgs_momx(k,i,j,zdir) = - 0.125_rp & ! 2/4/4
523  * ( dens(k,i,j)+dens(k+1,i,j)+dens(k,i+1,j)+dens(k+1,i+1,j) ) &
524  * ( km(k,i,j)+km(k+1,i,j)+km(k,i+1,j)+km(k+1,i+1,j) ) &
525  * s31_y(k,i,j)
526  enddo
527  enddo
528  enddo
529 #ifdef DEBUG
530  i = iundef; j = iundef; k = iundef
531 #endif
532  do j = jjs, jje
533  do i = iis, iie
534  qflx_sgs_momx(ks-1,i,j,zdir) = 0.0_rp ! bottom boundary
535  qflx_sgs_momx(ke ,i,j,zdir) = 0.0_rp ! top boundary
536  enddo
537  enddo
538 #ifdef DEBUG
539  i = iundef; j = iundef; k = iundef
540 #endif
541  ! (cell center)
542  do j = jjs, jje
543  do i = iis, iie+1
544  do k = ks, ke
545 #ifdef DEBUG
546  call check( __line__, dens(k,i,j) )
547  call check( __line__, km(k,i,j) )
548  call check( __line__, s11_c(k,i,j) )
549 #endif
550  qflx_sgs_momx(k,i,j,xdir) = dens(k,i,j) * ( - 2.0_rp * km(k,i,j) * s11_c(k,i,j) )
551  enddo
552  enddo
553  enddo
554 #ifdef DEBUG
555  i = iundef; j = iundef; k = iundef
556 #endif
557  ! (z edge)
558  do j = jjs-1, jje
559  do i = iis, iie
560  do k = ks, ke
561 #ifdef DEBUG
562  call check( __line__, dens(k,i,j) )
563  call check( __line__, dens(k,i+1,j) )
564  call check( __line__, dens(k,i,j+1) )
565  call check( __line__, dens(k,i+1,j+1) )
566  call check( __line__, km(k,i,j) )
567  call check( __line__, km(k,i+1,j) )
568  call check( __line__, km(k,i,j+1) )
569  call check( __line__, km(k,i+1,j+1) )
570  call check( __line__, s12_z(k,i,j) )
571 #endif
572  qflx_sgs_momx(k,i,j,ydir) = - 0.125_rp & ! 2/4/4
573  * ( dens(k,i,j)+dens(k,i+1,j)+dens(k,i,j+1)+dens(k,i+1,j+1) ) &
574  * ( km(k,i,j)+km(k,i+1,j)+km(k,i,j+1)+km(k,i+1,j+1) ) &
575  * s12_z(k,i,j)
576  enddo
577  enddo
578  enddo
579 #ifdef DEBUG
580  i = iundef; j = iundef; k = iundef
581 #endif
582 
583  if ( atmos_phy_tb_d1980_implicit ) then
584  call calc_tend_momx( tend, & ! (out)
585  qflx_sgs_momx, & ! (in)
586  gsqrt, j13g, j23g, j33g, mapf, & ! (in)
587  iis, iie, jjs, jje ) ! (in)
588 
589  do j = jjs, jje
590  do i = iis, iie
591 
592  ap = - dt * 0.25_rp * ( dens(ks ,i ,j)*km(ks ,i ,j) &
593  + dens(ks+1,i ,j)*km(ks+1,i ,j) &
594  + dens(ks ,i+1,j)*km(ks ,i+1,j) &
595  + dens(ks+1,i+1,j)*km(ks+1,i+1,j) ) &
596  * rfdz(ks) / gsqrt(ks,i,j,i_uyw)
597  a(ks,i,j) = ap * rcdz(ks) / gsqrt(ks,i,j,i_uyz)
598  c(ks,i,j) = 0.0_rp
599  b(ks,i,j) = - a(ks,i,j) + 0.5_rp * ( dens(ks,i,j)+dens(ks,i+1,j) )
600  do k = ks+1, ke-1
601  c(k,i,j) = ap * rcdz(k) / gsqrt(k,i,j,i_uyz)
602  ap = - dt * 0.25_rp * ( dens(k ,i ,j)*km(k ,i ,j) &
603  + dens(k+1,i ,j)*km(k+1,i ,j) &
604  + dens(k ,i+1,j)*km(k ,i+1,j) &
605  + dens(k+1,i+1,j)*km(k+1,i+1,j) ) &
606  * rfdz(k) / gsqrt(k,i,j,i_uyw)
607  a(k,i,j) = ap * rcdz(k) / gsqrt(k,i,j,i_uyz)
608  b(k,i,j) = - a(k,i,j) - c(k,i,j) + 0.5_rp * ( dens(k,i,j)+dens(k,i+1,j) )
609  end do
610  a(ke,i,j) = 0.0_rp
611  c(ke,i,j) = ap * rcdz(ke) / gsqrt(ke,i,j,i_uyz)
612  b(ke,i,j) = - c(ke,i,j) + 0.5_rp * ( dens(ke,i,j)+dens(ke,i+1,j) )
613 
614  do k = ks, ke
615  d(k) = tend(k,i,j)
616  end do
617 
618  call diffusion_solver( &
619  tend(:,i,j), & ! (out)
620  a(:,i,j), b(:,i,j), c(:,i,j), d, & ! (in)
621  ke ) ! (in)
622 
623  do k = ks, ke-1
624  qflx_sgs_momx(k,i,j,zdir) = qflx_sgs_momx(k,i,j,zdir) &
625  - 0.25_rp * ( dens(k ,i ,j)*km(k ,i ,j) &
626  + dens(k+1,i ,j)*km(k+1,i ,j) &
627  + dens(k ,i+1,j)*km(k ,i+1,j) &
628  + dens(k+1,i+1,j)*km(k+1,i+1,j) ) &
629  * dt * ( tend(k+1,i,j) - tend(k,i,j) ) * rfdz(k) / gsqrt(k,i,j,i_uyw)
630  end do
631 
632  end do
633  end do
634 
635  end if
636 
637  !##### momentum equation (y) #####
638  ! (x edge)
639  do j = jjs, jje
640  do i = iis, iie
641  do k = ks, ke-1
642 #ifdef DEBUG
643  call check( __line__, dens(k,i,j) )
644  call check( __line__, dens(k,i,j+1) )
645  call check( __line__, dens(k+1,i,j) )
646  call check( __line__, dens(k+1,i,j+1) )
647  call check( __line__, km(k,i,j) )
648  call check( __line__, km(k,i,j+1) )
649  call check( __line__, km(k+1,i,j) )
650  call check( __line__, km(k+1,i,j+1) )
651  call check( __line__, s23_x(k,i,j) )
652 #endif
653  qflx_sgs_momy(k,i,j,zdir) = - 0.125_rp & ! 2/4/4
654  * ( dens(k,i,j)+dens(k+1,i,j)+dens(k,i,j+1)+dens(k+1,i,j+1) ) &
655  * ( km(k,i,j)+km(k+1,i,j)+km(k,i,j+1)+km(k+1,i,j+1) ) &
656  * s23_x(k,i,j)
657  enddo
658  enddo
659  enddo
660 #ifdef DEBUG
661  i = iundef; j = iundef; k = iundef
662 #endif
663  do j = jjs, jje
664  do i = iis, iie
665  qflx_sgs_momy(ks-1,i,j,zdir) = 0.0_rp ! bottom boundary
666  qflx_sgs_momy(ke ,i,j,zdir) = 0.0_rp ! top boundary
667  enddo
668  enddo
669 #ifdef DEBUG
670  i = iundef; j = iundef; k = iundef
671 #endif
672 
673  ! (z edge)
674  do j = jjs, jje
675  do i = iis-1, iie
676  do k = ks, ke
677 #ifdef DEBUG
678  call check( __line__, dens(k,i,j) )
679  call check( __line__, dens(k,i+1,j) )
680  call check( __line__, dens(k,i,j+1) )
681  call check( __line__, dens(k,i+1,j+1) )
682  call check( __line__, km(k,i,j) )
683  call check( __line__, km(k,i+1,j) )
684  call check( __line__, km(k,i,j+1) )
685  call check( __line__, km(k,i+1,j+1) )
686  call check( __line__, s12_z(k,i,j) )
687 #endif
688  qflx_sgs_momy(k,i,j,xdir) = - 0.125_rp & !
689  * ( dens(k,i,j)+dens(k,i+1,j)+dens(k,i,j+1)+dens(k,i+1,j+1) ) &
690  * ( km(k,i,j)+km(k,i+1,j)+km(k,i,j+1)+km(k,i+1,j+1) ) &
691  * s12_z(k,i,j)
692  enddo
693  enddo
694  enddo
695 #ifdef DEBUG
696  i = iundef; j = iundef; k = iundef
697 #endif
698 
699  ! (z-x plane)
700  do j = jjs, jje+1
701  do i = iis, iie
702  do k = ks, ke
703 #ifdef DEBUG
704  call check( __line__, dens(k,i,j) )
705  call check( __line__, km(k,i,j) )
706  call check( __line__, s22_c(k,i,j) )
707 #endif
708  qflx_sgs_momy(k,i,j,ydir) = dens(k,i,j) * ( - 2.0_rp * km(k,i,j) * s22_c(k,i,j) )
709  enddo
710  enddo
711  enddo
712 #ifdef DEBUG
713  i = iundef; j = iundef; k = iundef
714 #endif
715 
716  if ( atmos_phy_tb_d1980_implicit ) then
717  call calc_tend_momy( tend, & ! (out)
718  qflx_sgs_momy, & ! (in)
719  gsqrt, j13g, j23g, j33g, mapf, & ! (in)
720  iis, iie, jjs, jje ) ! (in)
721 
722  do j = jjs, jje
723  do i = iis, iie
724 
725  ap = - dt * 0.25_rp * ( dens(ks ,i,j )*km(ks ,i,j ) &
726  + dens(ks+1,i,j )*km(ks+1,i,j ) &
727  + dens(ks ,i,j+1)*km(ks ,i,j+1) &
728  + dens(ks+1,i,j+1)*km(ks+1,i,j+1) ) &
729  * rfdz(ks) / gsqrt(ks,i,j,i_xvw)
730  a(ks,i,j) = ap * rcdz(ks) / gsqrt(ks,i,j,i_xvz)
731  c(ks,i,j) = 0.0_rp
732  b(ks,i,j) = - a(ks,i,j) + 0.5_rp * ( dens(ks,i,j)+dens(ks,i,j+1) )
733  do k = ks+1, ke-1
734  c(k,i,j) = ap * rcdz(k) / gsqrt(k,i,j,i_xvz)
735  ap = - dt * 0.25_rp * ( dens(k ,i,j )*km(k ,i,j ) &
736  + dens(k+1,i,j )*km(k+1,i,j ) &
737  + dens(k ,i,j+1)*km(k ,i,j+1) &
738  + dens(k+1,i,j+1)*km(k+1,i,j+1) ) &
739  * rfdz(k) / gsqrt(k,i,j,i_xvw)
740  a(k,i,j) = ap * rcdz(k) / gsqrt(k,i,j,i_xvz)
741  b(k,i,j) = - a(k,i,j) - c(k,i,j) + 0.5_rp * ( dens(k,i,j)+dens(k,i,j+1) )
742  end do
743  a(ke,i,j) = 0.0_rp
744  c(ke,i,j) = ap * rcdz(ke) / gsqrt(ke,i,j,i_xvz)
745  b(ke,i,j) = - c(ke,i,j) + 0.5_rp * ( dens(ke,i,j)+dens(ke,i,j+1) )
746 
747  do k = ks, ke
748  d(k) = tend(k,i,j)
749  end do
750 
751  call diffusion_solver( &
752  tend(:,i,j), & ! (out)
753  a(:,i,j), b(:,i,j), c(:,i,j), d, & ! (in)
754  ke ) ! (in)
755 
756  do k = ks, ke-1
757  qflx_sgs_momy(k,i,j,zdir) = qflx_sgs_momy(k,i,j,zdir) &
758  - 0.25_rp * ( dens(k ,i,j )*km(k ,i,j ) &
759  + dens(k+1,i,j )*km(k+1,i,j ) &
760  + dens(k ,i,j+1)*km(k ,i,j+1) &
761  + dens(k+1,i,j+1)*km(k+1,i,j+1) ) &
762  * dt * ( tend(k+1,i,j) - tend(k,i,j) ) * rfdz(k) / gsqrt(k,i,j,i_xvw)
763  end do
764 
765  end do
766  end do
767 
768  end if
769 
770  !##### Thermodynamic Equation #####
771 
772  if ( atmos_phy_tb_d1980_implicit ) then
773 
774  do j = jjs, jje
775  do i = iis, iie
776 
777  ap = - dt * 0.25_rp * ( dens(ks,i,j)+dens(ks+1,i,j) ) &
778  * ( kh(ks,i,j)+kh(ks+1,i,j) ) &
779  * rfdz(ks) / gsqrt(ks,i,j,i_xyw)
780  a(ks,i,j) = ap * rcdz(ks) / gsqrt(ks,i,j,i_xyz)
781  c(ks,i,j) = 0.0_rp
782  b(ks,i,j) = - a(ks,i,j) + dens(ks,i,j)
783  do k = ks+1, ke-1
784  c(k,i,j) = ap * rcdz(k) / gsqrt(k,i,j,i_xyz)
785  ap = - dt * 0.25_rp * ( dens(k,i,j)+dens(k+1,i,j) ) &
786  * ( kh(k,i,j)+kh(k+1,i,j) ) &
787  * rfdz(k) / gsqrt(k,i,j,i_xyw)
788  a(k,i,j) = ap * rcdz(k) / gsqrt(k,i,j,i_xyz)
789  b(k,i,j) = - a(k,i,j) - c(k,i,j) + dens(k,i,j)
790  end do
791  a(ke,i,j) = 0.0_rp
792  c(ke,i,j) = ap * rcdz(ke) / gsqrt(ke,i,j,i_xyz)
793  b(ke,i,j) = - c(ke,i,j) + dens(ke,i,j)
794 
795  end do
796  end do
797 
798  end if
799 
800  call calc_flux_phi( &
801  qflx_sgs_rhot, &
802  dens, pott, kh, 1.0_rp, &
803  gsqrt, j13g, j23g, j33g, mapf, &
804  a, b, c, dt, &
805  atmos_phy_tb_d1980_implicit, &
806  iis, iie, jjs, jje )
807 
808  enddo
809  enddo
810 
811 
812  !##### Tracers #####
813  do iq = 1, qa
814 
815  do jjs = js, je, jblock
816  jje = jjs+jblock-1
817  do iis = is, ie, iblock
818  iie = iis+iblock-1
819 
820  call calc_flux_phi( &
821  qflx_sgs_rhoq(:,:,:,:,iq), &
822  dens, qtrc(:,:,:,iq), kh, 1.0_rp, &
823  gsqrt, j13g, j23g, j33g, mapf, &
824  a, b, c, dt, &
825  atmos_phy_tb_d1980_implicit, &
826  iis, iie, jjs, jje )
827 
828  enddo
829  enddo
830 #ifdef DEBUG
831  iis = iundef; iie = iundef; jjs = iundef; jje = iundef
832 #endif
833 
834  enddo ! scalar quantities loop
835 #ifdef DEBUG
836  iq = iundef
837 #endif
838 
839  ! TKE
840  do jjs = js, je, jblock
841  jje = jjs+jblock-1
842  do iis = is, ie, iblock
843  iie = iis+iblock-1
844 
845  if ( atmos_phy_tb_d1980_implicit ) then
846 
847  do j = jjs, jje
848  do i = iis, iie
849 
850  ap = - dt * 0.25_rp * ( dens(ks,i,j)+dens(ks+1,i,j) ) &
851  * 2.0_rp * ( km(ks,i,j)+km(ks+1,i,j) ) &
852  * rfdz(ks) / gsqrt(ks,i,j,i_xyw)
853  a(ks,i,j) = ap * rcdz(ks) / gsqrt(ks,i,j,i_xyz)
854  c(ks,i,j) = 0.0_rp
855  b(ks,i,j) = - a(ks,i,j) + dens(ks,i,j)
856  do k = ks+1, ke-1
857  c(k,i,j) = ap * rcdz(k) / gsqrt(k,i,j,i_xyz)
858  ap = - dt * 0.25_rp * ( dens(k,i,j)+dens(k+1,i,j) ) &
859  * 2.0_rp * ( km(k,i,j)+km(k+1,i,j) ) &
860  * rfdz(k) / gsqrt(k,i,j,i_xyw)
861  a(k,i,j) = ap * rcdz(k) / gsqrt(k,i,j,i_xyz)
862  b(k,i,j) = - a(k,i,j) - c(k,i,j) + dens(k,i,j)
863  end do
864  a(ke,i,j) = 0.0_rp
865  c(ke,i,j) = ap * rcdz(ke) / gsqrt(ke,i,j,i_xyz)
866  b(ke,i,j) = - c(ke,i,j) + dens(ke,i,j)
867 
868  end do
869  end do
870 
871  end if
872 
873  call calc_flux_phi( &
874  qflx_tke, &
875  dens, tke, km, 2.0_rp, &
876  gsqrt, j13g, j23g, j33g, mapf, &
877  a, b, c, dt, &
878  atmos_phy_tb_d1980_implicit, &
879  iis, iie, jjs, jje )
880 
881  do j = jjs, jje
882  do i = iis, iie
883  do k = ks, ke-1
884  qflx_tke(k,i,j,zdir) = qflx_tke(k,i,j,zdir) &
885  + 0.5_rp * ( momz(k,i,j) * ( tke(k+1,i,j)+tke(k,i,j) ) &
886  - abs(momz(k,i,j)) * ( tke(k+1,i,j)-tke(k,i,j) ) )
887  end do
888  end do
889  end do
890 
891  do j = jjs, jje
892  do i = iis-1, iie
893  do k = ks, ke
894  qflx_tke(k,i,j,xdir) = qflx_tke(k,i,j,xdir) &
895  + 0.5_rp * ( momx(k,i,j) * ( tke(k,i+1,j)+tke(k,i,j) ) &
896  - abs(momx(k,i,j)) * ( tke(k,i+1,j)-tke(k,i,j) ) )
897  end do
898  end do
899  end do
900 
901  do j = jjs-1, jje
902  do i = iis, iie
903  do k = ks, ke
904  qflx_tke(k,i,j,ydir) = qflx_tke(k,i,j,ydir) &
905  + 0.5_rp * ( momy(k,i,j) * ( tke(k,i,j+1)+tke(k,i,j) ) &
906  - abs(momy(k,i,j)) * ( tke(k,i,j+1)-tke(k,i,j) ) )
907  end do
908  end do
909  end do
910 
911  call calc_tend_phi( tke_t, & ! (out)
912  qflx_tke, & ! (in)
913  gsqrt, j13g, j23g, j33g, mapf, & ! (in)
914  iis, iie, jjs, jje ) ! (in)
915  do j = jjs, jje
916  do i = iis, iie
917  do k = ks, ke
918  tke_t(k,i,j) = min( &
919  tke_t(k,i,j) &
920  + km(k,i,j) * s2(k,i,j) &
921  - kh(k,i,j) * n2(k,i,j) &
922  - ( c1 + c2*l(k,i,j)/delta(k,i,j) ) * tke(k,i,j)**(1.5_rp) / l(k,i,j), &
923  tke(k,i,j)/dt )
924  enddo
925  enddo
926  enddo
927 #ifdef DEBUG
928  i = iundef; j = iundef; k = iundef
929 #endif
930 
931  end do
932  end do
933 
934 
935  return
936  end subroutine atmos_phy_tb_d1980
937 
938 end module scale_atmos_phy_tb_d1980
integer, public is
start point of inner domain: x, local
module DEBUG
Definition: scale_debug.F90:13
integer, public i_xvz
subroutine, public atmos_phy_tb_calc_tend_momz(MOMZ_t_TB, QFLX_MOMZ, GSQRT, J13G, J23G, J33G, MAPF, IIS, IIE, JJS, JJE)
integer, public je
end point of inner domain: y, local
real(rp), dimension(:), allocatable, public grid_rcdy
reciprocal of center-dy
subroutine, public prc_mpistop
Abort MPI.
subroutine, public atmos_phy_tb_diffusion_solver(phi, a, b, c, d, KE_TB)
subroutine, public atmos_phy_tb_calc_flux_phi(qflx_phi, DENS, PHI, Kh, FACT, GSQRT, J13G, J23G, J33G, MAPF, a, b, c, dt, implicit, IIS, IIE, JJS, JJE)
subroutine, public atmos_phy_tb_d1980_setup(TYPE_TB, CDZ, CDX, CDY, CZ)
real(rp), dimension(:), allocatable, public grid_fdy
y-length of grid(j+1) to grid(j) [m]
integer, public iblock
block size for cache blocking: x
subroutine, public atmos_phy_tb_calc_tend_momy(MOMY_t_TB, QFLX_MOMY, GSQRT, J13G, J23G, J33G, MAPF, IIS, IIE, JJS, JJE)
logical, public io_l
output log or not? (this process)
Definition: scale_stdio.F90:59
real(rp), dimension(:), allocatable, public grid_rcdx
reciprocal of center-dx
integer, parameter, public zdir
module STDIO
Definition: scale_stdio.F90:12
integer, parameter, public ydir
integer, public ke
end point of inner domain: z, local
integer, parameter, public xdir
integer, public qa
integer, public i_xy
real(rp), dimension(:), allocatable, public grid_rfdy
reciprocal of face-dy
subroutine, public atmos_phy_tb_d1980(qflx_sgs_momz, qflx_sgs_momx, qflx_sgs_momy, qflx_sgs_rhot, qflx_sgs_rhoq, tke, tke_t, Km, 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, public check(current_line, v)
Undefined value checker.
Definition: scale_debug.F90:58
real(rp), dimension(:), allocatable, public grid_rcdz
reciprocal of center-dz
integer, public i_xvw
real(rp), public const_undef
Definition: scale_const.F90:43
module grid index
subroutine, public atmos_phy_tb_calc_tend_phi(phi_t_TB, QFLX_phi, GSQRT, J13G, J23G, J33G, MAPF, IIS, IIE, JJS, JJE)
module TRACER
integer, public ia
of x whole cells (local, with HALO)
subroutine, public atmos_phy_tb_calc_tend_momx(MOMX_t_TB, QFLX_MOMX, GSQRT, J13G, J23G, J33G, MAPF, IIS, IIE, JJS, JJE)
module ATMOSPHERE / Physics Turbulence
module GRIDTRANS
real(rp), dimension(:), allocatable, public grid_fdz
z-length of grid(k+1) to grid(k) [m]
integer, public ka
of z whole cells (local, with HALO)
integer, public i_uy
integer, public jblock
block size for cache blocking: y
integer, public i_xyw
real(rp), public const_grav
standard acceleration of gravity [m/s2]
Definition: scale_const.F90:48
integer, public js
start point of inner domain: y, local
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)
integer, parameter, public const_undef2
undefined value (INT2)
Definition: scale_const.F90:40
module PROCESS
integer, public i_uyw
module CONSTANT
Definition: scale_const.F90:14
integer, public ks
start point of inner domain: z, local
module ATMOSPHERE / Physics Turbulence
module GRID (cartesian)
integer, public i_xv
module profiler
Definition: scale_prof.F90:10
integer, public ie
end point of inner domain: x, local
real(rp), public const_eps
small number
Definition: scale_const.F90:36
logical, public io_lnml
output log or not? (for namelist, this process)
Definition: scale_stdio.F90:60
integer, public i_uyz
module PRECISION
real(rp), dimension(:), allocatable, public grid_fdx
x-length of grid(i+1) to grid(i) [m]
real(rp), dimension(:), allocatable, public grid_rfdx
reciprocal of face-dx
real(rp), dimension(:), allocatable, public grid_rfdz
reciprocal of face-dz
integer, public io_fid_conf
Config file ID.
Definition: scale_stdio.F90:55
integer, public i_uvz
integer, public i_xyz
integer, public io_fid_log
Log file ID.
Definition: scale_stdio.F90:56
integer, public i_uv
integer, public ja
of y whole cells (local, with HALO)