SCALE-RM
Functions/Subroutines
scale_atmos_dyn_tstep_short_fvm_hivi Module Reference

module Atmosphere / Dynamics RK More...

Functions/Subroutines

subroutine, public atmos_dyn_tstep_short_fvm_hivi_regist (ATMOS_DYN_TYPE, VA_out, VAR_NAME, VAR_DESC, VAR_UNIT)
 Register. More...
 
subroutine, public atmos_dyn_tstep_short_fvm_hivi_setup
 Setup. More...
 
subroutine, public atmos_dyn_tstep_short_fvm_hivi (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, TwoD, dtrk, last)
 
subroutine solve_bicgstab (DPRES_N, DPRES, M, B, TwoD)
 
subroutine make_matrix (M, POTT, RCs2T, GRAV, G, J33G, RCDZ, RFDZ, RCDX, RFDX, RCDY, RFDY, FDZ, rdt, FACT_N, FACT_F, TwoD, I_XYZ, I_XYW, IIS, IIE, JJS, JJE)
 

Detailed Description

module Atmosphere / Dynamics RK

Description
HIVI FVM scheme for Atmospheric dynamical process
Author
Team SCALE
History
  • 2014-04-15 (S.Nishizawa) [new] newly impremented
NAMELIST
  • PARAM_ATMOS_DYN_TSTEP_FVM_HIVI
    nametypedefault valuecomment
    ITMAX integer
    EPSILON real(RP)

History Output
No history output

Function/Subroutine Documentation

◆ atmos_dyn_tstep_short_fvm_hivi_regist()

subroutine, public scale_atmos_dyn_tstep_short_fvm_hivi::atmos_dyn_tstep_short_fvm_hivi_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 73 of file scale_atmos_dyn_tstep_short_fvm_hivi.F90.

79  use scale_prc, only: &
80  prc_abort
81  implicit none
82 
83  character(len=*), intent(in) :: ATMOS_DYN_TYPE
84  integer, intent(out) :: VA_out
85  character(len=H_SHORT), intent(out) :: VAR_NAME(:)
86  character(len=H_MID), intent(out) :: VAR_DESC(:)
87  character(len=H_SHORT), intent(out) :: VAR_UNIT(:)
88  !---------------------------------------------------------------------------
89 
90  if ( atmos_dyn_type /= 'FVM-HIVI' .AND. atmos_dyn_type /= 'HIVI' ) then
91  log_error("ATMOS_DYN_Tstep_short_fvm_hivi_regist",*) 'ATMOS_DYN_TYPE is not FVM-HIVI. Check!'
92  call prc_abort
93  endif
94 
95  va_out = va_fvm_hivi
96  var_name(:) = ""
97  var_desc(:) = ""
98  var_unit(:) = ""
99 
100  log_newline
101  log_info("ATMOS_DYN_Tstep_short_fvm_hivi_regist",*) 'Register additional prognostic variables (HIVI)'
102  if ( va_out < 1 ) then
103  log_info_cont(*) '=> nothing.'
104  endif
105 
106  return
scale_prc
module PROCESS
Definition: scale_prc.F90:11
scale_prc::prc_abort
subroutine, public prc_abort
Abort Process.
Definition: scale_prc.F90:350

References scale_prc::prc_abort().

Referenced by scale_atmos_dyn_tstep_short::atmos_dyn_tstep_short_regist().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ atmos_dyn_tstep_short_fvm_hivi_setup()

subroutine, public scale_atmos_dyn_tstep_short_fvm_hivi::atmos_dyn_tstep_short_fvm_hivi_setup

Setup.

Definition at line 111 of file scale_atmos_dyn_tstep_short_fvm_hivi.F90.

112  use scale_prc, only: &
113  prc_abort
114  implicit none
115 
116  namelist / param_atmos_dyn_tstep_fvm_hivi / &
117  itmax, &
118  epsilon
119 
120  integer :: ierr
121  !---------------------------------------------------------------------------
122 
123  log_info("ATMOS_DYN_Tstep_short_fvm_hivi_setup",*) 'HIVI Setup'
124 #ifdef HIVI_BICGSTAB
125  log_info("ATMOS_DYN_Tstep_short_fvm_hivi_setup",*) 'USING Bi-CGSTAB'
126 #else
127  log_info("ATMOS_DYN_Tstep_short_fvm_hivi_setup",*) 'USING Multi-Grid'
128  log_error("ATMOS_DYN_Tstep_short_fvm_hivi_setup",*) 'Not Implemented yet'
129  call prc_abort
130 #endif
131 
132  ! currently, vertical difference scheme for potential temperature is the CD4
133  itmax = 100
134  epsilon = 0.1_rp ** (rp*2)
135 
136  !--- read namelist
137  rewind(io_fid_conf)
138  read(io_fid_conf,nml=param_atmos_dyn_tstep_fvm_hivi,iostat=ierr)
139 
140  if( ierr < 0 ) then !--- missing
141  log_info("ATMOS_DYN_Tstep_short_fvm_hivi_setup",*) 'Not found namelist. Default used.'
142  elseif( ierr > 0 ) then !--- fatal error
143  log_error("ATMOS_DYN_Tstep_short_fvm_hivi_setup",*) 'Not appropriate names in namelist PARAM_ATMOS_DYN_TSTEP_FVM_HIVI. Check!'
144  call prc_abort
145  endif
146  log_nml(param_atmos_dyn_tstep_fvm_hivi)
147 
148  if ( rp == dp ) then
149  mtype = mpi_double_precision
150  elseif( rp == sp ) then
151  mtype = mpi_real
152  else
153  log_error("ATMOS_DYN_Tstep_short_fvm_hivi_setup",*) 'Unsupported precision'
154  call prc_abort
155  endif
156 
157  return

References scale_precision::dp, scale_io::io_fid_conf, scale_prc::prc_abort(), and scale_precision::sp.

Referenced by scale_atmos_dyn_tstep_short::atmos_dyn_tstep_short_regist().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ atmos_dyn_tstep_short_fvm_hivi()

