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, dtrk, last)
 
subroutine solve_bicgstab (DPRES_N, DPRES, M, B)
 
subroutine make_matrix (M, POTT, RCs2T, GRAV, G, J33G, RCDZ, RFDZ, RCDX, RFDX, RCDY, RFDY, FDZ, rdt, FACT_N, FACT_F, 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 79 of file scale_atmos_dyn_tstep_short_fvm_hivi.F90.

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

Referenced by scale_atmos_dyn_tstep_short::atmos_dyn_tstep_short_regist().

79  use scale_process, only: &
80  prc_mpistop
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  write(*,*) 'xxx ATMOS_DYN_TYPE is not FVM-HIVI. Check!'
92  call prc_mpistop
93  endif
94 
95  va_out = va_fvm_hivi
96  var_name(:) = ""
97  var_desc(:) = ""
98  var_unit(:) = ""
99 
100  if( io_l ) write(io_fid_log,*)
101  if( io_l ) write(io_fid_log,*) '*** Register additional prognostic variables (HIVI)'
102  if ( va_out < 1 ) then
103  if( io_l ) write(io_fid_log,*) '*** => nothing.'
104  endif
105 
106  return
scale_process::prc_mpistop
subroutine, public prc_mpistop
Abort MPI.
Definition: scale_process.F90:267
scale_process
module PROCESS
Definition: scale_process.F90:14
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 112 of file scale_atmos_dyn_tstep_short_fvm_hivi.F90.

References scale_precision::dp, scale_stdio::io_fid_conf, scale_stdio::io_fid_log, scale_stdio::io_fid_nml, scale_stdio::io_l, scale_stdio::io_nml, scale_process::prc_mpistop(), scale_precision::rp, and scale_precision::sp.

Referenced by scale_atmos_dyn_tstep_short::atmos_dyn_tstep_short_regist().

112  use scale_process, only: &
113  prc_mpistop
114  implicit none
115 
116  namelist / param_atmos_dyn_tstep_fvm_hivi / &
117  itmax, &
118  epsilon
119 
120  integer :: ierr
121  !---------------------------------------------------------------------------
122 
123  if( io_l ) write(io_fid_log,*) '*** HIVI Setup'
124 #ifdef HIVI_BICGSTAB
125  if( io_l ) write(io_fid_log,*) '*** USING Bi-CGSTAB'
126 #else
127  if( io_l ) write(io_fid_log,*) '*** USING Multi-Grid'
128  write(*,*) 'xxx Not Implemented yet'
129  call prc_mpistop
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  if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
142  elseif( ierr > 0 ) then !--- fatal error
143  write(*,*) 'xxx Not appropriate names in namelist PARAM_ATMOS_DYN_TSTEP_FVM_HIVI. Check!'
144  call prc_mpistop
145  endif
146  if( io_nml ) write(io_fid_nml,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  write(*,*) 'xxx Unsupported precision'
154  call prc_mpistop
155  endif
156 
157  return
scale_process::prc_mpistop
subroutine, public prc_mpistop
Abort MPI.
Definition: scale_process.F90:267
scale_process
module PROCESS
Definition: scale_process.F90:14
scale_precision::sp
integer, parameter, public sp
Definition: scale_precision.F90:33
scale_precision::rp
integer, parameter, public rp
Definition: scale_precision.F90:43
scale_precision::dp
integer, parameter, public dp
Definition: scale_precision.F90:34
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(1,ia,ja), intent(in)  CORIOLI,
real(rp), dimension(ka,ia,ja,5,3), intent(in)  num_diff,
real(rp), dimension(ka), intent(in)  wdamp_coef,
real(rp), intent(in)  divdmp_coef,
real(rp), dimension(ka,ia,ja), intent(in)  DDIV,
logical, intent(in)  FLAG_FCT_MOMENTUM,
logical, intent(in)  FLAG_FCT_T,
logical, intent(in)  FLAG_FCT_ALONG_STREAM,
real(rp), dimension(ka), intent(in)  CDZ,
real(rp), dimension(ka-1), intent(in)  FDZ,
real(rp), dimension(ia-1), intent(in)  FDX,
real(rp), dimension(ja-1), intent(in)  FDY,
real(rp), dimension(ka), intent(in)  RCDZ,
real(rp), dimension(ia), intent(in)  RCDX,
real(rp), dimension(ja), intent(in)  RCDY,
real(rp), dimension(ka-1), intent(in)  RFDZ,
real(rp), dimension(ia-1), intent(in)  RFDX,
real(rp), dimension(ja-1), intent(in)  RFDY,
real(rp), dimension (ka,ia,ja), intent(in)  PHI,
real(rp), dimension (ka,ia,ja,7), intent(in)  GSQRT,
real(rp), dimension (ka,ia,ja,7), intent(in)  J13G,
real(rp), dimension (ka,ia,ja,7), intent(in)  J23G,
real(rp), intent(in)  J33G,
real(rp), dimension (ia,ja,2,4), intent(in)  MAPF,
real(rp), dimension(ka,ia,ja), intent(in)  REF_dens,
real(rp), dimension(ka,ia,ja), intent(in)  REF_rhot,
logical, intent(in)  BND_W,
logical, intent(in)  BND_E,
logical, intent(in)  BND_S,
logical, intent(in)  BND_N,
real(rp), intent(in)  dtrk,
logical, intent(in)  last 
)
Parameters
[in]phigeopotential
[in]gsqrtvertical metrics {G}^1/2
[in]j13g(1,3) element of Jacobian matrix
[in]j23g(2,3) element of Jacobian matrix
[in]j33g(3,3) element of Jacobian matrix
[in]mapfmap factor
[in]ref_densreference density

Definition at line 179 of file scale_atmos_dyn_tstep_short_fvm_hivi.F90.

References scale_atmos_dyn_common::atmos_dyn_fct(), scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_flux_valuew_z, scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxj13_uyz, scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxj13_xvz, scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxj13_xyw, scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxj23_uyz, scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxj23_xvz, scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxj23_xyw, scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxx_uyz, scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxx_xvz, scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxx_xyw, scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxx_xyz, scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxy_uyz, scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxy_xvz, scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxy_xyw, scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxy_xyz, scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxz_uyz, scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxz_xvz, scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxz_xyw, scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxz_xyz, scale_debug::check(), scale_const::const_grav, scale_const::const_pre00, scale_index::i_dens, scale_index::i_momx, scale_index::i_momy, scale_index::i_momz, scale_index::i_rhot, scale_gridtrans::i_uv, scale_gridtrans::i_uvz, scale_gridtrans::i_uy, scale_gridtrans::i_uyw, scale_gridtrans::i_uyz, scale_gridtrans::i_xv, scale_gridtrans::i_xvw, scale_gridtrans::i_xvz, scale_gridtrans::i_xy, scale_gridtrans::i_xyw, scale_gridtrans::i_xyz, scale_grid_index::iblock, scale_grid_index::ie, scale_grid_index::is, scale_grid_index::jblock, scale_grid_index::je, scale_grid_index::js, scale_grid_index::ka, scale_grid_index::ke, scale_grid_index::ks, make_matrix(), solve_bicgstab(), scale_grid_index::xdir, scale_grid_index::ydir, and scale_grid_index::zdir.

Referenced by scale_atmos_dyn_tstep_short::atmos_dyn_tstep_short_regist().

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

Definition at line 1477 of file scale_atmos_dyn_tstep_short_fvm_hivi.F90.

References scale_debug::check(), scale_comm::comm_world, scale_grid_index::iblock, scale_grid_index::ie, scale_stdio::io_l, scale_grid_index::is, scale_grid_index::jblock, scale_grid_index::je, scale_grid_index::js, scale_grid_index::ke, scale_grid_index::ks, and scale_process::prc_mpistop().

Referenced by atmos_dyn_tstep_short_fvm_hivi().

1477  use scale_process, only: &
1478  prc_mpistop
1479  use scale_comm, only: &
1480  comm_world, &
1481  comm_vars8, &
1482  comm_wait
1483  implicit none
1484  real(RP), intent(out) :: DPRES_N(KA,IA,JA)
1485  real(RP), intent(in) :: DPRES(KA,IA,JA)
1486  real(RP), intent(in) :: M(7,KA,IA,JA)
1487  real(RP), intent(in) :: B(KA,IA,JA)
1488 
1489  real(RP) :: r0(KA,IA,JA)
1490 
1491  real(RP) :: p(KA,IA,JA)
1492  real(RP) :: Mp(KA,IA,JA)
1493  real(RP) :: s(KA,IA,JA)
1494  real(RP) :: Ms(KA,IA,JA)
1495  real(RP) :: al, be, w
1496 
1497  real(RP), pointer :: r(:,:,:)
1498  real(RP), pointer :: rn(:,:,:)
1499  real(RP), pointer :: swap(:,:,:)
1500  real(RP), target :: v0(KA,IA,JA)
1501  real(RP), target :: v1(KA,IA,JA)
1502  real(RP) :: r0r
1503  real(RP) :: norm, error, error2
1504 
1505  real(RP) :: iprod(2)
1506  real(RP) :: buf(2)
1507 
1508  integer :: k, i, j
1509  integer :: iis, iie, jjs, jje
1510  integer :: iter
1511  integer :: ierror
1512 
1513 #ifdef DEBUG
1514  r0(:,:,:) = undef
1515  p(:,:,:) = undef
1516  mp(:,:,:) = undef
1517  s(:,:,:) = undef
1518  ms(:,:,:) = undef
1519 
1520  v0(:,:,:) = undef
1521  v1(:,:,:) = undef
1522 #endif
1523 
1524  r => v0
1525  rn => v1
1526 
1527  call mul_matrix( v1, m, dpres ) ! v1 = M x0
1528 
1529  norm = 0.0_rp
1530  r0r = 0.0_rp
1531 
1532  do jjs = js, je, jblock
1533  jje = jjs+jblock-1
1534  do iis = is, ie, iblock
1535  iie = iis+iblock-1
1536 
1537  do j = jjs, jje
1538  do i = iis, iie
1539  do k = ks, ke
1540  norm = norm + b(k,i,j)**2
1541  enddo
1542  enddo
1543  enddo
1544 #ifdef DEBUG
1545  k = iundef; i = iundef; j = iundef
1546 #endif
1547 
1548  ! r = b - M x0
1549  do j = jjs, jje
1550  do i = iis, iie
1551  do k = ks, ke
1552 #ifdef DEBUG
1553  call check( __line__, b(k,i,j) )
1554  call check( __line__, v1(k,i,j) )
1555 #endif
1556  r(k,i,j) = b(k,i,j) - v1(k,i,j)
1557  enddo
1558  enddo
1559  enddo
1560 #ifdef DEBUG
1561  k = iundef; i = iundef; j = iundef
1562 #endif
1563 
1564  do j = jjs, jje
1565  do i = iis, iie
1566  do k = ks, ke
1567  r0(k,i,j) = r(k,i,j)
1568  p(k,i,j) = r(k,i,j)
1569  enddo
1570  enddo
1571  enddo
1572 #ifdef DEBUG
1573  k = iundef; i = iundef; j = iundef
1574 #endif
1575 
1576  do j = jjs, jje
1577  do i = iis, iie
1578  do k = ks, ke
1579 #ifdef DEBUG
1580  call check( __line__, r(k,i,j) )
1581  call check( __line__, r0(k,i,j) )
1582 #endif
1583  r0r = r0r + r0(k,i,j) * r(k,i,j)
1584  enddo
1585  enddo
1586  enddo
1587 #ifdef DEBUG
1588  k = iundef; i = iundef; j = iundef
1589 #endif
1590 
1591  enddo
1592  enddo
1593 
1594  do j = js-1, je+1
1595  do i = is-1, ie+1
1596  do k = ks, ke
1597  dpres_n(k,i,j) = dpres(k,i,j)
1598  end do
1599  end do
1600  end do
1601 
1602 
1603  iprod(1) = r0r
1604  iprod(2) = norm
1605  call mpi_allreduce(iprod, buf, 2, mtype, mpi_sum, comm_world, ierror)
1606  r0r = buf(1)
1607  norm = buf(2)
1608 
1609  error2 = norm
1610 
1611  do iter = 1, itmax
1612 
1613  error = 0.0_rp
1614  do j = js, je
1615  do i = is, ie
1616  do k = ks, ke
1617 #ifdef DEBUG
1618  call check( __line__, r(k,i,j) )
1619 #endif
1620  error = error + r(k,i,j)**2
1621  enddo
1622  enddo
1623  enddo
1624 #ifdef DEBUG
1625  k = iundef; i = iundef; j = iundef
1626 #endif
1627  call mpi_allreduce(error, buf, 1, mtype, mpi_sum, comm_world, ierror)
1628  error = buf(1)
1629 
1630 #ifdef DEBUG
1631  if( io_l ) write(*,*) iter, error/norm
1632 #endif
1633  if ( sqrt(error/norm) < epsilon .OR. error > error2 ) then
1634 #ifdef DEBUG
1635  if( io_l ) write(*,*) "Bi-CGSTAB converged:", iter
1636 #endif
1637  exit
1638  endif
1639  error2 = error
1640 
1641  call comm_vars8( p, 1 )
1642  call comm_wait ( p, 1 )
1643  call mul_matrix( mp, m, p )
1644 
1645  iprod(1) = 0.0_rp
1646  do j = js, je
1647  do i = is, ie
1648  do k = ks, ke
1649 #ifdef DEBUG
1650  call check( __line__, r0(k,i,j) )
1651  call check( __line__, mp(k,i,j) )
1652 #endif
1653  iprod(1) = iprod(1) + r0(k,i,j) * mp(k,i,j)
1654  enddo
1655  enddo
1656  enddo
1657 #ifdef DEBUG
1658  k = iundef; i = iundef; j = iundef
1659 #endif
1660  call mpi_allreduce(iprod, buf, 1, mtype, mpi_sum, comm_world, ierror)
1661  al = r0r / buf(1) ! (r0,r) / (r0,Mp)
1662 
1663  do j = js, je
1664  do i = is, ie
1665  do k = ks, ke
1666 #ifdef DEBUG
1667  call check( __line__, r(k,i,j) )
1668  call check( __line__, mp(k,i,j) )
1669 #endif
1670  s(k,i,j) = r(k,i,j) - al*mp(k,i,j)
1671  enddo
1672  enddo
1673  enddo
1674 #ifdef DEBUG
1675  k = iundef; i = iundef; j = iundef
1676 #endif
1677 
1678  call comm_vars8( s, 1 )
1679  call comm_wait ( s, 1 )
1680  call mul_matrix( ms, m, s )
1681  iprod(1) = 0.0_rp
1682  iprod(2) = 0.0_rp
1683  do j = js, je
1684  do i = is, ie
1685  do k = ks, ke
1686 #ifdef DEBUG
1687  call check( __line__, ms(k,i,j) )
1688  call check( __line__, s(k,i,j) )
1689 #endif
1690  iprod(1) = iprod(1) + ms(k,i,j) * s(k,i,j)
1691  iprod(2) = iprod(2) + ms(k,i,j) * ms(k,i,j)
1692  enddo
1693  enddo
1694  enddo
1695 #ifdef DEBUG
1696  k = iundef; i = iundef; j = iundef
1697 #endif
1698  call mpi_allreduce(iprod, buf, 2, mtype, mpi_sum, comm_world, ierror)
1699  w = buf(1) / buf(2) ! (Ms,s) / (Ms,Ms)
1700 
1701  iprod(1) = 0.0_rp
1702 
1703  do jjs = js, je, jblock
1704  jje = jjs+jblock-1
1705  do iis = is, ie, iblock
1706  iie = iis+iblock-1
1707 
1708  do j = jjs, jje
1709  do i = iis, iie
1710  do k = ks, ke
1711 #ifdef DEBUG
1712  call check( __line__, dpres_n(k,i,j) )
1713  call check( __line__, p(k,i,j) )
1714  call check( __line__, s(k,i,j) )
1715 #endif
1716  dpres_n(k,i,j) = dpres_n(k,i,j) + al*p(k,i,j) + w*s(k,i,j)
1717  enddo
1718  enddo
1719  enddo
1720 #ifdef DEBUG
1721  k = iundef; i = iundef; j = iundef
1722 #endif
1723 
1724  do j = jjs, jje
1725  do i = iis, iie
1726  do k = ks, ke
1727 #ifdef DEBUG
1728  call check( __line__, s(k,i,j) )
1729  call check( __line__, ms(k,i,j) )
1730 #endif
1731  rn(k,i,j) = s(k,i,j) - w*ms(k,i,j)
1732  enddo
1733  enddo
1734  enddo
1735 #ifdef DEBUG
1736  k = iundef; i = iundef; j = iundef
1737 #endif
1738 
1739  do j = jjs, jje
1740  do i = iis, iie
1741  do k = ks, ke
1742 #ifdef DEBUG
1743  call check( __line__, r0(k,i,j) )
1744  call check( __line__, rn(k,i,j) )
1745 #endif
1746  iprod(1) = iprod(1) + r0(k,i,j) * rn(k,i,j)
1747  enddo
1748  enddo
1749  enddo
1750 #ifdef DEBUG
1751  k = iundef; i = iundef; j = iundef
1752 #endif
1753 
1754  enddo
1755  enddo
1756 
1757  be = al/w / r0r
1758 
1759  call mpi_allreduce(iprod, r0r, 1, mtype, mpi_sum, comm_world, ierror)
1760 
1761  be = be * r0r ! al/w * (r0,rn)/(r0,r)
1762  do j = js, je
1763  do i = is, ie
1764  do k = ks, ke
1765 #ifdef DEBUG
1766  call check( __line__, rn(k,i,j) )
1767  call check( __line__, p(k,i,j) )
1768  call check( __line__, mp(k,i,j) )
1769 #endif
1770  p(k,i,j) = rn(k,i,j) + be * ( p(k,i,j) - w*mp(k,i,j) )
1771  enddo
1772  enddo
1773  enddo
1774 #ifdef DEBUG
1775  k = iundef; i = iundef; j = iundef
1776 #endif
1777 
1778  swap => rn
1779  rn => r
1780  r => swap
1781 #ifdef DEBUG
1782  rn(:,:,:) = undef
1783 #endif
1784  enddo
1785 
1786  if ( iter >= itmax ) then
1787  write(*,*) 'xxx [atmos_dyn_hivi] Bi-CGSTAB'
1788  write(*,*) 'xxx not converged', error, norm
1789  call prc_mpistop
1790  endif
1791 
1792  return
scale_grid_index::is
integer, public is
start point of inner domain: x, local
Definition: scale_grid_index.F90:60
scale_grid_index::je
integer, public je
end point of inner domain: y, local
Definition: scale_grid_index.F90:63
scale_process::prc_mpistop
subroutine, public prc_mpistop
Abort MPI.
Definition: scale_process.F90:267
scale_grid_index::iblock
integer, public iblock
block size for cache blocking: x
Definition: scale_grid_index.F90:47
scale_grid_index::ke
integer, public ke
end point of inner domain: z, local
Definition: scale_grid_index.F90:59
scale_comm::comm_world
integer, public comm_world
communication world ID
Definition: scale_comm.F90:118
scale_grid_index::jblock
integer, public jblock
block size for cache blocking: y
Definition: scale_grid_index.F90:48
scale_comm
module COMMUNICATION
Definition: scale_comm.F90:23
scale_grid_index::js
integer, public js
start point of inner domain: y, local
Definition: scale_grid_index.F90:62
scale_process
module PROCESS
Definition: scale_process.F90:14
scale_grid_index::ks
integer, public ks
start point of inner domain: z, local
Definition: scale_grid_index.F90:58
scale_grid_index::ie
integer, public ie
end point of inner domain: x, local
Definition: scale_grid_index.F90:61
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,
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 1804 of file scale_atmos_dyn_tstep_short_fvm_hivi.F90.

References scale_debug::check(), scale_grid_index::ie, scale_grid_index::is, scale_grid_index::je, scale_grid_index::js, scale_grid_index::ke, and scale_grid_index::ks.

Referenced by atmos_dyn_tstep_short_fvm_hivi().

1804  implicit none
1805  real(RP), intent(inout) :: M(7,KA,IA,JA)
1806  real(RP), intent(in) :: POTT(KA,IA,JA)
1807  real(RP), intent(in) :: RCs2T(KA,IA,JA)
1808  real(RP), intent(in) :: GRAV
1809  real(RP), intent(in) :: G(KA,IA,JA,7)
1810  real(RP), intent(in) :: J33G
1811  real(RP), intent(in) :: RCDZ(KA)
1812  real(RP), intent(in) :: RFDZ(KA-1)
1813  real(RP), intent(in) :: RFDX(IA-1)
1814  real(RP), intent(in) :: RCDX(IA)
1815  real(RP), intent(in) :: RCDY(JA)
1816  real(RP), intent(in) :: RFDY(JA-1)
1817  real(RP), intent(in) :: FDZ(KA-1)
1818  real(RP), intent(in) :: rdt
1819  real(RP), intent(in) :: FACT_N
1820  real(RP), intent(in) :: FACT_F
1821  integer, intent(in) :: I_XYZ
1822  integer, intent(in) :: I_XYW
1823  integer, intent(in) :: IIS
1824  integer, intent(in) :: IIE
1825  integer, intent(in) :: JJS
1826  integer, intent(in) :: JJE
1827 
1828  integer :: k, i, j
1829 
1830  do j = jjs, jje
1831  do i = iis, iie
1832  do k = ks+2, ke-2
1833 #ifdef DEBUG
1834  call check( __line__, pott(k-2,i,j) )
1835  call check( __line__, pott(k-1,i,j) )
1836  call check( __line__, pott(k ,i,j) )
1837  call check( __line__, pott(k+1,i,j) )
1838  call check( __line__, pott(k+2,i,j) )
1839  call check( __line__, pott(k,i-2,j) )
1840  call check( __line__, pott(k,i-1,j) )
1841  call check( __line__, pott(k,i ,j) )
1842  call check( __line__, pott(k,i+1,j) )
1843  call check( __line__, pott(k,i+2,j) )
1844  call check( __line__, pott(k,i,j-2) )
1845  call check( __line__, pott(k,i,j-1) )
1846  call check( __line__, pott(k,i,j ) )
1847  call check( __line__, pott(k,i,j+1) )
1848  call check( __line__, pott(k,i,j+2) )
1849  call check( __line__, g(k-1,i,j,i_xyw) )
1850  call check( __line__, g(k+1,i,j,i_xyw) )
1851  call check( __line__, g(k,i,j,i_xyz) )
1852  call check( __line__, rcs2t(k-1,i,j) )
1853  call check( __line__, rcs2t(k ,i,j) )
1854  call check( __line__, rcs2t(k+1,i,j) )
1855 #endif
1856  ! k,i,j
1857  m(1,k,i,j) = &
1858  - ( ( fact_n * (pott(k+1,i,j)+pott(k ,i,j)) &
1859  + fact_f * (pott(k+2,i,j)+pott(k-1,i,j)) ) * rfdz(k ) / g(k ,i,j,i_xyw) &
1860  + ( fact_n * (pott(k ,i,j)+pott(k-1,i,j)) &
1861  + fact_f * (pott(k+1,i,j)+pott(k-2,i,j)) ) * rfdz(k-1) / g(k-1,i,j,i_xyw) &
1862  ) * j33g * j33g * rfdz(k) &
1863  - ( ( fact_n * (pott(k,i+1,j)+pott(k,i ,j)) &
1864  + fact_f * (pott(k,i+2,j)+pott(k,i-1,j)) ) * rfdx(i ) &
1865  + ( fact_n * (pott(k,i ,j)+pott(k,i-1,j)) &
1866  + fact_f * (pott(k,i+1,j)+pott(k,i-2,j)) ) * rfdx(i-1) &
1867  ) * g(k,i,j,i_xyz) * rfdx(i) &
1868  - ( ( fact_n * (pott(k,i,j+1)+pott(k,i,j )) &
1869  + fact_f * (pott(k,i,j+2)+pott(k,i,j-1)) ) * rfdy(j ) &
1870  + ( fact_n * (pott(k,i,j )+pott(k,i,j-1)) &
1871  + fact_f * (pott(k,i,j-1)+pott(k,i,j-2)) ) * rfdy(j-1) &
1872  ) * g(k,i,j,i_xyz) * rfdy(j) &
1873  - g(k,i,j,i_xyz) * rcs2t(k,i,j) * rdt * rdt
1874  ! k-1
1875  m(2,k,i,j) = j33g * j33g / g(k-1,i,j,i_xyw) &
1876  * ( fact_n * (pott(k ,i,j)+pott(k-1,i,j)) &
1877  + fact_f * (pott(k+1,i,j)+pott(k-2,i,j)) ) &
1878  * rfdz(k-1) * rcdz(k) &
1879  - grav * j33g * rcs2t(k-1,i,j) / ( fdz(k)+fdz(k-1) )
1880  ! k+1
1881  m(3,k,i,j) = j33g * j33g / g(k+1,i,j,i_xyw) &
1882  * ( fact_n * (pott(k+1,i,j)+pott(k ,i,j)) &
1883  + fact_f * (pott(k+2,i,j)+pott(k-1,i,j)) ) &
1884  * rfdz(k ) * rcdz(k) &
1885  + grav * j33g * rcs2t(k+1,i,j) / ( fdz(k)+fdz(k-1) )
1886  enddo
1887  enddo
1888  enddo
1889 #ifdef DEBUG
1890  k = iundef; i = iundef; j = iundef
1891 #endif
1892  do j = jjs, jje
1893  do i = iis, iie
1894 #ifdef DEBUG
1895  call check( __line__, pott(ks ,i,j) )
1896  call check( __line__, pott(ks+1,i,j) )
1897  call check( __line__, pott(ks,i-2,j) )
1898  call check( __line__, pott(ks,i-1,j) )
1899  call check( __line__, pott(ks,i ,j) )
1900  call check( __line__, pott(ks,i+1,j) )
1901  call check( __line__, pott(ks,i+2,j) )
1902  call check( __line__, pott(ks,i,j-2) )
1903  call check( __line__, pott(ks,i,j-1) )
1904  call check( __line__, pott(ks,i,j ) )
1905  call check( __line__, pott(ks,i,j+1) )
1906  call check( __line__, pott(ks,i,j+2) )
1907  call check( __line__, g(ks,i,j,i_xyz) )
1908  call check( __line__, rcs2t(ks,i,j) )
1909  call check( __line__, pott(ks ,i,j) )
1910  call check( __line__, pott(ks+1,i,j) )
1911  call check( __line__, g(ks+1,i,j,i_xyw) )
1912  call check( __line__, rcs2t(ks+1,i,j) )
1913 #endif
1914  ! k,i,j
1915  m(1,ks,i,j) = &
1916  - ( 0.5_rp * (pott(ks+1,i,j)+pott(ks ,i,j)) * rfdz(ks ) &
1917  ) * j33g * j33g / g(ks ,i,j,i_xyw) * rfdz(ks) &
1918  - ( ( fact_n * (pott(ks,i+1,j)+pott(ks,i ,j)) &
1919  + fact_f * (pott(ks,i+2,j)+pott(ks,i-1,j)) ) * rfdx(i ) &
1920  + ( fact_n * (pott(ks,i ,j)+pott(ks,i-1,j)) &
1921  + fact_f * (pott(ks,i+1,j)+pott(ks,i-2,j)) ) * rfdx(i-1) &
1922  ) * g(ks,i,j,i_xyz) * rfdx(i) &
1923  - ( ( fact_n * (pott(ks,i,j+1)+pott(ks,i,j )) &
1924  + fact_f * (pott(ks,i,j+2)+pott(ks,i,j-1)) ) * rfdy(j ) &
1925  + ( fact_n * (pott(ks,i,j )+pott(ks,i,j-1)) &
1926  + fact_f * (pott(ks,i,j-1)+pott(ks,i,j-2)) ) * rfdy(j-1) &
1927  ) * g(ks,i,j,i_xyz) * rfdy(j) &
1928  - g(ks,i,j,i_xyz) * rcs2t(ks,i,j) * rdt * rdt &
1929  + grav * j33g * 0.5_rp * rcs2t(ks,i,j) * rcdz(ks)
1930  ! k+1
1931  m(3,ks,i,j) = j33g * j33g / g(ks+1,i,j,i_xyw) &
1932  * 0.5_rp * (pott(ks+1,i,j)+pott(ks ,i,j)) &
1933  * rfdz(ks ) * rcdz(ks) &
1934  + grav * j33g * 0.5_rp * rcs2t(ks+1,i,j) * rcdz(ks)
1935 #ifdef DEBUG
1936  call check( __line__, pott(ks ,i,j) )
1937  call check( __line__, pott(ks+1,i,j) )
1938  call check( __line__, pott(ks+2,i,j) )
1939  call check( __line__, pott(ks+3,i,j) )
1940  call check( __line__, pott(ks+1,i-2,j) )
1941  call check( __line__, pott(ks+1,i-1,j) )
1942  call check( __line__, pott(ks+1,i ,j) )
1943  call check( __line__, pott(ks+1,i+1,j) )
1944  call check( __line__, pott(ks+1,i+2,j) )
1945  call check( __line__, pott(ks+1,i,j-2) )
1946  call check( __line__, pott(ks+1,i,j-1) )
1947  call check( __line__, pott(ks+1,i,j ) )
1948  call check( __line__, pott(ks+1,i,j+1) )
1949  call check( __line__, pott(ks+1,i,j+2) )
1950  call check( __line__, g(ks ,i,j,i_xyw) )
1951  call check( __line__, g(ks+2,i,j,i_xyw) )
1952  call check( __line__, g(ks+1,i,j,i_xyz) )
1953  call check( __line__, rcs2t(ks ,i,j) )
1954  call check( __line__, rcs2t(ks+1 ,i,j) )
1955  call check( __line__, rcs2t(ks+2,i,j) )
1956 #endif
1957  ! k,i,j
1958  m(1,ks+1,i,j) = &
1959  - ( ( fact_n * (pott(ks+2,i,j)+pott(ks+1,i,j)) &
1960  + fact_f * (pott(ks+3,i,j)+pott(ks ,i,j)) ) * rfdz(ks+1) / g(ks+1,i,j,i_xyw) &
1961  + 0.5_rp * (pott(ks+1,i,j)+pott(ks ,i,j)) * rfdz(ks ) / g(ks ,i,j,i_xyw) &
1962  ) * j33g * j33g * rfdz(ks+1) &
1963  - ( ( fact_n * (pott(ks+1,i+1,j)+pott(ks+1,i ,j)) &
1964  + fact_f * (pott(ks+1,i+2,j)+pott(ks+1,i-1,j)) ) * rfdx(i ) &
1965  + ( fact_n * (pott(ks+1,i ,j)+pott(ks+1,i-1,j)) &
1966  + fact_f * (pott(ks+1,i+1,j)+pott(ks+1,i-2,j)) ) * rfdx(i-1) &
1967  ) * g(ks+1,i,j,i_xyz) * rfdx(i) &
1968  - ( ( fact_n * (pott(ks+1,i,j+1)+pott(ks+1,i,j )) &
1969  + fact_f * (pott(ks+1,i,j+2)+pott(ks+1,i,j-1)) ) * rfdy(j ) &
1970  + ( fact_n * (pott(ks+1,i,j )+pott(ks+1,i,j-1)) &
1971  + fact_f * (pott(ks+1,i,j-1)+pott(ks+1,i,j-2)) ) * rfdy(j-1) &
1972  ) * g(ks+1,i,j,i_xyz) * rfdy(j) &
1973  - g(ks+1,i,j,i_xyz) * rcs2t(ks+1,i,j) * rdt * rdt
1974  ! k-1
1975  m(2,ks+1,i,j) = j33g * j33g / g(ks,i,j,i_xyw) &
1976  * 0.5_rp * (pott(ks+1,i,j)+pott(ks ,i,j)) &
1977  * rfdz(ks ) * rcdz(ks+1) &
1978  - grav * j33g * rcs2t(ks ,i,j) / ( fdz(ks+1)+fdz(ks) )
1979  ! k+1
1980  m(3,ks+1,i,j) = j33g * j33g / g(ks+2,i,j,i_xyw) &
1981  * ( fact_n * (pott(ks+2,i,j)+pott(ks+1,i,j)) &
1982  + fact_f * (pott(ks+3,i,j)+pott(ks ,i,j)) ) &
1983  * rfdz(ks+1) * rcdz(ks+1) &
1984  + grav * j33g * rcs2t(ks+2,i,j) / ( fdz(ks+1)+fdz(ks) )
1985 #ifdef DEBUG
1986  call check( __line__, pott(ke-3,i,j) )
1987  call check( __line__, pott(ke-2,i,j) )
1988  call check( __line__, pott(ke-1,i,j) )
1989  call check( __line__, pott(ke ,i,j) )
1990  call check( __line__, pott(ke-1,i-2,j) )
1991  call check( __line__, pott(ke-1,i-1,j) )
1992  call check( __line__, pott(ke-1,i ,j) )
1993  call check( __line__, pott(ke-1,i+1,j) )
1994  call check( __line__, pott(ke-1,i+2,j) )
1995  call check( __line__, pott(ke-1,i,j-2) )
1996  call check( __line__, pott(ke-1,i,j-1) )
1997  call check( __line__, pott(ke-1,i,j ) )
1998  call check( __line__, pott(ke-1,i,j+1) )
1999  call check( __line__, pott(ke-1,i,j+2) )
2000  call check( __line__, g(ke-2,i,j,i_xyw) )
2001  call check( __line__, g(ke ,i,j,i_xyw) )
2002  call check( __line__, g(ke-1,i,j,i_xyz) )
2003  call check( __line__, rcs2t(ke-2,i,j) )
2004  call check( __line__, rcs2t(ke-1,i,j) )
2005  call check( __line__, rcs2t(ke ,i,j) )
2006 #endif
2007  ! k,i,j
2008  m(1,ke-1,i,j) = &
2009  - ( 0.5_rp * (pott(ke ,i,j)+pott(ke-1,i,j)) * rfdz(ke-1) / g(ke-1,i,j,i_xyw) &
2010  + ( fact_n * (pott(ke-1,i,j)+pott(ke-2,i,j)) &
2011  + fact_f * (pott(ke ,i,j)+pott(ke-3,i,j)) ) * rfdz(ke-2) / g(ke-2,i,j,i_xyw) &
2012  ) * j33g * j33g * rfdz(ke-1) &
2013  - ( ( fact_n * (pott(ke-1,i+1,j)+pott(ke-1,i ,j)) &
2014  + fact_f * (pott(ke-1,i+2,j)+pott(ke-1,i-1,j)) ) * rfdx(i ) &
2015  + ( fact_n * (pott(ke-1,i ,j)+pott(ke-1,i-1,j)) &
2016  + fact_f * (pott(ke-1,i+1,j)+pott(ke-1,i-2,j)) ) * rfdx(i-1) &
2017  ) * g(ke-1,i,j,i_xyz) * rfdx(i) &
2018  - ( ( fact_n * (pott(ke-1,i,j+1)+pott(ke-1,i,j )) &
2019  + fact_f * (pott(ke-1,i,j+2)+pott(ke-1,i,j-1)) ) * rfdy(j ) &
2020  + ( fact_n * (pott(ke-1,i,j )+pott(ke-1,i,j-1)) &
2021  + fact_f * (pott(ke-1,i,j-1)+pott(ke-1,i,j-2)) ) * rfdy(j-1) &
2022  ) * g(ke-1,i,j,i_xyz) * rfdy(j) &
2023  - g(ke-1,i,j,i_xyz) * rcs2t(ke-1,i,j) * rdt * rdt
2024  ! k-1
2025  m(2,ke-1,i,j) = j33g * j33g / g(ke-2,i,j,i_xyw) &
2026  * ( fact_n * (pott(ke-1,i,j)+pott(ke-2,i,j)) &
2027  + fact_f * (pott(ke ,i,j)+pott(ke-3,i,j)) ) &
2028  * rfdz(ke-2) * rcdz(ke-1) &
2029  - grav * j33g * rcs2t(ke-2,i,j) / ( fdz(ke-1)+fdz(ke-2) )
2030  ! k+1
2031  m(3,ke-1,i,j) = j33g * j33g / g(ke ,i,j,i_xyw) &
2032  * 0.5_rp * (pott(ke ,i,j)+pott(ke-1,i,j)) &
2033  * rfdz(ke-1) * rcdz(ke-1) &
2034  + grav * j33g * rcs2t(ke ,i,j) / ( fdz(ke-1)+fdz(ke-2) )
2035 #ifdef DEBUG
2036  call check( __line__, pott(ke-1,i,j) )
2037  call check( __line__, pott(ke ,i,j) )
2038  call check( __line__, pott(ke,i-2,j) )
2039  call check( __line__, pott(ke,i-1,j) )
2040  call check( __line__, pott(ke,i ,j) )
2041  call check( __line__, pott(ke,i+1,j) )
2042  call check( __line__, pott(ke,i+2,j) )
2043  call check( __line__, pott(ke,i,j-2) )
2044  call check( __line__, pott(ke,i,j-1) )
2045  call check( __line__, pott(ke,i,j ) )
2046  call check( __line__, pott(ke,i,j+1) )
2047  call check( __line__, pott(ke,i,j+2) )
2048  call check( __line__, g(ke-1,i,j,i_xyw) )
2049  call check( __line__, g(ke,i,j,i_xyz) )
2050  call check( __line__, rcs2t(ke-1,i,j) )
2051  call check( __line__, rcs2t(ke,i,j) )
2052 #endif
2053  ! k,i,j
2054  m(1,ke,i,j) = &
2055  - ( &
2056  + 0.5_rp * (pott(ke ,i,j)+pott(ke-1,i,j)) * rfdz(ke-1) / g(ke-1,i,j,i_xyw) &
2057  ) * j33g * j33g * rfdz(ke) &
2058  - ( ( fact_n * (pott(ke,i+1,j)+pott(ke,i ,j)) &
2059  + fact_f * (pott(ke,i+2,j)+pott(ke,i-1,j)) ) * rfdx(i ) &
2060  + ( fact_n * (pott(ke,i ,j)+pott(ke,i-1,j)) &
2061  + fact_f * (pott(ke,i+1,j)+pott(ke,i-2,j)) ) * rfdx(i-1) &
2062  ) * g(ke,i,j,i_xyz) * rfdx(i) &
2063  - ( ( fact_n * (pott(ke,i,j+1)+pott(ke,i,j )) &
2064  + fact_f * (pott(ke,i,j+2)+pott(ke,i,j-1)) ) * rfdy(j ) &
2065  + ( fact_n * (pott(ke,i,j )+pott(ke,i,j-1)) &
2066  + fact_f * (pott(ke,i,j-1)+pott(ke,i,j-2)) ) * rfdy(j-1) &
2067  ) * g(ke,i,j,i_xyz) * rfdy(j) &
2068  - g(ke,i,j,i_xyz) * rcs2t(ke,i,j) * rdt * rdt &
2069  - grav * j33g * 0.5_rp * rcs2t(ke,i,j) * rcdz(ke)
2070  ! k-1
2071  m(2,ke,i,j) = j33g * j33g / g(ke-1,i,j,i_xyw) &
2072  * 0.5_rp * (pott(ke ,i,j)+pott(ke-1,i,j)) &
2073  * rfdz(ke-1) * rcdz(ke) &
2074  - grav * j33g * 0.5_rp * rcs2t(ke,i,j) * rcdz(ke)
2075  enddo
2076  enddo
2077 #ifdef DEBUG
2078  k = iundef; i = iundef; j = iundef
2079 #endif
2080 
2081  do j = jjs, jje
2082  do i = iis, iie
2083  do k = ks, ke
2084 #ifdef DEBUG
2085  call check( __line__, g(k,i-1,j,i_xyz) )
2086  call check( __line__, g(k,i+1,j,i_xyz) )
2087  call check( __line__, g(k,i,j-1,i_xyz) )
2088  call check( __line__, g(k,i,j+1,i_xyz) )
2089  call check( __line__, pott(k,i-2,j ) )
2090  call check( __line__, pott(k,i-1,j ) )
2091  call check( __line__, pott(k,i ,j ) )
2092  call check( __line__, pott(k,i+1,j ) )
2093  call check( __line__, pott(k,i+2,j ) )
2094  call check( __line__, pott(k,i ,j-2) )
2095  call check( __line__, pott(k,i ,j-1) )
2096  call check( __line__, pott(k,i ,j ) )
2097  call check( __line__, pott(k,i ,j+1) )
2098  call check( __line__, pott(k,i ,j+2) )
2099 #endif
2100  ! i-1
2101  m(4,k,i,j) = g(k,i-1,j,i_xyz) &
2102  * ( fact_n * (pott(k,i ,j)+pott(k,i-1,j)) &
2103  + fact_f * (pott(k,i+1,j)+pott(k,i-2,j)) ) &
2104  * rfdx(i-1) * rcdx(i)
2105  ! i+1
2106  m(5,k,i,j) = g(k,i+1,j,i_xyz) &
2107  * ( fact_n * (pott(k,i+1,j)+pott(k,i ,j)) &
2108  + fact_f * (pott(k,i+2,j)+pott(k,i-1,j)) ) &
2109  * rfdx(i ) * rcdx(i)
2110  ! j-1
2111  m(6,k,i,j) = g(k,i,j-1,i_xyz) &
2112  * ( fact_n * (pott(k,i,j )+pott(k,i,j-1)) &
2113  + fact_f * (pott(k,i,j+1)+pott(k,i,j-2)) ) &
2114  * rfdy(j-1) * rcdy(j)
2115  ! j+1
2116  m(7,k,i,j) = g(k,i,j+1,i_xyz) &
2117  * ( fact_n * (pott(k,i,j+1)+pott(k,i,j )) &
2118  + fact_f * (pott(k,i,j+2)+pott(k,i,j-1)) ) &
2119  * rfdy(j ) * rcdy(j)
2120  enddo
2121  enddo
2122  enddo
2123 #ifdef DEBUG
2124  k = iundef; i = iundef; j = iundef
2125 #endif
2126 
2127  return
scale_grid_index::ke
integer, public ke
end point of inner domain: z, local
Definition: scale_grid_index.F90:59
scale_gridtrans::i_xyw
integer, public i_xyw
Definition: scale_gridtrans.F90:50
scale_grid_index::ks
integer, public ks
start point of inner domain: z, local
Definition: scale_grid_index.F90:58
scale_gridtrans::i_xyz
integer, public i_xyz
Definition: scale_gridtrans.F90:49
Here is the call graph for this function:
Here is the caller graph for this function:
Generated by   doxygen 1.8.14