SCALE-RM
Functions/Subroutines
scale_atmos_phy_tb_d1980 Module Reference

module ATMOSPHERE / Physics Turbulence More...

Functions/Subroutines

subroutine, public atmos_phy_tb_d1980_config (TYPE_TB, I_TKE_out)
 Config. More...
 
subroutine, public atmos_phy_tb_d1980_setup (CDZ, CDX, CDY, CZ)
 Setup. More...
 
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)
 

Detailed Description

module ATMOSPHERE / Physics Turbulence

Description
Sub-grid scale turbulence process 1.5th TKE model Deardorff (1980)
Author
Team SCALE
Reference
  • Deardorff, 1980: Stratocumulus-capped mixed layers derived from a three-dimensional model. Bound.-Layer Meteor., 18, 495-527
NAMELIST
  • PARAM_ATMOS_PHY_TB_D1980
    nametypedefault valuecomment
    ATMOS_PHY_TB_D1980_KM_MAX real(RP) 1000.0_RP
    ATMOS_PHY_TB_D1980_IMPLICIT logical .false.

History Output
No history output

Function/Subroutine Documentation

◆ atmos_phy_tb_d1980_config()

subroutine, public scale_atmos_phy_tb_d1980::atmos_phy_tb_d1980_config ( character(len=*), intent(in)  TYPE_TB,
integer, intent(out)  I_TKE_out 
)

Config.

Definition at line 80 of file scale_atmos_phy_tb_d1980.F90.

References scale_prc::prc_abort(), and scale_tracer::tracer_regist().

Referenced by mod_atmos_phy_tb_driver::atmos_phy_tb_driver_tracer_setup().

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
module TRACER
module PROCESS
Definition: scale_prc.F90:11
subroutine, public prc_abort
Abort Process.
Definition: scale_prc.F90:338
subroutine, public tracer_regist(QS, NQ, NAME, DESC, UNIT, CV, CP, R, ADVC, MASS)
Regist tracer.
Here is the call graph for this function:
Here is the caller graph for this function:

◆ atmos_phy_tb_d1980_setup()

subroutine, public scale_atmos_phy_tb_d1980::atmos_phy_tb_d1980_setup ( real(rp), dimension(ka), intent(in)  CDZ,
real(rp), dimension(ia), intent(in)  CDX,
real(rp), dimension(ja), intent(in)  CDY,
real(rp), dimension (ka,ia,ja), intent(in)  CZ 
)

Setup.

Definition at line 114 of file scale_atmos_phy_tb_d1980.F90.

References scale_atmos_grid_cartesc_index::ia, scale_atmos_grid_cartesc_index::ie, scale_io::io_fid_conf, scale_atmos_grid_cartesc_index::is, scale_atmos_grid_cartesc_index::ja, scale_atmos_grid_cartesc_index::je, scale_atmos_grid_cartesc_index::js, scale_atmos_grid_cartesc_index::ka, scale_atmos_grid_cartesc_index::ke, scale_atmos_grid_cartesc_index::ks, and scale_prc::prc_abort().

Referenced by mod_atmos_phy_tb_driver::atmos_phy_tb_driver_setup().

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
integer, public ia
of whole cells: x, local, with HALO
integer, public ja
of whole cells: y, local, with HALO
integer, public io_fid_conf
Config file ID.
Definition: scale_io.F90:55
integer, public is
start point of inner domain: x, local
integer, public ie
end point of inner domain: x, local
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
integer, public ks
start point of inner domain: z, local
subroutine, public prc_abort
Abort Process.
Definition: scale_prc.F90:338
integer, public js
start point of inner domain: y, local
integer, public ka
of whole cells: z, local, with HALO
Here is the call graph for this function:
Here is the caller graph for this function:

◆ atmos_phy_tb_d1980()

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

Definition at line 172 of file scale_atmos_phy_tb_d1980.F90.

References scale_atmos_grid_cartesc::atmos_grid_cartesc_fdz, scale_atmos_grid_cartesc::atmos_grid_cartesc_rcdx, scale_atmos_grid_cartesc::atmos_grid_cartesc_rcdy, scale_atmos_grid_cartesc::atmos_grid_cartesc_rcdz, scale_atmos_grid_cartesc::atmos_grid_cartesc_rfdx, scale_atmos_grid_cartesc::atmos_grid_cartesc_rfdy, scale_atmos_grid_cartesc::atmos_grid_cartesc_rfdz, scale_atmos_phy_tb_common::atmos_phy_tb_calc_flux_phi(), scale_atmos_phy_tb_common::atmos_phy_tb_calc_strain_tensor(), scale_atmos_phy_tb_common::atmos_phy_tb_calc_tend_momx(), scale_atmos_phy_tb_common::atmos_phy_tb_calc_tend_momy(), scale_atmos_phy_tb_common::atmos_phy_tb_calc_tend_momz(), scale_atmos_phy_tb_common::atmos_phy_tb_calc_tend_phi(), scale_atmos_phy_tb_common::atmos_phy_tb_diffusion_solver(), scale_debug::check(), scale_const::const_grav, scale_atmos_grid_cartesc_index::i_uyw, scale_atmos_grid_cartesc_index::i_uyz, scale_atmos_grid_cartesc_index::i_xvw, scale_atmos_grid_cartesc_index::i_xvz, scale_atmos_grid_cartesc_index::i_xyw, scale_atmos_grid_cartesc_index::i_xyz, scale_atmos_grid_cartesc_index::iblock, scale_atmos_grid_cartesc_index::ie, scale_atmos_grid_cartesc_index::is, scale_atmos_grid_cartesc_index::jblock, scale_atmos_grid_cartesc_index::je, scale_atmos_grid_cartesc_index::js, scale_atmos_grid_cartesc_index::ke, scale_atmos_grid_cartesc_index::ks, scale_tracer::qa, scale_tracer::tracer_advc, scale_atmos_grid_cartesc_index::xdir, scale_atmos_grid_cartesc_index::ydir, and scale_atmos_grid_cartesc_index::zdir.

Referenced by mod_atmos_phy_tb_driver::atmos_phy_tb_driver_calc_tendency().

172  use scale_precision
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
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_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
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
real(rp), dimension(:), allocatable, public atmos_grid_cartesc_fdz
z-length of grid(i+1) to grid(i) [m]
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 / grid / cartesC index
integer, public ke
end point of inner domain: z, local
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)
module atmosphere / grid / cartesC
integer, public ks
start point of inner domain: z, local
integer, public jblock
block size for cache blocking: y
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)
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
real(rp), dimension(:), allocatable, public atmos_grid_cartesc_rfdx
reciprocal of face-dx
Here is the call graph for this function:
Here is the caller graph for this function: