SCALE-RM
Functions/Subroutines
scale_atmos_dyn_tstep_short_fvm_hevi Module Reference

module Atmosphere / Dynamics RK More...

Functions/Subroutines

subroutine, public atmos_dyn_tstep_short_fvm_hevi_regist (ATMOS_DYN_TYPE, VA_out, VAR_NAME, VAR_DESC, VAR_UNIT)
 Register. More...
 
subroutine, public atmos_dyn_tstep_short_fvm_hevi_setup
 Setup. More...
 
subroutine, public atmos_dyn_tstep_short_fvm_hevi (DENS_RK, MOMZ_RK, MOMX_RK, MOMY_RK, RHOT_RK, PROG_RK, mflx_hi, tflx_hi, DENS0, MOMZ0, MOMX0, MOMY0, RHOT0, DENS, MOMZ, MOMX, MOMY, RHOT, DENS_t, MOMZ_t, MOMX_t, MOMY_t, RHOT_t, PROG0, PROG, DPRES0, RT2P, CORIOLI, num_diff, divdmp_coef, DDIV, FLAG_FCT_MOMENTUM, FLAG_FCT_T, FLAG_FCT_ALONG_STREAM, CDZ, FDZ, FDX, FDY, RCDZ, RCDX, RCDY, RFDZ, RFDX, RFDY, PHI, GSQRT, J13G, J23G, J33G, MAPF, REF_dens, REF_rhot, BND_W, BND_E, BND_S, BND_N, dtrk, dt)
 
subroutine solve_direct (C, F1, F2, F3)
 

Detailed Description

module Atmosphere / Dynamics RK

Description
HEVI FVM scheme for Atmospheric dynamical process
Author
Team SCALE
History
  • 2013-06-18 (S.Nishizawa) [new] newly impremented
  • 2014-04-04 (S.Nishizawa) [mod] support terrain-following coordinate
NAMELIST
  • No namelist group
History Output
namedescriptionunitvariable
DENS_t_advch tendency of density (horiz. advection) kg/m3/s advch_t
DENS_t_advcv tendency of density (vert. advection) kg/m3/s advcv_t
MOMX_t_advch tendency of momentum x (horiz. advection) kg/m2/s2 advch_t
MOMX_t_advcv tendency of momentum x (vert. advection) kg/m2/s2 advcv_t
MOMX_t_cf tendency of momentum x (coliolis force) kg/m2/s2 cf_t
MOMX_t_ddiv tendency of momentum x (divergence damping) kg/m2/s2 ddiv_t
MOMX_t_pg tendency of momentum x (pressure gradient) kg/m2/s2 pg_t
MOMY_t_advch tendency of momentum y (horiz. advection) kg/m2/s2 advch_t
MOMY_t_advcv tendency of momentum y (vert. advection) kg/m2/s2 advcv_t
MOMY_t_cf tendency of momentum y (coliolis force) kg/m2/s2 cf_t
MOMY_t_ddiv tendency of momentum y (divergence damping) kg/m2/s2 ddiv_t
MOMY_t_pg tendency of momentum y (pressure gradient) kg/m2/s2 pg_t
MOMZ_t_advch tendency of momentum z (horiz. advection) kg/m2/s2 advch_t
MOMZ_t_advcv tendency of momentum z (vert. advection) kg/m2/s2 advcv_t
MOMZ_t_ddiv tendency of momentum z (divergence damping) kg/m2/s2 ddiv_t
MOMZ_t_pg tendency of momentum z (pressure gradient) kg/m2/s2 pg_t
RHOT_t_advch tendency of rho*theta (horiz. advection) K kg/m3/s advch_t
RHOT_t_advcv tendency of rho*theta (vert. advection) K kg/m3/s advcv_t

Function/Subroutine Documentation

◆ atmos_dyn_tstep_short_fvm_hevi_regist()

subroutine, public scale_atmos_dyn_tstep_short_fvm_hevi::atmos_dyn_tstep_short_fvm_hevi_regist ( character(len=*), intent(in)  ATMOS_DYN_TYPE,
integer, intent(out)  VA_out,
character(len=h_short), dimension(:), intent(out)  VAR_NAME,
character(len=h_mid), dimension(:), intent(out)  VAR_DESC,
character(len=h_short), dimension(:), intent(out)  VAR_UNIT 
)

Register.

Parameters
[out]va_outnumber of prognostic variables
[out]var_namename of the variables
[out]var_descdesc. of the variables
[out]var_unitunit of the variables

Definition at line 97 of file scale_atmos_dyn_tstep_short_fvm_hevi.F90.

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

Referenced by scale_atmos_dyn_tstep_short::atmos_dyn_tstep_short_regist().

97  use scale_process, only: &
99  implicit none
100 
101  character(len=*), intent(in) :: atmos_dyn_type
102  integer, intent(out) :: va_out
103  character(len=H_SHORT), intent(out) :: var_name(:)
104  character(len=H_MID), intent(out) :: var_desc(:)
105  character(len=H_SHORT), intent(out) :: var_unit(:)
106  !---------------------------------------------------------------------------
107 
108  if( io_l ) write(io_fid_log,*) '*** HEVI Register'
109 
110  if ( atmos_dyn_type /= 'FVM-HEVI' .AND. atmos_dyn_type /= 'HEVI' ) then
111  write(*,*) 'xxx ATMOS_DYN_TYPE is not FVM-HEVI. Check!'
112  call prc_mpistop
113  endif
114 
115  va_out = va_fvm_hevi
116  var_name(:) = ""
117  var_desc(:) = ""
118  var_unit(:) = ""
119 
120  return
subroutine, public prc_mpistop
Abort MPI.
module PROCESS
Here is the call graph for this function:
Here is the caller graph for this function:

◆ atmos_dyn_tstep_short_fvm_hevi_setup()

subroutine, public scale_atmos_dyn_tstep_short_fvm_hevi::atmos_dyn_tstep_short_fvm_hevi_setup ( )

Setup.

Definition at line 126 of file scale_atmos_dyn_tstep_short_fvm_hevi.F90.

References scale_grid_index::ia, scale_stdio::io_fid_log, scale_stdio::io_l, scale_grid_index::ja, and scale_grid_index::ka.

Referenced by scale_atmos_dyn_tstep_short::atmos_dyn_tstep_short_regist().

126  implicit none
127  !---------------------------------------------------------------------------
128 
129  if( io_l ) write(io_fid_log,*) '*** HEVI Setup'
130 #ifdef HEVI_BICGSTAB
131  if ( io_l ) write(io_fid_log,*) '*** USING Bi-CGSTAB'
132 #elif defined(HEVI_LAPACK)
133  if ( io_l ) write(io_fid_log,*) '*** USING LAPACK'
134 #else
135  if ( io_l ) write(io_fid_log,*) '*** USING DIRECT'
136 #endif
137 
138 #ifdef HIST_TEND
139  allocate( advch_t(ka,ia,ja,5) )
140  allocate( advcv_t(ka,ia,ja,5) )
141  allocate( ddiv_t(ka,ia,ja,3) )
142  allocate( pg_t(ka,ia,ja,3) )
143  allocate( cf_t(ka,ia,ja,2) )
144 
145  advch_t = 0.0_rp
146  advcv_t = 0.0_rp
147  ddiv_t = 0.0_rp
148  pg_t = 0.0_rp
149  cf_t = 0.0_rp
150 #endif
151 
152  return
integer, public ia
of x whole cells (local, with HALO)
integer, public ka
of z whole cells (local, with HALO)
integer, public ja
of y whole cells (local, with HALO)
Here is the caller graph for this function:

◆ atmos_dyn_tstep_short_fvm_hevi()

subroutine, public scale_atmos_dyn_tstep_short_fvm_hevi::atmos_dyn_tstep_short_fvm_hevi ( real(rp), dimension(ka,ia,ja), intent(out)  DENS_RK,
real(rp), dimension(ka,ia,ja), intent(out)  MOMZ_RK,
real(rp), dimension(ka,ia,ja), intent(out)  MOMX_RK,
real(rp), dimension(ka,ia,ja), intent(out)  MOMY_RK,
real(rp), dimension(ka,ia,ja), intent(out)  RHOT_RK,
real(rp), dimension(ka,ia,ja,va), intent(out)  PROG_RK,
real(rp), dimension(ka,ia,ja,3), intent(inout)  mflx_hi,
real(rp), dimension(ka,ia,ja,3), intent(out)  tflx_hi,
real(rp), dimension(ka,ia,ja), intent(in), target  DENS0,
real(rp), dimension(ka,ia,ja), intent(in), target  MOMZ0,
real(rp), dimension(ka,ia,ja), intent(in), target  MOMX0,
real(rp), dimension(ka,ia,ja), intent(in), target  MOMY0,
real(rp), dimension(ka,ia,ja), intent(in), target  RHOT0,
real(rp), dimension(ka,ia,ja), intent(in)  DENS,
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_t,
real(rp), dimension(ka,ia,ja), intent(in)  MOMZ_t,
real(rp), dimension(ka,ia,ja), intent(in)  MOMX_t,
real(rp), dimension(ka,ia,ja), intent(in)  MOMY_t,
real(rp), dimension(ka,ia,ja), intent(in)  RHOT_t,
real(rp), dimension(ka,ia,ja,va), intent(in)  PROG0,
real(rp), dimension (ka,ia,ja,va), intent(in)  PROG,
real(rp), dimension(ka,ia,ja), intent(in)  DPRES0,
real(rp), dimension(ka,ia,ja), intent(in)  RT2P,
real(rp), dimension(1,ia,ja), intent(in)  CORIOLI,
real(rp), dimension(ka,ia,ja,5,3), intent(in)  num_diff,
real(rp), intent(in)  divdmp_coef,
real(rp), dimension(ka,ia,ja), intent(in)  DDIV,
logical, intent(in)  FLAG_FCT_MOMENTUM,
logical, intent(in)  FLAG_FCT_T,
logical, intent(in)  FLAG_FCT_ALONG_STREAM,
real(rp), dimension(ka), intent(in)  CDZ,
real(rp), dimension(ka-1), intent(in)  FDZ,
real(rp), dimension(ia-1), intent(in)  FDX,
real(rp), dimension(ja-1), intent(in)  FDY,
real(rp), dimension(ka), intent(in)  RCDZ,
real(rp), dimension(ia), intent(in)  RCDX,
real(rp), dimension(ja), intent(in)  RCDY,
real(rp), dimension(ka-1), intent(in)  RFDZ,
real(rp), dimension(ia-1), intent(in)  RFDX,
real(rp), dimension(ja-1), intent(in)  RFDY,
real(rp), dimension (ka,ia,ja), intent(in)  PHI,
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(rp), dimension(ka,ia,ja), intent(in)  REF_dens,
real(rp), dimension(ka,ia,ja), intent(in)  REF_rhot,
logical, intent(in)  BND_W,
logical, intent(in)  BND_E,
logical, intent(in)  BND_S,
logical, intent(in)  BND_N,
real(rp), intent(in)  dtrk,
real(rp), intent(in)  dt 
)
Parameters
[in]phigeopotential
[in]gsqrtvertical metrics {G}^1/2
[in]j13g(1,3) element of Jacobian matrix
[in]j23g(2,3) element of Jacobian matrix
[in]j33g(3,3) element of Jacobian matrix
[in]mapfmap factor
[in]ref_densreference density

Definition at line 174 of file scale_atmos_dyn_tstep_short_fvm_hevi.F90.

References scale_atmos_dyn_common::atmos_dyn_fct(), scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_flux_valuew_z, scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxj13_uyz, scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxj13_xvz, scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxj13_xyw, scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxj23_uyz, scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxj23_xvz, scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxj23_xyw, scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxx_uyz, scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxx_xvz, scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxx_xyw, scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxx_xyz, scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxy_uyz, scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxy_xvz, scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxy_xyw, scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxy_xyz, scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxz_uyz, scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxz_xvz, scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxz_xyw, scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxz_xyz, scale_debug::check(), dc_types::dp, scale_index::i_dens, scale_index::i_momx, scale_index::i_momy, scale_index::i_momz, scale_index::i_rhot, scale_gridtrans::i_uv, scale_gridtrans::i_uvz, scale_gridtrans::i_uy, scale_gridtrans::i_uyw, scale_gridtrans::i_uyz, scale_gridtrans::i_xv, scale_gridtrans::i_xvw, scale_gridtrans::i_xvz, scale_gridtrans::i_xy, scale_gridtrans::i_xyw, scale_gridtrans::i_xyz, scale_grid_index::iblock, scale_grid_index::ie, scale_grid_index::ieh, scale_grid_index::ihalo, scale_grid_index::is, scale_grid_index::jblock, scale_grid_index::je, scale_grid_index::jeh, scale_grid_index::jhalo, scale_grid_index::js, scale_grid_index::ke, scale_grid_index::kmax, scale_grid_index::ks, scale_process::prc_mpistop(), scale_precision::rp, solve_direct(), scale_grid_index::xdir, scale_grid_index::ydir, and scale_grid_index::zdir.

Referenced by scale_atmos_dyn_tstep_short::atmos_dyn_tstep_short_regist().

174  use scale_grid_index
175  use scale_const, only: &
176 #ifdef dry
177  rdry => const_rdry, &
178  cvdry => const_cvdry, &
179  cpdry => const_cpdry, &
180 #endif
181  grav => const_grav, &
182  p00 => const_pre00
183  use scale_atmos_dyn_common, only: &
185  use scale_atmos_dyn_fvm_flux, only: &
205  use scale_gridtrans, only: &
206  i_xyz, &
207  i_xyw, &
208  i_uyw, &
209  i_xvw, &
210  i_uyz, &
211  i_xvz, &
212  i_uvz, &
213  i_xy , &
214  i_uy , &
215  i_xv , &
216  i_uv
217 #ifdef HIST_TEND
218  use scale_history, only: &
219  hist_in
220 #endif
221  implicit none
222 
223  real(RP), intent(out) :: dens_rk(ka,ia,ja) ! prognostic variables
224  real(RP), intent(out) :: momz_rk(ka,ia,ja) !
225  real(RP), intent(out) :: momx_rk(ka,ia,ja) !
226  real(RP), intent(out) :: momy_rk(ka,ia,ja) !
227  real(RP), intent(out) :: rhot_rk(ka,ia,ja) !
228 
229  real(RP), intent(out) :: prog_rk(ka,ia,ja,va) !
230 
231  real(RP), intent(inout) :: mflx_hi(ka,ia,ja,3) ! rho * vel(x,y,z)
232  real(RP), intent(out) :: tflx_hi(ka,ia,ja,3) ! rho * theta * vel(x,y,z)
233 
234  real(RP), intent(in),target :: dens0(ka,ia,ja) ! prognostic variables
235  real(RP), intent(in),target :: momz0(ka,ia,ja) ! at previous dynamical time step
236  real(RP), intent(in),target :: momx0(ka,ia,ja) !
237  real(RP), intent(in),target :: momy0(ka,ia,ja) !
238  real(RP), intent(in),target :: rhot0(ka,ia,ja) !
239 
240  real(RP), intent(in) :: dens(ka,ia,ja) ! prognostic variables
241  real(RP), intent(in) :: momz(ka,ia,ja) ! at previous RK step
242  real(RP), intent(in) :: momx(ka,ia,ja) !
243  real(RP), intent(in) :: momy(ka,ia,ja) !
244  real(RP), intent(in) :: rhot(ka,ia,ja) !
245 
246  real(RP), intent(in) :: dens_t(ka,ia,ja)
247  real(RP), intent(in) :: momz_t(ka,ia,ja)
248  real(RP), intent(in) :: momx_t(ka,ia,ja)
249  real(RP), intent(in) :: momy_t(ka,ia,ja)
250  real(RP), intent(in) :: rhot_t(ka,ia,ja)
251 
252  real(RP), intent(in) :: prog0(ka,ia,ja,va)
253  real(RP), intent(in) :: prog (ka,ia,ja,va)
254 
255  real(RP), intent(in) :: dpres0(ka,ia,ja)
256  real(RP), intent(in) :: rt2p(ka,ia,ja)
257  real(RP), intent(in) :: corioli(1,ia,ja)
258  real(RP), intent(in) :: num_diff(ka,ia,ja,5,3)
259  real(RP), intent(in) :: divdmp_coef
260  real(RP), intent(in) :: ddiv(ka,ia,ja)
261 
262  logical, intent(in) :: flag_fct_momentum
263  logical, intent(in) :: flag_fct_t
264  logical, intent(in) :: flag_fct_along_stream
265 
266  real(RP), intent(in) :: cdz(ka)
267  real(RP), intent(in) :: fdz(ka-1)
268  real(RP), intent(in) :: fdx(ia-1)
269  real(RP), intent(in) :: fdy(ja-1)
270  real(RP), intent(in) :: rcdz(ka)
271  real(RP), intent(in) :: rcdx(ia)
272  real(RP), intent(in) :: rcdy(ja)
273  real(RP), intent(in) :: rfdz(ka-1)
274  real(RP), intent(in) :: rfdx(ia-1)
275  real(RP), intent(in) :: rfdy(ja-1)
276 
277  real(RP), intent(in) :: phi (ka,ia,ja)
278  real(RP), intent(in) :: gsqrt (ka,ia,ja,7)
279  real(RP), intent(in) :: j13g (ka,ia,ja,7)
280  real(RP), intent(in) :: j23g (ka,ia,ja,7)
281  real(RP), intent(in) :: j33g
282  real(RP), intent(in) :: mapf (ia,ja,2,4)
283  real(RP), intent(in) :: ref_dens(ka,ia,ja)
284  real(RP), intent(in) :: ref_rhot(ka,ia,ja)
285 
286  logical, intent(in) :: bnd_w
287  logical, intent(in) :: bnd_e
288  logical, intent(in) :: bnd_s
289  logical, intent(in) :: bnd_n
290 
291  real(RP), intent(in) :: dtrk
292  real(RP), intent(in) :: dt
293 
294 
295  ! diagnostic variables (work space)
296  real(RP) :: pott(ka,ia,ja) ! potential temperature [K]
297  real(RP) :: dpres(ka,ia,ja) ! pressure deviation from reference pressure
298 
299  real(RP) :: qflx_hi (ka,ia,ja,3)
300  real(RP) :: qflx_j13(ka,ia,ja)
301  real(RP) :: qflx_j23(ka,ia,ja)
302 
303  real(RP) :: advch ! horizontal advection
304  real(RP) :: advcv ! vertical advection
305  real(RP) :: div ! divergence damping
306  real(RP) :: pg ! pressure gradient force
307  real(RP) :: cf ! colioris force
308 #ifdef HIST_TEND
309  real(RP) :: advch_t(ka,ia,ja,5)
310  real(RP) :: advcv_t(ka,ia,ja,5)
311  real(RP) :: ddiv_t(ka,ia,ja,3)
312  real(RP) :: pg_t(ka,ia,ja,3)
313  real(RP) :: cf_t(ka,ia,ja,2)
314  logical :: lhist
315 #endif
316 
317  ! for implicit solver
318  real(RP) :: a(ka,ia,ja)
319  real(RP) :: b
320  real(RP) :: sr(ka,ia,ja)
321  real(RP) :: sw(ka,ia,ja)
322  real(RP) :: st(ka,ia,ja)
323  real(RP) :: pt(ka,ia,ja)
324  real(RP) :: c(kmax-1,ia,ja)
325 
326  real(RP) :: f1(ka,ia,ja)
327  real(RP) :: f2(ka,ia,ja)
328  real(RP) :: f3(ka,ia,ja)
329 
330  integer :: iis, iie, jjs, jje
331  integer :: k, i, j
332  integer :: iss, iee
333 
334 #ifdef DEBUG
335  pott(:,:,:) = undef
336  dpres(:,:,:) = undef
337 
338  pt(:,:,:) = undef
339 
340  qflx_hi(:,:,:,:) = undef
341  qflx_j13(:,:,:) = undef
342  qflx_j23(:,:,:) = undef
343 #endif
344 
345 #ifdef HIST_TEND
346  lhist = dt .eq. dtrk
347 #endif
348 
349 #ifdef PROFILE_FIPP
350  call fipp_start()
351 #endif
352 
353  ifs_off = 1
354  jfs_off = 1
355  if ( bnd_w ) ifs_off = 0
356  if ( bnd_s ) jfs_off = 0
357 
358 
359  do jjs = js, je, jblock
360  jje = jjs+jblock-1
361  do iis = is, ie, iblock
362  iie = iis+iblock-1
363 
364  profile_start("hevi_pres")
365  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
366  do j = jjs, jje+1
367  do i = iis, iie+1
368  do k = ks, ke
369 #ifdef DEBUG
370  call check( __line__, dpres0(k,i,j) )
371  call check( __line__, rt2p(k,i,j) )
372  call check( __line__, rhot(k,i,j) )
373  call check( __line__, ref_rhot(k,i,j) )
374 #endif
375  dpres(k,i,j) = dpres0(k,i,j) + rt2p(k,i,j) * ( rhot(k,i,j) - ref_rhot(k,i,j) )
376  enddo
377  dpres(ks-1,i,j) = dpres0(ks-1,i,j) - dens(ks,i,j) * ( phi(ks-1,i,j) - phi(ks+1,i,j) )
378  dpres(ke+1,i,j) = dpres0(ke+1,i,j) - dens(ke,i,j) * ( phi(ke+1,i,j) - phi(ke-1,i,j) )
379  enddo
380  enddo
381  profile_stop("hevi_pres")
382 
383  !##### continuity equation #####
384 
385  profile_start("hevi_mflx_z")
386  ! at (x, y, w)
387  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
388  do j = jjs, jje
389  do i = iis-1, iie
390  mflx_hi(ks-1,i,j,zdir) = 0.0_rp
391  do k = ks, ke-1
392 #ifdef DEBUG
393  call check( __line__, momx(k+1,i ,j) )
394  call check( __line__, momx(k+1,i-1,j) )
395  call check( __line__, momx(k ,i ,j) )
396  call check( __line__, momx(k ,i+1,j) )
397  call check( __line__, momy(k+1,i,j) )
398  call check( __line__, momy(k+1,i,j-1) )
399  call check( __line__, momy(k ,i,j) )
400  call check( __line__, momy(k ,i,j-1) )
401 #endif
402  mflx_hi(k,i,j,zdir) = j13g(k,i,j,i_xyw) * 0.25_rp / mapf(i,j,2,i_xy) &
403  * ( momx(k+1,i,j)+momx(k+1,i-1,j) &
404  + momx(k ,i,j)+momx(k ,i-1,j) ) &
405  + j23g(k,i,j,i_xyw) * 0.25_rp / mapf(i,j,1,i_xy) &
406  * ( momy(k+1,i,j)+momy(k+1,i,j-1) &
407  + momy(k ,i,j)+momy(k ,i,j-1) ) &
408  + gsqrt(k,i,j,i_xyw) / ( mapf(i,j,1,i_xy)*mapf(i,j,2,i_xy) ) * num_diff(k,i,j,i_dens,zdir)
409  enddo
410  mflx_hi(ke ,i,j,zdir) = 0.0_rp
411  enddo
412  enddo
413 #ifdef DEBUG
414  k = iundef; i = iundef; j = iundef
415 #endif
416  profile_stop("hevi_mflx_z")
417 
418  profile_start("hevi_mflx_x")
419  iss = max(iis-1,is-ifs_off)
420  iee = min(iie,ieh)
421  ! at (u, y, z)
422  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
423  do j = jjs, jje
424  do i = iss, iee
425  do k = ks, ke
426 #ifdef DEBUG
427  call check( __line__, gsqrt(k,i,j,i_uyz) )
428  call check( __line__, momx(k,i,j) )
429  call check( __line__, num_diff(k,i,j,i_dens,xdir) )
430 #endif
431  mflx_hi(k,i,j,xdir) = gsqrt(k,i,j,i_uyz) / mapf(i,j,2,i_uy) &
432  * ( momx(k,i,j) + num_diff(k,i,j,i_dens,xdir) )
433  enddo
434  enddo
435  enddo
436 #ifdef DEBUG
437  k = iundef; i = iundef; j = iundef
438 #endif
439  profile_stop("hevi_mflx_x")
440 
441  ! at (x, v, z)
442  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
443  do j = max(jjs-1,js-jfs_off), min(jje,jeh)
444  do i = iis, iie
445  do k = ks, ke
446 #ifdef DEBUG
447  call check( __line__, gsqrt(k,i,j,i_xvz) )
448  call check( __line__, momy(k,i,j) )
449  call check( __line__, num_diff(k,i,j,i_dens,ydir) )
450 #endif
451  mflx_hi(k,i,j,ydir) = gsqrt(k,i,j,i_xvz) / mapf(i,j,1,i_xv) &
452  * ( momy(k,i,j) + num_diff(k,i,j,i_dens,ydir) )
453  enddo
454  enddo
455  enddo
456 #ifdef DEBUG
457  k = iundef; i = iundef; j = iundef
458 #endif
459 
460  !--- update density
461  profile_start("hevi_sr")
462  !$omp parallel do private(i,j,k,advch) OMP_SCHEDULE_ collapse(2)
463  do j = jjs, jje
464  do i = iis, iie
465  do k = ks, ke
466 #ifdef DEBUG
467  call check( __line__, dens0(k,i,j) )
468  call check( __line__, mflx_hi(k ,i ,j ,xdir) )
469  call check( __line__, mflx_hi(k ,i-1,j ,xdir) )
470  call check( __line__, mflx_hi(k ,i ,j ,ydir) )
471  call check( __line__, mflx_hi(k ,i ,j-1,ydir) )
472  call check( __line__, dens_t(k,i,j) )
473 #endif
474  advch = - ( ( mflx_hi(k,i,j,zdir)-mflx_hi(k-1,i ,j, zdir) ) * rcdz(k) &
475  + ( mflx_hi(k,i,j,xdir)-mflx_hi(k ,i-1,j, xdir) ) * rcdx(i) &
476  + ( mflx_hi(k,i,j,ydir)-mflx_hi(k ,i, j-1,ydir) ) * rcdy(j) ) &
477  * mapf(i,j,1,i_xy) * mapf(i,j,2,i_xy) / gsqrt(k,i,j,i_xyz)
478  sr(k,i,j) = advch + dens_t(k,i,j)
479 #ifdef HIST_TEND
480  if ( lhist ) advch_t(k,i,j,i_dens) = advch
481 #endif
482  enddo
483  enddo
484  enddo
485 #ifdef DEBUG
486  k = iundef; i = iundef; j = iundef
487 #endif
488  profile_stop("hevi_sr")
489 
490  !##### momentum equation (z) #####
491 
492  ! at (x, y, z)
493  ! not that z-index is added by -1
494  profile_start("hevi_momz_qflxhi_z")
495  call atmos_dyn_fvm_fluxz_xyw( qflx_hi(:,:,:,zdir), & ! (out)
496  momz, momz, dens, & ! (in)
497  gsqrt(:,:,:,i_xyz), j33g, & ! (in)
498  num_diff(:,:,:,i_momz,zdir), & ! (in)
499  cdz, fdz, dtrk, &
500  iis, iie, jjs, jje ) ! (in)
501  profile_stop("hevi_momz_qflxhi_z")
502 
503  profile_start("hevi_momz_qflxj")
504  call atmos_dyn_fvm_fluxj13_xyw( qflx_j13, & ! (out)
505  momx, momz, dens, & ! (in)
506  gsqrt(:,:,:,i_xyz), j13g(:,:,:,i_xyz), mapf(:,:,:,i_xy), & ! (in)
507  cdz, &
508  iis, iie, jjs, jje ) ! (in)
509  call atmos_dyn_fvm_fluxj23_xyw( qflx_j23, & ! (out)
510  momy, momz, dens, & ! (in)
511  gsqrt(:,:,:,i_xyz), j23g(:,:,:,i_xyz), mapf(:,:,:,i_xy), & ! (in)
512  cdz, &
513  iis, iie, jjs, jje ) ! (in)
514  profile_stop("hevi_momz_qflxj")
515 
516  ! at (u, y, w)
517  profile_start("hevi_momz_qflxhi_x")
518  call atmos_dyn_fvm_fluxx_xyw( qflx_hi(:,:,:,xdir), & ! (out)
519  momx, momz, dens, & ! (in)
520  gsqrt(:,:,:,i_uyw), mapf(:,:,:,i_uy), & ! (in)
521  num_diff(:,:,:,i_momz,xdir), & ! (in)
522  cdz, & ! (in)
523  iis, iie, jjs, jje ) ! (in)
524  profile_stop("hevi_momz_qflxhi_x")
525 
526  ! at (x, v, w)
527  profile_start("hevi_momz_qflxhi_y")
528  call atmos_dyn_fvm_fluxy_xyw( qflx_hi(:,:,:,ydir), & ! (out)
529  momy, momz, dens, & ! (in)
530  gsqrt(:,:,:,i_xvw), mapf(:,:,:,i_xv), & ! (in)
531  num_diff(:,:,:,i_momz,ydir), & ! (in)
532  cdz, & ! (in)
533  iis, iie, jjs, jje ) ! (in)
534  profile_stop("hevi_momz_qflxhi_y")
535 
536  !--- update momentum(z)
537  profile_start("hevi_sw")
538  !$omp parallel do private(i,j,k,advcv,advch,cf,div) OMP_SCHEDULE_ collapse(2)
539  do j = jjs, jje
540  do i = iis, iie
541  do k = ks, ke-1
542 #ifdef DEBUG
543  call check( __line__, qflx_hi(k ,i ,j ,zdir) )
544  call check( __line__, qflx_hi(k-1,i ,j ,zdir) )
545  call check( __line__, qflx_j13(k ,i ,j) )
546  call check( __line__, qflx_j13(k-1,i ,j) )
547  call check( __line__, qflx_j23(k ,i ,j) )
548  call check( __line__, qflx_j23(k-1,i ,j) )
549  call check( __line__, qflx_hi(k ,i ,j ,xdir) )
550  call check( __line__, qflx_hi(k ,i-1,j ,xdir) )
551  call check( __line__, qflx_hi(k ,i ,j ,ydir) )
552  call check( __line__, qflx_hi(k ,i ,j-1,ydir) )
553  call check( __line__, ddiv(k ,i,j) )
554  call check( __line__, ddiv(k+1,i,j) )
555  call check( __line__, momz0(k,i,j) )
556  call check( __line__, momz_t(k,i,j) )
557 #endif
558  advcv = - ( qflx_hi(k,i,j,zdir) - qflx_hi(k-1,i ,j ,zdir) ) * rfdz(k)
559  advch = - ( ( qflx_j13(k,i,j) - qflx_j13(k-1,i ,j ) &
560  + qflx_j23(k,i,j) - qflx_j23(k-1,i ,j ) ) * rfdz(k) &
561  + ( qflx_hi(k,i,j,xdir) - qflx_hi(k,i-1,j ,xdir) ) * rcdx(i) &
562  + ( qflx_hi(k,i,j,ydir) - qflx_hi(k,i ,j-1,ydir) ) * rcdy(j) ) &
563  * mapf(i,j,1,i_xy) * mapf(i,j,2,i_xy)
564  div = divdmp_coef / dtrk * ( ddiv(k+1,i,j)-ddiv(k,i,j) ) * fdz(k) ! divergence damping
565  sw(k,i,j) = ( advcv + advch ) / gsqrt(k,i,j,i_xyw) &
566  + div + momz_t(k,i,j)
567 #ifdef HIST_TEND
568  if ( lhist ) then
569  advcv_t(k,i,j,i_momz) = advcv / gsqrt(k,i,j,i_xyw)
570  advch_t(k,i,j,i_momz) = advch / gsqrt(k,i,j,i_xyw)
571  ddiv_t(k,i,j,1) = div
572  endif
573 #endif
574  enddo
575  enddo
576  enddo
577  profile_stop("hevi_sw")
578 #ifdef DEBUG
579  k = iundef; i = iundef; j = iundef
580 #endif
581 
582 
583  !##### Thermodynamic Equation #####
584 
585  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
586  do j = jjs-jhalo, jje+jhalo
587  do i = iis-ihalo, iie+ihalo
588  do k = ks, ke
589 #ifdef DEBUG
590  call check( __line__, rhot(k,i,j) )
591  call check( __line__, dens(k,i,j) )
592 #endif
593  pott(k,i,j) = rhot(k,i,j) / dens(k,i,j)
594  enddo
595  enddo
596  enddo
597 #ifdef DEBUG
598  k = iundef; i = iundef; j = iundef
599 #endif
600 
601  ! at (x, y, w)
602  call atmos_dyn_fvm_fluxz_xyz( tflx_hi(:,:,:,zdir), & ! (out)
603  mflx_hi(:,:,:,zdir), pott, gsqrt(:,:,:,i_xyw), & ! (in)
604  num_diff(:,:,:,i_rhot,zdir), & ! (in)
605  cdz, & ! (in)
606  iis, iie, jjs, jje ) ! (in)
607 
608  ! at (u, y, z)
609  call atmos_dyn_fvm_fluxx_xyz( tflx_hi(:,:,:,xdir), & ! (out)
610  mflx_hi(:,:,:,xdir), pott, gsqrt(:,:,:,i_uyz), & ! (in)
611  num_diff(:,:,:,i_rhot,xdir), & ! (in)
612  cdz, & ! (in)
613  iis, iie, jjs, jje ) ! (in)
614 
615  ! at (x, v, z)
616  call atmos_dyn_fvm_fluxy_xyz( tflx_hi(:,:,:,ydir), & ! (out)
617  mflx_hi(:,:,:,ydir), pott, gsqrt(:,:,:,i_xvz), & ! (in)
618  num_diff(:,:,:,i_rhot,ydir), & ! (in)
619  cdz, & ! (in)
620  iis, iie, jjs, jje ) ! (in)
621 
622 
623  profile_start("hevi_st")
624  !$omp parallel do private(i,j,k,advch) OMP_SCHEDULE_ collapse(2)
625  do j = jjs, jje
626  do i = iis, iie
627  do k = ks, ke
628 #ifdef DEBUG
629  call check( __line__, tflx_hi(k ,i ,j ,zdir) )
630  call check( __line__, tflx_hi(k-1,i ,j ,zdir) )
631  call check( __line__, tflx_hi(k ,i ,j ,xdir) )
632  call check( __line__, tflx_hi(k ,i-1,j ,xdir) )
633  call check( __line__, tflx_hi(k ,i ,j ,ydir) )
634  call check( __line__, tflx_hi(k ,i ,j-1,ydir) )
635  call check( __line__, rhot_t(k,i,j) )
636 #endif
637  advch = - ( ( tflx_hi(k,i,j,zdir) - tflx_hi(k-1,i ,j ,zdir) ) * rcdz(k) &
638  + ( tflx_hi(k,i,j,xdir) - tflx_hi(k ,i-1,j ,xdir) ) * rcdx(i) &
639  + ( tflx_hi(k,i,j,ydir) - tflx_hi(k ,i ,j-1,ydir) ) * rcdy(j) ) &
640  * mapf(i,j,1,i_xy) * mapf(i,j,2,i_xy) / gsqrt(k,i,j,i_xyz)
641  st(k,i,j) = advch + rhot_t(k,i,j)
642 #ifdef HIST_TEND
643  if ( lhist ) then
644  advch_t(k,i,j,i_rhot) = advch
645  endif
646 #endif
647  enddo
648  enddo
649  enddo
650 #ifdef DEBUG
651  k = iundef; i = iundef; j = iundef
652 #endif
653  profile_stop("hevi_st")
654 
655  ! implicit solver
656 
657  profile_start("hevi_solver")
658 
659 !OCL INDEPENDENT
660 !OCL PREFETCH_SEQUENTIAL(SOFT)
661  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
662  do j = jjs, jje
663  do i = iis, iie
664 
665  call atmos_dyn_fvm_flux_valuew_z( pt(:,i,j), & ! (out)
666  momz(:,i,j), pott(:,i,j), gsqrt(:,i,j,i_xyz), & ! (in)
667  cdz )
668 
669  do k = ks, ke
670  a(k,i,j) = dtrk**2 * j33g * rcdz(k) * rt2p(k,i,j) * j33g / gsqrt(k,i,j,i_xyz)
671  enddo
672  b = grav * dtrk**2 * j33g / ( cdz(ks+1) + cdz(ks) )
673  f1(ks,i,j) = - ( pt(ks+1,i,j) * rfdz(ks) * a(ks+1,i,j) + b ) / gsqrt(ks,i,j,i_xyw)
674  f2(ks,i,j) = 1.0_rp + ( pt(ks ,i,j) * rfdz(ks) * ( a(ks+1,i,j)+a(ks,i,j) ) ) / gsqrt(ks,i,j,i_xyw)
675  do k = ks+1, ke-2
676  b = grav * dtrk**2 * j33g / ( cdz(k+1) + cdz(k) )
677  f1(k,i,j) = - ( pt(k+1,i,j) * rfdz(k) * a(k+1,i,j) + b ) / gsqrt(k,i,j,i_xyw)
678  f2(k,i,j) = 1.0_rp + ( pt(k ,i,j) * rfdz(k) * ( a(k+1,i,j)+a(k,i,j) ) ) / gsqrt(k,i,j,i_xyw)
679  f3(k,i,j) = - ( pt(k-1,i,j) * rfdz(k) * a(k,i,j) - b ) / gsqrt(k,i,j,i_xyw)
680  enddo
681  b = grav * dtrk**2 * j33g / ( cdz(ke) + cdz(ke-1) )
682  f2(ke-1,i,j) = 1.0_rp + ( pt(ke-1,i,j) * rfdz(ke-1) * ( a(ke,i,j)+a(ke-1,i,j) ) ) / gsqrt(ke-1,i,j,i_xyw)
683  f3(ke-1,i,j) = - ( pt(ke-2,i,j) * rfdz(ke-1) * a(ke-1,i,j) - b ) / gsqrt(ke-1,i,j,i_xyw)
684  do k = ks, ke-1
685  pg = - ( dpres(k+1,i,j) + rt2p(k+1,i,j)*dtrk*st(k+1,i,j) &
686  - dpres(k ,i,j) - rt2p(k ,i,j)*dtrk*st(k ,i,j) ) &
687  * rfdz(k) * j33g / gsqrt(k,i,j,i_xyw) &
688  - grav &
689  * ( f2h(k,1,i_xyz) * ( dens(k+1,i,j) - ref_dens(k+1,i,j) + sr(k+1,i,j) * dtrk ) &
690  + f2h(k,2,i_xyz) * ( dens(k ,i,j) - ref_dens(k ,i,j) + sr(k ,i,j) * dtrk ) )
691  c(k-ks+1,i,j) = momz(k,i,j) + dtrk * ( pg + sw(k,i,j) )
692 #ifdef HIST_TEND
693  if ( lhist ) pg_t(k,i,j,1) = pg
694 #endif
695  enddo
696 
697 #ifdef HEVI_BICGSTAB
698  call solve_bicgstab( &
699  c(:,i,j), & ! (inout)
700  f1(:,i,j), f2(:,i,j), f3(:,i,j) ) ! (in)
701 #elif defined(HEVI_LAPACK)
702  call solve_lapack( &
703  c(:,i,j), & ! (inout)
704 #ifdef DEBUG
705  i, j, & ! (in)
706 #endif
707  f1(:,i,j), f2(:,i,j), f3(:,i,j) ) ! (in)
708 #else
709  call solve_direct( &
710  c(:,i,j), & ! (inout)
711  f1(:,i,j), f2(:,i,j), f3(:,i,j) ) ! (in)
712 #endif
713 
714  do k = ks, ke-1
715 #ifdef DEBUG_HEVI2HEVE
716  ! for debug (change to explicit integration)
717  c(k-ks+1,i,j) = momz(k,i,j)
718  mflx_hi(k,i,j,zdir) = mflx_hi(k,i,j,zdir) &
719  + j33g * momz(k,i,j) / ( mapf(i,j,1,i_xy) * mapf(i,j,2,i_xy) )
720  momz_rk(k,i,j) = momz0(k,i,j) &
721  + dtrk*( &
722  - j33g * ( dpres(k+1,i,j)-dpres(k,i,j) ) * rfdz(k) / gsqrt(k,i,j,i_xyw) &
723  - grav * ( f2h(k,2,i_xyz)*(dens(k,i,j)-ref_dens(k,i,j))+f2h(k,1,i_xyz)*(dens(k+1,i,j)-ref_dens(k+1,i,j)) ) &
724  + sw(k,i,j) )
725 #else
726  ! z-momentum flux
727  mflx_hi(k,i,j,zdir) = mflx_hi(k,i,j,zdir) &
728  + j33g * c(k-ks+1,i,j) / ( mapf(i,j,1,i_xy) * mapf(i,j,2,i_xy) )
729  ! z-momentum
730  momz_rk(k,i,j) = momz0(k,i,j) &
731  + ( c(k-ks+1,i,j) - momz(k,i,j) )
732 #endif
733  enddo
734  momz_rk(ks-1,i,j) = 0.0_rp
735  momz_rk(ke ,i,j) = 0.0_rp
736 
737  ! density and rho*theta
738  advcv = - c(1,i,j) * j33g * rcdz(ks) / gsqrt(ks,i,j,i_xyz) ! C(0) = 0
739  dens_rk(ks,i,j) = dens0(ks,i,j) + dtrk * ( advcv + sr(ks,i,j) )
740 #ifdef HIST_TEND
741  if ( lhist ) advcv_t(ks,i,j,i_dens) = advcv
742 #endif
743  advcv = - c(1,i,j)*pt(ks,i,j) * j33g * rcdz(ks) / gsqrt(ks,i,j,i_xyz) ! C(0) = 0
744  rhot_rk(ks,i,j) = rhot0(ks,i,j) + dtrk * ( advcv + st(ks,i,j) )
745 #ifdef HIST_TEND
746  if ( lhist ) advcv_t(ks,i,j,i_rhot) = advcv
747 #endif
748  do k = ks+1, ke-1
749  advcv = - ( c(k-ks+1,i,j) - c(k-ks,i,j) ) &
750  * j33g * rcdz(k) / gsqrt(k,i,j,i_xyz)
751  dens_rk(k,i,j) = dens0(k,i,j) + dtrk * ( advcv + sr(k,i,j) )
752 #ifdef HIST_TEND
753  if ( lhist ) advcv_t(k,i,j,i_dens) = advcv
754 #endif
755  advcv = - ( c(k-ks+1,i,j)*pt(k,i,j) - c(k-ks,i,j)*pt(k-1,i,j) ) &
756  * j33g * rcdz(k) / gsqrt(k,i,j,i_xyz)
757  rhot_rk(k,i,j) = rhot0(k,i,j) + dtrk * ( advcv + st(k,i,j) )
758 #ifdef HIST_TEND
759  if ( lhist ) advcv_t(k,i,j,i_rhot) = advcv
760 #endif
761  enddo
762  advcv = c(ke-ks,i,j) * j33g * rcdz(ke) / gsqrt(ke,i,j,i_xyz) ! C(KE-KS+1) = 0
763  dens_rk(ke,i,j) = dens0(ke,i,j) + dtrk * ( advcv + sr(ke,i,j) )
764 #ifdef HIST_TEND
765  if ( lhist ) advcv_t(ke,i,j,i_dens) = advcv
766 #endif
767  advcv = c(ke-ks,i,j) * pt(ke-1,i,j) * j33g * rcdz(ke) / gsqrt(ke,i,j,i_xyz) ! C(KE-KS+1) = 0
768  rhot_rk(ke,i,j) = rhot0(ke,i,j) + dtrk * ( advcv + st(ke,i,j) )
769 #ifdef HIST_TEND
770  if ( lhist ) advcv_t(ke,i,j,i_rhot) = advcv
771 #endif
772 
773 #ifdef DEBUG
774  call check_equation( &
775  c(:,i,j), &
776  dens(:,i,j), momz(:,i,j), rhot(:,i,j), dpres(:,i,j), &
777  ref_dens(:,i,j), &
778  sr(:,i,j), sw(:,i,j), st(:,i,j), &
779  j33g, gsqrt(:,i,j,:), &
780  rt2p(:,i,j), &
781  dtrk, i, j )
782 #endif
783 
784  enddo
785  enddo
786 #ifdef DEBUG
787  k = iundef; i = iundef; j = iundef
788 #endif
789 
790  profile_stop("hevi_solver")
791 
792 
793  !##### momentum equation (x) #####
794 
795  profile_start("hevi_momx")
796  ! at (u, y, w)
797  call atmos_dyn_fvm_fluxz_uyz( qflx_hi(:,:,:,zdir), & ! (out)
798  momz, momx, dens, & ! (in)
799  gsqrt(:,:,:,i_uyw), j33g, & ! (in)
800  num_diff(:,:,:,i_momx,zdir), & ! (in)
801  cdz, & ! (in)
802  iis, iie, jjs, jje ) ! (in)
803  call atmos_dyn_fvm_fluxj13_uyz( qflx_j13, & ! (out)
804  momx, momx, dens, & ! (in)
805  gsqrt(:,:,:,i_uyz), j13g(:,:,:,i_uyw), mapf(:,:,:,i_uy), & ! (in)
806  cdz, & ! (in)
807  iis, iie, jjs, jje ) ! (in)
808  call atmos_dyn_fvm_fluxj23_uyz( qflx_j23, & ! (out)
809  momy, momx, dens, & ! (in)
810  gsqrt(:,:,:,i_uyz), j23g(:,:,:,i_uyw), mapf(:,:,:,i_uy), & ! (in)
811  cdz, & ! (in)
812  iis, iie, jjs, jje ) ! (in)
813 
814  ! at (x, y, z)
815  ! note that x-index is added by -1
816  call atmos_dyn_fvm_fluxx_uyz( qflx_hi(:,:,:,xdir), & ! (out)
817  momx, momx, dens, & ! (in)
818  gsqrt(:,:,:,i_xyz), mapf(:,:,:,i_xy), & ! (in)
819  num_diff(:,:,:,i_momx,xdir), & ! (in)
820  cdz, & ! (in)
821  iis, iie, jjs, jje ) ! (in)
822 
823  ! at (u, v, z)
824  call atmos_dyn_fvm_fluxy_uyz( qflx_hi(:,:,:,ydir), & ! (out)
825  momy, momx, dens, & ! (in)
826  gsqrt(:,:,:,i_uvz), mapf(:,:,1,i_uv), & ! (in)
827  num_diff(:,:,:,i_momx,ydir), & ! (in)
828  cdz, & ! (in)
829  iis, iie, jjs, jje ) ! (in)
830 
831  !--- update momentum(x)
832  iee = min(iie,ieh)
833  !$omp parallel do private(i,j,k,advch,advcv,pg,cf,div) OMP_SCHEDULE_ collapse(2)
834  do j = jjs, jje
835  do i = iis, iee
836  do k = ks, ke
837 #ifdef DEBUG
838  call check( __line__, qflx_hi(k ,i ,j ,zdir) )
839  call check( __line__, qflx_hi(k-1,i ,j ,zdir) )
840  call check( __line__, qflx_hi(k ,i ,j ,xdir) )
841  call check( __line__, qflx_hi(k ,i-1,j ,xdir) )
842  call check( __line__, qflx_hi(k ,i ,j ,ydir) )
843  call check( __line__, qflx_hi(k ,i ,j-1,ydir) )
844  call check( __line__, dpres(k,i+1,j) )
845  call check( __line__, dpres(k,i ,j) )
846  call check( __line__, corioli(1,i ,j) )
847  call check( __line__, corioli(1,i+1,j) )
848  call check( __line__, momy(k,i ,j ) )
849  call check( __line__, momy(k,i+1,j ) )
850  call check( __line__, momy(k,i ,j-1) )
851  call check( __line__, momy(k,i+1,j-1) )
852  call check( __line__, ddiv(k,i+1,j) )
853  call check( __line__, ddiv(k,i ,j) )
854  call check( __line__, momx0(k,i,j) )
855 #endif
856  advcv = - ( qflx_hi(k,i,j,zdir) - qflx_hi(k-1,i ,j ,zdir) ) * rcdz(k)
857  advch = - ( ( qflx_j13(k,i,j) - qflx_j13(k-1,i ,j ) &
858  + qflx_j23(k,i,j) - qflx_j23(k-1,i ,j ) ) * rfdz(k) &
859  + ( qflx_hi(k,i,j,xdir) - qflx_hi(k ,i-1,j ,xdir) ) * rfdx(i) &
860  + ( qflx_hi(k,i,j,ydir) - qflx_hi(k ,i ,j-1,ydir) ) * rcdy(j) ) &
861  * mapf(i,j,1,i_uy) * mapf(i,j,2,i_uy)
862  pg = ( ( gsqrt(k,i+1,j,i_xyz) * dpres(k,i+1,j) & ! [x,y,z]
863  - gsqrt(k,i ,j,i_xyz) * dpres(k,i ,j) & ! [x,y,z]
864  ) * rfdx(i) &
865  + ( j13g(k ,i,j,i_uyw) &
866  * 0.5_rp * ( f2h(k,1,i_uyz) * ( dpres(k+1,i+1,j)+dpres(k+1,i,j) ) &
867  + f2h(k,2,i_uyz) * ( dpres(k ,i+1,j)+dpres(k ,i,j) ) ) & ! [x,y,z->u,y,w]
868  - j13g(k-1,i,j,i_uyw) &
869  * 0.5_rp * ( f2h(k,1,i_uyz) * ( dpres(k ,i+1,j)+dpres(k ,i,j) ) &
870  + f2h(k,2,i_uyz) * ( dpres(k-1,i+1,j)+dpres(k-1,i,j) ) ) & ! [x,y,z->u,y,w]
871  ) * rcdz(k) ) &
872  * mapf(i,j,1,i_uy)
873  cf = 0.125_rp * ( corioli(1,i+1,j )+corioli(1,i,j ) ) & ! [x,y,z->u,y,z]
874  * ( momy(k,i+1,j )+momy(k,i,j ) &
875  + momy(k,i+1,j-1)+momy(k,i,j-1) ) & ! [x,v,z->u,y,z]
876  + 0.25_rp * mapf(i,j,1,i_uy) * mapf(i,j,2,i_uy) &
877  * ( momy(k,i,j) + momy(k,i,j-1) + momy(k,i+1,j) + momy(k,i+1,j-1) ) &
878  * ( ( momy(k,i,j) + momy(k,i,j-1) + momy(k,i+1,j) + momy(k,i+1,j-1) ) * 0.25_rp &
879  * ( 1.0_rp/mapf(i+1,j,2,i_xy) - 1.0_rp/mapf(i,j,2,i_xy) ) * rcdx(i) &
880  - momx(k,i,j) &
881  * ( 1.0_rp/mapf(i,j,1,i_uv) - 1.0_rp/mapf(i,j-1,1,i_uv) ) * rfdy(j) ) &
882  * 2.0_rp / ( dens(k,i+1,j) + dens(k,i,j) ) ! metric term
883  div = divdmp_coef / dtrk * ( ddiv(k,i+1,j)/mapf(i+1,j,2,i_xy) - ddiv(k,i,j)/mapf(i,j,1,i_xy) ) &
884  * mapf(i,j,1,i_uy) * mapf(i,j,2,i_uy) * fdx(i) ! divergence damping
885  momx_rk(k,i,j) = momx0(k,i,j) &
886  + dtrk * ( ( advcv + advch - pg ) / gsqrt(k,i,j,i_uyz) + cf + div + momx_t(k,i,j) )
887 #ifdef HIST_TEND
888  if ( lhist ) then
889  advcv_t(k,i,j,i_momx) = advcv / gsqrt(k,i,j,i_uyz)
890  advch_t(k,i,j,i_momx) = advch / gsqrt(k,i,j,i_uyz)
891  pg_t(k,i,j,2) = - pg / gsqrt(k,i,j,i_uyz)
892  cf_t(k,i,j,1) = cf
893  ddiv_t(k,i,j,2) = div
894  endif
895 #endif
896  enddo
897  enddo
898  enddo
899  profile_stop("hevi_momx")
900 #ifdef DEBUG
901  k = iundef; i = iundef; j = iundef
902 #endif
903 
904  !##### momentum equation (y) #####
905  profile_start("hevi_momy")
906  ! at (x, v, w)
907  call atmos_dyn_fvm_fluxz_xvz( qflx_hi(:,:,:,zdir), & ! (out)
908  momz, momy, dens, & ! (in)
909  gsqrt(:,:,:,i_xvw), j33g, & ! (in)
910  num_diff(:,:,:,i_momy,zdir), & ! (in)
911  cdz, & ! (in)
912  iis, iie, jjs, jje ) ! (in)
913  call atmos_dyn_fvm_fluxj13_xvz( qflx_j13, & ! (out)
914  momx, momy, dens, & ! (in)
915  gsqrt(:,:,:,i_xvz), j13g(:,:,:,i_xvw), mapf(:,:,:,i_xv), & ! (in)
916  cdz, & ! (in)
917  iis, iie, jjs, jje ) ! (in)
918  call atmos_dyn_fvm_fluxj23_xvz( qflx_j23, & ! (out)
919  momy, momy, dens, & ! (in)
920  gsqrt(:,:,:,i_xvz), j23g(:,:,:,i_xvw), mapf(:,:,:,i_xv), & ! (in)
921  cdz, & ! (in)
922  iis, iie, jjs, jje ) ! (in)
923 
924  ! at (u, v, z)
925  call atmos_dyn_fvm_fluxx_xvz( qflx_hi(:,:,:,xdir), & ! (out)
926  momx, momy, dens, & ! (in)
927  gsqrt(:,:,:,i_uvz), mapf(:,:,:,i_uv), & ! (in)
928  num_diff(:,:,:,i_momy,xdir), & ! (in)
929  cdz, & ! (in)
930  iis, iie, jjs, jje ) ! (in)
931 
932  ! at (x, y, z)
933  ! note that y-index is added by -1
934  call atmos_dyn_fvm_fluxy_xvz( qflx_hi(:,:,:,ydir), & ! (out)
935  momy, momy, dens, & ! (in)
936  gsqrt(:,:,:,i_xyz), mapf(:,:,:,i_xy), & ! (in)
937  num_diff(:,:,:,i_momy,ydir), & ! (in
938  cdz, & ! (in)
939  iis, iie, jjs, jje ) ! (in)
940 
941  !--- update momentum(y)
942  !$omp parallel do private(i,j,k,advch,advcv,pg,cf,div) OMP_SCHEDULE_ collapse(2)
943  do j = jjs, min(jje,jeh)
944  do i = iis, iie
945  do k = ks, ke
946 #ifdef DEBUG
947  call check( __line__, qflx_hi(k ,i ,j ,zdir) )
948  call check( __line__, qflx_hi(k-1,i ,j ,zdir) )
949  call check( __line__, qflx_hi(k ,i ,j ,xdir) )
950  call check( __line__, qflx_hi(k ,i-1,j ,xdir) )
951  call check( __line__, qflx_hi(k ,i ,j ,ydir) )
952  call check( __line__, qflx_hi(k ,i ,j-1,ydir) )
953  call check( __line__, dpres(k,i,j ) )
954  call check( __line__, dpres(k,i,j+1) )
955  call check( __line__, corioli(1,i,j ) )
956  call check( __line__, corioli(1,i,j+1) )
957  call check( __line__, momx(k,i ,j ) )
958  call check( __line__, momx(k,i ,j+1) )
959  call check( __line__, momx(k,i-1,j ) )
960  call check( __line__, momx(k,i-1,j+1) )
961  call check( __line__, ddiv(k,i,j+1) )
962  call check( __line__, ddiv(k,i,j ) )
963  call check( __line__, momy_t(k,i,j) )
964  call check( __line__, momy0(k,i,j) )
965 #endif
966  advcv = - ( qflx_hi(k,i,j,zdir) - qflx_hi(k-1,i ,j ,zdir) ) * rcdz(k)
967  advch = - ( ( qflx_j13(k,i,j) - qflx_j13(k-1,i ,j ) &
968  + qflx_j23(k,i,j) - qflx_j23(k-1,i ,j ) ) * rcdz(k) &
969  + ( qflx_hi(k,i,j,xdir) - qflx_hi(k ,i-1,j ,xdir) ) * rcdx(i) &
970  + ( qflx_hi(k,i,j,ydir) - qflx_hi(k ,i ,j-1,ydir) ) * rfdy(j) ) &
971  * mapf(i,j,1,i_xv) * mapf(i,j,2,i_xv)
972  pg = ( ( gsqrt(k,i,j+1,i_xyz) * dpres(k,i,j+1) & ! [x,y,z]
973  - gsqrt(k,i,j ,i_xyz) * dpres(k,i,j ) & ! [x,y,z]
974  ) * rfdy(j) &
975  + ( j23g(k ,i,j,i_xvw) &
976  * 0.5_rp * ( f2h(k ,1,i_xvz) * ( dpres(k+1,i,j+1)+dpres(k+1,i,j) ) &
977  + f2h(k ,2,i_xvz) * ( dpres(k ,i,j+1)+dpres(k ,i,j) ) ) & ! [x,y,z->x,v,w]
978  - j23g(k-1,i,j,i_xvw) &
979  * 0.5_rp * ( f2h(k-1,1,i_xvz) * ( dpres(k ,i,j+1)+dpres(k ,i,j) ) &
980  + f2h(k-1,2,i_xvz) * ( dpres(k-1,i,j+1)+dpres(k-1,i,j) ) ) & ! [x,y,z->x,v,w]
981  ) * rcdz(k) ) &
982  * mapf(i,j,2,i_xv)
983  cf = - 0.125_rp * ( corioli(1,i ,j+1)+corioli(1,i ,j) ) & ! [x,y,z->x,v,z]
984  * ( momx(k,i ,j+1)+momx(k,i ,j) &
985  + momx(k,i-1,j+1)+momx(k,i-1,j) ) & ! [u,y,z->x,v,z]
986  - 0.25_rp * mapf(i,j,1,i_xv) * mapf(i,j,2,i_xv) &
987  * ( momx(k,i,j) + momx(k,i-1,j) + momx(k,i,j+1) + momx(k,i-1,j+1) ) &
988  * ( momy(k,i,j) &
989  * ( 1.0_rp/mapf(i,j,2,i_uv) - 1.0_rp/mapf(i-1,j,2,i_uv) ) * rcdx(i) &
990  - 0.25_rp * ( momx(k,i,j)+momx(k,i-1,j)+momx(k,i,j+1)+momx(k,i-1,j+1) ) &
991  * ( 1.0_rp/mapf(i,j+1,1,i_xy) - 1.0_rp/mapf(i,j,1,i_xy) ) * rfdy(j) ) &
992  * 2.0_rp / ( dens(k,i,j+1) + dens(k,i,j) ) ! metoric term
993  div = divdmp_coef / dtrk * ( ddiv(k,i,j+1)/mapf(i,j+1,1,i_xy) - ddiv(k,i,j)/mapf(i,j,1,i_xy) ) &
994  * mapf(i,j,1,i_xv) * mapf(i,j,2,i_xv) * fdy(j) ! divergence damping
995  momy_rk(k,i,j) = momy0(k,i,j) &
996  + dtrk * ( ( advcv + advch - pg ) / gsqrt(k,i,j,i_xvz) + cf + div + momy_t(k,i,j) )
997 #ifdef HIST_TEND
998  if ( lhist ) then
999  advcv_t(k,i,j,i_momy) = advcv / gsqrt(k,i,j,i_xvz)
1000  advch_t(k,i,j,i_momy) = advch / gsqrt(k,i,j,i_xvz)
1001  pg_t(k,i,j,3) = - pg / gsqrt(k,i,j,i_xvz)
1002  cf_t(k,i,j,2) = cf
1003  ddiv_t(k,i,j,3) = div
1004  endif
1005 #endif
1006  enddo
1007  enddo
1008  enddo
1009  profile_stop("hevi_momy")
1010 #ifdef DEBUG
1011  k = iundef; i = iundef; j = iundef
1012 #endif
1013 
1014  enddo
1015  enddo
1016 
1017 #ifdef PROFILE_FIPP
1018  call fipp_stop()
1019 #endif
1020 
1021 #ifdef HIST_TEND
1022  if ( lhist ) then
1023  call hist_in(advcv_t(:,:,:,i_dens), 'DENS_t_advcv', 'tendency of density (vert. advection)', 'kg/m3/s' )
1024  call hist_in(advcv_t(:,:,:,i_momz), 'MOMZ_t_advcv', 'tendency of momentum z (vert. advection)', 'kg/m2/s2', zdim='half')
1025  call hist_in(advcv_t(:,:,:,i_momx), 'MOMX_t_advcv', 'tendency of momentum x (vert. advection)', 'kg/m2/s2', xdim='half')
1026  call hist_in(advcv_t(:,:,:,i_momy), 'MOMY_t_advcv', 'tendency of momentum y (vert. advection)', 'kg/m2/s2', ydim='half')
1027  call hist_in(advcv_t(:,:,:,i_rhot), 'RHOT_t_advcv', 'tendency of rho*theta (vert. advection)', 'K kg/m3/s' )
1028 
1029  call hist_in(advch_t(:,:,:,i_dens), 'DENS_t_advch', 'tendency of density (horiz. advection)', 'kg/m3/s' )
1030  call hist_in(advch_t(:,:,:,i_momz), 'MOMZ_t_advch', 'tendency of momentum z (horiz. advection)', 'kg/m2/s2', zdim='half')
1031  call hist_in(advch_t(:,:,:,i_momx), 'MOMX_t_advch', 'tendency of momentum x (horiz. advection)', 'kg/m2/s2', xdim='half')
1032  call hist_in(advch_t(:,:,:,i_momy), 'MOMY_t_advch', 'tendency of momentum y (horiz. advection)', 'kg/m2/s2', ydim='half')
1033  call hist_in(advch_t(:,:,:,i_rhot), 'RHOT_t_advch', 'tendency of rho*theta (horiz. advection)', 'K kg/m3/s' )
1034 
1035  call hist_in(pg_t(:,:,:,1), 'MOMZ_t_pg', 'tendency of momentum z (pressure gradient)', 'kg/m2/s2', zdim='half')
1036  call hist_in(pg_t(:,:,:,2), 'MOMX_t_pg', 'tendency of momentum x (pressure gradient)', 'kg/m2/s2', xdim='half')
1037  call hist_in(pg_t(:,:,:,3), 'MOMY_t_pg', 'tendency of momentum y (pressure gradient)', 'kg/m2/s2', ydim='half')
1038 
1039  call hist_in(ddiv_t(:,:,:,1), 'MOMZ_t_ddiv', 'tendency of momentum z (divergence damping)', 'kg/m2/s2', zdim='half')
1040  call hist_in(ddiv_t(:,:,:,2), 'MOMX_t_ddiv', 'tendency of momentum x (divergence damping)', 'kg/m2/s2', xdim='half')
1041  call hist_in(ddiv_t(:,:,:,3), 'MOMY_t_ddiv', 'tendency of momentum y (divergence damping)', 'kg/m2/s2', ydim='half')
1042 
1043  call hist_in(cf_t(:,:,:,1), 'MOMX_t_cf', 'tendency of momentum x (coliolis force)', 'kg/m2/s2', xdim='half')
1044  call hist_in(cf_t(:,:,:,2), 'MOMY_t_cf', 'tendency of momentum y (coliolis force)', 'kg/m2/s2', ydim='half')
1045  endif
1046 #endif
1047 
integer, parameter, public i_rhot
Definition: scale_index.F90:35
integer, public is
start point of inner domain: x, local
integer, public i_xvz
integer, public je
end point of inner domain: y, local
integer, public va
Definition: scale_index.F90:38
procedure(flux_mom), pointer, public atmos_dyn_fvm_fluxx_xvz
procedure(flux_j), pointer, public atmos_dyn_fvm_fluxj23_uyz
integer, public iblock
block size for cache blocking: x
integer, parameter, public i_momx
Definition: scale_index.F90:33
integer, parameter, public zdir
integer, parameter, public i_momz
Definition: scale_index.F90:32
procedure(flux_j), pointer, public atmos_dyn_fvm_fluxj23_xyw
procedure(flux_phi), pointer, public atmos_dyn_fvm_fluxx_xyz
integer, parameter, public ydir
integer, public ke
end point of inner domain: z, local
procedure(flux_j), pointer, public atmos_dyn_fvm_fluxj13_xyw
integer, parameter, public xdir
integer, public i_xy
integer, parameter, public i_dens
Definition: scale_index.F90:31
integer, parameter, public i_momy
Definition: scale_index.F90:34
integer, public i_xvw
procedure(flux_mom), pointer, public atmos_dyn_fvm_fluxx_xyw
module grid index
procedure(flux_mom), pointer, public atmos_dyn_fvm_fluxy_xyw
integer, public ia
of x whole cells (local, with HALO)
procedure(flux_z), pointer, public atmos_dyn_fvm_fluxz_xvz
module GRIDTRANS
integer, public ka
of z whole cells (local, with HALO)
integer, public i_uy
integer, public jblock
block size for cache blocking: y
integer, public i_xyw
integer, public kmax
of computational cells: z
integer, public jhalo
of halo cells: y
integer, public js
start point of inner domain: y, local
module Atmosphere / Dynamics common
integer, public i_uyw
procedure(valuew), pointer, public atmos_dyn_fvm_flux_valuew_z
module CONSTANT
Definition: scale_const.F90:14
integer, public ks
start point of inner domain: z, local
procedure(flux_j), pointer, public atmos_dyn_fvm_fluxj13_xvz
integer, public jeh
end point of inner domain: y, local (half level)
integer, public ieh
end point of inner domain: x, local (half level)
integer, public i_xv
integer, public ie
end point of inner domain: x, local
procedure(flux_z), pointer, public atmos_dyn_fvm_fluxz_uyz
module scale_atmos_dyn_fvm_flux
integer, public i_uyz
subroutine, public atmos_dyn_fct(qflx_anti, phi_in, DENS0, DENS, qflx_hi, qflx_lo, mflx_hi, rdz, rdx, rdy, GSQRT, MAPF, dt, flag_vect)
Flux Correction Transport Limiter.
module HISTORY
procedure(flux_j), pointer, public atmos_dyn_fvm_fluxj13_uyz
procedure(flux_j), pointer, public atmos_dyn_fvm_fluxj23_xvz
integer, public i_uvz
integer, public i_xyz
integer, public i_uv
procedure(flux_phi), pointer, public atmos_dyn_fvm_fluxz_xyz
procedure(flux_mom), pointer, public atmos_dyn_fvm_fluxx_uyz
procedure(flux_phi), pointer, public atmos_dyn_fvm_fluxy_xyz
procedure(flux_mom), pointer, public atmos_dyn_fvm_fluxy_uyz
procedure(flux_wz), pointer, public atmos_dyn_fvm_fluxz_xyw
integer, public ihalo
of halo cells: x
procedure(flux_mom), pointer, public atmos_dyn_fvm_fluxy_xvz
integer, public ja
of y whole cells (local, with HALO)
Here is the call graph for this function:
Here is the caller graph for this function:

◆ solve_direct()

subroutine scale_atmos_dyn_tstep_short_fvm_hevi::solve_direct ( real(rp), dimension(kmax-1), intent(inout)  C,
real(rp), dimension(ka), intent(in)  F1,
real(rp), dimension(ka), intent(in)  F2,
real(rp), dimension(ka), intent(in)  F3 
)

Definition at line 1280 of file scale_atmos_dyn_tstep_short_fvm_hevi.F90.

References scale_const::const_eps, scale_const::const_grav, scale_grid::grid_rcdz, scale_grid::grid_rfdz, scale_gridtrans::i_xyw, scale_gridtrans::i_xyz, scale_grid_index::ke, scale_grid_index::kmax, scale_grid_index::ks, and scale_process::prc_mpistop().

Referenced by atmos_dyn_tstep_short_fvm_hevi().

1280 
1281  use scale_process, only: &
1282  prc_mpistop
1283  implicit none
1284  real(RP), intent(inout) :: c(kmax-1)
1285  real(RP), intent(in) :: f1(ka)
1286  real(RP), intent(in) :: f2(ka)
1287  real(RP), intent(in) :: f3(ka)
1288 
1289  real(RP) :: e(kmax-2)
1290  real(RP) :: f(kmax-2)
1291 
1292  real(RP) :: rdenom
1293 
1294  integer :: k
1295 
1296  rdenom = 1.0_rp / f2(ks)
1297  e(1) = - f1(ks) * rdenom
1298  f(1) = c(1) * rdenom
1299  do k = 2, kmax-2
1300  rdenom = 1.0_rp / ( f2(k+ks-1) + f3(k+ks-1) * e(k-1) )
1301  e(k) = - f1(k+ks-1) * rdenom
1302  f(k) = ( c(k) - f3(k+ks-1) * f(k-1) ) * rdenom
1303  enddo
1304 
1305  ! C = \rho w
1306  c(kmax-1) = ( c(kmax-1) - f3(ke-1) * f(kmax-2) ) &
1307  / ( f2(ke-1) + f3(ke-1) * e(kmax-2) ) ! C(KMAX-1) = f(KMAX-1)
1308  do k = kmax-2, 1, -1
1309  c(k) = e(k) * c(k+1) + f(k)
1310  enddo
1311 
1312  return
subroutine, public prc_mpistop
Abort MPI.
integer, public ke
end point of inner domain: z, local
integer, public ka
of z whole cells (local, with HALO)
integer, public kmax
of computational cells: z
module PROCESS
integer, public ks
start point of inner domain: z, local
Here is the call graph for this function:
Here is the caller graph for this function: