SCALE-RM
scale_atmos_dyn_tstep_short_fvm_hivi.F90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
13 !-------------------------------------------------------------------------------
14 #define HIVI_BICGSTAB 1
15 #include "scalelib.h"
16 module scale_atmos_dyn_tstep_short_fvm_hivi
17  !-----------------------------------------------------------------------------
18  !
19  !++ used modules
20  !
21  use mpi
22  use scale_precision
23  use scale_io
24  use scale_prof
25  use scale_atmos_grid_cartesc_index
26  use scale_index
27  use scale_tracer
28 #if defined DEBUG || defined QUICKDEBUG
29  use scale_debug, only: &
30  check
31  use scale_const, only: &
32  undef => const_undef, &
33  iundef => const_undef2
34 #endif
35  !-----------------------------------------------------------------------------
36  implicit none
37  private
38  !-----------------------------------------------------------------------------
39  !
40  !++ Public procedure
41  !
42  public :: atmos_dyn_tstep_short_fvm_hivi_regist
43  public :: atmos_dyn_tstep_short_fvm_hivi_setup
44  public :: atmos_dyn_tstep_short_fvm_hivi
45 
46  !-----------------------------------------------------------------------------
47  !
48  !++ Public parameters & variables
49  !
50  !-----------------------------------------------------------------------------
51  !
52  !++ Private procedure
53  !
54  !-----------------------------------------------------------------------------
55  !
56  !++ Private parameters & variables
57  !
58  integer, private, parameter :: VA_FVM_HIVI = 0
59 
60  integer, private :: ITMAX
61  real(RP), private :: epsilon
62 
63  integer, private :: mtype ! MPI DATATYPE
64 
65  ! tentative
66  real(RP), private, parameter :: FACT_N = 7.0_rp / 12.0_rp
67  real(RP), private, parameter :: FACT_F = -1.0_rp / 12.0_rp
68 
69  !-----------------------------------------------------------------------------
70 contains
71  !-----------------------------------------------------------------------------
73  subroutine atmos_dyn_tstep_short_fvm_hivi_regist( &
74  ATMOS_DYN_TYPE, &
75  VA_out, &
76  VAR_NAME, &
77  VAR_DESC, &
78  VAR_UNIT )
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
107  end subroutine atmos_dyn_tstep_short_fvm_hivi_regist
108 
109  !-----------------------------------------------------------------------------
111  subroutine atmos_dyn_tstep_short_fvm_hivi_setup
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
158  end subroutine atmos_dyn_tstep_short_fvm_hivi_setup
159 
160  !-----------------------------------------------------------------------------
161  subroutine atmos_dyn_tstep_short_fvm_hivi( &
162  DENS_RK, MOMZ_RK, MOMX_RK, MOMY_RK, RHOT_RK, &
163  PROG_RK, &
164  mflx_hi, tflx_hi, &
165  DENS0, MOMZ0, MOMX0, MOMY0, RHOT0, &
166  DENS, MOMZ, MOMX, MOMY, RHOT, &
167  DENS_t, MOMZ_t, MOMX_t, MOMY_t, RHOT_t, &
168  PROG0, PROG, &
169  DPRES0, RT2P, CORIOLI, &
170  num_diff, wdamp_coef, divdmp_coef, DDIV, &
171  FLAG_FCT_MOMENTUM, FLAG_FCT_T, &
172  FLAG_FCT_ALONG_STREAM, &
173  CDZ, FDZ, FDX, FDY, &
174  RCDZ, RCDX, RCDY, RFDZ, RFDX, RFDY, &
175  PHI, GSQRT, J13G, J23G, J33G, MAPF, &
176  REF_dens, REF_rhot, &
177  BND_W, BND_E, BND_S, BND_N, TwoD, &
178  dtrk, last )
179  use scale_atmos_grid_cartesc_index
180  use scale_const, only: &
181  grav => const_grav, &
182  p00 => const_pre00
183  use scale_comm_cartesc, only: &
184  comm_vars8, &
185  comm_wait
186  use scale_atmos_dyn_fvm_fct, only: &
187  atmos_dyn_fvm_fct
188  use scale_atmos_dyn_fvm_flux, only: &
189  atmos_dyn_fvm_flux_valuew_z, &
190  atmos_dyn_fvm_fluxz_xyz, &
191  atmos_dyn_fvm_fluxx_xyz, &
192  atmos_dyn_fvm_fluxy_xyz, &
193  atmos_dyn_fvm_fluxz_xyw, &
194  atmos_dyn_fvm_fluxj13_xyw, &
195  atmos_dyn_fvm_fluxj23_xyw, &
196  atmos_dyn_fvm_fluxx_xyw, &
197  atmos_dyn_fvm_fluxy_xyw, &
198  atmos_dyn_fvm_fluxz_uyz, &
199  atmos_dyn_fvm_fluxj13_uyz, &
200  atmos_dyn_fvm_fluxj23_uyz, &
201  atmos_dyn_fvm_fluxx_uyz, &
202  atmos_dyn_fvm_fluxy_uyz, &
203  atmos_dyn_fvm_fluxz_xvz, &
204  atmos_dyn_fvm_fluxj13_xvz, &
205  atmos_dyn_fvm_fluxj23_xvz, &
206  atmos_dyn_fvm_fluxx_xvz, &
207  atmos_dyn_fvm_fluxy_xvz
208  implicit none
209 
210  real(rp), intent(out) :: dens_rk(ka,ia,ja) ! prognostic variables
211  real(rp), intent(out) :: momz_rk(ka,ia,ja) !
212  real(rp), intent(out) :: momx_rk(ka,ia,ja) !
213  real(rp), intent(out) :: momy_rk(ka,ia,ja) !
214  real(rp), intent(out) :: rhot_rk(ka,ia,ja) !
215 
216  real(rp), intent(out) :: prog_rk(ka,ia,ja,va) !
217 
218  real(rp), intent(inout) :: mflx_hi(ka,ia,ja,3) ! rho * vel(x,y,z)
219  real(rp), intent(out) :: tflx_hi(ka,ia,ja,3) ! rho * theta * vel(x,y,z)
220 
221  real(rp), intent(in),target :: dens0(ka,ia,ja) ! prognostic variables
222  real(rp), intent(in),target :: momz0(ka,ia,ja) ! at previous dynamical time step
223  real(rp), intent(in),target :: momx0(ka,ia,ja) !
224  real(rp), intent(in),target :: momy0(ka,ia,ja) !
225  real(rp), intent(in),target :: rhot0(ka,ia,ja) !
226 
227  real(rp), intent(in) :: dens(ka,ia,ja) ! prognostic variables
228  real(rp), intent(in) :: momz(ka,ia,ja) ! at previous RK step
229  real(rp), intent(in) :: momx(ka,ia,ja) !
230  real(rp), intent(in) :: momy(ka,ia,ja) !
231  real(rp), intent(in) :: rhot(ka,ia,ja) !
232 
233  real(rp), intent(in) :: dens_t(ka,ia,ja)
234  real(rp), intent(in) :: momz_t(ka,ia,ja)
235  real(rp), intent(in) :: momx_t(ka,ia,ja)
236  real(rp), intent(in) :: momy_t(ka,ia,ja)
237  real(rp), intent(in) :: rhot_t(ka,ia,ja)
238 
239  real(rp), intent(in) :: prog0(ka,ia,ja,va)
240  real(rp), intent(in) :: prog (ka,ia,ja,va)
241 
242  real(rp), intent(in) :: dpres0(ka,ia,ja)
243  real(rp), intent(in) :: rt2p(ka,ia,ja)
244  real(rp), intent(in) :: corioli(ia,ja)
245  real(rp), intent(in) :: num_diff(ka,ia,ja,5,3)
246  real(rp), intent(in) :: wdamp_coef(ka)
247  real(rp), intent(in) :: divdmp_coef
248  real(rp), intent(in) :: ddiv(ka,ia,ja)
249 
250  logical, intent(in) :: flag_fct_momentum
251  logical, intent(in) :: flag_fct_t
252  logical, intent(in) :: flag_fct_along_stream
253 
254  real(rp), intent(in) :: cdz(ka)
255  real(rp), intent(in) :: fdz(ka-1)
256  real(rp), intent(in) :: fdx(ia-1)
257  real(rp), intent(in) :: fdy(ja-1)
258  real(rp), intent(in) :: rcdz(ka)
259  real(rp), intent(in) :: rcdx(ia)
260  real(rp), intent(in) :: rcdy(ja)
261  real(rp), intent(in) :: rfdz(ka-1)
262  real(rp), intent(in) :: rfdx(ia-1)
263  real(rp), intent(in) :: rfdy(ja-1)
264 
265  real(rp), intent(in) :: phi (ka,ia,ja)
266  real(rp), intent(in) :: gsqrt (ka,ia,ja,7)
267  real(rp), intent(in) :: j13g (ka,ia,ja,7)
268  real(rp), intent(in) :: j23g (ka,ia,ja,7)
269  real(rp), intent(in) :: j33g
270  real(rp), intent(in) :: mapf (ia,ja,2,4)
271  real(rp), intent(in) :: ref_dens(ka,ia,ja)
272  real(rp), intent(in) :: ref_rhot(ka,ia,ja)
273 
274  logical, intent(in) :: bnd_w
275  logical, intent(in) :: bnd_e
276  logical, intent(in) :: bnd_s
277  logical, intent(in) :: bnd_n
278  logical, intent(in) :: twod
279 
280  real(rp), intent(in) :: dtrk
281  logical, intent(in) :: last
282 
283 
284  ! diagnostic variables (work space)
285  real(rp) :: pott(ka,ia,ja) ! potential temperature [K]
286  real(rp) :: ddens(ka,ia,ja) ! density deviation from reference
287  real(rp) :: dpres(ka,ia,ja) ! pressure deviation from reference
288  real(rp) :: dpres_n(ka,ia,ja) ! pressure deviation at t=n+1 [Pa]
289 
290  real(rp) :: sr(ka,ia,ja)
291  real(rp) :: sw(ka,ia,ja)
292  real(rp) :: su(ka,ia,ja)
293  real(rp) :: sv(ka,ia,ja)
294  real(rp) :: st(ka,ia,ja)
295  real(rp) :: rcs2t(ka,ia,ja)
296  real(rp) :: b(ka,ia,ja)
297 
298  real(rp) :: duvw
299 
300  real(rp) :: qflx_hi (ka,ia,ja,3)
301  real(rp) :: qflx_j13(ka,ia,ja)
302  real(rp) :: qflx_j23(ka,ia,ja)
303  real(rp) :: mflx_hi2(ka,ia,ja,3)
304 
305  ! for implicit solver
306  real(rp) :: m(7,ka,ia,ja)
307  real(rp) :: r(ka,ia,ja)
308  real(rp) :: p(ka,ia,ja)
309 
310  real(rp) :: zero(ka,ia,ja)
311 
312  real(rp) :: rdt
313 
314  integer :: iis, iie, jjs, jje
315  integer :: k, i, j
316 
317  zero = 0.0_rp
318 
319 #ifdef DEBUG
320  dpres(:,:,:) = undef
321  pott(:,:,:) = undef
322 
323  qflx_j13(:,:,:) = undef
324  qflx_j23(:,:,:) = undef
325  mflx_hi2(:,:,:,:) = undef
326 
327  sr(:,:,:) = undef
328  sw(:,:,:) = undef
329  su(:,:,:) = undef
330  sv(:,:,:) = undef
331  st(:,:,:) = undef
332  rcs2t(:,:,:) = undef
333 
334  b(:,:,:) = undef
335  r(:,:,:) = undef
336  p(:,:,:) = undef
337 #endif
338 
339 #if defined DEBUG || defined QUICKDEBUG
340  dens_rk( 1:ks-1,:,:) = undef
341  dens_rk(ke+1:ka ,:,:) = undef
342  momz_rk( 1:ks-1,:,:) = undef
343  momz_rk(ke+1:ka ,:,:) = undef
344  momx_rk( 1:ks-1,:,:) = undef
345  momx_rk(ke+1:ka ,:,:) = undef
346  momy_rk( 1:ks-1,:,:) = undef
347  momy_rk(ke+1:ka ,:,:) = undef
348  rhot_rk( 1:ks-1,:,:) = undef
349  rhot_rk(ke+1:ka ,:,:) = undef
350  prog_rk( 1:ks-1,:,:,:) = undef
351  prog_rk(ke+1:ka ,:,:,:) = undef
352 #endif
353 
354  rdt = 1.0_rp / dtrk
355 
356  do jjs = js, je, jblock
357  jje = jjs+jblock-1
358  do iis = is, ie, iblock
359  iie = iis+iblock-1
360 
361  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
362  do j = jjs-1, jje+1
363  do i = max(iis-1,1), min(iie+1,ia)
364  do k = ks, ke
365 #ifdef DEBUG
366  call check( __line__, dpres0(k,i,j) )
367  call check( __line__, rt2p(k,i,j) )
368  call check( __line__, rhot(k,i,j) )
369  call check( __line__, ref_rhot(k,i,j) )
370 #endif
371  dpres(k,i,j) = dpres0(k,i,j) + rt2p(k,i,j) * ( rhot(k,i,j) - ref_rhot(k,i,j) )
372  enddo
373  dpres(ks-1,i,j) = dpres0(ks-1,i,j) - dens(ks,i,j) * ( phi(ks-1,i,j) - phi(ks+1,i,j) )
374  dpres(ke+1,i,j) = dpres0(ke+1,i,j) - dens(ke,i,j) * ( phi(ke+1,i,j) - phi(ke-1,i,j) )
375  enddo
376  enddo
377 
378  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
379  do j = jjs-1, jje+1
380  do i = iis-1, iie+1
381  do k = ks, ke
382 #ifdef DEBUG
383  call check( __line__, dens(k,i,j) )
384  call check( __line__, ref_dens(k,i,j) )
385 #endif
386  ddens(k,i,j) = dens(k,i,j) - ref_dens(k,i,j)
387  enddo
388  enddo
389  enddo
390 
391  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
392  do j = jjs-2, jje+2
393  do i = iis-2, iie+2
394  do k = ks, ke
395 #ifdef DEBUG
396  call check( __line__, rhot(k,i,j) )
397  call check( __line__, dens(k,i,j) )
398 #endif
399  pott(k,i,j) = rhot(k,i,j) / dens(k,i,j)
400  enddo
401  enddo
402  enddo
403 #ifdef DEBUG
404  k = iundef; i = iundef; j = iundef
405 #endif
406 
407  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
408  do j = jjs, jje
409  do i = iis, iie
410  do k = ks, ke
411 #ifdef DEBUG
412  call check( __line__, rt2p(k,i,j) )
413 #endif
414  rcs2t(k,i,j) = 1.0_rp / rt2p(k,i,j)
415  enddo
416  enddo
417  enddo
418 #ifdef DEBUG
419  k = iundef; i = iundef; j = iundef
420 #endif
421 
422  !##### continuity equation #####
423 
424  ! at (x, y, w)
425  if ( twod ) then
426  !$omp parallel do private(j,k) OMP_SCHEDULE_
427  do j = jjs, jje
428  do k = ks, ke-1
429 #ifdef DEBUG
430  call check( __line__, momx(k+1,is,j) )
431  call check( __line__, momx(k ,is,j) )
432  call check( __line__, momy(k+1,is,j) )
433  call check( __line__, momy(k+1,is,j-1) )
434  call check( __line__, momy(k ,is,j) )
435  call check( __line__, momy(k ,is,j-1) )
436  call check( __line__, num_diff(k,is,j,i_dens,zdir) )
437  call check( __line__, gsqrt(k,is,j,i_xyw) )
438 #endif
439  mflx_hi2(k,is,j,zdir) = j23g(k,is,j,i_xyw) * 0.25_rp * ( momy(k+1,is,j)+momy(k+1,is,j-1) &
440  + momy(k ,is,j)+momy(k ,is,j-1) ) &
441  + gsqrt(k,is,j,i_xyw) * num_diff(k,is,j,i_dens,zdir)
442  enddo
443  enddo
444  else
445  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
446  do j = jjs, jje
447  do i = iis-1, iie
448  do k = ks, ke-1
449 #ifdef DEBUG
450  call check( __line__, momx(k+1,i ,j) )
451  call check( __line__, momx(k+1,i-1,j) )
452  call check( __line__, momx(k ,i ,j) )
453  call check( __line__, momx(k ,i+1,j) )
454  call check( __line__, momy(k+1,i,j) )
455  call check( __line__, momy(k+1,i,j-1) )
456  call check( __line__, momy(k ,i,j) )
457  call check( __line__, momy(k ,i,j-1) )
458  call check( __line__, num_diff(k,i,j,i_dens,zdir) )
459  call check( __line__, gsqrt(k,i,j,i_xyw) )
460 #endif
461  mflx_hi2(k,i,j,zdir) = j13g(k,i,j,i_xyw) * 0.25_rp * ( momx(k+1,i,j)+momx(k+1,i-1,j) &
462  + momx(k ,i,j)+momx(k ,i-1,j) ) &
463  + j23g(k,i,j,i_xyw) * 0.25_rp * ( momy(k+1,i,j)+momy(k+1,i,j-1) &
464  + momy(k ,i,j)+momy(k ,i,j-1) ) &
465  + gsqrt(k,i,j,i_xyw) * num_diff(k,i,j,i_dens,zdir)
466  enddo
467  enddo
468  enddo
469  end if
470 #ifdef DEBUG
471  k = iundef; i = iundef; j = iundef
472 #endif
473 
474  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
475  do j = jjs, jje
476  do i = iis, iie
477  mflx_hi2(ks-1,i,j,zdir) = 0.0_rp
478  mflx_hi2(ke ,i,j,zdir) = 0.0_rp
479  enddo
480  enddo
481 
482  ! at (u, y, z)
483  if ( .not. twod ) then
484  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
485  do j = jjs , jje
486  do i = iis-1, iie
487  do k = ks, ke
488 #ifdef DEBUG
489  call check( __line__, gsqrt(k,i,j,i_uyz) )
490  call check( __line__, num_diff(k,i,j,i_dens,xdir) )
491 #endif
492  mflx_hi2(k,i,j,xdir) = gsqrt(k,i,j,i_uyz) * num_diff(k,i,j,i_dens,xdir)
493  enddo
494  enddo
495  enddo
496  end if
497 #ifdef DEBUG
498  k = iundef; i = iundef; j = iundef
499 #endif
500  ! at (x, v, z)
501  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
502  do j = jjs-1, jje
503  do i = iis , iie
504  do k = ks, ke
505 #ifdef DEBUG
506  call check( __line__, gsqrt(k,i,j,i_xvz) )
507  call check( __line__, num_diff(k,i,j,i_dens,ydir) )
508 #endif
509  mflx_hi2(k,i,j,ydir) = gsqrt(k,i,j,i_xvz) * num_diff(k,i,j,i_dens,ydir)
510  enddo
511  enddo
512  enddo
513 #ifdef DEBUG
514  k = iundef; i = iundef; j = iundef
515 #endif
516 
517  !##### momentum equation (z) #####
518 
519  ! at (x, y, z)
520  ! note than z-index is added by -1
521  call atmos_dyn_fvm_fluxz_xyw( qflx_hi(:,:,:,zdir), & ! (out)
522  momz, momz, dens, & ! (in)
523  gsqrt(:,:,:,i_xyz), j33g, & ! (in)
524  num_diff(:,:,:,i_momz,zdir), & ! (in)
525  cdz, fdz, dtrk, &
526  iis, iie, jjs, jje ) ! (in)
527  if ( .not. twod ) &
528  call atmos_dyn_fvm_fluxj13_xyw( qflx_j13, & ! (out)
529  momx, momz, dens, & ! (in)
530  gsqrt(:,:,:,i_xyz), j13g(:,:,:,i_xyz), mapf(:,:,:,i_xy), & ! (in)
531  cdz, twod, &
532  iis, iie, jjs, jje ) ! (in)
533  call atmos_dyn_fvm_fluxj23_xyw( qflx_j23, & ! (out)
534  momy, momz, dens, & ! (in)
535  gsqrt(:,:,:,i_xyz), j23g(:,:,:,i_xyz), mapf(:,:,:,i_xy), & ! (in)
536  cdz, twod, &
537  iis, iie, jjs, jje ) ! (in)
538 
539  if ( .not. twod ) &
540  call atmos_dyn_fvm_fluxx_xyw( qflx_hi(:,:,:,xdir), & ! (out)
541  momx, momz, dens, & ! (in)
542  gsqrt(:,:,:,i_uyw), mapf(:,:,:,i_uy), & ! (in)
543  num_diff(:,:,:,i_momz,xdir), & ! (in)
544  cdz, twod, & ! (in)
545  iis, iie, jjs, jje ) ! (in)
546 
547  ! at (x, v, w)
548  call atmos_dyn_fvm_fluxy_xyw( qflx_hi(:,:,:,ydir), & ! (out)
549  momy, momz, dens, & ! (in)
550  gsqrt(:,:,:,i_xvw), mapf(:,:,:,i_xv), & ! (in)
551  num_diff(:,:,:,i_momz,ydir), & ! (in)
552  cdz, twod, & ! (in)
553  iis, iie, jjs, jje ) ! (in)
554 
555  !--- update momentum(z)
556  if ( twod ) then
557  !$omp parallel do private(j,k) OMP_SCHEDULE_
558  do j = jjs, jje
559  do k = ks, ke-1
560 #ifdef DEBUG
561  call check( __line__, qflx_hi(k ,is,j ,zdir) )
562  call check( __line__, qflx_hi(k-1,is,j ,zdir) )
563  call check( __line__, qflx_j23(k ,is,j ) )
564  call check( __line__, qflx_j23(k-1,is,j ) )
565  call check( __line__, qflx_hi(k ,is,j ,ydir) )
566  call check( __line__, qflx_hi(k ,is,j-1,ydir) )
567  call check( __line__, ddiv(k ,is,j) )
568  call check( __line__, ddiv(k+1,is,j) )
569  call check( __line__, momz0(k,is,j) )
570  call check( __line__, momz_t(k,is,j) )
571 #endif
572  sw(k,is,j) = &
573  - ( ( qflx_hi(k,is,j,zdir) - qflx_hi(k-1,is,j ,zdir) &
574  + qflx_j23(k,is,j) - qflx_j23(k-1,is,j ) ) * rfdz(k) &
575  + ( qflx_hi(k,is,j,ydir) - qflx_hi(k ,is,j-1,ydir) ) * rcdy(j) &
576  ) / gsqrt(k,is,j,i_xyw) &
577  - wdamp_coef(k) * momz0(k,is,j) & ! Rayleigh damping
578  + divdmp_coef * rdt * ( ddiv(k+1,is,j)-ddiv(k,is,j) ) * fdz(k) & ! divergence damping
579  + momz_t(k,is,j)
580  enddo
581  enddo
582  else
583  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
584  do j = jjs, jje
585  do i = iis, iie
586  do k = ks, ke-1
587 #ifdef DEBUG
588  call check( __line__, qflx_hi(k ,i ,j ,zdir) )
589  call check( __line__, qflx_hi(k-1,i ,j ,zdir) )
590  call check( __line__, qflx_j13(k ,i ,j ) )
591  call check( __line__, qflx_j13(k-1,i ,j ) )
592  call check( __line__, qflx_j23(k ,i ,j ) )
593  call check( __line__, qflx_j23(k-1,i ,j ) )
594  call check( __line__, qflx_hi(k ,i ,j ,xdir) )
595  call check( __line__, qflx_hi(k ,i-1,j ,xdir) )
596  call check( __line__, qflx_hi(k ,i ,j ,ydir) )
597  call check( __line__, qflx_hi(k ,i ,j-1,ydir) )
598  call check( __line__, ddiv(k ,i,j) )
599  call check( __line__, ddiv(k+1,i,j) )
600  call check( __line__, momz0(k,i,j) )
601  call check( __line__, momz_t(k,i,j) )
602 #endif
603  sw(k,i,j) = &
604  - ( ( qflx_hi(k,i,j,zdir) - qflx_hi(k-1,i ,j ,zdir) &
605  + qflx_j13(k,i,j) - qflx_j13(k-1,i ,j ) &
606  + qflx_j23(k,i,j) - qflx_j23(k-1,i ,j ) ) * rfdz(k) &
607  + ( qflx_hi(k,i,j,xdir) - qflx_hi(k ,i-1,j ,xdir) ) * rcdx(i) &
608  + ( qflx_hi(k,i,j,ydir) - qflx_hi(k ,i ,j-1,ydir) ) * rcdy(j) &
609  ) / gsqrt(k,i,j,i_xyw) &
610  - wdamp_coef(k) * momz0(k,i,j) & ! Rayleigh damping
611  + divdmp_coef * rdt * ( ddiv(k+1,i,j)-ddiv(k,i,j) ) * fdz(k) & ! divergence damping
612  + momz_t(k,i,j)
613  enddo
614  enddo
615  enddo
616  end if
617 #ifdef DEBUG
618  k = iundef; i = iundef; j = iundef
619 #endif
620  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
621  do j = jjs, jje
622  do i = iis, iie
623  sw(ks-1,i,j) = 0.0_rp
624  sw(ke ,i,j) = 0.0_rp
625  enddo
626  enddo
627 #ifdef DEBUG
628  k = iundef; i = iundef; j = iundef
629 #endif
630 
631 
632  !##### momentum equation (x) #####
633 
634  ! at (u, y, w)
635  call atmos_dyn_fvm_fluxz_uyz( qflx_hi(:,:,:,zdir), & ! (out)
636  momz, momx, dens, & ! (in)
637  gsqrt(:,:,:,i_uyw), j33g, & ! (in)
638  num_diff(:,:,:,i_momx,zdir), & ! (in)
639  cdz, twod, & ! (in)
640  iis, iie, jjs, jje ) ! (in)
641  if ( .not. twod ) &
642  call atmos_dyn_fvm_fluxj13_uyz( qflx_j13, & ! (out)
643  momx, momx, dens, & ! (in)
644  gsqrt(:,:,:,i_uyz), j13g(:,:,:,i_uyw), mapf(:,:,:,i_uy), & ! (in)
645  cdz, twod, & ! (in)
646  iis, iie, jjs, jje ) ! (in)
647  call atmos_dyn_fvm_fluxj23_uyz( qflx_j23, & ! (out)
648  momy, momx, dens, & ! (in)
649  gsqrt(:,:,:,i_uyz), j23g(:,:,:,i_uyw), mapf(:,:,:,i_uy), & ! (in)
650  cdz, twod, & ! (in)
651  iis, iie, jjs, jje ) ! (in)
652 
653  ! at (x, y, z)
654  ! note that x-index is added by -1
655  if ( .not. twod ) &
656  call atmos_dyn_fvm_fluxx_uyz( qflx_hi(:,:,:,xdir), & ! (out)
657  momx, momx, dens, & ! (in)
658  gsqrt(:,:,:,i_xyz), mapf(:,:,:,i_xy), & ! (in)
659  num_diff(:,:,:,i_momx,xdir), & ! (in)
660  cdz, twod, & ! (in)
661  iis, iie, jjs, jje ) ! (in)
662 
663  ! at (u, v, z)
664  call atmos_dyn_fvm_fluxy_uyz( qflx_hi(:,:,:,ydir), & ! (out)
665  momy, momx, dens, & ! (in)
666  gsqrt(:,:,:,i_uvz), mapf(:,:,:,i_uv), & ! (in)
667  num_diff(:,:,:,i_momx,ydir), & ! (in)
668  cdz, twod, & ! (in)
669  iis, iie, jjs, jje ) ! (in)
670 
671  !--- update momentum(x)
672  if ( twod ) then
673  !$omp parallel do private(j,k) OMP_SCHEDULE_
674  do j = jjs, jje
675  do k = ks, ke
676 #ifdef DEBUG
677  call check( __line__, qflx_hi(k ,is,j ,zdir) )
678  call check( __line__, qflx_hi(k-1,is,j ,zdir) )
679  call check( __line__, qflx_hi(k ,is,j ,ydir) )
680  call check( __line__, qflx_hi(k ,is,j-1,ydir) )
681  call check( __line__, corioli(is,j) )
682  call check( __line__, momy(k,is,j ) )
683  call check( __line__, momy(k,is,j-1) )
684  call check( __line__, momx0(k,is,j) )
685 #endif
686  su(k,is,j) = ( &
687  - ( ( qflx_hi(k,is,j,zdir) - qflx_hi(k-1,is,j ,zdir) &
688  + qflx_j23(k,is,j) - qflx_j23(k-1,is,j) ) * rcdz(k) &
689  + ( qflx_hi(k,is,j,ydir) - qflx_hi(k ,is,j-1,ydir) ) * rcdy(j) ) &
690  ) / gsqrt(k,is,j,i_uyz) &
691  + 0.5_rp * corioli(is,j) * ( momy(k,is,j)+momy(k,is,j-1) ) & ! coriolis force
692  + momx_t(k,is,j)
693  enddo
694  enddo
695  else
696  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
697  do j = jjs, jje
698  do i = iis, iie
699  do k = ks, ke
700 #ifdef DEBUG
701  call check( __line__, qflx_hi(k ,i ,j ,zdir) )
702  call check( __line__, qflx_hi(k-1,i ,j ,zdir) )
703  call check( __line__, qflx_hi(k ,i ,j ,xdir) )
704  call check( __line__, qflx_hi(k ,i-1,j ,xdir) )
705  call check( __line__, qflx_hi(k ,i ,j ,ydir) )
706  call check( __line__, qflx_hi(k ,i ,j-1,ydir) )
707  call check( __line__, dpres(k+1,i+1,j) )
708  call check( __line__, dpres(k+1,i ,j) )
709  call check( __line__, dpres(k-1,i+1,j) )
710  call check( __line__, dpres(k-1,i ,j) )
711  call check( __line__, corioli(i ,j) )
712  call check( __line__, corioli(i+1,j) )
713  call check( __line__, momy(k,i ,j ) )
714  call check( __line__, momy(k,i+1,j ) )
715  call check( __line__, momy(k,i ,j-1) )
716  call check( __line__, momy(k,i+1,j-1) )
717  call check( __line__, ddiv(k,i+1,j) )
718  call check( __line__, ddiv(k,i ,j) )
719  call check( __line__, momx0(k,i,j) )
720 #endif
721  su(k,i,j) = ( &
722  - ( ( qflx_hi(k,i,j,zdir) - qflx_hi(k-1,i ,j ,zdir) &
723  + qflx_j13(k,i,j) - qflx_j13(k-1,i ,j) &
724  + qflx_j23(k,i,j) - qflx_j23(k-1,i ,j) ) * rcdz(k) &
725  + ( qflx_hi(k,i,j,xdir) - qflx_hi(k ,i-1,j ,xdir) ) * rfdx(i) &
726  + ( qflx_hi(k,i,j,ydir) - qflx_hi(k ,i ,j-1,ydir) ) * rcdy(j) ) &
727  - ( j13g(k+1,i,j,i_uyz) * ( dpres(k+1,i+1,j)+dpres(k+1,i,j) ) &
728  - j13g(k-1,i,j,i_uyz) * ( dpres(k-1,i+1,j)+dpres(k-1,i,j) ) ) &
729  * 0.5_rp / ( fdz(k+1)+fdz(k) ) &
730  ) / gsqrt(k,i,j,i_uyz) &
731  + 0.125_rp * ( corioli(i,j)+corioli(i+1,j) ) &
732  * ( momy(k,i,j)+momy(k,i+1,j)+momy(k,i,j-1)+momy(k,i+1,j-1) ) & ! coriolis force
733  + divdmp_coef * rdt * ( ddiv(k,i+1,j)-ddiv(k,i,j) ) * fdx(i) & ! divergence damping
734  + momx_t(k,i,j)
735  enddo
736  enddo
737  enddo
738  end if
739 #ifdef DEBUG
740  k = iundef; i = iundef; j = iundef
741 #endif
742 
743  !##### momentum equation (y) #####
744 
745  ! at (x, v, w)
746  call atmos_dyn_fvm_fluxz_xvz( qflx_hi(:,:,:,zdir), & ! (out)
747  momz, momy, dens, & ! (in)
748  gsqrt(:,:,:,i_xvw), j33g, & ! (in)
749  num_diff(:,:,:,i_momy,zdir), & ! (in)
750  cdz, twod, & ! (in)
751  iis, iie, jjs, jje ) ! (in)
752  if ( .not. twod ) &
753  call atmos_dyn_fvm_fluxj13_xvz( qflx_j13, & ! (out)
754  momx, momy, dens, & ! (in)
755  gsqrt(:,:,:,i_xvz), j13g(:,:,:,i_xvw), mapf(:,:,:,i_xv), & ! (in)
756  cdz, twod, & ! (in)
757  iis, iie, jjs, jje ) ! (in)
758  call atmos_dyn_fvm_fluxj23_xvz( qflx_j23, & ! (out)
759  momy, momy, dens, & ! (in)
760  gsqrt(:,:,:,i_xvz), j23g(:,:,:,i_xvw), mapf(:,:,:,i_xv), & ! (in)
761  cdz, twod, & ! (in)
762  iis, iie, jjs, jje ) ! (in)
763 
764  ! at (u, v, z)
765  if ( .not. twod ) &
766  call atmos_dyn_fvm_fluxx_xvz( qflx_hi(:,:,:,xdir), & ! (out)
767  momx, momy, dens, & ! (in)
768  gsqrt(:,:,:,i_uvz), mapf(:,:,:,i_uv), & ! (in)
769  num_diff(:,:,:,i_momy,xdir), & ! (in)
770  cdz, twod, & ! (in)
771  iis, iie, jjs, jje ) ! (in)
772 
773  ! at (x, y, z)
774  ! note that y-index is added by -1
775  call atmos_dyn_fvm_fluxy_xvz( qflx_hi(:,:,:,ydir), & ! (out)
776  momy, momy, dens, & ! (in)
777  gsqrt(:,:,:,i_xyz), mapf(:,:,:,i_xy), & ! (in)
778  num_diff(:,:,:,i_momy,ydir), & ! (in
779  cdz, twod, & ! (in)
780  iis, iie, jjs, jje ) ! (in)
781 
782  !--- update momentum(y)
783  if ( twod ) then
784  !$omp parallel do private(j,k) OMP_SCHEDULE_
785  do j = jjs, jje
786  do k = ks, ke
787 #ifdef DEBUG
788  call check( __line__, qflx_hi(k ,is,j ,zdir) )
789  call check( __line__, qflx_hi(k-1,is,j ,zdir) )
790  call check( __line__, qflx_hi(k ,is,j ,ydir) )
791  call check( __line__, qflx_hi(k ,is,j-1,ydir) )
792  call check( __line__, dpres(k+1,is,j ) )
793  call check( __line__, dpres(k+1,is,j+1) )
794  call check( __line__, dpres(k-1,is,j ) )
795  call check( __line__, dpres(k-1,is,j+1) )
796  call check( __line__, corioli(is,j ) )
797  call check( __line__, corioli(is,j+1) )
798  call check( __line__, momx(k,is,j ) )
799  call check( __line__, momx(k,is,j+1) )
800  call check( __line__, ddiv(k,is,j+1) )
801  call check( __line__, ddiv(k,is,j ) )
802  call check( __line__, momy_t(k,is,j) )
803 #endif
804  sv(k,is,j) = &
805  ( - ( ( qflx_hi(k,is,j,zdir) - qflx_hi(k-1,is,j ,zdir) &
806  + qflx_j23(k,is,j) - qflx_j23(k-1,is,j ) ) * rcdz(k) &
807  + ( qflx_hi(k,is,j,ydir) - qflx_hi(k ,is,j-1,ydir) ) * rfdy(j) ) &
808  - ( j23g(k+1,is,j,i_xvz) * ( dpres(k+1,is,j+1)+dpres(k+1,is,j) ) &
809  - j23g(k-1,is,j,i_xvz) * ( dpres(k-1,is,j+1)+dpres(k-1,is,j) ) ) &
810  * 0.5_rp / ( fdz(k+1)+fdz(k) ) & ! pressure gradient force
811  ) / gsqrt(k,is,j,i_xvz) &
812  - 0.25_rp * ( corioli(is,j+1)+corioli(is,j) ) &
813  * ( momx(k,is,j+1)+momx(k,is,j) ) & ! coriolis force
814  + divdmp_coef * rdt * ( ddiv(k,is,j+1)-ddiv(k,is,j) ) * fdy(j) & ! divergence damping
815  + momy_t(k,is,j)
816  enddo
817  enddo
818  else
819  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
820  do j = jjs, jje
821  do i = iis, iie
822  do k = ks, ke
823 #ifdef DEBUG
824  call check( __line__, qflx_hi(k ,i ,j ,zdir) )
825  call check( __line__, qflx_hi(k-1,i ,j ,zdir) )
826  call check( __line__, qflx_hi(k ,i ,j ,xdir) )
827  call check( __line__, qflx_hi(k ,i-1,j ,xdir) )
828  call check( __line__, qflx_hi(k ,i ,j ,ydir) )
829  call check( __line__, qflx_hi(k ,i ,j-1,ydir) )
830  call check( __line__, dpres(k+1,i,j ) )
831  call check( __line__, dpres(k+1,i,j+1) )
832  call check( __line__, dpres(k-1,i,j ) )
833  call check( __line__, dpres(k-1,i,j+1) )
834  call check( __line__, corioli(i,j ) )
835  call check( __line__, corioli(i,j+1) )
836  call check( __line__, momx(k,i ,j ) )
837  call check( __line__, momx(k,i ,j+1) )
838  call check( __line__, momx(k,i-1,j ) )
839  call check( __line__, momx(k,i-1,j+1) )
840  call check( __line__, ddiv(k,i,j+1) )
841  call check( __line__, ddiv(k,i,j ) )
842  call check( __line__, momy_t(k,i,j) )
843 #endif
844  sv(k,i,j) = &
845  ( - ( ( qflx_hi(k,i,j,zdir) - qflx_hi(k-1,i ,j ,zdir) &
846  + qflx_j13(k,i,j) - qflx_j13(k-1,i ,j ) &
847  + qflx_j23(k,i,j) - qflx_j23(k-1,i ,j ) ) * rcdz(k) &
848  + ( qflx_hi(k,i,j,xdir) - qflx_hi(k ,i-1,j ,xdir) ) * rcdx(i) &
849  + ( qflx_hi(k,i,j,ydir) - qflx_hi(k ,i ,j-1,ydir) ) * rfdy(j) ) &
850  - ( j23g(k+1,i,j,i_xvz) * ( dpres(k+1,i,j+1)+dpres(k+1,i,j) ) &
851  - j23g(k-1,i,j,i_xvz) * ( dpres(k-1,i,j+1)+dpres(k-1,i,j) ) ) &
852  * 0.5_rp / ( fdz(k+1)+fdz(k) ) & ! pressure gradient force
853  ) / gsqrt(k,i,j,i_xvz) &
854  - 0.125_rp * ( corioli(i ,j+1)+corioli(i ,j) ) &
855  * ( momx(k,i,j+1)+momx(k,i,j)+momx(k,i-1,j+1)+momx(k,i-1,j) ) & ! coriolis force
856  + divdmp_coef * rdt * ( ddiv(k,i,j+1)-ddiv(k,i,j) ) * fdy(j) & ! divergence damping
857  + momy_t(k,i,j)
858  enddo
859  enddo
860  enddo
861  end if
862 #ifdef DEBUG
863  k = iundef; i = iundef; j = iundef
864 #endif
865 
866  !##### Thermodynamic Equation #####
867 
868 
869  ! at (x, y, w)
870  call atmos_dyn_fvm_fluxz_xyz( tflx_hi(:,:,:,zdir), & ! (out)
871  mflx_hi2(:,:,:,zdir), pott, gsqrt(:,:,:,i_xyw), & ! (in)
872  num_diff(:,:,:,i_rhot,zdir), & ! (in)
873  cdz, & ! (in)
874  iis, iie, jjs, jje ) ! (in)
875 
876  ! at (u, y, z)
877  if ( .not. twod ) &
878  call atmos_dyn_fvm_fluxx_xyz( tflx_hi(:,:,:,xdir), & ! (out)
879  mflx_hi2(:,:,:,xdir), pott, gsqrt(:,:,:,i_uyz), & ! (in)
880  num_diff(:,:,:,i_rhot,xdir), & ! (in)
881  cdz, & ! (in)
882  iis, iie, jjs, jje ) ! (in)
883 
884  ! at (x, v, z)
885  call atmos_dyn_fvm_fluxy_xyz( tflx_hi(:,:,:,ydir), & ! (out)
886  mflx_hi2(:,:,:,ydir), pott, gsqrt(:,:,:,i_xvz), & ! (in)
887  num_diff(:,:,:,i_rhot,ydir), & ! (in)
888  cdz, & ! (in)
889  iis, iie, jjs, jje ) ! (in)
890 
891  if ( twod ) then
892  !$omp parallel do private(j,k) OMP_SCHEDULE_
893  do j = jjs, jje
894  do k = ks, ke
895 #ifdef DEBUG
896  call check( __line__, tflx_hi(k ,is,j ,zdir) )
897  call check( __line__, tflx_hi(k-1,is,j ,zdir) )
898  call check( __line__, tflx_hi(k ,is,j ,ydir) )
899  call check( __line__, tflx_hi(k ,is,j-1,ydir) )
900  call check( __line__, rhot_t(k,is,j) )
901 #endif
902  st(k,is,j) = &
903  - ( ( tflx_hi(k,is,j,zdir) - tflx_hi(k-1,is,j ,zdir) ) * rcdz(k) &
904  + ( tflx_hi(k,is,j,ydir) - tflx_hi(k ,is,j-1,ydir) ) * rcdy(j) &
905  ) / gsqrt(k,is,j,i_xyz) &
906  + rhot_t(k,is,j)
907  enddo
908  enddo
909  else
910  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
911  do j = jjs, jje
912  do i = iis, iie
913  do k = ks, ke
914 #ifdef DEBUG
915  call check( __line__, tflx_hi(k ,i ,j ,zdir) )
916  call check( __line__, tflx_hi(k-1,i ,j ,zdir) )
917  call check( __line__, tflx_hi(k ,i ,j ,xdir) )
918  call check( __line__, tflx_hi(k ,i-1,j ,xdir) )
919  call check( __line__, tflx_hi(k ,i ,j ,ydir) )
920  call check( __line__, tflx_hi(k ,i ,j-1,ydir) )
921  call check( __line__, rhot_t(k,i,j) )
922 #endif
923  st(k,i,j) = &
924  - ( ( tflx_hi(k,i,j,zdir) - tflx_hi(k-1,i ,j ,zdir) ) * rcdz(k) &
925  + ( tflx_hi(k,i,j,xdir) - tflx_hi(k ,i-1,j ,xdir) ) * rcdx(i) &
926  + ( tflx_hi(k,i,j,ydir) - tflx_hi(k ,i ,j-1,ydir) ) * rcdy(j) &
927  ) / gsqrt(k,i,j,i_xyz) &
928  + rhot_t(k,i,j)
929  enddo
930  enddo
931  enddo
932  end if
933 #ifdef DEBUG
934  k = iundef; i = iundef; j = iundef
935 #endif
936 
937  enddo
938  enddo
939 
940  ! implicit solver
941 #ifdef HIVI_BICGSTAB
942 
943 #ifdef DEBUG
944  m(:,:,:,:) = undef
945 #endif
946 
947  call comm_vars8( su, 1 )
948  call comm_vars8( sv, 2 )
949  call comm_wait ( su, 1 )
950  call comm_wait ( sv, 2 )
951 
952  do jjs = js, je, jblock
953  jje = jjs+jblock-1
954  do iis = is, ie, iblock
955  iie = iis+iblock-1
956 
957  if ( twod ) then
958 
959  i = is
960  ! r = b - M x0
961  do j = jjs, jje
962  do k = ks+2, ke-2
963 #ifdef DEBUG
964  call check( __line__, momz(k-1,i,j) )
965  call check( __line__, momz(k ,i,j) )
966  call check( __line__, momx(k,i ,j) )
967  call check( __line__, momy(k,i,j-1) )
968  call check( __line__, momy(k,i,j ) )
969  call check( __line__, pott(k-2,i,j) )
970  call check( __line__, pott(k-1,i,j) )
971  call check( __line__, pott(k ,i,j) )
972  call check( __line__, pott(k+1,i,j) )
973  call check( __line__, pott(k+2,i,j) )
974  call check( __line__, pott(k,i ,j) )
975  call check( __line__, pott(k,i,j-2) )
976  call check( __line__, pott(k,i,j-1) )
977  call check( __line__, pott(k,i,j ) )
978  call check( __line__, pott(k,i,j+1) )
979  call check( __line__, pott(k,i,j+2) )
980  call check( __line__, sw(k-1,i,j) )
981  call check( __line__, sw(k ,i,j) )
982  call check( __line__, su(k,i-1,j) )
983  call check( __line__, su(k,i ,j) )
984  call check( __line__, sv(k,i,j-1) )
985  call check( __line__, sv(k,i,j ) )
986  call check( __line__, gsqrt(k,i ,j,i_uyz) )
987  call check( __line__, gsqrt(k,i,j-1,i_xvz) )
988  call check( __line__, gsqrt(k,i,j ,i_xvz) )
989  call check( __line__, gsqrt(k,i,j,i_xyz) )
990  call check( __line__, st(k,i,j) )
991  call check( __line__, dpres(k-1,i,j) )
992  call check( __line__, dpres(k ,i,j) )
993  call check( __line__, dpres(k+1,i,j) )
994  call check( __line__, rt2p(k-1,i,j) )
995  call check( __line__, rt2p(k ,i,j) )
996  call check( __line__, rt2p(k+1,i,j) )
997  call check( __line__, rhot(k-1,i,j) )
998  call check( __line__, rhot(k+1,i,j) )
999  call check( __line__, ddens(k-1,i,j) )
1000  call check( __line__, ddens(k+1,i,j) )
1001 #endif
1002  b(k,i,j) = ( &
1003  ( j33g * ( momz(k ,i,j) + dtrk*sw(k ,i,j) ) &
1004  * ( fact_n*(pott(k+1,i,j)+pott(k ,i,j)) &
1005  + fact_f*(pott(k+2,i,j)+pott(k-1,i,j)) ) &
1006  - j33g * ( momz(k-1,i,j) + dtrk*sw(k-1,i,j) ) &
1007  * ( fact_n*(pott(k ,i,j)+pott(k-1,i,j)) &
1008  + fact_f*(pott(k+1,i,j)+pott(k-2,i,j)) ) ) * rcdz(k) &
1009  + ( gsqrt(k,i,j ,i_xvz) * ( momy(k,i,j ) + dtrk*sv(k,i,j ) ) &
1010  * ( fact_n*(pott(k,i,j+1)+pott(k,i,j )) &
1011  + fact_f*(pott(k,i,j+2)+pott(k,i,j-1)) ) &
1012  - gsqrt(k,i,j-1,i_xvz) * ( momy(k,i,j-1) + dtrk*sv(k,i,j-1) ) &
1013  * ( fact_n*(pott(k,i,j )+pott(k,i,j-1)) &
1014  + fact_f*(pott(k,i,j+1)+pott(k,i,j-2)) ) ) * rcdy(j) &
1015  + gsqrt(k,i,j,i_xyz) * ( st(k,i,j) - dpres(k,i,j) * rcs2t(k,i,j) * rdt ) &
1016  ) * rdt &
1017  + grav * j33g * ( ( dpres(k+1,i,j)*rcs2t(k+1,i,j) &
1018  - dpres(k-1,i,j)*rcs2t(k-1,i,j) ) &
1019  - ( rhot(k+1,i,j)*ddens(k+1,i,j)/dens(k+1,i,j) &
1020  - rhot(k-1,i,j)*ddens(k-1,i,j)/dens(k-1,i,j) ) &
1021  ) / ( fdz(k) + fdz(k-1) )
1022  enddo
1023  enddo
1024 #ifdef DEBUG
1025  k = iundef; i = iundef; j = iundef
1026 #endif
1027  do j = jjs, jje
1028  do i = iis, iie
1029 #ifdef DEBUG
1030  call check( __line__, momz(ks,i,j) )
1031  call check( __line__, momx(ks,i,j) )
1032  call check( __line__, momy(ks,i,j) )
1033  call check( __line__, sw(ks,i,j) )
1034  call check( __line__, su(ks,i-1,j) )
1035  call check( __line__, su(ks,i ,j) )
1036  call check( __line__, sv(ks,i,j-1) )
1037  call check( __line__, sv(ks,i,j ) )
1038  call check( __line__, st(ks,i,j) )
1039  call check( __line__, pott(ks ,i,j) )
1040  call check( __line__, pott(ks+1,i,j) )
1041  call check( __line__, pott(ks,i ,j) )
1042  call check( __line__, pott(ks,i,j-2) )
1043  call check( __line__, pott(ks,i,j-1) )
1044  call check( __line__, pott(ks,i,j ) )
1045  call check( __line__, pott(ks,i,j+1) )
1046  call check( __line__, pott(ks,i,j+2) )
1047  call check( __line__, gsqrt(ks,i,j,i_uyz) )
1048  call check( __line__, gsqrt(ks,i,j,i_xvz) )
1049  call check( __line__, gsqrt(ks,i,j,i_xyz) )
1050  call check( __line__, dpres(ks ,i,j) )
1051  call check( __line__, dpres(ks+1,i,j) )
1052  call check( __line__, rt2p(ks ,i,j) )
1053  call check( __line__, rt2p(ks+1,i,j) )
1054  call check( __line__, rhot(ks+1,i,j) )
1055  call check( __line__, ddens(ks+1,i,j) )
1056 #endif
1057  b(ks,i,j) = ( &
1058  ( j33g * ( momz(ks ,i,j) + dtrk*sw(ks ,i,j) ) &
1059  * 0.5_rp*(pott(ks+1,i,j)+pott(ks ,i,j)) ) * rcdz(ks) &
1060  + ( gsqrt(ks,i,j ,i_xvz) * ( momy(ks,i,j ) + dtrk*sv(ks,i,j ) ) &
1061  * ( fact_n*(pott(ks,i,j+1)+pott(ks,i,j )) &
1062  + fact_f*(pott(ks,i,j+2)+pott(ks,i,j-1)) ) &
1063  - gsqrt(ks,i,j-1,i_xvz) * ( momy(ks,i,j-1) + dtrk*sv(ks,i,j-1) ) &
1064  * ( fact_n*(pott(ks,i,j )+pott(ks,i,j-1)) &
1065  + fact_f*(pott(ks,i,j+1)+pott(ks,i,j-2)) ) ) * rcdy(j) &
1066  + gsqrt(ks,i,j,i_xyz) * ( st(ks,i,j) - dpres(ks,i,j) * rcs2t(ks,i,j) * rdt ) &
1067  ) * rdt &
1068  + grav * j33g * 0.5_rp * ( ( dpres(ks,i,j)+dpres(ks+1,i,j) ) * rcs2t(ks,i,j) &
1069  - ( ddens(ks,i,j)+ddens(ks+1,i,j) ) * pott(ks,i,j) ) * rcdz(ks)
1070 #ifdef DEBUG
1071  call check( __line__, momz(ks ,i,j) )
1072  call check( __line__, momz(ks+1,i,j) )
1073  call check( __line__, momx(ks+1,i,j) )
1074  call check( __line__, momy(ks+1,i,j-1) )
1075  call check( __line__, momy(ks+1,i,j ) )
1076  call check( __line__, pott(ks ,i,j) )
1077  call check( __line__, pott(ks+1 ,i,j) )
1078  call check( __line__, pott(ks+1+1,i,j) )
1079  call check( __line__, pott(ks+1+2,i,j) )
1080  call check( __line__, pott(ks+1,i,j) )
1081  call check( __line__, pott(ks+1,i,j-2) )
1082  call check( __line__, pott(ks+1,i,j-1) )
1083  call check( __line__, pott(ks+1,i,j ) )
1084  call check( __line__, pott(ks+1,i,j+1) )
1085  call check( __line__, pott(ks+1,i,j+2) )
1086  call check( __line__, sw(ks ,i,j) )
1087  call check( __line__, sw(ks+1,i,j) )
1088  call check( __line__, su(ks+1,i,j) )
1089  call check( __line__, sv(ks+1,i,j-1) )
1090  call check( __line__, sv(ks+1,i,j ) )
1091  call check( __line__, gsqrt(ks+1,i,j,i_uyz) )
1092  call check( __line__, gsqrt(ks+1,i,j-1,i_xvz) )
1093  call check( __line__, gsqrt(ks+1,i,j ,i_xvz) )
1094  call check( __line__, gsqrt(ks+1,i,j,i_xyz) )
1095  call check( __line__, st(ks+1,i,j) )
1096  call check( __line__, dpres(ks ,i,j) )
1097  call check( __line__, dpres(ks+1,i,j) )
1098  call check( __line__, dpres(ks+2,i,j) )
1099  call check( __line__, rt2p(ks ,i,j) )
1100  call check( __line__, rt2p(ks+1,i,j) )
1101  call check( __line__, rt2p(ks+2,i,j) )
1102  call check( __line__, rhot(ks ,i,j) )
1103  call check( __line__, rhot(ks+2,i,j) )
1104  call check( __line__, ddens(ks ,i,j) )
1105  call check( __line__, ddens(ks+2,i,j) )
1106 #endif
1107  b(ks+1,i,j) = ( &
1108  ( j33g * ( momz(ks+1,i,j) + dtrk*sw(ks+1,i,j) ) &
1109  * ( fact_n*(pott(ks+2,i,j)+pott(ks+1,i,j)) &
1110  + fact_f*(pott(ks+3,i,j)+pott(ks ,i,j)) ) &
1111  - j33g * ( momz(ks+1-1,i,j) + dtrk*sw(ks+1-1,i,j) ) &
1112  * ( 0.5_rp*(pott(ks+1,i,j)+pott(ks,i,j)) ) ) * rcdz(ks+1) &
1113  + ( gsqrt(ks+1,i,j ,i_xvz) * ( momy(ks+1,i,j ) + dtrk*sv(ks+1,i,j ) ) &
1114  * ( fact_n*(pott(ks+1,i,j+1)+pott(ks+1,i,j )) &
1115  + fact_f*(pott(ks+1,i,j+2)+pott(ks+1,i,j-1)) ) &
1116  - gsqrt(ks+1,i,j-1,i_xvz) * ( momy(ks+1,i,j-1) + dtrk*sv(ks+1,i,j-1) ) &
1117  * ( fact_n*(pott(ks+1,i,j )+pott(ks+1,i,j-1)) &
1118  + fact_f*(pott(ks+1,i,j+1)+pott(ks+1,i,j-2)) ) ) * rcdy(j) &
1119  + gsqrt(ks+1,i,j,i_xyz) * ( st(ks+1,i,j) - dpres(ks+1,i,j) * rcs2t(ks+1,i,j) * rdt ) &
1120  ) * rdt &
1121  + grav * j33g * ( ( dpres(ks+2,i,j)*rcs2t(ks+2,i,j) &
1122  - dpres(ks ,i,j)*rcs2t(ks ,i,j) ) &
1123  - ( rhot(ks+2,i,j)*ddens(ks+2,i,j)/dens(ks+2,i,j) &
1124  - rhot(ks ,i,j)*ddens(ks ,i,j)/dens(ks ,i,j) ) &
1125  ) / ( fdz(ks+1) + fdz(ks) )
1126 #ifdef DEBUG
1127  call check( __line__, momz(ke-2,i,j) )
1128  call check( __line__, momz(ke-1,i,j) )
1129  call check( __line__, momx(ke-1,i,j) )
1130  call check( __line__, momy(ke-1,i,j-1) )
1131  call check( __line__, momy(ke-1,i,j ) )
1132  call check( __line__, pott(ke-3,i,j) )
1133  call check( __line__, pott(ke-2,i,j) )
1134  call check( __line__, pott(ke-1,i,j) )
1135  call check( __line__, pott(ke ,i,j) )
1136  call check( __line__, pott(ke-1,i ,j) )
1137  call check( __line__, pott(ke-1,i,j-2) )
1138  call check( __line__, pott(ke-1,i,j-1) )
1139  call check( __line__, pott(ke-1,i,j ) )
1140  call check( __line__, pott(ke-1,i,j+1) )
1141  call check( __line__, pott(ke-1,i,j+2) )
1142  call check( __line__, sw(ke-2,i,j) )
1143  call check( __line__, sw(ke-1,i,j) )
1144  call check( __line__, su(ke-1,i,j) )
1145  call check( __line__, sv(ke-1,i,j-1) )
1146  call check( __line__, sv(ke-1,i,j ) )
1147  call check( __line__, gsqrt(ke-1,i,j,i_uyz) )
1148  call check( __line__, gsqrt(ke-1,i,j-1,i_xvz) )
1149  call check( __line__, gsqrt(ke-1,i,j ,i_xvz) )
1150  call check( __line__, gsqrt(ke-1,i,j,i_xyz) )
1151  call check( __line__, st(ke-1,i,j) )
1152  call check( __line__, dpres(ke-2,i,j) )
1153  call check( __line__, dpres(ke-1,i,j) )
1154  call check( __line__, dpres(ke ,i,j) )
1155  call check( __line__, rt2p(ke-2,i,j) )
1156  call check( __line__, rt2p(ke-1,i,j) )
1157  call check( __line__, rt2p(ke ,i,j) )
1158  call check( __line__, rhot(ke-2,i,j) )
1159  call check( __line__, rhot(ke,i,j) )
1160  call check( __line__, ddens(ke-2,i,j) )
1161  call check( __line__, ddens(ke,i,j) )
1162 #endif
1163  b(ke-1,i,j) = ( &
1164  ( j33g * ( momz(ke-1,i,j) + dtrk*sw(ke-1,i,j) ) &
1165  * ( 0.5_rp*(pott(ke ,i,j)+pott(ke-1,i,j)) ) &
1166  - j33g * ( momz(ke-2,i,j) + dtrk*sw(ke-2,i,j) ) &
1167  * ( fact_n*(pott(ke-1,i,j)+pott(ke-2,i,j)) &
1168  + fact_f*(pott(ke ,i,j)+pott(ke-3,i,j)) ) ) * rcdz(ke-1) &
1169  + ( gsqrt(ke-1,i,j ,i_xvz) * ( momy(ke-1,i,j ) + dtrk*sv(ke-1,i,j ) ) &
1170  * ( fact_n*(pott(ke-1,i,j+1)+pott(ke-1,i,j )) &
1171  + fact_f*(pott(ke-1,i,j+2)+pott(ke-1,i,j-1)) ) &
1172  - gsqrt(ke-1,i,j-1,i_xvz) * ( momy(ke-1,i,j-1) + dtrk*sv(ke-1,i,j-1) ) &
1173  * ( fact_n*(pott(ke-1,i,j )+pott(ke-1,i,j-1)) &
1174  + fact_f*(pott(ke-1,i,j+1)+pott(ke-1,i,j-2)) ) ) * rcdy(j) &
1175  + gsqrt(ke-1,i,j,i_xyz) * ( st(ke-1,i,j) - dpres(ke-1,i,j) * rcs2t(ke-1,i,j) * rdt ) &
1176  ) * rdt &
1177  + grav * j33g * ( ( dpres(ke ,i,j)*rcs2t(ke ,i,j) &
1178  - dpres(ke-2,i,j)*rcs2t(ke-2,i,j) ) &
1179  - ( rhot(ke ,i,j)*ddens(ke ,i,j)/dens(ke ,i,j) &
1180  - rhot(ke-2,i,j)*ddens(ke-2,i,j)/dens(ke-2,i,j) )&
1181  ) / ( fdz(ke-1) + fdz(ke-1-1) )
1182 #ifdef DEBUG
1183  call check( __line__, momz(ke-1,i,j) )
1184  call check( __line__, momx(ke,i,j) )
1185  call check( __line__, momy(ke,i,j-1) )
1186  call check( __line__, momy(ke,i,j ) )
1187  call check( __line__, sw(ke-1,i,j) )
1188  call check( __line__, su(ke,i,j) )
1189  call check( __line__, sv(ke,i,j-1) )
1190  call check( __line__, sv(ke,i,j ) )
1191  call check( __line__, pott(ke-1,i,j) )
1192  call check( __line__, pott(ke ,i,j) )
1193  call check( __line__, pott(ke,i,j) )
1194  call check( __line__, pott(ke,i,j-2) )
1195  call check( __line__, pott(ke,i,j-1) )
1196  call check( __line__, pott(ke,i,j ) )
1197  call check( __line__, pott(ke,i,j+1) )
1198  call check( __line__, pott(ke,i,j+2) )
1199  call check( __line__, gsqrt(ke,i,j,i_uyz) )
1200  call check( __line__, gsqrt(ke,i,j-1,i_xvz) )
1201  call check( __line__, gsqrt(ke,i,j ,i_xvz) )
1202  call check( __line__, gsqrt(ke,i,j,i_xyz) )
1203  call check( __line__, st(ke,i,j) )
1204  call check( __line__, dpres(ke-1,i,j) )
1205  call check( __line__, dpres(ke ,i,j) )
1206  call check( __line__, rt2p(ke-1,i,j) )
1207  call check( __line__, rt2p(ke ,i,j) )
1208  call check( __line__, rhot(ke-1,i,j) )
1209  call check( __line__, ddens(ke-1,i,j) )
1210 #endif
1211  b(ke,i,j) = ( &
1212  ( &
1213  - j33g * ( momz(ke-1,i,j) + dtrk*sw(ke-1,i,j) ) &
1214  * 0.5_rp*(pott(ke ,i,j)+pott(ke-1,i,j)) ) * rcdz(ke) &
1215  + ( gsqrt(ke,i,j ,i_xvz) * ( momy(ke,i,j ) + dtrk*sv(ke,i,j ) ) &
1216  * ( fact_n*(pott(ke,i,j+1)+pott(ke,i,j )) &
1217  + fact_f*(pott(ke,i,j+2)+pott(ke,i,j-1)) ) &
1218  - gsqrt(ke,i,j-1,i_xvz) * ( momy(ke,i,j-1) + dtrk*sv(ke,i,j-1) ) &
1219  * ( fact_n*(pott(ke,i,j )+pott(ke,i,j-1)) &
1220  + fact_f*(pott(ke,i,j+1)+pott(ke,i,j-2)) ) ) * rcdy(j) &
1221  + gsqrt(ke,i,j,i_xyz) * ( st(ke,i,j) - dpres(ke,i,j) * rcs2t(ke,i,j) * rdt ) &
1222  ) * rdt &
1223  + grav * j33g * 0.5_rp * ( - ( dpres(ke,i,j)+dpres(ke-1,i,j) ) * rcs2t(ke,i,j) &
1224  + ( ddens(ke,i,j)+ddens(ke-1,i,j) ) * pott(ke,i,j) ) * rcdz(ke)
1225  enddo
1226  enddo
1227 #ifdef DEBUG
1228  k = iundef; i = iundef; j = iundef
1229 #endif
1230 
1231  else
1232 
1233 
1234  ! r = b - M x0
1235  do j = jjs, jje
1236  do i = iis, iie
1237  do k = ks+2, ke-2
1238 #ifdef DEBUG
1239  call check( __line__, momz(k-1,i,j) )
1240  call check( __line__, momz(k ,i,j) )
1241  call check( __line__, momx(k,i-1,j) )
1242  call check( __line__, momx(k,i ,j) )
1243  call check( __line__, momy(k,i,j-1) )
1244  call check( __line__, momy(k,i,j ) )
1245  call check( __line__, pott(k-2,i,j) )
1246  call check( __line__, pott(k-1,i,j) )
1247  call check( __line__, pott(k ,i,j) )
1248  call check( __line__, pott(k+1,i,j) )
1249  call check( __line__, pott(k+2,i,j) )
1250  call check( __line__, pott(k,i-2,j) )
1251  call check( __line__, pott(k,i-1,j) )
1252  call check( __line__, pott(k,i ,j) )
1253  call check( __line__, pott(k,i+1,j) )
1254  call check( __line__, pott(k,i+2,j) )
1255  call check( __line__, pott(k,i,j-2) )
1256  call check( __line__, pott(k,i,j-1) )
1257  call check( __line__, pott(k,i,j ) )
1258  call check( __line__, pott(k,i,j+1) )
1259  call check( __line__, pott(k,i,j+2) )
1260  call check( __line__, sw(k-1,i,j) )
1261  call check( __line__, sw(k ,i,j) )
1262  call check( __line__, su(k,i-1,j) )
1263  call check( __line__, su(k,i ,j) )
1264  call check( __line__, sv(k,i,j-1) )
1265  call check( __line__, sv(k,i,j ) )
1266  call check( __line__, gsqrt(k,i ,j,i_uyz) )
1267  call check( __line__, gsqrt(k,i-1,j,i_uyz) )
1268  call check( __line__, gsqrt(k,i,j-1,i_xvz) )
1269  call check( __line__, gsqrt(k,i,j ,i_xvz) )
1270  call check( __line__, gsqrt(k,i,j,i_xyz) )
1271  call check( __line__, st(k,i,j) )
1272  call check( __line__, dpres(k-1,i,j) )
1273  call check( __line__, dpres(k ,i,j) )
1274  call check( __line__, dpres(k+1,i,j) )
1275  call check( __line__, rt2p(k-1,i,j) )
1276  call check( __line__, rt2p(k ,i,j) )
1277  call check( __line__, rt2p(k+1,i,j) )
1278  call check( __line__, rhot(k-1,i,j) )
1279  call check( __line__, rhot(k+1,i,j) )
1280  call check( __line__, ddens(k-1,i,j) )
1281  call check( __line__, ddens(k+1,i,j) )
1282 #endif
1283  b(k,i,j) = ( &
1284  ( j33g * ( momz(k ,i,j) + dtrk*sw(k ,i,j) ) &
1285  * ( fact_n*(pott(k+1,i,j)+pott(k ,i,j)) &
1286  + fact_f*(pott(k+2,i,j)+pott(k-1,i,j)) ) &
1287  - j33g * ( momz(k-1,i,j) + dtrk*sw(k-1,i,j) ) &
1288  * ( fact_n*(pott(k ,i,j)+pott(k-1,i,j)) &
1289  + fact_f*(pott(k+1,i,j)+pott(k-2,i,j)) ) ) * rcdz(k) &
1290  + ( gsqrt(k,i ,j,i_uyz) * ( momx(k,i ,j) + dtrk*su(k,i ,j) ) &
1291  * ( fact_n*(pott(k,i+1,j)+pott(k,i ,j)) &
1292  + fact_f*(pott(k,i+2,j)+pott(k,i-1,j)) ) &
1293  - gsqrt(k,i-1,j,i_uyz) * ( momx(k,i-1,j) + dtrk*su(k,i-1,j) ) &
1294  * ( fact_n*(pott(k,i ,j)+pott(k,i-1,j)) &
1295  + fact_f*(pott(k,i+1,j)+pott(k,i-2,j)) ) ) * rcdx(i) &
1296  + ( gsqrt(k,i,j ,i_xvz) * ( momy(k,i,j ) + dtrk*sv(k,i,j ) ) &
1297  * ( fact_n*(pott(k,i,j+1)+pott(k,i,j )) &
1298  + fact_f*(pott(k,i,j+2)+pott(k,i,j-1)) ) &
1299  - gsqrt(k,i,j-1,i_xvz) * ( momy(k,i,j-1) + dtrk*sv(k,i,j-1) ) &
1300  * ( fact_n*(pott(k,i,j )+pott(k,i,j-1)) &
1301  + fact_f*(pott(k,i,j+1)+pott(k,i,j-2)) ) ) * rcdy(j) &
1302  + gsqrt(k,i,j,i_xyz) * ( st(k,i,j) - dpres(k,i,j) * rcs2t(k,i,j) * rdt ) &
1303  ) * rdt &
1304  + grav * j33g * ( ( dpres(k+1,i,j)*rcs2t(k+1,i,j) &
1305  - dpres(k-1,i,j)*rcs2t(k-1,i,j) ) &
1306  - ( rhot(k+1,i,j)*ddens(k+1,i,j)/dens(k+1,i,j) &
1307  - rhot(k-1,i,j)*ddens(k-1,i,j)/dens(k-1,i,j) ) &
1308  ) / ( fdz(k) + fdz(k-1) )
1309  enddo
1310  enddo
1311  enddo
1312 #ifdef DEBUG
1313  k = iundef; i = iundef; j = iundef
1314 #endif
1315  do j = jjs, jje
1316  do i = iis, iie
1317 #ifdef DEBUG
1318  call check( __line__, momz(ks,i,j) )
1319  call check( __line__, momx(ks,i,j) )
1320  call check( __line__, momy(ks,i,j) )
1321  call check( __line__, sw(ks,i,j) )
1322  call check( __line__, su(ks,i-1,j) )
1323  call check( __line__, su(ks,i ,j) )
1324  call check( __line__, sv(ks,i,j-1) )
1325  call check( __line__, sv(ks,i,j ) )
1326  call check( __line__, st(ks,i,j) )
1327  call check( __line__, pott(ks ,i,j) )
1328  call check( __line__, pott(ks+1,i,j) )
1329  call check( __line__, pott(ks,i-2,j) )
1330  call check( __line__, pott(ks,i-1,j) )
1331  call check( __line__, pott(ks,i ,j) )
1332  call check( __line__, pott(ks,i+1,j) )
1333  call check( __line__, pott(ks,i+2,j) )
1334  call check( __line__, pott(ks,i,j-2) )
1335  call check( __line__, pott(ks,i,j-1) )
1336  call check( __line__, pott(ks,i,j ) )
1337  call check( __line__, pott(ks,i,j+1) )
1338  call check( __line__, pott(ks,i,j+2) )
1339  call check( __line__, gsqrt(ks,i,j,i_uyz) )
1340  call check( __line__, gsqrt(ks,i,j,i_xvz) )
1341  call check( __line__, gsqrt(ks,i,j,i_xyz) )
1342  call check( __line__, dpres(ks ,i,j) )
1343  call check( __line__, dpres(ks+1,i,j) )
1344  call check( __line__, rt2p(ks ,i,j) )
1345  call check( __line__, rt2p(ks+1,i,j) )
1346  call check( __line__, rhot(ks+1,i,j) )
1347  call check( __line__, ddens(ks+1,i,j) )
1348 #endif
1349  b(ks,i,j) = ( &
1350  ( j33g * ( momz(ks ,i,j) + dtrk*sw(ks ,i,j) ) &
1351  * 0.5_rp*(pott(ks+1,i,j)+pott(ks ,i,j)) ) * rcdz(ks) &
1352  + ( gsqrt(ks,i ,j,i_uyz) * ( momx(ks,i ,j) + dtrk*su(ks,i ,j) ) &
1353  * ( fact_n*(pott(ks,i+1,j)+pott(ks,i ,j)) &
1354  + fact_f*(pott(ks,i+2,j)+pott(ks,i-1,j)) ) &
1355  - gsqrt(ks,i-1,j,i_uyz) * ( momx(ks,i-1,j) + dtrk*su(ks,i-1,j) ) &
1356  * ( fact_n*(pott(ks,i ,j)+pott(ks,i-1,j)) &
1357  + fact_f*(pott(ks,i+1,j)+pott(ks,i-2,j)) ) ) * rcdx(i) &
1358  + ( gsqrt(ks,i,j ,i_xvz) * ( momy(ks,i,j ) + dtrk*sv(ks,i,j ) ) &
1359  * ( fact_n*(pott(ks,i,j+1)+pott(ks,i,j )) &
1360  + fact_f*(pott(ks,i,j+2)+pott(ks,i,j-1)) ) &
1361  - gsqrt(ks,i,j-1,i_xvz) * ( momy(ks,i,j-1) + dtrk*sv(ks,i,j-1) ) &
1362  * ( fact_n*(pott(ks,i,j )+pott(ks,i,j-1)) &
1363  + fact_f*(pott(ks,i,j+1)+pott(ks,i,j-2)) ) ) * rcdy(j) &
1364  + gsqrt(ks,i,j,i_xyz) * ( st(ks,i,j) - dpres(ks,i,j) * rcs2t(ks,i,j) * rdt ) &
1365  ) * rdt &
1366  + grav * j33g * 0.5_rp * ( ( dpres(ks,i,j)+dpres(ks+1,i,j) ) * rcs2t(ks,i,j) &
1367  - ( ddens(ks,i,j)+ddens(ks+1,i,j) ) * pott(ks,i,j) ) * rcdz(ks)
1368 #ifdef DEBUG
1369  call check( __line__, momz(ks ,i,j) )
1370  call check( __line__, momz(ks+1 ,i,j) )
1371  call check( __line__, momx(ks+1,i-1,j) )
1372  call check( __line__, momx(ks+1,i ,j) )
1373  call check( __line__, momy(ks+1,i,j-1) )
1374  call check( __line__, momy(ks+1,i,j ) )
1375  call check( __line__, pott(ks ,i,j) )
1376  call check( __line__, pott(ks+1 ,i,j) )
1377  call check( __line__, pott(ks+1+1,i,j) )
1378  call check( __line__, pott(ks+1+2,i,j) )
1379  call check( __line__, pott(ks+1,i-2,j) )
1380  call check( __line__, pott(ks+1,i-1,j) )
1381  call check( __line__, pott(ks+1,i ,j) )
1382  call check( __line__, pott(ks+1,i+1,j) )
1383  call check( __line__, pott(ks+1,i+2,j) )
1384  call check( __line__, pott(ks+1,i,j-2) )
1385  call check( __line__, pott(ks+1,i,j-1) )
1386  call check( __line__, pott(ks+1,i,j ) )
1387  call check( __line__, pott(ks+1,i,j+1) )
1388  call check( __line__, pott(ks+1,i,j+2) )
1389  call check( __line__, sw(ks ,i,j) )
1390  call check( __line__, sw(ks+1,i,j) )
1391  call check( __line__, su(ks+1,i-1,j) )
1392  call check( __line__, su(ks+1,i ,j) )
1393  call check( __line__, sv(ks+1,i,j-1) )
1394  call check( __line__, sv(ks+1,i,j ) )
1395  call check( __line__, gsqrt(ks+1,i ,j,i_uyz) )
1396  call check( __line__, gsqrt(ks+1,i-1,j,i_uyz) )
1397  call check( __line__, gsqrt(ks+1,i,j-1,i_xvz) )
1398  call check( __line__, gsqrt(ks+1,i,j ,i_xvz) )
1399  call check( __line__, gsqrt(ks+1,i,j,i_xyz) )
1400  call check( __line__, st(ks+1,i,j) )
1401  call check( __line__, dpres(ks ,i,j) )
1402  call check( __line__, dpres(ks+1,i,j) )
1403  call check( __line__, dpres(ks+2,i,j) )
1404  call check( __line__, rt2p(ks ,i,j) )
1405  call check( __line__, rt2p(ks+1,i,j) )
1406  call check( __line__, rt2p(ks+2,i,j) )
1407  call check( __line__, rhot(ks ,i,j) )
1408  call check( __line__, rhot(ks+2,i,j) )
1409  call check( __line__, ddens(ks ,i,j) )
1410  call check( __line__, ddens(ks+2,i,j) )
1411 #endif
1412  b(ks+1,i,j) = ( &
1413  ( j33g * ( momz(ks+1,i,j) + dtrk*sw(ks+1,i,j) ) &
1414  * ( fact_n*(pott(ks+2,i,j)+pott(ks+1,i,j)) &
1415  + fact_f*(pott(ks+3,i,j)+pott(ks ,i,j)) ) &
1416  - j33g * ( momz(ks+1-1,i,j) + dtrk*sw(ks+1-1,i,j) ) &
1417  * ( 0.5_rp*(pott(ks+1,i,j)+pott(ks,i,j)) ) ) * rcdz(ks+1) &
1418  + ( gsqrt(ks+1,i ,j,i_uyz) * ( momx(ks+1,i ,j) + dtrk*su(ks+1,i ,j) ) &
1419  * ( fact_n*(pott(ks+1,i+1,j)+pott(ks+1,i ,j)) &
1420  + fact_f*(pott(ks+1,i+2,j)+pott(ks+1,i-1,j)) ) &
1421  - gsqrt(ks+1,i-1,j,i_uyz) * ( momx(ks+1,i-1,j) + dtrk*su(ks+1,i-1,j) ) &
1422  * ( fact_n*(pott(ks+1,i ,j)+pott(ks+1,i-1,j)) &
1423  + fact_f*(pott(ks+1,i+1,j)+pott(ks+1,i-2,j)) ) ) * rcdx(i) &
1424  + ( gsqrt(ks+1,i,j ,i_xvz) * ( momy(ks+1,i,j ) + dtrk*sv(ks+1,i,j ) ) &
1425  * ( fact_n*(pott(ks+1,i,j+1)+pott(ks+1,i,j )) &
1426  + fact_f*(pott(ks+1,i,j+2)+pott(ks+1,i,j-1)) ) &
1427  - gsqrt(ks+1,i,j-1,i_xvz) * ( momy(ks+1,i,j-1) + dtrk*sv(ks+1,i,j-1) ) &
1428  * ( fact_n*(pott(ks+1,i,j )+pott(ks+1,i,j-1)) &
1429  + fact_f*(pott(ks+1,i,j+1)+pott(ks+1,i,j-2)) ) ) * rcdy(j) &
1430  + gsqrt(ks+1,i,j,i_xyz) * ( st(ks+1,i,j) - dpres(ks+1,i,j) * rcs2t(ks+1,i,j) * rdt ) &
1431  ) * rdt &
1432  + grav * j33g * ( ( dpres(ks+2,i,j)*rcs2t(ks+2,i,j) &
1433  - dpres(ks ,i,j)*rcs2t(ks ,i,j) ) &
1434  - ( rhot(ks+2,i,j)*ddens(ks+2,i,j)/dens(ks+2,i,j) &
1435  - rhot(ks ,i,j)*ddens(ks ,i,j)/dens(ks ,i,j) ) &
1436  ) / ( fdz(ks+1) + fdz(ks) )
1437 #ifdef DEBUG
1438  call check( __line__, momz(ke-2,i,j) )
1439  call check( __line__, momz(ke-1,i,j) )
1440  call check( __line__, momx(ke-1,i-1,j) )
1441  call check( __line__, momx(ke-1,i ,j) )
1442  call check( __line__, momy(ke-1,i,j-1) )
1443  call check( __line__, momy(ke-1,i,j ) )
1444  call check( __line__, pott(ke-3,i,j) )
1445  call check( __line__, pott(ke-2,i,j) )
1446  call check( __line__, pott(ke-1,i,j) )
1447  call check( __line__, pott(ke ,i,j) )
1448  call check( __line__, pott(ke-1,i-2,j) )
1449  call check( __line__, pott(ke-1,i-1,j) )
1450  call check( __line__, pott(ke-1,i ,j) )
1451  call check( __line__, pott(ke-1,i+1,j) )
1452  call check( __line__, pott(ke-1,i+2,j) )
1453  call check( __line__, pott(ke-1,i,j-2) )
1454  call check( __line__, pott(ke-1,i,j-1) )
1455  call check( __line__, pott(ke-1,i,j ) )
1456  call check( __line__, pott(ke-1,i,j+1) )
1457  call check( __line__, pott(ke-1,i,j+2) )
1458  call check( __line__, sw(ke-2,i,j) )
1459  call check( __line__, sw(ke-1,i,j) )
1460  call check( __line__, su(ke-1,i-1,j) )
1461  call check( __line__, su(ke-1,i ,j) )
1462  call check( __line__, sv(ke-1,i,j-1) )
1463  call check( __line__, sv(ke-1,i,j ) )
1464  call check( __line__, gsqrt(ke-1,i ,j,i_uyz) )
1465  call check( __line__, gsqrt(ke-1,i-1,j,i_uyz) )
1466  call check( __line__, gsqrt(ke-1,i,j-1,i_xvz) )
1467  call check( __line__, gsqrt(ke-1,i,j ,i_xvz) )
1468  call check( __line__, gsqrt(ke-1,i,j,i_xyz) )
1469  call check( __line__, st(ke-1,i,j) )
1470  call check( __line__, dpres(ke-2,i,j) )
1471  call check( __line__, dpres(ke-1,i,j) )
1472  call check( __line__, dpres(ke ,i,j) )
1473  call check( __line__, rt2p(ke-2,i,j) )
1474  call check( __line__, rt2p(ke-1,i,j) )
1475  call check( __line__, rt2p(ke ,i,j) )
1476  call check( __line__, rhot(ke-2,i,j) )
1477  call check( __line__, rhot(ke,i,j) )
1478  call check( __line__, ddens(ke-2,i,j) )
1479  call check( __line__, ddens(ke,i,j) )
1480 #endif
1481  b(ke-1,i,j) = ( &
1482  ( j33g * ( momz(ke-1,i,j) + dtrk*sw(ke-1,i,j) ) &
1483  * ( 0.5_rp*(pott(ke ,i,j)+pott(ke-1,i,j)) ) &
1484  - j33g * ( momz(ke-2,i,j) + dtrk*sw(ke-2,i,j) ) &
1485  * ( fact_n*(pott(ke-1,i,j)+pott(ke-2,i,j)) &
1486  + fact_f*(pott(ke ,i,j)+pott(ke-3,i,j)) ) ) * rcdz(ke-1) &
1487  + ( gsqrt(ke-1,i ,j,i_uyz) * ( momx(ke-1,i ,j) + dtrk*su(ke-1,i ,j) ) &
1488  * ( fact_n*(pott(ke-1,i+1,j)+pott(ke-1,i ,j)) &
1489  + fact_f*(pott(ke-1,i+2,j)+pott(ke-1,i-1,j)) ) &
1490  - gsqrt(ke-1,i-1,j,i_uyz) * ( momx(ke-1,i-1,j) + dtrk*su(ke-1,i-1,j) ) &
1491  * ( fact_n*(pott(ke-1,i ,j)+pott(ke-1,i-1,j)) &
1492  + fact_f*(pott(ke-1,i+1,j)+pott(ke-1,i-2,j)) ) ) * rcdx(i) &
1493  + ( gsqrt(ke-1,i,j ,i_xvz) * ( momy(ke-1,i,j ) + dtrk*sv(ke-1,i,j ) ) &
1494  * ( fact_n*(pott(ke-1,i,j+1)+pott(ke-1,i,j )) &
1495  + fact_f*(pott(ke-1,i,j+2)+pott(ke-1,i,j-1)) ) &
1496  - gsqrt(ke-1,i,j-1,i_xvz) * ( momy(ke-1,i,j-1) + dtrk*sv(ke-1,i,j-1) ) &
1497  * ( fact_n*(pott(ke-1,i,j )+pott(ke-1,i,j-1)) &
1498  + fact_f*(pott(ke-1,i,j+1)+pott(ke-1,i,j-2)) ) ) * rcdy(j) &
1499  + gsqrt(ke-1,i,j,i_xyz) * ( st(ke-1,i,j) - dpres(ke-1,i,j) * rcs2t(ke-1,i,j) * rdt ) &
1500  ) * rdt &
1501  + grav * j33g * ( ( dpres(ke ,i,j)*rcs2t(ke ,i,j) &
1502  - dpres(ke-2,i,j)*rcs2t(ke-2,i,j) ) &
1503  - ( rhot(ke ,i,j)*ddens(ke ,i,j)/dens(ke ,i,j) &
1504  - rhot(ke-2,i,j)*ddens(ke-2,i,j)/dens(ke-2,i,j) )&
1505  ) / ( fdz(ke-1) + fdz(ke-1-1) )
1506 #ifdef DEBUG
1507  call check( __line__, momz(ke-1,i,j) )
1508  call check( __line__, momx(ke,i-1,j) )
1509  call check( __line__, momx(ke,i ,j) )
1510  call check( __line__, momy(ke,i,j-1) )
1511  call check( __line__, momy(ke,i,j ) )
1512  call check( __line__, sw(ke-1,i,j) )
1513  call check( __line__, su(ke,i-1,j) )
1514  call check( __line__, su(ke,i ,j) )
1515  call check( __line__, sv(ke,i,j-1) )
1516  call check( __line__, sv(ke,i,j ) )
1517  call check( __line__, gsqrt(ke,i-1,j,i_uyz) )
1518  call check( __line__, pott(ke-1,i,j) )
1519  call check( __line__, pott(ke ,i,j) )
1520  call check( __line__, pott(ke,i-2,j) )
1521  call check( __line__, pott(ke,i-1,j) )
1522  call check( __line__, pott(ke,i ,j) )
1523  call check( __line__, pott(ke,i+1,j) )
1524  call check( __line__, pott(ke,i+2,j) )
1525  call check( __line__, pott(ke,i,j-2) )
1526  call check( __line__, pott(ke,i,j-1) )
1527  call check( __line__, pott(ke,i,j ) )
1528  call check( __line__, pott(ke,i,j+1) )
1529  call check( __line__, pott(ke,i,j+2) )
1530  call check( __line__, gsqrt(ke,i-1,j,i_uyz) )
1531  call check( __line__, gsqrt(ke,i ,j,i_uyz) )
1532  call check( __line__, gsqrt(ke,i,j-1,i_xvz) )
1533  call check( __line__, gsqrt(ke,i,j ,i_xvz) )
1534  call check( __line__, gsqrt(ke,i,j,i_xyz) )
1535  call check( __line__, st(ke,i,j) )
1536  call check( __line__, dpres(ke-1,i,j) )
1537  call check( __line__, dpres(ke ,i,j) )
1538  call check( __line__, rt2p(ke-1,i,j) )
1539  call check( __line__, rt2p(ke ,i,j) )
1540  call check( __line__, rhot(ke-1,i,j) )
1541  call check( __line__, ddens(ke-1,i,j) )
1542 #endif
1543  b(ke,i,j) = ( &
1544  ( &
1545  - j33g * ( momz(ke-1,i,j) + dtrk*sw(ke-1,i,j) ) &
1546  * 0.5_rp*(pott(ke ,i,j)+pott(ke-1,i,j)) ) * rcdz(ke) &
1547  + ( gsqrt(ke,i ,j,i_uyz) * ( momx(ke,i ,j) + dtrk*su(ke,i ,j) ) &
1548  * ( fact_n*(pott(ke,i+1,j)+pott(ke,i ,j)) &
1549  + fact_f*(pott(ke,i+2,j)+pott(ke,i-1,j)) ) &
1550  - gsqrt(ke,i-1,j,i_uyz) * ( momx(ke,i-1,j) + dtrk*su(ke,i-1,j) ) &
1551  * ( fact_n*(pott(ke,i ,j)+pott(ke,i-1,j)) &
1552  + fact_f*(pott(ke,i+1,j)+pott(ke,i-2,j)) ) ) * rcdx(i) &
1553  + ( gsqrt(ke,i,j ,i_xvz) * ( momy(ke,i,j ) + dtrk*sv(ke,i,j ) ) &
1554  * ( fact_n*(pott(ke,i,j+1)+pott(ke,i,j )) &
1555  + fact_f*(pott(ke,i,j+2)+pott(ke,i,j-1)) ) &
1556  - gsqrt(ke,i,j-1,i_xvz) * ( momy(ke,i,j-1) + dtrk*sv(ke,i,j-1) ) &
1557  * ( fact_n*(pott(ke,i,j )+pott(ke,i,j-1)) &
1558  + fact_f*(pott(ke,i,j+1)+pott(ke,i,j-2)) ) ) * rcdy(j) &
1559  + gsqrt(ke,i,j,i_xyz) * ( st(ke,i,j) - dpres(ke,i,j) * rcs2t(ke,i,j) * rdt ) &
1560  ) * rdt &
1561  + grav * j33g * 0.5_rp * ( - ( dpres(ke,i,j)+dpres(ke-1,i,j) ) * rcs2t(ke,i,j) &
1562  + ( ddens(ke,i,j)+ddens(ke-1,i,j) ) * pott(ke,i,j) ) * rcdz(ke)
1563  enddo
1564  enddo
1565 #ifdef DEBUG
1566  k = iundef; i = iundef; j = iundef
1567 #endif
1568 
1569  end if
1570 
1571  call make_matrix( &
1572  m, & ! (inout)
1573  pott, rcs2t, grav, &
1574  gsqrt, j33g, &
1575  rcdz, rfdz, rcdx, rfdx, rcdy,rfdy, fdz, &
1576  rdt, &
1577  fact_n, fact_f, &
1578  twod, &
1579  i_xyz, i_xyw, &
1580  iis, iie, jjs, jje )
1581 
1582  enddo
1583  enddo
1584 
1585  call solve_bicgstab( &
1586  dpres_n, & ! (out)
1587  dpres, & ! (in)
1588  m, b, & ! (in)
1589  twod ) ! (in)
1590 
1591  call comm_vars8( dpres_n, 1 )
1592  call comm_wait ( dpres_n, 1 )
1593 
1594 #ifdef DEBUG
1595  call check_solver( dpres_n, m, b, twod )
1596 #endif
1597 
1598 #else
1599 
1600  call solve_multigrid
1601 
1602 #endif
1603 
1604  do jjs = js, je, jblock
1605  jje = jjs+jblock-1
1606  do iis = is, ie, iblock
1607  iie = iis+iblock-1
1608 
1609  !##### momentum equation (z) #####
1610 
1611  !--- update momentum(z)
1612  !$omp parallel do private(i,j,k,duvw) OMP_SCHEDULE_ collapse(2)
1613  do j = jjs, jje
1614  do i = iis, iie
1615  do k = ks, ke-1
1616 #ifdef DEBUG
1617  call check( __line__, dpres_n(k+1,i,j) )
1618  call check( __line__, dpres_n(k ,i,j) )
1619  call check( __line__, dpres(k+1,i,j) )
1620  call check( __line__, dpres(k ,i,j) )
1621  call check( __line__, ddens(k+1,i,j) )
1622  call check( __line__, ddens(k ,i,j) )
1623  call check( __line__, ref_dens(k+1,i,j) )
1624  call check( __line__, ref_dens(k,i,j) )
1625  call check( __line__, momz0(k,i,j) )
1626 #endif
1627  duvw = dtrk * ( ( &
1628  - j33g * ( dpres_n(k+1,i,j) - dpres_n(k,i,j) ) * rfdz(k) &
1629  ) / gsqrt(k,i,j,i_uyz) &
1630  - 0.5_rp * grav &
1631  * ( ddens(k+1,i,j) &
1632  + ( dpres_n(k+1,i,j) - dpres(k+1,i,j) ) &
1633  / ( pott(k+1,i,j) * rt2p(k+1,i,j) ) &
1634  + ddens(k ,i,j) &
1635  + ( dpres_n(k ,i,j) - dpres(k ,i,j) ) &
1636  / ( pott(k ,i,j) * rt2p(k ,i,j) ) ) &
1637  + sw(k,i,j) )
1638  momz_rk(k,i,j) = momz0(k,i,j) + duvw
1639  mflx_hi(k,i,j,zdir) = j33g * ( momz(k,i,j) + duvw )
1640  enddo
1641  enddo
1642  enddo
1643 #ifdef DEBUG
1644  k = iundef; i = iundef; j = iundef
1645 #endif
1646  do j = jjs, jje
1647  do i = iis, iie
1648  mflx_hi(ks-1,i,j,zdir) = 0.0_rp
1649  mflx_hi(ke ,i,j,zdir) = 0.0_rp
1650  enddo
1651  enddo
1652 #ifdef DEBUG
1653  k = iundef; i = iundef; j = iundef
1654 #endif
1655 
1656 
1657  !##### momentum equation (x) #####
1658 
1659  !--- update momentum(x)
1660  if ( twod ) then
1661  !$omp parallel do private(j,k,duvw) OMP_SCHEDULE_
1662  do j = jjs , jje
1663  do k = ks, ke
1664 #ifdef DEBUG
1665  call check( __line__, su(k,is,j) )
1666  call check( __line__, momx0(k,is,j) )
1667 #endif
1668  duvw = dtrk * su(k,is,j)
1669 
1670  momx_rk(k,is,j) = momx0(k,is,j) + duvw
1671  enddo
1672  enddo
1673  else
1674  !$omp parallel do private(i,j,k,duvw) OMP_SCHEDULE_ collapse(2)
1675  do j = jjs , jje
1676  do i = iis-1, iie
1677  do k = ks, ke
1678 #ifdef DEBUG
1679  call check( __line__, dpres_n(k,i+1,j) )
1680  call check( __line__, dpres_n(k,i ,j) )
1681  call check( __line__, gsqrt(k,i+1,j,i_xyz) )
1682  call check( __line__, gsqrt(k,i ,j,i_xyz) )
1683  call check( __line__, gsqrt(k,i,j,i_uyz) )
1684  call check( __line__, su(k,i,j) )
1685  call check( __line__, momx0(k,i,j) )
1686 #endif
1687  duvw = dtrk * ( ( &
1688  - ( gsqrt(k,i+1,j,i_xyz) * dpres_n(k,i+1,j) &
1689  - gsqrt(k,i ,j,i_xyz) * dpres_n(k,i ,j) ) * rfdx(i) & ! pressure gradient force
1690  ) / gsqrt(k,i,j,i_uyz) &
1691  + su(k,i,j) )
1692 
1693  momx_rk(k,i,j) = momx0(k,i,j) + duvw
1694  mflx_hi(k,i,j,xdir) = gsqrt(k,i,j,i_uyz) * ( momx(k,i,j) + duvw )
1695  enddo
1696  enddo
1697  enddo
1698  end if
1699 #ifdef DEBUG
1700  k = iundef; i = iundef; j = iundef
1701 #endif
1702 
1703  !##### momentum equation (y) #####
1704 
1705  !--- update momentum(y)
1706  !$omp parallel do private(i,j,k,duvw) OMP_SCHEDULE_ collapse(2)
1707  do j = jjs-1, jje
1708  do i = iis, iie
1709  do k = ks, ke
1710 #ifdef DEBUG
1711  call check( __line__, dpres_n(k,i,j ) )
1712  call check( __line__, dpres_n(k,i,j+1) )
1713  call check( __line__, gsqrt(k,i,j ,i_xyz) )
1714  call check( __line__, gsqrt(k,i,j+1,i_xyz) )
1715  call check( __line__, gsqrt(k,i,j,i_xvz) )
1716  call check( __line__, sv(k,i,j) )
1717  call check( __line__, momy0(k,i,j) )
1718 #endif
1719  duvw = dtrk * ( ( &
1720  - ( gsqrt(k,i,j+1,i_xyz) * dpres_n(k,i,j+1) &
1721  - gsqrt(k,i,j ,i_xyz) * dpres_n(k,i,j ) ) * rfdy(j) &
1722  ) / gsqrt(k,i,j,i_xvz) &
1723  + sv(k,i,j) )
1724  momy_rk(k,i,j) = momy0(k,i,j) + duvw
1725  mflx_hi(k,i,j,ydir) = gsqrt(k,i,j,i_xvz) * ( momy(k,i,j) + duvw )
1726  enddo
1727  enddo
1728  enddo
1729 #ifdef DEBUG
1730  k = iundef; i = iundef; j = iundef
1731 #endif
1732 
1733 
1734  !##### Thermodynamic Equation #####
1735 
1736  ! at (x, y, w)
1737  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
1738  do j = jjs, jje
1739  do i = iis, iie
1740  do k = ks+1, ke-2
1741 #ifdef DEBUG
1742  call check( __line__, mflx_hi(k,i,j,zdir) )
1743  call check( __line__, pott(k,i-1,j) )
1744  call check( __line__, pott(k,i ,j) )
1745  call check( __line__, pott(k,i+1,j) )
1746  call check( __line__, pott(k,i+1,j) )
1747  call check( __line__, num_diff(k,i,j,i_rhot,xdir) )
1748 #endif
1749  tflx_hi(k,i,j,zdir) = mflx_hi(k,i,j,zdir) &
1750  * ( fact_n * ( pott(k+1,i,j) + pott(k ,i,j) ) &
1751  + fact_f * ( pott(k+2,i,j) + pott(k-1,i,j) ) )
1752  enddo
1753  enddo
1754  enddo
1755 #ifdef DEBUG
1756  k = iundef; i = iundef; j = iundef
1757 #endif
1758  !$omp parallel do private(i,j) OMP_SCHEDULE_ collapse(2)
1759  do j = jjs, jje
1760  do i = iis, iie
1761  tflx_hi(ks-1,i,j,zdir) = 0.0_rp
1762  tflx_hi(ks ,i,j,zdir) = mflx_hi(ks ,i,j,zdir) * 0.5_rp * ( pott(ks+1,i,j) + pott(ks ,i,j) )
1763  tflx_hi(ke-1,i,j,zdir) = mflx_hi(ke-1,i,j,zdir) * 0.5_rp * ( pott(ke ,i,j) + pott(ke-1,i,j) )
1764  tflx_hi(ke ,i,j,zdir) = 0.0_rp
1765  enddo
1766  enddo
1767 #ifdef DEBUG
1768  k = iundef; i = iundef; j = iundef
1769 #endif
1770  ! at (u, y, z)
1771  if ( .not. twod ) then
1772  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
1773  do j = jjs, jje
1774  do i = iis-1, iie
1775  do k = ks, ke
1776 #ifdef DEBUG
1777  call check( __line__, mflx_hi(k,i,j,xdir) )
1778  call check( __line__, pott(k,i-1,j) )
1779  call check( __line__, pott(k,i ,j) )
1780  call check( __line__, pott(k,i+1,j) )
1781  call check( __line__, pott(k,i+1,j) )
1782  call check( __line__, num_diff(k,i,j,i_rhot,xdir) )
1783 #endif
1784  tflx_hi(k,i,j,xdir) = mflx_hi(k,i,j,xdir) &
1785  * ( fact_n * ( pott(k,i+1,j)+pott(k,i ,j) ) &
1786  + fact_f * ( pott(k,i+2,j)+pott(k,i-1,j) ) )
1787  enddo
1788  enddo
1789  enddo
1790  end if
1791 #ifdef DEBUG
1792  k = iundef; i = iundef; j = iundef
1793 #endif
1794  ! at (x, v, z)
1795  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
1796  do j = jjs-1, jje
1797  do i = iis, iie
1798  do k = ks, ke
1799 #ifdef DEBUG
1800  call check( __line__, mflx_hi(k,i,j,ydir) )
1801  call check( __line__, pott(k,i,j-1) )
1802  call check( __line__, pott(k,i,j ) )
1803  call check( __line__, pott(k,i,j+1) )
1804  call check( __line__, pott(k,i,j+2) )
1805  call check( __line__, num_diff(k,i,j,i_rhot,ydir) )
1806 #endif
1807  tflx_hi(k,i,j,ydir) = mflx_hi(k,i,j,ydir) &
1808  * ( fact_n * ( pott(k,i,j+1)+pott(k,i,j ) ) &
1809  + fact_f * ( pott(k,i,j+2)+pott(k,i,j-1) ) )
1810  enddo
1811  enddo
1812  enddo
1813 #ifdef DEBUG
1814  k = iundef; i = iundef; j = iundef
1815 #endif
1816 
1817  ! rho*theta
1818  if ( twod ) then
1819  !$omp parallel do
1820  do j = jjs, jje
1821  do k = ks, ke
1822  rhot_rk(k,is,j) = rhot0(k,is,j) &
1823  + dtrk * ( - ( ( tflx_hi(k,is,j,zdir) - tflx_hi(k-1,is,j ,zdir) ) * rcdz(k) &
1824  + ( tflx_hi(k,is,j,ydir) - tflx_hi(k ,is,j-1,ydir) ) * rcdy(j) ) &
1825  / gsqrt(k,is,j,i_xyz) &
1826  + st(k,is,j) )
1827  enddo
1828  enddo
1829  else
1830  !$omp parallel do
1831  do j = jjs, jje
1832  do i = iis, iie
1833  do k = ks, ke
1834  rhot_rk(k,i,j) = rhot0(k,i,j) &
1835  + dtrk * ( - ( ( tflx_hi(k,i,j,zdir) - tflx_hi(k-1,i ,j ,zdir) ) * rcdz(k) &
1836  + ( tflx_hi(k,i,j,xdir) - tflx_hi(k ,i-1,j ,xdir) ) * rcdx(i) &
1837  + ( tflx_hi(k,i,j,ydir) - tflx_hi(k ,i ,j-1,ydir) ) * rcdy(j) ) &
1838  / gsqrt(k,i,j,i_xyz) &
1839  + st(k,i,j) )
1840  enddo
1841  enddo
1842  enddo
1843  end if
1844 
1845  !##### continuous equation #####
1846 
1847  ! total momentum flux
1848  !$omp parallel do
1849  do j = jjs, jje
1850  do i = iis, iie
1851  do k = ks, ke-1
1852  mflx_hi(k,i,j,zdir) = mflx_hi(k,i,j,zdir) + mflx_hi2(k,i,j,zdir)
1853  enddo
1854  enddo
1855  enddo
1856  if ( .not. twod ) then
1857  !$omp parallel do
1858  do j = jjs , jje
1859  do i = iis-1, iie
1860  do k = ks, ke
1861  mflx_hi(k,i,j,xdir) = mflx_hi(k,i,j,xdir) + mflx_hi2(k,i,j,xdir)
1862  enddo
1863  enddo
1864  enddo
1865  end if
1866  !$omp parallel do
1867  do j = jjs-1, jje
1868  do i = iis , iie
1869  do k = ks, ke
1870  mflx_hi(k,i,j,ydir) = mflx_hi(k,i,j,ydir) + mflx_hi2(k,i,j,ydir)
1871  enddo
1872  enddo
1873  enddo
1874 
1875  !--- update density
1876  if ( twod ) then
1877  !$omp parallel do private(j,k) OMP_SCHEDULE_
1878  do j = jjs, jje
1879  do k = ks, ke
1880 #ifdef DEBUG
1881  call check( __line__, dens0(k,is,j) )
1882  call check( __line__, mflx_hi(k ,is,j ,zdir) )
1883  call check( __line__, mflx_hi(k-1,is,j ,zdir) )
1884  call check( __line__, mflx_hi(k ,is,j ,ydir) )
1885  call check( __line__, mflx_hi(k ,is,j-1,ydir) )
1886  call check( __line__, dens_t(k,is,j) )
1887 #endif
1888  dens_rk(k,is,j) = dens0(k,is,j) &
1889  + dtrk * ( - ( ( mflx_hi(k,is,j,zdir)-mflx_hi(k-1,is,j, zdir) ) * rcdz(k) &
1890  + ( mflx_hi(k,is,j,ydir)-mflx_hi(k ,is,j-1,ydir) ) * rcdy(j) ) &
1891  / gsqrt(k,is,j,i_xyz) & ! divergence
1892  + dens_t(k,is,j) )
1893  enddo
1894  enddo
1895  else
1896  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
1897  do j = jjs, jje
1898  do i = iis, iie
1899  do k = ks, ke
1900 #ifdef DEBUG
1901  call check( __line__, dens0(k,i,j) )
1902  call check( __line__, mflx_hi(k ,i ,j ,zdir) )
1903  call check( __line__, mflx_hi(k-1,i ,j ,zdir) )
1904  call check( __line__, mflx_hi(k ,i ,j ,xdir) )
1905  call check( __line__, mflx_hi(k ,i-1,j ,xdir) )
1906  call check( __line__, mflx_hi(k ,i ,j ,ydir) )
1907  call check( __line__, mflx_hi(k ,i ,j-1,ydir) )
1908  call check( __line__, dens_t(k,i,j) )
1909 #endif
1910  dens_rk(k,i,j) = dens0(k,i,j) &
1911  + dtrk * ( - ( ( mflx_hi(k,i,j,zdir)-mflx_hi(k-1,i ,j, zdir) ) * rcdz(k) &
1912  + ( mflx_hi(k,i,j,xdir)-mflx_hi(k ,i-1,j, xdir) ) * rcdx(i) &
1913  + ( mflx_hi(k,i,j,ydir)-mflx_hi(k ,i, j-1,ydir) ) * rcdy(j) ) &
1914  / gsqrt(k,i,j,i_xyz) & ! divergence
1915  + dens_t(k,i,j) )
1916  enddo
1917  enddo
1918  enddo
1919  end if
1920 #ifdef DEBUG
1921  k = iundef; i = iundef; j = iundef
1922 #endif
1923 
1924 
1925 #ifdef DEBUG
1926  call check_pres( &
1927  dpres_n, dpres, rhot_rk, rhot, dens_rk, dens, b,&
1928  rt2p )
1929 #endif
1930 
1931  enddo
1932  enddo
1933 
1934  return
1935  end subroutine atmos_dyn_tstep_short_fvm_hivi
1936 
1937 #ifdef HIVI_BICGSTAB
1938  subroutine solve_bicgstab( &
1939  DPRES_N, &
1940  DPRES, &
1941  M, B, &
1942  TwoD )
1943  use scale_prc, only: &
1944  prc_abort
1945  use scale_comm_cartesc, only: &
1946  comm_world, &
1947  comm_vars8, &
1948  comm_wait
1949  implicit none
1950  real(RP), intent(out) :: DPRES_N(KA,IA,JA)
1951  real(RP), intent(in) :: DPRES(KA,IA,JA)
1952  real(RP), intent(in) :: M(7,KA,IA,JA)
1953  real(RP), intent(in) :: B(KA,IA,JA)
1954  logical, intent(in) :: TwoD
1955 
1956  real(RP) :: r0(KA,IA,JA)
1957 
1958  real(RP) :: p(KA,IA,JA)
1959  real(RP) :: Mp(KA,IA,JA)
1960  real(RP) :: s(KA,IA,JA)
1961  real(RP) :: Ms(KA,IA,JA)
1962  real(RP) :: al, be, w
1963 
1964  real(RP), pointer :: r(:,:,:)
1965  real(RP), pointer :: rn(:,:,:)
1966  real(RP), pointer :: swap(:,:,:)
1967  real(RP), target :: v0(KA,IA,JA)
1968  real(RP), target :: v1(KA,IA,JA)
1969  real(RP) :: r0r
1970  real(RP) :: norm, error, error2
1971 
1972  real(RP) :: iprod(2)
1973  real(RP) :: buf(2)
1974 
1975  integer :: k, i, j
1976  integer :: iis, iie, jjs, jje
1977  integer :: iter
1978  integer :: ierror
1979 
1980 #ifdef DEBUG
1981  r0(:,:,:) = undef
1982  p(:,:,:) = undef
1983  mp(:,:,:) = undef
1984  s(:,:,:) = undef
1985  ms(:,:,:) = undef
1986 
1987  v0(:,:,:) = undef
1988  v1(:,:,:) = undef
1989 #endif
1990 
1991  r => v0
1992  rn => v1
1993 
1994  call mul_matrix( v1, m, dpres, twod ) ! v1 = M x0
1995 
1996  norm = 0.0_rp
1997  r0r = 0.0_rp
1998 
1999  do jjs = js, je, jblock
2000  jje = jjs+jblock-1
2001  do iis = is, ie, iblock
2002  iie = iis+iblock-1
2003 
2004  do j = jjs, jje
2005  do i = iis, iie
2006  do k = ks, ke
2007  norm = norm + b(k,i,j)**2
2008  enddo
2009  enddo
2010  enddo
2011 #ifdef DEBUG
2012  k = iundef; i = iundef; j = iundef
2013 #endif
2014 
2015  ! r = b - M x0
2016  do j = jjs, jje
2017  do i = iis, iie
2018  do k = ks, ke
2019 #ifdef DEBUG
2020  call check( __line__, b(k,i,j) )
2021  call check( __line__, v1(k,i,j) )
2022 #endif
2023  r(k,i,j) = b(k,i,j) - v1(k,i,j)
2024  enddo
2025  enddo
2026  enddo
2027 #ifdef DEBUG
2028  k = iundef; i = iundef; j = iundef
2029 #endif
2030 
2031  do j = jjs, jje
2032  do i = iis, iie
2033  do k = ks, ke
2034  r0(k,i,j) = r(k,i,j)
2035  p(k,i,j) = r(k,i,j)
2036  enddo
2037  enddo
2038  enddo
2039 #ifdef DEBUG
2040  k = iundef; i = iundef; j = iundef
2041 #endif
2042 
2043  do j = jjs, jje
2044  do i = iis, iie
2045  do k = ks, ke
2046 #ifdef DEBUG
2047  call check( __line__, r(k,i,j) )
2048  call check( __line__, r0(k,i,j) )
2049 #endif
2050  r0r = r0r + r0(k,i,j) * r(k,i,j)
2051  enddo
2052  enddo
2053  enddo
2054 #ifdef DEBUG
2055  k = iundef; i = iundef; j = iundef
2056 #endif
2057 
2058  enddo
2059  enddo
2060 
2061  do j = js-1, je+1
2062  do i = is-1, ie+1
2063  do k = ks, ke
2064  dpres_n(k,i,j) = dpres(k,i,j)
2065  end do
2066  end do
2067  end do
2068 
2069 
2070  iprod(1) = r0r
2071  iprod(2) = norm
2072  call mpi_allreduce(iprod, buf, 2, mtype, mpi_sum, comm_world, ierror)
2073  r0r = buf(1)
2074  norm = buf(2)
2075 
2076  error2 = norm
2077 
2078  do iter = 1, itmax
2079 
2080  error = 0.0_rp
2081  do j = js, je
2082  do i = is, ie
2083  do k = ks, ke
2084 #ifdef DEBUG
2085  call check( __line__, r(k,i,j) )
2086 #endif
2087  error = error + r(k,i,j)**2
2088  enddo
2089  enddo
2090  enddo
2091 #ifdef DEBUG
2092  k = iundef; i = iundef; j = iundef
2093 #endif
2094  call mpi_allreduce(error, buf, 1, mtype, mpi_sum, comm_world, ierror)
2095  error = buf(1)
2096 
2097 #ifdef DEBUG
2098  log_info("solve_bicgstab",*) iter, error/norm
2099 #endif
2100  if ( sqrt(error/norm) < epsilon .OR. error > error2 ) then
2101 #ifdef DEBUG
2102  log_info("solve_bicgstab",*) "Bi-CGSTAB converged:", iter
2103 #endif
2104  exit
2105  endif
2106  error2 = error
2107 
2108  call comm_vars8( p, 1 )
2109  call comm_wait ( p, 1 )
2110  call mul_matrix( mp, m, p, twod )
2111 
2112  iprod(1) = 0.0_rp
2113  do j = js, je
2114  do i = is, ie
2115  do k = ks, ke
2116 #ifdef DEBUG
2117  call check( __line__, r0(k,i,j) )
2118  call check( __line__, mp(k,i,j) )
2119 #endif
2120  iprod(1) = iprod(1) + r0(k,i,j) * mp(k,i,j)
2121  enddo
2122  enddo
2123  enddo
2124 #ifdef DEBUG
2125  k = iundef; i = iundef; j = iundef
2126 #endif
2127  call mpi_allreduce(iprod, buf, 1, mtype, mpi_sum, comm_world, ierror)
2128  al = r0r / buf(1) ! (r0,r) / (r0,Mp)
2129 
2130  do j = js, je
2131  do i = is, ie
2132  do k = ks, ke
2133 #ifdef DEBUG
2134  call check( __line__, r(k,i,j) )
2135  call check( __line__, mp(k,i,j) )
2136 #endif
2137  s(k,i,j) = r(k,i,j) - al*mp(k,i,j)
2138  enddo
2139  enddo
2140  enddo
2141 #ifdef DEBUG
2142  k = iundef; i = iundef; j = iundef
2143 #endif
2144 
2145  call comm_vars8( s, 1 )
2146  call comm_wait ( s, 1 )
2147  call mul_matrix( ms, m, s, twod )
2148  iprod(1) = 0.0_rp
2149  iprod(2) = 0.0_rp
2150  do j = js, je
2151  do i = is, ie
2152  do k = ks, ke
2153 #ifdef DEBUG
2154  call check( __line__, ms(k,i,j) )
2155  call check( __line__, s(k,i,j) )
2156 #endif
2157  iprod(1) = iprod(1) + ms(k,i,j) * s(k,i,j)
2158  iprod(2) = iprod(2) + ms(k,i,j) * ms(k,i,j)
2159  enddo
2160  enddo
2161  enddo
2162 #ifdef DEBUG
2163  k = iundef; i = iundef; j = iundef
2164 #endif
2165  call mpi_allreduce(iprod, buf, 2, mtype, mpi_sum, comm_world, ierror)
2166  w = buf(1) / buf(2) ! (Ms,s) / (Ms,Ms)
2167 
2168  iprod(1) = 0.0_rp
2169 
2170  do jjs = js, je, jblock
2171  jje = jjs+jblock-1
2172  do iis = is, ie, iblock
2173  iie = iis+iblock-1
2174 
2175  do j = jjs, jje
2176  do i = iis, iie
2177  do k = ks, ke
2178 #ifdef DEBUG
2179  call check( __line__, dpres_n(k,i,j) )
2180  call check( __line__, p(k,i,j) )
2181  call check( __line__, s(k,i,j) )
2182 #endif
2183  dpres_n(k,i,j) = dpres_n(k,i,j) + al*p(k,i,j) + w*s(k,i,j)
2184  enddo
2185  enddo
2186  enddo
2187 #ifdef DEBUG
2188  k = iundef; i = iundef; j = iundef
2189 #endif
2190 
2191  do j = jjs, jje
2192  do i = iis, iie
2193  do k = ks, ke
2194 #ifdef DEBUG
2195  call check( __line__, s(k,i,j) )
2196  call check( __line__, ms(k,i,j) )
2197 #endif
2198  rn(k,i,j) = s(k,i,j) - w*ms(k,i,j)
2199  enddo
2200  enddo
2201  enddo
2202 #ifdef DEBUG
2203  k = iundef; i = iundef; j = iundef
2204 #endif
2205 
2206  do j = jjs, jje
2207  do i = iis, iie
2208  do k = ks, ke
2209 #ifdef DEBUG
2210  call check( __line__, r0(k,i,j) )
2211  call check( __line__, rn(k,i,j) )
2212 #endif
2213  iprod(1) = iprod(1) + r0(k,i,j) * rn(k,i,j)
2214  enddo
2215  enddo
2216  enddo
2217 #ifdef DEBUG
2218  k = iundef; i = iundef; j = iundef
2219 #endif
2220 
2221  enddo
2222  enddo
2223 
2224  be = al/w / r0r
2225 
2226  call mpi_allreduce(iprod, r0r, 1, mtype, mpi_sum, comm_world, ierror)
2227 
2228  be = be * r0r ! al/w * (r0,rn)/(r0,r)
2229  do j = js, je
2230  do i = is, ie
2231  do k = ks, ke
2232 #ifdef DEBUG
2233  call check( __line__, rn(k,i,j) )
2234  call check( __line__, p(k,i,j) )
2235  call check( __line__, mp(k,i,j) )
2236 #endif
2237  p(k,i,j) = rn(k,i,j) + be * ( p(k,i,j) - w*mp(k,i,j) )
2238  enddo
2239  enddo
2240  enddo
2241 #ifdef DEBUG
2242  k = iundef; i = iundef; j = iundef
2243 #endif
2244 
2245  swap => rn
2246  rn => r
2247  r => swap
2248 #ifdef DEBUG
2249  rn(:,:,:) = undef
2250 #endif
2251  enddo
2252 
2253  if ( iter >= itmax ) then
2254  log_error("solve_bicgstab",*) 'not converged', error, norm
2255  call prc_abort
2256  endif
2257 
2258  return
2259  end subroutine solve_bicgstab
2260 
2261  subroutine make_matrix(&
2262  M, &
2263  POTT, RCs2T, GRAV, &
2264  G, J33G, &
2265  RCDZ, RFDZ, RCDX, RFDX, RCDY, RFDY, FDZ, &
2266  rdt, &
2267  FACT_N, FACT_F, &
2268  TwoD, &
2269  I_XYZ, I_XYW, &
2270  IIS, IIE, JJS, JJE )
2271  implicit none
2272  real(RP), intent(inout) :: M(7,KA,IA,JA)
2273  real(RP), intent(in) :: POTT(KA,IA,JA)
2274  real(RP), intent(in) :: RCs2T(KA,IA,JA)
2275  real(RP), intent(in) :: GRAV
2276  real(RP), intent(in) :: G(KA,IA,JA,7)
2277  real(RP), intent(in) :: J33G
2278  real(RP), intent(in) :: RCDZ(KA)
2279  real(RP), intent(in) :: RFDZ(KA-1)
2280  real(RP), intent(in) :: RFDX(IA-1)
2281  real(RP), intent(in) :: RCDX(IA)
2282  real(RP), intent(in) :: RCDY(JA)
2283  real(RP), intent(in) :: RFDY(JA-1)
2284  real(RP), intent(in) :: FDZ(KA-1)
2285  real(RP), intent(in) :: rdt
2286  real(RP), intent(in) :: FACT_N
2287  real(RP), intent(in) :: FACT_F
2288  logical, intent(in) :: TwoD
2289  integer, intent(in) :: I_XYZ
2290  integer, intent(in) :: I_XYW
2291  integer, intent(in) :: IIS
2292  integer, intent(in) :: IIE
2293  integer, intent(in) :: JJS
2294  integer, intent(in) :: JJE
2295 
2296  integer :: k, i, j
2297 
2298  if ( twod ) then
2299  i = is
2300  !$omp parallel do
2301  do j = jjs, jje
2302  do k = ks+2, ke-2
2303 #ifdef DEBUG
2304  call check( __line__, pott(k-2,i,j) )
2305  call check( __line__, pott(k-1,i,j) )
2306  call check( __line__, pott(k ,i,j) )
2307  call check( __line__, pott(k+1,i,j) )
2308  call check( __line__, pott(k+2,i,j) )
2309  call check( __line__, pott(k,i ,j) )
2310  call check( __line__, pott(k,i,j-2) )
2311  call check( __line__, pott(k,i,j-1) )
2312  call check( __line__, pott(k,i,j ) )
2313  call check( __line__, pott(k,i,j+1) )
2314  call check( __line__, pott(k,i,j+2) )
2315  call check( __line__, g(k-1,i,j,i_xyw) )
2316  call check( __line__, g(k+1,i,j,i_xyw) )
2317  call check( __line__, g(k,i,j,i_xyz) )
2318  call check( __line__, rcs2t(k-1,i,j) )
2319  call check( __line__, rcs2t(k ,i,j) )
2320  call check( __line__, rcs2t(k+1,i,j) )
2321 #endif
2322  ! k,i,j
2323  m(1,k,i,j) = &
2324  - ( ( fact_n * (pott(k+1,i,j)+pott(k ,i,j)) &
2325  + fact_f * (pott(k+2,i,j)+pott(k-1,i,j)) ) * rfdz(k ) / g(k ,i,j,i_xyw) &
2326  + ( fact_n * (pott(k ,i,j)+pott(k-1,i,j)) &
2327  + fact_f * (pott(k+1,i,j)+pott(k-2,i,j)) ) * rfdz(k-1) / g(k-1,i,j,i_xyw) &
2328  ) * j33g * j33g * rfdz(k) &
2329  - ( ( fact_n * (pott(k,i,j+1)+pott(k,i,j )) &
2330  + fact_f * (pott(k,i,j+2)+pott(k,i,j-1)) ) * rfdy(j ) &
2331  + ( fact_n * (pott(k,i,j )+pott(k,i,j-1)) &
2332  + fact_f * (pott(k,i,j-1)+pott(k,i,j-2)) ) * rfdy(j-1) &
2333  ) * g(k,i,j,i_xyz) * rfdy(j) &
2334  - g(k,i,j,i_xyz) * rcs2t(k,i,j) * rdt * rdt
2335  ! k-1
2336  m(2,k,i,j) = j33g * j33g / g(k-1,i,j,i_xyw) &
2337  * ( fact_n * (pott(k ,i,j)+pott(k-1,i,j)) &
2338  + fact_f * (pott(k+1,i,j)+pott(k-2,i,j)) ) &
2339  * rfdz(k-1) * rcdz(k) &
2340  - grav * j33g * rcs2t(k-1,i,j) / ( fdz(k)+fdz(k-1) )
2341  ! k+1
2342  m(3,k,i,j) = j33g * j33g / g(k+1,i,j,i_xyw) &
2343  * ( fact_n * (pott(k+1,i,j)+pott(k ,i,j)) &
2344  + fact_f * (pott(k+2,i,j)+pott(k-1,i,j)) ) &
2345  * rfdz(k ) * rcdz(k) &
2346  + grav * j33g * rcs2t(k+1,i,j) / ( fdz(k)+fdz(k-1) )
2347  enddo
2348  enddo
2349 
2350  !$omp parallel do
2351  do j = jjs, jje
2352 #ifdef DEBUG
2353  call check( __line__, pott(ks ,i,j) )
2354  call check( __line__, pott(ks+1,i,j) )
2355  call check( __line__, pott(ks,i ,j) )
2356  call check( __line__, pott(ks,i,j-2) )
2357  call check( __line__, pott(ks,i,j-1) )
2358  call check( __line__, pott(ks,i,j ) )
2359  call check( __line__, pott(ks,i,j+1) )
2360  call check( __line__, pott(ks,i,j+2) )
2361  call check( __line__, g(ks,i,j,i_xyz) )
2362  call check( __line__, rcs2t(ks,i,j) )
2363  call check( __line__, pott(ks ,i,j) )
2364  call check( __line__, pott(ks+1,i,j) )
2365  call check( __line__, g(ks+1,i,j,i_xyw) )
2366  call check( __line__, rcs2t(ks+1,i,j) )
2367 #endif
2368  ! k,i,j
2369  m(1,ks,i,j) = &
2370  - ( 0.5_rp * (pott(ks+1,i,j)+pott(ks ,i,j)) * rfdz(ks ) &
2371  ) * j33g * j33g / g(ks ,i,j,i_xyw) * rfdz(ks) &
2372  - ( ( fact_n * (pott(ks,i,j+1)+pott(ks,i,j )) &
2373  + fact_f * (pott(ks,i,j+2)+pott(ks,i,j-1)) ) * rfdy(j ) &
2374  + ( fact_n * (pott(ks,i,j )+pott(ks,i,j-1)) &
2375  + fact_f * (pott(ks,i,j-1)+pott(ks,i,j-2)) ) * rfdy(j-1) &
2376  ) * g(ks,i,j,i_xyz) * rfdy(j) &
2377  - g(ks,i,j,i_xyz) * rcs2t(ks,i,j) * rdt * rdt &
2378  + grav * j33g * 0.5_rp * rcs2t(ks,i,j) * rcdz(ks)
2379  ! k+1
2380  m(3,ks,i,j) = j33g * j33g / g(ks+1,i,j,i_xyw) &
2381  * 0.5_rp * (pott(ks+1,i,j)+pott(ks ,i,j)) &
2382  * rfdz(ks ) * rcdz(ks) &
2383  + grav * j33g * 0.5_rp * rcs2t(ks+1,i,j) * rcdz(ks)
2384 #ifdef DEBUG
2385  call check( __line__, pott(ks ,i,j) )
2386  call check( __line__, pott(ks+1,i,j) )
2387  call check( __line__, pott(ks+2,i,j) )
2388  call check( __line__, pott(ks+3,i,j) )
2389  call check( __line__, pott(ks+1,i ,j) )
2390  call check( __line__, pott(ks+1,i,j-2) )
2391  call check( __line__, pott(ks+1,i,j-1) )
2392  call check( __line__, pott(ks+1,i,j ) )
2393  call check( __line__, pott(ks+1,i,j+1) )
2394  call check( __line__, pott(ks+1,i,j+2) )
2395  call check( __line__, g(ks ,i,j,i_xyw) )
2396  call check( __line__, g(ks+2,i,j,i_xyw) )
2397  call check( __line__, g(ks+1,i,j,i_xyz) )
2398  call check( __line__, rcs2t(ks ,i,j) )
2399  call check( __line__, rcs2t(ks+1 ,i,j) )
2400  call check( __line__, rcs2t(ks+2,i,j) )
2401 #endif
2402  ! k,i,j
2403  m(1,ks+1,i,j) = &
2404  - ( ( fact_n * (pott(ks+2,i,j)+pott(ks+1,i,j)) &
2405  + fact_f * (pott(ks+3,i,j)+pott(ks ,i,j)) ) * rfdz(ks+1) / g(ks+1,i,j,i_xyw) &
2406  + 0.5_rp * (pott(ks+1,i,j)+pott(ks ,i,j)) * rfdz(ks ) / g(ks ,i,j,i_xyw) &
2407  ) * j33g * j33g * rfdz(ks+1) &
2408  - ( ( fact_n * (pott(ks+1,i,j+1)+pott(ks+1,i,j )) &
2409  + fact_f * (pott(ks+1,i,j+2)+pott(ks+1,i,j-1)) ) * rfdy(j ) &
2410  + ( fact_n * (pott(ks+1,i,j )+pott(ks+1,i,j-1)) &
2411  + fact_f * (pott(ks+1,i,j-1)+pott(ks+1,i,j-2)) ) * rfdy(j-1) &
2412  ) * g(ks+1,i,j,i_xyz) * rfdy(j) &
2413  - g(ks+1,i,j,i_xyz) * rcs2t(ks+1,i,j) * rdt * rdt
2414  ! k-1
2415  m(2,ks+1,i,j) = j33g * j33g / g(ks,i,j,i_xyw) &
2416  * 0.5_rp * (pott(ks+1,i,j)+pott(ks ,i,j)) &
2417  * rfdz(ks ) * rcdz(ks+1) &
2418  - grav * j33g * rcs2t(ks ,i,j) / ( fdz(ks+1)+fdz(ks) )
2419  ! k+1
2420  m(3,ks+1,i,j) = j33g * j33g / g(ks+2,i,j,i_xyw) &
2421  * ( fact_n * (pott(ks+2,i,j)+pott(ks+1,i,j)) &
2422  + fact_f * (pott(ks+3,i,j)+pott(ks ,i,j)) ) &
2423  * rfdz(ks+1) * rcdz(ks+1) &
2424  + grav * j33g * rcs2t(ks+2,i,j) / ( fdz(ks+1)+fdz(ks) )
2425 #ifdef DEBUG
2426  call check( __line__, pott(ke-3,i,j) )
2427  call check( __line__, pott(ke-2,i,j) )
2428  call check( __line__, pott(ke-1,i,j) )
2429  call check( __line__, pott(ke ,i,j) )
2430  call check( __line__, pott(ke-1,i ,j) )
2431  call check( __line__, pott(ke-1,i,j-2) )
2432  call check( __line__, pott(ke-1,i,j-1) )
2433  call check( __line__, pott(ke-1,i,j ) )
2434  call check( __line__, pott(ke-1,i,j+1) )
2435  call check( __line__, pott(ke-1,i,j+2) )
2436  call check( __line__, g(ke-2,i,j,i_xyw) )
2437  call check( __line__, g(ke ,i,j,i_xyw) )
2438  call check( __line__, g(ke-1,i,j,i_xyz) )
2439  call check( __line__, rcs2t(ke-2,i,j) )
2440  call check( __line__, rcs2t(ke-1,i,j) )
2441  call check( __line__, rcs2t(ke ,i,j) )
2442 #endif
2443  ! k,i,j
2444  m(1,ke-1,i,j) = &
2445  - ( 0.5_rp * (pott(ke ,i,j)+pott(ke-1,i,j)) * rfdz(ke-1) / g(ke-1,i,j,i_xyw) &
2446  + ( fact_n * (pott(ke-1,i,j)+pott(ke-2,i,j)) &
2447  + fact_f * (pott(ke ,i,j)+pott(ke-3,i,j)) ) * rfdz(ke-2) / g(ke-2,i,j,i_xyw) &
2448  ) * j33g * j33g * rfdz(ke-1) &
2449  - ( ( fact_n * (pott(ke-1,i,j+1)+pott(ke-1,i,j )) &
2450  + fact_f * (pott(ke-1,i,j+2)+pott(ke-1,i,j-1)) ) * rfdy(j ) &
2451  + ( fact_n * (pott(ke-1,i,j )+pott(ke-1,i,j-1)) &
2452  + fact_f * (pott(ke-1,i,j-1)+pott(ke-1,i,j-2)) ) * rfdy(j-1) &
2453  ) * g(ke-1,i,j,i_xyz) * rfdy(j) &
2454  - g(ke-1,i,j,i_xyz) * rcs2t(ke-1,i,j) * rdt * rdt
2455  ! k-1
2456  m(2,ke-1,i,j) = j33g * j33g / g(ke-2,i,j,i_xyw) &
2457  * ( fact_n * (pott(ke-1,i,j)+pott(ke-2,i,j)) &
2458  + fact_f * (pott(ke ,i,j)+pott(ke-3,i,j)) ) &
2459  * rfdz(ke-2) * rcdz(ke-1) &
2460  - grav * j33g * rcs2t(ke-2,i,j) / ( fdz(ke-1)+fdz(ke-2) )
2461  ! k+1
2462  m(3,ke-1,i,j) = j33g * j33g / g(ke ,i,j,i_xyw) &
2463  * 0.5_rp * (pott(ke ,i,j)+pott(ke-1,i,j)) &
2464  * rfdz(ke-1) * rcdz(ke-1) &
2465  + grav * j33g * rcs2t(ke ,i,j) / ( fdz(ke-1)+fdz(ke-2) )
2466 #ifdef DEBUG
2467  call check( __line__, pott(ke-1,i,j) )
2468  call check( __line__, pott(ke ,i,j) )
2469  call check( __line__, pott(ke,i ,j) )
2470  call check( __line__, pott(ke,i,j-2) )
2471  call check( __line__, pott(ke,i,j-1) )
2472  call check( __line__, pott(ke,i,j ) )
2473  call check( __line__, pott(ke,i,j+1) )
2474  call check( __line__, pott(ke,i,j+2) )
2475  call check( __line__, g(ke-1,i,j,i_xyw) )
2476  call check( __line__, g(ke,i,j,i_xyz) )
2477  call check( __line__, rcs2t(ke-1,i,j) )
2478  call check( __line__, rcs2t(ke,i,j) )
2479 #endif
2480  ! k,i,j
2481  m(1,ke,i,j) = &
2482  - ( &
2483  + 0.5_rp * (pott(ke ,i,j)+pott(ke-1,i,j)) * rfdz(ke-1) / g(ke-1,i,j,i_xyw) &
2484  ) * j33g * j33g * rfdz(ke) &
2485  - ( ( fact_n * (pott(ke,i,j+1)+pott(ke,i,j )) &
2486  + fact_f * (pott(ke,i,j+2)+pott(ke,i,j-1)) ) * rfdy(j ) &
2487  + ( fact_n * (pott(ke,i,j )+pott(ke,i,j-1)) &
2488  + fact_f * (pott(ke,i,j-1)+pott(ke,i,j-2)) ) * rfdy(j-1) &
2489  ) * g(ke,i,j,i_xyz) * rfdy(j) &
2490  - g(ke,i,j,i_xyz) * rcs2t(ke,i,j) * rdt * rdt &
2491  - grav * j33g * 0.5_rp * rcs2t(ke,i,j) * rcdz(ke)
2492  ! k-1
2493  m(2,ke,i,j) = j33g * j33g / g(ke-1,i,j,i_xyw) &
2494  * 0.5_rp * (pott(ke ,i,j)+pott(ke-1,i,j)) &
2495  * rfdz(ke-1) * rcdz(ke) &
2496  - grav * j33g * 0.5_rp * rcs2t(ke,i,j) * rcdz(ke)
2497  enddo
2498 
2499  else
2500  !$omp parallel do
2501  do j = jjs, jje
2502  do i = iis, iie
2503  do k = ks+2, ke-2
2504 #ifdef DEBUG
2505  call check( __line__, pott(k-2,i,j) )
2506  call check( __line__, pott(k-1,i,j) )
2507  call check( __line__, pott(k ,i,j) )
2508  call check( __line__, pott(k+1,i,j) )
2509  call check( __line__, pott(k+2,i,j) )
2510  call check( __line__, pott(k,i-2,j) )
2511  call check( __line__, pott(k,i-1,j) )
2512  call check( __line__, pott(k,i ,j) )
2513  call check( __line__, pott(k,i+1,j) )
2514  call check( __line__, pott(k,i+2,j) )
2515  call check( __line__, pott(k,i,j-2) )
2516  call check( __line__, pott(k,i,j-1) )
2517  call check( __line__, pott(k,i,j ) )
2518  call check( __line__, pott(k,i,j+1) )
2519  call check( __line__, pott(k,i,j+2) )
2520  call check( __line__, g(k-1,i,j,i_xyw) )
2521  call check( __line__, g(k+1,i,j,i_xyw) )
2522  call check( __line__, g(k,i,j,i_xyz) )
2523  call check( __line__, rcs2t(k-1,i,j) )
2524  call check( __line__, rcs2t(k ,i,j) )
2525  call check( __line__, rcs2t(k+1,i,j) )
2526 #endif
2527  ! k,i,j
2528  m(1,k,i,j) = &
2529  - ( ( fact_n * (pott(k+1,i,j)+pott(k ,i,j)) &
2530  + fact_f * (pott(k+2,i,j)+pott(k-1,i,j)) ) * rfdz(k ) / g(k ,i,j,i_xyw) &
2531  + ( fact_n * (pott(k ,i,j)+pott(k-1,i,j)) &
2532  + fact_f * (pott(k+1,i,j)+pott(k-2,i,j)) ) * rfdz(k-1) / g(k-1,i,j,i_xyw) &
2533  ) * j33g * j33g * rfdz(k) &
2534  - ( ( fact_n * (pott(k,i+1,j)+pott(k,i ,j)) &
2535  + fact_f * (pott(k,i+2,j)+pott(k,i-1,j)) ) * rfdx(i ) &
2536  + ( fact_n * (pott(k,i ,j)+pott(k,i-1,j)) &
2537  + fact_f * (pott(k,i+1,j)+pott(k,i-2,j)) ) * rfdx(i-1) &
2538  ) * g(k,i,j,i_xyz) * rfdx(i) &
2539  - ( ( fact_n * (pott(k,i,j+1)+pott(k,i,j )) &
2540  + fact_f * (pott(k,i,j+2)+pott(k,i,j-1)) ) * rfdy(j ) &
2541  + ( fact_n * (pott(k,i,j )+pott(k,i,j-1)) &
2542  + fact_f * (pott(k,i,j-1)+pott(k,i,j-2)) ) * rfdy(j-1) &
2543  ) * g(k,i,j,i_xyz) * rfdy(j) &
2544  - g(k,i,j,i_xyz) * rcs2t(k,i,j) * rdt * rdt
2545  ! k-1
2546  m(2,k,i,j) = j33g * j33g / g(k-1,i,j,i_xyw) &
2547  * ( fact_n * (pott(k ,i,j)+pott(k-1,i,j)) &
2548  + fact_f * (pott(k+1,i,j)+pott(k-2,i,j)) ) &
2549  * rfdz(k-1) * rcdz(k) &
2550  - grav * j33g * rcs2t(k-1,i,j) / ( fdz(k)+fdz(k-1) )
2551  ! k+1
2552  m(3,k,i,j) = j33g * j33g / g(k+1,i,j,i_xyw) &
2553  * ( fact_n * (pott(k+1,i,j)+pott(k ,i,j)) &
2554  + fact_f * (pott(k+2,i,j)+pott(k-1,i,j)) ) &
2555  * rfdz(k ) * rcdz(k) &
2556  + grav * j33g * rcs2t(k+1,i,j) / ( fdz(k)+fdz(k-1) )
2557  enddo
2558  enddo
2559  enddo
2560 
2561  !$omp parallel do
2562  do j = jjs, jje
2563  do i = iis, iie
2564 #ifdef DEBUG
2565  call check( __line__, pott(ks ,i,j) )
2566  call check( __line__, pott(ks+1,i,j) )
2567  call check( __line__, pott(ks,i-2,j) )
2568  call check( __line__, pott(ks,i-1,j) )
2569  call check( __line__, pott(ks,i ,j) )
2570  call check( __line__, pott(ks,i+1,j) )
2571  call check( __line__, pott(ks,i+2,j) )
2572  call check( __line__, pott(ks,i,j-2) )
2573  call check( __line__, pott(ks,i,j-1) )
2574  call check( __line__, pott(ks,i,j ) )
2575  call check( __line__, pott(ks,i,j+1) )
2576  call check( __line__, pott(ks,i,j+2) )
2577  call check( __line__, g(ks,i,j,i_xyz) )
2578  call check( __line__, rcs2t(ks,i,j) )
2579  call check( __line__, pott(ks ,i,j) )
2580  call check( __line__, pott(ks+1,i,j) )
2581  call check( __line__, g(ks+1,i,j,i_xyw) )
2582  call check( __line__, rcs2t(ks+1,i,j) )
2583 #endif
2584  ! k,i,j
2585  m(1,ks,i,j) = &
2586  - ( 0.5_rp * (pott(ks+1,i,j)+pott(ks ,i,j)) * rfdz(ks ) &
2587  ) * j33g * j33g / g(ks ,i,j,i_xyw) * rfdz(ks) &
2588  - ( ( fact_n * (pott(ks,i+1,j)+pott(ks,i ,j)) &
2589  + fact_f * (pott(ks,i+2,j)+pott(ks,i-1,j)) ) * rfdx(i ) &
2590  + ( fact_n * (pott(ks,i ,j)+pott(ks,i-1,j)) &
2591  + fact_f * (pott(ks,i+1,j)+pott(ks,i-2,j)) ) * rfdx(i-1) &
2592  ) * g(ks,i,j,i_xyz) * rfdx(i) &
2593  - ( ( fact_n * (pott(ks,i,j+1)+pott(ks,i,j )) &
2594  + fact_f * (pott(ks,i,j+2)+pott(ks,i,j-1)) ) * rfdy(j ) &
2595  + ( fact_n * (pott(ks,i,j )+pott(ks,i,j-1)) &
2596  + fact_f * (pott(ks,i,j-1)+pott(ks,i,j-2)) ) * rfdy(j-1) &
2597  ) * g(ks,i,j,i_xyz) * rfdy(j) &
2598  - g(ks,i,j,i_xyz) * rcs2t(ks,i,j) * rdt * rdt &
2599  + grav * j33g * 0.5_rp * rcs2t(ks,i,j) * rcdz(ks)
2600  ! k+1
2601  m(3,ks,i,j) = j33g * j33g / g(ks+1,i,j,i_xyw) &
2602  * 0.5_rp * (pott(ks+1,i,j)+pott(ks ,i,j)) &
2603  * rfdz(ks ) * rcdz(ks) &
2604  + grav * j33g * 0.5_rp * rcs2t(ks+1,i,j) * rcdz(ks)
2605 #ifdef DEBUG
2606  call check( __line__, pott(ks ,i,j) )
2607  call check( __line__, pott(ks+1,i,j) )
2608  call check( __line__, pott(ks+2,i,j) )
2609  call check( __line__, pott(ks+3,i,j) )
2610  call check( __line__, pott(ks+1,i-2,j) )
2611  call check( __line__, pott(ks+1,i-1,j) )
2612  call check( __line__, pott(ks+1,i ,j) )
2613  call check( __line__, pott(ks+1,i+1,j) )
2614  call check( __line__, pott(ks+1,i+2,j) )
2615  call check( __line__, pott(ks+1,i,j-2) )
2616  call check( __line__, pott(ks+1,i,j-1) )
2617  call check( __line__, pott(ks+1,i,j ) )
2618  call check( __line__, pott(ks+1,i,j+1) )
2619  call check( __line__, pott(ks+1,i,j+2) )
2620  call check( __line__, g(ks ,i,j,i_xyw) )
2621  call check( __line__, g(ks+2,i,j,i_xyw) )
2622  call check( __line__, g(ks+1,i,j,i_xyz) )
2623  call check( __line__, rcs2t(ks ,i,j) )
2624  call check( __line__, rcs2t(ks+1 ,i,j) )
2625  call check( __line__, rcs2t(ks+2,i,j) )
2626 #endif
2627  ! k,i,j
2628  m(1,ks+1,i,j) = &
2629  - ( ( fact_n * (pott(ks+2,i,j)+pott(ks+1,i,j)) &
2630  + fact_f * (pott(ks+3,i,j)+pott(ks ,i,j)) ) * rfdz(ks+1) / g(ks+1,i,j,i_xyw) &
2631  + 0.5_rp * (pott(ks+1,i,j)+pott(ks ,i,j)) * rfdz(ks ) / g(ks ,i,j,i_xyw) &
2632  ) * j33g * j33g * rfdz(ks+1) &
2633  - ( ( fact_n * (pott(ks+1,i+1,j)+pott(ks+1,i ,j)) &
2634  + fact_f * (pott(ks+1,i+2,j)+pott(ks+1,i-1,j)) ) * rfdx(i ) &
2635  + ( fact_n * (pott(ks+1,i ,j)+pott(ks+1,i-1,j)) &
2636  + fact_f * (pott(ks+1,i+1,j)+pott(ks+1,i-2,j)) ) * rfdx(i-1) &
2637  ) * g(ks+1,i,j,i_xyz) * rfdx(i) &
2638  - ( ( fact_n * (pott(ks+1,i,j+1)+pott(ks+1,i,j )) &
2639  + fact_f * (pott(ks+1,i,j+2)+pott(ks+1,i,j-1)) ) * rfdy(j ) &
2640  + ( fact_n * (pott(ks+1,i,j )+pott(ks+1,i,j-1)) &
2641  + fact_f * (pott(ks+1,i,j-1)+pott(ks+1,i,j-2)) ) * rfdy(j-1) &
2642  ) * g(ks+1,i,j,i_xyz) * rfdy(j) &
2643  - g(ks+1,i,j,i_xyz) * rcs2t(ks+1,i,j) * rdt * rdt
2644  ! k-1
2645  m(2,ks+1,i,j) = j33g * j33g / g(ks,i,j,i_xyw) &
2646  * 0.5_rp * (pott(ks+1,i,j)+pott(ks ,i,j)) &
2647  * rfdz(ks ) * rcdz(ks+1) &
2648  - grav * j33g * rcs2t(ks ,i,j) / ( fdz(ks+1)+fdz(ks) )
2649  ! k+1
2650  m(3,ks+1,i,j) = j33g * j33g / g(ks+2,i,j,i_xyw) &
2651  * ( fact_n * (pott(ks+2,i,j)+pott(ks+1,i,j)) &
2652  + fact_f * (pott(ks+3,i,j)+pott(ks ,i,j)) ) &
2653  * rfdz(ks+1) * rcdz(ks+1) &
2654  + grav * j33g * rcs2t(ks+2,i,j) / ( fdz(ks+1)+fdz(ks) )
2655 #ifdef DEBUG
2656  call check( __line__, pott(ke-3,i,j) )
2657  call check( __line__, pott(ke-2,i,j) )
2658  call check( __line__, pott(ke-1,i,j) )
2659  call check( __line__, pott(ke ,i,j) )
2660  call check( __line__, pott(ke-1,i-2,j) )
2661  call check( __line__, pott(ke-1,i-1,j) )
2662  call check( __line__, pott(ke-1,i ,j) )
2663  call check( __line__, pott(ke-1,i+1,j) )
2664  call check( __line__, pott(ke-1,i+2,j) )
2665  call check( __line__, pott(ke-1,i,j-2) )
2666  call check( __line__, pott(ke-1,i,j-1) )
2667  call check( __line__, pott(ke-1,i,j ) )
2668  call check( __line__, pott(ke-1,i,j+1) )
2669  call check( __line__, pott(ke-1,i,j+2) )
2670  call check( __line__, g(ke-2,i,j,i_xyw) )
2671  call check( __line__, g(ke ,i,j,i_xyw) )
2672  call check( __line__, g(ke-1,i,j,i_xyz) )
2673  call check( __line__, rcs2t(ke-2,i,j) )
2674  call check( __line__, rcs2t(ke-1,i,j) )
2675  call check( __line__, rcs2t(ke ,i,j) )
2676 #endif
2677  ! k,i,j
2678  m(1,ke-1,i,j) = &
2679  - ( 0.5_rp * (pott(ke ,i,j)+pott(ke-1,i,j)) * rfdz(ke-1) / g(ke-1,i,j,i_xyw) &
2680  + ( fact_n * (pott(ke-1,i,j)+pott(ke-2,i,j)) &
2681  + fact_f * (pott(ke ,i,j)+pott(ke-3,i,j)) ) * rfdz(ke-2) / g(ke-2,i,j,i_xyw) &
2682  ) * j33g * j33g * rfdz(ke-1) &
2683  - ( ( fact_n * (pott(ke-1,i+1,j)+pott(ke-1,i ,j)) &
2684  + fact_f * (pott(ke-1,i+2,j)+pott(ke-1,i-1,j)) ) * rfdx(i ) &
2685  + ( fact_n * (pott(ke-1,i ,j)+pott(ke-1,i-1,j)) &
2686  + fact_f * (pott(ke-1,i+1,j)+pott(ke-1,i-2,j)) ) * rfdx(i-1) &
2687  ) * g(ke-1,i,j,i_xyz) * rfdx(i) &
2688  - ( ( fact_n * (pott(ke-1,i,j+1)+pott(ke-1,i,j )) &
2689  + fact_f * (pott(ke-1,i,j+2)+pott(ke-1,i,j-1)) ) * rfdy(j ) &
2690  + ( fact_n * (pott(ke-1,i,j )+pott(ke-1,i,j-1)) &
2691  + fact_f * (pott(ke-1,i,j-1)+pott(ke-1,i,j-2)) ) * rfdy(j-1) &
2692  ) * g(ke-1,i,j,i_xyz) * rfdy(j) &
2693  - g(ke-1,i,j,i_xyz) * rcs2t(ke-1,i,j) * rdt * rdt
2694  ! k-1
2695  m(2,ke-1,i,j) = j33g * j33g / g(ke-2,i,j,i_xyw) &
2696  * ( fact_n * (pott(ke-1,i,j)+pott(ke-2,i,j)) &
2697  + fact_f * (pott(ke ,i,j)+pott(ke-3,i,j)) ) &
2698  * rfdz(ke-2) * rcdz(ke-1) &
2699  - grav * j33g * rcs2t(ke-2,i,j) / ( fdz(ke-1)+fdz(ke-2) )
2700  ! k+1
2701  m(3,ke-1,i,j) = j33g * j33g / g(ke ,i,j,i_xyw) &
2702  * 0.5_rp * (pott(ke ,i,j)+pott(ke-1,i,j)) &
2703  * rfdz(ke-1) * rcdz(ke-1) &
2704  + grav * j33g * rcs2t(ke ,i,j) / ( fdz(ke-1)+fdz(ke-2) )
2705 #ifdef DEBUG
2706  call check( __line__, pott(ke-1,i,j) )
2707  call check( __line__, pott(ke ,i,j) )
2708  call check( __line__, pott(ke,i-2,j) )
2709  call check( __line__, pott(ke,i-1,j) )
2710  call check( __line__, pott(ke,i ,j) )
2711  call check( __line__, pott(ke,i+1,j) )
2712  call check( __line__, pott(ke,i+2,j) )
2713  call check( __line__, pott(ke,i,j-2) )
2714  call check( __line__, pott(ke,i,j-1) )
2715  call check( __line__, pott(ke,i,j ) )
2716  call check( __line__, pott(ke,i,j+1) )
2717  call check( __line__, pott(ke,i,j+2) )
2718  call check( __line__, g(ke-1,i,j,i_xyw) )
2719  call check( __line__, g(ke,i,j,i_xyz) )
2720  call check( __line__, rcs2t(ke-1,i,j) )
2721  call check( __line__, rcs2t(ke,i,j) )
2722 #endif
2723  ! k,i,j
2724  m(1,ke,i,j) = &
2725  - ( &
2726  + 0.5_rp * (pott(ke ,i,j)+pott(ke-1,i,j)) * rfdz(ke-1) / g(ke-1,i,j,i_xyw) &
2727  ) * j33g * j33g * rfdz(ke) &
2728  - ( ( fact_n * (pott(ke,i+1,j)+pott(ke,i ,j)) &
2729  + fact_f * (pott(ke,i+2,j)+pott(ke,i-1,j)) ) * rfdx(i ) &
2730  + ( fact_n * (pott(ke,i ,j)+pott(ke,i-1,j)) &
2731  + fact_f * (pott(ke,i+1,j)+pott(ke,i-2,j)) ) * rfdx(i-1) &
2732  ) * g(ke,i,j,i_xyz) * rfdx(i) &
2733  - ( ( fact_n * (pott(ke,i,j+1)+pott(ke,i,j )) &
2734  + fact_f * (pott(ke,i,j+2)+pott(ke,i,j-1)) ) * rfdy(j ) &
2735  + ( fact_n * (pott(ke,i,j )+pott(ke,i,j-1)) &
2736  + fact_f * (pott(ke,i,j-1)+pott(ke,i,j-2)) ) * rfdy(j-1) &
2737  ) * g(ke,i,j,i_xyz) * rfdy(j) &
2738  - g(ke,i,j,i_xyz) * rcs2t(ke,i,j) * rdt * rdt &
2739  - grav * j33g * 0.5_rp * rcs2t(ke,i,j) * rcdz(ke)
2740  ! k-1
2741  m(2,ke,i,j) = j33g * j33g / g(ke-1,i,j,i_xyw) &
2742  * 0.5_rp * (pott(ke ,i,j)+pott(ke-1,i,j)) &
2743  * rfdz(ke-1) * rcdz(ke) &
2744  - grav * j33g * 0.5_rp * rcs2t(ke,i,j) * rcdz(ke)
2745  enddo
2746  enddo
2747 
2748  end if
2749 #ifdef DEBUG
2750  k = iundef; i = iundef; j = iundef
2751 #endif
2752 
2753  if ( twod ) then
2754  i = is
2755  !$omp parallel do
2756  do j = jjs, jje
2757  do k = ks, ke
2758 #ifdef DEBUG
2759  call check( __line__, g(k,i,j-1,i_xyz) )
2760  call check( __line__, g(k,i,j+1,i_xyz) )
2761  call check( __line__, pott(k,i ,j ) )
2762  call check( __line__, pott(k,i ,j-2) )
2763  call check( __line__, pott(k,i ,j-1) )
2764  call check( __line__, pott(k,i ,j ) )
2765  call check( __line__, pott(k,i ,j+1) )
2766  call check( __line__, pott(k,i ,j+2) )
2767 #endif
2768  ! i-1
2769  !M(4,k,i,j) = 0.0_RP
2770  ! i+1
2771  !M(5,k,i,j) = 0.0_RP
2772  ! j-1
2773  m(6,k,i,j) = g(k,i,j-1,i_xyz) &
2774  * ( fact_n * (pott(k,i,j )+pott(k,i,j-1)) &
2775  + fact_f * (pott(k,i,j+1)+pott(k,i,j-2)) ) &
2776  * rfdy(j-1) * rcdy(j)
2777  ! j+1
2778  m(7,k,i,j) = g(k,i,j+1,i_xyz) &
2779  * ( fact_n * (pott(k,i,j+1)+pott(k,i,j )) &
2780  + fact_f * (pott(k,i,j+2)+pott(k,i,j-1)) ) &
2781  * rfdy(j ) * rcdy(j)
2782  enddo
2783  enddo
2784  else
2785  !$omp parallel do
2786  do j = jjs, jje
2787  do i = iis, iie
2788  do k = ks, ke
2789 #ifdef DEBUG
2790  call check( __line__, g(k,i-1,j,i_xyz) )
2791  call check( __line__, g(k,i+1,j,i_xyz) )
2792  call check( __line__, g(k,i,j-1,i_xyz) )
2793  call check( __line__, g(k,i,j+1,i_xyz) )
2794  call check( __line__, pott(k,i-2,j ) )
2795  call check( __line__, pott(k,i-1,j ) )
2796  call check( __line__, pott(k,i ,j ) )
2797  call check( __line__, pott(k,i+1,j ) )
2798  call check( __line__, pott(k,i+2,j ) )
2799  call check( __line__, pott(k,i ,j-2) )
2800  call check( __line__, pott(k,i ,j-1) )
2801  call check( __line__, pott(k,i ,j ) )
2802  call check( __line__, pott(k,i ,j+1) )
2803  call check( __line__, pott(k,i ,j+2) )
2804 #endif
2805  ! i-1
2806  m(4,k,i,j) = g(k,i-1,j,i_xyz) &
2807  * ( fact_n * (pott(k,i ,j)+pott(k,i-1,j)) &
2808  + fact_f * (pott(k,i+1,j)+pott(k,i-2,j)) ) &
2809  * rfdx(i-1) * rcdx(i)
2810  ! i+1
2811  m(5,k,i,j) = g(k,i+1,j,i_xyz) &
2812  * ( fact_n * (pott(k,i+1,j)+pott(k,i ,j)) &
2813  + fact_f * (pott(k,i+2,j)+pott(k,i-1,j)) ) &
2814  * rfdx(i ) * rcdx(i)
2815  ! j-1
2816  m(6,k,i,j) = g(k,i,j-1,i_xyz) &
2817  * ( fact_n * (pott(k,i,j )+pott(k,i,j-1)) &
2818  + fact_f * (pott(k,i,j+1)+pott(k,i,j-2)) ) &
2819  * rfdy(j-1) * rcdy(j)
2820  ! j+1
2821  m(7,k,i,j) = g(k,i,j+1,i_xyz) &
2822  * ( fact_n * (pott(k,i,j+1)+pott(k,i,j )) &
2823  + fact_f * (pott(k,i,j+2)+pott(k,i,j-1)) ) &
2824  * rfdy(j ) * rcdy(j)
2825  enddo
2826  enddo
2827  enddo
2828  end if
2829 #ifdef DEBUG
2830  k = iundef; i = iundef; j = iundef
2831 #endif
2832 
2833  return
2834  end subroutine make_matrix
2835 
2836  subroutine mul_matrix(V, M, C, TwoD)
2837  implicit none
2838  real(RP), intent(out) :: V(KA,IA,JA)
2839  real(RP), intent(in) :: M(7,KA,IA,JA)
2840  real(RP), intent(in) :: C(KA,IA,JA)
2841  logical, intent(in) :: TwoD
2842 
2843  integer :: k, i, j
2844 
2845 #ifdef DEBUG
2846  v(:,:,:) = undef
2847 #endif
2848 
2849  if ( twod ) then
2850  i = is
2851  !$omp parallel do
2852  do j = js, je
2853  do k = ks+1, ke-1
2854 #ifdef DEBUG
2855  call check( __line__, m(1,k,i,j) )
2856  call check( __line__, m(2,k,i,j) )
2857  call check( __line__, m(3,k,i,j) )
2858  call check( __line__, m(6,k,i,j) )
2859  call check( __line__, m(7,k,i,j) )
2860  call check( __line__, c(k ,i ,j ) )
2861  call check( __line__, c(k-1,i ,j ) )
2862  call check( __line__, c(k+1,i ,j ) )
2863  call check( __line__, c(k ,i ,j-1) )
2864  call check( __line__, c(k ,i ,j+1) )
2865 #endif
2866  v(k,i,j) = m(1,k,i,j) * c(k ,i ,j ) &
2867  + m(2,k,i,j) * c(k-1,i ,j ) &
2868  + m(3,k,i,j) * c(k+1,i ,j ) &
2869  + m(6,k,i,j) * c(k ,i ,j-1) &
2870  + m(7,k,i,j) * c(k ,i ,j+1)
2871  enddo
2872  enddo
2873  !$omp parallel do
2874  do j = js, je
2875  do i = is, ie
2876 #ifdef DEBUG
2877  call check( __line__, m(1,ks,i,j) )
2878  call check( __line__, m(3,ks,i,j) )
2879  call check( __line__, m(6,ks,i,j) )
2880  call check( __line__, m(7,ks,i,j) )
2881  call check( __line__, c(ks ,i ,j ) )
2882  call check( __line__, c(ks+1,i ,j ) )
2883  call check( __line__, c(ks ,i ,j-1) )
2884  call check( __line__, c(ks ,i ,j+1) )
2885  call check( __line__, m(1,ke,i,j) )
2886  call check( __line__, m(2,ke,i,j) )
2887  call check( __line__, m(6,ke,i,j) )
2888  call check( __line__, m(7,ke,i,j) )
2889  call check( __line__, c(ke ,i ,j ) )
2890  call check( __line__, c(ke-1,i ,j ) )
2891  call check( __line__, c(ke ,i ,j-1) )
2892  call check( __line__, c(ke ,i ,j+1) )
2893 #endif
2894  v(ks,i,j) = m(1,ks,i,j) * c(ks ,i ,j ) &
2895  + m(3,ks,i,j) * c(ks+1,i ,j ) &
2896  + m(6,ks,i,j) * c(ks ,i ,j-1) &
2897  + m(7,ks,i,j) * c(ks ,i ,j+1)
2898  v(ke,i,j) = m(1,ke,i,j) * c(ke ,i ,j ) &
2899  + m(2,ke,i,j) * c(ke-1,i ,j ) &
2900  + m(6,ke,i,j) * c(ke ,i ,j-1) &
2901  + m(7,ke,i,j) * c(ke ,i ,j+1)
2902  enddo
2903  enddo
2904 
2905  else
2906  !$omp parallel do
2907  do j = js, je
2908  do i = is, ie
2909  do k = ks+1, ke-1
2910 #ifdef DEBUG
2911  call check( __line__, m(1,k,i,j) )
2912  call check( __line__, m(2,k,i,j) )
2913  call check( __line__, m(3,k,i,j) )
2914  call check( __line__, m(4,k,i,j) )
2915  call check( __line__, m(5,k,i,j) )
2916  call check( __line__, m(6,k,i,j) )
2917  call check( __line__, m(7,k,i,j) )
2918  call check( __line__, c(k ,i ,j ) )
2919  call check( __line__, c(k-1,i ,j ) )
2920  call check( __line__, c(k+1,i ,j ) )
2921  call check( __line__, c(k ,i-1,j ) )
2922  call check( __line__, c(k ,i+1,j ) )
2923  call check( __line__, c(k ,i ,j-1) )
2924  call check( __line__, c(k ,i ,j+1) )
2925 #endif
2926  v(k,i,j) = m(1,k,i,j) * c(k ,i ,j ) &
2927  + m(2,k,i,j) * c(k-1,i ,j ) &
2928  + m(3,k,i,j) * c(k+1,i ,j ) &
2929  + m(4,k,i,j) * c(k ,i-1,j ) &
2930  + m(5,k,i,j) * c(k ,i+1,j ) &
2931  + m(6,k,i,j) * c(k ,i ,j-1) &
2932  + m(7,k,i,j) * c(k ,i ,j+1)
2933  enddo
2934  enddo
2935  enddo
2936 
2937  !$omp parallel do
2938  do j = js, je
2939  do i = is, ie
2940 #ifdef DEBUG
2941  call check( __line__, m(1,ks,i,j) )
2942  call check( __line__, m(3,ks,i,j) )
2943  call check( __line__, m(4,ks,i,j) )
2944  call check( __line__, m(5,ks,i,j) )
2945  call check( __line__, m(6,ks,i,j) )
2946  call check( __line__, m(7,ks,i,j) )
2947  call check( __line__, c(ks ,i ,j ) )
2948  call check( __line__, c(ks+1,i ,j ) )
2949  call check( __line__, c(ks ,i-1,j ) )
2950  call check( __line__, c(ks ,i+1,j ) )
2951  call check( __line__, c(ks ,i ,j-1) )
2952  call check( __line__, c(ks ,i ,j+1) )
2953  call check( __line__, m(1,ke,i,j) )
2954  call check( __line__, m(2,ke,i,j) )
2955  call check( __line__, m(4,ke,i,j) )
2956  call check( __line__, m(5,ke,i,j) )
2957  call check( __line__, m(6,ke,i,j) )
2958  call check( __line__, m(7,ke,i,j) )
2959  call check( __line__, c(ke ,i ,j ) )
2960  call check( __line__, c(ke-1,i ,j ) )
2961  call check( __line__, c(ke ,i-1,j ) )
2962  call check( __line__, c(ke ,i+1,j ) )
2963  call check( __line__, c(ke ,i ,j-1) )
2964  call check( __line__, c(ke ,i ,j+1) )
2965 #endif
2966  v(ks,i,j) = m(1,ks,i,j) * c(ks ,i ,j ) &
2967  + m(3,ks,i,j) * c(ks+1,i ,j ) &
2968  + m(4,ks,i,j) * c(ks ,i-1,j ) &
2969  + m(5,ks,i,j) * c(ks ,i+1,j ) &
2970  + m(6,ks,i,j) * c(ks ,i ,j-1) &
2971  + m(7,ks,i,j) * c(ks ,i ,j+1)
2972  v(ke,i,j) = m(1,ke,i,j) * c(ke ,i ,j ) &
2973  + m(2,ke,i,j) * c(ke-1,i ,j ) &
2974  + m(4,ke,i,j) * c(ke ,i-1,j ) &
2975  + m(5,ke,i,j) * c(ke ,i+1,j ) &
2976  + m(6,ke,i,j) * c(ke ,i ,j-1) &
2977  + m(7,ke,i,j) * c(ke ,i ,j+1)
2978  enddo
2979  enddo
2980 
2981  end if
2982 #ifdef DEBUG
2983  k = iundef; i = iundef; j = iundef
2984 #endif
2985 
2986  return
2987  end subroutine mul_matrix
2988 
2989 #else
2990  subroutine solve_multigrid
2991  end subroutine solve_multigrid
2992 #endif
2993 
2994 #ifdef DEBUG
2995  subroutine check_solver( &
2996  DPRES, M, B, TwoD )
2997  use scale_prc, only: &
2998  prc_abort
2999  real(RP), intent(in) :: DPRES(KA,IA,JA)
3000  real(RP), intent(in) :: M(7,KA,IA,JA)
3001  real(RP), intent(in) :: B(KA,IA,JA)
3002  logical, intent(in) :: TwoD
3003 
3004  real(RP) :: B2(KA,IA,JA)
3005  integer :: k, i, j
3006  real(RP) :: err
3007 
3008  call mul_matrix(b2, m, dpres, twod)
3009 
3010  do j = js, je
3011  do i = is, ie
3012  do k = ks, ke
3013  err = abs( b2(k,i,j) - b(k,i,j) )
3014  if ( err > 1.e-5_rp .and. abs( err / b(k,i,j) ) > 1.e-5_rp ) then
3015  log_error("check_solver",*) "solver error is too large: ", k,i,j, b(k,i,j), b2(k,i,j)
3016  call prc_abort
3017  endif
3018  enddo
3019  enddo
3020  enddo
3021 
3022  end subroutine check_solver
3023 
3024  subroutine check_pres( &
3025  DPRES_N, DPRES, &
3026  RHOT_RK, RHOT, &
3027  DENS_RK, DENS, &
3028  B, &
3029  RT2P )
3030  use scale_prc, only: &
3031  prc_abort
3032  real(RP), intent(in) :: DPRES_N(KA,IA,JA)
3033  real(RP), intent(in) :: DPRES(KA,IA,JA)
3034  real(RP), intent(in) :: RHOT_RK(KA,IA,JA)
3035  real(RP), intent(in) :: RHOT(KA,IA,JA)
3036  real(RP), intent(in) :: DENS_RK(KA,IA,JA)
3037  real(RP), intent(in) :: DENS(KA,IA,JA)
3038  real(RP), intent(in) :: B(KA,IA,JA)
3039  real(RP), intent(in) :: RT2P(KA,IA,JA)
3040 
3041  real(RP) :: lhs, rhs
3042  integer :: k,i,j
3043 
3044  do j = js, je
3045  do i = is, ie
3046  do k = ks, ke
3047  lhs = dpres_n(k,i,j) - dpres(k,i,j)
3048  rhs = rt2p(k,i,j) * ( rhot_rk(k,i,j) - rhot(k,i,j) )
3049  if ( abs( (lhs - rhs) / lhs ) > 1e-15 ) then
3050  log_error("check_pres",*) "error is too large: ", k,i,j, lhs, rhs, &
3051  dpres_n(k,i,j),dpres(k,i,j),rhot_rk(k,i,j),rhot(k,i,j),dens_rk(k,i,j),dens(k,i,j),b(k,i,j)
3052  call prc_abort
3053  endif
3054  enddo
3055  enddo
3056  enddo
3057  end subroutine check_pres
3058 #endif
3059 
3060 end module scale_atmos_dyn_tstep_short_fvm_hivi
scale_const::const_grav
real(rp), public const_grav
standard acceleration of gravity [m/s2]
Definition: scale_const.F90:46
scale_precision::sp
integer, parameter, public sp
Definition: scale_precision.F90:31
scale_atmos_grid_cartesc_index::i_uy
integer, public i_uy
Definition: scale_atmos_grid_cartesC_index.F90:99
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::i_xv
integer, public i_xv
Definition: scale_atmos_grid_cartesC_index.F90:100
scale_prc::prc_abort
subroutine, public prc_abort
Abort Process.
Definition: scale_prc.F90:342
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:212
scale_index::i_momz
integer, parameter, public i_momz
Definition: scale_index.F90:29
scale_atmos_dyn_tstep_short_fvm_hivi::make_matrix
subroutine make_matrix(M, POTT, RCs2T, GRAV, G, J33G, RCDZ, RFDZ, RCDX, RFDX, RCDY, RFDY, FDZ, rdt, FACT_N, FACT_F, TwoD, I_XYZ, I_XYW, IIS, IIE, JJS, JJE)
Definition: scale_atmos_dyn_tstep_short_fvm_hivi.F90:2271
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:173
scale_atmos_dyn_tstep_short_fvm_hivi::solve_bicgstab
subroutine solve_bicgstab(DPRES_N, DPRES, M, B, TwoD)
Definition: scale_atmos_dyn_tstep_short_fvm_hivi.F90:1943
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:201
scale_atmos_dyn_fvm_fct::atmos_dyn_fvm_fct
subroutine, public atmos_dyn_fvm_fct(qflx_anti, phi_in, DENS0, DENS, qflx_hi, qflx_lo, mflx_hi, rdz, rdx, rdy, GSQRT, MAPF, TwoD, dt, flag_vect)
Setup.
Definition: scale_atmos_dyn_fvm_fct.F90:72
scale_index::i_momx
integer, parameter, public i_momx
Definition: scale_index.F90:30
scale_index
module Index
Definition: scale_index.F90:11
scale_const::const_undef2
integer, parameter, public const_undef2
undefined value (INT2)
Definition: scale_const.F90:38
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:228
scale_precision
module PRECISION
Definition: scale_precision.F90:14
scale_atmos_grid_cartesc_index::ka
integer, public ka
Definition: scale_atmos_grid_cartesC_index.F90:47
scale_atmos_dyn_fvm_flux
module scale_atmos_dyn_fvm_flux
Definition: scale_atmos_dyn_fvm_flux.F90:13
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_uv
integer, public i_uv
Definition: scale_atmos_grid_cartesC_index.F90:101
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:196
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:168
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_fluxx_xvz
procedure(flux_mom), pointer, public atmos_dyn_fvm_fluxx_xvz
Definition: scale_atmos_dyn_fvm_flux.F90:239
scale_index::i_dens
integer, parameter, public i_dens
Definition: scale_index.F90:28
scale_atmos_dyn_tstep_short_fvm_hivi::atmos_dyn_tstep_short_fvm_hivi_regist
subroutine, public atmos_dyn_tstep_short_fvm_hivi_regist(ATMOS_DYN_TYPE, VA_out, VAR_NAME, VAR_DESC, VAR_UNIT)
Register.
Definition: scale_atmos_dyn_tstep_short_fvm_hivi.F90:79
scale_atmos_dyn_tstep_short_fvm_hivi::atmos_dyn_tstep_short_fvm_hivi
subroutine, public atmos_dyn_tstep_short_fvm_hivi(DENS_RK, MOMZ_RK, MOMX_RK, MOMY_RK, RHOT_RK, PROG_RK, mflx_hi, tflx_hi, DENS0, MOMZ0, MOMX0, MOMY0, RHOT0, DENS, MOMZ, MOMX, MOMY, RHOT, DENS_t, MOMZ_t, MOMX_t, MOMY_t, RHOT_t, PROG0, PROG, DPRES0, RT2P, CORIOLI, num_diff, wdamp_coef, divdmp_coef, DDIV, FLAG_FCT_MOMENTUM, FLAG_FCT_T, FLAG_FCT_ALONG_STREAM, CDZ, FDZ, FDX, FDY, RCDZ, RCDX, RCDY, RFDZ, RFDX, RFDY, PHI, GSQRT, J13G, J23G, J33G, MAPF, REF_dens, REF_rhot, BND_W, BND_E, BND_S, BND_N, TwoD, dtrk, last)
Definition: scale_atmos_dyn_tstep_short_fvm_hivi.F90:179
scale_index::i_momy
integer, parameter, public i_momy
Definition: scale_index.F90:31
scale_atmos_grid_cartesc_index::i_xyz
integer, public i_xyz
Definition: scale_atmos_grid_cartesC_index.F90:90
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:190
scale_index::va
integer, public va
Definition: scale_index.F90:35
scale_prc
module PROCESS
Definition: scale_prc.F90:11
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:185
scale_atmos_grid_cartesc_index::i_xy
integer, public i_xy
Definition: scale_atmos_grid_cartesC_index.F90:98
scale_atmos_grid_cartesc_index::ie
integer, public ie
end point of inner domain: x, local
Definition: scale_atmos_grid_cartesC_index.F90:54
scale_atmos_grid_cartesc_index::i_uyz
integer, public i_uyz
Definition: scale_atmos_grid_cartesC_index.F90:94
scale_io
module STDIO
Definition: scale_io.F90:10
scale_atmos_grid_cartesc_index::i_xvw
integer, public i_xvw
Definition: scale_atmos_grid_cartesC_index.F90:93
scale_tracer::k
real(rp), public k
Definition: scale_tracer.F90:44
scale_atmos_grid_cartesc_index
module atmosphere / grid / cartesC index
Definition: scale_atmos_grid_cartesC_index.F90:12
scale_const
module CONSTANT
Definition: scale_const.F90:11
scale_atmos_grid_cartesc_index::ia
integer, public ia
Definition: scale_atmos_grid_cartesC_index.F90:48
scale_debug::check
subroutine, public check(current_line, v)
Undefined value checker.
Definition: scale_debug.F90:56
scale_atmos_grid_cartesc_index::zdir
integer, parameter, public zdir
Definition: scale_atmos_grid_cartesC_index.F90:32
scale_atmos_dyn_tstep_short_fvm_hivi::atmos_dyn_tstep_short_fvm_hivi_setup
subroutine, public atmos_dyn_tstep_short_fvm_hivi_setup
Setup.
Definition: scale_atmos_dyn_tstep_short_fvm_hivi.F90:112
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_fluxz_xvz
procedure(flux_z), pointer, public atmos_dyn_fvm_fluxz_xvz
Definition: scale_atmos_dyn_fvm_flux.F90:234
scale_prof
module profiler
Definition: scale_prof.F90:11
scale_index::i_rhot
integer, parameter, public i_rhot
Definition: scale_index.F90:32
scale_atmos_grid_cartesc_index::i_uyw
integer, public i_uyw
Definition: scale_atmos_grid_cartesC_index.F90:92
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_precision::dp
integer, parameter, public dp
Definition: scale_precision.F90:32
scale_atmos_grid_cartesc_index::ja
integer, public ja
Definition: scale_atmos_grid_cartesC_index.F90:49
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:159
scale_tracer
module TRACER
Definition: scale_tracer.F90:12
scale_comm_cartesc::comm_world
integer, public comm_world
communication world ID
Definition: scale_comm_cartesC.F90:99
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:207
scale_atmos_grid_cartesc_index::xdir
integer, parameter, public xdir
Definition: scale_atmos_grid_cartesC_index.F90:33
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_debug
module DEBUG
Definition: scale_debug.F90:11
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:255
scale_comm_cartesc
module COMMUNICATION
Definition: scale_comm_cartesC.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::ydir
integer, parameter, public ydir
Definition: scale_atmos_grid_cartesC_index.F90:34
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:250
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:163
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:223
scale_atmos_dyn_tstep_short_fvm_hivi
module Atmosphere / Dynamics RK
Definition: scale_atmos_dyn_tstep_short_fvm_hivi.F90:16
scale_atmos_grid_cartesc_index::i_uvz
integer, public i_uvz
Definition: scale_atmos_grid_cartesC_index.F90:96
scale_const::const_undef
real(rp), public const_undef
Definition: scale_const.F90:41
scale_atmos_dyn_fvm_fct
module Atmosphere / Dynamics common
Definition: scale_atmos_dyn_fvm_fct.F90:12
scale_const::const_pre00
real(rp), public const_pre00
pressure reference [Pa]
Definition: scale_const.F90:88
scale_io::io_fid_conf
integer, public io_fid_conf
Config file ID.
Definition: scale_io.F90:56
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_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:180
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:244
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:217
scale_atmos_grid_cartesc_index::i_xyw
integer, public i_xyw
Definition: scale_atmos_grid_cartesC_index.F90:91
Generated by   doxygen 1.8.17