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
History
  • 2015-10-29 (S.Nishizawa) [new]
  • 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 82 of file scale_atmos_phy_tb_d1980.F90.

References scale_stdio::io_fid_log, scale_stdio::io_l, scale_process::prc_mpistop(), and scale_tracer::tracer_regist().

Referenced by scale_atmos_phy_tb::atmos_phy_tb_config().

82  use scale_process, only: &
84  use scale_tracer, only: &
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
subroutine, public prc_mpistop
Abort MPI.
module TRACER
module PROCESS
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 116 of file scale_atmos_phy_tb_d1980.F90.

References scale_grid_index::ia, scale_grid_index::ie, scale_stdio::io_fid_conf, scale_stdio::io_fid_log, scale_stdio::io_fid_nml, scale_stdio::io_l, scale_stdio::io_nml, scale_grid_index::is, scale_grid_index::ja, scale_grid_index::je, scale_grid_index::js, scale_grid_index::ka, scale_grid_index::ke, scale_grid_index::ks, and scale_process::prc_mpistop().

Referenced by scale_atmos_phy_tb::atmos_phy_tb_config().

116  use scale_process, only: &
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
integer, public is
start point of inner domain: x, local
integer, public je
end point of inner domain: y, local
subroutine, public prc_mpistop
Abort MPI.
integer, public js
start point of inner domain: y, local
module PROCESS
integer, public ie
end point of inner domain: x, local
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 173 of file scale_atmos_phy_tb_d1980.F90.

References 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_grid::grid_fdz, scale_grid::grid_rcdx, scale_grid::grid_rcdy, scale_grid::grid_rcdz, scale_grid::grid_rfdx, scale_grid::grid_rfdy, scale_grid::grid_rfdz, scale_gridtrans::i_uyw, scale_gridtrans::i_uyz, scale_gridtrans::i_xvw, scale_gridtrans::i_xvz, scale_gridtrans::i_xyw, scale_gridtrans::i_xyz, scale_grid_index::iblock, scale_grid_index::ie, scale_stdio::io_fid_log, scale_stdio::io_l, scale_grid_index::is, scale_grid_index::jblock, scale_grid_index::je, scale_grid_index::js, scale_grid_index::ke, scale_grid_index::ks, scale_tracer::qa, scale_tracer::tracer_advc, scale_grid_index::xdir, scale_grid_index::ydir, and scale_grid_index::zdir.

Referenced by scale_atmos_phy_tb::atmos_phy_tb_config().

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
integer, public is
start point of inner domain: x, local
integer, public i_xvz
subroutine, public atmos_phy_tb_calc_tend_momz(MOMZ_t_TB, QFLX_MOMZ, GSQRT, J13G, J23G, J33G, MAPF, IIS, IIE, JJS, JJE)
integer, public je
end point of inner domain: y, local
real(rp), dimension(:), allocatable, public grid_rcdy
reciprocal of center-dy
subroutine, public atmos_phy_tb_diffusion_solver(phi, a, b, c, d, KE_TB)
subroutine, public atmos_phy_tb_calc_flux_phi(qflx_phi, DENS, PHI, Kh, FACT, GSQRT, J13G, J23G, J33G, MAPF, a, b, c, dt, implicit, IIS, IIE, JJS, JJE)
integer, public iblock
block size for cache blocking: x
subroutine, public atmos_phy_tb_calc_tend_momy(MOMY_t_TB, QFLX_MOMY, GSQRT, J13G, J23G, J33G, MAPF, IIS, IIE, JJS, JJE)
real(rp), dimension(:), allocatable, public grid_rcdx
reciprocal of center-dx
integer, parameter, public zdir
integer, parameter, public ydir
integer, parameter, public xdir
logical, dimension(qa_max), public tracer_advc
real(rp), dimension(:), allocatable, public grid_rfdy
reciprocal of face-dy
real(rp), dimension(:), allocatable, public grid_rcdz
reciprocal of center-dz
integer, public i_xvw
module grid index
subroutine, public atmos_phy_tb_calc_tend_phi(phi_t_TB, QFLX_phi, GSQRT, J13G, J23G, J33G, MAPF, IIS, IIE, JJS, JJE)
module TRACER
subroutine, public atmos_phy_tb_calc_tend_momx(MOMX_t_TB, QFLX_MOMX, GSQRT, J13G, J23G, J33G, MAPF, IIS, IIE, JJS, JJE)
module GRIDTRANS
real(rp), dimension(:), allocatable, public grid_fdz
z-length of grid(k+1) to grid(k) [m]
integer, public jblock
block size for cache blocking: y
integer, public i_xyw
real(rp), public const_grav
standard acceleration of gravity [m/s2]
Definition: scale_const.F90:48
integer, public js
start point of inner domain: y, local
subroutine, public atmos_phy_tb_calc_strain_tensor(S33_C, S11_C, S22_C, S31_C, S12_C, S23_C, S12_Z, S23_X, S31_Y, S2, DENS, MOMZ, MOMX, MOMY, GSQRT, J13G, J23G, J33G, MAPF)
integer, public i_uyw
module CONSTANT
Definition: scale_const.F90:14
module ATMOSPHERE / Physics Turbulence
module GRID (cartesian)
integer, public ie
end point of inner domain: x, local
integer, public i_uyz
module PRECISION
real(rp), dimension(:), allocatable, public grid_rfdx
reciprocal of face-dx
real(rp), dimension(:), allocatable, public grid_rfdz
reciprocal of face-dz
integer, public i_xyz
Here is the call graph for this function:
Here is the caller graph for this function: