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_prc::prc_abort().

Referenced by scale_atmos_dyn_tstep_short::atmos_dyn_tstep_short_regist().

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:338
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_io::io_fid_conf, scale_prc::prc_abort(), scale_precision::rp, and scale_precision::sp.

Referenced by scale_atmos_dyn_tstep_short::atmos_dyn_tstep_short_regist().

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
scale_io::io_fid_conf
integer, public io_fid_conf
Config file ID.
Definition: scale_io.F90:55
scale_prc
module PROCESS
Definition: scale_prc.F90:11
scale_prc::prc_abort
subroutine, public prc_abort
Abort Process.
Definition: scale_prc.F90:338
scale_precision::sp
integer, parameter, public sp
Definition: scale_precision.F90:31
scale_precision::rp
integer, parameter, public rp
Definition: scale_precision.F90:41
scale_precision::dp
integer, parameter, public dp
Definition: scale_precision.F90:32
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_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::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::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().

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

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().

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

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().

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