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, 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)
 

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
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( io_l ) write(io_fid_log,*) '*** HEVE Register'
87 
88  if ( atmos_dyn_type /= 'FVM-HEVE' .AND. atmos_dyn_type /= 'HEVE' ) then
89  write(*,*) 'xxx ATMOS_DYN_TYPE is not FVM-HEVE. Check!'
90  call prc_mpistop
91  endif
92 
93  va_out = va_fvm_heve
94  var_name(:) = ""
95  var_desc(:) = ""
96  var_unit(:) = ""
97 
98  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 104 of file scale_atmos_dyn_tstep_short_fvm_heve.F90.

Referenced by scale_atmos_dyn_tstep_short::atmos_dyn_tstep_short_regist().

104  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), 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 126 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_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::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().

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