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