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