subroutine, public scale_atmos_dyn_tstep_short_fvm_hivi::atmos_dyn_tstep_short_fvm_hivi ( 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(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,
logical, intent(in)  TwoD,
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 161 of file scale_atmos_dyn_tstep_short_fvm_hivi.F90.

179  use scale_atmos_grid_cartesc_index
180  use scale_const, only: &
181  grav => const_grav, &
182  p00 => const_pre00
183  use scale_comm_cartesc, only: &
184  comm_vars8, &
185  comm_wait
186  use scale_atmos_dyn_fvm_fct, only: &
187  atmos_dyn_fvm_fct
188  use scale_atmos_dyn_fvm_flux, only: &
189  atmos_dyn_fvm_flux_valuew_z, &
190  atmos_dyn_fvm_fluxz_xyz, &
191  atmos_dyn_fvm_fluxx_xyz, &
192  atmos_dyn_fvm_fluxy_xyz, &
193  atmos_dyn_fvm_fluxz_xyw, &
194  atmos_dyn_fvm_fluxj13_xyw, &
195  atmos_dyn_fvm_fluxj23_xyw, &
196  atmos_dyn_fvm_fluxx_xyw, &
197  atmos_dyn_fvm_fluxy_xyw, &
198  atmos_dyn_fvm_fluxz_uyz, &
199  atmos_dyn_fvm_fluxj13_uyz, &
200  atmos_dyn_fvm_fluxj23_uyz, &
201  atmos_dyn_fvm_fluxx_uyz, &
202  atmos_dyn_fvm_fluxy_uyz, &
203  atmos_dyn_fvm_fluxz_xvz, &
204  atmos_dyn_fvm_fluxj13_xvz, &
205  atmos_dyn_fvm_fluxj23_xvz, &
206  atmos_dyn_fvm_fluxx_xvz, &
207  atmos_dyn_fvm_fluxy_xvz
208  implicit none
209 
210  real(RP), intent(out) :: DENS_RK(KA,IA,JA) ! prognostic variables
211  real(RP), intent(out) :: MOMZ_RK(KA,IA,JA) !
212  real(RP), intent(out) :: MOMX_RK(KA,IA,JA) !
213  real(RP), intent(out) :: MOMY_RK(KA,IA,JA) !
214  real(RP), intent(out) :: RHOT_RK(KA,IA,JA) !
215 
216  real(RP), intent(out) :: PROG_RK(KA,IA,JA,VA) !
217 
218  real(RP), intent(inout) :: mflx_hi(KA,IA,JA,3) ! rho * vel(x,y,z)
219  real(RP), intent(out) :: tflx_hi(KA,IA,JA,3) ! rho * theta * vel(x,y,z)
220 
221  real(RP), intent(in),target :: DENS0(KA,IA,JA) ! prognostic variables
222  real(RP), intent(in),target :: MOMZ0(KA,IA,JA) ! at previous dynamical time step
223  real(RP), intent(in),target :: MOMX0(KA,IA,JA) !
224  real(RP), intent(in),target :: MOMY0(KA,IA,JA) !
225  real(RP), intent(in),target :: RHOT0(KA,IA,JA) !
226 
227  real(RP), intent(in) :: DENS(KA,IA,JA) ! prognostic variables
228  real(RP), intent(in) :: MOMZ(KA,IA,JA) ! at previous RK step
229  real(RP), intent(in) :: MOMX(KA,IA,JA) !
230  real(RP), intent(in) :: MOMY(KA,IA,JA) !
231  real(RP), intent(in) :: RHOT(KA,IA,JA) !
232 
233  real(RP), intent(in) :: DENS_t(KA,IA,JA)
234  real(RP), intent(in) :: MOMZ_t(KA,IA,JA)
235  real(RP), intent(in) :: MOMX_t(KA,IA,JA)
236  real(RP), intent(in) :: MOMY_t(KA,IA,JA)
237  real(RP), intent(in) :: RHOT_t(KA,IA,JA)
238 
239  real(RP), intent(in) :: PROG0(KA,IA,JA,VA)
240  real(RP), intent(in) :: PROG (KA,IA,JA,VA)
241 
242  real(RP), intent(in) :: DPRES0(KA,IA,JA)
243  real(RP), intent(in) :: RT2P(KA,IA,JA)
244  real(RP), intent(in) :: CORIOLI(IA,JA)
245  real(RP), intent(in) :: num_diff(KA,IA,JA,5,3)
246  real(RP), intent(in) :: wdamp_coef(KA)
247  real(RP), intent(in) :: divdmp_coef
248  real(RP), intent(in) :: DDIV(KA,IA,JA)
249 
250  logical, intent(in) :: FLAG_FCT_MOMENTUM
251  logical, intent(in) :: FLAG_FCT_T
252  logical, intent(in) :: FLAG_FCT_ALONG_STREAM
253 
254  real(RP), intent(in) :: CDZ(KA)
255  real(RP), intent(in) :: FDZ(KA-1)
256  real(RP), intent(in) :: FDX(IA-1)
257  real(RP), intent(in) :: FDY(JA-1)
258  real(RP), intent(in) :: RCDZ(KA)
259  real(RP), intent(in) :: RCDX(IA)
260  real(RP), intent(in) :: RCDY(JA)
261  real(RP), intent(in) :: RFDZ(KA-1)
262  real(RP), intent(in) :: RFDX(IA-1)
263  real(RP), intent(in) :: RFDY(JA-1)
264 
265  real(RP), intent(in) :: PHI (KA,IA,JA)
266  real(RP), intent(in) :: GSQRT (KA,IA,JA,7)
267  real(RP), intent(in) :: J13G (KA,IA,JA,7)
268  real(RP), intent(in) :: J23G (KA,IA,JA,7)
269  real(RP), intent(in) :: J33G
270  real(RP), intent(in) :: MAPF (IA,JA,2,4)
271  real(RP), intent(in) :: REF_dens(KA,IA,JA)
272  real(RP), intent(in) :: REF_rhot(KA,IA,JA)
273 
274  logical, intent(in) :: BND_W
275  logical, intent(in) :: BND_E
276  logical, intent(in) :: BND_S
277  logical, intent(in) :: BND_N
278  logical, intent(in) :: TwoD
279 
280  real(RP), intent(in) :: dtrk
281  logical, intent(in) :: last
282 
283 
284  ! diagnostic variables (work space)
285  real(RP) :: POTT(KA,IA,JA) ! potential temperature [K]
286  real(RP) :: DDENS(KA,IA,JA) ! density deviation from reference
287  real(RP) :: DPRES(KA,IA,JA) ! pressure deviation from reference
288  real(RP) :: DPRES_N(KA,IA,JA) ! pressure deviation at t=n+1 [Pa]
289 
290  real(RP) :: Sr(KA,IA,JA)
291  real(RP) :: Sw(KA,IA,JA)
292  real(RP) :: Su(KA,IA,JA)
293  real(RP) :: Sv(KA,IA,JA)
294  real(RP) :: St(KA,IA,JA)
295  real(RP) :: RCs2T(KA,IA,JA)
296  real(RP) :: B(KA,IA,JA)
297 
298  real(RP) :: duvw
299 
300  real(RP) :: qflx_hi (KA,IA,JA,3)
301  real(RP) :: qflx_J13(KA,IA,JA)
302  real(RP) :: qflx_J23(KA,IA,JA)
303  real(RP) :: mflx_hi2(KA,IA,JA,3)
304 
305  ! for implicit solver
306  real(RP) :: M(7,KA,IA,JA)
307  real(RP) :: r(KA,IA,JA)
308  real(RP) :: p(KA,IA,JA)
309 
310  real(RP) :: zero(KA,IA,JA)
311 
312  real(RP) :: rdt
313 
314  integer :: IIS, IIE, JJS, JJE
315  integer :: k, i, j
316 
317  zero = 0.0_rp
318 
319 #ifdef DEBUG
320  dpres(:,:,:) = undef
321  pott(:,:,:) = undef
322 
323  qflx_j13(:,:,:) = undef
324  qflx_j23(:,:,:) = undef
325  mflx_hi2(:,:,:,:) = undef
326 
327  sr(:,:,:) = undef
328  sw(:,:,:) = undef
329  su(:,:,:) = undef
330  sv(:,:,:) = undef
331  st(:,:,:) = undef
332  rcs2t(:,:,:) = undef
333 
334  b(:,:,:) = undef
335  r(:,:,:) = undef
336  p(:,:,:) = undef
337 #endif
338 
339 #if defined DEBUG || defined QUICKDEBUG
340  dens_rk( 1:ks-1,:,:) = undef
341  dens_rk(ke+1:ka ,:,:) = undef
342  momz_rk( 1:ks-1,:,:) = undef
343  momz_rk(ke+1:ka ,:,:) = undef
344  momx_rk( 1:ks-1,:,:) = undef
345  momx_rk(ke+1:ka ,:,:) = undef
346  momy_rk( 1:ks-1,:,:) = undef
347  momy_rk(ke+1:ka ,:,:) = undef
348  rhot_rk( 1:ks-1,:,:) = undef
349  rhot_rk(ke+1:ka ,:,:) = undef
350  prog_rk( 1:ks-1,:,:,:) = undef
351  prog_rk(ke+1:ka ,:,:,:) = undef
352 #endif
353 
354  rdt = 1.0_rp / dtrk
355 
356  do jjs = js, je, jblock
357  jje = jjs+jblock-1
358  do iis = is, ie, iblock
359  iie = iis+iblock-1
360 
361  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
362  do j = jjs-1, jje+1
363  do i = max(iis-1,1), min(iie+1,ia)
364  do k = ks, ke
365 #ifdef DEBUG
366  call check( __line__, dpres0(k,i,j) )
367  call check( __line__, rt2p(k,i,j) )
368  call check( __line__, rhot(k,i,j) )
369  call check( __line__, ref_rhot(k,i,j) )
370 #endif
371  dpres(k,i,j) = dpres0(k,i,j) + rt2p(k,i,j) * ( rhot(k,i,j) - ref_rhot(k,i,j) )
372  enddo
373  dpres(ks-1,i,j) = dpres0(ks-1,i,j) - dens(ks,i,j) * ( phi(ks-1,i,j) - phi(ks+1,i,j) )
374  dpres(ke+1,i,j) = dpres0(ke+1,i,j) - dens(ke,i,j) * ( phi(ke+1,i,j) - phi(ke-1,i,j) )
375  enddo
376  enddo
377 
378  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
379  do j = jjs-1, jje+1
380  do i = iis-1, iie+1
381  do k = ks, ke
382 #ifdef DEBUG
383  call check( __line__, dens(k,i,j) )
384  call check( __line__, ref_dens(k,i,j) )
385 #endif
386  ddens(k,i,j) = dens(k,i,j) - ref_dens(k,i,j)
387  enddo
388  enddo
389  enddo
390 
391  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
392  do j = jjs-2, jje+2
393  do i = iis-2, iie+2
394  do k = ks, ke
395 #ifdef DEBUG
396  call check( __line__, rhot(k,i,j) )
397  call check( __line__, dens(k,i,j) )
398 #endif
399  pott(k,i,j) = rhot(k,i,j) / dens(k,i,j)
400  enddo
401  enddo
402  enddo
403 #ifdef DEBUG
404  k = iundef; i = iundef; j = iundef
405 #endif
406 
407  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
408  do j = jjs, jje
409  do i = iis, iie
410  do k = ks, ke
411 #ifdef DEBUG
412  call check( __line__, rt2p(k,i,j) )
413 #endif
414  rcs2t(k,i,j) = 1.0_rp / rt2p(k,i,j)
415  enddo
416  enddo
417  enddo
418 #ifdef DEBUG
419  k = iundef; i = iundef; j = iundef
420 #endif
421 
422  !##### continuity equation #####
423 
424  ! at (x, y, w)
425  if ( twod ) then
426  !$omp parallel do private(j,k) OMP_SCHEDULE_
427  do j = jjs, jje
428  do k = ks, ke-1
429 #ifdef DEBUG
430  call check( __line__, momx(k+1,is,j) )
431  call check( __line__, momx(k ,is,j) )
432  call check( __line__, momy(k+1,is,j) )
433  call check( __line__, momy(k+1,is,j-1) )
434  call check( __line__, momy(k ,is,j) )
435  call check( __line__, momy(k ,is,j-1) )
436  call check( __line__, num_diff(k,is,j,i_dens,zdir) )
437  call check( __line__, gsqrt(k,is,j,i_xyw) )
438 #endif
439  mflx_hi2(k,is,j,zdir) = j23g(k,is,j,i_xyw) * 0.25_rp * ( momy(k+1,is,j)+momy(k+1,is,j-1) &
440  + momy(k ,is,j)+momy(k ,is,j-1) ) &
441  + gsqrt(k,is,j,i_xyw) * num_diff(k,is,j,i_dens,zdir)
442  enddo
443  enddo
444  else
445  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
446  do j = jjs, jje
447  do i = iis-1, iie
448  do k = ks, ke-1
449 #ifdef DEBUG
450  call check( __line__, momx(k+1,i ,j) )
451  call check( __line__, momx(k+1,i-1,j) )
452  call check( __line__, momx(k ,i ,j) )
453  call check( __line__, momx(k ,i+1,j) )
454  call check( __line__, momy(k+1,i,j) )
455  call check( __line__, momy(k+1,i,j-1) )
456  call check( __line__, momy(k ,i,j) )
457  call check( __line__, momy(k ,i,j-1) )
458  call check( __line__, num_diff(k,i,j,i_dens,zdir) )
459  call check( __line__, gsqrt(k,i,j,i_xyw) )
460 #endif
461  mflx_hi2(k,i,j,zdir) = j13g(k,i,j,i_xyw) * 0.25_rp * ( momx(k+1,i,j)+momx(k+1,i-1,j) &
462  + momx(k ,i,j)+momx(k ,i-1,j) ) &
463  + j23g(k,i,j,i_xyw) * 0.25_rp * ( momy(k+1,i,j)+momy(k+1,i,j-1) &
464  + momy(k ,i,j)+momy(k ,i,j-1) ) &
465  + gsqrt(k,i,j,i_xyw) * num_diff(k,i,j,i_dens,zdir)
466  enddo
467  enddo
468  enddo
469  end if
470 #ifdef DEBUG
471  k = iundef; i = iundef; j = iundef
472 #endif
473 
474  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
475  do j = jjs, jje
476  do i = iis, iie
477  mflx_hi2(ks-1,i,j,zdir) = 0.0_rp
478  mflx_hi2(ke ,i,j,zdir) = 0.0_rp
479  enddo
480  enddo
481 
482  ! at (u, y, z)
483  if ( .not. twod ) then
484  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
485  do j = jjs , jje
486  do i = iis-1, iie
487  do k = ks, ke
488 #ifdef DEBUG
489  call check( __line__, gsqrt(k,i,j,i_uyz) )
490  call check( __line__, num_diff(k,i,j,i_dens,xdir) )
491 #endif
492  mflx_hi2(k,i,j,xdir) = gsqrt(k,i,j,i_uyz) * num_diff(k,i,j,i_dens,xdir)
493  enddo
494  enddo
495  enddo
496  end if
497 #ifdef DEBUG
498  k = iundef; i = iundef; j = iundef
499 #endif
500  ! at (x, v, z)
501  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
502  do j = jjs-1, jje
503  do i = iis , iie
504  do k = ks, ke
505 #ifdef DEBUG
506  call check( __line__, gsqrt(k,i,j,i_xvz) )
507  call check( __line__, num_diff(k,i,j,i_dens,ydir) )
508 #endif
509  mflx_hi2(k,i,j,ydir) = gsqrt(k,i,j,i_xvz) * num_diff(k,i,j,i_dens,ydir)
510  enddo
511  enddo
512  enddo
513 #ifdef DEBUG
514  k = iundef; i = iundef; j = iundef
515 #endif
516 
517  !##### momentum equation (z) #####
518 
519  ! at (x, y, z)
520  ! note than z-index is added by -1
521  call atmos_dyn_fvm_fluxz_xyw( qflx_hi(:,:,:,zdir), & ! (out)
522  momz, momz, dens, & ! (in)
523  gsqrt(:,:,:,i_xyz), j33g, & ! (in)
524  num_diff(:,:,:,i_momz,zdir), & ! (in)
525  cdz, fdz, dtrk, &
526  iis, iie, jjs, jje ) ! (in)
527  if ( .not. twod ) &
528  call atmos_dyn_fvm_fluxj13_xyw( qflx_j13, & ! (out)
529  momx, momz, dens, & ! (in)
530  gsqrt(:,:,:,i_xyz), j13g(:,:,:,i_xyz), mapf(:,:,:,i_xy), & ! (in)
531  cdz, twod, &
532  iis, iie, jjs, jje ) ! (in)
533  call atmos_dyn_fvm_fluxj23_xyw( qflx_j23, & ! (out)
534  momy, momz, dens, & ! (in)
535  gsqrt(:,:,:,i_xyz), j23g(:,:,:,i_xyz), mapf(:,:,:,i_xy), & ! (in)
536  cdz, twod, &
537  iis, iie, jjs, jje ) ! (in)
538 
539  if ( .not. twod ) &
540  call atmos_dyn_fvm_fluxx_xyw( qflx_hi(:,:,:,xdir), & ! (out)
541  momx, momz, dens, & ! (in)
542  gsqrt(:,:,:,i_uyw), mapf(:,:,:,i_uy), & ! (in)
543  num_diff(:,:,:,i_momz,xdir), & ! (in)
544  cdz, twod, & ! (in)
545  iis, iie, jjs, jje ) ! (in)
546 
547  ! at (x, v, w)
548  call atmos_dyn_fvm_fluxy_xyw( qflx_hi(:,:,:,ydir), & ! (out)
549  momy, momz, dens, & ! (in)
550  gsqrt(:,:,:,i_xvw), mapf(:,:,:,i_xv), & ! (in)
551  num_diff(:,:,:,i_momz,ydir), & ! (in)
552  cdz, twod, & ! (in)
553  iis, iie, jjs, jje ) ! (in)
554 
555  !--- update momentum(z)
556  if ( twod ) then
557  !$omp parallel do private(j,k) OMP_SCHEDULE_
558  do j = jjs, jje
559  do k = ks, ke-1
560 #ifdef DEBUG
561  call check( __line__, qflx_hi(k ,is,j ,zdir) )
562  call check( __line__, qflx_hi(k-1,is,j ,zdir) )
563  call check( __line__, qflx_j23(k ,is,j ) )
564  call check( __line__, qflx_j23(k-1,is,j ) )
565  call check( __line__, qflx_hi(k ,is,j ,ydir) )
566  call check( __line__, qflx_hi(k ,is,j-1,ydir) )
567  call check( __line__, ddiv(k ,is,j) )
568  call check( __line__, ddiv(k+1,is,j) )
569  call check( __line__, momz0(k,is,j) )
570  call check( __line__, momz_t(k,is,j) )
571 #endif
572  sw(k,is,j) = &
573  - ( ( qflx_hi(k,is,j,zdir) - qflx_hi(k-1,is,j ,zdir) &
574  + qflx_j23(k,is,j) - qflx_j23(k-1,is,j ) ) * rfdz(k) &
575  + ( qflx_hi(k,is,j,ydir) - qflx_hi(k ,is,j-1,ydir) ) * rcdy(j) &
576  ) / gsqrt(k,is,j,i_xyw) &
577  - wdamp_coef(k) * momz0(k,is,j) & ! Rayleigh damping
578  + divdmp_coef * rdt * ( ddiv(k+1,is,j)-ddiv(k,is,j) ) * fdz(k) & ! divergence damping
579  + momz_t(k,is,j)
580  enddo
581  enddo
582  else
583  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
584  do j = jjs, jje
585  do i = iis, iie
586  do k = ks, ke-1
587 #ifdef DEBUG
588  call check( __line__, qflx_hi(k ,i ,j ,zdir) )
589  call check( __line__, qflx_hi(k-1,i ,j ,zdir) )
590  call check( __line__, qflx_j13(k ,i ,j ) )
591  call check( __line__, qflx_j13(k-1,i ,j ) )
592  call check( __line__, qflx_j23(k ,i ,j ) )
593  call check( __line__, qflx_j23(k-1,i ,j ) )
594  call check( __line__, qflx_hi(k ,i ,j ,xdir) )
595  call check( __line__, qflx_hi(k ,i-1,j ,xdir) )
596  call check( __line__, qflx_hi(k ,i ,j ,ydir) )
597  call check( __line__, qflx_hi(k ,i ,j-1,ydir) )
598  call check( __line__, ddiv(k ,i,j) )
599  call check( __line__, ddiv(k+1,i,j) )
600  call check( __line__, momz0(k,i,j) )
601  call check( __line__, momz_t(k,i,j) )
602 #endif
603  sw(k,i,j) = &
604  - ( ( qflx_hi(k,i,j,zdir) - qflx_hi(k-1,i ,j ,zdir) &
605  + qflx_j13(k,i,j) - qflx_j13(k-1,i ,j ) &
606  + qflx_j23(k,i,j) - qflx_j23(k-1,i ,j ) ) * rfdz(k) &
607  + ( qflx_hi(k,i,j,xdir) - qflx_hi(k ,i-1,j ,xdir) ) * rcdx(i) &
608  + ( qflx_hi(k,i,j,ydir) - qflx_hi(k ,i ,j-1,ydir) ) * rcdy(j) &
609  ) / gsqrt(k,i,j,i_xyw) &
610  - wdamp_coef(k) * momz0(k,i,j) & ! Rayleigh damping
611  + divdmp_coef * rdt * ( ddiv(k+1,i,j)-ddiv(k,i,j) ) * fdz(k) & ! divergence damping
612  + momz_t(k,i,j)
613  enddo
614  enddo
615  enddo
616  end if
617 #ifdef DEBUG
618  k = iundef; i = iundef; j = iundef
619 #endif
620  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
621  do j = jjs, jje
622  do i = iis, iie
623  sw(ks-1,i,j) = 0.0_rp
624  sw(ke ,i,j) = 0.0_rp
625  enddo
626  enddo
627 #ifdef DEBUG
628  k = iundef; i = iundef; j = iundef
629 #endif
630 
631 
632  !##### momentum equation (x) #####
633 
634  ! at (u, y, w)
635  call atmos_dyn_fvm_fluxz_uyz( qflx_hi(:,:,:,zdir), & ! (out)
636  momz, momx, dens, & ! (in)
637  gsqrt(:,:,:,i_uyw), j33g, & ! (in)
638  num_diff(:,:,:,i_momx,zdir), & ! (in)
639  cdz, twod, & ! (in)
640  iis, iie, jjs, jje ) ! (in)
641  if ( .not. twod ) &
642  call atmos_dyn_fvm_fluxj13_uyz( qflx_j13, & ! (out)
643  momx, momx, dens, & ! (in)
644  gsqrt(:,:,:,i_uyz), j13g(:,:,:,i_uyw), mapf(:,:,:,i_uy), & ! (in)
645  cdz, twod, & ! (in)
646  iis, iie, jjs, jje ) ! (in)
647  call atmos_dyn_fvm_fluxj23_uyz( qflx_j23, & ! (out)
648  momy, momx, dens, & ! (in)
649  gsqrt(:,:,:,i_uyz), j23g(:,:,:,i_uyw), mapf(:,:,:,i_uy), & ! (in)
650  cdz, twod, & ! (in)
651  iis, iie, jjs, jje ) ! (in)
652 
653  ! at (x, y, z)
654  ! note that x-index is added by -1
655  if ( .not. twod ) &
656  call atmos_dyn_fvm_fluxx_uyz( qflx_hi(:,:,:,xdir), & ! (out)
657  momx, momx, dens, & ! (in)
658  gsqrt(:,:,:,i_xyz), mapf(:,:,:,i_xy), & ! (in)
659  num_diff(:,:,:,i_momx,xdir), & ! (in)
660  cdz, twod, & ! (in)
661  iis, iie, jjs, jje ) ! (in)
662 
663  ! at (u, v, z)
664  call atmos_dyn_fvm_fluxy_uyz( qflx_hi(:,:,:,ydir), & ! (out)
665  momy, momx, dens, & ! (in)
666  gsqrt(:,:,:,i_uvz), mapf(:,:,:,i_uv), & ! (in)
667  num_diff(:,:,:,i_momx,ydir), & ! (in)
668  cdz, twod, & ! (in)
669  iis, iie, jjs, jje ) ! (in)
670 
671  !--- update momentum(x)
672  if ( twod ) then
673  !$omp parallel do private(j,k) OMP_SCHEDULE_
674  do j = jjs, jje
675  do k = ks, ke
676 #ifdef DEBUG
677  call check( __line__, qflx_hi(k ,is,j ,zdir) )
678  call check( __line__, qflx_hi(k-1,is,j ,zdir) )
679  call check( __line__, qflx_hi(k ,is,j ,ydir) )
680  call check( __line__, qflx_hi(k ,is,j-1,ydir) )
681  call check( __line__, corioli(is,j) )
682  call check( __line__, momy(k,is,j ) )
683  call check( __line__, momy(k,is,j-1) )
684  call check( __line__, momx0(k,is,j) )
685 #endif
686  su(k,is,j) = ( &
687  - ( ( qflx_hi(k,is,j,zdir) - qflx_hi(k-1,is,j ,zdir) &
688  + qflx_j23(k,is,j) - qflx_j23(k-1,is,j) ) * rcdz(k) &
689  + ( qflx_hi(k,is,j,ydir) - qflx_hi(k ,is,j-1,ydir) ) * rcdy(j) ) &
690  ) / gsqrt(k,is,j,i_uyz) &
691  + 0.5_rp * corioli(is,j) * ( momy(k,is,j)+momy(k,is,j-1) ) & ! coriolis force
692  + momx_t(k,is,j)
693  enddo
694  enddo
695  else
696  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
697  do j = jjs, jje
698  do i = iis, iie
699  do k = ks, ke
700 #ifdef DEBUG
701  call check( __line__, qflx_hi(k ,i ,j ,zdir) )
702  call check( __line__, qflx_hi(k-1,i ,j ,zdir) )
703  call check( __line__, qflx_hi(k ,i ,j ,xdir) )
704  call check( __line__, qflx_hi(k ,i-1,j ,xdir) )
705  call check( __line__, qflx_hi(k ,i ,j ,ydir) )
706  call check( __line__, qflx_hi(k ,i ,j-1,ydir) )
707  call check( __line__, dpres(k+1,i+1,j) )
708  call check( __line__, dpres(k+1,i ,j) )
709  call check( __line__, dpres(k-1,i+1,j) )
710  call check( __line__, dpres(k-1,i ,j) )
711  call check( __line__, corioli(i ,j) )
712  call check( __line__, corioli(i+1,j) )
713  call check( __line__, momy(k,i ,j ) )
714  call check( __line__, momy(k,i+1,j ) )
715  call check( __line__, momy(k,i ,j-1) )
716  call check( __line__, momy(k,i+1,j-1) )
717  call check( __line__, ddiv(k,i+1,j) )
718  call check( __line__, ddiv(k,i ,j) )
719  call check( __line__, momx0(k,i,j) )
720 #endif
721  su(k,i,j) = ( &
722  - ( ( qflx_hi(k,i,j,zdir) - qflx_hi(k-1,i ,j ,zdir) &
723  + qflx_j13(k,i,j) - qflx_j13(k-1,i ,j) &
724  + qflx_j23(k,i,j) - qflx_j23(k-1,i ,j) ) * rcdz(k) &
725  + ( qflx_hi(k,i,j,xdir) - qflx_hi(k ,i-1,j ,xdir) ) * rfdx(i) &
726  + ( qflx_hi(k,i,j,ydir) - qflx_hi(k ,i ,j-1,ydir) ) * rcdy(j) ) &
727  - ( j13g(k+1,i,j,i_uyz) * ( dpres(k+1,i+1,j)+dpres(k+1,i,j) ) &
728  - j13g(k-1,i,j,i_uyz) * ( dpres(k-1,i+1,j)+dpres(k-1,i,j) ) ) &
729  * 0.5_rp / ( fdz(k+1)+fdz(k) ) &
730  ) / gsqrt(k,i,j,i_uyz) &
731  + 0.125_rp * ( corioli(i,j)+corioli(i+1,j) ) &
732  * ( momy(k,i,j)+momy(k,i+1,j)+momy(k,i,j-1)+momy(k,i+1,j-1) ) & ! coriolis force
733  + divdmp_coef * rdt * ( ddiv(k,i+1,j)-ddiv(k,i,j) ) * fdx(i) & ! divergence damping
734  + momx_t(k,i,j)
735  enddo
736  enddo
737  enddo
738  end if
739 #ifdef DEBUG
740  k = iundef; i = iundef; j = iundef
741 #endif
742 
743  !##### momentum equation (y) #####
744 
745  ! at (x, v, w)
746  call atmos_dyn_fvm_fluxz_xvz( qflx_hi(:,:,:,zdir), & ! (out)
747  momz, momy, dens, & ! (in)
748  gsqrt(:,:,:,i_xvw), j33g, & ! (in)
749  num_diff(:,:,:,i_momy,zdir), & ! (in)
750  cdz, twod, & ! (in)
751  iis, iie, jjs, jje ) ! (in)
752  if ( .not. twod ) &
753  call atmos_dyn_fvm_fluxj13_xvz( qflx_j13, & ! (out)
754  momx, momy, dens, & ! (in)
755  gsqrt(:,:,:,i_xvz), j13g(:,:,:,i_xvw), mapf(:,:,:,i_xv), & ! (in)
756  cdz, twod, & ! (in)
757  iis, iie, jjs, jje ) ! (in)
758  call atmos_dyn_fvm_fluxj23_xvz( qflx_j23, & ! (out)
759  momy, momy, dens, & ! (in)
760  gsqrt(:,:,:,i_xvz), j23g(:,:,:,i_xvw), mapf(:,:,:,i_xv), & ! (in)
761  cdz, twod, & ! (in)
762  iis, iie, jjs, jje ) ! (in)
763 
764  ! at (u, v, z)
765  if ( .not. twod ) &
766  call atmos_dyn_fvm_fluxx_xvz( qflx_hi(:,:,:,xdir), & ! (out)
767  momx, momy, dens, & ! (in)
768  gsqrt(:,:,:,i_uvz), mapf(:,:,:,i_uv), & ! (in)
769  num_diff(:,:,:,i_momy,xdir), & ! (in)
770  cdz, twod, & ! (in)
771  iis, iie, jjs, jje ) ! (in)
772 
773  ! at (x, y, z)
774  ! note that y-index is added by -1
775  call atmos_dyn_fvm_fluxy_xvz( qflx_hi(:,:,:,ydir), & ! (out)
776  momy, momy, dens, & ! (in)
777  gsqrt(:,:,:,i_xyz), mapf(:,:,:,i_xy), & ! (in)
778  num_diff(:,:,:,i_momy,ydir), & ! (in
779  cdz, twod, & ! (in)
780  iis, iie, jjs, jje ) ! (in)
781 
782  !--- update momentum(y)
783  if ( twod ) then
784  !$omp parallel do private(j,k) OMP_SCHEDULE_
785  do j = jjs, jje
786  do k = ks, ke
787 #ifdef DEBUG
788  call check( __line__, qflx_hi(k ,is,j ,zdir) )
789  call check( __line__, qflx_hi(k-1,is,j ,zdir) )
790  call check( __line__, qflx_hi(k ,is,j ,ydir) )
791  call check( __line__, qflx_hi(k ,is,j-1,ydir) )
792  call check( __line__, dpres(k+1,is,j ) )
793  call check( __line__, dpres(k+1,is,j+1) )
794  call check( __line__, dpres(k-1,is,j ) )
795  call check( __line__, dpres(k-1,is,j+1) )
796  call check( __line__, corioli(is,j ) )
797  call check( __line__, corioli(is,j+1) )
798  call check( __line__, momx(k,is,j ) )
799  call check( __line__, momx(k,is,j+1) )
800  call check( __line__, ddiv(k,is,j+1) )
801  call check( __line__, ddiv(k,is,j ) )
802  call check( __line__, momy_t(k,is,j) )
803 #endif
804  sv(k,is,j) = &
805  ( - ( ( qflx_hi(k,is,j,zdir) - qflx_hi(k-1,is,j ,zdir) &
806  + qflx_j23(k,is,j) - qflx_j23(k-1,is,j ) ) * rcdz(k) &
807  + ( qflx_hi(k,is,j,ydir) - qflx_hi(k ,is,j-1,ydir) ) * rfdy(j) ) &
808  - ( j23g(k+1,is,j,i_xvz) * ( dpres(k+1,is,j+1)+dpres(k+1,is,j) ) &
809  - j23g(k-1,is,j,i_xvz) * ( dpres(k-1,is,j+1)+dpres(k-1,is,j) ) ) &
810  * 0.5_rp / ( fdz(k+1)+fdz(k) ) & ! pressure gradient force
811  ) / gsqrt(k,is,j,i_xvz) &
812  - 0.25_rp * ( corioli(is,j+1)+corioli(is,j) ) &
813  * ( momx(k,is,j+1)+momx(k,is,j) ) & ! coriolis force
814  + divdmp_coef * rdt * ( ddiv(k,is,j+1)-ddiv(k,is,j) ) * fdy(j) & ! divergence damping
815  + momy_t(k,is,j)
816  enddo
817  enddo
818  else
819  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
820  do j = jjs, jje
821  do i = iis, iie
822  do k = ks, ke
823 #ifdef DEBUG
824  call check( __line__, qflx_hi(k ,i ,j ,zdir) )
825  call check( __line__, qflx_hi(k-1,i ,j ,zdir) )
826  call check( __line__, qflx_hi(k ,i ,j ,xdir) )
827  call check( __line__, qflx_hi(k ,i-1,j ,xdir) )
828  call check( __line__, qflx_hi(k ,i ,j ,ydir) )
829  call check( __line__, qflx_hi(k ,i ,j-1,ydir) )
830  call check( __line__, dpres(k+1,i,j ) )
831  call check( __line__, dpres(k+1,i,j+1) )
832  call check( __line__, dpres(k-1,i,j ) )
833  call check( __line__, dpres(k-1,i,j+1) )
834  call check( __line__, corioli(i,j ) )
835  call check( __line__, corioli(i,j+1) )
836  call check( __line__, momx(k,i ,j ) )
837  call check( __line__, momx(k,i ,j+1) )
838  call check( __line__, momx(k,i-1,j ) )
839  call check( __line__, momx(k,i-1,j+1) )
840  call check( __line__, ddiv(k,i,j+1) )
841  call check( __line__, ddiv(k,i,j ) )
842  call check( __line__, momy_t(k,i,j) )
843 #endif
844  sv(k,i,j) = &
845  ( - ( ( qflx_hi(k,i,j,zdir) - qflx_hi(k-1,i ,j ,zdir) &
846  + qflx_j13(k,i,j) - qflx_j13(k-1,i ,j ) &
847  + qflx_j23(k,i,j) - qflx_j23(k-1,i ,j ) ) * rcdz(k) &
848  + ( qflx_hi(k,i,j,xdir) - qflx_hi(k ,i-1,j ,xdir) ) * rcdx(i) &
849  + ( qflx_hi(k,i,j,ydir) - qflx_hi(k ,i ,j-1,ydir) ) * rfdy(j) ) &
850  - ( j23g(k+1,i,j,i_xvz) * ( dpres(k+1,i,j+1)+dpres(k+1,i,j) ) &
851  - j23g(k-1,i,j,i_xvz) * ( dpres(k-1,i,j+1)+dpres(k-1,i,j) ) ) &
852  * 0.5_rp / ( fdz(k+1)+fdz(k) ) & ! pressure gradient force
853  ) / gsqrt(k,i,j,i_xvz) &
854  - 0.125_rp * ( corioli(i ,j+1)+corioli(i ,j) ) &
855  * ( momx(k,i,j+1)+momx(k,i,j)+momx(k,i-1,j+1)+momx(k,i-1,j) ) & ! coriolis force
856  + divdmp_coef * rdt * ( ddiv(k,i,j+1)-ddiv(k,i,j) ) * fdy(j) & ! divergence damping
857  + momy_t(k,i,j)
858  enddo
859  enddo
860  enddo
861  end if
862 #ifdef DEBUG
863  k = iundef; i = iundef; j = iundef
864 #endif
865 
866  !##### Thermodynamic Equation #####
867 
868 
869  ! at (x, y, w)
870  call atmos_dyn_fvm_fluxz_xyz( tflx_hi(:,:,:,zdir), & ! (out)
871  mflx_hi2(:,:,:,zdir), pott, gsqrt(:,:,:,i_xyw), & ! (in)
872  num_diff(:,:,:,i_rhot,zdir), & ! (in)
873  cdz, & ! (in)
874  iis, iie, jjs, jje ) ! (in)
875 
876  ! at (u, y, z)
877  if ( .not. twod ) &
878  call atmos_dyn_fvm_fluxx_xyz( tflx_hi(:,:,:,xdir), & ! (out)
879  mflx_hi2(:,:,:,xdir), pott, gsqrt(:,:,:,i_uyz), & ! (in)
880  num_diff(:,:,:,i_rhot,xdir), & ! (in)
881  cdz, & ! (in)
882  iis, iie, jjs, jje ) ! (in)
883 
884  ! at (x, v, z)
885  call atmos_dyn_fvm_fluxy_xyz( tflx_hi(:,:,:,ydir), & ! (out)
886  mflx_hi2(:,:,:,ydir), pott, gsqrt(:,:,:,i_xvz), & ! (in)
887  num_diff(:,:,:,i_rhot,ydir), & ! (in)
888  cdz, & ! (in)
889  iis, iie, jjs, jje ) ! (in)
890 
891  if ( twod ) then
892  !$omp parallel do private(j,k) OMP_SCHEDULE_
893  do j = jjs, jje
894  do k = ks, ke
895 #ifdef DEBUG
896  call check( __line__, tflx_hi(k ,is,j ,zdir) )
897  call check( __line__, tflx_hi(k-1,is,j ,zdir) )
898  call check( __line__, tflx_hi(k ,is,j ,ydir) )
899  call check( __line__, tflx_hi(k ,is,j-1,ydir) )
900  call check( __line__, rhot_t(k,is,j) )
901 #endif
902  st(k,is,j) = &
903  - ( ( tflx_hi(k,is,j,zdir) - tflx_hi(k-1,is,j ,zdir) ) * rcdz(k) &
904  + ( tflx_hi(k,is,j,ydir) - tflx_hi(k ,is,j-1,ydir) ) * rcdy(j) &
905  ) / gsqrt(k,is,j,i_xyz) &
906  + rhot_t(k,is,j)
907  enddo
908  enddo
909  else
910  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
911  do j = jjs, jje
912  do i = iis, iie
913  do k = ks, ke
914 #ifdef DEBUG
915  call check( __line__, tflx_hi(k ,i ,j ,zdir) )
916  call check( __line__, tflx_hi(k-1,i ,j ,zdir) )
917  call check( __line__, tflx_hi(k ,i ,j ,xdir) )
918  call check( __line__, tflx_hi(k ,i-1,j ,xdir) )
919  call check( __line__, tflx_hi(k ,i ,j ,ydir) )
920  call check( __line__, tflx_hi(k ,i ,j-1,ydir) )
921  call check( __line__, rhot_t(k,i,j) )
922 #endif
923  st(k,i,j) = &
924  - ( ( tflx_hi(k,i,j,zdir) - tflx_hi(k-1,i ,j ,zdir) ) * rcdz(k) &
925  + ( tflx_hi(k,i,j,xdir) - tflx_hi(k ,i-1,j ,xdir) ) * rcdx(i) &
926  + ( tflx_hi(k,i,j,ydir) - tflx_hi(k ,i ,j-1,ydir) ) * rcdy(j) &
927  ) / gsqrt(k,i,j,i_xyz) &
928  + rhot_t(k,i,j)
929  enddo
930  enddo
931  enddo
932  end if
933 #ifdef DEBUG
934  k = iundef; i = iundef; j = iundef
935 #endif
936 
937  enddo
938  enddo
939 
940  ! implicit solver
941 #ifdef HIVI_BICGSTAB
942 
943 #ifdef DEBUG
944  m(:,:,:,:) = undef
945 #endif
946 
947  call comm_vars8( su, 1 )
948  call comm_vars8( sv, 2 )
949  call comm_wait ( su, 1 )
950  call comm_wait ( sv, 2 )
951 
952  do jjs = js, je, jblock
953  jje = jjs+jblock-1
954  do iis = is, ie, iblock
955  iie = iis+iblock-1
956 
957  if ( twod ) then
958 
959  i = is
960  ! r = b - M x0
961  do j = jjs, jje
962  do k = ks+2, ke-2
963 #ifdef DEBUG
964  call check( __line__, momz(k-1,i,j) )
965  call check( __line__, momz(k ,i,j) )
966  call check( __line__, momx(k,i ,j) )
967  call check( __line__, momy(k,i,j-1) )
968  call check( __line__, momy(k,i,j ) )
969  call check( __line__, pott(k-2,i,j) )
970  call check( __line__, pott(k-1,i,j) )
971  call check( __line__, pott(k ,i,j) )
972  call check( __line__, pott(k+1,i,j) )
973  call check( __line__, pott(k+2,i,j) )
974  call check( __line__, pott(k,i ,j) )
975  call check( __line__, pott(k,i,j-2) )
976  call check( __line__, pott(k,i,j-1) )
977  call check( __line__, pott(k,i,j ) )
978  call check( __line__, pott(k,i,j+1) )
979  call check( __line__, pott(k,i,j+2) )
980  call check( __line__, sw(k-1,i,j) )
981  call check( __line__, sw(k ,i,j) )
982  call check( __line__, su(k,i-1,j) )
983  call check( __line__, su(k,i ,j) )
984  call check( __line__, sv(k,i,j-1) )
985  call check( __line__, sv(k,i,j ) )
986  call check( __line__, gsqrt(k,i ,j,i_uyz) )
987  call check( __line__, gsqrt(k,i,j-1,i_xvz) )
988  call check( __line__, gsqrt(k,i,j ,i_xvz) )
989  call check( __line__, gsqrt(k,i,j,i_xyz) )
990  call check( __line__, st(k,i,j) )
991  call check( __line__, dpres(k-1,i,j) )
992  call check( __line__, dpres(k ,i,j) )
993  call check( __line__, dpres(k+1,i,j) )
994  call check( __line__, rt2p(k-1,i,j) )
995  call check( __line__, rt2p(k ,i,j) )
996  call check( __line__, rt2p(k+1,i,j) )
997  call check( __line__, rhot(k-1,i,j) )
998  call check( __line__, rhot(k+1,i,j) )
999  call check( __line__, ddens(k-1,i,j) )
1000  call check( __line__, ddens(k+1,i,j) )
1001 #endif
1002  b(k,i,j) = ( &
1003  ( j33g * ( momz(k ,i,j) + dtrk*sw(k ,i,j) ) &
1004  * ( fact_n*(pott(k+1,i,j)+pott(k ,i,j)) &
1005  + fact_f*(pott(k+2,i,j)+pott(k-1,i,j)) ) &
1006  - j33g * ( momz(k-1,i,j) + dtrk*sw(k-1,i,j) ) &
1007  * ( fact_n*(pott(k ,i,j)+pott(k-1,i,j)) &
1008  + fact_f*(pott(k+1,i,j)+pott(k-2,i,j)) ) ) * rcdz(k) &
1009  + ( gsqrt(k,i,j ,i_xvz) * ( momy(k,i,j ) + dtrk*sv(k,i,j ) ) &
1010  * ( fact_n*(pott(k,i,j+1)+pott(k,i,j )) &
1011  + fact_f*(pott(k,i,j+2)+pott(k,i,j-1)) ) &
1012  - gsqrt(k,i,j-1,i_xvz) * ( momy(k,i,j-1) + dtrk*sv(k,i,j-1) ) &
1013  * ( fact_n*(pott(k,i,j )+pott(k,i,j-1)) &
1014  + fact_f*(pott(k,i,j+1)+pott(k,i,j-2)) ) ) * rcdy(j) &
1015  + gsqrt(k,i,j,i_xyz) * ( st(k,i,j) - dpres(k,i,j) * rcs2t(k,i,j) * rdt ) &
1016  ) * rdt &
1017  + grav * j33g * ( ( dpres(k+1,i,j)*rcs2t(k+1,i,j) &
1018  - dpres(k-1,i,j)*rcs2t(k-1,i,j) ) &
1019  - ( rhot(k+1,i,j)*ddens(k+1,i,j)/dens(k+1,i,j) &
1020  - rhot(k-1,i,j)*ddens(k-1,i,j)/dens(k-1,i,j) ) &
1021  ) / ( fdz(k) + fdz(k-1) )
1022  enddo
1023  enddo
1024 #ifdef DEBUG
1025  k = iundef; i = iundef; j = iundef
1026 #endif
1027  do j = jjs, jje
1028  do i = iis, iie
1029 #ifdef DEBUG
1030  call check( __line__, momz(ks,i,j) )
1031  call check( __line__, momx(ks,i,j) )
1032  call check( __line__, momy(ks,i,j) )
1033  call check( __line__, sw(ks,i,j) )
1034  call check( __line__, su(ks,i-1,j) )
1035  call check( __line__, su(ks,i ,j) )
1036  call check( __line__, sv(ks,i,j-1) )
1037  call check( __line__, sv(ks,i,j ) )
1038  call check( __line__, st(ks,i,j) )
1039  call check( __line__, pott(ks ,i,j) )
1040  call check( __line__, pott(ks+1,i,j) )
1041  call check( __line__, pott(ks,i ,j) )
1042  call check( __line__, pott(ks,i,j-2) )
1043  call check( __line__, pott(ks,i,j-1) )
1044  call check( __line__, pott(ks,i,j ) )
1045  call check( __line__, pott(ks,i,j+1) )
1046  call check( __line__, pott(ks,i,j+2) )
1047  call check( __line__, gsqrt(ks,i,j,i_uyz) )
1048  call check( __line__, gsqrt(ks,i,j,i_xvz) )
1049  call check( __line__, gsqrt(ks,i,j,i_xyz) )
1050  call check( __line__, dpres(ks ,i,j) )
1051  call check( __line__, dpres(ks+1,i,j) )
1052  call check( __line__, rt2p(ks ,i,j) )
1053  call check( __line__, rt2p(ks+1,i,j) )
1054  call check( __line__, rhot(ks+1,i,j) )
1055  call check( __line__, ddens(ks+1,i,j) )
1056 #endif
1057  b(ks,i,j) = ( &
1058  ( j33g * ( momz(ks ,i,j) + dtrk*sw(ks ,i,j) ) &
1059  * 0.5_rp*(pott(ks+1,i,j)+pott(ks ,i,j)) ) * rcdz(ks) &
1060  + ( gsqrt(ks,i,j ,i_xvz) * ( momy(ks,i,j ) + dtrk*sv(ks,i,j ) ) &
1061  * ( fact_n*(pott(ks,i,j+1)+pott(ks,i,j )) &
1062  + fact_f*(pott(ks,i,j+2)+pott(ks,i,j-1)) ) &
1063  - gsqrt(ks,i,j-1,i_xvz) * ( momy(ks,i,j-1) + dtrk*sv(ks,i,j-1) ) &
1064  * ( fact_n*(pott(ks,i,j )+pott(ks,i,j-1)) &
1065  + fact_f*(pott(ks,i,j+1)+pott(ks,i,j-2)) ) ) * rcdy(j) &
1066  + gsqrt(ks,i,j,i_xyz) * ( st(ks,i,j) - dpres(ks,i,j) * rcs2t(ks,i,j) * rdt ) &
1067  ) * rdt &
1068  + grav * j33g * 0.5_rp * ( ( dpres(ks,i,j)+dpres(ks+1,i,j) ) * rcs2t(ks,i,j) &
1069  - ( ddens(ks,i,j)+ddens(ks+1,i,j) ) * pott(ks,i,j) ) * rcdz(ks)
1070 #ifdef DEBUG
1071  call check( __line__, momz(ks ,i,j) )
1072  call check( __line__, momz(ks+1,i,j) )
1073  call check( __line__, momx(ks+1,i,j) )
1074  call check( __line__, momy(ks+1,i,j-1) )
1075  call check( __line__, momy(ks+1,i,j ) )
1076  call check( __line__, pott(ks ,i,j) )
1077  call check( __line__, pott(ks+1 ,i,j) )
1078  call check( __line__, pott(ks+1+1,i,j) )
1079  call check( __line__, pott(ks+1+2,i,j) )
1080  call check( __line__, pott(ks+1,i,j) )
1081  call check( __line__, pott(ks+1,i,j-2) )
1082  call check( __line__, pott(ks+1,i,j-1) )
1083  call check( __line__, pott(ks+1,i,j ) )
1084  call check( __line__, pott(ks+1,i,j+1) )
1085  call check( __line__, pott(ks+1,i,j+2) )
1086  call check( __line__, sw(ks ,i,j) )
1087  call check( __line__, sw(ks+1,i,j) )
1088  call check( __line__, su(ks+1,i,j) )
1089  call check( __line__, sv(ks+1,i,j-1) )
1090  call check( __line__, sv(ks+1,i,j ) )
1091  call check( __line__, gsqrt(ks+1,i,j,i_uyz) )
1092  call check( __line__, gsqrt(ks+1,i,j-1,i_xvz) )
1093  call check( __line__, gsqrt(ks+1,i,j ,i_xvz) )
1094  call check( __line__, gsqrt(ks+1,i,j,i_xyz) )
1095  call check( __line__, st(ks+1,i,j) )
1096  call check( __line__, dpres(ks ,i,j) )
1097  call check( __line__, dpres(ks+1,i,j) )
1098  call check( __line__, dpres(ks+2,i,j) )
1099  call check( __line__, rt2p(ks ,i,j) )
1100  call check( __line__, rt2p(ks+1,i,j) )
1101  call check( __line__, rt2p(ks+2,i,j) )
1102  call check( __line__, rhot(ks ,i,j) )
1103  call check( __line__, rhot(ks+2,i,j) )
1104  call check( __line__, ddens(ks ,i,j) )
1105  call check( __line__, ddens(ks+2,i,j) )
1106 #endif
1107  b(ks+1,i,j) = ( &
1108  ( j33g * ( momz(ks+1,i,j) + dtrk*sw(ks+1,i,j) ) &
1109  * ( fact_n*(pott(ks+2,i,j)+pott(ks+1,i,j)) &
1110  + fact_f*(pott(ks+3,i,j)+pott(ks ,i,j)) ) &
1111  - j33g * ( momz(ks+1-1,i,j) + dtrk*sw(ks+1-1,i,j) ) &
1112  * ( 0.5_rp*(pott(ks+1,i,j)+pott(ks,i,j)) ) ) * rcdz(ks+1) &
1113  + ( gsqrt(ks+1,i,j ,i_xvz) * ( momy(ks+1,i,j ) + dtrk*sv(ks+1,i,j ) ) &
1114  * ( fact_n*(pott(ks+1,i,j+1)+pott(ks+1,i,j )) &
1115  + fact_f*(pott(ks+1,i,j+2)+pott(ks+1,i,j-1)) ) &
1116  - gsqrt(ks+1,i,j-1,i_xvz) * ( momy(ks+1,i,j-1) + dtrk*sv(ks+1,i,j-1) ) &
1117  * ( fact_n*(pott(ks+1,i,j )+pott(ks+1,i,j-1)) &
1118  + fact_f*(pott(ks+1,i,j+1)+pott(ks+1,i,j-2)) ) ) * rcdy(j) &
1119  + gsqrt(ks+1,i,j,i_xyz) * ( st(ks+1,i,j) - dpres(ks+1,i,j) * rcs2t(ks+1,i,j) * rdt ) &
1120  ) * rdt &
1121  + grav * j33g * ( ( dpres(ks+2,i,j)*rcs2t(ks+2,i,j) &
1122  - dpres(ks ,i,j)*rcs2t(ks ,i,j) ) &
1123  - ( rhot(ks+2,i,j)*ddens(ks+2,i,j)/dens(ks+2,i,j) &
1124  - rhot(ks ,i,j)*ddens(ks ,i,j)/dens(ks ,i,j) ) &
1125  ) / ( fdz(ks+1) + fdz(ks) )
1126 #ifdef DEBUG
1127  call check( __line__, momz(ke-2,i,j) )
1128  call check( __line__, momz(ke-1,i,j) )
1129  call check( __line__, momx(ke-1,i,j) )
1130  call check( __line__, momy(ke-1,i,j-1) )
1131  call check( __line__, momy(ke-1,i,j ) )
1132  call check( __line__, pott(ke-3,i,j) )
1133  call check( __line__, pott(ke-2,i,j) )
1134  call check( __line__, pott(ke-1,i,j) )
1135  call check( __line__, pott(ke ,i,j) )
1136  call check( __line__, pott(ke-1,i ,j) )
1137  call check( __line__, pott(ke-1,i,j-2) )
1138  call check( __line__, pott(ke-1,i,j-1) )
1139  call check( __line__, pott(ke-1,i,j ) )
1140  call check( __line__, pott(ke-1,i,j+1) )
1141  call check( __line__, pott(ke-1,i,j+2) )
1142  call check( __line__, sw(ke-2,i,j) )
1143  call check( __line__, sw(ke-1,i,j) )
1144  call check( __line__, su(ke-1,i,j) )
1145  call check( __line__, sv(ke-1,i,j-1) )
1146  call check( __line__, sv(ke-1,i,j ) )
1147  call check( __line__, gsqrt(ke-1,i,j,i_uyz) )
1148  call check( __line__, gsqrt(ke-1,i,j-1,i_xvz) )
1149  call check( __line__, gsqrt(ke-1,i,j ,i_xvz) )
1150  call check( __line__, gsqrt(ke-1,i,j,i_xyz) )
1151  call check( __line__, st(ke-1,i,j) )
1152  call check( __line__, dpres(ke-2,i,j) )
1153  call check( __line__, dpres(ke-1,i,j) )
1154  call check( __line__, dpres(ke ,i,j) )
1155  call check( __line__, rt2p(ke-2,i,j) )
1156  call check( __line__, rt2p(ke-1,i,j) )
1157  call check( __line__, rt2p(ke ,i,j) )
1158  call check( __line__, rhot(ke-2,i,j) )
1159  call check( __line__, rhot(ke,i,j) )
1160  call check( __line__, ddens(ke-2,i,j) )
1161  call check( __line__, ddens(ke,i,j) )
1162 #endif
1163  b(ke-1,i,j) = ( &
1164  ( j33g * ( momz(ke-1,i,j) + dtrk*sw(ke-1,i,j) ) &
1165  * ( 0.5_rp*(pott(ke ,i,j)+pott(ke-1,i,j)) ) &
1166  - j33g * ( momz(ke-2,i,j) + dtrk*sw(ke-2,i,j) ) &
1167  * ( fact_n*(pott(ke-1,i,j)+pott(ke-2,i,j)) &
1168  + fact_f*(pott(ke ,i,j)+pott(ke-3,i,j)) ) ) * rcdz(ke-1) &
1169  + ( gsqrt(ke-1,i,j ,i_xvz) * ( momy(ke-1,i,j ) + dtrk*sv(ke-1,i,j ) ) &
1170  * ( fact_n*(pott(ke-1,i,j+1)+pott(ke-1,i,j )) &
1171  + fact_f*(pott(ke-1,i,j+2)+pott(ke-1,i,j-1)) ) &
1172  - gsqrt(ke-1,i,j-1,i_xvz) * ( momy(ke-1,i,j-1) + dtrk*sv(ke-1,i,j-1) ) &
1173  * ( fact_n*(pott(ke-1,i,j )+pott(ke-1,i,j-1)) &
1174  + fact_f*(pott(ke-1,i,j+1)+pott(ke-1,i,j-2)) ) ) * rcdy(j) &
1175  + gsqrt(ke-1,i,j,i_xyz) * ( st(ke-1,i,j) - dpres(ke-1,i,j) * rcs2t(ke-1,i,j) * rdt ) &
1176  ) * rdt &
1177  + grav * j33g * ( ( dpres(ke ,i,j)*rcs2t(ke ,i,j) &
1178  - dpres(ke-2,i,j)*rcs2t(ke-2,i,j) ) &
1179  - ( rhot(ke ,i,j)*ddens(ke ,i,j)/dens(ke ,i,j) &
1180  - rhot(ke-2,i,j)*ddens(ke-2,i,j)/dens(ke-2,i,j) )&
1181  ) / ( fdz(ke-1) + fdz(ke-1-1) )
1182 #ifdef DEBUG
1183  call check( __line__, momz(ke-1,i,j) )
1184  call check( __line__, momx(ke,i,j) )
1185  call check( __line__, momy(ke,i,j-1) )
1186  call check( __line__, momy(ke,i,j ) )
1187  call check( __line__, sw(ke-1,i,j) )
1188  call check( __line__, su(ke,i,j) )
1189  call check( __line__, sv(ke,i,j-1) )
1190  call check( __line__, sv(ke,i,j ) )
1191  call check( __line__, pott(ke-1,i,j) )
1192  call check( __line__, pott(ke ,i,j) )
1193  call check( __line__, pott(ke,i,j) )
1194  call check( __line__, pott(ke,i,j-2) )
1195  call check( __line__, pott(ke,i,j-1) )
1196  call check( __line__, pott(ke,i,j ) )
1197  call check( __line__, pott(ke,i,j+1) )
1198  call check( __line__, pott(ke,i,j+2) )
1199  call check( __line__, gsqrt(ke,i,j,i_uyz) )
1200  call check( __line__, gsqrt(ke,i,j-1,i_xvz) )
1201  call check( __line__, gsqrt(ke,i,j ,i_xvz) )
1202  call check( __line__, gsqrt(ke,i,j,i_xyz) )
1203  call check( __line__, st(ke,i,j) )
1204  call check( __line__, dpres(ke-1,i,j) )
1205  call check( __line__, dpres(ke ,i,j) )
1206  call check( __line__, rt2p(ke-1,i,j) )
1207  call check( __line__, rt2p(ke ,i,j) )
1208  call check( __line__, rhot(ke-1,i,j) )
1209  call check( __line__, ddens(ke-1,i,j) )
1210 #endif
1211  b(ke,i,j) = ( &
1212  ( &
1213  - j33g * ( momz(ke-1,i,j) + dtrk*sw(ke-1,i,j) ) &
1214  * 0.5_rp*(pott(ke ,i,j)+pott(ke-1,i,j)) ) * rcdz(ke) &
1215  + ( gsqrt(ke,i,j ,i_xvz) * ( momy(ke,i,j ) + dtrk*sv(ke,i,j ) ) &
1216  * ( fact_n*(pott(ke,i,j+1)+pott(ke,i,j )) &
1217  + fact_f*(pott(ke,i,j+2)+pott(ke,i,j-1)) ) &
1218  - gsqrt(ke,i,j-1,i_xvz) * ( momy(ke,i,j-1) + dtrk*sv(ke,i,j-1) ) &
1219  * ( fact_n*(pott(ke,i,j )+pott(ke,i,j-1)) &
1220  + fact_f*(pott(ke,i,j+1)+pott(ke,i,j-2)) ) ) * rcdy(j) &
1221  + gsqrt(ke,i,j,i_xyz) * ( st(ke,i,j) - dpres(ke,i,j) * rcs2t(ke,i,j) * rdt ) &
1222  ) * rdt &
1223  + grav * j33g * 0.5_rp * ( - ( dpres(ke,i,j)+dpres(ke-1,i,j) ) * rcs2t(ke,i,j) &
1224  + ( ddens(ke,i,j)+ddens(ke-1,i,j) ) * pott(ke,i,j) ) * rcdz(ke)
1225  enddo
1226  enddo
1227 #ifdef DEBUG
1228  k = iundef; i = iundef; j = iundef
1229 #endif
1230 
1231  else
1232 
1233 
1234  ! r = b - M x0
1235  do j = jjs, jje
1236  do i = iis, iie
1237  do k = ks+2, ke-2
1238 #ifdef DEBUG
1239  call check( __line__, momz(k-1,i,j) )
1240  call check( __line__, momz(k ,i,j) )
1241  call check( __line__, momx(k,i-1,j) )
1242  call check( __line__, momx(k,i ,j) )
1243  call check( __line__, momy(k,i,j-1) )
1244  call check( __line__, momy(k,i,j ) )
1245  call check( __line__, pott(k-2,i,j) )
1246  call check( __line__, pott(k-1,i,j) )
1247  call check( __line__, pott(k ,i,j) )
1248  call check( __line__, pott(k+1,i,j) )
1249  call check( __line__, pott(k+2,i,j) )
1250  call check( __line__, pott(k,i-2,j) )
1251  call check( __line__, pott(k,i-1,j) )
1252  call check( __line__, pott(k,i ,j) )
1253  call check( __line__, pott(k,i+1,j) )
1254  call check( __line__, pott(k,i+2,j) )
1255  call check( __line__, pott(k,i,j-2) )
1256  call check( __line__, pott(k,i,j-1) )
1257  call check( __line__, pott(k,i,j ) )
1258  call check( __line__, pott(k,i,j+1) )
1259  call check( __line__, pott(k,i,j+2) )
1260  call check( __line__, sw(k-1,i,j) )
1261  call check( __line__, sw(k ,i,j) )
1262  call check( __line__, su(k,i-1,j) )
1263  call check( __line__, su(k,i ,j) )
1264  call check( __line__, sv(k,i,j-1) )
1265  call check( __line__, sv(k,i,j ) )
1266  call check( __line__, gsqrt(k,i ,j,i_uyz) )
1267  call check( __line__, gsqrt(k,i-1,j,i_uyz) )
1268  call check( __line__, gsqrt(k,i,j-1,i_xvz) )
1269  call check( __line__, gsqrt(k,i,j ,i_xvz) )
1270  call check( __line__, gsqrt(k,i,j,i_xyz) )
1271  call check( __line__, st(k,i,j) )
1272  call check( __line__, dpres(k-1,i,j) )
1273  call check( __line__, dpres(k ,i,j) )
1274  call check( __line__, dpres(k+1,i,j) )
1275  call check( __line__, rt2p(k-1,i,j) )
1276  call check( __line__, rt2p(k ,i,j) )
1277  call check( __line__, rt2p(k+1,i,j) )
1278  call check( __line__, rhot(k-1,i,j) )
1279  call check( __line__, rhot(k+1,i,j) )
1280  call check( __line__, ddens(k-1,i,j) )
1281  call check( __line__, ddens(k+1,i,j) )
1282 #endif
1283  b(k,i,j) = ( &
1284  ( j33g * ( momz(k ,i,j) + dtrk*sw(k ,i,j) ) &
1285  * ( fact_n*(pott(k+1,i,j)+pott(k ,i,j)) &
1286  + fact_f*(pott(k+2,i,j)+pott(k-1,i,j)) ) &
1287  - j33g * ( momz(k-1,i,j) + dtrk*sw(k-1,i,j) ) &
1288  * ( fact_n*(pott(k ,i,j)+pott(k-1,i,j)) &
1289  + fact_f*(pott(k+1,i,j)+pott(k-2,i,j)) ) ) * rcdz(k) &
1290  + ( gsqrt(k,i ,j,i_uyz) * ( momx(k,i ,j) + dtrk*su(k,i ,j) ) &
1291  * ( fact_n*(pott(k,i+1,j)+pott(k,i ,j)) &
1292  + fact_f*(pott(k,i+2,j)+pott(k,i-1,j)) ) &
1293  - gsqrt(k,i-1,j,i_uyz) * ( momx(k,i-1,j) + dtrk*su(k,i-1,j) ) &
1294  * ( fact_n*(pott(k,i ,j)+pott(k,i-1,j)) &
1295  + fact_f*(pott(k,i+1,j)+pott(k,i-2,j)) ) ) * rcdx(i) &
1296  + ( gsqrt(k,i,j ,i_xvz) * ( momy(k,i,j ) + dtrk*sv(k,i,j ) ) &
1297  * ( fact_n*(pott(k,i,j+1)+pott(k,i,j )) &
1298  + fact_f*(pott(k,i,j+2)+pott(k,i,j-1)) ) &
1299  - gsqrt(k,i,j-1,i_xvz) * ( momy(k,i,j-1) + dtrk*sv(k,i,j-1) ) &
1300  * ( fact_n*(pott(k,i,j )+pott(k,i,j-1)) &
1301  + fact_f*(pott(k,i,j+1)+pott(k,i,j-2)) ) ) * rcdy(j) &
1302  + gsqrt(k,i,j,i_xyz) * ( st(k,i,j) - dpres(k,i,j) * rcs2t(k,i,j) * rdt ) &
1303  ) * rdt &
1304  + grav * j33g * ( ( dpres(k+1,i,j)*rcs2t(k+1,i,j) &
1305  - dpres(k-1,i,j)*rcs2t(k-1,i,j) ) &
1306  - ( rhot(k+1,i,j)*ddens(k+1,i,j)/dens(k+1,i,j) &
1307  - rhot(k-1,i,j)*ddens(k-1,i,j)/dens(k-1,i,j) ) &
1308  ) / ( fdz(k) + fdz(k-1) )
1309  enddo
1310  enddo
1311  enddo
1312 #ifdef DEBUG
1313  k = iundef; i = iundef; j = iundef
1314 #endif
1315  do j = jjs, jje
1316  do i = iis, iie
1317 #ifdef DEBUG
1318  call check( __line__, momz(ks,i,j) )
1319  call check( __line__, momx(ks,i,j) )
1320  call check( __line__, momy(ks,i,j) )
1321  call check( __line__, sw(ks,i,j) )
1322  call check( __line__, su(ks,i-1,j) )
1323  call check( __line__, su(ks,i ,j) )
1324  call check( __line__, sv(ks,i,j-1) )
1325  call check( __line__, sv(ks,i,j ) )
1326  call check( __line__, st(ks,i,j) )
1327  call check( __line__, pott(ks ,i,j) )
1328  call check( __line__, pott(ks+1,i,j) )
1329  call check( __line__, pott(ks,i-2,j) )
1330  call check( __line__, pott(ks,i-1,j) )
1331  call check( __line__, pott(ks,i ,j) )
1332  call check( __line__, pott(ks,i+1,j) )
1333  call check( __line__, pott(ks,i+2,j) )
1334  call check( __line__, pott(ks,i,j-2) )
1335  call check( __line__, pott(ks,i,j-1) )
1336  call check( __line__, pott(ks,i,j ) )
1337  call check( __line__, pott(ks,i,j+1) )
1338  call check( __line__, pott(ks,i,j+2) )
1339  call check( __line__, gsqrt(ks,i,j,i_uyz) )
1340  call check( __line__, gsqrt(ks,i,j,i_xvz) )
1341  call check( __line__, gsqrt(ks,i,j,i_xyz) )
1342  call check( __line__, dpres(ks ,i,j) )
1343  call check( __line__, dpres(ks+1,i,j) )
1344  call check( __line__, rt2p(ks ,i,j) )
1345  call check( __line__, rt2p(ks+1,i,j) )
1346  call check( __line__, rhot(ks+1,i,j) )
1347  call check( __line__, ddens(ks+1,i,j) )
1348 #endif
1349  b(ks,i,j) = ( &
1350  ( j33g * ( momz(ks ,i,j) + dtrk*sw(ks ,i,j) ) &
1351  * 0.5_rp*(pott(ks+1,i,j)+pott(ks ,i,j)) ) * rcdz(ks) &
1352  + ( gsqrt(ks,i ,j,i_uyz) * ( momx(ks,i ,j) + dtrk*su(ks,i ,j) ) &
1353  * ( fact_n*(pott(ks,i+1,j)+pott(ks,i ,j)) &
1354  + fact_f*(pott(ks,i+2,j)+pott(ks,i-1,j)) ) &
1355  - gsqrt(ks,i-1,j,i_uyz) * ( momx(ks,i-1,j) + dtrk*su(ks,i-1,j) ) &
1356  * ( fact_n*(pott(ks,i ,j)+pott(ks,i-1,j)) &
1357  + fact_f*(pott(ks,i+1,j)+pott(ks,i-2,j)) ) ) * rcdx(i) &
1358  + ( gsqrt(ks,i,j ,i_xvz) * ( momy(ks,i,j ) + dtrk*sv(ks,i,j ) ) &
1359  * ( fact_n*(pott(ks,i,j+1)+pott(ks,i,j )) &
1360  + fact_f*(pott(ks,i,j+2)+pott(ks,i,j-1)) ) &
1361  - gsqrt(ks,i,j-1,i_xvz) * ( momy(ks,i,j-1) + dtrk*sv(ks,i,j-1) ) &
1362  * ( fact_n*(pott(ks,i,j )+pott(ks,i,j-1)) &
1363  + fact_f*(pott(ks,i,j+1)+pott(ks,i,j-2)) ) ) * rcdy(j) &
1364  + gsqrt(ks,i,j,i_xyz) * ( st(ks,i,j) - dpres(ks,i,j) * rcs2t(ks,i,j) * rdt ) &
1365  ) * rdt &
1366  + grav * j33g * 0.5_rp * ( ( dpres(ks,i,j)+dpres(ks+1,i,j) ) * rcs2t(ks,i,j) &
1367  - ( ddens(ks,i,j)+ddens(ks+1,i,j) ) * pott(ks,i,j) ) * rcdz(ks)
1368 #ifdef DEBUG
1369  call check( __line__, momz(ks ,i,j) )
1370  call check( __line__, momz(ks+1 ,i,j) )
1371  call check( __line__, momx(ks+1,i-1,j) )
1372  call check( __line__, momx(ks+1,i ,j) )
1373  call check( __line__, momy(ks+1,i,j-1) )
1374  call check( __line__, momy(ks+1,i,j ) )
1375  call check( __line__, pott(ks ,i,j) )
1376  call check( __line__, pott(ks+1 ,i,j) )
1377  call check( __line__, pott(ks+1+1,i,j) )
1378  call check( __line__, pott(ks+1+2,i,j) )
1379  call check( __line__, pott(ks+1,i-2,j) )
1380  call check( __line__, pott(ks+1,i-1,j) )
1381  call check( __line__, pott(ks+1,i ,j) )
1382  call check( __line__, pott(ks+1,i+1,j) )
1383  call check( __line__, pott(ks+1,i+2,j) )
1384  call check( __line__, pott(ks+1,i,j-2) )
1385  call check( __line__, pott(ks+1,i,j-1) )
1386  call check( __line__, pott(ks+1,i,j ) )
1387  call check( __line__, pott(ks+1,i,j+1) )
1388  call check( __line__, pott(ks+1,i,j+2) )
1389  call check( __line__, sw(ks ,i,j) )
1390  call check( __line__, sw(ks+1,i,j) )
1391  call check( __line__, su(ks+1,i-1,j) )
1392  call check( __line__, su(ks+1,i ,j) )
1393  call check( __line__, sv(ks+1,i,j-1) )
1394  call check( __line__, sv(ks+1,i,j ) )
1395  call check( __line__, gsqrt(ks+1,i ,j,i_uyz) )
1396  call check( __line__, gsqrt(ks+1,i-1,j,i_uyz) )
1397  call check( __line__, gsqrt(ks+1,i,j-1,i_xvz) )
1398  call check( __line__, gsqrt(ks+1,i,j ,i_xvz) )
1399  call check( __line__, gsqrt(ks+1,i,j,i_xyz) )
1400  call check( __line__, st(ks+1,i,j) )
1401  call check( __line__, dpres(ks ,i,j) )
1402  call check( __line__, dpres(ks+1,i,j) )
1403  call check( __line__, dpres(ks+2,i,j) )
1404  call check( __line__, rt2p(ks ,i,j) )
1405  call check( __line__, rt2p(ks+1,i,j) )
1406  call check( __line__, rt2p(ks+2,i,j) )
1407  call check( __line__, rhot(ks ,i,j) )
1408  call check( __line__, rhot(ks+2,i,j) )
1409  call check( __line__, ddens(ks ,i,j) )
1410  call check( __line__, ddens(ks+2,i,j) )
1411 #endif
1412  b(ks+1,i,j) = ( &
1413  ( j33g * ( momz(ks+1,i,j) + dtrk*sw(ks+1,i,j) ) &
1414  * ( fact_n*(pott(ks+2,i,j)+pott(ks+1,i,j)) &
1415  + fact_f*(pott(ks+3,i,j)+pott(ks ,i,j)) ) &
1416  - j33g * ( momz(ks+1-1,i,j) + dtrk*sw(ks+1-1,i,j) ) &
1417  * ( 0.5_rp*(pott(ks+1,i,j)+pott(ks,i,j)) ) ) * rcdz(ks+1) &
1418  + ( gsqrt(ks+1,i ,j,i_uyz) * ( momx(ks+1,i ,j) + dtrk*su(ks+1,i ,j) ) &
1419  * ( fact_n*(pott(ks+1,i+1,j)+pott(ks+1,i ,j)) &
1420  + fact_f*(pott(ks+1,i+2,j)+pott(ks+1,i-1,j)) ) &
1421  - gsqrt(ks+1,i-1,j,i_uyz) * ( momx(ks+1,i-1,j) + dtrk*su(ks+1,i-1,j) ) &
1422  * ( fact_n*(pott(ks+1,i ,j)+pott(ks+1,i-1,j)) &
1423  + fact_f*(pott(ks+1,i+1,j)+pott(ks+1,i-2,j)) ) ) * rcdx(i) &
1424  + ( gsqrt(ks+1,i,j ,i_xvz) * ( momy(ks+1,i,j ) + dtrk*sv(ks+1,i,j ) ) &
1425  * ( fact_n*(pott(ks+1,i,j+1)+pott(ks+1,i,j )) &
1426  + fact_f*(pott(ks+1,i,j+2)+pott(ks+1,i,j-1)) ) &
1427  - gsqrt(ks+1,i,j-1,i_xvz) * ( momy(ks+1,i,j-1) + dtrk*sv(ks+1,i,j-1) ) &
1428  * ( fact_n*(pott(ks+1,i,j )+pott(ks+1,i,j-1)) &
1429  + fact_f*(pott(ks+1,i,j+1)+pott(ks+1,i,j-2)) ) ) * rcdy(j) &
1430  + gsqrt(ks+1,i,j,i_xyz) * ( st(ks+1,i,j) - dpres(ks+1,i,j) * rcs2t(ks+1,i,j) * rdt ) &
1431  ) * rdt &
1432  + grav * j33g * ( ( dpres(ks+2,i,j)*rcs2t(ks+2,i,j) &
1433  - dpres(ks ,i,j)*rcs2t(ks ,i,j) ) &
1434  - ( rhot(ks+2,i,j)*ddens(ks+2,i,j)/dens(ks+2,i,j) &
1435  - rhot(ks ,i,j)*ddens(ks ,i,j)/dens(ks ,i,j) ) &
1436  ) / ( fdz(ks+1) + fdz(ks) )
1437 #ifdef DEBUG
1438  call check( __line__, momz(ke-2,i,j) )
1439  call check( __line__, momz(ke-1,i,j) )
1440  call check( __line__, momx(ke-1,i-1,j) )
1441  call check( __line__, momx(ke-1,i ,j) )
1442  call check( __line__, momy(ke-1,i,j-1) )
1443  call check( __line__, momy(ke-1,i,j ) )
1444  call check( __line__, pott(ke-3,i,j) )
1445  call check( __line__, pott(ke-2,i,j) )
1446  call check( __line__, pott(ke-1,i,j) )
1447  call check( __line__, pott(ke ,i,j) )
1448  call check( __line__, pott(ke-1,i-2,j) )
1449  call check( __line__, pott(ke-1,i-1,j) )
1450  call check( __line__, pott(ke-1,i ,j) )
1451  call check( __line__, pott(ke-1,i+1,j) )
1452  call check( __line__, pott(ke-1,i+2,j) )
1453  call check( __line__, pott(ke-1,i,j-2) )
1454  call check( __line__, pott(ke-1,i,j-1) )
1455  call check( __line__, pott(ke-1,i,j ) )
1456  call check( __line__, pott(ke-1,i,j+1) )
1457  call check( __line__, pott(ke-1,i,j+2) )
1458  call check( __line__, sw(ke-2,i,j) )
1459  call check( __line__, sw(ke-1,i,j) )
1460  call check( __line__, su(ke-1,i-1,j) )
1461  call check( __line__, su(ke-1,i ,j) )
1462  call check( __line__, sv(ke-1,i,j-1) )
1463  call check( __line__, sv(ke-1,i,j ) )
1464  call check( __line__, gsqrt(ke-1,i ,j,i_uyz) )
1465  call check( __line__, gsqrt(ke-1,i-1,j,i_uyz) )
1466  call check( __line__, gsqrt(ke-1,i,j-1,i_xvz) )
1467  call check( __line__, gsqrt(ke-1,i,j ,i_xvz) )
1468  call check( __line__, gsqrt(ke-1,i,j,i_xyz) )
1469  call check( __line__, st(ke-1,i,j) )
1470  call check( __line__, dpres(ke-2,i,j) )
1471  call check( __line__, dpres(ke-1,i,j) )
1472  call check( __line__, dpres(ke ,i,j) )
1473  call check( __line__, rt2p(ke-2,i,j) )
1474  call check( __line__, rt2p(ke-1,i,j) )
1475  call check( __line__, rt2p(ke ,i,j) )
1476  call check( __line__, rhot(ke-2,i,j) )
1477  call check( __line__, rhot(ke,i,j) )
1478  call check( __line__, ddens(ke-2,i,j) )
1479  call check( __line__, ddens(ke,i,j) )
1480 #endif
1481  b(ke-1,i,j) = ( &
1482  ( j33g * ( momz(ke-1,i,j) + dtrk*sw(ke-1,i,j) ) &
1483  * ( 0.5_rp*(pott(ke ,i,j)+pott(ke-1,i,j)) ) &
1484  - j33g * ( momz(ke-2,i,j) + dtrk*sw(ke-2,i,j) ) &
1485  * ( fact_n*(pott(ke-1,i,j)+pott(ke-2,i,j)) &
1486  + fact_f*(pott(ke ,i,j)+pott(ke-3,i,j)) ) ) * rcdz(ke-1) &
1487  + ( gsqrt(ke-1,i ,j,i_uyz) * ( momx(ke-1,i ,j) + dtrk*su(ke-1,i ,j) ) &
1488  * ( fact_n*(pott(ke-1,i+1,j)+pott(ke-1,i ,j)) &
1489  + fact_f*(pott(ke-1,i+2,j)+pott(ke-1,i-1,j)) ) &
1490  - gsqrt(ke-1,i-1,j,i_uyz) * ( momx(ke-1,i-1,j) + dtrk*su(ke-1,i-1,j) ) &
1491  * ( fact_n*(pott(ke-1,i ,j)+pott(ke-1,i-1,j)) &
1492  + fact_f*(pott(ke-1,i+1,j)+pott(ke-1,i-2,j)) ) ) * rcdx(i) &
1493  + ( gsqrt(ke-1,i,j ,i_xvz) * ( momy(ke-1,i,j ) + dtrk*sv(ke-1,i,j ) ) &
1494  * ( fact_n*(pott(ke-1,i,j+1)+pott(ke-1,i,j )) &
1495  + fact_f*(pott(ke-1,i,j+2)+pott(ke-1,i,j-1)) ) &
1496  - gsqrt(ke-1,i,j-1,i_xvz) * ( momy(ke-1,i,j-1) + dtrk*sv(ke-1,i,j-1) ) &
1497  * ( fact_n*(pott(ke-1,i,j )+pott(ke-1,i,j-1)) &
1498  + fact_f*(pott(ke-1,i,j+1)+pott(ke-1,i,j-2)) ) ) * rcdy(j) &
1499  + gsqrt(ke-1,i,j,i_xyz) * ( st(ke-1,i,j) - dpres(ke-1,i,j) * rcs2t(ke-1,i,j) * rdt ) &
1500  ) * rdt &
1501  + grav * j33g * ( ( dpres(ke ,i,j)*rcs2t(ke ,i,j) &
1502  - dpres(ke-2,i,j)*rcs2t(ke-2,i,j) ) &
1503  - ( rhot(ke ,i,j)*ddens(ke ,i,j)/dens(ke ,i,j) &
1504  - rhot(ke-2,i,j)*ddens(ke-2,i,j)/dens(ke-2,i,j) )&
1505  ) / ( fdz(ke-1) + fdz(ke-1-1) )
1506 #ifdef DEBUG
1507  call check( __line__, momz(ke-1,i,j) )
1508  call check( __line__, momx(ke,i-1,j) )
1509  call check( __line__, momx(ke,i ,j) )
1510  call check( __line__, momy(ke,i,j-1) )
1511  call check( __line__, momy(ke,i,j ) )
1512  call check( __line__, sw(ke-1,i,j) )
1513  call check( __line__, su(ke,i-1,j) )
1514  call check( __line__, su(ke,i ,j) )
1515  call check( __line__, sv(ke,i,j-1) )
1516  call check( __line__, sv(ke,i,j ) )
1517  call check( __line__, gsqrt(ke,i-1,j,i_uyz) )
1518  call check( __line__, pott(ke-1,i,j) )
1519  call check( __line__, pott(ke ,i,j) )
1520  call check( __line__, pott(ke,i-2,j) )
1521  call check( __line__, pott(ke,i-1,j) )
1522  call check( __line__, pott(ke,i ,j) )
1523  call check( __line__, pott(ke,i+1,j) )
1524  call check( __line__, pott(ke,i+2,j) )
1525  call check( __line__, pott(ke,i,j-2) )
1526  call check( __line__, pott(ke,i,j-1) )
1527  call check( __line__, pott(ke,i,j ) )
1528  call check( __line__, pott(ke,i,j+1) )
1529  call check( __line__, pott(ke,i,j+2) )
1530  call check( __line__, gsqrt(ke,i-1,j,i_uyz) )
1531  call check( __line__, gsqrt(ke,i ,j,i_uyz) )
1532  call check( __line__, gsqrt(ke,i,j-1,i_xvz) )
1533  call check( __line__, gsqrt(ke,i,j ,i_xvz) )
1534  call check( __line__, gsqrt(ke,i,j,i_xyz) )
1535  call check( __line__, st(ke,i,j) )
1536  call check( __line__, dpres(ke-1,i,j) )
1537  call check( __line__, dpres(ke ,i,j) )
1538  call check( __line__, rt2p(ke-1,i,j) )
1539  call check( __line__, rt2p(ke ,i,j) )
1540  call check( __line__, rhot(ke-1,i,j) )
1541  call check( __line__, ddens(ke-1,i,j) )
1542 #endif
1543  b(ke,i,j) = ( &
1544  ( &
1545  - j33g * ( momz(ke-1,i,j) + dtrk*sw(ke-1,i,j) ) &
1546  * 0.5_rp*(pott(ke ,i,j)+pott(ke-1,i,j)) ) * rcdz(ke) &
1547  + ( gsqrt(ke,i ,j,i_uyz) * ( momx(ke,i ,j) + dtrk*su(ke,i ,j) ) &
1548  * ( fact_n*(pott(ke,i+1,j)+pott(ke,i ,j)) &
1549  + fact_f*(pott(ke,i+2,j)+pott(ke,i-1,j)) ) &
1550  - gsqrt(ke,i-1,j,i_uyz) * ( momx(ke,i-1,j) + dtrk*su(ke,i-1,j) ) &
1551  * ( fact_n*(pott(ke,i ,j)+pott(ke,i-1,j)) &
1552  + fact_f*(pott(ke,i+1,j)+pott(ke,i-2,j)) ) ) * rcdx(i) &
1553  + ( gsqrt(ke,i,j ,i_xvz) * ( momy(ke,i,j ) + dtrk*sv(ke,i,j ) ) &
1554  * ( fact_n*(pott(ke,i,j+1)+pott(ke,i,j )) &
1555  + fact_f*(pott(ke,i,j+2)+pott(ke,i,j-1)) ) &
1556  - gsqrt(ke,i,j-1,i_xvz) * ( momy(ke,i,j-1) + dtrk*sv(ke,i,j-1) ) &
1557  * ( fact_n*(pott(ke,i,j )+pott(ke,i,j-1)) &
1558  + fact_f*(pott(ke,i,j+1)+pott(ke,i,j-2)) ) ) * rcdy(j) &
1559  + gsqrt(ke,i,j,i_xyz) * ( st(ke,i,j) - dpres(ke,i,j) * rcs2t(ke,i,j) * rdt ) &
1560  ) * rdt &
1561  + grav * j33g * 0.5_rp * ( - ( dpres(ke,i,j)+dpres(ke-1,i,j) ) * rcs2t(ke,i,j) &
1562  + ( ddens(ke,i,j)+ddens(ke-1,i,j) ) * pott(ke,i,j) ) * rcdz(ke)
1563  enddo
1564  enddo
1565 #ifdef DEBUG
1566  k = iundef; i = iundef; j = iundef
1567 #endif
1568 
1569  end if
1570 
1571  call make_matrix( &
1572  m, & ! (inout)
1573  pott, rcs2t, grav, &
1574  gsqrt, j33g, &
1575  rcdz, rfdz, rcdx, rfdx, rcdy,rfdy, fdz, &
1576  rdt, &
1577  fact_n, fact_f, &
1578  twod, &
1579  i_xyz, i_xyw, &
1580  iis, iie, jjs, jje )
1581 
1582  enddo
1583  enddo
1584 
1585  call solve_bicgstab( &
1586  dpres_n, & ! (out)
1587  dpres, & ! (in)
1588  m, b, & ! (in)
1589  twod ) ! (in)
1590 
1591  call comm_vars8( dpres_n, 1 )
1592  call comm_wait ( dpres_n, 1 )
1593 
1594 #ifdef DEBUG
1595  call check_solver( dpres_n, m, b, twod )
1596 #endif
1597 
1598 #else
1599 
1600  call solve_multigrid
1601 
1602 #endif
1603 
1604  do jjs = js, je, jblock
1605  jje = jjs+jblock-1
1606  do iis = is, ie, iblock
1607  iie = iis+iblock-1
1608 
1609  !##### momentum equation (z) #####
1610 
1611  !--- update momentum(z)
1612  !$omp parallel do private(i,j,k,duvw) OMP_SCHEDULE_ collapse(2)
1613  do j = jjs, jje
1614  do i = iis, iie
1615  do k = ks, ke-1
1616 #ifdef DEBUG
1617  call check( __line__, dpres_n(k+1,i,j) )
1618  call check( __line__, dpres_n(k ,i,j) )
1619  call check( __line__, dpres(k+1,i,j) )
1620  call check( __line__, dpres(k ,i,j) )
1621  call check( __line__, ddens(k+1,i,j) )
1622  call check( __line__, ddens(k ,i,j) )
1623  call check( __line__, ref_dens(k+1,i,j) )
1624  call check( __line__, ref_dens(k,i,j) )
1625  call check( __line__, momz0(k,i,j) )
1626 #endif
1627  duvw = dtrk * ( ( &
1628  - j33g * ( dpres_n(k+1,i,j) - dpres_n(k,i,j) ) * rfdz(k) &
1629  ) / gsqrt(k,i,j,i_uyz) &
1630  - 0.5_rp * grav &
1631  * ( ddens(k+1,i,j) &
1632  + ( dpres_n(k+1,i,j) - dpres(k+1,i,j) ) &
1633  / ( pott(k+1,i,j) * rt2p(k+1,i,j) ) &
1634  + ddens(k ,i,j) &
1635  + ( dpres_n(k ,i,j) - dpres(k ,i,j) ) &
1636  / ( pott(k ,i,j) * rt2p(k ,i,j) ) ) &
1637  + sw(k,i,j) )
1638  momz_rk(k,i,j) = momz0(k,i,j) + duvw
1639  mflx_hi(k,i,j,zdir) = j33g * ( momz(k,i,j) + duvw )
1640  enddo
1641  enddo
1642  enddo
1643 #ifdef DEBUG
1644  k = iundef; i = iundef; j = iundef
1645 #endif
1646  do j = jjs, jje
1647  do i = iis, iie
1648  mflx_hi(ks-1,i,j,zdir) = 0.0_rp
1649  mflx_hi(ke ,i,j,zdir) = 0.0_rp
1650  enddo
1651  enddo
1652 #ifdef DEBUG
1653  k = iundef; i = iundef; j = iundef
1654 #endif
1655 
1656 
1657  !##### momentum equation (x) #####
1658 
1659  !--- update momentum(x)
1660  if ( twod ) then
1661  !$omp parallel do private(j,k,duvw) OMP_SCHEDULE_
1662  do j = jjs , jje
1663  do k = ks, ke
1664 #ifdef DEBUG
1665  call check( __line__, su(k,is,j) )
1666  call check( __line__, momx0(k,is,j) )
1667 #endif
1668  duvw = dtrk * su(k,is,j)
1669 
1670  momx_rk(k,is,j) = momx0(k,is,j) + duvw
1671  enddo
1672  enddo
1673  else
1674  !$omp parallel do private(i,j,k,duvw) OMP_SCHEDULE_ collapse(2)
1675  do j = jjs , jje
1676  do i = iis-1, iie
1677  do k = ks, ke
1678 #ifdef DEBUG
1679  call check( __line__, dpres_n(k,i+1,j) )
1680  call check( __line__, dpres_n(k,i ,j) )
1681  call check( __line__, gsqrt(k,i+1,j,i_xyz) )
1682  call check( __line__, gsqrt(k,i ,j,i_xyz) )
1683  call check( __line__, gsqrt(k,i,j,i_uyz) )
1684  call check( __line__, su(k,i,j) )
1685  call check( __line__, momx0(k,i,j) )
1686 #endif
1687  duvw = dtrk * ( ( &
1688  - ( gsqrt(k,i+1,j,i_xyz) * dpres_n(k,i+1,j) &
1689  - gsqrt(k,i ,j,i_xyz) * dpres_n(k,i ,j) ) * rfdx(i) & ! pressure gradient force
1690  ) / gsqrt(k,i,j,i_uyz) &
1691  + su(k,i,j) )
1692 
1693  momx_rk(k,i,j) = momx0(k,i,j) + duvw
1694  mflx_hi(k,i,j,xdir) = gsqrt(k,i,j,i_uyz) * ( momx(k,i,j) + duvw )
1695  enddo
1696  enddo
1697  enddo
1698  end if
1699 #ifdef DEBUG
1700  k = iundef; i = iundef; j = iundef
1701 #endif
1702 
1703  !##### momentum equation (y) #####
1704 
1705  !--- update momentum(y)
1706  !$omp parallel do private(i,j,k,duvw) OMP_SCHEDULE_ collapse(2)
1707  do j = jjs-1, jje
1708  do i = iis, iie
1709  do k = ks, ke
1710 #ifdef DEBUG
1711  call check( __line__, dpres_n(k,i,j ) )
1712  call check( __line__, dpres_n(k,i,j+1) )
1713  call check( __line__, gsqrt(k,i,j ,i_xyz) )
1714  call check( __line__, gsqrt(k,i,j+1,i_xyz) )
1715  call check( __line__, gsqrt(k,i,j,i_xvz) )
1716  call check( __line__, sv(k,i,j) )
1717  call check( __line__, momy0(k,i,j) )
1718 #endif
1719  duvw = dtrk * ( ( &
1720  - ( gsqrt(k,i,j+1,i_xyz) * dpres_n(k,i,j+1) &
1721  - gsqrt(k,i,j ,i_xyz) * dpres_n(k,i,j ) ) * rfdy(j) &
1722  ) / gsqrt(k,i,j,i_xvz) &
1723  + sv(k,i,j) )
1724  momy_rk(k,i,j) = momy0(k,i,j) + duvw
1725  mflx_hi(k,i,j,ydir) = gsqrt(k,i,j,i_xvz) * ( momy(k,i,j) + duvw )
1726  enddo
1727  enddo
1728  enddo
1729 #ifdef DEBUG
1730  k = iundef; i = iundef; j = iundef
1731 #endif
1732 
1733 
1734  !##### Thermodynamic Equation #####
1735 
1736  ! at (x, y, w)
1737  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
1738  do j = jjs, jje
1739  do i = iis, iie
1740  do k = ks+1, ke-2
1741 #ifdef DEBUG
1742  call check( __line__, mflx_hi(k,i,j,zdir) )
1743  call check( __line__, pott(k,i-1,j) )
1744  call check( __line__, pott(k,i ,j) )
1745  call check( __line__, pott(k,i+1,j) )
1746  call check( __line__, pott(k,i+1,j) )
1747  call check( __line__, num_diff(k,i,j,i_rhot,xdir) )
1748 #endif
1749  tflx_hi(k,i,j,zdir) = mflx_hi(k,i,j,zdir) &
1750  * ( fact_n * ( pott(k+1,i,j) + pott(k ,i,j) ) &
1751  + fact_f * ( pott(k+2,i,j) + pott(k-1,i,j) ) )
1752  enddo
1753  enddo
1754  enddo
1755 #ifdef DEBUG
1756  k = iundef; i = iundef; j = iundef
1757 #endif
1758  !$omp parallel do private(i,j) OMP_SCHEDULE_ collapse(2)
1759  do j = jjs, jje
1760  do i = iis, iie
1761  tflx_hi(ks-1,i,j,zdir) = 0.0_rp
1762  tflx_hi(ks ,i,j,zdir) = mflx_hi(ks ,i,j,zdir) * 0.5_rp * ( pott(ks+1,i,j) + pott(ks ,i,j) )
1763  tflx_hi(ke-1,i,j,zdir) = mflx_hi(ke-1,i,j,zdir) * 0.5_rp * ( pott(ke ,i,j) + pott(ke-1,i,j) )
1764  tflx_hi(ke ,i,j,zdir) = 0.0_rp
1765  enddo
1766  enddo
1767 #ifdef DEBUG
1768  k = iundef; i = iundef; j = iundef
1769 #endif
1770  ! at (u, y, z)
1771  if ( .not. twod ) then
1772  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
1773  do j = jjs, jje
1774  do i = iis-1, iie
1775  do k = ks, ke
1776 #ifdef DEBUG
1777  call check( __line__, mflx_hi(k,i,j,xdir) )
1778  call check( __line__, pott(k,i-1,j) )
1779  call check( __line__, pott(k,i ,j) )
1780  call check( __line__, pott(k,i+1,j) )
1781  call check( __line__, pott(k,i+1,j) )
1782  call check( __line__, num_diff(k,i,j,i_rhot,xdir) )
1783 #endif
1784  tflx_hi(k,i,j,xdir) = mflx_hi(k,i,j,xdir) &
1785  * ( fact_n * ( pott(k,i+1,j)+pott(k,i ,j) ) &
1786  + fact_f * ( pott(k,i+2,j)+pott(k,i-1,j) ) )
1787  enddo
1788  enddo
1789  enddo
1790  end if
1791 #ifdef DEBUG
1792  k = iundef; i = iundef; j = iundef
1793 #endif
1794  ! at (x, v, z)
1795  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
1796  do j = jjs-1, jje
1797  do i = iis, iie
1798  do k = ks, ke
1799 #ifdef DEBUG
1800  call check( __line__, mflx_hi(k,i,j,ydir) )
1801  call check( __line__, pott(k,i,j-1) )
1802  call check( __line__, pott(k,i,j ) )
1803  call check( __line__, pott(k,i,j+1) )
1804  call check( __line__, pott(k,i,j+2) )
1805  call check( __line__, num_diff(k,i,j,i_rhot,ydir) )
1806 #endif
1807  tflx_hi(k,i,j,ydir) = mflx_hi(k,i,j,ydir) &
1808  * ( fact_n * ( pott(k,i,j+1)+pott(k,i,j ) ) &
1809  + fact_f * ( pott(k,i,j+2)+pott(k,i,j-1) ) )
1810  enddo
1811  enddo
1812  enddo
1813 #ifdef DEBUG
1814  k = iundef; i = iundef; j = iundef
1815 #endif
1816 
1817  ! rho*theta
1818  if ( twod ) then
1819  !$omp parallel do
1820  do j = jjs, jje
1821  do k = ks, ke
1822  rhot_rk(k,is,j) = rhot0(k,is,j) &
1823  + dtrk * ( - ( ( tflx_hi(k,is,j,zdir) - tflx_hi(k-1,is,j ,zdir) ) * rcdz(k) &
1824  + ( tflx_hi(k,is,j,ydir) - tflx_hi(k ,is,j-1,ydir) ) * rcdy(j) ) &
1825  / gsqrt(k,is,j,i_xyz) &
1826  + st(k,is,j) )
1827  enddo
1828  enddo
1829  else
1830  !$omp parallel do
1831  do j = jjs, jje
1832  do i = iis, iie
1833  do k = ks, ke
1834  rhot_rk(k,i,j) = rhot0(k,i,j) &
1835  + dtrk * ( - ( ( tflx_hi(k,i,j,zdir) - tflx_hi(k-1,i ,j ,zdir) ) * rcdz(k) &
1836  + ( tflx_hi(k,i,j,xdir) - tflx_hi(k ,i-1,j ,xdir) ) * rcdx(i) &
1837  + ( tflx_hi(k,i,j,ydir) - tflx_hi(k ,i ,j-1,ydir) ) * rcdy(j) ) &
1838  / gsqrt(k,i,j,i_xyz) &
1839  + st(k,i,j) )
1840  enddo
1841  enddo
1842  enddo
1843  end if
1844 
1845  !##### continuous equation #####
1846 
1847  ! total momentum flux
1848  !$omp parallel do
1849  do j = jjs, jje
1850  do i = iis, iie
1851  do k = ks, ke-1
1852  mflx_hi(k,i,j,zdir) = mflx_hi(k,i,j,zdir) + mflx_hi2(k,i,j,zdir)
1853  enddo
1854  enddo
1855  enddo
1856  if ( .not. twod ) then
1857  !$omp parallel do
1858  do j = jjs , jje
1859  do i = iis-1, iie
1860  do k = ks, ke
1861  mflx_hi(k,i,j,xdir) = mflx_hi(k,i,j,xdir) + mflx_hi2(k,i,j,xdir)
1862  enddo
1863  enddo
1864  enddo
1865  end if
1866  !$omp parallel do
1867  do j = jjs-1, jje
1868  do i = iis , iie
1869  do k = ks, ke
1870  mflx_hi(k,i,j,ydir) = mflx_hi(k,i,j,ydir) + mflx_hi2(k,i,j,ydir)
1871  enddo
1872  enddo
1873  enddo
1874 
1875  !--- update density
1876  if ( twod ) then
1877  !$omp parallel do private(j,k) OMP_SCHEDULE_
1878  do j = jjs, jje
1879  do k = ks, ke
1880 #ifdef DEBUG
1881  call check( __line__, dens0(k,is,j) )
1882  call check( __line__, mflx_hi(k ,is,j ,zdir) )
1883  call check( __line__, mflx_hi(k-1,is,j ,zdir) )
1884  call check( __line__, mflx_hi(k ,is,j ,ydir) )
1885  call check( __line__, mflx_hi(k ,is,j-1,ydir) )
1886  call check( __line__, dens_t(k,is,j) )
1887 #endif
1888  dens_rk(k,is,j) = dens0(k,is,j) &
1889  + dtrk * ( - ( ( mflx_hi(k,is,j,zdir)-mflx_hi(k-1,is,j, zdir) ) * rcdz(k) &
1890  + ( mflx_hi(k,is,j,ydir)-mflx_hi(k ,is,j-1,ydir) ) * rcdy(j) ) &
1891  / gsqrt(k,is,j,i_xyz) & ! divergence
1892  + dens_t(k,is,j) )
1893  enddo
1894  enddo
1895  else
1896  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
1897  do j = jjs, jje
1898  do i = iis, iie
1899  do k = ks, ke
1900 #ifdef DEBUG
1901  call check( __line__, dens0(k,i,j) )
1902  call check( __line__, mflx_hi(k ,i ,j ,zdir) )
1903  call check( __line__, mflx_hi(k-1,i ,j ,zdir) )
1904  call check( __line__, mflx_hi(k ,i ,j ,xdir) )
1905  call check( __line__, mflx_hi(k ,i-1,j ,xdir) )
1906  call check( __line__, mflx_hi(k ,i ,j ,ydir) )
1907  call check( __line__, mflx_hi(k ,i ,j-1,ydir) )
1908  call check( __line__, dens_t(k,i,j) )
1909 #endif
1910  dens_rk(k,i,j) = dens0(k,i,j) &
1911  + dtrk * ( - ( ( mflx_hi(k,i,j,zdir)-mflx_hi(k-1,i ,j, zdir) ) * rcdz(k) &
1912  + ( mflx_hi(k,i,j,xdir)-mflx_hi(k ,i-1,j, xdir) ) * rcdx(i) &
1913  + ( mflx_hi(k,i,j,ydir)-mflx_hi(k ,i, j-1,ydir) ) * rcdy(j) ) &
1914  / gsqrt(k,i,j,i_xyz) & ! divergence
1915  + dens_t(k,i,j) )
1916  enddo
1917  enddo
1918  enddo
1919  end if
1920 #ifdef DEBUG
1921  k = iundef; i = iundef; j = iundef
1922 #endif
1923 
1924 
1925 #ifdef DEBUG
1926  call check_pres( &
1927  dpres_n, dpres, rhot_rk, rhot, dens_rk, dens, b,&
1928  rt2p )
1929 #endif
1930 
1931  enddo
1932  enddo
1933 
1934  return
scale_atmos_dyn_fvm_fct
module Atmosphere / Dynamics common
Definition: scale_atmos_dyn_fvm_fct.F90:12
scale_atmos_dyn_fvm_fct::atmos_dyn_fvm_fct
subroutine, public atmos_dyn_fvm_fct(qflx_anti, phi_in, DENS0, DENS, qflx_hi, qflx_lo, mflx_hi, rdz, rdx, rdy, GSQRT, MAPF, TwoD, dt, flag_vect)
Setup.
Definition: scale_atmos_dyn_fvm_fct.F90:72
scale_atmos_dyn_fvm_flux
module scale_atmos_dyn_fvm_flux
Definition: scale_atmos_dyn_fvm_flux.F90:13
scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxj13_uyz
procedure(flux_j), pointer, public atmos_dyn_fvm_fluxj13_uyz
Definition: scale_atmos_dyn_fvm_flux.F90:225
scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxx_uyz
procedure(flux_mom), pointer, public atmos_dyn_fvm_fluxx_uyz
Definition: scale_atmos_dyn_fvm_flux.F90:214
scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxx_xyw
procedure(flux_mom), pointer, public atmos_dyn_fvm_fluxx_xyw
Definition: scale_atmos_dyn_fvm_flux.F90:187
scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxy_xvz
procedure(flux_mom), pointer, public atmos_dyn_fvm_fluxy_xvz
Definition: scale_atmos_dyn_fvm_flux.F90:246
scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxx_xvz
procedure(flux_mom), pointer, public atmos_dyn_fvm_fluxx_xvz
Definition: scale_atmos_dyn_fvm_flux.F90:241
scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxj23_xvz
procedure(flux_j), pointer, public atmos_dyn_fvm_fluxj23_xvz
Definition: scale_atmos_dyn_fvm_flux.F90:257
scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxz_xyw
procedure(flux_wz), pointer, public atmos_dyn_fvm_fluxz_xyw
Definition: scale_atmos_dyn_fvm_flux.F90:182
scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxz_uyz
procedure(flux_z), pointer, public atmos_dyn_fvm_fluxz_uyz
Definition: scale_atmos_dyn_fvm_flux.F90:209
scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxj13_xyw
procedure(flux_j), pointer, public atmos_dyn_fvm_fluxj13_xyw
Definition: scale_atmos_dyn_fvm_flux.F90:198
scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxj23_uyz
procedure(flux_j), pointer, public atmos_dyn_fvm_fluxj23_uyz
Definition: scale_atmos_dyn_fvm_flux.F90:230
scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxx_xyz
procedure(flux_phi), pointer, public atmos_dyn_fvm_fluxx_xyz
Definition: scale_atmos_dyn_fvm_flux.F90:170
scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxy_xyw
procedure(flux_mom), pointer, public atmos_dyn_fvm_fluxy_xyw
Definition: scale_atmos_dyn_fvm_flux.F90:192
scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxz_xyz
procedure(flux_phi), pointer, public atmos_dyn_fvm_fluxz_xyz
Definition: scale_atmos_dyn_fvm_flux.F90:165
scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxj13_xvz
procedure(flux_j), pointer, public atmos_dyn_fvm_fluxj13_xvz
Definition: scale_atmos_dyn_fvm_flux.F90:252
scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxz_xvz
procedure(flux_z), pointer, public atmos_dyn_fvm_fluxz_xvz
Definition: scale_atmos_dyn_fvm_flux.F90:236
scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxy_xyz
procedure(flux_phi), pointer, public atmos_dyn_fvm_fluxy_xyz
Definition: scale_atmos_dyn_fvm_flux.F90:175
scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_flux_valuew_z
procedure(valuew), pointer, public atmos_dyn_fvm_flux_valuew_z
Definition: scale_atmos_dyn_fvm_flux.F90:160
scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxy_uyz
procedure(flux_mom), pointer, public atmos_dyn_fvm_fluxy_uyz
Definition: scale_atmos_dyn_fvm_flux.F90:219
scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxj23_xyw
procedure(flux_j), pointer, public atmos_dyn_fvm_fluxj23_xyw
Definition: scale_atmos_dyn_fvm_flux.F90:203
scale_atmos_grid_cartesc_index
module atmosphere / grid / cartesC index
Definition: scale_atmos_grid_cartesC_index.F90:12
scale_atmos_grid_cartesc_index::i_xvw
integer, public i_xvw
Definition: scale_atmos_grid_cartesC_index.F90:94
scale_atmos_grid_cartesc_index::xdir
integer, parameter, public xdir
Definition: scale_atmos_grid_cartesC_index.F90:33
scale_atmos_grid_cartesc_index::i_uyz
integer, public i_uyz
Definition: scale_atmos_grid_cartesC_index.F90:95
scale_atmos_grid_cartesc_index::i_uy
integer, public i_uy
Definition: scale_atmos_grid_cartesC_index.F90:100
scale_atmos_grid_cartesc_index::i_xy
integer, public i_xy
Definition: scale_atmos_grid_cartesC_index.F90:99
scale_atmos_grid_cartesc_index::ie
integer, public ie
end point of inner domain: x, local
Definition: scale_atmos_grid_cartesC_index.F90:54
scale_atmos_grid_cartesc_index::i_uvz
integer, public i_uvz
Definition: scale_atmos_grid_cartesC_index.F90:97
scale_atmos_grid_cartesc_index::je
integer, public je
end point of inner domain: y, local
Definition: scale_atmos_grid_cartesC_index.F90:56
scale_atmos_grid_cartesc_index::is
integer, public is
start point of inner domain: x, local
Definition: scale_atmos_grid_cartesC_index.F90:53
scale_atmos_grid_cartesc_index::ia
integer, public ia
Definition: scale_atmos_grid_cartesC_index.F90:48
scale_atmos_grid_cartesc_index::i_xv
integer, public i_xv
Definition: scale_atmos_grid_cartesC_index.F90:101
scale_atmos_grid_cartesc_index::i_xvz
integer, public i_xvz
Definition: scale_atmos_grid_cartesC_index.F90:96
scale_atmos_grid_cartesc_index::i_xyw
integer, public i_xyw
Definition: scale_atmos_grid_cartesC_index.F90:92
scale_atmos_grid_cartesc_index::iblock
integer, public iblock
block size for cache blocking: x
Definition: scale_atmos_grid_cartesC_index.F90:40
scale_atmos_grid_cartesc_index::i_xyz
integer, public i_xyz
Definition: scale_atmos_grid_cartesC_index.F90:91
scale_atmos_grid_cartesc_index::zdir
integer, parameter, public zdir
Definition: scale_atmos_grid_cartesC_index.F90:32
scale_atmos_grid_cartesc_index::jblock
integer, public jblock
block size for cache blocking: y
Definition: scale_atmos_grid_cartesC_index.F90:41
scale_atmos_grid_cartesc_index::ka
integer, public ka
Definition: scale_atmos_grid_cartesC_index.F90:47
scale_atmos_grid_cartesc_index::js
integer, public js
start point of inner domain: y, local
Definition: scale_atmos_grid_cartesC_index.F90:55
scale_atmos_grid_cartesc_index::ydir
integer, parameter, public ydir
Definition: scale_atmos_grid_cartesC_index.F90:34
scale_atmos_grid_cartesc_index::ke
integer, public ke
end point of inner domain: z, local
Definition: scale_atmos_grid_cartesC_index.F90:52
scale_atmos_grid_cartesc_index::ks
integer, public ks
start point of inner domain: z, local
Definition: scale_atmos_grid_cartesC_index.F90:51
scale_atmos_grid_cartesc_index::i_uyw
integer, public i_uyw
Definition: scale_atmos_grid_cartesC_index.F90:93
scale_atmos_grid_cartesc_index::i_uv
integer, public i_uv
Definition: scale_atmos_grid_cartesC_index.F90:102
scale_comm_cartesc
module COMMUNICATION
Definition: scale_comm_cartesC.F90:11
scale_const
module CONSTANT
Definition: scale_const.F90:11
scale_const::const_pre00
real(rp), public const_pre00
pressure reference [Pa]
Definition: scale_const.F90:97
scale_const::const_grav
real(rp), public const_grav
standard acceleration of gravity [m/s2]
Definition: scale_const.F90:49

References scale_atmos_dyn_fvm_fct::atmos_dyn_fvm_fct(), scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_flux_valuew_z, scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxj13_uyz, scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxj13_xvz, scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxj13_xyw, scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxj23_uyz, scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxj23_xvz, scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxj23_xyw, scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxx_uyz, scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxx_xvz, scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxx_xyw, scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxx_xyz, scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxy_uyz, scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxy_xvz, scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxy_xyw, scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxy_xyz, scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxz_uyz, scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxz_xvz, scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxz_xyw, scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxz_xyz, scale_debug::check(), 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_atmos_grid_cartesc_index::i_uv, scale_atmos_grid_cartesc_index::i_uvz, scale_atmos_grid_cartesc_index::i_uy, scale_atmos_grid_cartesc_index::i_uyw, scale_atmos_grid_cartesc_index::i_uyz, scale_atmos_grid_cartesc_index::i_xv, scale_atmos_grid_cartesc_index::i_xvw, scale_atmos_grid_cartesc_index::i_xvz, scale_atmos_grid_cartesc_index::i_xy, scale_atmos_grid_cartesc_index::i_xyw, scale_atmos_grid_cartesc_index::i_xyz, scale_atmos_grid_cartesc_index::ia, scale_atmos_grid_cartesc_index::iblock, scale_atmos_grid_cartesc_index::ie, scale_atmos_grid_cartesc_index::is, scale_atmos_grid_cartesc_index::jblock, scale_atmos_grid_cartesc_index::je, scale_atmos_grid_cartesc_index::js, scale_tracer::k, scale_atmos_grid_cartesc_index::ka, scale_atmos_grid_cartesc_index::ke, scale_atmos_grid_cartesc_index::ks, make_matrix(), solve_bicgstab(), scale_atmos_grid_cartesc_index::xdir, scale_atmos_grid_cartesc_index::ydir, and scale_atmos_grid_cartesc_index::zdir.

Referenced by scale_atmos_dyn_tstep_short::atmos_dyn_tstep_short_regist().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ solve_bicgstab()

subroutine scale_atmos_dyn_tstep_short_fvm_hivi::solve_bicgstab ( real(rp), dimension(ka,ia,ja), intent(out)  DPRES_N,
real(rp), dimension(ka,ia,ja), intent(in)  DPRES,
real(rp), dimension(7,ka,ia,ja), intent(in)  M,
real(rp), dimension(ka,ia,ja), intent(in)  B,
logical, intent(in)  TwoD 
)

Definition at line 1938 of file scale_atmos_dyn_tstep_short_fvm_hivi.F90.

1943  use scale_prc, only: &
1944  prc_abort
1945  use scale_comm_cartesc, only: &
1946  comm_world, &
1947  comm_vars8, &
1948  comm_wait
1949  implicit none
1950  real(RP), intent(out) :: DPRES_N(KA,IA,JA)
1951  real(RP), intent(in) :: DPRES(KA,IA,JA)
1952  real(RP), intent(in) :: M(7,KA,IA,JA)
1953  real(RP), intent(in) :: B(KA,IA,JA)
1954  logical, intent(in) :: TwoD
1955 
1956  real(RP) :: r0(KA,IA,JA)
1957 
1958  real(RP) :: p(KA,IA,JA)
1959  real(RP) :: Mp(KA,IA,JA)
1960  real(RP) :: s(KA,IA,JA)
1961  real(RP) :: Ms(KA,IA,JA)
1962  real(RP) :: al, be, w
1963 
1964  real(RP), pointer :: r(:,:,:)
1965  real(RP), pointer :: rn(:,:,:)
1966  real(RP), pointer :: swap(:,:,:)
1967  real(RP), target :: v0(KA,IA,JA)
1968  real(RP), target :: v1(KA,IA,JA)
1969  real(RP) :: r0r
1970  real(RP) :: norm, error, error2
1971 
1972  real(RP) :: iprod(2)
1973  real(RP) :: buf(2)
1974 
1975  integer :: k, i, j
1976  integer :: iis, iie, jjs, jje
1977  integer :: iter
1978  integer :: ierror
1979 
1980 #ifdef DEBUG
1981  r0(:,:,:) = undef
1982  p(:,:,:) = undef
1983  mp(:,:,:) = undef
1984  s(:,:,:) = undef
1985  ms(:,:,:) = undef
1986 
1987  v0(:,:,:) = undef
1988  v1(:,:,:) = undef
1989 #endif
1990 
1991  r => v0
1992  rn => v1
1993 
1994  call mul_matrix( v1, m, dpres, twod ) ! v1 = M x0
1995 
1996  norm = 0.0_rp
1997  r0r = 0.0_rp
1998 
1999  do jjs = js, je, jblock
2000  jje = jjs+jblock-1
2001  do iis = is, ie, iblock
2002  iie = iis+iblock-1
2003 
2004  do j = jjs, jje
2005  do i = iis, iie
2006  do k = ks, ke
2007  norm = norm + b(k,i,j)**2
2008  enddo
2009  enddo
2010  enddo
2011 #ifdef DEBUG
2012  k = iundef; i = iundef; j = iundef
2013 #endif
2014 
2015  ! r = b - M x0
2016  do j = jjs, jje
2017  do i = iis, iie
2018  do k = ks, ke
2019 #ifdef DEBUG
2020  call check( __line__, b(k,i,j) )
2021  call check( __line__, v1(k,i,j) )
2022 #endif
2023  r(k,i,j) = b(k,i,j) - v1(k,i,j)
2024  enddo
2025  enddo
2026  enddo
2027 #ifdef DEBUG
2028  k = iundef; i = iundef; j = iundef
2029 #endif
2030 
2031  do j = jjs, jje
2032  do i = iis, iie
2033  do k = ks, ke
2034  r0(k,i,j) = r(k,i,j)
2035  p(k,i,j) = r(k,i,j)
2036  enddo
2037  enddo
2038  enddo
2039 #ifdef DEBUG
2040  k = iundef; i = iundef; j = iundef
2041 #endif
2042 
2043  do j = jjs, jje
2044  do i = iis, iie
2045  do k = ks, ke
2046 #ifdef DEBUG
2047  call check( __line__, r(k,i,j) )
2048  call check( __line__, r0(k,i,j) )
2049 #endif
2050  r0r = r0r + r0(k,i,j) * r(k,i,j)
2051  enddo
2052  enddo
2053  enddo
2054 #ifdef DEBUG
2055  k = iundef; i = iundef; j = iundef
2056 #endif
2057 
2058  enddo
2059  enddo
2060 
2061  do j = js-1, je+1
2062  do i = is-1, ie+1
2063  do k = ks, ke
2064  dpres_n(k,i,j) = dpres(k,i,j)
2065  end do
2066  end do
2067  end do
2068 
2069 
2070  iprod(1) = r0r
2071  iprod(2) = norm
2072  call mpi_allreduce(iprod, buf, 2, mtype, mpi_sum, comm_world, ierror)
2073  r0r = buf(1)
2074  norm = buf(2)
2075 
2076  error2 = norm
2077 
2078  do iter = 1, itmax
2079 
2080  error = 0.0_rp
2081  do j = js, je
2082  do i = is, ie
2083  do k = ks, ke
2084 #ifdef DEBUG
2085  call check( __line__, r(k,i,j) )
2086 #endif
2087  error = error + r(k,i,j)**2
2088  enddo
2089  enddo
2090  enddo
2091 #ifdef DEBUG
2092  k = iundef; i = iundef; j = iundef
2093 #endif
2094  call mpi_allreduce(error, buf, 1, mtype, mpi_sum, comm_world, ierror)
2095  error = buf(1)
2096 
2097 #ifdef DEBUG
2098  log_info("solve_bicgstab",*) iter, error/norm
2099 #endif
2100  if ( sqrt(error/norm) < epsilon .OR. error > error2 ) then
2101 #ifdef DEBUG
2102  log_info("solve_bicgstab",*) "Bi-CGSTAB converged:", iter
2103 #endif
2104  exit
2105  endif
2106  error2 = error
2107 
2108  call comm_vars8( p, 1 )
2109  call comm_wait ( p, 1 )
2110  call mul_matrix( mp, m, p, twod )
2111 
2112  iprod(1) = 0.0_rp
2113  do j = js, je
2114  do i = is, ie
2115  do k = ks, ke
2116 #ifdef DEBUG
2117  call check( __line__, r0(k,i,j) )
2118  call check( __line__, mp(k,i,j) )
2119 #endif
2120  iprod(1) = iprod(1) + r0(k,i,j) * mp(k,i,j)
2121  enddo
2122  enddo
2123  enddo
2124 #ifdef DEBUG
2125  k = iundef; i = iundef; j = iundef
2126 #endif
2127  call mpi_allreduce(iprod, buf, 1, mtype, mpi_sum, comm_world, ierror)
2128  al = r0r / buf(1) ! (r0,r) / (r0,Mp)
2129 
2130  do j = js, je
2131  do i = is, ie
2132  do k = ks, ke
2133 #ifdef DEBUG
2134  call check( __line__, r(k,i,j) )
2135  call check( __line__, mp(k,i,j) )
2136 #endif
2137  s(k,i,j) = r(k,i,j) - al*mp(k,i,j)
2138  enddo
2139  enddo
2140  enddo
2141 #ifdef DEBUG
2142  k = iundef; i = iundef; j = iundef
2143 #endif
2144 
2145  call comm_vars8( s, 1 )
2146  call comm_wait ( s, 1 )
2147  call mul_matrix( ms, m, s, twod )
2148  iprod(1) = 0.0_rp
2149  iprod(2) = 0.0_rp
2150  do j = js, je
2151  do i = is, ie
2152  do k = ks, ke
2153 #ifdef DEBUG
2154  call check( __line__, ms(k,i,j) )
2155  call check( __line__, s(k,i,j) )
2156 #endif
2157  iprod(1) = iprod(1) + ms(k,i,j) * s(k,i,j)
2158  iprod(2) = iprod(2) + ms(k,i,j) * ms(k,i,j)
2159  enddo
2160  enddo
2161  enddo
2162 #ifdef DEBUG
2163  k = iundef; i = iundef; j = iundef
2164 #endif
2165  call mpi_allreduce(iprod, buf, 2, mtype, mpi_sum, comm_world, ierror)
2166  w = buf(1) / buf(2) ! (Ms,s) / (Ms,Ms)
2167 
2168  iprod(1) = 0.0_rp
2169 
2170  do jjs = js, je, jblock
2171  jje = jjs+jblock-1
2172  do iis = is, ie, iblock
2173  iie = iis+iblock-1
2174 
2175  do j = jjs, jje
2176  do i = iis, iie
2177  do k = ks, ke
2178 #ifdef DEBUG
2179  call check( __line__, dpres_n(k,i,j) )
2180  call check( __line__, p(k,i,j) )
2181  call check( __line__, s(k,i,j) )
2182 #endif
2183  dpres_n(k,i,j) = dpres_n(k,i,j) + al*p(k,i,j) + w*s(k,i,j)
2184  enddo
2185  enddo
2186  enddo
2187 #ifdef DEBUG
2188  k = iundef; i = iundef; j = iundef
2189 #endif
2190 
2191  do j = jjs, jje
2192  do i = iis, iie
2193  do k = ks, ke
2194 #ifdef DEBUG
2195  call check( __line__, s(k,i,j) )
2196  call check( __line__, ms(k,i,j) )
2197 #endif
2198  rn(k,i,j) = s(k,i,j) - w*ms(k,i,j)
2199  enddo
2200  enddo
2201  enddo
2202 #ifdef DEBUG
2203  k = iundef; i = iundef; j = iundef
2204 #endif
2205 
2206  do j = jjs, jje
2207  do i = iis, iie
2208  do k = ks, ke
2209 #ifdef DEBUG
2210  call check( __line__, r0(k,i,j) )
2211  call check( __line__, rn(k,i,j) )
2212 #endif
2213  iprod(1) = iprod(1) + r0(k,i,j) * rn(k,i,j)
2214  enddo
2215  enddo
2216  enddo
2217 #ifdef DEBUG
2218  k = iundef; i = iundef; j = iundef
2219 #endif
2220 
2221  enddo
2222  enddo
2223 
2224  be = al/w / r0r
2225 
2226  call mpi_allreduce(iprod, r0r, 1, mtype, mpi_sum, comm_world, ierror)
2227 
2228  be = be * r0r ! al/w * (r0,rn)/(r0,r)
2229  do j = js, je
2230  do i = is, ie
2231  do k = ks, ke
2232 #ifdef DEBUG
2233  call check( __line__, rn(k,i,j) )
2234  call check( __line__, p(k,i,j) )
2235  call check( __line__, mp(k,i,j) )
2236 #endif
2237  p(k,i,j) = rn(k,i,j) + be * ( p(k,i,j) - w*mp(k,i,j) )
2238  enddo
2239  enddo
2240  enddo
2241 #ifdef DEBUG
2242  k = iundef; i = iundef; j = iundef
2243 #endif
2244 
2245  swap => rn
2246  rn => r
2247  r => swap
2248 #ifdef DEBUG
2249  rn(:,:,:) = undef
2250 #endif
2251  enddo
2252 
2253  if ( iter >= itmax ) then
2254  log_error("solve_bicgstab",*) 'not converged', error, norm
2255  call prc_abort
2256  endif
2257 
2258  return
scale_comm_cartesc::comm_world
integer, public comm_world
communication world ID
Definition: scale_comm_cartesC.F90:106

References scale_debug::check(), scale_comm_cartesc::comm_world, scale_atmos_grid_cartesc_index::iblock, scale_atmos_grid_cartesc_index::ie, scale_atmos_grid_cartesc_index::is, scale_atmos_grid_cartesc_index::jblock, scale_atmos_grid_cartesc_index::je, scale_atmos_grid_cartesc_index::js, scale_atmos_grid_cartesc_index::ke, scale_atmos_grid_cartesc_index::ks, and scale_prc::prc_abort().

Referenced by atmos_dyn_tstep_short_fvm_hivi().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ make_matrix()

subroutine scale_atmos_dyn_tstep_short_fvm_hivi::make_matrix ( real(rp), dimension(7,ka,ia,ja), intent(inout)  M,
real(rp), dimension(ka,ia,ja), intent(in)  POTT,
real(rp), dimension(ka,ia,ja), intent(in)  RCs2T,
real(rp), intent(in)  GRAV,
real(rp), dimension(ka,ia,ja,7), intent(in)  G,
real(rp), intent(in)  J33G,
real(rp), dimension(ka), intent(in)  RCDZ,
real(rp), dimension(ka-1), intent(in)  RFDZ,
real(rp), dimension(ia), intent(in)  RCDX,
real(rp), dimension(ia-1), intent(in)  RFDX,
real(rp), dimension(ja), intent(in)  RCDY,
real(rp), dimension(ja-1), intent(in)  RFDY,
real(rp), dimension(ka-1), intent(in)  FDZ,
real(rp), intent(in)  rdt,
real(rp), intent(in)  FACT_N,
real(rp), intent(in)  FACT_F,
logical, intent(in)  TwoD,
integer, intent(in)  I_XYZ,
integer, intent(in)  I_XYW,
integer, intent(in)  IIS,
integer, intent(in)  IIE,
integer, intent(in)  JJS,
integer, intent(in)  JJE 
)

Definition at line 2261 of file scale_atmos_dyn_tstep_short_fvm_hivi.F90.

2271  implicit none
2272  real(RP), intent(inout) :: M(7,KA,IA,JA)
2273  real(RP), intent(in) :: POTT(KA,IA,JA)
2274  real(RP), intent(in) :: RCs2T(KA,IA,JA)
2275  real(RP), intent(in) :: GRAV
2276  real(RP), intent(in) :: G(KA,IA,JA,7)
2277  real(RP), intent(in) :: J33G
2278  real(RP), intent(in) :: RCDZ(KA)
2279  real(RP), intent(in) :: RFDZ(KA-1)
2280  real(RP), intent(in) :: RFDX(IA-1)
2281  real(RP), intent(in) :: RCDX(IA)
2282  real(RP), intent(in) :: RCDY(JA)
2283  real(RP), intent(in) :: RFDY(JA-1)
2284  real(RP), intent(in) :: FDZ(KA-1)
2285  real(RP), intent(in) :: rdt
2286  real(RP), intent(in) :: FACT_N
2287  real(RP), intent(in) :: FACT_F
2288  logical, intent(in) :: TwoD
2289  integer, intent(in) :: I_XYZ
2290  integer, intent(in) :: I_XYW
2291  integer, intent(in) :: IIS
2292  integer, intent(in) :: IIE
2293  integer, intent(in) :: JJS
2294  integer, intent(in) :: JJE
2295 
2296  integer :: k, i, j
2297 
2298  if ( twod ) then
2299  i = is
2300  !$omp parallel do
2301  do j = jjs, jje
2302  do k = ks+2, ke-2
2303 #ifdef DEBUG
2304  call check( __line__, pott(k-2,i,j) )
2305  call check( __line__, pott(k-1,i,j) )
2306  call check( __line__, pott(k ,i,j) )
2307  call check( __line__, pott(k+1,i,j) )
2308  call check( __line__, pott(k+2,i,j) )
2309  call check( __line__, pott(k,i ,j) )
2310  call check( __line__, pott(k,i,j-2) )
2311  call check( __line__, pott(k,i,j-1) )
2312  call check( __line__, pott(k,i,j ) )
2313  call check( __line__, pott(k,i,j+1) )
2314  call check( __line__, pott(k,i,j+2) )
2315  call check( __line__, g(k-1,i,j,i_xyw) )
2316  call check( __line__, g(k+1,i,j,i_xyw) )
2317  call check( __line__, g(k,i,j,i_xyz) )
2318  call check( __line__, rcs2t(k-1,i,j) )
2319  call check( __line__, rcs2t(k ,i,j) )
2320  call check( __line__, rcs2t(k+1,i,j) )
2321 #endif
2322  ! k,i,j
2323  m(1,k,i,j) = &
2324  - ( ( fact_n * (pott(k+1,i,j)+pott(k ,i,j)) &
2325  + fact_f * (pott(k+2,i,j)+pott(k-1,i,j)) ) * rfdz(k ) / g(k ,i,j,i_xyw) &
2326  + ( fact_n * (pott(k ,i,j)+pott(k-1,i,j)) &
2327  + fact_f * (pott(k+1,i,j)+pott(k-2,i,j)) ) * rfdz(k-1) / g(k-1,i,j,i_xyw) &
2328  ) * j33g * j33g * rfdz(k) &
2329  - ( ( fact_n * (pott(k,i,j+1)+pott(k,i,j )) &
2330  + fact_f * (pott(k,i,j+2)+pott(k,i,j-1)) ) * rfdy(j ) &
2331  + ( fact_n * (pott(k,i,j )+pott(k,i,j-1)) &
2332  + fact_f * (pott(k,i,j-1)+pott(k,i,j-2)) ) * rfdy(j-1) &
2333  ) * g(k,i,j,i_xyz) * rfdy(j) &
2334  - g(k,i,j,i_xyz) * rcs2t(k,i,j) * rdt * rdt
2335  ! k-1
2336  m(2,k,i,j) = j33g * j33g / g(k-1,i,j,i_xyw) &
2337  * ( fact_n * (pott(k ,i,j)+pott(k-1,i,j)) &
2338  + fact_f * (pott(k+1,i,j)+pott(k-2,i,j)) ) &
2339  * rfdz(k-1) * rcdz(k) &
2340  - grav * j33g * rcs2t(k-1,i,j) / ( fdz(k)+fdz(k-1) )
2341  ! k+1
2342  m(3,k,i,j) = j33g * j33g / g(k+1,i,j,i_xyw) &
2343  * ( fact_n * (pott(k+1,i,j)+pott(k ,i,j)) &
2344  + fact_f * (pott(k+2,i,j)+pott(k-1,i,j)) ) &
2345  * rfdz(k ) * rcdz(k) &
2346  + grav * j33g * rcs2t(k+1,i,j) / ( fdz(k)+fdz(k-1) )
2347  enddo
2348  enddo
2349 
2350  !$omp parallel do
2351  do j = jjs, jje
2352 #ifdef DEBUG
2353  call check( __line__, pott(ks ,i,j) )
2354  call check( __line__, pott(ks+1,i,j) )
2355  call check( __line__, pott(ks,i ,j) )
2356  call check( __line__, pott(ks,i,j-2) )
2357  call check( __line__, pott(ks,i,j-1) )
2358  call check( __line__, pott(ks,i,j ) )
2359  call check( __line__, pott(ks,i,j+1) )
2360  call check( __line__, pott(ks,i,j+2) )
2361  call check( __line__, g(ks,i,j,i_xyz) )
2362  call check( __line__, rcs2t(ks,i,j) )
2363  call check( __line__, pott(ks ,i,j) )
2364  call check( __line__, pott(ks+1,i,j) )
2365  call check( __line__, g(ks+1,i,j,i_xyw) )
2366  call check( __line__, rcs2t(ks+1,i,j) )
2367 #endif
2368  ! k,i,j
2369  m(1,ks,i,j) = &
2370  - ( 0.5_rp * (pott(ks+1,i,j)+pott(ks ,i,j)) * rfdz(ks ) &
2371  ) * j33g * j33g / g(ks ,i,j,i_xyw) * rfdz(ks) &
2372  - ( ( fact_n * (pott(ks,i,j+1)+pott(ks,i,j )) &
2373  + fact_f * (pott(ks,i,j+2)+pott(ks,i,j-1)) ) * rfdy(j ) &
2374  + ( fact_n * (pott(ks,i,j )+pott(ks,i,j-1)) &
2375  + fact_f * (pott(ks,i,j-1)+pott(ks,i,j-2)) ) * rfdy(j-1) &
2376  ) * g(ks,i,j,i_xyz) * rfdy(j) &
2377  - g(ks,i,j,i_xyz) * rcs2t(ks,i,j) * rdt * rdt &
2378  + grav * j33g * 0.5_rp * rcs2t(ks,i,j) * rcdz(ks)
2379  ! k+1
2380  m(3,ks,i,j) = j33g * j33g / g(ks+1,i,j,i_xyw) &
2381  * 0.5_rp * (pott(ks+1,i,j)+pott(ks ,i,j)) &
2382  * rfdz(ks ) * rcdz(ks) &
2383  + grav * j33g * 0.5_rp * rcs2t(ks+1,i,j) * rcdz(ks)
2384 #ifdef DEBUG
2385  call check( __line__, pott(ks ,i,j) )
2386  call check( __line__, pott(ks+1,i,j) )
2387  call check( __line__, pott(ks+2,i,j) )
2388  call check( __line__, pott(ks+3,i,j) )
2389  call check( __line__, pott(ks+1,i ,j) )
2390  call check( __line__, pott(ks+1,i,j-2) )
2391  call check( __line__, pott(ks+1,i,j-1) )
2392  call check( __line__, pott(ks+1,i,j ) )
2393  call check( __line__, pott(ks+1,i,j+1) )
2394  call check( __line__, pott(ks+1,i,j+2) )
2395  call check( __line__, g(ks ,i,j,i_xyw) )
2396  call check( __line__, g(ks+2,i,j,i_xyw) )
2397  call check( __line__, g(ks+1,i,j,i_xyz) )
2398  call check( __line__, rcs2t(ks ,i,j) )
2399  call check( __line__, rcs2t(ks+1 ,i,j) )
2400  call check( __line__, rcs2t(ks+2,i,j) )
2401 #endif
2402  ! k,i,j
2403  m(1,ks+1,i,j) = &
2404  - ( ( fact_n * (pott(ks+2,i,j)+pott(ks+1,i,j)) &
2405  + fact_f * (pott(ks+3,i,j)+pott(ks ,i,j)) ) * rfdz(ks+1) / g(ks+1,i,j,i_xyw) &
2406  + 0.5_rp * (pott(ks+1,i,j)+pott(ks ,i,j)) * rfdz(ks ) / g(ks ,i,j,i_xyw) &
2407  ) * j33g * j33g * rfdz(ks+1) &
2408  - ( ( fact_n * (pott(ks+1,i,j+1)+pott(ks+1,i,j )) &
2409  + fact_f * (pott(ks+1,i,j+2)+pott(ks+1,i,j-1)) ) * rfdy(j ) &
2410  + ( fact_n * (pott(ks+1,i,j )+pott(ks+1,i,j-1)) &
2411  + fact_f * (pott(ks+1,i,j-1)+pott(ks+1,i,j-2)) ) * rfdy(j-1) &
2412  ) * g(ks+1,i,j,i_xyz) * rfdy(j) &
2413  - g(ks+1,i,j,i_xyz) * rcs2t(ks+1,i,j) * rdt * rdt
2414  ! k-1
2415  m(2,ks+1,i,j) = j33g * j33g / g(ks,i,j,i_xyw) &
2416  * 0.5_rp * (pott(ks+1,i,j)+pott(ks ,i,j)) &
2417  * rfdz(ks ) * rcdz(ks+1) &
2418  - grav * j33g * rcs2t(ks ,i,j) / ( fdz(ks+1)+fdz(ks) )
2419  ! k+1
2420  m(3,ks+1,i,j) = j33g * j33g / g(ks+2,i,j,i_xyw) &
2421  * ( fact_n * (pott(ks+2,i,j)+pott(ks+1,i,j)) &
2422  + fact_f * (pott(ks+3,i,j)+pott(ks ,i,j)) ) &
2423  * rfdz(ks+1) * rcdz(ks+1) &
2424  + grav * j33g * rcs2t(ks+2,i,j) / ( fdz(ks+1)+fdz(ks) )
2425 #ifdef DEBUG
2426  call check( __line__, pott(ke-3,i,j) )
2427  call check( __line__, pott(ke-2,i,j) )
2428  call check( __line__, pott(ke-1,i,j) )
2429  call check( __line__, pott(ke ,i,j) )
2430  call check( __line__, pott(ke-1,i ,j) )
2431  call check( __line__, pott(ke-1,i,j-2) )
2432  call check( __line__, pott(ke-1,i,j-1) )
2433  call check( __line__, pott(ke-1,i,j ) )
2434  call check( __line__, pott(ke-1,i,j+1) )
2435  call check( __line__, pott(ke-1,i,j+2) )
2436  call check( __line__, g(ke-2,i,j,i_xyw) )
2437  call check( __line__, g(ke ,i,j,i_xyw) )
2438  call check( __line__, g(ke-1,i,j,i_xyz) )
2439  call check( __line__, rcs2t(ke-2,i,j) )
2440  call check( __line__, rcs2t(ke-1,i,j) )
2441  call check( __line__, rcs2t(ke ,i,j) )
2442 #endif
2443  ! k,i,j
2444  m(1,ke-1,i,j) = &
2445  - ( 0.5_rp * (pott(ke ,i,j)+pott(ke-1,i,j)) * rfdz(ke-1) / g(ke-1,i,j,i_xyw) &
2446  + ( fact_n * (pott(ke-1,i,j)+pott(ke-2,i,j)) &
2447  + fact_f * (pott(ke ,i,j)+pott(ke-3,i,j)) ) * rfdz(ke-2) / g(ke-2,i,j,i_xyw) &
2448  ) * j33g * j33g * rfdz(ke-1) &
2449  - ( ( fact_n * (pott(ke-1,i,j+1)+pott(ke-1,i,j )) &
2450  + fact_f * (pott(ke-1,i,j+2)+pott(ke-1,i,j-1)) ) * rfdy(j ) &
2451  + ( fact_n * (pott(ke-1,i,j )+pott(ke-1,i,j-1)) &
2452  + fact_f * (pott(ke-1,i,j-1)+pott(ke-1,i,j-2)) ) * rfdy(j-1) &
2453  ) * g(ke-1,i,j,i_xyz) * rfdy(j) &
2454  - g(ke-1,i,j,i_xyz) * rcs2t(ke-1,i,j) * rdt * rdt
2455  ! k-1
2456  m(2,ke-1,i,j) = j33g * j33g / g(ke-2,i,j,i_xyw) &
2457  * ( fact_n * (pott(ke-1,i,j)+pott(ke-2,i,j)) &
2458  + fact_f * (pott(ke ,i,j)+pott(ke-3,i,j)) ) &
2459  * rfdz(ke-2) * rcdz(ke-1) &
2460  - grav * j33g * rcs2t(ke-2,i,j) / ( fdz(ke-1)+fdz(ke-2) )
2461  ! k+1
2462  m(3,ke-1,i,j) = j33g * j33g / g(ke ,i,j,i_xyw) &
2463  * 0.5_rp * (pott(ke ,i,j)+pott(ke-1,i,j)) &
2464  * rfdz(ke-1) * rcdz(ke-1) &
2465  + grav * j33g * rcs2t(ke ,i,j) / ( fdz(ke-1)+fdz(ke-2) )
2466 #ifdef DEBUG
2467  call check( __line__, pott(ke-1,i,j) )
2468  call check( __line__, pott(ke ,i,j) )
2469  call check( __line__, pott(ke,i ,j) )
2470  call check( __line__, pott(ke,i,j-2) )
2471  call check( __line__, pott(ke,i,j-1) )
2472  call check( __line__, pott(ke,i,j ) )
2473  call check( __line__, pott(ke,i,j+1) )
2474  call check( __line__, pott(ke,i,j+2) )
2475  call check( __line__, g(ke-1,i,j,i_xyw) )
2476  call check( __line__, g(ke,i,j,i_xyz) )
2477  call check( __line__, rcs2t(ke-1,i,j) )
2478  call check( __line__, rcs2t(ke,i,j) )
2479 #endif
2480  ! k,i,j
2481  m(1,ke,i,j) = &
2482  - ( &
2483  + 0.5_rp * (pott(ke ,i,j)+pott(ke-1,i,j)) * rfdz(ke-1) / g(ke-1,i,j,i_xyw) &
2484  ) * j33g * j33g * rfdz(ke) &
2485  - ( ( fact_n * (pott(ke,i,j+1)+pott(ke,i,j )) &
2486  + fact_f * (pott(ke,i,j+2)+pott(ke,i,j-1)) ) * rfdy(j ) &
2487  + ( fact_n * (pott(ke,i,j )+pott(ke,i,j-1)) &
2488  + fact_f * (pott(ke,i,j-1)+pott(ke,i,j-2)) ) * rfdy(j-1) &
2489  ) * g(ke,i,j,i_xyz) * rfdy(j) &
2490  - g(ke,i,j,i_xyz) * rcs2t(ke,i,j) * rdt * rdt &
2491  - grav * j33g * 0.5_rp * rcs2t(ke,i,j) * rcdz(ke)
2492  ! k-1
2493  m(2,ke,i,j) = j33g * j33g / g(ke-1,i,j,i_xyw) &
2494  * 0.5_rp * (pott(ke ,i,j)+pott(ke-1,i,j)) &
2495  * rfdz(ke-1) * rcdz(ke) &
2496  - grav * j33g * 0.5_rp * rcs2t(ke,i,j) * rcdz(ke)
2497  enddo
2498 
2499  else
2500  !$omp parallel do
2501  do j = jjs, jje
2502  do i = iis, iie
2503  do k = ks+2, ke-2
2504 #ifdef DEBUG
2505  call check( __line__, pott(k-2,i,j) )
2506  call check( __line__, pott(k-1,i,j) )
2507  call check( __line__, pott(k ,i,j) )
2508  call check( __line__, pott(k+1,i,j) )
2509  call check( __line__, pott(k+2,i,j) )
2510  call check( __line__, pott(k,i-2,j) )
2511  call check( __line__, pott(k,i-1,j) )
2512  call check( __line__, pott(k,i ,j) )
2513  call check( __line__, pott(k,i+1,j) )
2514  call check( __line__, pott(k,i+2,j) )
2515  call check( __line__, pott(k,i,j-2) )
2516  call check( __line__, pott(k,i,j-1) )
2517  call check( __line__, pott(k,i,j ) )
2518  call check( __line__, pott(k,i,j+1) )
2519  call check( __line__, pott(k,i,j+2) )
2520  call check( __line__, g(k-1,i,j,i_xyw) )
2521  call check( __line__, g(k+1,i,j,i_xyw) )
2522  call check( __line__, g(k,i,j,i_xyz) )
2523  call check( __line__, rcs2t(k-1,i,j) )
2524  call check( __line__, rcs2t(k ,i,j) )
2525  call check( __line__, rcs2t(k+1,i,j) )
2526 #endif
2527  ! k,i,j
2528  m(1,k,i,j) = &
2529  - ( ( fact_n * (pott(k+1,i,j)+pott(k ,i,j)) &
2530  + fact_f * (pott(k+2,i,j)+pott(k-1,i,j)) ) * rfdz(k ) / g(k ,i,j,i_xyw) &
2531  + ( fact_n * (pott(k ,i,j)+pott(k-1,i,j)) &
2532  + fact_f * (pott(k+1,i,j)+pott(k-2,i,j)) ) * rfdz(k-1) / g(k-1,i,j,i_xyw) &
2533  ) * j33g * j33g * rfdz(k) &
2534  - ( ( fact_n * (pott(k,i+1,j)+pott(k,i ,j)) &
2535  + fact_f * (pott(k,i+2,j)+pott(k,i-1,j)) ) * rfdx(i ) &
2536  + ( fact_n * (pott(k,i ,j)+pott(k,i-1,j)) &
2537  + fact_f * (pott(k,i+1,j)+pott(k,i-2,j)) ) * rfdx(i-1) &
2538  ) * g(k,i,j,i_xyz) * rfdx(i) &
2539  - ( ( fact_n * (pott(k,i,j+1)+pott(k,i,j )) &
2540  + fact_f * (pott(k,i,j+2)+pott(k,i,j-1)) ) * rfdy(j ) &
2541  + ( fact_n * (pott(k,i,j )+pott(k,i,j-1)) &
2542  + fact_f * (pott(k,i,j-1)+pott(k,i,j-2)) ) * rfdy(j-1) &
2543  ) * g(k,i,j,i_xyz) * rfdy(j) &
2544  - g(k,i,j,i_xyz) * rcs2t(k,i,j) * rdt * rdt
2545  ! k-1
2546  m(2,k,i,j) = j33g * j33g / g(k-1,i,j,i_xyw) &
2547  * ( fact_n * (pott(k ,i,j)+pott(k-1,i,j)) &
2548  + fact_f * (pott(k+1,i,j)+pott(k-2,i,j)) ) &
2549  * rfdz(k-1) * rcdz(k) &
2550  - grav * j33g * rcs2t(k-1,i,j) / ( fdz(k)+fdz(k-1) )
2551  ! k+1
2552  m(3,k,i,j) = j33g * j33g / g(k+1,i,j,i_xyw) &
2553  * ( fact_n * (pott(k+1,i,j)+pott(k ,i,j)) &
2554  + fact_f * (pott(k+2,i,j)+pott(k-1,i,j)) ) &
2555  * rfdz(k ) * rcdz(k) &
2556  + grav * j33g * rcs2t(k+1,i,j) / ( fdz(k)+fdz(k-1) )
2557  enddo
2558  enddo
2559  enddo
2560 
2561  !$omp parallel do
2562  do j = jjs, jje
2563  do i = iis, iie
2564 #ifdef DEBUG
2565  call check( __line__, pott(ks ,i,j) )
2566  call check( __line__, pott(ks+1,i,j) )
2567  call check( __line__, pott(ks,i-2,j) )
2568  call check( __line__, pott(ks,i-1,j) )
2569  call check( __line__, pott(ks,i ,j) )
2570  call check( __line__, pott(ks,i+1,j) )
2571  call check( __line__, pott(ks,i+2,j) )
2572  call check( __line__, pott(ks,i,j-2) )
2573  call check( __line__, pott(ks,i,j-1) )
2574  call check( __line__, pott(ks,i,j ) )
2575  call check( __line__, pott(ks,i,j+1) )
2576  call check( __line__, pott(ks,i,j+2) )
2577  call check( __line__, g(ks,i,j,i_xyz) )
2578  call check( __line__, rcs2t(ks,i,j) )
2579  call check( __line__, pott(ks ,i,j) )
2580  call check( __line__, pott(ks+1,i,j) )
2581  call check( __line__, g(ks+1,i,j,i_xyw) )
2582  call check( __line__, rcs2t(ks+1,i,j) )
2583 #endif
2584  ! k,i,j
2585  m(1,ks,i,j) = &
2586  - ( 0.5_rp * (pott(ks+1,i,j)+pott(ks ,i,j)) * rfdz(ks ) &
2587  ) * j33g * j33g / g(ks ,i,j,i_xyw) * rfdz(ks) &
2588  - ( ( fact_n * (pott(ks,i+1,j)+pott(ks,i ,j)) &
2589  + fact_f * (pott(ks,i+2,j)+pott(ks,i-1,j)) ) * rfdx(i ) &
2590  + ( fact_n * (pott(ks,i ,j)+pott(ks,i-1,j)) &
2591  + fact_f * (pott(ks,i+1,j)+pott(ks,i-2,j)) ) * rfdx(i-1) &
2592  ) * g(ks,i,j,i_xyz) * rfdx(i) &
2593  - ( ( fact_n * (pott(ks,i,j+1)+pott(ks,i,j )) &
2594  + fact_f * (pott(ks,i,j+2)+pott(ks,i,j-1)) ) * rfdy(j ) &
2595  + ( fact_n * (pott(ks,i,j )+pott(ks,i,j-1)) &
2596  + fact_f * (pott(ks,i,j-1)+pott(ks,i,j-2)) ) * rfdy(j-1) &
2597  ) * g(ks,i,j,i_xyz) * rfdy(j) &
2598  - g(ks,i,j,i_xyz) * rcs2t(ks,i,j) * rdt * rdt &
2599  + grav * j33g * 0.5_rp * rcs2t(ks,i,j) * rcdz(ks)
2600  ! k+1
2601  m(3,ks,i,j) = j33g * j33g / g(ks+1,i,j,i_xyw) &
2602  * 0.5_rp * (pott(ks+1,i,j)+pott(ks ,i,j)) &
2603  * rfdz(ks ) * rcdz(ks) &
2604  + grav * j33g * 0.5_rp * rcs2t(ks+1,i,j) * rcdz(ks)
2605 #ifdef DEBUG
2606  call check( __line__, pott(ks ,i,j) )
2607  call check( __line__, pott(ks+1,i,j) )
2608  call check( __line__, pott(ks+2,i,j) )
2609  call check( __line__, pott(ks+3,i,j) )
2610  call check( __line__, pott(ks+1,i-2,j) )
2611  call check( __line__, pott(ks+1,i-1,j) )
2612  call check( __line__, pott(ks+1,i ,j) )
2613  call check( __line__, pott(ks+1,i+1,j) )
2614  call check( __line__, pott(ks+1,i+2,j) )
2615  call check( __line__, pott(ks+1,i,j-2) )
2616  call check( __line__, pott(ks+1,i,j-1) )
2617  call check( __line__, pott(ks+1,i,j ) )
2618  call check( __line__, pott(ks+1,i,j+1) )
2619  call check( __line__, pott(ks+1,i,j+2) )
2620  call check( __line__, g(ks ,i,j,i_xyw) )
2621  call check( __line__, g(ks+2,i,j,i_xyw) )
2622  call check( __line__, g(ks+1,i,j,i_xyz) )
2623  call check( __line__, rcs2t(ks ,i,j) )
2624  call check( __line__, rcs2t(ks+1 ,i,j) )
2625  call check( __line__, rcs2t(ks+2,i,j) )
2626 #endif
2627  ! k,i,j
2628  m(1,ks+1,i,j) = &
2629  - ( ( fact_n * (pott(ks+2,i,j)+pott(ks+1,i,j)) &
2630  + fact_f * (pott(ks+3,i,j)+pott(ks ,i,j)) ) * rfdz(ks+1) / g(ks+1,i,j,i_xyw) &
2631  + 0.5_rp * (pott(ks+1,i,j)+pott(ks ,i,j)) * rfdz(ks ) / g(ks ,i,j,i_xyw) &
2632  ) * j33g * j33g * rfdz(ks+1) &
2633  - ( ( fact_n * (pott(ks+1,i+1,j)+pott(ks+1,i ,j)) &
2634  + fact_f * (pott(ks+1,i+2,j)+pott(ks+1,i-1,j)) ) * rfdx(i ) &
2635  + ( fact_n * (pott(ks+1,i ,j)+pott(ks+1,i-1,j)) &
2636  + fact_f * (pott(ks+1,i+1,j)+pott(ks+1,i-2,j)) ) * rfdx(i-1) &
2637  ) * g(ks+1,i,j,i_xyz) * rfdx(i) &
2638  - ( ( fact_n * (pott(ks+1,i,j+1)+pott(ks+1,i,j )) &
2639  + fact_f * (pott(ks+1,i,j+2)+pott(ks+1,i,j-1)) ) * rfdy(j ) &
2640  + ( fact_n * (pott(ks+1,i,j )+pott(ks+1,i,j-1)) &
2641  + fact_f * (pott(ks+1,i,j-1)+pott(ks+1,i,j-2)) ) * rfdy(j-1) &
2642  ) * g(ks+1,i,j,i_xyz) * rfdy(j) &
2643  - g(ks+1,i,j,i_xyz) * rcs2t(ks+1,i,j) * rdt * rdt
2644  ! k-1
2645  m(2,ks+1,i,j) = j33g * j33g / g(ks,i,j,i_xyw) &
2646  * 0.5_rp * (pott(ks+1,i,j)+pott(ks ,i,j)) &
2647  * rfdz(ks ) * rcdz(ks+1) &
2648  - grav * j33g * rcs2t(ks ,i,j) / ( fdz(ks+1)+fdz(ks) )
2649  ! k+1
2650  m(3,ks+1,i,j) = j33g * j33g / g(ks+2,i,j,i_xyw) &
2651  * ( fact_n * (pott(ks+2,i,j)+pott(ks+1,i,j)) &
2652  + fact_f * (pott(ks+3,i,j)+pott(ks ,i,j)) ) &
2653  * rfdz(ks+1) * rcdz(ks+1) &
2654  + grav * j33g * rcs2t(ks+2,i,j) / ( fdz(ks+1)+fdz(ks) )
2655 #ifdef DEBUG
2656  call check( __line__, pott(ke-3,i,j) )
2657  call check( __line__, pott(ke-2,i,j) )
2658  call check( __line__, pott(ke-1,i,j) )
2659  call check( __line__, pott(ke ,i,j) )
2660  call check( __line__, pott(ke-1,i-2,j) )
2661  call check( __line__, pott(ke-1,i-1,j) )
2662  call check( __line__, pott(ke-1,i ,j) )
2663  call check( __line__, pott(ke-1,i+1,j) )
2664  call check( __line__, pott(ke-1,i+2,j) )
2665  call check( __line__, pott(ke-1,i,j-2) )
2666  call check( __line__, pott(ke-1,i,j-1) )
2667  call check( __line__, pott(ke-1,i,j ) )
2668  call check( __line__, pott(ke-1,i,j+1) )
2669  call check( __line__, pott(ke-1,i,j+2) )
2670  call check( __line__, g(ke-2,i,j,i_xyw) )
2671  call check( __line__, g(ke ,i,j,i_xyw) )
2672  call check( __line__, g(ke-1,i,j,i_xyz) )
2673  call check( __line__, rcs2t(ke-2,i,j) )
2674  call check( __line__, rcs2t(ke-1,i,j) )
2675  call check( __line__, rcs2t(ke ,i,j) )
2676 #endif
2677  ! k,i,j
2678  m(1,ke-1,i,j) = &
2679  - ( 0.5_rp * (pott(ke ,i,j)+pott(ke-1,i,j)) * rfdz(ke-1) / g(ke-1,i,j,i_xyw) &
2680  + ( fact_n * (pott(ke-1,i,j)+pott(ke-2,i,j)) &
2681  + fact_f * (pott(ke ,i,j)+pott(ke-3,i,j)) ) * rfdz(ke-2) / g(ke-2,i,j,i_xyw) &
2682  ) * j33g * j33g * rfdz(ke-1) &
2683  - ( ( fact_n * (pott(ke-1,i+1,j)+pott(ke-1,i ,j)) &
2684  + fact_f * (pott(ke-1,i+2,j)+pott(ke-1,i-1,j)) ) * rfdx(i ) &
2685  + ( fact_n * (pott(ke-1,i ,j)+pott(ke-1,i-1,j)) &
2686  + fact_f * (pott(ke-1,i+1,j)+pott(ke-1,i-2,j)) ) * rfdx(i-1) &
2687  ) * g(ke-1,i,j,i_xyz) * rfdx(i) &
2688  - ( ( fact_n * (pott(ke-1,i,j+1)+pott(ke-1,i,j )) &
2689  + fact_f * (pott(ke-1,i,j+2)+pott(ke-1,i,j-1)) ) * rfdy(j ) &
2690  + ( fact_n * (pott(ke-1,i,j )+pott(ke-1,i,j-1)) &
2691  + fact_f * (pott(ke-1,i,j-1)+pott(ke-1,i,j-2)) ) * rfdy(j-1) &
2692  ) * g(ke-1,i,j,i_xyz) * rfdy(j) &
2693  - g(ke-1,i,j,i_xyz) * rcs2t(ke-1,i,j) * rdt * rdt
2694  ! k-1
2695  m(2,ke-1,i,j) = j33g * j33g / g(ke-2,i,j,i_xyw) &
2696  * ( fact_n * (pott(ke-1,i,j)+pott(ke-2,i,j)) &
2697  + fact_f * (pott(ke ,i,j)+pott(ke-3,i,j)) ) &
2698  * rfdz(ke-2) * rcdz(ke-1) &
2699  - grav * j33g * rcs2t(ke-2,i,j) / ( fdz(ke-1)+fdz(ke-2) )
2700  ! k+1
2701  m(3,ke-1,i,j) = j33g * j33g / g(ke ,i,j,i_xyw) &
2702  * 0.5_rp * (pott(ke ,i,j)+pott(ke-1,i,j)) &
2703  * rfdz(ke-1) * rcdz(ke-1) &
2704  + grav * j33g * rcs2t(ke ,i,j) / ( fdz(ke-1)+fdz(ke-2) )
2705 #ifdef DEBUG
2706  call check( __line__, pott(ke-1,i,j) )
2707  call check( __line__, pott(ke ,i,j) )
2708  call check( __line__, pott(ke,i-2,j) )
2709  call check( __line__, pott(ke,i-1,j) )
2710  call check( __line__, pott(ke,i ,j) )
2711  call check( __line__, pott(ke,i+1,j) )
2712  call check( __line__, pott(ke,i+2,j) )
2713  call check( __line__, pott(ke,i,j-2) )
2714  call check( __line__, pott(ke,i,j-1) )
2715  call check( __line__, pott(ke,i,j ) )
2716  call check( __line__, pott(ke,i,j+1) )
2717  call check( __line__, pott(ke,i,j+2) )
2718  call check( __line__, g(ke-1,i,j,i_xyw) )
2719  call check( __line__, g(ke,i,j,i_xyz) )
2720  call check( __line__, rcs2t(ke-1,i,j) )
2721  call check( __line__, rcs2t(ke,i,j) )
2722 #endif
2723  ! k,i,j
2724  m(1,ke,i,j) = &
2725  - ( &
2726  + 0.5_rp * (pott(ke ,i,j)+pott(ke-1,i,j)) * rfdz(ke-1) / g(ke-1,i,j,i_xyw) &
2727  ) * j33g * j33g * rfdz(ke) &
2728  - ( ( fact_n * (pott(ke,i+1,j)+pott(ke,i ,j)) &
2729  + fact_f * (pott(ke,i+2,j)+pott(ke,i-1,j)) ) * rfdx(i ) &
2730  + ( fact_n * (pott(ke,i ,j)+pott(ke,i-1,j)) &
2731  + fact_f * (pott(ke,i+1,j)+pott(ke,i-2,j)) ) * rfdx(i-1) &
2732  ) * g(ke,i,j,i_xyz) * rfdx(i) &
2733  - ( ( fact_n * (pott(ke,i,j+1)+pott(ke,i,j )) &
2734  + fact_f * (pott(ke,i,j+2)+pott(ke,i,j-1)) ) * rfdy(j ) &
2735  + ( fact_n * (pott(ke,i,j )+pott(ke,i,j-1)) &
2736  + fact_f * (pott(ke,i,j-1)+pott(ke,i,j-2)) ) * rfdy(j-1) &
2737  ) * g(ke,i,j,i_xyz) * rfdy(j) &
2738  - g(ke,i,j,i_xyz) * rcs2t(ke,i,j) * rdt * rdt &
2739  - grav * j33g * 0.5_rp * rcs2t(ke,i,j) * rcdz(ke)
2740  ! k-1
2741  m(2,ke,i,j) = j33g * j33g / g(ke-1,i,j,i_xyw) &
2742  * 0.5_rp * (pott(ke ,i,j)+pott(ke-1,i,j)) &
2743  * rfdz(ke-1) * rcdz(ke) &
2744  - grav * j33g * 0.5_rp * rcs2t(ke,i,j) * rcdz(ke)
2745  enddo
2746  enddo
2747 
2748  end if
2749 #ifdef DEBUG
2750  k = iundef; i = iundef; j = iundef
2751 #endif
2752 
2753  if ( twod ) then
2754  i = is
2755  !$omp parallel do
2756  do j = jjs, jje
2757  do k = ks, ke
2758 #ifdef DEBUG
2759  call check( __line__, g(k,i,j-1,i_xyz) )
2760  call check( __line__, g(k,i,j+1,i_xyz) )
2761  call check( __line__, pott(k,i ,j ) )
2762  call check( __line__, pott(k,i ,j-2) )
2763  call check( __line__, pott(k,i ,j-1) )
2764  call check( __line__, pott(k,i ,j ) )
2765  call check( __line__, pott(k,i ,j+1) )
2766  call check( __line__, pott(k,i ,j+2) )
2767 #endif
2768  ! i-1
2769  !M(4,k,i,j) = 0.0_RP
2770  ! i+1
2771  !M(5,k,i,j) = 0.0_RP
2772  ! j-1
2773  m(6,k,i,j) = g(k,i,j-1,i_xyz) &
2774  * ( fact_n * (pott(k,i,j )+pott(k,i,j-1)) &
2775  + fact_f * (pott(k,i,j+1)+pott(k,i,j-2)) ) &
2776  * rfdy(j-1) * rcdy(j)
2777  ! j+1
2778  m(7,k,i,j) = g(k,i,j+1,i_xyz) &
2779  * ( fact_n * (pott(k,i,j+1)+pott(k,i,j )) &
2780  + fact_f * (pott(k,i,j+2)+pott(k,i,j-1)) ) &
2781  * rfdy(j ) * rcdy(j)
2782  enddo
2783  enddo
2784  else
2785  !$omp parallel do
2786  do j = jjs, jje
2787  do i = iis, iie
2788  do k = ks, ke
2789 #ifdef DEBUG
2790  call check( __line__, g(k,i-1,j,i_xyz) )
2791  call check( __line__, g(k,i+1,j,i_xyz) )
2792  call check( __line__, g(k,i,j-1,i_xyz) )
2793  call check( __line__, g(k,i,j+1,i_xyz) )
2794  call check( __line__, pott(k,i-2,j ) )
2795  call check( __line__, pott(k,i-1,j ) )
2796  call check( __line__, pott(k,i ,j ) )
2797  call check( __line__, pott(k,i+1,j ) )
2798  call check( __line__, pott(k,i+2,j ) )
2799  call check( __line__, pott(k,i ,j-2) )
2800  call check( __line__, pott(k,i ,j-1) )
2801  call check( __line__, pott(k,i ,j ) )
2802  call check( __line__, pott(k,i ,j+1) )
2803  call check( __line__, pott(k,i ,j+2) )
2804 #endif
2805  ! i-1
2806  m(4,k,i,j) = g(k,i-1,j,i_xyz) &
2807  * ( fact_n * (pott(k,i ,j)+pott(k,i-1,j)) &
2808  + fact_f * (pott(k,i+1,j)+pott(k,i-2,j)) ) &
2809  * rfdx(i-1) * rcdx(i)
2810  ! i+1
2811  m(5,k,i,j) = g(k,i+1,j,i_xyz) &
2812  * ( fact_n * (pott(k,i+1,j)+pott(k,i ,j)) &
2813  + fact_f * (pott(k,i+2,j)+pott(k,i-1,j)) ) &
2814  * rfdx(i ) * rcdx(i)
2815  ! j-1
2816  m(6,k,i,j) = g(k,i,j-1,i_xyz) &
2817  * ( fact_n * (pott(k,i,j )+pott(k,i,j-1)) &
2818  + fact_f * (pott(k,i,j+1)+pott(k,i,j-2)) ) &
2819  * rfdy(j-1) * rcdy(j)
2820  ! j+1
2821  m(7,k,i,j) = g(k,i,j+1,i_xyz) &
2822  * ( fact_n * (pott(k,i,j+1)+pott(k,i,j )) &
2823  + fact_f * (pott(k,i,j+2)+pott(k,i,j-1)) ) &
2824  * rfdy(j ) * rcdy(j)
2825  enddo
2826  enddo
2827  enddo
2828  end if
2829 #ifdef DEBUG
2830  k = iundef; i = iundef; j = iundef
2831 #endif
2832 
2833  return

References scale_debug::check(), scale_atmos_grid_cartesc_index::ie, scale_atmos_grid_cartesc_index::is, scale_atmos_grid_cartesc_index::je, scale_atmos_grid_cartesc_index::js, scale_atmos_grid_cartesc_index::ke, scale_atmos_grid_cartesc_index::ks, and scale_prc::prc_abort().

Referenced by atmos_dyn_tstep_short_fvm_hivi().

Here is the call graph for this function:
Here is the caller graph for this function:
Generated by doxygen 1.9.1