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