SCALE-RM
Functions/Subroutines
scale_atmos_dyn_tstep_short_fvm_heve Module Reference

module Atmosphere / Dynamics RK More...

Functions/Subroutines

subroutine, public atmos_dyn_tstep_short_fvm_heve_regist (ATMOS_DYN_TYPE, VA_out, VAR_NAME, VAR_DESC, VAR_UNIT)
 Register. More...
 
subroutine, public atmos_dyn_tstep_short_fvm_heve_setup
 Setup. More...
 
subroutine, public atmos_dyn_tstep_short_fvm_heve (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, wdamp_coef, 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, last)
 

Detailed Description

module Atmosphere / Dynamics RK

Description
HEVE FVM scheme for Atmospheric dynamical process
Author
Team SCALE
History
  • 2013-06-18 (S.Nishizawa) [new] splited from dynamical core
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
MOMZ_t_wdamp tendency of momentum z (Raileight damping) kg/m2/s2 wdmp_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_heve_regist()

subroutine, public scale_atmos_dyn_tstep_short_fvm_heve::atmos_dyn_tstep_short_fvm_heve_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 75 of file scale_atmos_dyn_tstep_short_fvm_heve.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().

75  use scale_process, only: &
77  implicit none
78 
79  character(len=*), intent(in) :: ATMOS_DYN_TYPE
80  integer, intent(out) :: VA_out
81  character(len=H_SHORT), intent(out) :: VAR_NAME(:)
82  character(len=H_MID), intent(out) :: VAR_DESC(:)
83  character(len=H_SHORT), intent(out) :: VAR_UNIT(:)
84  !---------------------------------------------------------------------------
85 
86  if ( atmos_dyn_type /= 'FVM-HEVE' .AND. atmos_dyn_type /= 'HEVE' ) then
87  write(*,*) 'xxx ATMOS_DYN_TYPE is not FVM-HEVE. Check!'
88  call prc_mpistop
89  endif
90 
91  va_out = va_fvm_heve
92  var_name(:) = ""
93  var_desc(:) = ""
94  var_unit(:) = ""
95 
96  if( io_l ) write(io_fid_log,*)
97  if( io_l ) write(io_fid_log,*) '*** Register additional prognostic variables (HEVE)'
98  if ( va_out < 1 ) then
99  if( io_l ) write(io_fid_log,*) '*** => nothing.'
100  endif
101 
102  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_heve_setup()

subroutine, public scale_atmos_dyn_tstep_short_fvm_heve::atmos_dyn_tstep_short_fvm_heve_setup ( )

Setup.

Definition at line 108 of file scale_atmos_dyn_tstep_short_fvm_heve.F90.

Referenced by scale_atmos_dyn_tstep_short::atmos_dyn_tstep_short_regist().

108  return
Here is the caller graph for this function:

◆ atmos_dyn_tstep_short_fvm_heve()

subroutine, public scale_atmos_dyn_tstep_short_fvm_heve::atmos_dyn_tstep_short_fvm_heve ( 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), dimension(ka), intent(in)  wdamp_coef,
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,
logical, intent(in)  last 
)
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 130 of file scale_atmos_dyn_tstep_short_fvm_heve.F90.

References scale_atmos_dyn_common::atmos_dyn_fct(), 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_ud1::atmos_dyn_fvm_fluxx_uyz_ud1(), scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxx_xvz, scale_atmos_dyn_fvm_flux_ud1::atmos_dyn_fvm_fluxx_xvz_ud1(), scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxx_xyw, scale_atmos_dyn_fvm_flux_ud1::atmos_dyn_fvm_fluxx_xyw_ud1(), scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxx_xyz, scale_atmos_dyn_fvm_flux_ud1::atmos_dyn_fvm_fluxx_xyz_ud1(), scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxy_uyz, scale_atmos_dyn_fvm_flux_ud1::atmos_dyn_fvm_fluxy_uyz_ud1(), scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxy_xvz, scale_atmos_dyn_fvm_flux_ud1::atmos_dyn_fvm_fluxy_xvz_ud1(), scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxy_xyw, scale_atmos_dyn_fvm_flux_ud1::atmos_dyn_fvm_fluxy_xyw_ud1(), scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxy_xyz, scale_atmos_dyn_fvm_flux_ud1::atmos_dyn_fvm_fluxy_xyz_ud1(), scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxz_uyz, scale_atmos_dyn_fvm_flux_ud1::atmos_dyn_fvm_fluxz_uyz_ud1(), scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxz_xvz, scale_atmos_dyn_fvm_flux_ud1::atmos_dyn_fvm_fluxz_xvz_ud1(), scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxz_xyw, scale_atmos_dyn_fvm_flux_ud1::atmos_dyn_fvm_fluxz_xyw_ud1(), scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxz_xyz, scale_atmos_dyn_fvm_flux_ud1::atmos_dyn_fvm_fluxz_xyz_ud1(), scale_debug::check(), scale_const::const_eps, scale_const::const_grav, scale_const::const_pre00, 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::ka, scale_grid_index::ke, scale_grid_index::ks, 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().

130  use scale_grid_index
131  use scale_const, only: &
132  eps => const_eps, &
133  grav => const_grav, &
134  p00 => const_pre00
135  use scale_comm, only: &
136  comm_vars8, &
137  comm_wait
138  use scale_atmos_dyn_common, only: &
140  use scale_atmos_dyn_fvm_flux_ud1, only: &
153  use scale_atmos_dyn_fvm_flux, only: &
172  use scale_gridtrans, only: &
173  i_xyz, &
174  i_xyw, &
175  i_uyw, &
176  i_xvw, &
177  i_uyz, &
178  i_xvz, &
179  i_uvz, &
180  i_xy, &
181  i_uy, &
182  i_xv, &
183  i_uv
184 #ifdef HIST_TEND
185  use scale_history, only: &
186  hist_in
187 #endif
188  implicit none
189 
190  real(RP), intent(out) :: DENS_RK (KA,IA,JA) ! prognostic variables
191  real(RP), intent(out) :: MOMZ_RK (KA,IA,JA) !
192  real(RP), intent(out) :: MOMX_RK (KA,IA,JA) !
193  real(RP), intent(out) :: MOMY_RK (KA,IA,JA) !
194  real(RP), intent(out) :: RHOT_RK (KA,IA,JA) !
195  real(RP), intent(out) :: PROG_RK (KA,IA,JA,VA) !
196 
197  real(RP), intent(inout) :: mflx_hi (KA,IA,JA,3) ! mass flux
198  real(RP), intent(out) :: tflx_hi (KA,IA,JA,3) ! internal energy flux
199 
200  real(RP), intent(in), target :: DENS0 (KA,IA,JA) ! prognostic variables at previous dynamical time step
201  real(RP), intent(in), target :: MOMZ0 (KA,IA,JA) !
202  real(RP), intent(in), target :: MOMX0 (KA,IA,JA) !
203  real(RP), intent(in), target :: MOMY0 (KA,IA,JA) !
204  real(RP), intent(in), target :: RHOT0 (KA,IA,JA) !
205  real(RP), intent(in) :: PROG0 (KA,IA,JA,VA)
206 
207  real(RP), intent(in) :: DENS (KA,IA,JA) ! prognostic variables at previous RK step
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) :: PROG (KA,IA,JA,VA)
213 
214  real(RP), intent(in) :: DENS_t (KA,IA,JA) ! tendency
215  real(RP), intent(in) :: MOMZ_t (KA,IA,JA) !
216  real(RP), intent(in) :: MOMX_t (KA,IA,JA) !
217  real(RP), intent(in) :: MOMY_t (KA,IA,JA) !
218  real(RP), intent(in) :: RHOT_t (KA,IA,JA) !
219 
220  real(RP), intent(in) :: DPRES0 (KA,IA,JA)
221  real(RP), intent(in) :: RT2P (KA,IA,JA)
222  real(RP), intent(in) :: CORIOLI (1, IA,JA)
223  real(RP), intent(in) :: num_diff(KA,IA,JA,5,3)
224  real(RP), intent(in) :: wdamp_coef(KA)
225  real(RP), intent(in) :: divdmp_coef
226  real(RP), intent(in) :: DDIV (KA,IA,JA)
227 
228  logical, intent(in) :: FLAG_FCT_MOMENTUM
229  logical, intent(in) :: FLAG_FCT_T
230  logical, intent(in) :: FLAG_FCT_ALONG_STREAM
231 
232  real(RP), intent(in) :: CDZ (KA)
233  real(RP), intent(in) :: FDZ (KA-1)
234  real(RP), intent(in) :: FDX (IA-1)
235  real(RP), intent(in) :: FDY (JA-1)
236  real(RP), intent(in) :: RCDZ(KA)
237  real(RP), intent(in) :: RCDX(IA)
238  real(RP), intent(in) :: RCDY(JA)
239  real(RP), intent(in) :: RFDZ(KA-1)
240  real(RP), intent(in) :: RFDX(IA-1)
241  real(RP), intent(in) :: RFDY(JA-1)
242 
243  real(RP), intent(in) :: PHI (KA,IA,JA)
244  real(RP), intent(in) :: GSQRT (KA,IA,JA,7)
245  real(RP), intent(in) :: J13G (KA,IA,JA,7)
246  real(RP), intent(in) :: J23G (KA,IA,JA,7)
247  real(RP), intent(in) :: J33G
248  real(RP), intent(in) :: MAPF (IA,JA,2,4)
249  real(RP), intent(in) :: REF_dens(KA,IA,JA)
250  real(RP), intent(in) :: REF_rhot(KA,IA,JA)
251 
252  logical, intent(in) :: BND_W
253  logical, intent(in) :: BND_E
254  logical, intent(in) :: BND_S
255  logical, intent(in) :: BND_N
256 
257  real(RP), intent(in) :: dtrk
258  logical, intent(in) :: last
259 
260  ! diagnostic variables
261  real(RP) :: VELZ (KA,IA,JA) ! velocity w [m/s]
262  real(RP) :: VELX (KA,IA,JA) ! velocity u [m/s]
263  real(RP) :: VELY (KA,IA,JA) ! velocity v [m/s]
264  real(RP) :: POTT (KA,IA,JA) ! potential temperature [K]
265  real(RP) :: DPRES(KA,IA,JA) ! pressure - reference pressure
266 
267  real(RP) :: qflx_J13(KA,IA,JA)
268  real(RP) :: qflx_J23(KA,IA,JA)
269  real(RP) :: pgf (KA,IA,JA) ! pressure gradient force
270  real(RP) :: buoy (KA,IA,JA) ! buoyancy force
271  real(RP) :: cor (KA,IA,JA) ! Coriolis force
272 
273  ! flux
274  real(RP) :: qflx_hi (KA,IA,JA,3)
275 #ifndef NO_FCT_DYN
276  real(RP) :: qflx_lo (KA,IA,JA,3)
277  real(RP) :: qflx_anti(KA,IA,JA,3)
278  real(RP) :: tflx_lo (KA,IA,JA,3)
279  real(RP) :: tflx_anti(KA,IA,JA,3)
280  real(RP) :: DENS0_uvw(KA,IA,JA)
281  real(RP) :: DENS_uvw (KA,IA,JA)
282 #endif
283  real(RP) :: advch ! horizontal advection
284  real(RP) :: advcv ! vertical advection
285  real(RP) :: wdamp ! rayleight damping for W
286  real(RP) :: div ! divergence damping
287 #ifdef HIST_TEND
288  real(RP) :: advch_t(KA,IA,JA,5)
289  real(RP) :: advcv_t(KA,IA,JA,5)
290  real(RP) :: wdmp_t(KA,IA,JA)
291  real(RP) :: ddiv_t(KA,IA,JA,3)
292  real(RP) :: pg_t(KA,IA,JA,3)
293  real(RP) :: cf_t(KA,IA,JA,2)
294  logical :: lhist
295 #endif
296 
297  integer :: IIS, IIE
298  integer :: JJS, JJE
299  integer :: IFS_OFF, JFS_OFF
300  integer :: k, i, j
301  !---------------------------------------------------------------------------
302 
303 #ifdef DEBUG
304  velz(:,:,:) = undef
305  velx(:,:,:) = undef
306  vely(:,:,:) = undef
307  pott(:,:,:) = undef
308 
309  dpres(:,:,:) = undef
310 
311  mflx_hi(:,:,:,:) = undef
312  tflx_hi(:,:,:,:) = undef
313  qflx_hi(:,:,:,:) = undef
314 
315 #ifndef NO_FCT_DYN
316  qflx_lo(:,:,:,:) = undef
317  qflx_anti(:,:,:,:) = undef
318  tflx_lo(:,:,:,:) = undef
319  tflx_anti(:,:,:,:) = undef
320 #endif
321 #endif
322 
323 #if defined DEBUG || defined QUICKDEBUG
324  dens_rk( 1:ks-1,:,:) = undef
325  dens_rk(ke+1:ka ,:,:) = undef
326  momz_rk( 1:ks-1,:,:) = undef
327  momz_rk(ke+1:ka ,:,:) = undef
328  momx_rk( 1:ks-1,:,:) = undef
329  momx_rk(ke+1:ka ,:,:) = undef
330  momy_rk( 1:ks-1,:,:) = undef
331  momy_rk(ke+1:ka ,:,:) = undef
332  rhot_rk( 1:ks-1,:,:) = undef
333  rhot_rk(ke+1:ka ,:,:) = undef
334  prog_rk( 1:ks-1,:,:,:) = undef
335  prog_rk(ke+1:ka ,:,:,:) = undef
336 #endif
337 
338 #ifdef HIST_TEND
339  advch_t = 0.0_rp
340  advcv_t = 0.0_rp
341  wdmp_t = 0.0_rp
342  ddiv_t = 0.0_rp
343  pg_t = 0.0_rp
344  cf_t = 0.0_rp
345 
346  lhist = last
347 #endif
348 
349  ifs_off = 1
350  jfs_off = 1
351  if ( bnd_w ) ifs_off = 0
352  if ( bnd_s ) jfs_off = 0
353 
354 
355  do jjs = js, je, jblock
356  jje = jjs+jblock-1
357  do iis = is, ie, iblock
358  iie = iis+iblock-1
359 
360  ! pressure, pott. temp.
361 
362  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
363  do j = jjs-2, jje+2
364  do i = iis-2, iie+2
365  do k = ks, ke
366 #ifdef DEBUG
367  call check( __line__, dpres0(k,i,j) )
368  call check( __line__, rt2p(k,i,j) )
369  call check( __line__, rhot(k,i,j) )
370  call check( __line__, ref_rhot(k,i,j) )
371 #endif
372  dpres(k,i,j) = dpres0(k,i,j) + rt2p(k,i,j) * ( rhot(k,i,j) - ref_rhot(k,i,j) )
373  enddo
374  dpres(ks-1,i,j) = dpres0(ks-1,i,j) - dens(ks,i,j) * ( phi(ks-1,i,j) - phi(ks+1,i,j) )
375  dpres(ke+1,i,j) = dpres0(ke+1,i,j) - dens(ke,i,j) * ( phi(ke+1,i,j) - phi(ke-1,i,j) )
376  enddo
377  enddo
378 #ifdef DEBUG
379  k = iundef; i = iundef; j = iundef
380 #endif
381 
382  do j = jjs-jhalo, jje+jhalo
383  do i = iis-ihalo, iie+ihalo
384  do k = ks, ke
385 #ifdef DEBUG
386  call check( __line__, rhot(k,i,j) )
387  call check( __line__, dens(k,i,j) )
388 #endif
389  pott(k,i,j) = rhot(k,i,j) / dens(k,i,j)
390  enddo
391  enddo
392  enddo
393 #ifdef DEBUG
394  k = iundef; i = iundef; j = iundef
395 #endif
396 
397 #ifndef NO_FCT_DYN
398  enddo
399  enddo
400 
401  if ( flag_fct_momentum ) then
402  ! momentum -> velocity
403  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
404  do j = js-1, je+2
405  do i = is-1, ie+2
406  do k = ks, ke-1
407 #ifdef DEBUG
408  call check( __line__, momz0(k,i,j) )
409  call check( __line__, dens0(k ,i,j) )
410  call check( __line__, dens0(k+1,i,j) )
411 #endif
412  velz(k,i,j) = 2.0_rp * momz0(k,i,j) / ( dens0(k+1,i,j)+dens0(k,i,j) )
413  enddo
414  enddo
415  enddo
416 #ifdef DEBUG
417  k = iundef; i = iundef; j = iundef
418 #endif
419  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
420  do j = js-1, je+2
421  do i = is-1, ie+2
422  velz(ke,i,j) = 0.0_rp
423  enddo
424  enddo
425 #ifdef DEBUG
426  k = iundef; i = iundef; j = iundef
427 #endif
428 
429  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
430  do j = js-1, je+2
431  do i = is-2, ie+1
432  do k = ks, ke
433 #ifdef DEBUG
434  call check( __line__, momx0(k,i,j) )
435  call check( __line__, dens0(k,i ,j) )
436  call check( __line__, dens0(k,i+1,j) )
437 #endif
438  velx(k,i,j) = 2.0_rp * momx0(k,i,j) / ( dens0(k,i+1,j)+dens0(k,i,j) )
439  enddo
440  enddo
441  enddo
442 #ifdef DEBUG
443  k = iundef; i = iundef; j = iundef
444 #endif
445 
446  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
447  do j = js-2, je+1
448  do i = is-1, ie+2
449  do k = ks, ke
450 #ifdef DEBUG
451  call check( __line__, momy0(k,i,j) )
452  call check( __line__, dens0(k,i,j ) )
453  call check( __line__, dens0(k,i,j+1) )
454 #endif
455  vely(k,i,j) = 2.0_rp * momy0(k,i,j) / ( dens0(k,i,j+1)+dens0(k,i,j) )
456  enddo
457  enddo
458  enddo
459 #ifdef DEBUG
460  k = iundef; i = iundef; j = iundef
461 #endif
462  call comm_vars8( velz(:,:,:), 4 )
463  call comm_vars8( velx(:,:,:), 5 )
464  call comm_vars8( vely(:,:,:), 6 )
465  endif
466 
467  do jjs = js, je, jblock
468  jje = jjs+jblock-1
469  do iis = is, ie, iblock
470  iie = iis+iblock-1
471 #endif
472 
473  !########################################################################
474  ! continuity equation (total rho)
475  !########################################################################
476 
477  !-----< high order flux >-----
478 
479  ! at (x, y, w)
480 
481  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
482  do j = jjs, jje
483  do i = iis, iie
484  do k = ks, ke-1
485 #ifdef DEBUG
486  call check( __line__, momz(k+1,i,j) )
487  call check( __line__, momz(k ,i,j) )
488  call check( __line__, momz(k-1,i,j) )
489  call check( __line__, num_diff(k,i,j,i_dens,zdir) )
490 #endif
491  mflx_hi(k,i,j,zdir) = j33g * momz(k,i,j) / ( mapf(i,j,1,i_xy)*mapf(i,j,2,i_xy) ) &
492  + j13g(k,i,j,i_xyw) * 0.25_rp * ( momx(k+1,i,j)+momx(k+1,i-1,j) &
493  + momx(k ,i,j)+momx(k ,i-1,j) ) &
494  / mapf(i,j,2,i_xy) & ! [{u,y,z->x,y,w}]
495  + j23g(k,i,j,i_xyw) * 0.25_rp * ( momy(k+1,i,j)+momy(k+1,i,j-1) &
496  + momy(k ,i,j)+momy(k ,i,j-1) ) &
497  / mapf(i,j,1,i_xy) & ! [{x,v,z->x,y,w}]
498  + gsqrt(k,i,j,i_xyw) * num_diff(k,i,j,i_dens,zdir) / ( mapf(i,j,1,i_xy)*mapf(i,j,2,i_xy) )
499  enddo
500  enddo
501  enddo
502 #ifdef DEBUG
503  k = iundef; i = iundef; j = iundef
504 #endif
505 
506  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
507  do j = jjs, jje
508  do i = iis, iie
509  mflx_hi(ks-1,i,j,zdir) = 0.0_rp
510  mflx_hi(ke ,i,j,zdir) = 0.0_rp
511  enddo
512  enddo
513 #ifdef DEBUG
514  k = iundef; i = iundef; j = iundef
515 #endif
516 
517  ! at (u, y, z)
518 
519  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
520  do j = jjs , jje
521  do i = iis-ifs_off, min(iie,ieh)
522  do k = ks, ke
523 #ifdef DEBUG
524  call check( __line__, momx(k,i+1,j) )
525  call check( __line__, momx(k,i ,j) )
526  call check( __line__, momx(k,i-1,j) )
527  call check( __line__, num_diff(k,i,j,i_dens,xdir) )
528 #endif
529  mflx_hi(k,i,j,xdir) = gsqrt(k,i,j,i_uyz) / mapf(i,j,2,i_uy) &
530  * ( momx(k,i,j) + num_diff(k,i,j,i_dens,xdir) )
531  enddo
532  enddo
533  enddo
534 #ifdef DEBUG
535  k = iundef; i = iundef; j = iundef
536 #endif
537 
538  ! at (x, v, z)
539 
540  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
541  do j = jjs-jfs_off, min(jje,jeh)
542  do i = iis , iie
543  do k = ks, ke
544 #ifdef DEBUG
545  call check( __line__, momy(k,i,j+1) )
546  call check( __line__, momy(k,i,j ) )
547  call check( __line__, momy(k,i,j-1) )
548  call check( __line__, num_diff(k,i,j,i_dens,ydir) )
549 #endif
550  mflx_hi(k,i,j,ydir) = gsqrt(k,i,j,i_xvz) / mapf(i,j,1,i_xv) &
551  * ( momy(k,i,j) + num_diff(k,i,j,i_dens,ydir) )
552  enddo
553  enddo
554  enddo
555 #ifdef DEBUG
556  k = iundef; i = iundef; j = iundef
557 #endif
558 
559  !-----< update density >-----
560 
561  !$omp parallel do private(i,j,k,advcv,advch) OMP_SCHEDULE_ collapse(2)
562  do j = jjs, jje
563  do i = iis, iie
564  do k = ks, ke
565 #ifdef DEBUG
566  call check( __line__, dens0(k,i,j) )
567  call check( __line__, mflx_hi(k ,i ,j ,zdir) )
568  call check( __line__, mflx_hi(k-1,i ,j ,zdir) )
569  call check( __line__, mflx_hi(k ,i ,j ,xdir) )
570  call check( __line__, mflx_hi(k ,i-1,j ,xdir) )
571  call check( __line__, mflx_hi(k ,i ,j ,ydir) )
572  call check( __line__, mflx_hi(k ,i ,j-1,ydir) )
573  call check( __line__, dens_t(k,i,j) )
574 #endif
575  advcv = - ( mflx_hi(k,i,j,zdir)-mflx_hi(k-1,i, j, zdir) ) * rcdz(k)
576  advch = - ( mflx_hi(k,i,j,xdir)-mflx_hi(k ,i-1,j, xdir) ) * rcdx(i) &
577  - ( mflx_hi(k,i,j,ydir)-mflx_hi(k ,i, j-1,ydir) ) * rcdy(j)
578  dens_rk(k,i,j) = dens0(k,i,j) &
579  + dtrk * ( ( advcv + advch ) * mapf(i,j,1,i_xy) * mapf(i,j,2,i_xy) / gsqrt(k,i,j,i_xyz) &
580  + dens_t(k,i,j) )
581 #ifdef HIST_TEND
582  if ( lhist ) then
583  advcv_t(k,i,j,i_dens) = advcv * mapf(i,j,1,i_xy) * mapf(i,j,2,i_xy) / gsqrt(k,i,j,i_xyz)
584  advch_t(k,i,j,i_dens) = advch * mapf(i,j,1,i_xy) * mapf(i,j,2,i_xy) / gsqrt(k,i,j,i_xyz)
585  endif
586 #endif
587  enddo
588  enddo
589  enddo
590 #ifdef DEBUG
591  k = iundef; i = iundef; j = iundef
592 #endif
593 
594 
595  !########################################################################
596  ! momentum equation (z)
597  !########################################################################
598 
599  !-----< high order flux >-----
600 
601  ! at (x, y, z)
602  ! note than z-index is added by -1
603  call atmos_dyn_fvm_fluxz_xyw( qflx_hi(:,:,:,zdir), & ! (out)
604  momz, momz, dens, & ! (in)
605  gsqrt(:,:,:,i_xyz), j33g, & ! (in)
606  num_diff(:,:,:,i_momz,zdir), & ! (in)
607  cdz, fdz, dtrk, &
608  iis, iie, jjs, jje ) ! (in)
609  call atmos_dyn_fvm_fluxj13_xyw( qflx_j13, & ! (out)
610  momx, momz, dens, & ! (in)
611  gsqrt(:,:,:,i_xyz), j13g(:,:,:,i_xyz), mapf(:,:,:,i_xy), & ! (in)
612  cdz, &
613  iis, iie, jjs, jje ) ! (in)
614  call atmos_dyn_fvm_fluxj23_xyw( qflx_j23, & ! (out)
615  momy, momz, dens, & ! (in)
616  gsqrt(:,:,:,i_xyz), j23g(:,:,:,i_xyz), mapf(:,:,:,i_xy), & ! (in)
617  cdz, &
618  iis, iie, jjs, jje ) ! (in)
619 
620  ! at (u, y, w)
621  call atmos_dyn_fvm_fluxx_xyw( qflx_hi(:,:,:,xdir), & ! (out)
622  momx, momz, dens, & ! (in)
623  gsqrt(:,:,:,i_uyw), mapf(:,:,:,i_uy), & ! (in)
624  num_diff(:,:,:,i_momz,xdir), & ! (in)
625  cdz, & ! (in)
626  iis, iie, jjs, jje ) ! (in)
627 
628  ! at (x, v, w)
629  call atmos_dyn_fvm_fluxy_xyw( qflx_hi(:,:,:,ydir), & ! (out)
630  momy, momz, dens, & ! (in)
631  gsqrt(:,:,:,i_xvw), mapf(:,:,:,i_xv), & ! (in)
632  num_diff(:,:,:,i_momz,ydir), & ! (in)
633  cdz, & ! (in)
634  iis, iie, jjs, jje ) ! (in)
635 
636  ! pressure gradient force at (x, y, w)
637 
638  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
639  do j = jjs, jje
640  do i = iis, iie
641  do k = ks, ke-1
642  pgf(k,i,j) = j33g * ( dpres(k+1,i,j)-dpres(k,i,j) ) * rfdz(k) ! [x,y,z]
643  enddo
644  enddo
645  enddo
646 
647  ! buoyancy force at (x, y, w)
648 
649  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
650  do j = jjs, jje
651  do i = iis, iie
652  do k = ks, ke-1
653  buoy(k,i,j) = grav * gsqrt(k,i,j,i_xyw) &
654  * ( f2h(k,1,i_xyz) * ( dens(k+1,i,j)-ref_dens(k+1,i,j) ) &
655  + f2h(k,2,i_xyz) * ( dens(k ,i,j)-ref_dens(k ,i,j) ) )
656  enddo
657  enddo
658  enddo
659 
660  !-----< update momentum (z) -----
661 
662  !$omp parallel do private(i,j,k,advcv,advch,wdamp,div) OMP_SCHEDULE_ collapse(2)
663  do j = jjs, jje
664  do i = iis, iie
665  do k = ks, ke-1
666 #ifdef DEBUG
667  call check( __line__, qflx_hi(k ,i ,j ,zdir) )
668  call check( __line__, qflx_hi(k-1,i ,j ,zdir) )
669  call check( __line__, qflx_hi(k ,i ,j ,xdir) )
670  call check( __line__, qflx_hi(k ,i-1,j ,xdir) )
671  call check( __line__, qflx_hi(k ,i ,j ,ydir) )
672  call check( __line__, qflx_hi(k ,i ,j-1,ydir) )
673  call check( __line__, ddiv(k ,i,j) )
674  call check( __line__, ddiv(k+1,i,j) )
675  call check( __line__, momz0(k,i,j) )
676  call check( __line__, momz_t(k,i,j) )
677 #endif
678  advcv = - ( qflx_hi(k,i,j,zdir) - qflx_hi(k-1,i ,j ,zdir) &
679  + qflx_j13(k,i,j) - qflx_j13(k-1,i,j) &
680  + qflx_j23(k,i,j) - qflx_j23(k-1,i,j) ) * rfdz(k)
681  advch = - ( ( qflx_hi(k,i,j,xdir) - qflx_hi(k,i-1,j,xdir) ) * rcdx(i) &
682  + ( qflx_hi(k,i,j,ydir) - qflx_hi(k,i,j-1,ydir) ) * rcdy(j) ) &
683  * mapf(i,j,1,i_xy) * mapf(i,j,2,i_xy)
684  wdamp = - wdamp_coef(k) * momz0(k,i,j)
685  div = divdmp_coef / dtrk * fdz(k) * ( ddiv(k+1,i,j)-ddiv(k,i,j) ) ! divergence damping
686  momz_rk(k,i,j) = momz0(k,i,j) &
687  + dtrk * ( ( advcv + advch &
688  - pgf(k,i,j) & ! pressure gradient force
689  - buoy(k,i,j) & ! buoyancy force
690  ) / gsqrt(k,i,j,i_xyw) &
691  + wdamp & ! Rayleigh damping
692  + div &
693  + momz_t(k,i,j) ) ! physics tendency
694 #ifdef HIST_TEND
695  if ( lhist ) then
696  advcv_t(k,i,j,i_momz) = advcv / gsqrt(k,i,j,i_xyw)
697  advch_t(k,i,j,i_momz) = advch / gsqrt(k,i,j,i_xyw)
698  pg_t(k,i,j,1) = ( - pgf(k,i,j) - buoy(k,i,j) ) / gsqrt(k,i,j,i_xyw)
699  wdmp_t(k,i,j) = wdamp
700  ddiv_t(k,i,j,1) = div
701  endif
702 #endif
703  enddo
704  enddo
705  enddo
706 #ifdef DEBUG
707  k = iundef; i = iundef; j = iundef
708 #endif
709 
710  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
711  do j = jjs, jje
712  do i = iis, iie
713  momz_rk(ks-1,i,j) = 0.0_rp
714  momz_rk(ke ,i,j) = 0.0_rp
715 #ifdef HIST_TEND
716  if ( lhist ) then
717  advcv_t(ke,i,j,i_momz) = 0.0_rp
718  advch_t(ke,i,j,i_momz) = 0.0_rp
719  pg_t(ke,i,j,1) = 0.0_rp
720  wdmp_t(ke,i,j) = 0.0_rp
721  ddiv_t(ke,i,j,1) = 0.0_rp
722  endif
723 #endif
724  enddo
725  enddo
726 #ifdef DEBUG
727  k = iundef; i = iundef; j = iundef
728 #endif
729 
730 #ifndef NO_FCT_DYN
731  if ( flag_fct_momentum ) then
732 
733  ! monotonic flux
734  ! at (x, y, layer)
735  ! note than z-index is added by -1
736  call atmos_dyn_fvm_fluxz_xyw_ud1( qflx_lo(:,:,:,zdir), & ! (out)
737  momz, momz0, dens, & ! (in)
738  gsqrt(:,:,:,i_xyz), j33g, & ! (in)
739  num_diff(:,:,:,i_momz,zdir), & ! (in)
740  cdz, fdz, dtrk, & ! (in)
741  iis-1, iie+1, jjs-1, jje+1 ) ! (in)
742 
743  call atmos_dyn_fvm_fluxx_xyw_ud1( qflx_lo(:,:,:,xdir), & ! (out)
744  momx, momz0, dens, & ! (in)
745  gsqrt(:,:,:,i_uyz), mapf(:,:,:,i_uy), & ! (in)
746  num_diff(:,:,:,i_momz,xdir), & ! (in)
747  cdz, & ! (in)
748  iis-1, iie+1, jjs-1, jje+1 ) ! (in)
749 
750  call atmos_dyn_fvm_fluxy_xyw_ud1( qflx_lo(:,:,:,ydir), & ! (out)
751  momy, momz0, dens, & ! (in)
752  gsqrt(:,:,:,i_xvz), mapf(:,:,:,i_xv), & ! (in)
753  num_diff(:,:,:,i_momz,ydir), & ! (in)
754  cdz, & ! (in)
755  iis-1, iie+1, jjs-1, jje+1 ) ! (in)
756 
757  endif
758  enddo
759  enddo
760 
761  if ( flag_fct_momentum ) then
762 
763  call comm_vars8( dens_rk, 1 )
764  call comm_wait ( dens_rk, 1, .false. )
765 
766  do j = js, je
767  do i = is, ie
768  do k = ks, ke
769  qflx_hi(k,i,j,zdir) = qflx_hi(k,i,j,zdir) / ( mapf(i,j,1,i_xy) * mapf(i,j,2,i_xy) ) &
770  + qflx_j13(k,i,j) + qflx_j23(k,i,j)
771  enddo
772  enddo
773  enddo
774 
775  do j = js-1, je+1
776  do i = is-1, ie+1
777  do k = ks, ke-1
778  dens0_uvw(k,i,j) = 0.5_rp * ( dens0(k,i,j) + dens0(k+1,i,j) )
779  dens_uvw(k,i,j) = 0.5_rp * ( dens_rk(k,i,j) + dens_rk(k+1,i,j) )
780  enddo
781  enddo
782  enddo
783 
784  do j = js-1, je+1
785  do i = is-1, ie+1
786  dens_uvw(ke,i,j) = dens_uvw(ke-1,i,j)
787  dens0_uvw(ke,i,j) = dens0_uvw(ke-1,i,j)
788  enddo
789  enddo
790 
791  call comm_wait ( velz(:,:,:), 4 )
792 
793  call atmos_dyn_fct( qflx_anti, & ! (out)
794  velz, dens0_uvw, dens_uvw, & ! (in)
795  qflx_hi, qflx_lo, & ! (in)
796  mflx_hi, & ! (in)
797  rfdz, rcdx, rcdy, & ! (in)
798  gsqrt(:,:,:,i_xyw), & ! (in)
799  mapf(:,:,:,i_xy), dtrk, & ! (in)
800  flag_fct_along_stream ) ! (in)
801 
802  do jjs = js, je, jblock
803  jje = jjs+jblock-1
804  do iis = is, ie, iblock
805  iie = iis+iblock-1
806 
807  !--- update momentum(z)
808  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
809  do j = jjs, jje
810  do i = iis, iie
811  do k = ks, ke-1
812  momz_rk(k,i,j) = momz_rk(k,i,j) &
813  + dtrk * ( ( qflx_anti(k,i,j,zdir) - qflx_anti(k-1,i ,j ,zdir) ) * rfdz(k) &
814  + ( ( qflx_anti(k,i,j,xdir) - qflx_anti(k ,i-1,j ,xdir) ) * rcdx(i) &
815  + ( qflx_anti(k,i,j,ydir) - qflx_anti(k ,i ,j-1,ydir) ) * rcdy(j) ) &
816  * mapf(i,j,1,i_xy) * mapf(i,j,2,i_xy) ) &
817  / gsqrt(k,i,j,i_xyw)
818  enddo
819  enddo
820  enddo
821 
822  enddo
823  enddo
824 
825  endif ! FLAG_FCT_MOMENTUM
826 
827 #ifdef DEBUG
828  qflx_hi(:,:,:,:) = undef
829 #endif
830 
831  do jjs = js, je, jblock
832  jje = jjs+jblock-1
833  do iis = is, ie, iblock
834  iie = iis+iblock-1
835 #endif
836 
837  !########################################################################
838  ! momentum equation (x)
839  !########################################################################
840 
841  !-----< high order flux >-----
842 
843  ! at (u, y, w)
844  call atmos_dyn_fvm_fluxz_uyz( qflx_hi(:,:,:,zdir), & ! (out)
845  momz, momx, dens, & ! (in)
846  gsqrt(:,:,:,i_uyw), j33g, & ! (in)
847  num_diff(:,:,:,i_momx,zdir), & ! (in)
848  cdz, & ! (in)
849  iis, iie, jjs, jje ) ! (in)
850  call atmos_dyn_fvm_fluxj13_uyz( qflx_j13, & ! (out)
851  momx, momx, dens, & ! (in)
852  gsqrt(:,:,:,i_uyz), j13g(:,:,:,i_uyw), mapf(:,:,:,i_uy), & ! (in)
853  cdz, & ! (in)
854  iis, iie, jjs, jje ) ! (in)
855  call atmos_dyn_fvm_fluxj23_uyz( qflx_j23, & ! (out)
856  momy, momx, dens, & ! (in)
857  gsqrt(:,:,:,i_uyz), j23g(:,:,:,i_uyw), mapf(:,:,:,i_uy), & ! (in)
858  cdz, & ! (in)
859  iis, iie, jjs, jje ) ! (in)
860 
861  ! at (x, y, z)
862  ! note that x-index is added by -1
863  call atmos_dyn_fvm_fluxx_uyz( qflx_hi(:,:,:,xdir), & ! (out)
864  momx, momx, dens, & ! (in)
865  gsqrt(:,:,:,i_xyz), mapf(:,:,:,i_xy), & ! (in)
866  num_diff(:,:,:,i_momx,xdir), & ! (in)
867  cdz, & ! (in)
868  iis, iie, jjs, jje ) ! (in)
869 
870  ! at (u, v, z)
871  call atmos_dyn_fvm_fluxy_uyz( qflx_hi(:,:,:,ydir), & ! (out)
872  momy, momx, dens, & ! (in)
873  gsqrt(:,:,:,i_uvz), mapf(:,:,1,i_uv), & ! (in)
874  num_diff(:,:,:,i_momx,ydir), & ! (in)
875  cdz, & ! (in)
876  iis, iie, jjs, jje ) ! (in)
877 
878  ! pressure gradient force at (u, y, z)
879  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
880  do j = jjs, jje
881  do i = iis, iie
882  do k = ks, ke
883  pgf(k,i,j) = ( ( gsqrt(k,i+1,j,i_xyz) * dpres(k,i+1,j) & ! [x,y,z]
884  - gsqrt(k,i ,j,i_xyz) * dpres(k,i ,j) & ! [x,y,z]
885  ) * rfdx(i) &
886  + ( j13g(k ,i,j,i_uyw) &
887  * 0.5_rp * ( f2h(k,1,i_uyz) * ( dpres(k+1,i+1,j)+dpres(k+1,i,j) ) &
888  + f2h(k,2,i_uyz) * ( dpres(k ,i+1,j)+dpres(k ,i,j) ) ) & ! [x,y,z->u,y,w]
889  - j13g(k-1,i,j,i_uyw) &
890  * 0.5_rp * ( f2h(k,1,i_uyz) * ( dpres(k ,i+1,j)+dpres(k ,i,j) ) &
891  + f2h(k,2,i_uyz) * ( dpres(k-1,i+1,j)+dpres(k-1,i,j) ) ) & ! [x,y,z->u,y,w]
892  ) * rcdz(k) ) &
893  * mapf(i,j,1,i_uy)
894  enddo
895  enddo
896  enddo
897 
898  ! coriolis force at (u, y, z)
899 
900  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
901  do j = jjs, jje
902  do i = iis, iie
903  do k = ks, ke
904 #ifdef DEBUG
905  call check( __line__, momy(k,i ,j ) )
906  call check( __line__, momy(k,i+1,j ) )
907  call check( __line__, momy(k,i ,j-1) )
908  call check( __line__, momy(k,i+1,j-1) )
909 #endif
910  cor(k,i,j) = 0.125_rp * ( corioli(1,i+1,j )+corioli(1,i,j ) ) & ! [x,y,z->u,y,z]
911  * ( momy(k,i+1,j )+momy(k,i,j ) &
912  + momy(k,i+1,j-1)+momy(k,i,j-1) ) & ! [x,v,z->u,y,z]
913  + 0.25_rp * mapf(i,j,1,i_uy) * mapf(i,j,2,i_uy) &
914  * ( momy(k,i,j) + momy(k,i,j-1) + momy(k,i+1,j) + momy(k,i+1,j-1) ) &
915  * ( ( momy(k,i,j) + momy(k,i,j-1) + momy(k,i+1,j) + momy(k,i+1,j-1) ) * 0.25_rp &
916  * ( 1.0_rp/mapf(i+1,j,2,i_xy) - 1.0_rp/mapf(i,j,2,i_xy) ) * rfdx(i) &
917  - momx(k,i,j) &
918  * ( 1.0_rp/mapf(i,j,1,i_uv) - 1.0_rp/mapf(i,j-1,1,i_uv) ) * rcdy(j) ) &
919  * 2.0_rp / ( dens(k,i+1,j) + dens(k,i,j) ) ! metric term
920  enddo
921  enddo
922  enddo
923 
924  !-----< update momentum (x) >-----
925 
926  !$omp parallel do private(i,j,k,advcv,advch,div) OMP_SCHEDULE_ collapse(2)
927  do j = jjs, jje
928  do i = iis, min(iie, ieh)
929  do k = ks, ke
930 #ifdef DEBUG
931  call check( __line__, qflx_hi(k ,i ,j ,zdir) )
932  call check( __line__, qflx_hi(k-1,i ,j ,zdir) )
933  call check( __line__, qflx_hi(k ,i ,j ,xdir) )
934  call check( __line__, qflx_hi(k ,i-1,j ,xdir) )
935  call check( __line__, qflx_hi(k ,i ,j ,ydir) )
936  call check( __line__, qflx_hi(k ,i ,j-1,ydir) )
937  call check( __line__, ddiv(k,i+1,j) )
938  call check( __line__, ddiv(k,i ,j) )
939  call check( __line__, momx0(k,i,j) )
940 #endif
941  ! advection
942  advcv = - ( qflx_hi(k,i,j,zdir) - qflx_hi(k-1,i ,j ,zdir) &
943  + qflx_j13(k,i,j) - qflx_j13(k-1,i,j) &
944  + qflx_j23(k,i,j) - qflx_j23(k-1,i,j) ) * rcdz(k)
945  advch = - ( ( qflx_hi(k,i,j,xdir) - qflx_hi(k ,i-1,j ,xdir) ) * rfdx(i) &
946  + ( qflx_hi(k,i,j,ydir) - qflx_hi(k ,i ,j-1,ydir) ) * rcdy(j) ) &
947  * mapf(i,j,1,i_uy) * mapf(i,j,2,i_uy)
948  div = divdmp_coef / dtrk * fdx(i) * ( ddiv(k,i+1,j)-ddiv(k,i,j) ) ! divergence damping
949  momx_rk(k,i,j) = momx0(k,i,j) &
950  + dtrk * ( ( advcv + advch & ! advection
951  - pgf(k,i,j) & ! pressure gradient force
952  ) / gsqrt(k,i,j,i_uyz) &
953  + cor(k,i,j) & ! coriolis force
954  + div & ! divergence damping
955  + momx_t(k,i,j) ) ! physics tendency
956 #ifdef HIST_TEND
957  if ( lhist ) then
958  advcv_t(k,i,j,i_momx) = advcv / gsqrt(k,i,j,i_uyz)
959  advch_t(k,i,j,i_momx) = advch / gsqrt(k,i,j,i_uyz)
960  pg_t(k,i,j,2) = - pgf(k,i,j) / gsqrt(k,i,j,i_uyz)
961  cf_t(k,i,j,1) = cor(k,i,j)
962  ddiv_t(k,i,j,2) = div
963  endif
964 #endif
965  enddo
966  enddo
967  enddo
968 #ifdef DEBUG
969  k = iundef; i = iundef; j = iundef
970 #endif
971 
972 #ifndef NO_FCT_DYN
973  if ( flag_fct_momentum ) then
974 
975  call atmos_dyn_fvm_fluxz_uyz_ud1( qflx_lo(:,:,:,zdir), & ! (out)
976  momz, momx0, dens, & ! (in)
977  gsqrt(:,:,:,i_uyw), j33g, & ! (in)
978  num_diff(:,:,:,i_momx,zdir), & ! (in)
979  cdz, & ! (in)
980  iis-1, iie+1, jjs-1, jje+1 ) ! (in)
981 
982  ! note that x-index is added by -1
983  call atmos_dyn_fvm_fluxx_uyz_ud1( qflx_lo(:,:,:,xdir), & ! (out)
984  momx, momx0, dens, & ! (in)
985  gsqrt(:,:,:,i_xyz), mapf(:,:,:,i_uy), & ! (in)
986  num_diff(:,:,:,i_momx,xdir), & ! (in)
987  cdz, & ! (in)
988  iis-1, iie+1, jjs-1, jje+1 ) ! (in)
989 
990  call atmos_dyn_fvm_fluxy_uyz_ud1( qflx_lo(:,:,:,ydir), & ! (out)
991  momy, momx0, dens, & ! (in)
992  gsqrt(:,:,:,i_uvz), mapf(:,:,:,i_xv), & ! (in)
993  num_diff(:,:,:,i_momx,ydir), & ! (in)
994  cdz, & ! (in)
995  iis-1, iie+1, jjs-1, jje+1 ) ! (in)
996 
997  endif
998  enddo
999  enddo
1000 
1001  if ( flag_fct_momentum ) then
1002 
1003  do j = js, je
1004  do i = is, ie
1005  do k = ks, ke
1006  qflx_hi(k,i,j,zdir) = qflx_hi(k,i,j,zdir) / ( mapf(i,j,1,i_uy) * mapf(i,j,2,i_uy) )&
1007  + qflx_j13(k,i,j) + qflx_j23(k,i,j)
1008  enddo
1009  enddo
1010  enddo
1011 
1012  do j = js-1, je+1
1013  do i = is-1, ie+1
1014  do k = ks, ke
1015  dens0_uvw(k,i,j) = 0.5_rp * ( dens0(k,i,j) + dens0(k,i+1,j) )
1016  dens_uvw(k,i,j) = 0.5_rp * ( dens_rk(k,i,j) + dens_rk(k,i+1,j) )
1017  enddo
1018  enddo
1019  enddo
1020 
1021  call comm_wait ( velx(:,:,:), 5 )
1022 
1023  call atmos_dyn_fct( qflx_anti, & ! (out)
1024  velx, dens0_uvw, dens_uvw, & ! (in)
1025  qflx_hi, qflx_lo, & ! (in)
1026  mflx_hi, & ! (in)
1027  rcdz, rfdx, rcdy, & ! (in)
1028  gsqrt(:,:,:,i_uyz), & ! (in)
1029  mapf(:,:,:,i_uy), dtrk, & ! (in)
1030  flag_fct_along_stream ) ! (in)
1031 
1032  do jjs = js, je, jblock
1033  jje = jjs+jblock-1
1034  do iis = is, ie, iblock
1035  iie = iis+iblock-1
1036 
1037  !--- update momentum(x)
1038  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
1039  do j = jjs, jje
1040  do i = iis, min(iie,ieh)
1041  do k = ks, ke
1042 #ifdef DEBUG
1043  call check( __line__, momx_rk(k,i,j) )
1044  call check( __line__, qflx_anti(k ,i ,j ,zdir) )
1045  call check( __line__, qflx_anti(k-1,i ,j ,zdir) )
1046  call check( __line__, qflx_anti(k ,i ,j ,xdir) )
1047  call check( __line__, qflx_anti(k ,i-1,j ,xdir) )
1048  call check( __line__, qflx_anti(k ,i ,j ,ydir) )
1049  call check( __line__, qflx_anti(k ,i ,j-1,ydir) )
1050 #endif
1051  momx_rk(k,i,j) = momx_rk(k,i,j) &
1052  + dtrk * ( ( ( qflx_anti(k,i,j,zdir) - qflx_anti(k-1,i ,j ,zdir) ) * rcdz(k) &
1053  + ( qflx_anti(k,i,j,xdir) - qflx_anti(k ,i-1,j ,xdir) ) * rfdx(i) &
1054  + ( qflx_anti(k,i,j,ydir) - qflx_anti(k ,i ,j-1,ydir) ) * rcdy(j) ) ) &
1055  * mapf(i,j,1,i_uy) * mapf(i,j,2,i_uy) &
1056  / gsqrt(k,i,j,i_uyz)
1057  enddo
1058  enddo
1059  enddo
1060 #ifdef DEBUG
1061  k = iundef; i = iundef; j = iundef
1062 #endif
1063 
1064  enddo
1065  enddo
1066 
1067 #ifdef DEBUG
1068  qflx_lo(:,:,:,:) = undef
1069  qflx_anti(:,:,:,:) = undef
1070 #endif
1071 
1072  endif ! FLAG_FCT_MOMENTUM
1073 
1074 #ifdef DEBUG
1075  qflx_hi(:,:,:,:) = undef
1076 #endif
1077 
1078  do jjs = js, je, jblock
1079  jje = jjs+jblock-1
1080  do iis = is, ie, iblock
1081  iie = iis+iblock-1
1082 #endif
1083 
1084  !########################################################################
1085  ! momentum equation (y)
1086  !########################################################################
1087 
1088  !-----< high order flux >-----
1089 
1090  ! at (x, v, w)
1091  call atmos_dyn_fvm_fluxz_xvz( qflx_hi(:,:,:,zdir), & ! (out)
1092  momz, momy, dens, & ! (in)
1093  gsqrt(:,:,:,i_xvw), j33g, & ! (in)
1094  num_diff(:,:,:,i_momy,zdir), & ! (in)
1095  cdz, & ! (in)
1096  iis, iie, jjs, jje ) ! (in)
1097  call atmos_dyn_fvm_fluxj13_xvz( qflx_j13, & ! (out)
1098  momx, momy, dens, & ! (in)
1099  gsqrt(:,:,:,i_xvz), j13g(:,:,:,i_xvw), mapf(:,:,:,i_xv), & ! (in)
1100  cdz, & ! (in)
1101  iis, iie, jjs, jje ) ! (in)
1102  call atmos_dyn_fvm_fluxj23_xvz( qflx_j23, & ! (out)
1103  momy, momy, dens, & ! (in)
1104  gsqrt(:,:,:,i_xvz), j23g(:,:,:,i_xvw), mapf(:,:,:,i_xv), & ! (in)
1105  cdz, & ! (in)
1106  iis, iie, jjs, jje ) ! (in)
1107 
1108  ! at (u, v, z)
1109  call atmos_dyn_fvm_fluxx_xvz( qflx_hi(:,:,:,xdir), & ! (out)
1110  momx, momy, dens, & ! (in)
1111  gsqrt(:,:,:,i_uvz), mapf(:,:,:,i_uv), & ! (in)
1112  num_diff(:,:,:,i_momy,xdir), & ! (in)
1113  cdz, & ! (in)
1114  iis, iie, jjs, jje ) ! (in)
1115 
1116  ! at (x, y, z)
1117  ! note that y-index is added by -1
1118  call atmos_dyn_fvm_fluxy_xvz( qflx_hi(:,:,:,ydir), & ! (out)
1119  momy, momy, dens, & ! (in)
1120  gsqrt(:,:,:,i_xyz), mapf(:,:,:,i_xy), & ! (in)
1121  num_diff(:,:,:,i_momy,ydir), & ! (in
1122  cdz, & ! (in)
1123  iis, iie, jjs, jje ) ! (in)
1124 
1125  ! pressure gradient force at (x, v, z)
1126 
1127  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
1128  do j = jjs, jje
1129  do i = iis, iie
1130  do k = ks, ke
1131  pgf(k,i,j) = ( ( gsqrt(k,i,j+1,i_xyz) * dpres(k,i,j+1) & ! [x,y,z]
1132  - gsqrt(k,i,j ,i_xyz) * dpres(k,i,j ) & ! [x,y,z]
1133  ) * rfdy(j) &
1134  + ( j23g(k ,i,j,i_xvw) &
1135  * 0.5_rp * ( f2h(k ,1,i_xvz) * ( dpres(k+1,i,j+1)+dpres(k+1,i,j) ) &
1136  + f2h(k ,2,i_xvz) * ( dpres(k ,i,j+1)+dpres(k ,i,j) ) ) & ! [x,y,z->x,v,w]
1137  - j23g(k-1,i,j,i_xvw) &
1138  * 0.5_rp * ( f2h(k-1,1,i_xvz) * ( dpres(k ,i,j+1)+dpres(k ,i,j) ) &
1139  + f2h(k-1,2,i_xvz) * ( dpres(k-1,i,j+1)+dpres(k-1,i,j) ) ) & ! [x,y,z->x,v,w]
1140  ) * rcdz(k) ) &
1141  * mapf(i,j,2,i_xv)
1142  enddo
1143  enddo
1144  enddo
1145 
1146  ! coriolis force at (x, v, z)
1147 
1148  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
1149  do j = jjs, jje
1150  do i = iis, iie
1151  do k = ks, ke
1152 #ifdef DEBUG
1153  call check( __line__, momx(k,i ,j ) )
1154  call check( __line__, momx(k,i ,j+1) )
1155  call check( __line__, momx(k,i-1,j ) )
1156  call check( __line__, momx(k,i-1,j+1) )
1157 #endif
1158  cor(k,i,j) = - 0.125_rp * ( corioli(1,i ,j+1)+corioli(1,i ,j) ) & ! [x,y,z->x,v,z]
1159  * ( momx(k,i ,j+1)+momx(k,i ,j) &
1160  + momx(k,i-1,j+1)+momx(k,i-1,j) ) & ! [u,y,z->x,v,z]
1161  - 0.25_rp * mapf(i,j,1,i_xv) * mapf(i,j,2,i_xv) &
1162  * ( momx(k,i,j) + momx(k,i-1,j) + momx(k,i,j+1) + momx(k,i-1,j+1) )&
1163  * ( momy(k,i,j) &
1164  * ( 1.0_rp/mapf(i,j,2,i_uv) - 1.0_rp/mapf(i-1,j,2,i_uv) ) * rcdx(i) &
1165  - 0.25_rp * ( momx(k,i,j)+momx(k,i-1,j)+momx(k,i,j+1)+momx(k,i-1,j+1) ) &
1166  * ( 1.0_rp/mapf(i,j+1,1,i_xy) - 1.0_rp/mapf(i,j,1,i_xy) ) * rfdy(j) ) &
1167  * 2.0_rp / ( dens(k,i,j) + dens(k,i,j+1) ) ! metoric term
1168  enddo
1169  enddo
1170  enddo
1171 
1172  !-----< update momentum (y) >-----
1173 
1174  !$omp parallel do private(i,j,k,advcv,advch,div) OMP_SCHEDULE_ collapse(2)
1175  do j = jjs, min(jje, jeh)
1176  do i = iis, iie
1177  do k = ks, ke
1178 #ifdef DEBUG
1179  call check( __line__, qflx_hi(k ,i ,j ,zdir) )
1180  call check( __line__, qflx_hi(k-1,i ,j ,zdir) )
1181  call check( __line__, qflx_hi(k ,i ,j ,xdir) )
1182  call check( __line__, qflx_hi(k ,i-1,j ,xdir) )
1183  call check( __line__, qflx_hi(k ,i ,j ,ydir) )
1184  call check( __line__, qflx_hi(k ,i ,j-1,ydir) )
1185  call check( __line__, ddiv(k,i,j+1) )
1186  call check( __line__, ddiv(k,i,j ) )
1187  call check( __line__, momy_t(k,i,j) )
1188  call check( __line__, momy0(k,i,j) )
1189 #endif
1190 
1191  advcv = - ( qflx_hi(k,i,j,zdir) - qflx_hi(k-1,i ,j ,zdir) &
1192  + qflx_j13(k,i,j) - qflx_j13(k-1,i,j) &
1193  + qflx_j23(k,i,j) - qflx_j23(k-1,i,j) ) * rcdz(k)
1194  advch = - ( ( qflx_hi(k,i,j,xdir) - qflx_hi(k ,i-1,j ,xdir) ) * rcdx(i) &
1195  + ( qflx_hi(k,i,j,ydir) - qflx_hi(k ,i ,j-1,ydir) ) * rfdy(j) ) &
1196  * mapf(i,j,1,i_xv) * mapf(i,j,2,i_xv)
1197  div = divdmp_coef / dtrk * fdy(j) * ( ddiv(k,i,j+1)-ddiv(k,i,j) )
1198  momy_rk(k,i,j) = momy0(k,i,j) &
1199  + dtrk * ( ( advcv + advch & ! advection
1200  - pgf(k,i,j) & ! pressure gradient force
1201  ) / gsqrt(k,i,j,i_xvz) &
1202  + cor(k,i,j) & ! coriolis force
1203  + div & ! divergence damping
1204  + momy_t(k,i,j) ) ! physics tendency
1205 #ifdef HIST_TEND
1206  if ( lhist ) then
1207  advcv_t(k,i,j,i_momy) = advcv / gsqrt(k,i,j,i_uyz)
1208  advch_t(k,i,j,i_momy) = advch / gsqrt(k,i,j,i_uyz)
1209  pg_t(k,i,j,3) = - pgf(k,i,j) / gsqrt(k,i,j,i_uyz)
1210  cf_t(k,i,j,2) = cor(k,i,j)
1211  ddiv_t(k,i,j,3) = div
1212  endif
1213 #endif
1214  enddo
1215  enddo
1216  enddo
1217 #ifdef DEBUG
1218  k = iundef; i = iundef; j = iundef
1219 #endif
1220 
1221 #ifndef NO_FCT_DYN
1222  if ( flag_fct_momentum ) then
1223 
1224  ! monotonic flux
1225  ! at (x, v, interface)
1226  call atmos_dyn_fvm_fluxz_xvz_ud1( qflx_lo(:,:,:,zdir), & ! (out)
1227  momz, momy0, dens, & ! (in)
1228  gsqrt(:,:,:,i_xvz), j33g, & ! (in)
1229  num_diff(:,:,:,i_momy,zdir), & ! (in)
1230  cdz, & ! (in)
1231  iis-1, iie+1, jjs-1, jje+1 ) ! (in)
1232 
1233  ! at (u, v, layer)
1234  call atmos_dyn_fvm_fluxx_xvz_ud1( qflx_lo(:,:,:,xdir), & ! (out)
1235  momx, momy0, dens, & ! (in)
1236  gsqrt(:,:,:,i_uvz), mapf(:,:,:,i_xy), & ! (in)
1237  num_diff(:,:,:,i_momy,xdir), & ! (in)
1238  cdz, & ! (in)
1239  iis-1, iie+1, jjs-1, jje+1 ) ! (in)
1240 
1241  ! at (x, y, layer)
1242  ! note that y-index is added by -1
1243  call atmos_dyn_fvm_fluxy_xvz_ud1( qflx_lo(:,:,:,ydir), & ! (out)
1244  momy, momy0, dens, & ! (in)
1245  gsqrt(:,:,:,i_xyz), mapf(:,:,:,i_xy), & ! (in)
1246  num_diff(:,:,:,i_momy,ydir), & ! (in)
1247  cdz, & ! (in)
1248  iis-1, iie+1, jjs-1, jje+1 ) ! (in)
1249 
1250  endif
1251  enddo
1252  enddo
1253 
1254  if ( flag_fct_momentum ) then
1255 
1256  do j = js, je
1257  do i = is, ie
1258  do k = ks, ke
1259  qflx_hi(k,i,j,zdir) = qflx_hi(k,i,j,zdir) / ( mapf(i,j,1,i_xv) * mapf(i,j,2,i_xv) ) &
1260  + qflx_j13(k,i,j) + qflx_j23(k,i,j)
1261  enddo
1262  enddo
1263  enddo
1264 
1265  do j = js-1, je+1
1266  do i = is-1, ie+1
1267  do k = ks, ke
1268  dens0_uvw(k,i,j) = 0.5_rp * ( dens0(k,i,j) + dens0(k,i,j+1) )
1269  dens_uvw(k,i,j) = 0.5_rp * ( dens_rk(k,i,j) + dens_rk(k,i,j+1) )
1270  enddo
1271  enddo
1272  enddo
1273 
1274  call comm_wait ( vely(:,:,:), 6 )
1275 
1276  call atmos_dyn_fct( qflx_anti, & ! (out)
1277  vely, dens0_uvw, dens_uvw, & ! (in)
1278  qflx_hi, qflx_lo, & ! (in)
1279  mflx_hi, & ! (in)
1280  rcdz, rcdx, rfdy, & ! (in)
1281  gsqrt(:,:,:,i_xvz), & ! (in)
1282  mapf(:,:,:,i_xv), dtrk, & ! (in)
1283  flag_fct_along_stream ) ! (in)
1284 
1285  do jjs = js, je, jblock
1286  jje = jjs+jblock-1
1287  do iis = is, ie, iblock
1288  iie = iis+iblock-1
1289 
1290  !--- update momentum(y)
1291  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
1292  do j = jjs, min(jje, jeh)
1293  do i = iis, iie
1294  do k = ks, ke
1295 #ifdef DEBUG
1296  call check( __line__, momy_rk(k,i,j) )
1297  call check( __line__, qflx_anti(k ,i ,j ,zdir) )
1298  call check( __line__, qflx_anti(k-1,i ,j ,zdir) )
1299  call check( __line__, qflx_anti(k ,i ,j ,xdir) )
1300  call check( __line__, qflx_anti(k ,i-1,j ,xdir) )
1301  call check( __line__, qflx_anti(k ,i ,j ,ydir) )
1302  call check( __line__, qflx_anti(k ,i ,j-1,ydir) )
1303 #endif
1304  momy_rk(k,i,j) = momy_rk(k,i,j) &
1305  + dtrk * ( ( ( qflx_anti(k,i,j,zdir) - qflx_anti(k-1,i ,j ,zdir) ) * rcdz(k) &
1306  + ( qflx_anti(k,i,j,xdir) - qflx_anti(k ,i-1,j ,xdir) ) * rcdx(i) &
1307  + ( qflx_anti(k,i,j,ydir) - qflx_anti(k ,i ,j-1,ydir) ) * rfdy(j) ) ) &
1308  * mapf(i,j,1,i_xv) * mapf(i,j,2,i_xv) &
1309  / gsqrt(k,i,j,i_xvz)
1310  enddo
1311  enddo
1312  enddo
1313 #ifdef DEBUG
1314  k = iundef; i = iundef; j = iundef
1315 #endif
1316 
1317  enddo
1318  enddo
1319 
1320 #ifdef DEBUG
1321  qflx_lo(:,:,:,:) = undef
1322  qflx_anti(:,:,:,:) = undef
1323 #endif
1324  endif ! FLAG_FCT_MOMENTUM
1325 
1326 #ifdef DEBUG
1327  qflx_hi(ks:,:,:,:) = undef
1328 #endif
1329 
1330  do jjs = js, je, jblock
1331  jje = jjs+jblock-1
1332  do iis = is, ie, iblock
1333  iie = iis+iblock-1
1334 #endif
1335 
1336  !########################################################################
1337  ! Thermodynamic equation
1338  !########################################################################
1339 
1340  !-----< high order flux >-----
1341 
1342  ! at (x, y, w)
1343  call atmos_dyn_fvm_fluxz_xyz( tflx_hi(:,:,:,zdir), & ! (out)
1344  mflx_hi(:,:,:,zdir), pott, gsqrt(:,:,:,i_xyw), & ! (in)
1345  num_diff(:,:,:,i_rhot,zdir), & ! (in)
1346  cdz, & ! (in)
1347  iis, iie, jjs, jje ) ! (in)
1348 
1349  ! at (u, y, z)
1350  call atmos_dyn_fvm_fluxx_xyz( tflx_hi(:,:,:,xdir), & ! (out)
1351  mflx_hi(:,:,:,xdir), pott, gsqrt(:,:,:,i_uyz), & ! (in)
1352  num_diff(:,:,:,i_rhot,xdir), & ! (in)
1353  cdz, & ! (in)
1354  iis, iie, jjs, jje ) ! (in)
1355 
1356  ! at (x, v, z)
1357  call atmos_dyn_fvm_fluxy_xyz( tflx_hi(:,:,:,ydir), & ! (out)
1358  mflx_hi(:,:,:,ydir), pott, gsqrt(:,:,:,i_xvz), & ! (in)
1359  num_diff(:,:,:,i_rhot,ydir), & ! (in)
1360  cdz, & ! (in)
1361  iis, iie, jjs, jje ) ! (in)
1362 
1363  !-----< update rho*theta >-----
1364 
1365  !$omp parallel do private(i,j,k,advcv,advch) OMP_SCHEDULE_ collapse(2)
1366  do j = jjs, jje
1367  do i = iis, iie
1368  do k = ks, ke
1369 #ifdef DEBUG
1370  call check( __line__, tflx_hi(k ,i ,j ,zdir) )
1371  call check( __line__, tflx_hi(k-1,i ,j ,zdir) )
1372  call check( __line__, tflx_hi(k ,i ,j ,xdir) )
1373  call check( __line__, tflx_hi(k ,i-1,j ,xdir) )
1374  call check( __line__, tflx_hi(k ,i ,j ,ydir) )
1375  call check( __line__, tflx_hi(k ,i ,j-1,ydir) )
1376  call check( __line__, rhot_t(k,i,j) )
1377  call check( __line__, rhot0(k,i,j) )
1378 #endif
1379  advcv = - ( tflx_hi(k,i,j,zdir) - tflx_hi(k-1,i ,j ,zdir) ) * rcdz(k)
1380  advch = - ( tflx_hi(k,i,j,xdir) - tflx_hi(k ,i-1,j ,xdir) ) * rcdx(i) &
1381  - ( tflx_hi(k,i,j,ydir) - tflx_hi(k ,i ,j-1,ydir) ) * rcdy(j)
1382  rhot_rk(k,i,j) = rhot0(k,i,j) &
1383  + dtrk * ( ( advcv + advch ) * mapf(i,j,1,i_xy) * mapf(i,j,2,i_xy) / gsqrt(k,i,j,i_xyz) &
1384  + rhot_t(k,i,j) )
1385 #ifdef HIST_TEND
1386  if ( lhist ) then
1387  advcv_t(k,i,j,i_rhot) = advcv * mapf(i,j,1,i_xy) * mapf(i,j,2,i_xy)/ gsqrt(k,i,j,i_xyz)
1388  advch_t(k,i,j,i_rhot) = advch * mapf(i,j,1,i_xy) * mapf(i,j,2,i_xy)/ gsqrt(k,i,j,i_xyz)
1389  endif
1390 #endif
1391  enddo
1392  enddo
1393  enddo
1394 #ifdef DEBUG
1395  k = iundef; i = iundef; j = iundef
1396 #endif
1397 
1398  enddo
1399  enddo
1400 
1401 #ifndef NO_FCT_DYN
1402 
1403  if ( flag_fct_t ) then
1404 
1405  call comm_vars8( mflx_hi(:,:,:,zdir), 1 )
1406  call comm_vars8( mflx_hi(:,:,:,xdir), 2 )
1407  call comm_vars8( mflx_hi(:,:,:,ydir), 3 )
1408  call comm_wait ( mflx_hi(:,:,:,zdir), 1, .false. )
1409  call comm_wait ( mflx_hi(:,:,:,xdir), 2, .false. )
1410  call comm_wait ( mflx_hi(:,:,:,ydir), 3, .false. )
1411 
1412  if ( .NOT. flag_fct_momentum ) then
1413  call comm_vars8( dens_rk, 1 )
1414  call comm_wait ( dens_rk, 1, .false. )
1415  endif
1416 
1417  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
1418  do j = js-1, je+2
1419  do i = is-1, ie+2
1420  do k = ks, ke
1421  pott(k,i,j) = rhot0(k,i,j) / dens0(k,i,j)
1422  enddo
1423  enddo
1424  enddo
1425 
1426  do jjs = js, je, jblock
1427  jje = jjs+jblock-1
1428  do iis = is, ie, iblock
1429  iie = iis+iblock-1
1430 
1431  ! monotonic flux
1432  ! at (x, y, interface)
1433  call atmos_dyn_fvm_fluxz_xyz_ud1( tflx_lo(:,:,:,zdir), & ! (out)
1434  mflx_hi(:,:,:,zdir), pott, gsqrt(:,:,:,i_xyz), & ! (in)
1435  num_diff(:,:,:,i_rhot,zdir), & ! (in)
1436  cdz, & ! (in)
1437  iis-1, iie+1, jjs-1, jje+1 ) ! (in)
1438 
1439  ! at (u, y, layer)
1440  call atmos_dyn_fvm_fluxx_xyz_ud1( tflx_lo(:,:,:,xdir), & ! (out)
1441  mflx_hi(:,:,:,xdir), pott, gsqrt(:,:,:,i_uyz), & ! (in)
1442  num_diff(:,:,:,i_rhot,xdir), & ! (in)
1443  cdz, & ! (in)
1444  iis-1, iie+1, jjs-1, jje+1 ) ! (in)
1445 
1446  ! at (x, v, layer)
1447  call atmos_dyn_fvm_fluxy_xyz_ud1( tflx_lo(:,:,:,ydir), & ! (out)
1448  mflx_hi(:,:,:,ydir), pott, gsqrt(:,:,:,i_xvz), & ! (in)
1449  num_diff(:,:,:,i_rhot,ydir), & ! (in)
1450  cdz, & ! (in)
1451  iis-1, iie+1, jjs-1, jje+1 ) ! (in)
1452 
1453  enddo
1454  enddo
1455 
1456  call atmos_dyn_fct( tflx_anti, & ! (out)
1457  pott, dens0, dens_rk, & ! (out)
1458  tflx_hi, tflx_lo, & ! (in)
1459  mflx_hi, & ! (in)
1460  rcdz, rcdx, rcdy, & ! (in)
1461  gsqrt(:,:,:,i_xyz), & ! (in)
1462  mapf(:,:,:,i_xy), dtrk, & ! (in)
1463  flag_fct_along_stream ) ! (in)
1464 
1465  do jjs = js, je, jblock
1466  jje = jjs+jblock-1
1467  do iis = is, ie, iblock
1468  iie = iis+iblock-1
1469 
1470  !--- update rho*theta
1471  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
1472  do j = jjs, jje
1473  do i = iis, iie
1474  do k = ks, ke
1475 #ifdef DEBUG
1476  call check( __line__, rhot_rk(k,i,j) )
1477  call check( __line__, tflx_anti(k ,i ,j ,zdir) )
1478  call check( __line__, tflx_anti(k-1,i ,j ,zdir) )
1479  call check( __line__, tflx_anti(k ,i ,j ,xdir) )
1480  call check( __line__, tflx_anti(k ,i-1,j ,xdir) )
1481  call check( __line__, tflx_anti(k ,i ,j ,ydir) )
1482  call check( __line__, tflx_anti(k ,i ,j-1,ydir) )
1483 #endif
1484  rhot_rk(k,i,j) = rhot_rk(k,i,j) &
1485  + dtrk * ( ( tflx_anti(k,i,j,zdir) - tflx_anti(k-1,i ,j ,zdir) ) * rcdz(k) &
1486  + ( tflx_anti(k,i,j,xdir) - tflx_anti(k ,i-1,j ,xdir) ) * rcdx(i) &
1487  + ( tflx_anti(k,i,j,ydir) - tflx_anti(k ,i ,j-1,ydir) ) * rcdy(j) ) &
1488  * mapf(i,j,1,i_xy) * mapf(i,j,2,i_xy) &
1489  / gsqrt(k,i,j,i_xyz)
1490  enddo
1491  enddo
1492  enddo
1493 #ifdef DEBUG
1494  k = iundef; i = iundef; j = iundef
1495 #endif
1496  enddo
1497  enddo
1498  endif
1499 #endif
1500 
1501 #ifdef HIST_TEND
1502  if ( lhist ) then
1503  call hist_in(advcv_t(:,:,:,i_dens), 'DENS_t_advcv', 'tendency of density (vert. advection)', 'kg/m3/s' )
1504  call hist_in(advcv_t(:,:,:,i_momz), 'MOMZ_t_advcv', 'tendency of momentum z (vert. advection)', 'kg/m2/s2', zdim='half')
1505  call hist_in(advcv_t(:,:,:,i_momx), 'MOMX_t_advcv', 'tendency of momentum x (vert. advection)', 'kg/m2/s2', xdim='half')
1506  call hist_in(advcv_t(:,:,:,i_momy), 'MOMY_t_advcv', 'tendency of momentum y (vert. advection)', 'kg/m2/s2', ydim='half')
1507  call hist_in(advcv_t(:,:,:,i_rhot), 'RHOT_t_advcv', 'tendency of rho*theta (vert. advection)', 'K kg/m3/s' )
1508 
1509  call hist_in(advch_t(:,:,:,i_dens), 'DENS_t_advch', 'tendency of density (horiz. advection)', 'kg/m3/s' )
1510  call hist_in(advch_t(:,:,:,i_momz), 'MOMZ_t_advch', 'tendency of momentum z (horiz. advection)', 'kg/m2/s2', zdim='half')
1511  call hist_in(advch_t(:,:,:,i_momx), 'MOMX_t_advch', 'tendency of momentum x (horiz. advection)', 'kg/m2/s2', xdim='half')
1512  call hist_in(advch_t(:,:,:,i_momy), 'MOMY_t_advch', 'tendency of momentum y (horiz. advection)', 'kg/m2/s2', ydim='half')
1513  call hist_in(advch_t(:,:,:,i_rhot), 'RHOT_t_advch', 'tendency of rho*theta (horiz. advection)', 'K kg/m3/s' )
1514 
1515  call hist_in(pg_t(:,:,:,1), 'MOMZ_t_pg', 'tendency of momentum z (pressure gradient)', 'kg/m2/s2', zdim='half')
1516  call hist_in(pg_t(:,:,:,2), 'MOMX_t_pg', 'tendency of momentum x (pressure gradient)', 'kg/m2/s2', xdim='half')
1517  call hist_in(pg_t(:,:,:,3), 'MOMY_t_pg', 'tendency of momentum y (pressure gradient)', 'kg/m2/s2', ydim='half')
1518 
1519  call hist_in(wdmp_t(:,:,:), 'MOMZ_t_wdamp', 'tendency of momentum z (Raileight damping)', 'kg/m2/s2', zdim='half')
1520 
1521  call hist_in(ddiv_t(:,:,:,1), 'MOMZ_t_ddiv', 'tendency of momentum z (divergence damping)', 'kg/m2/s2', zdim='half')
1522  call hist_in(ddiv_t(:,:,:,2), 'MOMX_t_ddiv', 'tendency of momentum x (divergence damping)', 'kg/m2/s2', xdim='half')
1523  call hist_in(ddiv_t(:,:,:,3), 'MOMY_t_ddiv', 'tendency of momentum y (divergence damping)', 'kg/m2/s2', ydim='half')
1524 
1525  call hist_in(cf_t(:,:,:,1), 'MOMX_t_cf', 'tendency of momentum x (coliolis force)', 'kg/m2/s2', xdim='half')
1526  call hist_in(cf_t(:,:,:,2), 'MOMY_t_cf', 'tendency of momentum y (coliolis force)', 'kg/m2/s2', ydim='half')
1527  endif
1528 #endif
1529 
1530  return
integer, parameter, public i_rhot
Definition: scale_index.F90:35
subroutine, public atmos_dyn_fvm_fluxx_xyw_ud1(flux, mom, val, DENS, GSQRT, MAPF, num_diff, CDZ, IIS, IIE, JJS, JJE)
calculation X-flux at XYW
integer, public is
start point of inner domain: x, local
integer, public i_xvz
integer, public je
end point of inner domain: y, local
subroutine, public atmos_dyn_fvm_fluxy_xyz_ud1(flux, mflx, val, GSQRT, num_diff, CDZ, IIS, IIE, JJS, JJE)
calculation Y-flux at XYZ
subroutine, public atmos_dyn_fvm_fluxy_xyw_ud1(flux, mom, val, DENS, GSQRT, MAPF, num_diff, CDZ, IIS, IIE, JJS, JJE)
calculation Y-flux at XYW
subroutine, public atmos_dyn_fvm_fluxz_xyz_ud1(flux, mflx, val, GSQRT, num_diff, CDZ, IIS, IIE, JJS, JJE)
calculation z-flux at XYZ
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
subroutine, public atmos_dyn_fvm_fluxy_xvz_ud1(flux, mom, val, DENS, GSQRT, MAPF, num_diff, CDZ, IIS, IIE, JJS, JJE)
calculation Y-flux at XV
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
subroutine, public atmos_dyn_fvm_fluxx_uyz_ud1(flux, mom, val, DENS, GSQRT, MAPF, num_diff, CDZ, IIS, IIE, JJS, JJE)
calculation X-flux at UY
procedure(flux_z), pointer, public atmos_dyn_fvm_fluxz_xvz
module GRIDTRANS
subroutine, public atmos_dyn_fvm_fluxz_xvz_ud1(flux, mom, val, DENS, GSQRT, J33G, num_diff, CDZ, IIS, IIE, JJS, JJE)
calculation z-flux at XV
integer, public ka
of whole cells: z, local, with HALO
integer, public i_uy
real(rp), public const_pre00
pressure reference [Pa]
Definition: scale_const.F90:90
integer, public jblock
block size for cache blocking: y
integer, public i_xyw
integer, public jhalo
of halo cells: y
module COMMUNICATION
Definition: scale_comm.F90:23
subroutine, public atmos_dyn_fvm_fluxz_uyz_ud1(flux, mom, val, DENS, GSQRT, J33G, num_diff, CDZ, IIS, IIE, JJS, JJE)
calculation z-flux at UY
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
module Atmosphere / Dynamics common
integer, public i_uyw
subroutine, public atmos_dyn_fvm_fluxy_uyz_ud1(flux, mom, val, DENS, GSQRT, MAPF, num_diff, CDZ, IIS, IIE, JJS, JJE)
calculation Y-flux at UY
module CONSTANT
Definition: scale_const.F90:14
subroutine, public atmos_dyn_fvm_fluxx_xvz_ud1(flux, mom, val, DENS, GSQRT, MAPF, num_diff, CDZ, IIS, IIE, JJS, JJE)
calculation X-flux at XV
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
real(rp), public const_eps
small number
Definition: scale_const.F90:36
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
subroutine, public atmos_dyn_fvm_fluxz_xyw_ud1(flux, mom, val, DENS, GSQRT, J33G, num_diff, CDZ, FDZ, dtrk, IIS, IIE, JJS, JJE)
calculation z-flux at XYW
procedure(flux_phi), pointer, public atmos_dyn_fvm_fluxy_xyz
procedure(flux_mom), pointer, public atmos_dyn_fvm_fluxy_uyz
subroutine, public atmos_dyn_fvm_fluxx_xyz_ud1(flux, mflx, val, GSQRT, num_diff, CDZ, IIS, IIE, JJS, JJE)
calculation X-flux at XYZ
procedure(flux_wz), pointer, public atmos_dyn_fvm_fluxz_xyw
integer, public ihalo
of halo cells: x
module scale_atmos_dyn_fvm_flux_ud1
procedure(flux_mom), pointer, public atmos_dyn_fvm_fluxy_xvz
Here is the call graph for this function:
Here is the caller graph for this function: