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, divdmp_coef, DDIV, FLAG_FCT_MOMENTUM, FLAG_FCT_T, FLAG_FCT_ALONG_STREAM, CDZ, FDZ, FDX, FDY, RCDZ, RCDX, RCDY, RFDZ, RFDX, RFDY, PHI, GSQRT, J13G, J23G, J33G, MAPF, REF_dens, REF_rhot, BND_W, BND_E, BND_S, BND_N, dtrk, dt)
 
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

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( io_l ) write(io_fid_log,*) '*** HIVI Register'
91 
92  if ( atmos_dyn_type /= 'FVM-HIVI' .AND. atmos_dyn_type /= 'HIVI' ) then
93  write(*,*) 'xxx ATMOS_DYN_TYPE is not FVM-HIVI. Check!'
94  call prc_mpistop
95  endif
96 
97  va_out = va_fvm_hivi
98  var_name(:) = ""
99  var_desc(:) = ""
100  var_unit(:) = ""
101 
102  return
scale_process::prc_mpistop
subroutine, public prc_mpistop
Abort MPI.
Definition: scale_process.F90:234
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 108 of file scale_atmos_dyn_tstep_short_fvm_hivi.F90.

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

Referenced by scale_atmos_dyn_tstep_short::atmos_dyn_tstep_short_regist().

108  use scale_process, only: &
109  prc_mpistop
110  implicit none
111 
112  namelist / param_atmos_dyn_tstep_fvm_hivi / &
113  itmax, &
114  epsilon
115 
116  integer :: ierr
117  !---------------------------------------------------------------------------
118 
119  if( io_l ) write(io_fid_log,*) '*** HIVI Setup'
120 #ifdef HIVI_BICGSTAB
121  if ( io_l ) write(io_fid_log,*) '*** USING Bi-CGSTAB'
122 #else
123  if ( io_l ) write(io_fid_log,*) '*** USING Multi-Grid'
124  write(*,*) 'xxx Not Implemented yet'
125  call prc_mpistop
126 #endif
127 
128  ! currently, vertical difference scheme for potential temperature is the CD4
129  itmax = 100
130  epsilon = 0.1_rp ** (rp*2)
131 
132  !--- read namelist
133  rewind(io_fid_conf)
134  read(io_fid_conf,nml=param_atmos_dyn_tstep_fvm_hivi,iostat=ierr)
135 
136  if( ierr < 0 ) then !--- missing
137  if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
138  elseif( ierr > 0 ) then !--- fatal error
139  write(*,*) 'xxx Not appropriate names in namelist PARAM_ATMOS_DYN_TSTEP_FVM_HIVI. Check!'
140  call prc_mpistop
141  endif
142  if( io_l ) write(io_fid_log,nml=param_atmos_dyn_tstep_fvm_hivi)
143 
144  if ( rp == dp ) then
145  mtype = mpi_double_precision
146  elseif( rp == sp ) then
147  mtype = mpi_real
148  else
149  write(*,*) 'xxx Unsupported precision'
150  call prc_mpistop
151  endif
152 
153  return
scale_process::prc_mpistop
subroutine, public prc_mpistop
Abort MPI.
Definition: scale_process.F90:234
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), intent(in)  divdmp_coef,
real(rp), dimension(ka,ia,ja), intent(in)  DDIV,
logical, intent(in)  FLAG_FCT_MOMENTUM,
logical, intent(in)  FLAG_FCT_T,
logical, intent(in)  FLAG_FCT_ALONG_STREAM,
real(rp), dimension(ka), intent(in)  CDZ,
real(rp), dimension(ka-1), intent(in)  FDZ,
real(rp), dimension(ia-1), intent(in)  FDX,
real(rp), dimension(ja-1), intent(in)  FDY,
real(rp), dimension(ka), intent(in)  RCDZ,
real(rp), dimension(ia), intent(in)  RCDX,
real(rp), dimension(ja), intent(in)  RCDY,
real(rp), dimension(ka-1), intent(in)  RFDZ,
real(rp), dimension(ia-1), intent(in)  RFDX,
real(rp), dimension(ja-1), intent(in)  RFDY,
real(rp), dimension (ka,ia,ja), intent(in)  PHI,
real(rp), dimension (ka,ia,ja,7), intent(in)  GSQRT,
real(rp), dimension (ka,ia,ja,7), intent(in)  J13G,
real(rp), dimension (ka,ia,ja,7), intent(in)  J23G,
real(rp), intent(in)  J33G,
real(rp), dimension (ia,ja,2,4), intent(in)  MAPF,
real(rp), dimension(ka,ia,ja), intent(in)  REF_dens,
real(rp), dimension(ka,ia,ja), intent(in)  REF_rhot,
logical, intent(in)  BND_W,
logical, intent(in)  BND_E,
logical, intent(in)  BND_S,
logical, intent(in)  BND_N,
real(rp), intent(in)  dtrk,
real(rp), intent(in)  dt 
)
Parameters
[in]phigeopotential
[in]gsqrtvertical metrics {G}^1/2
[in]j13g(1,3) element of Jacobian matrix
[in]j23g(2,3) element of Jacobian matrix
[in]j33g(3,3) element of Jacobian matrix
[in]mapfmap factor
[in]ref_densreference density

Definition at line 175 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::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().

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

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

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