SCALE-RM
scale_atmos_dyn_tstep_short_fvm_heve.F90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
13 !-------------------------------------------------------------------------------
14 #include "inc_openmp.h"
15 module scale_atmos_dyn_tstep_short_fvm_heve
16  !-----------------------------------------------------------------------------
17  !
18  !++ used modules
19  !
20  use scale_precision
21  use scale_stdio
22  use scale_prof
23  use scale_grid_index
24  use scale_index
25  use scale_tracer
26  use scale_process
27 #if defined DEBUG || defined QUICKDEBUG
28  use scale_debug, only: &
29  check
30  use scale_const, only: &
31  undef => const_undef, &
32  iundef => const_undef2
33 #endif
34  !-----------------------------------------------------------------------------
35  implicit none
36  private
37  !-----------------------------------------------------------------------------
38  !
39  !++ Public procedure
40  !
41  public :: atmos_dyn_tstep_short_fvm_heve_regist
42  public :: atmos_dyn_tstep_short_fvm_heve_setup
43  public :: atmos_dyn_tstep_short_fvm_heve
44 
45  !-----------------------------------------------------------------------------
46  !
47  !++ Public parameters & variables
48  !
49  !-----------------------------------------------------------------------------
50  !
51  !++ Private procedure
52  !
53 #if 1
54 #define F2H(k,p,idx) (CDZ(k+p-1)*GSQRT(k+p-1,i,j,idx)/(CDZ(k)*GSQRT(k,i,j,idx)+CDZ(k+1)*GSQRT(k+1,i,j,idx)))
55 #else
56 #define F2H(k,p,idx) 0.5_RP
57 #endif
58 
59  !-----------------------------------------------------------------------------
60  !
61  !++ Private parameters & variables
62  !
63  integer, parameter :: VA_FVM_HEVE = 0
64 
65  !-----------------------------------------------------------------------------
66 contains
67  !-----------------------------------------------------------------------------
69  subroutine atmos_dyn_tstep_short_fvm_heve_regist( &
70  ATMOS_DYN_TYPE, &
71  VA_out, &
72  VAR_NAME, &
73  VAR_DESC, &
74  VAR_UNIT )
75  use scale_process, only: &
76  prc_mpistop
77  implicit none
78 
79  character(len=*), intent(in) :: atmos_dyn_type
80  integer, intent(out) :: va_out
81  character(len=H_SHORT), intent(out) :: var_name(:)
82  character(len=H_MID), intent(out) :: var_desc(:)
83  character(len=H_SHORT), intent(out) :: var_unit(:)
84  !---------------------------------------------------------------------------
85 
86  if ( atmos_dyn_type /= 'FVM-HEVE' .AND. atmos_dyn_type /= 'HEVE' ) then
87  write(*,*) 'xxx ATMOS_DYN_TYPE is not FVM-HEVE. Check!'
88  call prc_mpistop
89  endif
90 
91  va_out = va_fvm_heve
92  var_name(:) = ""
93  var_desc(:) = ""
94  var_unit(:) = ""
95 
96  if( io_l ) write(io_fid_log,*)
97  if( io_l ) write(io_fid_log,*) '*** Register additional prognostic variables (HEVE)'
98  if ( va_out < 1 ) then
99  if( io_l ) write(io_fid_log,*) '*** => nothing.'
100  endif
101 
102  return
103  end subroutine atmos_dyn_tstep_short_fvm_heve_regist
104 
105  !-----------------------------------------------------------------------------
107  subroutine atmos_dyn_tstep_short_fvm_heve_setup
108  return
109  end subroutine atmos_dyn_tstep_short_fvm_heve_setup
110 
111  !-----------------------------------------------------------------------------
112  subroutine atmos_dyn_tstep_short_fvm_heve( &
113  DENS_RK, MOMZ_RK, MOMX_RK, MOMY_RK, RHOT_RK, &
114  PROG_RK, &
115  mflx_hi, tflx_hi, &
116  DENS0, MOMZ0, MOMX0, MOMY0, RHOT0, &
117  DENS, MOMZ, MOMX, MOMY, RHOT, &
118  DENS_t, MOMZ_t, MOMX_t, MOMY_t, RHOT_t, &
119  PROG0, PROG, &
120  DPRES0, RT2P, CORIOLI, &
121  num_diff, wdamp_coef, divdmp_coef, DDIV, &
122  FLAG_FCT_MOMENTUM, FLAG_FCT_T, &
123  FLAG_FCT_ALONG_STREAM, &
124  CDZ, FDZ, FDX, FDY, &
125  RCDZ, RCDX, RCDY, RFDZ, RFDX, RFDY, &
126  PHI, GSQRT, J13G, J23G, J33G, MAPF, &
127  REF_dens, REF_rhot, &
128  BND_W, BND_E, BND_S, BND_N, &
129  dtrk, last )
130  use scale_grid_index
131  use scale_const, only: &
132  eps => const_eps, &
133  grav => const_grav, &
134  p00 => const_pre00
135  use scale_comm, only: &
136  comm_vars8, &
137  comm_wait
138  use scale_atmos_dyn_common, only: &
139  atmos_dyn_fct
140  use scale_atmos_dyn_fvm_flux_ud1, only: &
141  atmos_dyn_fvm_fluxz_xyw_ud1, &
142  atmos_dyn_fvm_fluxx_xyw_ud1, &
143  atmos_dyn_fvm_fluxy_xyw_ud1, &
144  atmos_dyn_fvm_fluxz_uyz_ud1, &
145  atmos_dyn_fvm_fluxx_uyz_ud1, &
146  atmos_dyn_fvm_fluxy_uyz_ud1, &
147  atmos_dyn_fvm_fluxz_xvz_ud1, &
148  atmos_dyn_fvm_fluxx_xvz_ud1, &
149  atmos_dyn_fvm_fluxy_xvz_ud1, &
150  atmos_dyn_fvm_fluxz_xyz_ud1, &
151  atmos_dyn_fvm_fluxx_xyz_ud1, &
152  atmos_dyn_fvm_fluxy_xyz_ud1
153  use scale_atmos_dyn_fvm_flux, only: &
154  atmos_dyn_fvm_fluxz_xyz, &
155  atmos_dyn_fvm_fluxx_xyz, &
156  atmos_dyn_fvm_fluxy_xyz, &
157  atmos_dyn_fvm_fluxz_xyw, &
158  atmos_dyn_fvm_fluxj13_xyw, &
159  atmos_dyn_fvm_fluxj23_xyw, &
160  atmos_dyn_fvm_fluxx_xyw, &
161  atmos_dyn_fvm_fluxy_xyw, &
162  atmos_dyn_fvm_fluxz_uyz, &
163  atmos_dyn_fvm_fluxj13_uyz, &
164  atmos_dyn_fvm_fluxj23_uyz, &
165  atmos_dyn_fvm_fluxx_uyz, &
166  atmos_dyn_fvm_fluxy_uyz, &
167  atmos_dyn_fvm_fluxz_xvz, &
168  atmos_dyn_fvm_fluxj13_xvz, &
169  atmos_dyn_fvm_fluxj23_xvz, &
170  atmos_dyn_fvm_fluxx_xvz, &
171  atmos_dyn_fvm_fluxy_xvz
172  use scale_gridtrans, only: &
173  i_xyz, &
174  i_xyw, &
175  i_uyw, &
176  i_xvw, &
177  i_uyz, &
178  i_xvz, &
179  i_uvz, &
180  i_xy, &
181  i_uy, &
182  i_xv, &
183  i_uv
184 #ifdef HIST_TEND
185  use scale_history, only: &
186  hist_in
187 #endif
188  implicit none
189 
190  real(RP), intent(out) :: dens_rk (ka,ia,ja) ! prognostic variables
191  real(RP), intent(out) :: momz_rk (ka,ia,ja) !
192  real(RP), intent(out) :: momx_rk (ka,ia,ja) !
193  real(RP), intent(out) :: momy_rk (ka,ia,ja) !
194  real(RP), intent(out) :: rhot_rk (ka,ia,ja) !
195  real(RP), intent(out) :: prog_rk (ka,ia,ja,va) !
196 
197  real(RP), intent(inout) :: mflx_hi (ka,ia,ja,3) ! mass flux
198  real(RP), intent(out) :: tflx_hi (ka,ia,ja,3) ! internal energy flux
199 
200  real(RP), intent(in), target :: dens0 (ka,ia,ja) ! prognostic variables at previous dynamical time step
201  real(RP), intent(in), target :: momz0 (ka,ia,ja) !
202  real(RP), intent(in), target :: momx0 (ka,ia,ja) !
203  real(RP), intent(in), target :: momy0 (ka,ia,ja) !
204  real(RP), intent(in), target :: rhot0 (ka,ia,ja) !
205  real(RP), intent(in) :: prog0 (ka,ia,ja,va)
206 
207  real(RP), intent(in) :: dens (ka,ia,ja) ! prognostic variables at previous RK step
208  real(RP), intent(in) :: momz (ka,ia,ja) !
209  real(RP), intent(in) :: momx (ka,ia,ja) !
210  real(RP), intent(in) :: momy (ka,ia,ja) !
211  real(RP), intent(in) :: rhot (ka,ia,ja) !
212  real(RP), intent(in) :: prog (ka,ia,ja,va)
213 
214  real(RP), intent(in) :: dens_t (ka,ia,ja) ! tendency
215  real(RP), intent(in) :: momz_t (ka,ia,ja) !
216  real(RP), intent(in) :: momx_t (ka,ia,ja) !
217  real(RP), intent(in) :: momy_t (ka,ia,ja) !
218  real(RP), intent(in) :: rhot_t (ka,ia,ja) !
219 
220  real(RP), intent(in) :: dpres0 (ka,ia,ja)
221  real(RP), intent(in) :: rt2p (ka,ia,ja)
222  real(RP), intent(in) :: corioli (1, ia,ja)
223  real(RP), intent(in) :: num_diff(ka,ia,ja,5,3)
224  real(RP), intent(in) :: wdamp_coef(ka)
225  real(RP), intent(in) :: divdmp_coef
226  real(RP), intent(in) :: ddiv (ka,ia,ja)
227 
228  logical, intent(in) :: flag_fct_momentum
229  logical, intent(in) :: flag_fct_t
230  logical, intent(in) :: flag_fct_along_stream
231 
232  real(RP), intent(in) :: cdz (ka)
233  real(RP), intent(in) :: fdz (ka-1)
234  real(RP), intent(in) :: fdx (ia-1)
235  real(RP), intent(in) :: fdy (ja-1)
236  real(RP), intent(in) :: rcdz(ka)
237  real(RP), intent(in) :: rcdx(ia)
238  real(RP), intent(in) :: rcdy(ja)
239  real(RP), intent(in) :: rfdz(ka-1)
240  real(RP), intent(in) :: rfdx(ia-1)
241  real(RP), intent(in) :: rfdy(ja-1)
242 
243  real(RP), intent(in) :: phi (ka,ia,ja)
244  real(RP), intent(in) :: gsqrt (ka,ia,ja,7)
245  real(RP), intent(in) :: j13g (ka,ia,ja,7)
246  real(RP), intent(in) :: j23g (ka,ia,ja,7)
247  real(RP), intent(in) :: j33g
248  real(RP), intent(in) :: mapf (ia,ja,2,4)
249  real(RP), intent(in) :: ref_dens(ka,ia,ja)
250  real(RP), intent(in) :: ref_rhot(ka,ia,ja)
251 
252  logical, intent(in) :: bnd_w
253  logical, intent(in) :: bnd_e
254  logical, intent(in) :: bnd_s
255  logical, intent(in) :: bnd_n
256 
257  real(RP), intent(in) :: dtrk
258  logical, intent(in) :: last
259 
260  ! diagnostic variables
261  real(RP) :: velz (ka,ia,ja) ! velocity w [m/s]
262  real(RP) :: velx (ka,ia,ja) ! velocity u [m/s]
263  real(RP) :: vely (ka,ia,ja) ! velocity v [m/s]
264  real(RP) :: pott (ka,ia,ja) ! potential temperature [K]
265  real(RP) :: dpres(ka,ia,ja) ! pressure - reference pressure
266 
267  real(RP) :: qflx_j13(ka,ia,ja)
268  real(RP) :: qflx_j23(ka,ia,ja)
269  real(RP) :: pgf (ka,ia,ja) ! pressure gradient force
270  real(RP) :: buoy (ka,ia,ja) ! buoyancy force
271  real(RP) :: cor (ka,ia,ja) ! Coriolis force
272 
273  ! flux
274  real(RP) :: qflx_hi (ka,ia,ja,3)
275 #ifndef NO_FCT_DYN
276  real(RP) :: qflx_lo (ka,ia,ja,3)
277  real(RP) :: qflx_anti(ka,ia,ja,3)
278  real(RP) :: tflx_lo (ka,ia,ja,3)
279  real(RP) :: tflx_anti(ka,ia,ja,3)
280  real(RP) :: dens0_uvw(ka,ia,ja)
281  real(RP) :: dens_uvw (ka,ia,ja)
282 #endif
283  real(RP) :: advch ! horizontal advection
284  real(RP) :: advcv ! vertical advection
285  real(RP) :: wdamp ! rayleight damping for W
286  real(RP) :: div ! divergence damping
287 #ifdef HIST_TEND
288  real(RP) :: advch_t(ka,ia,ja,5)
289  real(RP) :: advcv_t(ka,ia,ja,5)
290  real(RP) :: wdmp_t(ka,ia,ja)
291  real(RP) :: ddiv_t(ka,ia,ja,3)
292  real(RP) :: pg_t(ka,ia,ja,3)
293  real(RP) :: cf_t(ka,ia,ja,2)
294  logical :: lhist
295 #endif
296 
297  integer :: iis, iie
298  integer :: jjs, jje
299  integer :: ifs_off, jfs_off
300  integer :: k, i, j
301  !---------------------------------------------------------------------------
302 
303 #ifdef DEBUG
304  velz(:,:,:) = undef
305  velx(:,:,:) = undef
306  vely(:,:,:) = undef
307  pott(:,:,:) = undef
308 
309  dpres(:,:,:) = undef
310 
311  mflx_hi(:,:,:,:) = undef
312  tflx_hi(:,:,:,:) = undef
313  qflx_hi(:,:,:,:) = undef
314 
315 #ifndef NO_FCT_DYN
316  qflx_lo(:,:,:,:) = undef
317  qflx_anti(:,:,:,:) = undef
318  tflx_lo(:,:,:,:) = undef
319  tflx_anti(:,:,:,:) = undef
320 #endif
321 #endif
322 
323 #if defined DEBUG || defined QUICKDEBUG
324  dens_rk( 1:ks-1,:,:) = undef
325  dens_rk(ke+1:ka ,:,:) = undef
326  momz_rk( 1:ks-1,:,:) = undef
327  momz_rk(ke+1:ka ,:,:) = undef
328  momx_rk( 1:ks-1,:,:) = undef
329  momx_rk(ke+1:ka ,:,:) = undef
330  momy_rk( 1:ks-1,:,:) = undef
331  momy_rk(ke+1:ka ,:,:) = undef
332  rhot_rk( 1:ks-1,:,:) = undef
333  rhot_rk(ke+1:ka ,:,:) = undef
334  prog_rk( 1:ks-1,:,:,:) = undef
335  prog_rk(ke+1:ka ,:,:,:) = undef
336 #endif
337 
338 #ifdef HIST_TEND
339  advch_t = 0.0_rp
340  advcv_t = 0.0_rp
341  wdmp_t = 0.0_rp
342  ddiv_t = 0.0_rp
343  pg_t = 0.0_rp
344  cf_t = 0.0_rp
345 
346  lhist = last
347 #endif
348 
349  ifs_off = 1
350  jfs_off = 1
351  if ( bnd_w ) ifs_off = 0
352  if ( bnd_s ) jfs_off = 0
353 
354 
355  do jjs = js, je, jblock
356  jje = jjs+jblock-1
357  do iis = is, ie, iblock
358  iie = iis+iblock-1
359 
360  ! pressure, pott. temp.
361 
362  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
363  do j = jjs-2, jje+2
364  do i = iis-2, iie+2
365  do k = ks, ke
366 #ifdef DEBUG
367  call check( __line__, dpres0(k,i,j) )
368  call check( __line__, rt2p(k,i,j) )
369  call check( __line__, rhot(k,i,j) )
370  call check( __line__, ref_rhot(k,i,j) )
371 #endif
372  dpres(k,i,j) = dpres0(k,i,j) + rt2p(k,i,j) * ( rhot(k,i,j) - ref_rhot(k,i,j) )
373  enddo
374  dpres(ks-1,i,j) = dpres0(ks-1,i,j) - dens(ks,i,j) * ( phi(ks-1,i,j) - phi(ks+1,i,j) )
375  dpres(ke+1,i,j) = dpres0(ke+1,i,j) - dens(ke,i,j) * ( phi(ke+1,i,j) - phi(ke-1,i,j) )
376  enddo
377  enddo
378 #ifdef DEBUG
379  k = iundef; i = iundef; j = iundef
380 #endif
381 
382  do j = jjs-jhalo, jje+jhalo
383  do i = iis-ihalo, iie+ihalo
384  do k = ks, ke
385 #ifdef DEBUG
386  call check( __line__, rhot(k,i,j) )
387  call check( __line__, dens(k,i,j) )
388 #endif
389  pott(k,i,j) = rhot(k,i,j) / dens(k,i,j)
390  enddo
391  enddo
392  enddo
393 #ifdef DEBUG
394  k = iundef; i = iundef; j = iundef
395 #endif
396 
397 #ifndef NO_FCT_DYN
398  enddo
399  enddo
400 
401  if ( flag_fct_momentum ) then
402  ! momentum -> velocity
403  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
404  do j = js-1, je+2
405  do i = is-1, ie+2
406  do k = ks, ke-1
407 #ifdef DEBUG
408  call check( __line__, momz0(k,i,j) )
409  call check( __line__, dens0(k ,i,j) )
410  call check( __line__, dens0(k+1,i,j) )
411 #endif
412  velz(k,i,j) = 2.0_rp * momz0(k,i,j) / ( dens0(k+1,i,j)+dens0(k,i,j) )
413  enddo
414  enddo
415  enddo
416 #ifdef DEBUG
417  k = iundef; i = iundef; j = iundef
418 #endif
419  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
420  do j = js-1, je+2
421  do i = is-1, ie+2
422  velz(ke,i,j) = 0.0_rp
423  enddo
424  enddo
425 #ifdef DEBUG
426  k = iundef; i = iundef; j = iundef
427 #endif
428 
429  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
430  do j = js-1, je+2
431  do i = is-2, ie+1
432  do k = ks, ke
433 #ifdef DEBUG
434  call check( __line__, momx0(k,i,j) )
435  call check( __line__, dens0(k,i ,j) )
436  call check( __line__, dens0(k,i+1,j) )
437 #endif
438  velx(k,i,j) = 2.0_rp * momx0(k,i,j) / ( dens0(k,i+1,j)+dens0(k,i,j) )
439  enddo
440  enddo
441  enddo
442 #ifdef DEBUG
443  k = iundef; i = iundef; j = iundef
444 #endif
445 
446  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
447  do j = js-2, je+1
448  do i = is-1, ie+2
449  do k = ks, ke
450 #ifdef DEBUG
451  call check( __line__, momy0(k,i,j) )
452  call check( __line__, dens0(k,i,j ) )
453  call check( __line__, dens0(k,i,j+1) )
454 #endif
455  vely(k,i,j) = 2.0_rp * momy0(k,i,j) / ( dens0(k,i,j+1)+dens0(k,i,j) )
456  enddo
457  enddo
458  enddo
459 #ifdef DEBUG
460  k = iundef; i = iundef; j = iundef
461 #endif
462  call comm_vars8( velz(:,:,:), 4 )
463  call comm_vars8( velx(:,:,:), 5 )
464  call comm_vars8( vely(:,:,:), 6 )
465  endif
466 
467  do jjs = js, je, jblock
468  jje = jjs+jblock-1
469  do iis = is, ie, iblock
470  iie = iis+iblock-1
471 #endif
472 
473  !########################################################################
474  ! continuity equation (total rho)
475  !########################################################################
476 
477  !-----< high order flux >-----
478 
479  ! at (x, y, w)
480 
481  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
482  do j = jjs, jje
483  do i = iis, iie
484  do k = ks, ke-1
485 #ifdef DEBUG
486  call check( __line__, momz(k+1,i,j) )
487  call check( __line__, momz(k ,i,j) )
488  call check( __line__, momz(k-1,i,j) )
489  call check( __line__, num_diff(k,i,j,i_dens,zdir) )
490 #endif
491  mflx_hi(k,i,j,zdir) = j33g * momz(k,i,j) / ( mapf(i,j,1,i_xy)*mapf(i,j,2,i_xy) ) &
492  + j13g(k,i,j,i_xyw) * 0.25_rp * ( momx(k+1,i,j)+momx(k+1,i-1,j) &
493  + momx(k ,i,j)+momx(k ,i-1,j) ) &
494  / mapf(i,j,2,i_xy) & ! [{u,y,z->x,y,w}]
495  + j23g(k,i,j,i_xyw) * 0.25_rp * ( momy(k+1,i,j)+momy(k+1,i,j-1) &
496  + momy(k ,i,j)+momy(k ,i,j-1) ) &
497  / mapf(i,j,1,i_xy) & ! [{x,v,z->x,y,w}]
498  + gsqrt(k,i,j,i_xyw) * num_diff(k,i,j,i_dens,zdir) / ( mapf(i,j,1,i_xy)*mapf(i,j,2,i_xy) )
499  enddo
500  enddo
501  enddo
502 #ifdef DEBUG
503  k = iundef; i = iundef; j = iundef
504 #endif
505 
506  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
507  do j = jjs, jje
508  do i = iis, iie
509  mflx_hi(ks-1,i,j,zdir) = 0.0_rp
510  mflx_hi(ke ,i,j,zdir) = 0.0_rp
511  enddo
512  enddo
513 #ifdef DEBUG
514  k = iundef; i = iundef; j = iundef
515 #endif
516 
517  ! at (u, y, z)
518 
519  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
520  do j = jjs , jje
521  do i = iis-ifs_off, min(iie,ieh)
522  do k = ks, ke
523 #ifdef DEBUG
524  call check( __line__, momx(k,i+1,j) )
525  call check( __line__, momx(k,i ,j) )
526  call check( __line__, momx(k,i-1,j) )
527  call check( __line__, num_diff(k,i,j,i_dens,xdir) )
528 #endif
529  mflx_hi(k,i,j,xdir) = gsqrt(k,i,j,i_uyz) / mapf(i,j,2,i_uy) &
530  * ( momx(k,i,j) + num_diff(k,i,j,i_dens,xdir) )
531  enddo
532  enddo
533  enddo
534 #ifdef DEBUG
535  k = iundef; i = iundef; j = iundef
536 #endif
537 
538  ! at (x, v, z)
539 
540  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
541  do j = jjs-jfs_off, min(jje,jeh)
542  do i = iis , iie
543  do k = ks, ke
544 #ifdef DEBUG
545  call check( __line__, momy(k,i,j+1) )
546  call check( __line__, momy(k,i,j ) )
547  call check( __line__, momy(k,i,j-1) )
548  call check( __line__, num_diff(k,i,j,i_dens,ydir) )
549 #endif
550  mflx_hi(k,i,j,ydir) = gsqrt(k,i,j,i_xvz) / mapf(i,j,1,i_xv) &
551  * ( momy(k,i,j) + num_diff(k,i,j,i_dens,ydir) )
552  enddo
553  enddo
554  enddo
555 #ifdef DEBUG
556  k = iundef; i = iundef; j = iundef
557 #endif
558 
559  !-----< update density >-----
560 
561  !$omp parallel do private(i,j,k,advcv,advch) OMP_SCHEDULE_ collapse(2)
562  do j = jjs, jje
563  do i = iis, iie
564  do k = ks, ke
565 #ifdef DEBUG
566  call check( __line__, dens0(k,i,j) )
567  call check( __line__, mflx_hi(k ,i ,j ,zdir) )
568  call check( __line__, mflx_hi(k-1,i ,j ,zdir) )
569  call check( __line__, mflx_hi(k ,i ,j ,xdir) )
570  call check( __line__, mflx_hi(k ,i-1,j ,xdir) )
571  call check( __line__, mflx_hi(k ,i ,j ,ydir) )
572  call check( __line__, mflx_hi(k ,i ,j-1,ydir) )
573  call check( __line__, dens_t(k,i,j) )
574 #endif
575  advcv = - ( mflx_hi(k,i,j,zdir)-mflx_hi(k-1,i, j, zdir) ) * rcdz(k)
576  advch = - ( mflx_hi(k,i,j,xdir)-mflx_hi(k ,i-1,j, xdir) ) * rcdx(i) &
577  - ( mflx_hi(k,i,j,ydir)-mflx_hi(k ,i, j-1,ydir) ) * rcdy(j)
578  dens_rk(k,i,j) = dens0(k,i,j) &
579  + dtrk * ( ( advcv + advch ) * mapf(i,j,1,i_xy) * mapf(i,j,2,i_xy) / gsqrt(k,i,j,i_xyz) &
580  + dens_t(k,i,j) )
581 #ifdef HIST_TEND
582  if ( lhist ) then
583  advcv_t(k,i,j,i_dens) = advcv * mapf(i,j,1,i_xy) * mapf(i,j,2,i_xy) / gsqrt(k,i,j,i_xyz)
584  advch_t(k,i,j,i_dens) = advch * mapf(i,j,1,i_xy) * mapf(i,j,2,i_xy) / gsqrt(k,i,j,i_xyz)
585  endif
586 #endif
587  enddo
588  enddo
589  enddo
590 #ifdef DEBUG
591  k = iundef; i = iundef; j = iundef
592 #endif
593 
594 
595  !########################################################################
596  ! momentum equation (z)
597  !########################################################################
598 
599  !-----< high order flux >-----
600 
601  ! at (x, y, z)
602  ! note than z-index is added by -1
603  call atmos_dyn_fvm_fluxz_xyw( qflx_hi(:,:,:,zdir), & ! (out)
604  momz, momz, dens, & ! (in)
605  gsqrt(:,:,:,i_xyz), j33g, & ! (in)
606  num_diff(:,:,:,i_momz,zdir), & ! (in)
607  cdz, fdz, dtrk, &
608  iis, iie, jjs, jje ) ! (in)
609  call atmos_dyn_fvm_fluxj13_xyw( qflx_j13, & ! (out)
610  momx, momz, dens, & ! (in)
611  gsqrt(:,:,:,i_xyz), j13g(:,:,:,i_xyz), mapf(:,:,:,i_xy), & ! (in)
612  cdz, &
613  iis, iie, jjs, jje ) ! (in)
614  call atmos_dyn_fvm_fluxj23_xyw( qflx_j23, & ! (out)
615  momy, momz, dens, & ! (in)
616  gsqrt(:,:,:,i_xyz), j23g(:,:,:,i_xyz), mapf(:,:,:,i_xy), & ! (in)
617  cdz, &
618  iis, iie, jjs, jje ) ! (in)
619 
620  ! at (u, y, w)
621  call atmos_dyn_fvm_fluxx_xyw( qflx_hi(:,:,:,xdir), & ! (out)
622  momx, momz, dens, & ! (in)
623  gsqrt(:,:,:,i_uyw), mapf(:,:,:,i_uy), & ! (in)
624  num_diff(:,:,:,i_momz,xdir), & ! (in)
625  cdz, & ! (in)
626  iis, iie, jjs, jje ) ! (in)
627 
628  ! at (x, v, w)
629  call atmos_dyn_fvm_fluxy_xyw( qflx_hi(:,:,:,ydir), & ! (out)
630  momy, momz, dens, & ! (in)
631  gsqrt(:,:,:,i_xvw), mapf(:,:,:,i_xv), & ! (in)
632  num_diff(:,:,:,i_momz,ydir), & ! (in)
633  cdz, & ! (in)
634  iis, iie, jjs, jje ) ! (in)
635 
636  ! pressure gradient force at (x, y, w)
637 
638  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
639  do j = jjs, jje
640  do i = iis, iie
641  do k = ks, ke-1
642  pgf(k,i,j) = j33g * ( dpres(k+1,i,j)-dpres(k,i,j) ) * rfdz(k) ! [x,y,z]
643  enddo
644  enddo
645  enddo
646 
647  ! buoyancy force at (x, y, w)
648 
649  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
650  do j = jjs, jje
651  do i = iis, iie
652  do k = ks, ke-1
653  buoy(k,i,j) = grav * gsqrt(k,i,j,i_xyw) &
654  * ( f2h(k,1,i_xyz) * ( dens(k+1,i,j)-ref_dens(k+1,i,j) ) &
655  + f2h(k,2,i_xyz) * ( dens(k ,i,j)-ref_dens(k ,i,j) ) )
656  enddo
657  enddo
658  enddo
659 
660  !-----< update momentum (z) -----
661 
662  !$omp parallel do private(i,j,k,advcv,advch,wdamp,div) OMP_SCHEDULE_ collapse(2)
663  do j = jjs, jje
664  do i = iis, iie
665  do k = ks, ke-1
666 #ifdef DEBUG
667  call check( __line__, qflx_hi(k ,i ,j ,zdir) )
668  call check( __line__, qflx_hi(k-1,i ,j ,zdir) )
669  call check( __line__, qflx_hi(k ,i ,j ,xdir) )
670  call check( __line__, qflx_hi(k ,i-1,j ,xdir) )
671  call check( __line__, qflx_hi(k ,i ,j ,ydir) )
672  call check( __line__, qflx_hi(k ,i ,j-1,ydir) )
673  call check( __line__, ddiv(k ,i,j) )
674  call check( __line__, ddiv(k+1,i,j) )
675  call check( __line__, momz0(k,i,j) )
676  call check( __line__, momz_t(k,i,j) )
677 #endif
678  advcv = - ( qflx_hi(k,i,j,zdir) - qflx_hi(k-1,i ,j ,zdir) &
679  + qflx_j13(k,i,j) - qflx_j13(k-1,i,j) &
680  + qflx_j23(k,i,j) - qflx_j23(k-1,i,j) ) * rfdz(k)
681  advch = - ( ( qflx_hi(k,i,j,xdir) - qflx_hi(k,i-1,j,xdir) ) * rcdx(i) &
682  + ( qflx_hi(k,i,j,ydir) - qflx_hi(k,i,j-1,ydir) ) * rcdy(j) ) &
683  * mapf(i,j,1,i_xy) * mapf(i,j,2,i_xy)
684  wdamp = - wdamp_coef(k) * momz0(k,i,j)
685  div = divdmp_coef / dtrk * fdz(k) * ( ddiv(k+1,i,j)-ddiv(k,i,j) ) ! divergence damping
686  momz_rk(k,i,j) = momz0(k,i,j) &
687  + dtrk * ( ( advcv + advch &
688  - pgf(k,i,j) & ! pressure gradient force
689  - buoy(k,i,j) & ! buoyancy force
690  ) / gsqrt(k,i,j,i_xyw) &
691  + wdamp & ! Rayleigh damping
692  + div &
693  + momz_t(k,i,j) ) ! physics tendency
694 #ifdef HIST_TEND
695  if ( lhist ) then
696  advcv_t(k,i,j,i_momz) = advcv / gsqrt(k,i,j,i_xyw)
697  advch_t(k,i,j,i_momz) = advch / gsqrt(k,i,j,i_xyw)
698  pg_t(k,i,j,1) = ( - pgf(k,i,j) - buoy(k,i,j) ) / gsqrt(k,i,j,i_xyw)
699  wdmp_t(k,i,j) = wdamp
700  ddiv_t(k,i,j,1) = div
701  endif
702 #endif
703  enddo
704  enddo
705  enddo
706 #ifdef DEBUG
707  k = iundef; i = iundef; j = iundef
708 #endif
709 
710  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
711  do j = jjs, jje
712  do i = iis, iie
713  momz_rk(ks-1,i,j) = 0.0_rp
714  momz_rk(ke ,i,j) = 0.0_rp
715 #ifdef HIST_TEND
716  if ( lhist ) then
717  advcv_t(ke,i,j,i_momz) = 0.0_rp
718  advch_t(ke,i,j,i_momz) = 0.0_rp
719  pg_t(ke,i,j,1) = 0.0_rp
720  wdmp_t(ke,i,j) = 0.0_rp
721  ddiv_t(ke,i,j,1) = 0.0_rp
722  endif
723 #endif
724  enddo
725  enddo
726 #ifdef DEBUG
727  k = iundef; i = iundef; j = iundef
728 #endif
729 
730 #ifndef NO_FCT_DYN
731  if ( flag_fct_momentum ) then
732 
733  ! monotonic flux
734  ! at (x, y, layer)
735  ! note than z-index is added by -1
736  call atmos_dyn_fvm_fluxz_xyw_ud1( qflx_lo(:,:,:,zdir), & ! (out)
737  momz, momz0, dens, & ! (in)
738  gsqrt(:,:,:,i_xyz), j33g, & ! (in)
739  num_diff(:,:,:,i_momz,zdir), & ! (in)
740  cdz, fdz, dtrk, & ! (in)
741  iis-1, iie+1, jjs-1, jje+1 ) ! (in)
742 
743  call atmos_dyn_fvm_fluxx_xyw_ud1( qflx_lo(:,:,:,xdir), & ! (out)
744  momx, momz0, dens, & ! (in)
745  gsqrt(:,:,:,i_uyz), mapf(:,:,:,i_uy), & ! (in)
746  num_diff(:,:,:,i_momz,xdir), & ! (in)
747  cdz, & ! (in)
748  iis-1, iie+1, jjs-1, jje+1 ) ! (in)
749 
750  call atmos_dyn_fvm_fluxy_xyw_ud1( qflx_lo(:,:,:,ydir), & ! (out)
751  momy, momz0, dens, & ! (in)
752  gsqrt(:,:,:,i_xvz), mapf(:,:,:,i_xv), & ! (in)
753  num_diff(:,:,:,i_momz,ydir), & ! (in)
754  cdz, & ! (in)
755  iis-1, iie+1, jjs-1, jje+1 ) ! (in)
756 
757  endif
758  enddo
759  enddo
760 
761  if ( flag_fct_momentum ) then
762 
763  call comm_vars8( dens_rk, 1 )
764  call comm_wait ( dens_rk, 1, .false. )
765 
766  do j = js, je
767  do i = is, ie
768  do k = ks, ke
769  qflx_hi(k,i,j,zdir) = qflx_hi(k,i,j,zdir) / ( mapf(i,j,1,i_xy) * mapf(i,j,2,i_xy) ) &
770  + qflx_j13(k,i,j) + qflx_j23(k,i,j)
771  enddo
772  enddo
773  enddo
774 
775  do j = js-1, je+1
776  do i = is-1, ie+1
777  do k = ks, ke-1
778  dens0_uvw(k,i,j) = 0.5_rp * ( dens0(k,i,j) + dens0(k+1,i,j) )
779  dens_uvw(k,i,j) = 0.5_rp * ( dens_rk(k,i,j) + dens_rk(k+1,i,j) )
780  enddo
781  enddo
782  enddo
783 
784  do j = js-1, je+1
785  do i = is-1, ie+1
786  dens_uvw(ke,i,j) = dens_uvw(ke-1,i,j)
787  dens0_uvw(ke,i,j) = dens0_uvw(ke-1,i,j)
788  enddo
789  enddo
790 
791  call comm_wait ( velz(:,:,:), 4 )
792 
793  call atmos_dyn_fct( qflx_anti, & ! (out)
794  velz, dens0_uvw, dens_uvw, & ! (in)
795  qflx_hi, qflx_lo, & ! (in)
796  mflx_hi, & ! (in)
797  rfdz, rcdx, rcdy, & ! (in)
798  gsqrt(:,:,:,i_xyw), & ! (in)
799  mapf(:,:,:,i_xy), dtrk, & ! (in)
800  flag_fct_along_stream ) ! (in)
801 
802  do jjs = js, je, jblock
803  jje = jjs+jblock-1
804  do iis = is, ie, iblock
805  iie = iis+iblock-1
806 
807  !--- update momentum(z)
808  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
809  do j = jjs, jje
810  do i = iis, iie
811  do k = ks, ke-1
812  momz_rk(k,i,j) = momz_rk(k,i,j) &
813  + dtrk * ( ( qflx_anti(k,i,j,zdir) - qflx_anti(k-1,i ,j ,zdir) ) * rfdz(k) &
814  + ( ( qflx_anti(k,i,j,xdir) - qflx_anti(k ,i-1,j ,xdir) ) * rcdx(i) &
815  + ( qflx_anti(k,i,j,ydir) - qflx_anti(k ,i ,j-1,ydir) ) * rcdy(j) ) &
816  * mapf(i,j,1,i_xy) * mapf(i,j,2,i_xy) ) &
817  / gsqrt(k,i,j,i_xyw)
818  enddo
819  enddo
820  enddo
821 
822  enddo
823  enddo
824 
825  endif ! FLAG_FCT_MOMENTUM
826 
827 #ifdef DEBUG
828  qflx_hi(:,:,:,:) = undef
829 #endif
830 
831  do jjs = js, je, jblock
832  jje = jjs+jblock-1
833  do iis = is, ie, iblock
834  iie = iis+iblock-1
835 #endif
836 
837  !########################################################################
838  ! momentum equation (x)
839  !########################################################################
840 
841  !-----< high order flux >-----
842 
843  ! at (u, y, w)
844  call atmos_dyn_fvm_fluxz_uyz( qflx_hi(:,:,:,zdir), & ! (out)
845  momz, momx, dens, & ! (in)
846  gsqrt(:,:,:,i_uyw), j33g, & ! (in)
847  num_diff(:,:,:,i_momx,zdir), & ! (in)
848  cdz, & ! (in)
849  iis, iie, jjs, jje ) ! (in)
850  call atmos_dyn_fvm_fluxj13_uyz( qflx_j13, & ! (out)
851  momx, momx, dens, & ! (in)
852  gsqrt(:,:,:,i_uyz), j13g(:,:,:,i_uyw), mapf(:,:,:,i_uy), & ! (in)
853  cdz, & ! (in)
854  iis, iie, jjs, jje ) ! (in)
855  call atmos_dyn_fvm_fluxj23_uyz( qflx_j23, & ! (out)
856  momy, momx, dens, & ! (in)
857  gsqrt(:,:,:,i_uyz), j23g(:,:,:,i_uyw), mapf(:,:,:,i_uy), & ! (in)
858  cdz, & ! (in)
859  iis, iie, jjs, jje ) ! (in)
860 
861  ! at (x, y, z)
862  ! note that x-index is added by -1
863  call atmos_dyn_fvm_fluxx_uyz( qflx_hi(:,:,:,xdir), & ! (out)
864  momx, momx, dens, & ! (in)
865  gsqrt(:,:,:,i_xyz), mapf(:,:,:,i_xy), & ! (in)
866  num_diff(:,:,:,i_momx,xdir), & ! (in)
867  cdz, & ! (in)
868  iis, iie, jjs, jje ) ! (in)
869 
870  ! at (u, v, z)
871  call atmos_dyn_fvm_fluxy_uyz( qflx_hi(:,:,:,ydir), & ! (out)
872  momy, momx, dens, & ! (in)
873  gsqrt(:,:,:,i_uvz), mapf(:,:,1,i_uv), & ! (in)
874  num_diff(:,:,:,i_momx,ydir), & ! (in)
875  cdz, & ! (in)
876  iis, iie, jjs, jje ) ! (in)
877 
878  ! pressure gradient force at (u, y, z)
879  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
880  do j = jjs, jje
881  do i = iis, iie
882  do k = ks, ke
883  pgf(k,i,j) = ( ( gsqrt(k,i+1,j,i_xyz) * dpres(k,i+1,j) & ! [x,y,z]
884  - gsqrt(k,i ,j,i_xyz) * dpres(k,i ,j) & ! [x,y,z]
885  ) * rfdx(i) &
886  + ( j13g(k ,i,j,i_uyw) &
887  * 0.5_rp * ( f2h(k,1,i_uyz) * ( dpres(k+1,i+1,j)+dpres(k+1,i,j) ) &
888  + f2h(k,2,i_uyz) * ( dpres(k ,i+1,j)+dpres(k ,i,j) ) ) & ! [x,y,z->u,y,w]
889  - j13g(k-1,i,j,i_uyw) &
890  * 0.5_rp * ( f2h(k,1,i_uyz) * ( dpres(k ,i+1,j)+dpres(k ,i,j) ) &
891  + f2h(k,2,i_uyz) * ( dpres(k-1,i+1,j)+dpres(k-1,i,j) ) ) & ! [x,y,z->u,y,w]
892  ) * rcdz(k) ) &
893  * mapf(i,j,1,i_uy)
894  enddo
895  enddo
896  enddo
897 
898  ! coriolis force at (u, y, z)
899 
900  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
901  do j = jjs, jje
902  do i = iis, iie
903  do k = ks, ke
904 #ifdef DEBUG
905  call check( __line__, momy(k,i ,j ) )
906  call check( __line__, momy(k,i+1,j ) )
907  call check( __line__, momy(k,i ,j-1) )
908  call check( __line__, momy(k,i+1,j-1) )
909 #endif
910  cor(k,i,j) = 0.125_rp * ( corioli(1,i+1,j )+corioli(1,i,j ) ) & ! [x,y,z->u,y,z]
911  * ( momy(k,i+1,j )+momy(k,i,j ) &
912  + momy(k,i+1,j-1)+momy(k,i,j-1) ) & ! [x,v,z->u,y,z]
913  + 0.25_rp * mapf(i,j,1,i_uy) * mapf(i,j,2,i_uy) &
914  * ( momy(k,i,j) + momy(k,i,j-1) + momy(k,i+1,j) + momy(k,i+1,j-1) ) &
915  * ( ( momy(k,i,j) + momy(k,i,j-1) + momy(k,i+1,j) + momy(k,i+1,j-1) ) * 0.25_rp &
916  * ( 1.0_rp/mapf(i+1,j,2,i_xy) - 1.0_rp/mapf(i,j,2,i_xy) ) * rfdx(i) &
917  - momx(k,i,j) &
918  * ( 1.0_rp/mapf(i,j,1,i_uv) - 1.0_rp/mapf(i,j-1,1,i_uv) ) * rcdy(j) ) &
919  * 2.0_rp / ( dens(k,i+1,j) + dens(k,i,j) ) ! metric term
920  enddo
921  enddo
922  enddo
923 
924  !-----< update momentum (x) >-----
925 
926  !$omp parallel do private(i,j,k,advcv,advch,div) OMP_SCHEDULE_ collapse(2)
927  do j = jjs, jje
928  do i = iis, min(iie, ieh)
929  do k = ks, ke
930 #ifdef DEBUG
931  call check( __line__, qflx_hi(k ,i ,j ,zdir) )
932  call check( __line__, qflx_hi(k-1,i ,j ,zdir) )
933  call check( __line__, qflx_hi(k ,i ,j ,xdir) )
934  call check( __line__, qflx_hi(k ,i-1,j ,xdir) )
935  call check( __line__, qflx_hi(k ,i ,j ,ydir) )
936  call check( __line__, qflx_hi(k ,i ,j-1,ydir) )
937  call check( __line__, ddiv(k,i+1,j) )
938  call check( __line__, ddiv(k,i ,j) )
939  call check( __line__, momx0(k,i,j) )
940 #endif
941  ! advection
942  advcv = - ( qflx_hi(k,i,j,zdir) - qflx_hi(k-1,i ,j ,zdir) &
943  + qflx_j13(k,i,j) - qflx_j13(k-1,i,j) &
944  + qflx_j23(k,i,j) - qflx_j23(k-1,i,j) ) * rcdz(k)
945  advch = - ( ( qflx_hi(k,i,j,xdir) - qflx_hi(k ,i-1,j ,xdir) ) * rfdx(i) &
946  + ( qflx_hi(k,i,j,ydir) - qflx_hi(k ,i ,j-1,ydir) ) * rcdy(j) ) &
947  * mapf(i,j,1,i_uy) * mapf(i,j,2,i_uy)
948  div = divdmp_coef / dtrk * fdx(i) * ( ddiv(k,i+1,j)-ddiv(k,i,j) ) ! divergence damping
949  momx_rk(k,i,j) = momx0(k,i,j) &
950  + dtrk * ( ( advcv + advch & ! advection
951  - pgf(k,i,j) & ! pressure gradient force
952  ) / gsqrt(k,i,j,i_uyz) &
953  + cor(k,i,j) & ! coriolis force
954  + div & ! divergence damping
955  + momx_t(k,i,j) ) ! physics tendency
956 #ifdef HIST_TEND
957  if ( lhist ) then
958  advcv_t(k,i,j,i_momx) = advcv / gsqrt(k,i,j,i_uyz)
959  advch_t(k,i,j,i_momx) = advch / gsqrt(k,i,j,i_uyz)
960  pg_t(k,i,j,2) = - pgf(k,i,j) / gsqrt(k,i,j,i_uyz)
961  cf_t(k,i,j,1) = cor(k,i,j)
962  ddiv_t(k,i,j,2) = div
963  endif
964 #endif
965  enddo
966  enddo
967  enddo
968 #ifdef DEBUG
969  k = iundef; i = iundef; j = iundef
970 #endif
971 
972 #ifndef NO_FCT_DYN
973  if ( flag_fct_momentum ) then
974 
975  call atmos_dyn_fvm_fluxz_uyz_ud1( qflx_lo(:,:,:,zdir), & ! (out)
976  momz, momx0, dens, & ! (in)
977  gsqrt(:,:,:,i_uyw), j33g, & ! (in)
978  num_diff(:,:,:,i_momx,zdir), & ! (in)
979  cdz, & ! (in)
980  iis-1, iie+1, jjs-1, jje+1 ) ! (in)
981 
982  ! note that x-index is added by -1
983  call atmos_dyn_fvm_fluxx_uyz_ud1( qflx_lo(:,:,:,xdir), & ! (out)
984  momx, momx0, dens, & ! (in)
985  gsqrt(:,:,:,i_xyz), mapf(:,:,:,i_uy), & ! (in)
986  num_diff(:,:,:,i_momx,xdir), & ! (in)
987  cdz, & ! (in)
988  iis-1, iie+1, jjs-1, jje+1 ) ! (in)
989 
990  call atmos_dyn_fvm_fluxy_uyz_ud1( qflx_lo(:,:,:,ydir), & ! (out)
991  momy, momx0, dens, & ! (in)
992  gsqrt(:,:,:,i_uvz), mapf(:,:,:,i_xv), & ! (in)
993  num_diff(:,:,:,i_momx,ydir), & ! (in)
994  cdz, & ! (in)
995  iis-1, iie+1, jjs-1, jje+1 ) ! (in)
996 
997  endif
998  enddo
999  enddo
1000 
1001  if ( flag_fct_momentum ) then
1002 
1003  do j = js, je
1004  do i = is, ie
1005  do k = ks, ke
1006  qflx_hi(k,i,j,zdir) = qflx_hi(k,i,j,zdir) / ( mapf(i,j,1,i_uy) * mapf(i,j,2,i_uy) )&
1007  + qflx_j13(k,i,j) + qflx_j23(k,i,j)
1008  enddo
1009  enddo
1010  enddo
1011 
1012  do j = js-1, je+1
1013  do i = is-1, ie+1
1014  do k = ks, ke
1015  dens0_uvw(k,i,j) = 0.5_rp * ( dens0(k,i,j) + dens0(k,i+1,j) )
1016  dens_uvw(k,i,j) = 0.5_rp * ( dens_rk(k,i,j) + dens_rk(k,i+1,j) )
1017  enddo
1018  enddo
1019  enddo
1020 
1021  call comm_wait ( velx(:,:,:), 5 )
1022 
1023  call atmos_dyn_fct( qflx_anti, & ! (out)
1024  velx, dens0_uvw, dens_uvw, & ! (in)
1025  qflx_hi, qflx_lo, & ! (in)
1026  mflx_hi, & ! (in)
1027  rcdz, rfdx, rcdy, & ! (in)
1028  gsqrt(:,:,:,i_uyz), & ! (in)
1029  mapf(:,:,:,i_uy), dtrk, & ! (in)
1030  flag_fct_along_stream ) ! (in)
1031 
1032  do jjs = js, je, jblock
1033  jje = jjs+jblock-1
1034  do iis = is, ie, iblock
1035  iie = iis+iblock-1
1036 
1037  !--- update momentum(x)
1038  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
1039  do j = jjs, jje
1040  do i = iis, min(iie,ieh)
1041  do k = ks, ke
1042 #ifdef DEBUG
1043  call check( __line__, momx_rk(k,i,j) )
1044  call check( __line__, qflx_anti(k ,i ,j ,zdir) )
1045  call check( __line__, qflx_anti(k-1,i ,j ,zdir) )
1046  call check( __line__, qflx_anti(k ,i ,j ,xdir) )
1047  call check( __line__, qflx_anti(k ,i-1,j ,xdir) )
1048  call check( __line__, qflx_anti(k ,i ,j ,ydir) )
1049  call check( __line__, qflx_anti(k ,i ,j-1,ydir) )
1050 #endif
1051  momx_rk(k,i,j) = momx_rk(k,i,j) &
1052  + dtrk * ( ( ( qflx_anti(k,i,j,zdir) - qflx_anti(k-1,i ,j ,zdir) ) * rcdz(k) &
1053  + ( qflx_anti(k,i,j,xdir) - qflx_anti(k ,i-1,j ,xdir) ) * rfdx(i) &
1054  + ( qflx_anti(k,i,j,ydir) - qflx_anti(k ,i ,j-1,ydir) ) * rcdy(j) ) ) &
1055  * mapf(i,j,1,i_uy) * mapf(i,j,2,i_uy) &
1056  / gsqrt(k,i,j,i_uyz)
1057  enddo
1058  enddo
1059  enddo
1060 #ifdef DEBUG
1061  k = iundef; i = iundef; j = iundef
1062 #endif
1063 
1064  enddo
1065  enddo
1066 
1067 #ifdef DEBUG
1068  qflx_lo(:,:,:,:) = undef
1069  qflx_anti(:,:,:,:) = undef
1070 #endif
1071 
1072  endif ! FLAG_FCT_MOMENTUM
1073 
1074 #ifdef DEBUG
1075  qflx_hi(:,:,:,:) = undef
1076 #endif
1077 
1078  do jjs = js, je, jblock
1079  jje = jjs+jblock-1
1080  do iis = is, ie, iblock
1081  iie = iis+iblock-1
1082 #endif
1083 
1084  !########################################################################
1085  ! momentum equation (y)
1086  !########################################################################
1087 
1088  !-----< high order flux >-----
1089 
1090  ! at (x, v, w)
1091  call atmos_dyn_fvm_fluxz_xvz( qflx_hi(:,:,:,zdir), & ! (out)
1092  momz, momy, dens, & ! (in)
1093  gsqrt(:,:,:,i_xvw), j33g, & ! (in)
1094  num_diff(:,:,:,i_momy,zdir), & ! (in)
1095  cdz, & ! (in)
1096  iis, iie, jjs, jje ) ! (in)
1097  call atmos_dyn_fvm_fluxj13_xvz( qflx_j13, & ! (out)
1098  momx, momy, dens, & ! (in)
1099  gsqrt(:,:,:,i_xvz), j13g(:,:,:,i_xvw), mapf(:,:,:,i_xv), & ! (in)
1100  cdz, & ! (in)
1101  iis, iie, jjs, jje ) ! (in)
1102  call atmos_dyn_fvm_fluxj23_xvz( qflx_j23, & ! (out)
1103  momy, momy, dens, & ! (in)
1104  gsqrt(:,:,:,i_xvz), j23g(:,:,:,i_xvw), mapf(:,:,:,i_xv), & ! (in)
1105  cdz, & ! (in)
1106  iis, iie, jjs, jje ) ! (in)
1107 
1108  ! at (u, v, z)
1109  call atmos_dyn_fvm_fluxx_xvz( qflx_hi(:,:,:,xdir), & ! (out)
1110  momx, momy, dens, & ! (in)
1111  gsqrt(:,:,:,i_uvz), mapf(:,:,:,i_uv), & ! (in)
1112  num_diff(:,:,:,i_momy,xdir), & ! (in)
1113  cdz, & ! (in)
1114  iis, iie, jjs, jje ) ! (in)
1115 
1116  ! at (x, y, z)
1117  ! note that y-index is added by -1
1118  call atmos_dyn_fvm_fluxy_xvz( qflx_hi(:,:,:,ydir), & ! (out)
1119  momy, momy, dens, & ! (in)
1120  gsqrt(:,:,:,i_xyz), mapf(:,:,:,i_xy), & ! (in)
1121  num_diff(:,:,:,i_momy,ydir), & ! (in
1122  cdz, & ! (in)
1123  iis, iie, jjs, jje ) ! (in)
1124 
1125  ! pressure gradient force at (x, v, z)
1126 
1127  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
1128  do j = jjs, jje
1129  do i = iis, iie
1130  do k = ks, ke
1131  pgf(k,i,j) = ( ( gsqrt(k,i,j+1,i_xyz) * dpres(k,i,j+1) & ! [x,y,z]
1132  - gsqrt(k,i,j ,i_xyz) * dpres(k,i,j ) & ! [x,y,z]
1133  ) * rfdy(j) &
1134  + ( j23g(k ,i,j,i_xvw) &
1135  * 0.5_rp * ( f2h(k ,1,i_xvz) * ( dpres(k+1,i,j+1)+dpres(k+1,i,j) ) &
1136  + f2h(k ,2,i_xvz) * ( dpres(k ,i,j+1)+dpres(k ,i,j) ) ) & ! [x,y,z->x,v,w]
1137  - j23g(k-1,i,j,i_xvw) &
1138  * 0.5_rp * ( f2h(k-1,1,i_xvz) * ( dpres(k ,i,j+1)+dpres(k ,i,j) ) &
1139  + f2h(k-1,2,i_xvz) * ( dpres(k-1,i,j+1)+dpres(k-1,i,j) ) ) & ! [x,y,z->x,v,w]
1140  ) * rcdz(k) ) &
1141  * mapf(i,j,2,i_xv)
1142  enddo
1143  enddo
1144  enddo
1145 
1146  ! coriolis force at (x, v, z)
1147 
1148  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
1149  do j = jjs, jje
1150  do i = iis, iie
1151  do k = ks, ke
1152 #ifdef DEBUG
1153  call check( __line__, momx(k,i ,j ) )
1154  call check( __line__, momx(k,i ,j+1) )
1155  call check( __line__, momx(k,i-1,j ) )
1156  call check( __line__, momx(k,i-1,j+1) )
1157 #endif
1158  cor(k,i,j) = - 0.125_rp * ( corioli(1,i ,j+1)+corioli(1,i ,j) ) & ! [x,y,z->x,v,z]
1159  * ( momx(k,i ,j+1)+momx(k,i ,j) &
1160  + momx(k,i-1,j+1)+momx(k,i-1,j) ) & ! [u,y,z->x,v,z]
1161  - 0.25_rp * mapf(i,j,1,i_xv) * mapf(i,j,2,i_xv) &
1162  * ( momx(k,i,j) + momx(k,i-1,j) + momx(k,i,j+1) + momx(k,i-1,j+1) )&
1163  * ( momy(k,i,j) &
1164  * ( 1.0_rp/mapf(i,j,2,i_uv) - 1.0_rp/mapf(i-1,j,2,i_uv) ) * rcdx(i) &
1165  - 0.25_rp * ( momx(k,i,j)+momx(k,i-1,j)+momx(k,i,j+1)+momx(k,i-1,j+1) ) &
1166  * ( 1.0_rp/mapf(i,j+1,1,i_xy) - 1.0_rp/mapf(i,j,1,i_xy) ) * rfdy(j) ) &
1167  * 2.0_rp / ( dens(k,i,j) + dens(k,i,j+1) ) ! metoric term
1168  enddo
1169  enddo
1170  enddo
1171 
1172  !-----< update momentum (y) >-----
1173 
1174  !$omp parallel do private(i,j,k,advcv,advch,div) OMP_SCHEDULE_ collapse(2)
1175  do j = jjs, min(jje, jeh)
1176  do i = iis, iie
1177  do k = ks, ke
1178 #ifdef DEBUG
1179  call check( __line__, qflx_hi(k ,i ,j ,zdir) )
1180  call check( __line__, qflx_hi(k-1,i ,j ,zdir) )
1181  call check( __line__, qflx_hi(k ,i ,j ,xdir) )
1182  call check( __line__, qflx_hi(k ,i-1,j ,xdir) )
1183  call check( __line__, qflx_hi(k ,i ,j ,ydir) )
1184  call check( __line__, qflx_hi(k ,i ,j-1,ydir) )
1185  call check( __line__, ddiv(k,i,j+1) )
1186  call check( __line__, ddiv(k,i,j ) )
1187  call check( __line__, momy_t(k,i,j) )
1188  call check( __line__, momy0(k,i,j) )
1189 #endif
1190 
1191  advcv = - ( qflx_hi(k,i,j,zdir) - qflx_hi(k-1,i ,j ,zdir) &
1192  + qflx_j13(k,i,j) - qflx_j13(k-1,i,j) &
1193  + qflx_j23(k,i,j) - qflx_j23(k-1,i,j) ) * rcdz(k)
1194  advch = - ( ( qflx_hi(k,i,j,xdir) - qflx_hi(k ,i-1,j ,xdir) ) * rcdx(i) &
1195  + ( qflx_hi(k,i,j,ydir) - qflx_hi(k ,i ,j-1,ydir) ) * rfdy(j) ) &
1196  * mapf(i,j,1,i_xv) * mapf(i,j,2,i_xv)
1197  div = divdmp_coef / dtrk * fdy(j) * ( ddiv(k,i,j+1)-ddiv(k,i,j) )
1198  momy_rk(k,i,j) = momy0(k,i,j) &
1199  + dtrk * ( ( advcv + advch & ! advection
1200  - pgf(k,i,j) & ! pressure gradient force
1201  ) / gsqrt(k,i,j,i_xvz) &
1202  + cor(k,i,j) & ! coriolis force
1203  + div & ! divergence damping
1204  + momy_t(k,i,j) ) ! physics tendency
1205 #ifdef HIST_TEND
1206  if ( lhist ) then
1207  advcv_t(k,i,j,i_momy) = advcv / gsqrt(k,i,j,i_uyz)
1208  advch_t(k,i,j,i_momy) = advch / gsqrt(k,i,j,i_uyz)
1209  pg_t(k,i,j,3) = - pgf(k,i,j) / gsqrt(k,i,j,i_uyz)
1210  cf_t(k,i,j,2) = cor(k,i,j)
1211  ddiv_t(k,i,j,3) = div
1212  endif
1213 #endif
1214  enddo
1215  enddo
1216  enddo
1217 #ifdef DEBUG
1218  k = iundef; i = iundef; j = iundef
1219 #endif
1220 
1221 #ifndef NO_FCT_DYN
1222  if ( flag_fct_momentum ) then
1223 
1224  ! monotonic flux
1225  ! at (x, v, interface)
1226  call atmos_dyn_fvm_fluxz_xvz_ud1( qflx_lo(:,:,:,zdir), & ! (out)
1227  momz, momy0, dens, & ! (in)
1228  gsqrt(:,:,:,i_xvz), j33g, & ! (in)
1229  num_diff(:,:,:,i_momy,zdir), & ! (in)
1230  cdz, & ! (in)
1231  iis-1, iie+1, jjs-1, jje+1 ) ! (in)
1232 
1233  ! at (u, v, layer)
1234  call atmos_dyn_fvm_fluxx_xvz_ud1( qflx_lo(:,:,:,xdir), & ! (out)
1235  momx, momy0, dens, & ! (in)
1236  gsqrt(:,:,:,i_uvz), mapf(:,:,:,i_xy), & ! (in)
1237  num_diff(:,:,:,i_momy,xdir), & ! (in)
1238  cdz, & ! (in)
1239  iis-1, iie+1, jjs-1, jje+1 ) ! (in)
1240 
1241  ! at (x, y, layer)
1242  ! note that y-index is added by -1
1243  call atmos_dyn_fvm_fluxy_xvz_ud1( qflx_lo(:,:,:,ydir), & ! (out)
1244  momy, momy0, dens, & ! (in)
1245  gsqrt(:,:,:,i_xyz), mapf(:,:,:,i_xy), & ! (in)
1246  num_diff(:,:,:,i_momy,ydir), & ! (in)
1247  cdz, & ! (in)
1248  iis-1, iie+1, jjs-1, jje+1 ) ! (in)
1249 
1250  endif
1251  enddo
1252  enddo
1253 
1254  if ( flag_fct_momentum ) then
1255 
1256  do j = js, je
1257  do i = is, ie
1258  do k = ks, ke
1259  qflx_hi(k,i,j,zdir) = qflx_hi(k,i,j,zdir) / ( mapf(i,j,1,i_xv) * mapf(i,j,2,i_xv) ) &
1260  + qflx_j13(k,i,j) + qflx_j23(k,i,j)
1261  enddo
1262  enddo
1263  enddo
1264 
1265  do j = js-1, je+1
1266  do i = is-1, ie+1
1267  do k = ks, ke
1268  dens0_uvw(k,i,j) = 0.5_rp * ( dens0(k,i,j) + dens0(k,i,j+1) )
1269  dens_uvw(k,i,j) = 0.5_rp * ( dens_rk(k,i,j) + dens_rk(k,i,j+1) )
1270  enddo
1271  enddo
1272  enddo
1273 
1274  call comm_wait ( vely(:,:,:), 6 )
1275 
1276  call atmos_dyn_fct( qflx_anti, & ! (out)
1277  vely, dens0_uvw, dens_uvw, & ! (in)
1278  qflx_hi, qflx_lo, & ! (in)
1279  mflx_hi, & ! (in)
1280  rcdz, rcdx, rfdy, & ! (in)
1281  gsqrt(:,:,:,i_xvz), & ! (in)
1282  mapf(:,:,:,i_xv), dtrk, & ! (in)
1283  flag_fct_along_stream ) ! (in)
1284 
1285  do jjs = js, je, jblock
1286  jje = jjs+jblock-1
1287  do iis = is, ie, iblock
1288  iie = iis+iblock-1
1289 
1290  !--- update momentum(y)
1291  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
1292  do j = jjs, min(jje, jeh)
1293  do i = iis, iie
1294  do k = ks, ke
1295 #ifdef DEBUG
1296  call check( __line__, momy_rk(k,i,j) )
1297  call check( __line__, qflx_anti(k ,i ,j ,zdir) )
1298  call check( __line__, qflx_anti(k-1,i ,j ,zdir) )
1299  call check( __line__, qflx_anti(k ,i ,j ,xdir) )
1300  call check( __line__, qflx_anti(k ,i-1,j ,xdir) )
1301  call check( __line__, qflx_anti(k ,i ,j ,ydir) )
1302  call check( __line__, qflx_anti(k ,i ,j-1,ydir) )
1303 #endif
1304  momy_rk(k,i,j) = momy_rk(k,i,j) &
1305  + dtrk * ( ( ( qflx_anti(k,i,j,zdir) - qflx_anti(k-1,i ,j ,zdir) ) * rcdz(k) &
1306  + ( qflx_anti(k,i,j,xdir) - qflx_anti(k ,i-1,j ,xdir) ) * rcdx(i) &
1307  + ( qflx_anti(k,i,j,ydir) - qflx_anti(k ,i ,j-1,ydir) ) * rfdy(j) ) ) &
1308  * mapf(i,j,1,i_xv) * mapf(i,j,2,i_xv) &
1309  / gsqrt(k,i,j,i_xvz)
1310  enddo
1311  enddo
1312  enddo
1313 #ifdef DEBUG
1314  k = iundef; i = iundef; j = iundef
1315 #endif
1316 
1317  enddo
1318  enddo
1319 
1320 #ifdef DEBUG
1321  qflx_lo(:,:,:,:) = undef
1322  qflx_anti(:,:,:,:) = undef
1323 #endif
1324  endif ! FLAG_FCT_MOMENTUM
1325 
1326 #ifdef DEBUG
1327  qflx_hi(ks:,:,:,:) = undef
1328 #endif
1329 
1330  do jjs = js, je, jblock
1331  jje = jjs+jblock-1
1332  do iis = is, ie, iblock
1333  iie = iis+iblock-1
1334 #endif
1335 
1336  !########################################################################
1337  ! Thermodynamic equation
1338  !########################################################################
1339 
1340  !-----< high order flux >-----
1341 
1342  ! at (x, y, w)
1343  call atmos_dyn_fvm_fluxz_xyz( tflx_hi(:,:,:,zdir), & ! (out)
1344  mflx_hi(:,:,:,zdir), pott, gsqrt(:,:,:,i_xyw), & ! (in)
1345  num_diff(:,:,:,i_rhot,zdir), & ! (in)
1346  cdz, & ! (in)
1347  iis, iie, jjs, jje ) ! (in)
1348 
1349  ! at (u, y, z)
1350  call atmos_dyn_fvm_fluxx_xyz( tflx_hi(:,:,:,xdir), & ! (out)
1351  mflx_hi(:,:,:,xdir), pott, gsqrt(:,:,:,i_uyz), & ! (in)
1352  num_diff(:,:,:,i_rhot,xdir), & ! (in)
1353  cdz, & ! (in)
1354  iis, iie, jjs, jje ) ! (in)
1355 
1356  ! at (x, v, z)
1357  call atmos_dyn_fvm_fluxy_xyz( tflx_hi(:,:,:,ydir), & ! (out)
1358  mflx_hi(:,:,:,ydir), pott, gsqrt(:,:,:,i_xvz), & ! (in)
1359  num_diff(:,:,:,i_rhot,ydir), & ! (in)
1360  cdz, & ! (in)
1361  iis, iie, jjs, jje ) ! (in)
1362 
1363  !-----< update rho*theta >-----
1364 
1365  !$omp parallel do private(i,j,k,advcv,advch) OMP_SCHEDULE_ collapse(2)
1366  do j = jjs, jje
1367  do i = iis, iie
1368  do k = ks, ke
1369 #ifdef DEBUG
1370  call check( __line__, tflx_hi(k ,i ,j ,zdir) )
1371  call check( __line__, tflx_hi(k-1,i ,j ,zdir) )
1372  call check( __line__, tflx_hi(k ,i ,j ,xdir) )
1373  call check( __line__, tflx_hi(k ,i-1,j ,xdir) )
1374  call check( __line__, tflx_hi(k ,i ,j ,ydir) )
1375  call check( __line__, tflx_hi(k ,i ,j-1,ydir) )
1376  call check( __line__, rhot_t(k,i,j) )
1377  call check( __line__, rhot0(k,i,j) )
1378 #endif
1379  advcv = - ( tflx_hi(k,i,j,zdir) - tflx_hi(k-1,i ,j ,zdir) ) * rcdz(k)
1380  advch = - ( tflx_hi(k,i,j,xdir) - tflx_hi(k ,i-1,j ,xdir) ) * rcdx(i) &
1381  - ( tflx_hi(k,i,j,ydir) - tflx_hi(k ,i ,j-1,ydir) ) * rcdy(j)
1382  rhot_rk(k,i,j) = rhot0(k,i,j) &
1383  + dtrk * ( ( advcv + advch ) * mapf(i,j,1,i_xy) * mapf(i,j,2,i_xy) / gsqrt(k,i,j,i_xyz) &
1384  + rhot_t(k,i,j) )
1385 #ifdef HIST_TEND
1386  if ( lhist ) then
1387  advcv_t(k,i,j,i_rhot) = advcv * mapf(i,j,1,i_xy) * mapf(i,j,2,i_xy)/ gsqrt(k,i,j,i_xyz)
1388  advch_t(k,i,j,i_rhot) = advch * mapf(i,j,1,i_xy) * mapf(i,j,2,i_xy)/ gsqrt(k,i,j,i_xyz)
1389  endif
1390 #endif
1391  enddo
1392  enddo
1393  enddo
1394 #ifdef DEBUG
1395  k = iundef; i = iundef; j = iundef
1396 #endif
1397 
1398  enddo
1399  enddo
1400 
1401 #ifndef NO_FCT_DYN
1402 
1403  if ( flag_fct_t ) then
1404 
1405  call comm_vars8( mflx_hi(:,:,:,zdir), 1 )
1406  call comm_vars8( mflx_hi(:,:,:,xdir), 2 )
1407  call comm_vars8( mflx_hi(:,:,:,ydir), 3 )
1408  call comm_wait ( mflx_hi(:,:,:,zdir), 1, .false. )
1409  call comm_wait ( mflx_hi(:,:,:,xdir), 2, .false. )
1410  call comm_wait ( mflx_hi(:,:,:,ydir), 3, .false. )
1411 
1412  if ( .NOT. flag_fct_momentum ) then
1413  call comm_vars8( dens_rk, 1 )
1414  call comm_wait ( dens_rk, 1, .false. )
1415  endif
1416 
1417  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
1418  do j = js-1, je+2
1419  do i = is-1, ie+2
1420  do k = ks, ke
1421  pott(k,i,j) = rhot0(k,i,j) / dens0(k,i,j)
1422  enddo
1423  enddo
1424  enddo
1425 
1426  do jjs = js, je, jblock
1427  jje = jjs+jblock-1
1428  do iis = is, ie, iblock
1429  iie = iis+iblock-1
1430 
1431  ! monotonic flux
1432  ! at (x, y, interface)
1433  call atmos_dyn_fvm_fluxz_xyz_ud1( tflx_lo(:,:,:,zdir), & ! (out)
1434  mflx_hi(:,:,:,zdir), pott, gsqrt(:,:,:,i_xyz), & ! (in)
1435  num_diff(:,:,:,i_rhot,zdir), & ! (in)
1436  cdz, & ! (in)
1437  iis-1, iie+1, jjs-1, jje+1 ) ! (in)
1438 
1439  ! at (u, y, layer)
1440  call atmos_dyn_fvm_fluxx_xyz_ud1( tflx_lo(:,:,:,xdir), & ! (out)
1441  mflx_hi(:,:,:,xdir), pott, gsqrt(:,:,:,i_uyz), & ! (in)
1442  num_diff(:,:,:,i_rhot,xdir), & ! (in)
1443  cdz, & ! (in)
1444  iis-1, iie+1, jjs-1, jje+1 ) ! (in)
1445 
1446  ! at (x, v, layer)
1447  call atmos_dyn_fvm_fluxy_xyz_ud1( tflx_lo(:,:,:,ydir), & ! (out)
1448  mflx_hi(:,:,:,ydir), pott, gsqrt(:,:,:,i_xvz), & ! (in)
1449  num_diff(:,:,:,i_rhot,ydir), & ! (in)
1450  cdz, & ! (in)
1451  iis-1, iie+1, jjs-1, jje+1 ) ! (in)
1452 
1453  enddo
1454  enddo
1455 
1456  call atmos_dyn_fct( tflx_anti, & ! (out)
1457  pott, dens0, dens_rk, & ! (out)
1458  tflx_hi, tflx_lo, & ! (in)
1459  mflx_hi, & ! (in)
1460  rcdz, rcdx, rcdy, & ! (in)
1461  gsqrt(:,:,:,i_xyz), & ! (in)
1462  mapf(:,:,:,i_xy), dtrk, & ! (in)
1463  flag_fct_along_stream ) ! (in)
1464 
1465  do jjs = js, je, jblock
1466  jje = jjs+jblock-1
1467  do iis = is, ie, iblock
1468  iie = iis+iblock-1
1469 
1470  !--- update rho*theta
1471  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
1472  do j = jjs, jje
1473  do i = iis, iie
1474  do k = ks, ke
1475 #ifdef DEBUG
1476  call check( __line__, rhot_rk(k,i,j) )
1477  call check( __line__, tflx_anti(k ,i ,j ,zdir) )
1478  call check( __line__, tflx_anti(k-1,i ,j ,zdir) )
1479  call check( __line__, tflx_anti(k ,i ,j ,xdir) )
1480  call check( __line__, tflx_anti(k ,i-1,j ,xdir) )
1481  call check( __line__, tflx_anti(k ,i ,j ,ydir) )
1482  call check( __line__, tflx_anti(k ,i ,j-1,ydir) )
1483 #endif
1484  rhot_rk(k,i,j) = rhot_rk(k,i,j) &
1485  + dtrk * ( ( tflx_anti(k,i,j,zdir) - tflx_anti(k-1,i ,j ,zdir) ) * rcdz(k) &
1486  + ( tflx_anti(k,i,j,xdir) - tflx_anti(k ,i-1,j ,xdir) ) * rcdx(i) &
1487  + ( tflx_anti(k,i,j,ydir) - tflx_anti(k ,i ,j-1,ydir) ) * rcdy(j) ) &
1488  * mapf(i,j,1,i_xy) * mapf(i,j,2,i_xy) &
1489  / gsqrt(k,i,j,i_xyz)
1490  enddo
1491  enddo
1492  enddo
1493 #ifdef DEBUG
1494  k = iundef; i = iundef; j = iundef
1495 #endif
1496  enddo
1497  enddo
1498  endif
1499 #endif
1500 
1501 #ifdef HIST_TEND
1502  if ( lhist ) then
1503  call hist_in(advcv_t(:,:,:,i_dens), 'DENS_t_advcv', 'tendency of density (vert. advection)', 'kg/m3/s' )
1504  call hist_in(advcv_t(:,:,:,i_momz), 'MOMZ_t_advcv', 'tendency of momentum z (vert. advection)', 'kg/m2/s2', zdim='half')
1505  call hist_in(advcv_t(:,:,:,i_momx), 'MOMX_t_advcv', 'tendency of momentum x (vert. advection)', 'kg/m2/s2', xdim='half')
1506  call hist_in(advcv_t(:,:,:,i_momy), 'MOMY_t_advcv', 'tendency of momentum y (vert. advection)', 'kg/m2/s2', ydim='half')
1507  call hist_in(advcv_t(:,:,:,i_rhot), 'RHOT_t_advcv', 'tendency of rho*theta (vert. advection)', 'K kg/m3/s' )
1508 
1509  call hist_in(advch_t(:,:,:,i_dens), 'DENS_t_advch', 'tendency of density (horiz. advection)', 'kg/m3/s' )
1510  call hist_in(advch_t(:,:,:,i_momz), 'MOMZ_t_advch', 'tendency of momentum z (horiz. advection)', 'kg/m2/s2', zdim='half')
1511  call hist_in(advch_t(:,:,:,i_momx), 'MOMX_t_advch', 'tendency of momentum x (horiz. advection)', 'kg/m2/s2', xdim='half')
1512  call hist_in(advch_t(:,:,:,i_momy), 'MOMY_t_advch', 'tendency of momentum y (horiz. advection)', 'kg/m2/s2', ydim='half')
1513  call hist_in(advch_t(:,:,:,i_rhot), 'RHOT_t_advch', 'tendency of rho*theta (horiz. advection)', 'K kg/m3/s' )
1514 
1515  call hist_in(pg_t(:,:,:,1), 'MOMZ_t_pg', 'tendency of momentum z (pressure gradient)', 'kg/m2/s2', zdim='half')
1516  call hist_in(pg_t(:,:,:,2), 'MOMX_t_pg', 'tendency of momentum x (pressure gradient)', 'kg/m2/s2', xdim='half')
1517  call hist_in(pg_t(:,:,:,3), 'MOMY_t_pg', 'tendency of momentum y (pressure gradient)', 'kg/m2/s2', ydim='half')
1518 
1519  call hist_in(wdmp_t(:,:,:), 'MOMZ_t_wdamp', 'tendency of momentum z (Raileight damping)', 'kg/m2/s2', zdim='half')
1520 
1521  call hist_in(ddiv_t(:,:,:,1), 'MOMZ_t_ddiv', 'tendency of momentum z (divergence damping)', 'kg/m2/s2', zdim='half')
1522  call hist_in(ddiv_t(:,:,:,2), 'MOMX_t_ddiv', 'tendency of momentum x (divergence damping)', 'kg/m2/s2', xdim='half')
1523  call hist_in(ddiv_t(:,:,:,3), 'MOMY_t_ddiv', 'tendency of momentum y (divergence damping)', 'kg/m2/s2', ydim='half')
1524 
1525  call hist_in(cf_t(:,:,:,1), 'MOMX_t_cf', 'tendency of momentum x (coliolis force)', 'kg/m2/s2', xdim='half')
1526  call hist_in(cf_t(:,:,:,2), 'MOMY_t_cf', 'tendency of momentum y (coliolis force)', 'kg/m2/s2', ydim='half')
1527  endif
1528 #endif
1529 
1530  return
1531  end subroutine atmos_dyn_tstep_short_fvm_heve
1532 
1533 end module scale_atmos_dyn_tstep_short_fvm_heve
scale_index::i_rhot
integer, parameter, public i_rhot
Definition: scale_index.F90:35
scale_atmos_dyn_fvm_flux_ud1::atmos_dyn_fvm_fluxx_xyw_ud1
subroutine, public atmos_dyn_fvm_fluxx_xyw_ud1(flux, mom, val, DENS, GSQRT, MAPF, num_diff, CDZ, IIS, IIE, JJS, JJE)
calculation X-flux at XYW
Definition: scale_atmos_dyn_fvm_flux_ud1.F90:522
scale_grid_index::is
integer, public is
start point of inner domain: x, local
Definition: scale_grid_index.F90:60
scale_debug
module DEBUG
Definition: scale_debug.F90:13
scale_gridtrans::i_xvz
integer, public i_xvz
Definition: scale_gridtrans.F90:54
scale_grid_index::je
integer, public je
end point of inner domain: y, local
Definition: scale_grid_index.F90:63
scale_atmos_dyn_fvm_flux_ud1::atmos_dyn_fvm_fluxy_xyz_ud1
subroutine, public atmos_dyn_fvm_fluxy_xyz_ud1(flux, mflx, val, GSQRT, num_diff, CDZ, IIS, IIE, JJS, JJE)
calculation Y-flux at XYZ
Definition: scale_atmos_dyn_fvm_flux_ud1.F90:262
scale_atmos_dyn_fvm_flux_ud1::atmos_dyn_fvm_fluxy_xyw_ud1
subroutine, public atmos_dyn_fvm_fluxy_xyw_ud1(flux, mom, val, DENS, GSQRT, MAPF, num_diff, CDZ, IIS, IIE, JJS, JJE)
calculation Y-flux at XYW
Definition: scale_atmos_dyn_fvm_flux_ud1.F90:594
scale_process::prc_mpistop
subroutine, public prc_mpistop
Abort MPI.
Definition: scale_process.F90:267
scale_index::va
integer, public va
Definition: scale_index.F90:38
scale_atmos_dyn_fvm_flux_ud1::atmos_dyn_fvm_fluxz_xyz_ud1
subroutine, public atmos_dyn_fvm_fluxz_xyz_ud1(flux, mflx, val, GSQRT, num_diff, CDZ, IIS, IIE, JJS, JJE)
calculation z-flux at XYZ
Definition: scale_atmos_dyn_fvm_flux_ud1.F90:149
scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxx_xvz
procedure(flux_mom), pointer, public atmos_dyn_fvm_fluxx_xvz
Definition: scale_atmos_dyn_fvm_flux.F90:238
scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxj23_uyz
procedure(flux_j), pointer, public atmos_dyn_fvm_fluxj23_uyz
Definition: scale_atmos_dyn_fvm_flux.F90:227
scale_grid_index::iblock
integer, public iblock
block size for cache blocking: x
Definition: scale_grid_index.F90:47
scale_index::i_momx
integer, parameter, public i_momx
Definition: scale_index.F90:33
scale_stdio::io_l
logical, public io_l
output log or not? (this process)
Definition: scale_stdio.F90:61
scale_grid_index::zdir
integer, parameter, public zdir
Definition: scale_grid_index.F90:35
scale_atmos_dyn_fvm_flux_ud1::atmos_dyn_fvm_fluxy_xvz_ud1
subroutine, public atmos_dyn_fvm_fluxy_xvz_ud1(flux, mom, val, DENS, GSQRT, MAPF, num_diff, CDZ, IIS, IIE, JJS, JJE)
calculation Y-flux at XV
Definition: scale_atmos_dyn_fvm_flux_ud1.F90:1236
scale_index::i_momz
integer, parameter, public i_momz
Definition: scale_index.F90:32
scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxj23_xyw
procedure(flux_j), pointer, public atmos_dyn_fvm_fluxj23_xyw
Definition: scale_atmos_dyn_fvm_flux.F90:200
scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxx_xyz
procedure(flux_phi), pointer, public atmos_dyn_fvm_fluxx_xyz
Definition: scale_atmos_dyn_fvm_flux.F90:167
scale_stdio
module STDIO
Definition: scale_stdio.F90:12
scale_grid_index::ydir
integer, parameter, public ydir
Definition: scale_grid_index.F90:37
scale_grid_index::ke
integer, public ke
end point of inner domain: z, local
Definition: scale_grid_index.F90:59
scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxj13_xyw
procedure(flux_j), pointer, public atmos_dyn_fvm_fluxj13_xyw
Definition: scale_atmos_dyn_fvm_flux.F90:195
scale_grid_index::xdir
integer, parameter, public xdir
Definition: scale_grid_index.F90:36
scale_gridtrans::i_xy
integer, public i_xy
Definition: scale_gridtrans.F90:57
scale_atmos_dyn_tstep_short_fvm_heve::atmos_dyn_tstep_short_fvm_heve
subroutine, public atmos_dyn_tstep_short_fvm_heve(DENS_RK, MOMZ_RK, MOMX_RK, MOMY_RK, RHOT_RK, PROG_RK, mflx_hi, tflx_hi, DENS0, MOMZ0, MOMX0, MOMY0, RHOT0, DENS, MOMZ, MOMX, MOMY, RHOT, DENS_t, MOMZ_t, MOMX_t, MOMY_t, RHOT_t, PROG0, PROG, DPRES0, RT2P, CORIOLI, num_diff, wdamp_coef, divdmp_coef, DDIV, FLAG_FCT_MOMENTUM, FLAG_FCT_T, FLAG_FCT_ALONG_STREAM, CDZ, FDZ, FDX, FDY, RCDZ, RCDX, RCDY, RFDZ, RFDX, RFDY, PHI, GSQRT, J13G, J23G, J33G, MAPF, REF_dens, REF_rhot, BND_W, BND_E, BND_S, BND_N, dtrk, last)
Definition: scale_atmos_dyn_tstep_short_fvm_heve.F90:130
scale_debug::check
subroutine, public check(current_line, v)
Undefined value checker.
Definition: scale_debug.F90:58
scale_index::i_dens
integer, parameter, public i_dens
Definition: scale_index.F90:31
scale_index::i_momy
integer, parameter, public i_momy
Definition: scale_index.F90:34
scale_gridtrans::i_xvw
integer, public i_xvw
Definition: scale_gridtrans.F90:52
scale_const::const_undef
real(rp), public const_undef
Definition: scale_const.F90:43
scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxx_xyw
procedure(flux_mom), pointer, public atmos_dyn_fvm_fluxx_xyw
Definition: scale_atmos_dyn_fvm_flux.F90:184
scale_grid_index
module grid index
Definition: scale_grid_index.F90:14
scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxy_xyw
procedure(flux_mom), pointer, public atmos_dyn_fvm_fluxy_xyw
Definition: scale_atmos_dyn_fvm_flux.F90:189
scale_tracer
module TRACER
Definition: scale_tracer.F90:15
scale_index
module Index
Definition: scale_index.F90:14
scale_atmos_dyn_fvm_flux_ud1::atmos_dyn_fvm_fluxx_uyz_ud1
subroutine, public atmos_dyn_fvm_fluxx_uyz_ud1(flux, mom, val, DENS, GSQRT, MAPF, num_diff, CDZ, IIS, IIE, JJS, JJE)
calculation X-flux at UY
Definition: scale_atmos_dyn_fvm_flux_ud1.F90:870
scale_grid_index::ia
integer, public ia
of whole cells: x, local, with HALO
Definition: scale_grid_index.F90:55
scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxz_xvz
procedure(flux_z), pointer, public atmos_dyn_fvm_fluxz_xvz
Definition: scale_atmos_dyn_fvm_flux.F90:233
scale_gridtrans
module GRIDTRANS
Definition: scale_gridtrans.F90:15
scale_atmos_dyn_tstep_short_fvm_heve
module Atmosphere / Dynamics RK
Definition: scale_atmos_dyn_tstep_short_fvm_heve.F90:15
scale_atmos_dyn_fvm_flux_ud1::atmos_dyn_fvm_fluxz_xvz_ud1
subroutine, public atmos_dyn_fvm_fluxz_xvz_ud1(flux, mom, val, DENS, GSQRT, J33G, num_diff, CDZ, IIS, IIE, JJS, JJE)
calculation z-flux at XV
Definition: scale_atmos_dyn_fvm_flux_ud1.F90:980
scale_grid_index::ka
integer, public ka
of whole cells: z, local, with HALO
Definition: scale_grid_index.F90:54
scale_gridtrans::i_uy
integer, public i_uy
Definition: scale_gridtrans.F90:58
scale_atmos_dyn_tstep_short_fvm_heve::atmos_dyn_tstep_short_fvm_heve_setup
subroutine, public atmos_dyn_tstep_short_fvm_heve_setup
Setup.
Definition: scale_atmos_dyn_tstep_short_fvm_heve.F90:108
scale_const::const_pre00
real(rp), public const_pre00
pressure reference [Pa]
Definition: scale_const.F90:90
scale_grid_index::jblock
integer, public jblock
block size for cache blocking: y
Definition: scale_grid_index.F90:48
scale_gridtrans::i_xyw
integer, public i_xyw
Definition: scale_gridtrans.F90:50
scale_grid_index::jhalo
integer, public jhalo
of halo cells: y
Definition: scale_grid_index.F90:52
scale_comm
module COMMUNICATION
Definition: scale_comm.F90:23
scale_atmos_dyn_fvm_flux_ud1::atmos_dyn_fvm_fluxz_uyz_ud1
subroutine, public atmos_dyn_fvm_fluxz_uyz_ud1(flux, mom, val, DENS, GSQRT, J33G, num_diff, CDZ, IIS, IIE, JJS, JJE)
calculation z-flux at UY
Definition: scale_atmos_dyn_fvm_flux_ud1.F90:667
scale_const::const_grav
real(rp), public const_grav
standard acceleration of gravity [m/s2]
Definition: scale_const.F90:48
scale_grid_index::js
integer, public js
start point of inner domain: y, local
Definition: scale_grid_index.F90:62
scale_const::const_undef2
integer, parameter, public const_undef2
undefined value (INT2)
Definition: scale_const.F90:40
scale_atmos_dyn_common
module Atmosphere / Dynamics common
Definition: scale_atmos_dyn_common.F90:15
scale_process
module PROCESS
Definition: scale_process.F90:14
scale_gridtrans::i_uyw
integer, public i_uyw
Definition: scale_gridtrans.F90:51
scale_atmos_dyn_fvm_flux_ud1::atmos_dyn_fvm_fluxy_uyz_ud1
subroutine, public atmos_dyn_fvm_fluxy_uyz_ud1(flux, mom, val, DENS, GSQRT, MAPF, num_diff, CDZ, IIS, IIE, JJS, JJE)
calculation Y-flux at UY
Definition: scale_atmos_dyn_fvm_flux_ud1.F90:925
scale_const
module CONSTANT
Definition: scale_const.F90:14
scale_atmos_dyn_fvm_flux_ud1::atmos_dyn_fvm_fluxx_xvz_ud1
subroutine, public atmos_dyn_fvm_fluxx_xvz_ud1(flux, mom, val, DENS, GSQRT, MAPF, num_diff, CDZ, IIS, IIE, JJS, JJE)
calculation X-flux at XV
Definition: scale_atmos_dyn_fvm_flux_ud1.F90:1183
scale_grid_index::ks
integer, public ks
start point of inner domain: z, local
Definition: scale_grid_index.F90:58
scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxj13_xvz
procedure(flux_j), pointer, public atmos_dyn_fvm_fluxj13_xvz
Definition: scale_atmos_dyn_fvm_flux.F90:249
scale_grid_index::jeh
integer, public jeh
end point of inner domain: y, local (half level)
Definition: scale_grid_index.F90:76
scale_grid_index::ieh
integer, public ieh
end point of inner domain: x, local (half level)
Definition: scale_grid_index.F90:75
scale_gridtrans::i_xv
integer, public i_xv
Definition: scale_gridtrans.F90:59
scale_prof
module profiler
Definition: scale_prof.F90:10
scale_grid_index::ie
integer, public ie
end point of inner domain: x, local
Definition: scale_grid_index.F90:61
scale_const::const_eps
real(rp), public const_eps
small number
Definition: scale_const.F90:36
scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxz_uyz
procedure(flux_z), pointer, public atmos_dyn_fvm_fluxz_uyz
Definition: scale_atmos_dyn_fvm_flux.F90:206
scale_atmos_dyn_tstep_short_fvm_heve::atmos_dyn_tstep_short_fvm_heve_regist
subroutine, public atmos_dyn_tstep_short_fvm_heve_regist(ATMOS_DYN_TYPE, VA_out, VAR_NAME, VAR_DESC, VAR_UNIT)
Register.
Definition: scale_atmos_dyn_tstep_short_fvm_heve.F90:75
scale_atmos_dyn_fvm_flux
module scale_atmos_dyn_fvm_flux
Definition: scale_atmos_dyn_fvm_flux.F90:15
scale_gridtrans::i_uyz
integer, public i_uyz
Definition: scale_gridtrans.F90:53
scale_precision
module PRECISION
Definition: scale_precision.F90:13
scale_atmos_dyn_common::atmos_dyn_fct
subroutine, public atmos_dyn_fct(qflx_anti, phi_in, DENS0, DENS, qflx_hi, qflx_lo, mflx_hi, rdz, rdx, rdy, GSQRT, MAPF, dt, flag_vect)
Flux Correction Transport Limiter.
Definition: scale_atmos_dyn_common.F90:1942
scale_history
module HISTORY
Definition: scale_history.F90:16
scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxj13_uyz
procedure(flux_j), pointer, public atmos_dyn_fvm_fluxj13_uyz
Definition: scale_atmos_dyn_fvm_flux.F90:222
scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxj23_xvz
procedure(flux_j), pointer, public atmos_dyn_fvm_fluxj23_xvz
Definition: scale_atmos_dyn_fvm_flux.F90:254
scale_gridtrans::i_uvz
integer, public i_uvz
Definition: scale_gridtrans.F90:55
scale_gridtrans::i_xyz
integer, public i_xyz
Definition: scale_gridtrans.F90:49
scale_stdio::io_fid_log
integer, public io_fid_log
Log file ID.
Definition: scale_stdio.F90:56
scale_gridtrans::i_uv
integer, public i_uv
Definition: scale_gridtrans.F90:60
scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxz_xyz
procedure(flux_phi), pointer, public atmos_dyn_fvm_fluxz_xyz
Definition: scale_atmos_dyn_fvm_flux.F90:162
scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxx_uyz
procedure(flux_mom), pointer, public atmos_dyn_fvm_fluxx_uyz
Definition: scale_atmos_dyn_fvm_flux.F90:211
scale_atmos_dyn_fvm_flux_ud1::atmos_dyn_fvm_fluxz_xyw_ud1
subroutine, public atmos_dyn_fvm_fluxz_xyw_ud1(flux, mom, val, DENS, GSQRT, J33G, num_diff, CDZ, FDZ, dtrk, IIS, IIE, JJS, JJE)
calculation z-flux at XYW
Definition: scale_atmos_dyn_fvm_flux_ud1.F90:313
scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxy_xyz
procedure(flux_phi), pointer, public atmos_dyn_fvm_fluxy_xyz
Definition: scale_atmos_dyn_fvm_flux.F90:172
scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxy_uyz
procedure(flux_mom), pointer, public atmos_dyn_fvm_fluxy_uyz
Definition: scale_atmos_dyn_fvm_flux.F90:216
scale_atmos_dyn_fvm_flux_ud1::atmos_dyn_fvm_fluxx_xyz_ud1
subroutine, public atmos_dyn_fvm_fluxx_xyz_ud1(flux, mflx, val, GSQRT, num_diff, CDZ, IIS, IIE, JJS, JJE)
calculation X-flux at XYZ
Definition: scale_atmos_dyn_fvm_flux_ud1.F90:214
scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxz_xyw
procedure(flux_wz), pointer, public atmos_dyn_fvm_fluxz_xyw
Definition: scale_atmos_dyn_fvm_flux.F90:179
scale_grid_index::ihalo
integer, public ihalo
of halo cells: x
Definition: scale_grid_index.F90:51
scale_atmos_dyn_fvm_flux_ud1
module scale_atmos_dyn_fvm_flux_ud1
Definition: scale_atmos_dyn_fvm_flux_ud1.F90:19
scale_atmos_dyn_fvm_flux::atmos_dyn_fvm_fluxy_xvz
procedure(flux_mom), pointer, public atmos_dyn_fvm_fluxy_xvz
Definition: scale_atmos_dyn_fvm_flux.F90:243
scale_grid_index::ja
integer, public ja
of whole cells: y, local, with HALO
Definition: scale_grid_index.F90:56
Generated by   doxygen 1.8.14