SCALE-RM
scale_atmos_dyn_fvm_flux_ud5.F90
Go to the documentation of this file.
1 
2 !-------------------------------------------------------------------------------
11 !-------------------------------------------------------------------------------
12 ! Warning: This file was generated from atmos-rm/dynamics/scale_atmos_dyn_fvm_flux_udcd.F90.erb.
13 ! Do not edit this file.
14 !-------------------------------------------------------------------------------
15 #include "scalelib.h"
16 module scale_atmos_dyn_fvm_flux_ud5
17  !-----------------------------------------------------------------------------
18  !
19  !++ used modules
20  !
21  use scale_precision
22  use scale_io
23  use scale_prof
24  use scale_atmos_grid_cartesc_index
25  use scale_index
26  use scale_tracer
27  use scale_prc
28 #ifdef DEBUG
29  use scale_debug, only: &
30  check
31  use scale_const, only: &
32  undef => const_undef, &
33  iundef => const_undef2
34 #endif
35  !-----------------------------------------------------------------------------
36  implicit none
37  private
38  !-----------------------------------------------------------------------------
39  !
40  !++ Public procedure
41  !
42  public :: atmos_dyn_fvm_flux_valuew_z_ud5
43 
44  public :: atmos_dyn_fvm_fluxz_xyz_ud5
45  public :: atmos_dyn_fvm_fluxx_xyz_ud5
46  public :: atmos_dyn_fvm_fluxy_xyz_ud5
47 
48  public :: atmos_dyn_fvm_fluxz_xyw_ud5
49  public :: atmos_dyn_fvm_fluxj13_xyw_ud5
50  public :: atmos_dyn_fvm_fluxj23_xyw_ud5
51  public :: atmos_dyn_fvm_fluxx_xyw_ud5
52  public :: atmos_dyn_fvm_fluxy_xyw_ud5
53 
54  public :: atmos_dyn_fvm_fluxz_uyz_ud5
55  public :: atmos_dyn_fvm_fluxj13_uyz_ud5
56  public :: atmos_dyn_fvm_fluxj23_uyz_ud5
57  public :: atmos_dyn_fvm_fluxx_uyz_ud5
58  public :: atmos_dyn_fvm_fluxy_uyz_ud5
59 
60  public :: atmos_dyn_fvm_fluxz_xvz_ud5
61  public :: atmos_dyn_fvm_fluxj13_xvz_ud5
62  public :: atmos_dyn_fvm_fluxj23_xvz_ud5
63  public :: atmos_dyn_fvm_fluxx_xvz_ud5
64  public :: atmos_dyn_fvm_fluxy_xvz_ud5
65 
66  !-----------------------------------------------------------------------------
67  !
68  !++ Public parameters & variables
69  !
70  !-----------------------------------------------------------------------------
71  !
72  !++ Private procedure
73  !
74 #if 1
75 #define F2H(k,p,q) (CDZ(k+p-1)*GSQRT(k+p-1,i,j)/(CDZ(k)*GSQRT(k,i,j)+CDZ(k+1)*GSQRT(k+1,i,j)))
76 #else
77 #define F2H(k,p,q) 0.5_RP
78 #endif
79  !-----------------------------------------------------------------------------
80  !
81  !++ Private parameters & variables
82  !
83 
84  real(RP), parameter :: f2 = 0.5_rp ! F2 is always used to calculate flux near boundary.
85 
86 
87  real(RP), parameter :: f31 = -1.0_rp/12.0_rp
88  real(RP), parameter :: f32 = 7.0_rp/12.0_rp
89  real(RP), parameter :: f33 = 3.0_rp/12.0_rp
90 
91 
92  real(RP), parameter :: f51 = 1.0_rp/60.0_rp
93  real(RP), parameter :: f52 = -8.0_rp/60.0_rp
94  real(RP), parameter :: f53 = 37.0_rp/60.0_rp
95  real(RP), parameter :: f54 = -5.0_rp/60.0_rp
96  real(RP), parameter :: f55 = 10.0_rp/60.0_rp
97 
98 
99 
100 contains
101  !-----------------------------------------------------------------------------
103 !OCL SERIAL
104  subroutine atmos_dyn_fvm_flux_valuew_z_ud5( &
105  valW, &
106  mflx, val, GSQRT, &
107  CDZ )
108  implicit none
109 
110  real(RP), intent(out) :: valW (ka)
111  real(RP), intent(in) :: mflx (ka)
112  real(RP), intent(in) :: val (ka)
113  real(RP), intent(in) :: GSQRT(ka)
114  real(RP), intent(in) :: CDZ (ka)
115 
116  integer :: k
117  !---------------------------------------------------------------------------
118 
119  do k = ks+2, ke-3
120 #ifdef DEBUG
121  call check( __line__, mflx(k) )
122 
123  call check( __line__, val(k) )
124  call check( __line__, val(k+1) )
125 
126  call check( __line__, val(k-1) )
127  call check( __line__, val(k+2) )
128 
129  call check( __line__, val(k-2) )
130  call check( __line__, val(k+3) )
131 
132 #endif
133  valw(k) = ( f51 * ( val(k+3)+val(k-2) ) &
134  + f52 * ( val(k+2)+val(k-1) ) &
135  + f53 * ( val(k+1)+val(k) ) ) &
136  - ( f51 * ( val(k+3)-val(k-2) ) &
137  + f54 * ( val(k+2)-val(k-1) ) &
138  + f55 * ( val(k+1)-val(k) ) ) * sign(1.0_rp,mflx(k))
139  enddo
140 #ifdef DEBUG
141  k = iundef
142 #endif
143 
144 #ifdef DEBUG
145 
146  call check( __line__, mflx(ks) )
147  call check( __line__, val(ks ) )
148  call check( __line__, val(ks+1) )
149  call check( __line__, mflx(ke-1) )
150  call check( __line__, val(ke ) )
151  call check( __line__, val(ke-1) )
152 
153  call check( __line__, mflx(ks+1) )
154  call check( __line__, val(ks+2 ) )
155  call check( __line__, val(ks+3) )
156  call check( __line__, mflx(ke-2) )
157  call check( __line__, val(ke-2 ) )
158  call check( __line__, val(ke-3) )
159 
160 #endif
161 
162  valw(ks) = f2 * ( val(ks+1)+val(ks) ) &
163  * ( 0.5_rp + sign(0.5_rp,mflx(ks)) ) &
164  + ( 2.0_rp * val(ks) + 5.0_rp * val(ks+1) - val(ks+2) ) / 6.0_rp &
165  * ( 0.5_rp - sign(0.5_rp,mflx(ks)) )
166  valw(ke-1) = ( 2.0_rp * val(ke) + 5.0_rp * val(ke-1) - val(ke-2) ) / 6.0_rp &
167  * ( 0.5_rp + sign(0.5_rp,mflx(ke-1)) ) &
168  + f2 * ( val(ke)+val(ke-1) ) &
169  * ( 0.5_rp - sign(0.5_rp,mflx(ke-1)) )
170 
171  valw(ks+1) = ( 2.0_rp * val(ks+2) + 5.0_rp * val(ks+1) - val(ks) ) / 6.0_rp &
172  * ( 0.5_rp + sign(0.5_rp,mflx(ks+1)) ) &
173  + ( - 3.0_rp * val(ks) &
174  + 27.0_rp * val(ks+1) &
175  + 47.0_rp * val(ks+2) &
176  - 13.0_rp * val(ks+3) &
177  + 2.0_rp * val(ks+4) ) / 60.0_rp &
178  * ( 0.5_rp - sign(0.5_rp,mflx(ks+1)) )
179  valw(ke-2) = ( - 3.0_rp * val(ke) &
180  + 27.0_rp * val(ke-1) &
181  + 47.0_rp * val(ke-2) &
182  - 13.0_rp * val(ke-3) &
183  + 2.0_rp * val(ke-4) ) / 60.0_rp &
184  * ( 0.5_rp + sign(0.5_rp,mflx(ke-2)) ) &
185  + ( 2.0_rp * val(ke-2) + 5.0_rp * val(ke-1) - val(ke) ) / 6.0_rp &
186  * ( 0.5_rp - sign(0.5_rp,mflx(ke-2)) )
187 
188 
189  return
190  end subroutine atmos_dyn_fvm_flux_valuew_z_ud5
191 
192  !-----------------------------------------------------------------------------
194  subroutine atmos_dyn_fvm_fluxz_xyz_ud5( &
195  flux, &
196  mflx, val, GSQRT, &
197  num_diff, &
198  CDZ, &
199  IIS, IIE, JJS, JJE )
200  use scale_const, only: &
201  eps => const_eps
202  implicit none
203 
204  real(RP), intent(inout) :: flux (ka,ia,ja)
205  real(RP), intent(in) :: mflx (ka,ia,ja)
206  real(RP), intent(in) :: val (ka,ia,ja)
207  real(RP), intent(in) :: GSQRT (ka,ia,ja)
208  real(RP), intent(in) :: num_diff(ka,ia,ja)
209  real(RP), intent(in) :: CDZ (ka)
210  integer, intent(in) :: IIS, IIE, JJS, JJE
211 
212  real(RP) :: vel
213  integer :: k, i, j
214  !---------------------------------------------------------------------------
215 
216  !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
217  !$omp private(vel) &
218  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mflx,val,flux,GSQRT,num_diff)
219  do j = jjs, jje
220  do i = iis, iie
221  do k = ks+2, ke-3
222 #ifdef DEBUG
223  call check( __line__, mflx(k,i,j) )
224 
225  call check( __line__, val(k,i,j) )
226  call check( __line__, val(k+1,i,j) )
227 
228  call check( __line__, val(k-1,i,j) )
229  call check( __line__, val(k+2,i,j) )
230 
231  call check( __line__, val(k-2,i,j) )
232  call check( __line__, val(k+3,i,j) )
233 
234 #endif
235  vel = mflx(k,i,j)
236  flux(k,i,j) = vel &
237  * ( ( f51 * ( val(k+3,i,j)+val(k-2,i,j) ) &
238  + f52 * ( val(k+2,i,j)+val(k-1,i,j) ) &
239  + f53 * ( val(k+1,i,j)+val(k,i,j) ) ) &
240  - ( f51 * ( val(k+3,i,j)-val(k-2,i,j) ) &
241  + f54 * ( val(k+2,i,j)-val(k-1,i,j) ) &
242  + f55 * ( val(k+1,i,j)-val(k,i,j) ) ) * sign(1.0_rp,vel) ) &
243  + gsqrt(k,i,j) * num_diff(k,i,j)
244  enddo
245  enddo
246  enddo
247 #ifdef DEBUG
248  k = iundef; i = iundef; j = iundef
249 #endif
250 
251  !$omp parallel do default(none) private(i,j) OMP_SCHEDULE_ collapse(2) &
252  !$omp private(vel) &
253  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mflx,val,flux,GSQRT,num_diff,EPS)
254  do j = jjs, jje
255  do i = iis, iie
256 #ifdef DEBUG
257 
258  call check( __line__, mflx(ks,i,j) )
259  call check( __line__, val(ks ,i,j) )
260  call check( __line__, val(ks+1,i,j) )
261  call check( __line__, mflx(ke-1,i,j) )
262  call check( __line__, val(ke ,i,j) )
263  call check( __line__, val(ke-1,i,j) )
264 
265  call check( __line__, mflx(ks+1,i,j) )
266  call check( __line__, val(ks+2 ,i,j) )
267  call check( __line__, val(ks+3,i,j) )
268  call check( __line__, mflx(ke-2,i,j) )
269  call check( __line__, val(ke-2 ,i,j) )
270  call check( __line__, val(ke-3,i,j) )
271 
272 #endif
273  flux(ks-1,i,j) = 0.0_rp
274 
275  vel = mflx(ks,i,j)
276  flux(ks,i,j) = vel &
277  * ( f2 * ( val(ks+1,i,j)+val(ks,i,j) ) &
278  * ( 0.5_rp + sign(0.5_rp,vel) ) &
279  + ( 2.0_rp * val(ks,i,j) + 5.0_rp * val(ks+1,i,j) - val(ks+2,i,j) ) / 6.0_rp &
280  * ( 0.5_rp - sign(0.5_rp,vel) ) ) &
281  + gsqrt(ks,i,j) * num_diff(ks,i,j)
282  vel = mflx(ke-1,i,j)
283  flux(ke-1,i,j) = vel &
284  * ( ( 2.0_rp * val(ke,i,j) + 5.0_rp * val(ke-1,i,j) - val(ke-2,i,j) ) / 6.0_rp &
285  * ( 0.5_rp + sign(0.5_rp,vel) ) &
286  + f2 * ( val(ke,i,j)+val(ke-1,i,j) ) &
287  * ( 0.5_rp - sign(0.5_rp,vel) ) ) &
288  + gsqrt(ke-1,i,j) * num_diff(ke-1,i,j)
289 
290  vel = mflx(ks+1,i,j)
291  flux(ks+1,i,j) = vel &
292  * ( ( 2.0_rp * val(ks+2,i,j) + 5.0_rp * val(ks+1,i,j) - val(ks,i,j) ) / 6.0_rp &
293  * ( 0.5_rp + sign(0.5_rp,vel) ) &
294  + ( - 3.0_rp * val(ks,i,j) &
295  + 27.0_rp * val(ks+1,i,j) &
296  + 47.0_rp * val(ks+2,i,j) &
297  - 13.0_rp * val(ks+3,i,j) &
298  + 2.0_rp * val(ks+4,i,j) ) / 60.0_rp &
299  * ( 0.5_rp - sign(0.5_rp,vel) ) ) &
300  + gsqrt(ks+1,i,j) * num_diff(ks+1,i,j)
301  vel = mflx(ke-2,i,j)
302  flux(ke-2,i,j) = vel &
303  * ( ( - 3.0_rp * val(ke,i,j) &
304  + 27.0_rp * val(ke-1,i,j) &
305  + 47.0_rp * val(ke-2,i,j) &
306  - 13.0_rp * val(ke-3,i,j) &
307  + 2.0_rp * val(ke-4,i,j) ) / 60.0_rp &
308  * ( 0.5_rp + sign(0.5_rp,vel) ) &
309  + ( 2.0_rp * val(ke-2,i,j) + 5.0_rp * val(ke-1,i,j) - val(ke,i,j) ) / 6.0_rp &
310  * ( 0.5_rp - sign(0.5_rp,vel) ) ) &
311  + gsqrt(ke-2,i,j) * num_diff(ke-2,i,j)
312 
313  flux(ke ,i,j) = 0.0_rp
314  enddo
315  enddo
316 #ifdef DEBUG
317  k = iundef; i = iundef; j = iundef
318 #endif
319 
320  return
321  end subroutine atmos_dyn_fvm_fluxz_xyz_ud5
322 
323  !-----------------------------------------------------------------------------
325  subroutine atmos_dyn_fvm_fluxx_xyz_ud5( &
326  flux, &
327  mflx, val, GSQRT, &
328  num_diff, &
329  CDZ, &
330  IIS, IIE, JJS, JJE )
331  implicit none
332 
333  real(RP), intent(inout) :: flux (ka,ia,ja)
334  real(RP), intent(in) :: mflx (ka,ia,ja)
335  real(RP), intent(in) :: val (ka,ia,ja)
336  real(RP), intent(in) :: GSQRT (ka,ia,ja)
337  real(RP), intent(in) :: num_diff(ka,ia,ja)
338  real(RP), intent(in) :: CDZ(ka)
339  integer, intent(in) :: IIS, IIE, JJS, JJE
340 
341  real(RP) :: vel
342  integer :: k, i, j
343  !---------------------------------------------------------------------------
344 
345  !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
346  !$omp private(vel) &
347  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mflx,val,flux,GSQRT,num_diff)
348  do j = jjs, jje
349  do i = iis-1, iie
350  do k = ks, ke
351 #ifdef DEBUG
352  call check( __line__, mflx(k,i,j) )
353 
354  call check( __line__, val(k,i,j) )
355  call check( __line__, val(k,i+1,j) )
356 
357  call check( __line__, val(k,i-1,j) )
358  call check( __line__, val(k,i+2,j) )
359 
360  call check( __line__, val(k,i-2,j) )
361  call check( __line__, val(k,i+3,j) )
362 
363 #endif
364  vel = mflx(k,i,j)
365  flux(k,i,j) = vel &
366  * ( ( f51 * ( val(k,i+3,j)+val(k,i-2,j) ) &
367  + f52 * ( val(k,i+2,j)+val(k,i-1,j) ) &
368  + f53 * ( val(k,i+1,j)+val(k,i,j) ) ) &
369  - ( f51 * ( val(k,i+3,j)-val(k,i-2,j) ) &
370  + f54 * ( val(k,i+2,j)-val(k,i-1,j) ) &
371  + f55 * ( val(k,i+1,j)-val(k,i,j) ) ) * sign(1.0_rp,vel) ) &
372  + gsqrt(k,i,j) * num_diff(k,i,j)
373  enddo
374  enddo
375  enddo
376 #ifdef DEBUG
377  k = iundef; i = iundef; j = iundef
378 #endif
379 
380  return
381  end subroutine atmos_dyn_fvm_fluxx_xyz_ud5
382 
383  !-----------------------------------------------------------------------------
385  subroutine atmos_dyn_fvm_fluxy_xyz_ud5( &
386  flux, &
387  mflx, val, GSQRT, &
388  num_diff, &
389  CDZ, &
390  IIS, IIE, JJS, JJE )
391  implicit none
392 
393  real(RP), intent(inout) :: flux (ka,ia,ja)
394  real(RP), intent(in) :: mflx (ka,ia,ja)
395  real(RP), intent(in) :: val (ka,ia,ja)
396  real(RP), intent(in) :: GSQRT (ka,ia,ja)
397  real(RP), intent(in) :: num_diff(ka,ia,ja)
398  real(RP), intent(in) :: CDZ(ka)
399  integer, intent(in) :: IIS, IIE, JJS, JJE
400 
401  real(RP) :: vel
402  integer :: k, i, j
403  !---------------------------------------------------------------------------
404 
405  !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
406  !$omp private(vel) &
407  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mflx,val,flux,GSQRT,num_diff)
408  do j = jjs-1, jje
409  do i = iis, iie
410  do k = ks, ke
411 #ifdef DEBUG
412  call check( __line__, mflx(k,i,j) )
413 
414  call check( __line__, val(k,i,j) )
415  call check( __line__, val(k,i,j+1) )
416 
417  call check( __line__, val(k,i,j-1) )
418  call check( __line__, val(k,i,j+2) )
419 
420  call check( __line__, val(k,i,j-2) )
421  call check( __line__, val(k,i,j+3) )
422 
423 #endif
424  vel = mflx(k,i,j)
425  flux(k,i,j) = vel &
426  * ( ( f51 * ( val(k,i,j+3)+val(k,i,j-2) ) &
427  + f52 * ( val(k,i,j+2)+val(k,i,j-1) ) &
428  + f53 * ( val(k,i,j+1)+val(k,i,j) ) ) &
429  - ( f51 * ( val(k,i,j+3)-val(k,i,j-2) ) &
430  + f54 * ( val(k,i,j+2)-val(k,i,j-1) ) &
431  + f55 * ( val(k,i,j+1)-val(k,i,j) ) ) * sign(1.0_rp,vel) ) &
432  + gsqrt(k,i,j) * num_diff(k,i,j)
433  enddo
434  enddo
435  enddo
436 #ifdef DEBUG
437  k = iundef; i = iundef; j = iundef
438 #endif
439 
440  return
441  end subroutine atmos_dyn_fvm_fluxy_xyz_ud5
442 
443 
444  !-----------------------------------------------------------------------------
446  subroutine atmos_dyn_fvm_fluxz_xyw_ud5( &
447  flux, &
448  mom, val, DENS, &
449  GSQRT, J33G, &
450  num_diff, &
451  CDZ, FDZ, &
452  dtrk, &
453  IIS, IIE, JJS, JJE )
454  implicit none
455 
456  real(RP), intent(inout) :: flux (ka,ia,ja)
457  real(RP), intent(in) :: mom (ka,ia,ja)
458  real(RP), intent(in) :: val (ka,ia,ja)
459  real(RP), intent(in) :: DENS (ka,ia,ja)
460  real(RP), intent(in) :: GSQRT (ka,ia,ja)
461  real(RP), intent(in) :: J33G
462  real(RP), intent(in) :: num_diff(ka,ia,ja)
463  real(RP), intent(in) :: CDZ (ka)
464  real(RP), intent(in) :: FDZ (ka-1)
465  real(RP), intent(in) :: dtrk
466  integer, intent(in) :: IIS, IIE, JJS, JJE
467 
468  real(RP) :: vel
469  integer :: k, i, j
470  !---------------------------------------------------------------------------
471 
472  ! note than z-index is added by -1
473 
474  !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
475  !$omp private(vel) &
476  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,flux,J33G,GSQRT,num_diff,DENS)
477  do j = jjs, jje
478  do i = iis, iie
479  do k = ks+3, ke-2
480 #ifdef DEBUG
481  call check( __line__, mom(k-1,i,j) )
482  call check( __line__, mom(k ,i,j) )
483 
484  call check( __line__, val(k-1,i,j) )
485  call check( __line__, val(k,i,j) )
486 
487  call check( __line__, val(k-2,i,j) )
488  call check( __line__, val(k+1,i,j) )
489 
490  call check( __line__, val(k-3,i,j) )
491  call check( __line__, val(k+2,i,j) )
492 
493 #endif
494  vel = ( 0.5_rp * ( mom(k-1,i,j) &
495  + mom(k,i,j) ) ) &
496  / dens(k,i,j)
497  flux(k-1,i,j) = j33g * vel &
498  * ( ( f51 * ( val(k+2,i,j)+val(k-3,i,j) ) &
499  + f52 * ( val(k+1,i,j)+val(k-2,i,j) ) &
500  + f53 * ( val(k,i,j)+val(k-1,i,j) ) ) &
501  - ( f51 * ( val(k+2,i,j)-val(k-3,i,j) ) &
502  + f54 * ( val(k+1,i,j)-val(k-2,i,j) ) &
503  + f55 * ( val(k,i,j)-val(k-1,i,j) ) ) * sign(1.0_rp,vel) ) &
504  + gsqrt(k,i,j) * num_diff(k,i,j)
505  enddo
506  enddo
507  enddo
508 #ifdef DEBUG
509  k = iundef; i = iundef; j = iundef
510 #endif
511 
512  !$omp parallel do default(none) private(i,j) OMP_SCHEDULE_ collapse(2) &
513  !$omp private(vel) &
514  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,J33G,GSQRT,num_diff,FDZ,dtrk)
515  do j = jjs, jje
516  do i = iis, iie
517 #ifdef DEBUG
518 
519  call check( __line__, val(ks ,i,j) )
520  call check( __line__, val(ks+1,i,j) )
521 
522 
523  call check( __line__, val(ke-2,i,j) )
524  call check( __line__, val(ke-1,i,j) )
525 
526 #endif
527  ! The boundary condition is qflx_hi + qflxJ13 + qfluxJ23 = 0 at KS.
528  ! The flux at KS can be non-zero.
529  ! To reduce calculations, all the fluxes are set to zero.
530  flux(ks-1,i,j) = 0.0_rp ! k = KS
531 
532  vel = ( 0.5_rp * ( mom(ks,i,j) &
533  + mom(ks+1,i,j) ) ) &
534  / dens(ks+1,i,j)
535  flux(ks,i,j) = j33g * vel &
536  * ( f2 * ( val(ks+1,i,j)+val(ks,i,j) ) &
537  * ( 0.5_rp + sign(0.5_rp,vel) ) &
538  + ( 2.0_rp * val(ks,i,j) + 5.0_rp * val(ks+1,i,j) - val(ks+2,i,j) ) / 6.0_rp &
539  * ( 0.5_rp - sign(0.5_rp,vel) ) ) &
540  + gsqrt(ks+1,i,j) * num_diff(ks+1,i,j) ! k = KS+1
541 
542  vel = ( 0.5_rp * ( mom(ks+1,i,j) &
543  + mom(ks+2,i,j) ) ) &
544  / dens(ks+2,i,j)
545  flux(ks+1,i,j) = j33g * vel &
546  * ( ( 2.0_rp * val(ks+2,i,j) + 5.0_rp * val(ks+1,i,j) - val(ks,i,j) ) / 6.0_rp &
547  * ( 0.5_rp + sign(0.5_rp,vel) ) &
548  + ( - 3.0_rp * val(ks,i,j) &
549  + 27.0_rp * val(ks+1,i,j) &
550  + 47.0_rp * val(ks+2,i,j) &
551  - 13.0_rp * val(ks+3,i,j) &
552  + 2.0_rp * val(ks+4,i,j) ) / 60.0_rp &
553  * ( 0.5_rp - sign(0.5_rp,vel) ) ) &
554  + gsqrt(ks+2,i,j) * num_diff(ks+2,i,j) ! k = KS+2
555 
556 
557 
558  vel = ( 0.5_rp * ( mom(ke-2,i,j) &
559  + mom(ke-1,i,j) ) ) &
560  / dens(ke-1,i,j)
561  flux(ke-2,i,j) = j33g * vel &
562  * ( ( - 3.0_rp * val(ke,i,j) &
563  + 27.0_rp * val(ke-1,i,j) &
564  + 47.0_rp * val(ke-2,i,j) &
565  - 13.0_rp * val(ke-3,i,j) &
566  + 2.0_rp * val(ke-4,i,j) ) / 60.0_rp &
567  * ( 0.5_rp + sign(0.5_rp,vel) ) &
568  + ( 2.0_rp * val(ke-2,i,j) + 5.0_rp * val(ke-1,i,j) - val(ke,i,j) ) / 6.0_rp &
569  * ( 0.5_rp - sign(0.5_rp,vel) ) ) &
570  + gsqrt(ke-1,i,j) * num_diff(ke-1,i,j) ! k = KE-1
571 
572  flux(ke-1,i,j) = 0.0_rp ! k = KE
573  flux(ke ,i,j) = 0.0_rp ! k = KE+1
574  enddo
575  enddo
576 
577  return
578  end subroutine atmos_dyn_fvm_fluxz_xyw_ud5
579 
580 
581  !-----------------------------------------------------------------------------
583  subroutine atmos_dyn_fvm_fluxj13_xyw_ud5( &
584  flux, &
585  mom, val, DENS, &
586  GSQRT, J13G, MAPF, &
587  CDZ, &
588  IIS, IIE, JJS, JJE )
589  implicit none
590 
591  real(RP), intent(inout) :: flux (ka,ia,ja)
592  real(RP), intent(in) :: mom (ka,ia,ja)
593  real(RP), intent(in) :: val (ka,ia,ja)
594  real(RP), intent(in) :: DENS (ka,ia,ja)
595  real(RP), intent(in) :: GSQRT (ka,ia,ja)
596  real(RP), intent(in) :: J13G (ka,ia,ja)
597  real(RP), intent(in) :: MAPF ( ia,ja,2)
598  real(RP), intent(in) :: CDZ (ka)
599  integer, intent(in) :: IIS, IIE, JJS, JJE
600 
601  real(RP) :: vel
602  integer :: k, i, j
603  !---------------------------------------------------------------------------
604 
605  !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
606  !$omp private(vel) &
607  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,J13G,MAPF)
608  do j = jjs, jje
609  do i = iis, iie
610  do k = ks+2, ke-2
611  vel = ( 0.5_rp * ( mom(k,i,j)+mom(k,i-1,j) ) ) &
612  / dens(k,i,j)
613  vel = vel * j13g(k,i,j)
614  flux(k-1,i,j) = vel / mapf(i,j,+2) &
615  * ( ( f51 * ( val(k+2,i,j)+val(k-3,i,j) ) &
616  + f52 * ( val(k+1,i,j)+val(k-2,i,j) ) &
617  + f53 * ( val(k,i,j)+val(k-1,i,j) ) ) &
618  - ( f51 * ( val(k+2,i,j)-val(k-3,i,j) ) &
619  + f54 * ( val(k+1,i,j)-val(k-2,i,j) ) &
620  + f55 * ( val(k,i,j)-val(k-1,i,j) ) ) * sign(1.0_rp,vel) )
621  enddo
622  enddo
623  enddo
624 
625  !$omp parallel do default(none) private(i,j) OMP_SCHEDULE_ collapse(2) &
626  !$omp private(vel) &
627  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,J13G,MAPF)
628  do j = jjs, jje
629  do i = iis, iie
630  ! The boundary condition is qflx_hi + qflxJ13 + qfluxJ23 = 0 at KS.
631  ! The flux at KS can be non-zero.
632  ! To reduce calculations, all the fluxes are set to zero.
633  flux(ks-1,i,j) = 0.0_rp ! k = KS
634 
635  ! physically incorrect but for numerical stability
636  vel = ( ( 0.5_rp * ( mom(ks+1,i,j)+mom(ks+1,i-1,j) ) ) / dens(ks+1,i,j) &
637  + ( 0.5_rp * ( mom(ks,i,j)+mom(ks,i-1,j) ) ) / dens(ks ,i,j) ) * 0.5_rp
638 ! vel = ( 0.5_RP * ( mom(KS+1,i,j)+mom(KS+1,i-1,j) ) ) &
639 ! / DENS(KS+1,i,j)
640  vel = vel * j13g(ks+1,i,j)
641  flux(ks,i,j) = vel / mapf(i,j,+2) &
642  * ( f2 * ( val(ks+1,i,j)+val(ks,i,j) ) &
643  * ( 0.5_rp + sign(0.5_rp,vel) ) &
644  + ( 2.0_rp * val(ks,i,j) + 5.0_rp * val(ks+1,i,j) - val(ks+2,i,j) ) / 6.0_rp &
645  * ( 0.5_rp - sign(0.5_rp,vel) ) ) ! k = KS+1
646 
647 
648  vel = ( 0.5_rp * ( mom(ke-1,i,j)+mom(ke-1,i-1,j) ) ) &
649  / dens(ke-1,i,j)
650  vel = vel * j13g(ke-1,i,j)
651  flux(ke-2,i,j) = vel / mapf(i,j,+2) &
652  * ( ( 2.0_rp * val(ke-1,i,j) + 5.0_rp * val(ke-2,i,j) - val(ke-3,i,j) ) / 6.0_rp &
653  * ( 0.5_rp + sign(0.5_rp,vel) ) &
654  + f2 * ( val(ke-1,i,j)+val(ke-2,i,j) ) &
655  * ( 0.5_rp - sign(0.5_rp,vel) ) ) ! k = KE-3
656 
657  flux(ke-1,i,j) = 0.0_rp
658  enddo
659  enddo
660 
661  return
662  end subroutine atmos_dyn_fvm_fluxj13_xyw_ud5
663 
664  !-----------------------------------------------------------------------------
666  subroutine atmos_dyn_fvm_fluxj23_xyw_ud5( &
667  flux, &
668  mom, val, DENS, &
669  GSQRT, J23G, MAPF, &
670  CDZ, &
671  IIS, IIE, JJS, JJE )
672  implicit none
673 
674  real(RP), intent(inout) :: flux (ka,ia,ja)
675  real(RP), intent(in) :: mom (ka,ia,ja)
676  real(RP), intent(in) :: val (ka,ia,ja)
677  real(RP), intent(in) :: DENS (ka,ia,ja)
678  real(RP), intent(in) :: GSQRT (ka,ia,ja)
679  real(RP), intent(in) :: J23G (ka,ia,ja)
680  real(RP), intent(in) :: MAPF ( ia,ja,2)
681  real(RP), intent(in) :: CDZ (ka)
682  integer, intent(in) :: IIS, IIE, JJS, JJE
683 
684  real(RP) :: vel
685  integer :: k, i, j
686  !---------------------------------------------------------------------------
687 
688  !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
689  !$omp private(vel) &
690  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,J23G,MAPF)
691  do j = jjs, jje
692  do i = iis, iie
693  do k = ks+2, ke-2
694  vel = ( 0.5_rp * ( mom(k,i,j)+mom(k,i,j-1) ) ) &
695  / dens(k,i,j)
696  vel = vel * j23g(k,i,j)
697  flux(k-1,i,j) = vel / mapf(i,j,+1) &
698  * ( ( f51 * ( val(k+2,i,j)+val(k-3,i,j) ) &
699  + f52 * ( val(k+1,i,j)+val(k-2,i,j) ) &
700  + f53 * ( val(k,i,j)+val(k-1,i,j) ) ) &
701  - ( f51 * ( val(k+2,i,j)-val(k-3,i,j) ) &
702  + f54 * ( val(k+1,i,j)-val(k-2,i,j) ) &
703  + f55 * ( val(k,i,j)-val(k-1,i,j) ) ) * sign(1.0_rp,vel) )
704  enddo
705  enddo
706  enddo
707 
708  !$omp parallel do default(none) private(i,j) OMP_SCHEDULE_ collapse(2) &
709  !$omp private(vel) &
710  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,J23G,MAPF)
711  do j = jjs, jje
712  do i = iis, iie
713  ! The boundary condition is qflx_hi + qflxJ13 + qfluxJ23 = 0 at KS.
714  ! The flux at KS can be non-zero.
715  ! To reduce calculations, all the fluxes are set to zero.
716  flux(ks-1,i,j) = 0.0_rp ! k = KS
717 
718  ! physically incorrect but for numerical stability
719  vel = ( ( 0.5_rp * ( mom(ks+1,i,j)+mom(ks+1,i,j-1) ) ) / dens(ks+1,i,j) &
720  + ( 0.5_rp * ( mom(ks,i,j)+mom(ks,i,j-1) ) ) / dens(ks ,i,j) ) * 0.5_rp
721 ! vel = ( 0.5_RP * ( mom(KS+1,i,j)+mom(KS+1,i,j-1) ) ) &
722 ! / DENS(KS+1,i,j)
723  vel = vel * j23g(ks+1,i,j)
724  flux(ks,i,j) = vel / mapf(i,j,+1) &
725  * ( f2 * ( val(ks+1,i,j)+val(ks,i,j) ) &
726  * ( 0.5_rp + sign(0.5_rp,vel) ) &
727  + ( 2.0_rp * val(ks,i,j) + 5.0_rp * val(ks+1,i,j) - val(ks+2,i,j) ) / 6.0_rp &
728  * ( 0.5_rp - sign(0.5_rp,vel) ) ) ! k = KS+1
729 
730 
731  vel = ( 0.5_rp * ( mom(ke-1,i,j)+mom(ke-1,i,j-1) ) ) &
732  / dens(ke-1,i,j)
733  vel = vel * j23g(ke-1,i,j)
734  flux(ke-2,i,j) = vel / mapf(i,j,+1) &
735  * ( ( 2.0_rp * val(ke-1,i,j) + 5.0_rp * val(ke-2,i,j) - val(ke-3,i,j) ) / 6.0_rp &
736  * ( 0.5_rp + sign(0.5_rp,vel) ) &
737  + f2 * ( val(ke-1,i,j)+val(ke-2,i,j) ) &
738  * ( 0.5_rp - sign(0.5_rp,vel) ) ) ! k = KE-3
739 
740  flux(ke-1,i,j) = 0.0_rp
741  enddo
742  enddo
743 
744  return
745  end subroutine atmos_dyn_fvm_fluxj23_xyw_ud5
746 
747 
748  !-----------------------------------------------------------------------------
750  subroutine atmos_dyn_fvm_fluxx_xyw_ud5( &
751  flux, &
752  mom, val, DENS, &
753  GSQRT, MAPF, &
754  num_diff, &
755  CDZ, &
756  IIS, IIE, JJS, JJE )
757  implicit none
758 
759  real(RP), intent(inout) :: flux (ka,ia,ja)
760  real(RP), intent(in) :: mom (ka,ia,ja)
761  real(RP), intent(in) :: val (ka,ia,ja)
762  real(RP), intent(in) :: DENS (ka,ia,ja)
763  real(RP), intent(in) :: GSQRT (ka,ia,ja)
764  real(RP), intent(in) :: MAPF ( ia,ja,2)
765  real(RP), intent(in) :: num_diff(ka,ia,ja)
766  real(RP), intent(in) :: CDZ (ka)
767  integer, intent(in) :: IIS, IIE, JJS, JJE
768 
769  real(RP) :: vel
770  integer :: k, i, j
771  !---------------------------------------------------------------------------
772 
773  !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
774  !$omp private(vel) &
775  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,GSQRT,MAPF,num_diff) &
776  !$omp shared(CDZ)
777  do j = jjs, jje
778  do i = iis-1, iie
779  do k = ks, ke-1
780 #ifdef DEBUG
781  call check( __line__, mom(k ,i,j) )
782  call check( __line__, mom(k+1,i,j) )
783 
784  call check( __line__, val(k,i,j) )
785  call check( __line__, val(k,i+1,j) )
786 
787  call check( __line__, val(k,i-1,j) )
788  call check( __line__, val(k,i+2,j) )
789 
790  call check( __line__, val(k,i-2,j) )
791  call check( __line__, val(k,i+3,j) )
792 
793 #endif
794  vel = ( f2h(k,1,i_uyz) &
795  * mom(k+1,i,j) &
796  + f2h(k,2,i_uyz) &
797  * mom(k,i,j) ) &
798  / ( f2h(k,1,i_uyz) &
799  * 0.5_rp * ( dens(k+1,i,j)+dens(k+1,i+1,j) ) &
800  + f2h(k,2,i_uyz) &
801  * 0.5_rp * ( dens(k,i,j)+dens(k,i+1,j) ) )
802  flux(k,i,j) = gsqrt(k,i,j) / mapf(i,j,+2) * vel &
803  * ( ( f51 * ( val(k,i+3,j)+val(k,i-2,j) ) &
804  + f52 * ( val(k,i+2,j)+val(k,i-1,j) ) &
805  + f53 * ( val(k,i+1,j)+val(k,i,j) ) ) &
806  - ( f51 * ( val(k,i+3,j)-val(k,i-2,j) ) &
807  + f54 * ( val(k,i+2,j)-val(k,i-1,j) ) &
808  + f55 * ( val(k,i+1,j)-val(k,i,j) ) ) * sign(1.0_rp,vel) ) &
809  + gsqrt(k,i,j) * num_diff(k,i,j)
810  enddo
811  enddo
812  enddo
813 #ifdef DEBUG
814  k = iundef; i = iundef; j = iundef
815 #endif
816 
817  !$omp parallel do default(none) private(i,j) OMP_SCHEDULE_ collapse(2) &
818  !$omp private(vel) &
819  !$omp shared(JJS,JJE,IIS,IIE,KE,flux)
820  do j = jjs, jje
821  do i = iis-1, iie
822  flux(ke,i,j) = 0.0_rp
823  enddo
824  enddo
825 #ifdef DEBUG
826  k = iundef; i = iundef; j = iundef
827 #endif
828 
829  return
830  end subroutine atmos_dyn_fvm_fluxx_xyw_ud5
831 
832  !-----------------------------------------------------------------------------
834  subroutine atmos_dyn_fvm_fluxy_xyw_ud5( &
835  flux, &
836  mom, val, DENS, &
837  GSQRT, MAPF, &
838  num_diff, &
839  CDZ, &
840  IIS, IIE, JJS, JJE )
841  implicit none
842 
843  real(RP), intent(inout) :: flux (ka,ia,ja)
844  real(RP), intent(in) :: mom (ka,ia,ja)
845  real(RP), intent(in) :: val (ka,ia,ja)
846  real(RP), intent(in) :: DENS (ka,ia,ja)
847  real(RP), intent(in) :: GSQRT (ka,ia,ja)
848  real(RP), intent(in) :: MAPF ( ia,ja,2)
849  real(RP), intent(in) :: num_diff(ka,ia,ja)
850  real(RP), intent(in) :: CDZ (ka)
851  integer, intent(in) :: IIS, IIE, JJS, JJE
852 
853  real(RP) :: vel
854  integer :: k, i, j
855  !---------------------------------------------------------------------------
856 
857  !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
858  !$omp private(vel) &
859  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,GSQRT,MAPF,num_diff) &
860  !$omp shared(CDZ)
861  do j = jjs-1, jje
862  do i = iis, iie
863  do k = ks, ke-1
864 #ifdef DEBUG
865  call check( __line__, mom(k ,i,j) )
866  call check( __line__, mom(k+1,i,j) )
867 
868  call check( __line__, val(k,i,j) )
869  call check( __line__, val(k,i,j+1) )
870 
871  call check( __line__, val(k,i,j-1) )
872  call check( __line__, val(k,i,j+2) )
873 
874  call check( __line__, val(k,i,j-2) )
875  call check( __line__, val(k,i,j+3) )
876 
877 #endif
878  vel = ( f2h(k,1,i_xvz) &
879  * mom(k+1,i,j) &
880  + f2h(k,2,i_xvz) &
881  * mom(k,i,j) ) &
882  / ( f2h(k,1,i_xvz) &
883  * 0.5_rp * ( dens(k+1,i,j)+dens(k+1,i,j+1) ) &
884  + f2h(k,2,i_xvz) &
885  * 0.5_rp * ( dens(k,i,j)+dens(k,i,j+1) ) )
886  flux(k,i,j) = gsqrt(k,i,j) / mapf(i,j,+1) * vel &
887  * ( ( f51 * ( val(k,i,j+3)+val(k,i,j-2) ) &
888  + f52 * ( val(k,i,j+2)+val(k,i,j-1) ) &
889  + f53 * ( val(k,i,j+1)+val(k,i,j) ) ) &
890  - ( f51 * ( val(k,i,j+3)-val(k,i,j-2) ) &
891  + f54 * ( val(k,i,j+2)-val(k,i,j-1) ) &
892  + f55 * ( val(k,i,j+1)-val(k,i,j) ) ) * sign(1.0_rp,vel) ) &
893  + gsqrt(k,i,j) * num_diff(k,i,j)
894  enddo
895  enddo
896  enddo
897 #ifdef DEBUG
898  k = iundef; i = iundef; j = iundef
899 #endif
900 
901  !$omp parallel do default(none) private(i,j) OMP_SCHEDULE_ collapse(2) &
902  !$omp private(vel) &
903  !$omp shared(JJS,JJE,IIS,IIE,KE,flux)
904  do j = jjs-1, jje
905  do i = iis, iie
906  flux(ke,i,j) = 0.0_rp
907  enddo
908  enddo
909 #ifdef DEBUG
910  k = iundef; i = iundef; j = iundef
911 #endif
912 
913  return
914  end subroutine atmos_dyn_fvm_fluxy_xyw_ud5
915 
916 
917  !-----------------------------------------------------------------------------
919  subroutine atmos_dyn_fvm_fluxz_uyz_ud5( &
920  flux, &
921  mom, val, DENS, &
922  GSQRT, J33G, &
923  num_diff, &
924  CDZ, &
925  IIS, IIE, JJS, JJE )
926  implicit none
927 
928  real(RP), intent(inout) :: flux (ka,ia,ja)
929  real(RP), intent(in) :: mom (ka,ia,ja)
930  real(RP), intent(in) :: val (ka,ia,ja)
931  real(RP), intent(in) :: DENS (ka,ia,ja)
932  real(RP), intent(in) :: GSQRT (ka,ia,ja)
933  real(RP), intent(in) :: J33G
934  real(RP), intent(in) :: num_diff(ka,ia,ja)
935  real(RP), intent(in) :: CDZ (ka)
936  integer, intent(in) :: IIS, IIE, JJS, JJE
937 
938  real(RP) :: vel
939  integer :: k, i, j
940  !---------------------------------------------------------------------------
941 
942  !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
943  !$omp private(vel) &
944  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,J33G,GSQRT,num_diff) &
945  !$omp shared(CDZ)
946  do j = jjs, jje
947  do i = iis, iie
948  do k = ks+2, ke-3
949 #ifdef DEBUG
950  call check( __line__, mom(k,i,j) )
951  call check( __line__, mom(k,i+1,j) )
952 
953  call check( __line__, val(k,i,j) )
954  call check( __line__, val(k+1,i,j) )
955 
956  call check( __line__, val(k-1,i,j) )
957  call check( __line__, val(k+2,i,j) )
958 
959  call check( __line__, val(k-2,i,j) )
960  call check( __line__, val(k+3,i,j) )
961 
962 #endif
963  vel = ( 0.5_rp * ( mom(k,i,j)+mom(k,i+1,j) ) ) &
964  / ( f2h(k,1,i_uyz) &
965  * 0.5_rp * ( dens(k+1,i,j)+dens(k+1,i+1,j) ) &
966  + f2h(k,2,i_uyz) &
967  * 0.5_rp * ( dens(k,i,j)+dens(k,i+1,j) ) )
968  flux(k,i,j) = j33g * vel &
969  * ( ( f51 * ( val(k+3,i,j)+val(k-2,i,j) ) &
970  + f52 * ( val(k+2,i,j)+val(k-1,i,j) ) &
971  + f53 * ( val(k+1,i,j)+val(k,i,j) ) ) &
972  - ( f51 * ( val(k+3,i,j)-val(k-2,i,j) ) &
973  + f54 * ( val(k+2,i,j)-val(k-1,i,j) ) &
974  + f55 * ( val(k+1,i,j)-val(k,i,j) ) ) * sign(1.0_rp,vel) ) &
975  + gsqrt(k,i,j) * num_diff(k,i,j)
976  enddo
977  enddo
978  enddo
979 #ifdef DEBUG
980  k = iundef; i = iundef; j = iundef
981 #endif
982 
983  !$omp parallel do default(none) private(i,j) OMP_SCHEDULE_ collapse(2) &
984  !$omp private(vel) &
985  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,flux,J33G,GSQRT,num_diff,DENS,CDZ)
986  do j = jjs, jje
987  do i = iis, iie
988 #ifdef DEBUG
989 
990  call check( __line__, mom(ks,i ,j) )
991  call check( __line__, mom(ks,i+1,j) )
992  call check( __line__, val(ks+1,i,j) )
993  call check( __line__, val(ks,i,j) )
994 
995  call check( __line__, mom(ks+1,i ,j) )
996  call check( __line__, mom(ks+1,i+1,j) )
997  call check( __line__, val(ks+3,i,j) )
998  call check( __line__, val(ks+2,i,j) )
999 
1000 #endif
1001  ! The boundary condition is qflx_hi + qflxJ13 + qfluxJ23 = 0 at KS-1.
1002  ! The flux at KS-1 can be non-zero.
1003  ! To reduce calculations, all the fluxes are set to zero.
1004  flux(ks-1,i,j) = 0.0_rp
1005 
1006  vel = ( 0.5_rp * ( mom(ks,i,j)+mom(ks,i+1,j) ) ) &
1007  / ( f2h(ks,1,i_uyz) &
1008  * 0.5_rp * ( dens(ks+1,i,j)+dens(ks+1,i+1,j) ) &
1009  + f2h(ks,2,i_uyz) &
1010  * 0.5_rp * ( dens(ks,i,j)+dens(ks,i+1,j) ) )
1011  flux(ks,i,j) = j33g * vel &
1012  * ( f2 * ( val(ks+1,i,j)+val(ks,i,j) ) &
1013  * ( 0.5_rp + sign(0.5_rp,vel) ) &
1014  + ( 2.0_rp * val(ks,i,j) + 5.0_rp * val(ks+1,i,j) - val(ks+2,i,j) ) / 6.0_rp &
1015  * ( 0.5_rp - sign(0.5_rp,vel) ) ) &
1016  + gsqrt(ks,i,j) * num_diff(ks,i,j)
1017  vel = ( 0.5_rp * ( mom(ke-1,i,j)+mom(ke-1,i+1,j) ) ) &
1018  / ( f2h(ke-1,1,i_uyz) &
1019  * 0.5_rp * ( dens(ke,i,j)+dens(ke,i+1,j) ) &
1020  + f2h(ke-1,2,i_uyz) &
1021  * 0.5_rp * ( dens(ke-1,i,j)+dens(ke-1,i+1,j) ) )
1022  flux(ke-1,i,j) = j33g * vel &
1023  * ( ( 2.0_rp * val(ke,i,j) + 5.0_rp * val(ke-1,i,j) - val(ke-2,i,j) ) / 6.0_rp &
1024  * ( 0.5_rp + sign(0.5_rp,vel) ) &
1025  + f2 * ( val(ke,i,j)+val(ke-1,i,j) ) &
1026  * ( 0.5_rp - sign(0.5_rp,vel) ) ) &
1027  + gsqrt(ke-1,i,j) * num_diff(ke-1,i,j)
1028 
1029  vel = ( 0.5_rp * ( mom(ks+1,i,j)+mom(ks+1,i+1,j) ) ) &
1030  / ( f2h(ks+1,1,i_uyz) &
1031  * 0.5_rp * ( dens(ks+2,i,j)+dens(ks+2,i+1,j) ) &
1032  + f2h(ks+1,2,i_uyz) &
1033  * 0.5_rp * ( dens(ks+1,i,j)+dens(ks+1,i+1,j) ) )
1034  flux(ks+1,i,j) = j33g * vel &
1035  * ( ( 2.0_rp * val(ks+2,i,j) + 5.0_rp * val(ks+1,i,j) - val(ks,i,j) ) / 6.0_rp &
1036  * ( 0.5_rp + sign(0.5_rp,vel) ) &
1037  + ( - 3.0_rp * val(ks,i,j) &
1038  + 27.0_rp * val(ks+1,i,j) &
1039  + 47.0_rp * val(ks+2,i,j) &
1040  - 13.0_rp * val(ks+3,i,j) &
1041  + 2.0_rp * val(ks+4,i,j) ) / 60.0_rp &
1042  * ( 0.5_rp - sign(0.5_rp,vel) ) ) &
1043  + gsqrt(ks+1,i,j) * num_diff(ks+1,i,j)
1044  vel = ( 0.5_rp * ( mom(ke-2,i,j)+mom(ke-2,i+1,j) ) ) &
1045  / ( f2h(ke-2,1,i_uyz) &
1046  * 0.5_rp * ( dens(ke-1,i,j)+dens(ke-1,i+1,j) ) &
1047  + f2h(ke-2,2,i_uyz) &
1048  * 0.5_rp * ( dens(ke-2,i,j)+dens(ke-2,i+1,j) ) )
1049  flux(ke-2,i,j) = j33g * vel &
1050  * ( ( - 3.0_rp * val(ke,i,j) &
1051  + 27.0_rp * val(ke-1,i,j) &
1052  + 47.0_rp * val(ke-2,i,j) &
1053  - 13.0_rp * val(ke-3,i,j) &
1054  + 2.0_rp * val(ke-4,i,j) ) / 60.0_rp &
1055  * ( 0.5_rp + sign(0.5_rp,vel) ) &
1056  + ( 2.0_rp * val(ke-2,i,j) + 5.0_rp * val(ke-1,i,j) - val(ke,i,j) ) / 6.0_rp &
1057  * ( 0.5_rp - sign(0.5_rp,vel) ) ) &
1058  + gsqrt(ke-2,i,j) * num_diff(ke-2,i,j)
1059 
1060  flux(ke,i,j) = 0.0_rp
1061  enddo
1062  enddo
1063 #ifdef DEBUG
1064  k = iundef; i = iundef; j = iundef
1065 #endif
1066 
1067  return
1068  end subroutine atmos_dyn_fvm_fluxz_uyz_ud5
1069 
1070  !-----------------------------------------------------------------------------
1072  subroutine atmos_dyn_fvm_fluxj13_uyz_ud5( &
1073  flux, &
1074  mom, val, DENS, &
1075  GSQRT, J13G, MAPF, &
1076  CDZ, &
1077  IIS, IIE, JJS, JJE )
1078  implicit none
1079 
1080  real(RP), intent(inout) :: flux (ka,ia,ja)
1081  real(RP), intent(in) :: mom (ka,ia,ja)
1082  real(RP), intent(in) :: val (ka,ia,ja)
1083  real(RP), intent(in) :: DENS (ka,ia,ja)
1084  real(RP), intent(in) :: GSQRT (ka,ia,ja)
1085  real(RP), intent(in) :: J13G (ka,ia,ja)
1086  real(RP), intent(in) :: MAPF ( ia,ja,2)
1087  real(RP), intent(in) :: CDZ (ka)
1088  integer, intent(in) :: IIS, IIE, JJS, JJE
1089 
1090  real(RP) :: vel
1091  integer :: k, i, j
1092  !---------------------------------------------------------------------------
1093 
1094  !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
1095  !$omp private(vel) &
1096  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,J13G,MAPF) &
1097  !$omp shared(GSQRT,CDZ)
1098  do j = jjs, jje
1099  do i = iis, iie
1100  do k = ks+2, ke-3
1101  vel = ( f2h(k,1,i_uyz) &
1102  * mom(k+1,i,j) &
1103  + f2h(k,2,i_uyz) &
1104  * mom(k,i,j) ) &
1105  / ( f2h(k,1,i_uyz) &
1106  * 0.5_rp * ( dens(k+1,i,j)+dens(k+1,i+1,j) ) &
1107  + f2h(k,2,i_uyz) &
1108  * 0.5_rp * ( dens(k,i,j)+dens(k,i+1,j) ) )
1109  vel = vel * j13g(k,i,j)
1110  flux(k,i,j) = vel / mapf(i,j,+2) &
1111  * ( ( f51 * ( val(k+3,i,j)+val(k-2,i,j) ) &
1112  + f52 * ( val(k+2,i,j)+val(k-1,i,j) ) &
1113  + f53 * ( val(k+1,i,j)+val(k,i,j) ) ) &
1114  - ( f51 * ( val(k+3,i,j)-val(k-2,i,j) ) &
1115  + f54 * ( val(k+2,i,j)-val(k-1,i,j) ) &
1116  + f55 * ( val(k+1,i,j)-val(k,i,j) ) ) * sign(1.0_rp,vel) )
1117  enddo
1118  enddo
1119  enddo
1120 
1121  !$omp parallel do default(none) private(i,j) OMP_SCHEDULE_ collapse(2) &
1122  !$omp private(vel) &
1123  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,J13G,MAPF) &
1124  !$omp shared(GSQRT,CDZ)
1125  do j = jjs, jje
1126  do i = iis, iie
1127  ! The boundary condition is qflx_hi + qflxJ13 + qfluxJ23 = 0 at KS-1.
1128  ! The flux at KS-1 can be non-zero.
1129  ! To reduce calculations, all the fluxes are set to zero.
1130  flux(ks-1,i,j) = 0.0_rp
1131 
1132  vel = ( f2h(ks,1,i_uyz) &
1133  * mom(ks+1,i,j) &
1134  + f2h(ks,2,i_uyz) &
1135  * mom(ks,i,j) ) &
1136  / ( f2h(ks,1,i_uyz) &
1137  * 0.5_rp * ( dens(ks+1,i,j)+dens(ks+1,i+1,j) ) &
1138  + f2h(ks,2,i_uyz) &
1139  * 0.5_rp * ( dens(ks,i,j)+dens(ks,i+1,j) ) )
1140  vel = vel * j13g(ks,i,j)
1141  flux(ks,i,j) = vel / mapf(i,j,+2) &
1142  * ( f2 * ( val(ks+1,i,j)+val(ks,i,j) ) &
1143  * ( 0.5_rp + sign(0.5_rp,vel) ) &
1144  + ( 2.0_rp * val(ks,i,j) + 5.0_rp * val(ks+1,i,j) - val(ks+2,i,j) ) / 6.0_rp &
1145  * ( 0.5_rp - sign(0.5_rp,vel) ) )
1146 
1147  vel = ( f2h(ke-1,1,i_uyz) &
1148  * mom(ke,i,j) &
1149  + f2h(ke-1,2,i_uyz) &
1150  * mom(ke-1,i,j) ) &
1151  / ( f2h(ke-1,1,i_uyz) &
1152  * 0.5_rp * ( dens(ke,i,j)+dens(ke,i+1,j) ) &
1153  + f2h(ke-1,2,i_uyz) &
1154  * 0.5_rp * ( dens(ke-1,i,j)+dens(ke-1,i+1,j) ) )
1155  vel = vel * j13g(ke-1,i,j)
1156  flux(ke-1,i,j) = vel / mapf(i,j,+2) &
1157  * ( ( 2.0_rp * val(ke,i,j) + 5.0_rp * val(ke-1,i,j) - val(ke-2,i,j) ) / 6.0_rp &
1158  * ( 0.5_rp + sign(0.5_rp,vel) ) &
1159  + f2 * ( val(ke,i,j)+val(ke-1,i,j) ) &
1160  * ( 0.5_rp - sign(0.5_rp,vel) ) )
1161 
1162  vel = ( f2h(ks+1,1,i_uyz) &
1163  * mom(ks+2,i,j) &
1164  + f2h(ks+1,2,i_uyz) &
1165  * mom(ks+1,i,j) ) &
1166  / ( f2h(ks+1,1,i_uyz) &
1167  * 0.5_rp * ( dens(ks+2,i,j)+dens(ks+2,i+1,j) ) &
1168  + f2h(ks+1,2,i_uyz) &
1169  * 0.5_rp * ( dens(ks+1,i,j)+dens(ks+1,i+1,j) ) )
1170  vel = vel * j13g(ks+1,i,j)
1171  flux(ks+1,i,j) = vel / mapf(i,j,+2) &
1172  * ( ( 2.0_rp * val(ks+2,i,j) + 5.0_rp * val(ks+1,i,j) - val(ks,i,j) ) / 6.0_rp &
1173  * ( 0.5_rp + sign(0.5_rp,vel) ) &
1174  + ( - 3.0_rp * val(ks,i,j) &
1175  + 27.0_rp * val(ks+1,i,j) &
1176  + 47.0_rp * val(ks+2,i,j) &
1177  - 13.0_rp * val(ks+3,i,j) &
1178  + 2.0_rp * val(ks+4,i,j) ) / 60.0_rp &
1179  * ( 0.5_rp - sign(0.5_rp,vel) ) )
1180 
1181  vel = ( f2h(ke-2,1,i_uyz) &
1182  * mom(ke-1,i,j) &
1183  + f2h(ke-2,2,i_uyz) &
1184  * mom(ke-2,i,j) ) &
1185  / ( f2h(ke-2,1,i_uyz) &
1186  * 0.5_rp * ( dens(ke-1,i,j)+dens(ke-1,i+1,j) ) &
1187  + f2h(ke-2,2,i_uyz) &
1188  * 0.5_rp * ( dens(ke-2,i,j)+dens(ke-2,i+1,j) ) )
1189  vel = vel * j13g(ke-2,i,j)
1190  flux(ke-2,i,j) = vel / mapf(i,j,+2) &
1191  * ( ( - 3.0_rp * val(ke,i,j) &
1192  + 27.0_rp * val(ke-1,i,j) &
1193  + 47.0_rp * val(ke-2,i,j) &
1194  - 13.0_rp * val(ke-3,i,j) &
1195  + 2.0_rp * val(ke-4,i,j) ) / 60.0_rp &
1196  * ( 0.5_rp + sign(0.5_rp,vel) ) &
1197  + ( 2.0_rp * val(ke-2,i,j) + 5.0_rp * val(ke-1,i,j) - val(ke,i,j) ) / 6.0_rp &
1198  * ( 0.5_rp - sign(0.5_rp,vel) ) )
1199 
1200  flux(ke ,i,j) = 0.0_rp
1201  enddo
1202  enddo
1203 
1204  return
1205  end subroutine atmos_dyn_fvm_fluxj13_uyz_ud5
1206 
1207  !-----------------------------------------------------------------------------
1209  subroutine atmos_dyn_fvm_fluxj23_uyz_ud5( &
1210  flux, &
1211  mom, val, DENS, &
1212  GSQRT, J23G, MAPF, &
1213  CDZ, &
1214  IIS, IIE, JJS, JJE )
1215  implicit none
1216 
1217  real(RP), intent(inout) :: flux (ka,ia,ja)
1218  real(RP), intent(in) :: mom (ka,ia,ja)
1219  real(RP), intent(in) :: val (ka,ia,ja)
1220  real(RP), intent(in) :: DENS (ka,ia,ja)
1221  real(RP), intent(in) :: GSQRT (ka,ia,ja)
1222  real(RP), intent(in) :: J23G (ka,ia,ja)
1223  real(RP), intent(in) :: MAPF ( ia,ja,2)
1224  real(RP), intent(in) :: CDZ (ka)
1225  integer, intent(in) :: IIS, IIE, JJS, JJE
1226 
1227  real(RP) :: vel
1228  integer :: k, i, j
1229  !---------------------------------------------------------------------------
1230 
1231  !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
1232  !$omp private(vel) &
1233  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,J23G,MAPF) &
1234  !$omp shared(GSQRT,CDZ)
1235  do j = jjs, jje
1236  do i = iis, iie
1237  do k = ks+2, ke-3
1238  vel = ( f2h(k,1,i_uyz) &
1239  * 0.25_rp * ( mom(k+1,i,j)+mom(k+1,i+1,j)+mom(k+1,i,j-1)+mom(k+1,i+1,j-1) ) &
1240  + f2h(k,2,i_uyz) &
1241  * 0.25_rp * ( mom(k,i,j)+mom(k,i+1,j)+mom(k,i,j-1)+mom(k,i+1,j-1) ) ) &
1242  / ( f2h(k,1,i_uyz) &
1243  * 0.5_rp * ( dens(k+1,i,j)+dens(k+1,i+1,j) ) &
1244  + f2h(k,2,i_uyz) &
1245  * 0.5_rp * ( dens(k,i,j)+dens(k,i+1,j) ) )
1246  vel = vel * j23g(k,i,j)
1247  flux(k,i,j) = vel / mapf(i,j,+1) &
1248  * ( ( f51 * ( val(k+3,i,j)+val(k-2,i,j) ) &
1249  + f52 * ( val(k+2,i,j)+val(k-1,i,j) ) &
1250  + f53 * ( val(k+1,i,j)+val(k,i,j) ) ) &
1251  - ( f51 * ( val(k+3,i,j)-val(k-2,i,j) ) &
1252  + f54 * ( val(k+2,i,j)-val(k-1,i,j) ) &
1253  + f55 * ( val(k+1,i,j)-val(k,i,j) ) ) * sign(1.0_rp,vel) )
1254  enddo
1255  enddo
1256  enddo
1257 
1258  !$omp parallel do default(none) private(i,j) OMP_SCHEDULE_ collapse(2) &
1259  !$omp private(vel) &
1260  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,J23G,MAPF) &
1261  !$omp shared(GSQRT,CDZ)
1262  do j = jjs, jje
1263  do i = iis, iie
1264  ! The boundary condition is qflx_hi + qflxJ13 + qfluxJ23 = 0 at KS-1.
1265  ! The flux at KS-1 can be non-zero.
1266  ! To reduce calculations, all the fluxes are set to zero.
1267  flux(ks-1,i,j) = 0.0_rp
1268 
1269  vel = ( f2h(ks,1,i_uyz) &
1270  * 0.25_rp * ( mom(ks+1,i,j)+mom(ks+1,i+1,j)+mom(ks+1,i,j-1)+mom(ks+1,i+1,j-1) ) &
1271  + f2h(ks,2,i_uyz) &
1272  * 0.25_rp * ( mom(ks,i,j)+mom(ks,i+1,j)+mom(ks,i,j-1)+mom(ks,i+1,j-1) ) ) &
1273  / ( f2h(ks,1,i_uyz) &
1274  * 0.5_rp * ( dens(ks+1,i,j)+dens(ks+1,i+1,j) ) &
1275  + f2h(ks,2,i_uyz) &
1276  * 0.5_rp * ( dens(ks,i,j)+dens(ks,i+1,j) ) )
1277  vel = vel * j23g(ks,i,j)
1278  flux(ks,i,j) = vel / mapf(i,j,+1) &
1279  * ( f2 * ( val(ks+1,i,j)+val(ks,i,j) ) &
1280  * ( 0.5_rp + sign(0.5_rp,vel) ) &
1281  + ( 2.0_rp * val(ks,i,j) + 5.0_rp * val(ks+1,i,j) - val(ks+2,i,j) ) / 6.0_rp &
1282  * ( 0.5_rp - sign(0.5_rp,vel) ) )
1283 
1284  vel = ( f2h(ke-1,1,i_uyz) &
1285  * 0.25_rp * ( mom(ke,i,j)+mom(ke,i+1,j)+mom(ke,i,j-1)+mom(ke,i+1,j-1) ) &
1286  + f2h(ke-1,2,i_uyz) &
1287  * 0.25_rp * ( mom(ke-1,i,j)+mom(ke-1,i+1,j)+mom(ke-1,i,j-1)+mom(ke-1,i+1,j-1) ) ) &
1288  / ( f2h(ke-1,1,i_uyz) &
1289  * 0.5_rp * ( dens(ke,i,j)+dens(ke,i+1,j) ) &
1290  + f2h(ke-1,2,i_uyz) &
1291  * 0.5_rp * ( dens(ke-1,i,j)+dens(ke-1,i+1,j) ) )
1292  vel = vel * j23g(ke-1,i,j)
1293  flux(ke-1,i,j) = vel / mapf(i,j,+1) &
1294  * ( ( 2.0_rp * val(ke,i,j) + 5.0_rp * val(ke-1,i,j) - val(ke-2,i,j) ) / 6.0_rp &
1295  * ( 0.5_rp + sign(0.5_rp,vel) ) &
1296  + f2 * ( val(ke,i,j)+val(ke-1,i,j) ) &
1297  * ( 0.5_rp - sign(0.5_rp,vel) ) )
1298 
1299  vel = ( f2h(ks+1,1,i_uyz) &
1300  * 0.25_rp * ( mom(ks+2,i,j)+mom(ks+2,i+1,j)+mom(ks+2,i,j-1)+mom(ks+2,i+1,j-1) ) &
1301  + f2h(ks+1,2,i_uyz) &
1302  * 0.25_rp * ( mom(ks+1,i,j)+mom(ks+1,i+1,j)+mom(ks+1,i,j-1)+mom(ks+1,i+1,j-1) ) ) &
1303  / ( f2h(ks+1,1,i_uyz) &
1304  * 0.5_rp * ( dens(ks+2,i,j)+dens(ks+2,i+1,j) ) &
1305  + f2h(ks+1,2,i_uyz) &
1306  * 0.5_rp * ( dens(ks+1,i,j)+dens(ks+1,i+1,j) ) )
1307  vel = vel * j23g(ks+1,i,j)
1308  flux(ks+1,i,j) = vel / mapf(i,j,+1) &
1309  * ( ( 2.0_rp * val(ks+2,i,j) + 5.0_rp * val(ks+1,i,j) - val(ks,i,j) ) / 6.0_rp &
1310  * ( 0.5_rp + sign(0.5_rp,vel) ) &
1311  + ( - 3.0_rp * val(ks,i,j) &
1312  + 27.0_rp * val(ks+1,i,j) &
1313  + 47.0_rp * val(ks+2,i,j) &
1314  - 13.0_rp * val(ks+3,i,j) &
1315  + 2.0_rp * val(ks+4,i,j) ) / 60.0_rp &
1316  * ( 0.5_rp - sign(0.5_rp,vel) ) )
1317 
1318  vel = ( f2h(ke-2,1,i_uyz) &
1319  * 0.25_rp * ( mom(ke-1,i,j)+mom(ke-1,i+1,j)+mom(ke-1,i,j-1)+mom(ke-1,i+1,j-1) ) &
1320  + f2h(ke-2,2,i_uyz) &
1321  * 0.25_rp * ( mom(ke-2,i,j)+mom(ke-2,i+1,j)+mom(ke-2,i,j-1)+mom(ke-2,i+1,j-1) ) ) &
1322  / ( f2h(ke-2,1,i_uyz) &
1323  * 0.5_rp * ( dens(ke-1,i,j)+dens(ke-1,i+1,j) ) &
1324  + f2h(ke-2,2,i_uyz) &
1325  * 0.5_rp * ( dens(ke-2,i,j)+dens(ke-2,i+1,j) ) )
1326  vel = vel * j23g(ke-2,i,j)
1327  flux(ke-2,i,j) = vel / mapf(i,j,+1) &
1328  * ( ( - 3.0_rp * val(ke,i,j) &
1329  + 27.0_rp * val(ke-1,i,j) &
1330  + 47.0_rp * val(ke-2,i,j) &
1331  - 13.0_rp * val(ke-3,i,j) &
1332  + 2.0_rp * val(ke-4,i,j) ) / 60.0_rp &
1333  * ( 0.5_rp + sign(0.5_rp,vel) ) &
1334  + ( 2.0_rp * val(ke-2,i,j) + 5.0_rp * val(ke-1,i,j) - val(ke,i,j) ) / 6.0_rp &
1335  * ( 0.5_rp - sign(0.5_rp,vel) ) )
1336 
1337  flux(ke ,i,j) = 0.0_rp
1338  enddo
1339  enddo
1340 
1341  return
1342  end subroutine atmos_dyn_fvm_fluxj23_uyz_ud5
1343 
1344  !-----------------------------------------------------------------------------
1346  subroutine atmos_dyn_fvm_fluxx_uyz_ud5( &
1347  flux, &
1348  mom, val, DENS, &
1349  GSQRT, MAPF, &
1350  num_diff, &
1351  CDZ, &
1352  IIS, IIE, JJS, JJE )
1353  implicit none
1354 
1355  real(RP), intent(inout) :: flux (ka,ia,ja)
1356  real(RP), intent(in) :: mom (ka,ia,ja)
1357  real(RP), intent(in) :: val (ka,ia,ja)
1358  real(RP), intent(in) :: DENS (ka,ia,ja)
1359  real(RP), intent(in) :: GSQRT (ka,ia,ja)
1360  real(RP), intent(in) :: MAPF ( ia,ja,2)
1361  real(RP), intent(in) :: num_diff(ka,ia,ja)
1362  real(RP), intent(in) :: CDZ (ka)
1363  integer, intent(in) :: IIS, IIE, JJS, JJE
1364 
1365  real(RP) :: vel
1366  integer :: k, i, j
1367  !---------------------------------------------------------------------------
1368 
1369  ! note that x-index is added by -1
1370 
1371  !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
1372  !$omp private(vel) &
1373  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,GSQRT,MAPF,num_diff)
1374  do j = jjs, jje
1375  do i = iis, iie+1
1376  do k = ks, ke
1377 #ifdef DEBUG
1378  call check( __line__, mom(k,i ,j) )
1379  call check( __line__, mom(k,i-1,j) )
1380 
1381  call check( __line__, val(k,i-1,j) )
1382  call check( __line__, val(k,i,j) )
1383 
1384  call check( __line__, val(k,i-2,j) )
1385  call check( __line__, val(k,i+1,j) )
1386 
1387  call check( __line__, val(k,i-3,j) )
1388  call check( __line__, val(k,i+2,j) )
1389 
1390 #endif
1391  vel = ( 0.5_rp * ( mom(k,i,j)+mom(k,i-1,j) ) ) &
1392  / ( dens(k,i,j) )
1393  flux(k,i-1,j) = gsqrt(k,i,j) / mapf(i,j,+2) * vel &
1394  * ( ( f51 * ( val(k,i+2,j)+val(k,i-3,j) ) &
1395  + f52 * ( val(k,i+1,j)+val(k,i-2,j) ) &
1396  + f53 * ( val(k,i,j)+val(k,i-1,j) ) ) &
1397  - ( f51 * ( val(k,i+2,j)-val(k,i-3,j) ) &
1398  + f54 * ( val(k,i+1,j)-val(k,i-2,j) ) &
1399  + f55 * ( val(k,i,j)-val(k,i-1,j) ) ) * sign(1.0_rp,vel) ) &
1400  + gsqrt(k,i,j) * num_diff(k,i,j)
1401  enddo
1402  enddo
1403  enddo
1404 #ifdef DEBUG
1405  k = iundef; i = iundef; j = iundef
1406 #endif
1407 
1408  return
1409  end subroutine atmos_dyn_fvm_fluxx_uyz_ud5
1410 
1411  !-----------------------------------------------------------------------------
1413  subroutine atmos_dyn_fvm_fluxy_uyz_ud5( &
1414  flux, &
1415  mom, val, DENS, &
1416  GSQRT, MAPF, &
1417  num_diff, &
1418  CDZ, &
1419  IIS, IIE, JJS, JJE )
1420  implicit none
1421 
1422  real(RP), intent(inout) :: flux (ka,ia,ja)
1423  real(RP), intent(in) :: mom (ka,ia,ja)
1424  real(RP), intent(in) :: val (ka,ia,ja)
1425  real(RP), intent(in) :: DENS (ka,ia,ja)
1426  real(RP), intent(in) :: GSQRT (ka,ia,ja)
1427  real(RP), intent(in) :: MAPF ( ia,ja,2)
1428  real(RP), intent(in) :: num_diff(ka,ia,ja)
1429  real(RP), intent(in) :: CDZ (ka)
1430  integer, intent(in) :: IIS, IIE, JJS, JJE
1431 
1432  real(RP) :: vel
1433  integer :: k, i, j
1434  !---------------------------------------------------------------------------
1435 
1436  !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
1437  !$omp private(vel) &
1438  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,GSQRT,MAPF,num_diff)
1439  do j = jjs-1, jje
1440  do i = iis, iie
1441  do k = ks, ke
1442 #ifdef DEBUG
1443  call check( __line__, mom(k,i ,j) )
1444  call check( __line__, mom(k,i-1,j) )
1445 
1446  call check( __line__, val(k,i,j) )
1447  call check( __line__, val(k,i,j+1) )
1448 
1449  call check( __line__, val(k,i,j-1) )
1450  call check( __line__, val(k,i,j+2) )
1451 
1452  call check( __line__, val(k,i,j-2) )
1453  call check( __line__, val(k,i,j+3) )
1454 
1455 #endif
1456  vel = ( 0.5_rp * ( mom(k,i,j)+mom(k,i+1,j) ) ) &
1457  / ( 0.25_rp * ( dens(k,i,j)+dens(k,i+1,j)+dens(k,i,j+1)+dens(k,i+1,j+1) ) )
1458  flux(k,i,j) = gsqrt(k,i,j) / mapf(i,j,+1) * vel &
1459  * ( ( f51 * ( val(k,i,j+3)+val(k,i,j-2) ) &
1460  + f52 * ( val(k,i,j+2)+val(k,i,j-1) ) &
1461  + f53 * ( val(k,i,j+1)+val(k,i,j) ) ) &
1462  - ( f51 * ( val(k,i,j+3)-val(k,i,j-2) ) &
1463  + f54 * ( val(k,i,j+2)-val(k,i,j-1) ) &
1464  + f55 * ( val(k,i,j+1)-val(k,i,j) ) ) * sign(1.0_rp,vel) ) &
1465  + gsqrt(k,i,j) * num_diff(k,i,j)
1466  enddo
1467  enddo
1468  enddo
1469 #ifdef DEBUG
1470  k = iundef; i = iundef; j = iundef
1471 #endif
1472 
1473  return
1474  end subroutine atmos_dyn_fvm_fluxy_uyz_ud5
1475 
1476 
1477 
1478  !-----------------------------------------------------------------------------
1480  subroutine atmos_dyn_fvm_fluxz_xvz_ud5( &
1481  flux, &
1482  mom, val, DENS, &
1483  GSQRT, J33G, &
1484  num_diff, &
1485  CDZ, &
1486  IIS, IIE, JJS, JJE )
1487  implicit none
1488 
1489  real(RP), intent(inout) :: flux (ka,ia,ja)
1490  real(RP), intent(in) :: mom (ka,ia,ja)
1491  real(RP), intent(in) :: val (ka,ia,ja)
1492  real(RP), intent(in) :: DENS (ka,ia,ja)
1493  real(RP), intent(in) :: GSQRT (ka,ia,ja)
1494  real(RP), intent(in) :: J33G
1495  real(RP), intent(in) :: num_diff(ka,ia,ja)
1496  real(RP), intent(in) :: CDZ (ka)
1497  integer, intent(in) :: IIS, IIE, JJS, JJE
1498 
1499  real(RP) :: vel
1500  integer :: k, i, j
1501  !---------------------------------------------------------------------------
1502 
1503  !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
1504  !$omp private(vel) &
1505  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,J33G,GSQRT,num_diff) &
1506  !$omp shared(CDZ)
1507  do j = jjs, jje
1508  do i = iis, iie
1509  do k = ks+2, ke-3
1510 #ifdef DEBUG
1511  call check( __line__, mom(k,i,j) )
1512  call check( __line__, mom(k,i,j+1) )
1513 
1514  call check( __line__, val(k,i,j) )
1515  call check( __line__, val(k+1,i,j) )
1516 
1517  call check( __line__, val(k-1,i,j) )
1518  call check( __line__, val(k+2,i,j) )
1519 
1520  call check( __line__, val(k-2,i,j) )
1521  call check( __line__, val(k+3,i,j) )
1522 
1523 #endif
1524  vel = ( 0.5_rp * ( mom(k,i,j)+mom(k,i,j+1) ) ) &
1525  / ( f2h(k,1,i_xvz) &
1526  * 0.5_rp * ( dens(k+1,i,j)+dens(k+1,i,j+1) ) &
1527  + f2h(k,2,i_xvz) &
1528  * 0.5_rp * ( dens(k,i,j)+dens(k,i,j+1) ) )
1529  flux(k,i,j) = j33g * vel &
1530  * ( ( f51 * ( val(k+3,i,j)+val(k-2,i,j) ) &
1531  + f52 * ( val(k+2,i,j)+val(k-1,i,j) ) &
1532  + f53 * ( val(k+1,i,j)+val(k,i,j) ) ) &
1533  - ( f51 * ( val(k+3,i,j)-val(k-2,i,j) ) &
1534  + f54 * ( val(k+2,i,j)-val(k-1,i,j) ) &
1535  + f55 * ( val(k+1,i,j)-val(k,i,j) ) ) * sign(1.0_rp,vel) ) &
1536  + gsqrt(k,i,j) * num_diff(k,i,j)
1537  enddo
1538  enddo
1539  enddo
1540 #ifdef DEBUG
1541  k = iundef; i = iundef; j = iundef
1542 #endif
1543 
1544  !$omp parallel do default(none) private(i,j) OMP_SCHEDULE_ collapse(2) &
1545  !$omp private(vel) &
1546  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,flux,J33G,GSQRT,num_diff,DENS,CDZ)
1547  do j = jjs, jje
1548  do i = iis, iie
1549 #ifdef DEBUG
1550 
1551  call check( __line__, mom(ks,i ,j) )
1552  call check( __line__, mom(ks,i,j+1) )
1553  call check( __line__, val(ks+1,i,j) )
1554  call check( __line__, val(ks,i,j) )
1555 
1556  call check( __line__, mom(ks+1,i ,j) )
1557  call check( __line__, mom(ks+1,i,j+1) )
1558  call check( __line__, val(ks+3,i,j) )
1559  call check( __line__, val(ks+2,i,j) )
1560 
1561 #endif
1562  ! The boundary condition is qflx_hi + qflxJ13 + qfluxJ23 = 0 at KS-1.
1563  ! The flux at KS-1 can be non-zero.
1564  ! To reduce calculations, all the fluxes are set to zero.
1565  flux(ks-1,i,j) = 0.0_rp
1566 
1567  vel = ( 0.5_rp * ( mom(ks,i,j)+mom(ks,i,j+1) ) ) &
1568  / ( f2h(ks,1,i_xvz) &
1569  * 0.5_rp * ( dens(ks+1,i,j)+dens(ks+1,i,j+1) ) &
1570  + f2h(ks,2,i_xvz) &
1571  * 0.5_rp * ( dens(ks,i,j)+dens(ks,i,j+1) ) )
1572  flux(ks,i,j) = j33g * vel &
1573  * ( f2 * ( val(ks+1,i,j)+val(ks,i,j) ) &
1574  * ( 0.5_rp + sign(0.5_rp,vel) ) &
1575  + ( 2.0_rp * val(ks,i,j) + 5.0_rp * val(ks+1,i,j) - val(ks+2,i,j) ) / 6.0_rp &
1576  * ( 0.5_rp - sign(0.5_rp,vel) ) ) &
1577  + gsqrt(ks,i,j) * num_diff(ks,i,j)
1578  vel = ( 0.5_rp * ( mom(ke-1,i,j)+mom(ke-1,i,j+1) ) ) &
1579  / ( f2h(ke-1,1,i_xvz) &
1580  * 0.5_rp * ( dens(ke,i,j)+dens(ke,i,j+1) ) &
1581  + f2h(ke-1,2,i_xvz) &
1582  * 0.5_rp * ( dens(ke-1,i,j)+dens(ke-1,i,j+1) ) )
1583  flux(ke-1,i,j) = j33g * vel &
1584  * ( ( 2.0_rp * val(ke,i,j) + 5.0_rp * val(ke-1,i,j) - val(ke-2,i,j) ) / 6.0_rp &
1585  * ( 0.5_rp + sign(0.5_rp,vel) ) &
1586  + f2 * ( val(ke,i,j)+val(ke-1,i,j) ) &
1587  * ( 0.5_rp - sign(0.5_rp,vel) ) ) &
1588  + gsqrt(ke-1,i,j) * num_diff(ke-1,i,j)
1589 
1590  vel = ( 0.5_rp * ( mom(ks+1,i,j)+mom(ks+1,i,j+1) ) ) &
1591  / ( f2h(ks+1,1,i_xvz) &
1592  * 0.5_rp * ( dens(ks+2,i,j)+dens(ks+2,i,j+1) ) &
1593  + f2h(ks+1,2,i_xvz) &
1594  * 0.5_rp * ( dens(ks+1,i,j)+dens(ks+1,i,j+1) ) )
1595  flux(ks+1,i,j) = j33g * vel &
1596  * ( ( 2.0_rp * val(ks+2,i,j) + 5.0_rp * val(ks+1,i,j) - val(ks,i,j) ) / 6.0_rp &
1597  * ( 0.5_rp + sign(0.5_rp,vel) ) &
1598  + ( - 3.0_rp * val(ks,i,j) &
1599  + 27.0_rp * val(ks+1,i,j) &
1600  + 47.0_rp * val(ks+2,i,j) &
1601  - 13.0_rp * val(ks+3,i,j) &
1602  + 2.0_rp * val(ks+4,i,j) ) / 60.0_rp &
1603  * ( 0.5_rp - sign(0.5_rp,vel) ) ) &
1604  + gsqrt(ks+1,i,j) * num_diff(ks+1,i,j)
1605  vel = ( 0.5_rp * ( mom(ke-2,i,j)+mom(ke-2,i,j+1) ) ) &
1606  / ( f2h(ke-2,1,i_xvz) &
1607  * 0.5_rp * ( dens(ke-1,i,j)+dens(ke-1,i,j+1) ) &
1608  + f2h(ke-2,2,i_xvz) &
1609  * 0.5_rp * ( dens(ke-2,i,j)+dens(ke-2,i,j+1) ) )
1610  flux(ke-2,i,j) = j33g * vel &
1611  * ( ( - 3.0_rp * val(ke,i,j) &
1612  + 27.0_rp * val(ke-1,i,j) &
1613  + 47.0_rp * val(ke-2,i,j) &
1614  - 13.0_rp * val(ke-3,i,j) &
1615  + 2.0_rp * val(ke-4,i,j) ) / 60.0_rp &
1616  * ( 0.5_rp + sign(0.5_rp,vel) ) &
1617  + ( 2.0_rp * val(ke-2,i,j) + 5.0_rp * val(ke-1,i,j) - val(ke,i,j) ) / 6.0_rp &
1618  * ( 0.5_rp - sign(0.5_rp,vel) ) ) &
1619  + gsqrt(ke-2,i,j) * num_diff(ke-2,i,j)
1620 
1621  flux(ke,i,j) = 0.0_rp
1622  enddo
1623  enddo
1624 #ifdef DEBUG
1625  k = iundef; i = iundef; j = iundef
1626 #endif
1627 
1628  return
1629  end subroutine atmos_dyn_fvm_fluxz_xvz_ud5
1630 
1631  !-----------------------------------------------------------------------------
1633  subroutine atmos_dyn_fvm_fluxj13_xvz_ud5( &
1634  flux, &
1635  mom, val, DENS, &
1636  GSQRT, J13G, MAPF, &
1637  CDZ, &
1638  IIS, IIE, JJS, JJE )
1639  implicit none
1640 
1641  real(RP), intent(inout) :: flux (ka,ia,ja)
1642  real(RP), intent(in) :: mom (ka,ia,ja)
1643  real(RP), intent(in) :: val (ka,ia,ja)
1644  real(RP), intent(in) :: DENS (ka,ia,ja)
1645  real(RP), intent(in) :: GSQRT (ka,ia,ja)
1646  real(RP), intent(in) :: J13G (ka,ia,ja)
1647  real(RP), intent(in) :: MAPF ( ia,ja,2)
1648  real(RP), intent(in) :: CDZ (ka)
1649  integer, intent(in) :: IIS, IIE, JJS, JJE
1650 
1651  real(RP) :: vel
1652  integer :: k, i, j
1653  !---------------------------------------------------------------------------
1654 
1655  !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
1656  !$omp private(vel) &
1657  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,J13G,MAPF) &
1658  !$omp shared(GSQRT,CDZ)
1659  do j = jjs, jje
1660  do i = iis, iie
1661  do k = ks+2, ke-3
1662  vel = ( f2h(k,1,i_xvz) &
1663  * 0.25_rp * ( mom(k+1,i,j)+mom(k+1,i-1,j)+mom(k+1,i,j+1)+mom(k+1,i-1,j+1) ) &
1664  + f2h(k,2,i_xvz) &
1665  * 0.25_rp * ( mom(k,i,j)+mom(k,i-1,j)+mom(k,i,j+1)+mom(k,i-1,j+1) ) ) &
1666  / ( f2h(k,1,i_xvz) &
1667  * 0.5_rp * ( dens(k+1,i,j)+dens(k+1,i,j+1) ) &
1668  + f2h(k,2,i_xvz) &
1669  * 0.5_rp * ( dens(k,i,j)+dens(k,i,j+1) ) )
1670  vel = vel * j13g(k,i,j)
1671  flux(k,i,j) = vel / mapf(i,j,+2) &
1672  * ( ( f51 * ( val(k+3,i,j)+val(k-2,i,j) ) &
1673  + f52 * ( val(k+2,i,j)+val(k-1,i,j) ) &
1674  + f53 * ( val(k+1,i,j)+val(k,i,j) ) ) &
1675  - ( f51 * ( val(k+3,i,j)-val(k-2,i,j) ) &
1676  + f54 * ( val(k+2,i,j)-val(k-1,i,j) ) &
1677  + f55 * ( val(k+1,i,j)-val(k,i,j) ) ) * sign(1.0_rp,vel) )
1678  enddo
1679  enddo
1680  enddo
1681 
1682  !$omp parallel do default(none) private(i,j) OMP_SCHEDULE_ collapse(2) &
1683  !$omp private(vel) &
1684  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,J13G,MAPF) &
1685  !$omp shared(GSQRT,CDZ)
1686  do j = jjs, jje
1687  do i = iis, iie
1688  ! The boundary condition is qflx_hi + qflxJ13 + qfluxJ23 = 0 at KS-1.
1689  ! The flux at KS-1 can be non-zero.
1690  ! To reduce calculations, all the fluxes are set to zero.
1691  flux(ks-1,i,j) = 0.0_rp
1692 
1693  vel = ( f2h(ks,1,i_xvz) &
1694  * 0.25_rp * ( mom(ks+1,i,j)+mom(ks+1,i-1,j)+mom(ks+1,i,j+1)+mom(ks+1,i-1,j+1) ) &
1695  + f2h(ks,2,i_xvz) &
1696  * 0.25_rp * ( mom(ks,i,j)+mom(ks,i-1,j)+mom(ks,i,j+1)+mom(ks,i-1,j+1) ) ) &
1697  / ( f2h(ks,1,i_xvz) &
1698  * 0.5_rp * ( dens(ks+1,i,j)+dens(ks+1,i,j+1) ) &
1699  + f2h(ks,2,i_xvz) &
1700  * 0.5_rp * ( dens(ks,i,j)+dens(ks,i,j+1) ) )
1701  vel = vel * j13g(ks,i,j)
1702  flux(ks,i,j) = vel / mapf(i,j,+2) &
1703  * ( f2 * ( val(ks+1,i,j)+val(ks,i,j) ) &
1704  * ( 0.5_rp + sign(0.5_rp,vel) ) &
1705  + ( 2.0_rp * val(ks,i,j) + 5.0_rp * val(ks+1,i,j) - val(ks+2,i,j) ) / 6.0_rp &
1706  * ( 0.5_rp - sign(0.5_rp,vel) ) )
1707 
1708  vel = ( f2h(ke-1,1,i_xvz) &
1709  * 0.25_rp * ( mom(ke,i,j)+mom(ke,i-1,j)+mom(ke,i,j+1)+mom(ke,i-1,j+1) ) &
1710  + f2h(ke-1,2,i_xvz) &
1711  * 0.25_rp * ( mom(ke-1,i,j)+mom(ke-1,i-1,j)+mom(ke-1,i,j+1)+mom(ke-1,i-1,j+1) ) ) &
1712  / ( f2h(ke-1,1,i_xvz) &
1713  * 0.5_rp * ( dens(ke,i,j)+dens(ke,i,j+1) ) &
1714  + f2h(ke-1,2,i_xvz) &
1715  * 0.5_rp * ( dens(ke-1,i,j)+dens(ke-1,i,j+1) ) )
1716  vel = vel * j13g(ke-1,i,j)
1717  flux(ke-1,i,j) = vel / mapf(i,j,+2) &
1718  * ( ( 2.0_rp * val(ke,i,j) + 5.0_rp * val(ke-1,i,j) - val(ke-2,i,j) ) / 6.0_rp &
1719  * ( 0.5_rp + sign(0.5_rp,vel) ) &
1720  + f2 * ( val(ke,i,j)+val(ke-1,i,j) ) &
1721  * ( 0.5_rp - sign(0.5_rp,vel) ) )
1722 
1723  vel = ( f2h(ks+1,1,i_xvz) &
1724  * 0.25_rp * ( mom(ks+2,i,j)+mom(ks+2,i-1,j)+mom(ks+2,i,j+1)+mom(ks+2,i-1,j+1) ) &
1725  + f2h(ks+1,2,i_xvz) &
1726  * 0.25_rp * ( mom(ks+1,i,j)+mom(ks+1,i-1,j)+mom(ks+1,i,j+1)+mom(ks+1,i-1,j+1) ) ) &
1727  / ( f2h(ks+1,1,i_xvz) &
1728  * 0.5_rp * ( dens(ks+2,i,j)+dens(ks+2,i,j+1) ) &
1729  + f2h(ks+1,2,i_xvz) &
1730  * 0.5_rp * ( dens(ks+1,i,j)+dens(ks+1,i,j+1) ) )
1731  vel = vel * j13g(ks+1,i,j)
1732  flux(ks+1,i,j) = vel / mapf(i,j,+2) &
1733  * ( ( 2.0_rp * val(ks+2,i,j) + 5.0_rp * val(ks+1,i,j) - val(ks,i,j) ) / 6.0_rp &
1734  * ( 0.5_rp + sign(0.5_rp,vel) ) &
1735  + ( - 3.0_rp * val(ks,i,j) &
1736  + 27.0_rp * val(ks+1,i,j) &
1737  + 47.0_rp * val(ks+2,i,j) &
1738  - 13.0_rp * val(ks+3,i,j) &
1739  + 2.0_rp * val(ks+4,i,j) ) / 60.0_rp &
1740  * ( 0.5_rp - sign(0.5_rp,vel) ) )
1741 
1742  vel = ( f2h(ke-2,1,i_xvz) &
1743  * 0.25_rp * ( mom(ke-1,i,j)+mom(ke-1,i-1,j)+mom(ke-1,i,j+1)+mom(ke-1,i-1,j+1) ) &
1744  + f2h(ke-2,2,i_xvz) &
1745  * 0.25_rp * ( mom(ke-2,i,j)+mom(ke-2,i-1,j)+mom(ke-2,i,j+1)+mom(ke-2,i-1,j+1) ) ) &
1746  / ( f2h(ke-2,1,i_xvz) &
1747  * 0.5_rp * ( dens(ke-1,i,j)+dens(ke-1,i,j+1) ) &
1748  + f2h(ke-2,2,i_xvz) &
1749  * 0.5_rp * ( dens(ke-2,i,j)+dens(ke-2,i,j+1) ) )
1750  vel = vel * j13g(ke-2,i,j)
1751  flux(ke-2,i,j) = vel / mapf(i,j,+2) &
1752  * ( ( - 3.0_rp * val(ke,i,j) &
1753  + 27.0_rp * val(ke-1,i,j) &
1754  + 47.0_rp * val(ke-2,i,j) &
1755  - 13.0_rp * val(ke-3,i,j) &
1756  + 2.0_rp * val(ke-4,i,j) ) / 60.0_rp &
1757  * ( 0.5_rp + sign(0.5_rp,vel) ) &
1758  + ( 2.0_rp * val(ke-2,i,j) + 5.0_rp * val(ke-1,i,j) - val(ke,i,j) ) / 6.0_rp &
1759  * ( 0.5_rp - sign(0.5_rp,vel) ) )
1760 
1761  flux(ke ,i,j) = 0.0_rp
1762  enddo
1763  enddo
1764 
1765  return
1766  end subroutine atmos_dyn_fvm_fluxj13_xvz_ud5
1767 
1768  !-----------------------------------------------------------------------------
1770  subroutine atmos_dyn_fvm_fluxj23_xvz_ud5( &
1771  flux, &
1772  mom, val, DENS, &
1773  GSQRT, J23G, MAPF, &
1774  CDZ, &
1775  IIS, IIE, JJS, JJE )
1776  implicit none
1777 
1778  real(RP), intent(inout) :: flux (ka,ia,ja)
1779  real(RP), intent(in) :: mom (ka,ia,ja)
1780  real(RP), intent(in) :: val (ka,ia,ja)
1781  real(RP), intent(in) :: DENS (ka,ia,ja)
1782  real(RP), intent(in) :: GSQRT (ka,ia,ja)
1783  real(RP), intent(in) :: J23G (ka,ia,ja)
1784  real(RP), intent(in) :: MAPF ( ia,ja,2)
1785  real(RP), intent(in) :: CDZ (ka)
1786  integer, intent(in) :: IIS, IIE, JJS, JJE
1787 
1788  real(RP) :: vel
1789  integer :: k, i, j
1790  !---------------------------------------------------------------------------
1791 
1792  !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
1793  !$omp private(vel) &
1794  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,J23G,MAPF) &
1795  !$omp shared(GSQRT,CDZ)
1796  do j = jjs, jje
1797  do i = iis, iie
1798  do k = ks+2, ke-3
1799  vel = ( f2h(k,1,i_xvz) &
1800  * mom(k+1,i,j) &
1801  + f2h(k,2,i_xvz) &
1802  * mom(k,i,j) ) &
1803  / ( f2h(k,1,i_xvz) &
1804  * 0.5_rp * ( dens(k+1,i,j)+dens(k+1,i,j+1) ) &
1805  + f2h(k,2,i_xvz) &
1806  * 0.5_rp * ( dens(k,i,j)+dens(k,i,j+1) ) )
1807  vel = vel * j23g(k,i,j)
1808  flux(k,i,j) = vel / mapf(i,j,+1) &
1809  * ( ( f51 * ( val(k+3,i,j)+val(k-2,i,j) ) &
1810  + f52 * ( val(k+2,i,j)+val(k-1,i,j) ) &
1811  + f53 * ( val(k+1,i,j)+val(k,i,j) ) ) &
1812  - ( f51 * ( val(k+3,i,j)-val(k-2,i,j) ) &
1813  + f54 * ( val(k+2,i,j)-val(k-1,i,j) ) &
1814  + f55 * ( val(k+1,i,j)-val(k,i,j) ) ) * sign(1.0_rp,vel) )
1815  enddo
1816  enddo
1817  enddo
1818 
1819  !$omp parallel do default(none) private(i,j) OMP_SCHEDULE_ collapse(2) &
1820  !$omp private(vel) &
1821  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,J23G,MAPF) &
1822  !$omp shared(GSQRT,CDZ)
1823  do j = jjs, jje
1824  do i = iis, iie
1825  ! The boundary condition is qflx_hi + qflxJ13 + qfluxJ23 = 0 at KS-1.
1826  ! The flux at KS-1 can be non-zero.
1827  ! To reduce calculations, all the fluxes are set to zero.
1828  flux(ks-1,i,j) = 0.0_rp
1829 
1830  vel = ( f2h(ks,1,i_xvz) &
1831  * mom(ks+1,i,j) &
1832  + f2h(ks,2,i_xvz) &
1833  * mom(ks,i,j) ) &
1834  / ( f2h(ks,1,i_xvz) &
1835  * 0.5_rp * ( dens(ks+1,i,j)+dens(ks+1,i,j+1) ) &
1836  + f2h(ks,2,i_xvz) &
1837  * 0.5_rp * ( dens(ks,i,j)+dens(ks,i,j+1) ) )
1838  vel = vel * j23g(ks,i,j)
1839  flux(ks,i,j) = vel / mapf(i,j,+1) &
1840  * ( f2 * ( val(ks+1,i,j)+val(ks,i,j) ) &
1841  * ( 0.5_rp + sign(0.5_rp,vel) ) &
1842  + ( 2.0_rp * val(ks,i,j) + 5.0_rp * val(ks+1,i,j) - val(ks+2,i,j) ) / 6.0_rp &
1843  * ( 0.5_rp - sign(0.5_rp,vel) ) )
1844 
1845  vel = ( f2h(ke-1,1,i_xvz) &
1846  * mom(ke,i,j) &
1847  + f2h(ke-1,2,i_xvz) &
1848  * mom(ke-1,i,j) ) &
1849  / ( f2h(ke-1,1,i_xvz) &
1850  * 0.5_rp * ( dens(ke,i,j)+dens(ke,i,j+1) ) &
1851  + f2h(ke-1,2,i_xvz) &
1852  * 0.5_rp * ( dens(ke-1,i,j)+dens(ke-1,i,j+1) ) )
1853  vel = vel * j23g(ke-1,i,j)
1854  flux(ke-1,i,j) = vel / mapf(i,j,+1) &
1855  * ( ( 2.0_rp * val(ke,i,j) + 5.0_rp * val(ke-1,i,j) - val(ke-2,i,j) ) / 6.0_rp &
1856  * ( 0.5_rp + sign(0.5_rp,vel) ) &
1857  + f2 * ( val(ke,i,j)+val(ke-1,i,j) ) &
1858  * ( 0.5_rp - sign(0.5_rp,vel) ) )
1859 
1860  vel = ( f2h(ks+1,1,i_xvz) &
1861  * mom(ks+2,i,j) &
1862  + f2h(ks+1,2,i_xvz) &
1863  * mom(ks+1,i,j) ) &
1864  / ( f2h(ks+1,1,i_xvz) &
1865  * 0.5_rp * ( dens(ks+2,i,j)+dens(ks+2,i,j+1) ) &
1866  + f2h(ks+1,2,i_xvz) &
1867  * 0.5_rp * ( dens(ks+1,i,j)+dens(ks+1,i,j+1) ) )
1868  vel = vel * j23g(ks+1,i,j)
1869  flux(ks+1,i,j) = vel / mapf(i,j,+1) &
1870  * ( ( 2.0_rp * val(ks+2,i,j) + 5.0_rp * val(ks+1,i,j) - val(ks,i,j) ) / 6.0_rp &
1871  * ( 0.5_rp + sign(0.5_rp,vel) ) &
1872  + ( - 3.0_rp * val(ks,i,j) &
1873  + 27.0_rp * val(ks+1,i,j) &
1874  + 47.0_rp * val(ks+2,i,j) &
1875  - 13.0_rp * val(ks+3,i,j) &
1876  + 2.0_rp * val(ks+4,i,j) ) / 60.0_rp &
1877  * ( 0.5_rp - sign(0.5_rp,vel) ) )
1878 
1879  vel = ( f2h(ke-2,1,i_xvz) &
1880  * mom(ke-1,i,j) &
1881  + f2h(ke-2,2,i_xvz) &
1882  * mom(ke-2,i,j) ) &
1883  / ( f2h(ke-2,1,i_xvz) &
1884  * 0.5_rp * ( dens(ke-1,i,j)+dens(ke-1,i,j+1) ) &
1885  + f2h(ke-2,2,i_xvz) &
1886  * 0.5_rp * ( dens(ke-2,i,j)+dens(ke-2,i,j+1) ) )
1887  vel = vel * j23g(ke-2,i,j)
1888  flux(ke-2,i,j) = vel / mapf(i,j,+1) &
1889  * ( ( - 3.0_rp * val(ke,i,j) &
1890  + 27.0_rp * val(ke-1,i,j) &
1891  + 47.0_rp * val(ke-2,i,j) &
1892  - 13.0_rp * val(ke-3,i,j) &
1893  + 2.0_rp * val(ke-4,i,j) ) / 60.0_rp &
1894  * ( 0.5_rp + sign(0.5_rp,vel) ) &
1895  + ( 2.0_rp * val(ke-2,i,j) + 5.0_rp * val(ke-1,i,j) - val(ke,i,j) ) / 6.0_rp &
1896  * ( 0.5_rp - sign(0.5_rp,vel) ) )
1897 
1898  flux(ke ,i,j) = 0.0_rp
1899  enddo
1900  enddo
1901 
1902  return
1903  end subroutine atmos_dyn_fvm_fluxj23_xvz_ud5
1904 
1905  !-----------------------------------------------------------------------------
1907  subroutine atmos_dyn_fvm_fluxx_xvz_ud5( &
1908  flux, &
1909  mom, val, DENS, &
1910  GSQRT, MAPF, &
1911  num_diff, &
1912  CDZ, &
1913  IIS, IIE, JJS, JJE )
1914  implicit none
1915 
1916  real(RP), intent(inout) :: flux (ka,ia,ja)
1917  real(RP), intent(in) :: mom (ka,ia,ja)
1918  real(RP), intent(in) :: val (ka,ia,ja)
1919  real(RP), intent(in) :: DENS (ka,ia,ja)
1920  real(RP), intent(in) :: GSQRT (ka,ia,ja)
1921  real(RP), intent(in) :: MAPF ( ia,ja,2)
1922  real(RP), intent(in) :: num_diff(ka,ia,ja)
1923  real(RP), intent(in) :: CDZ (ka)
1924  integer, intent(in) :: IIS, IIE, JJS, JJE
1925 
1926  real(RP) :: vel
1927  integer :: k, i, j
1928  !---------------------------------------------------------------------------
1929 
1930  !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
1931  !$omp private(vel) &
1932  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,GSQRT,MAPF,num_diff)
1933  do j = jjs, jje
1934  do i = iis-1, iie
1935  do k = ks, ke
1936 #ifdef DEBUG
1937  call check( __line__, mom(k,i ,j) )
1938  call check( __line__, mom(k,i,j-1) )
1939 
1940  call check( __line__, val(k,i,j) )
1941  call check( __line__, val(k,i+1,j) )
1942 
1943  call check( __line__, val(k,i-1,j) )
1944  call check( __line__, val(k,i+2,j) )
1945 
1946  call check( __line__, val(k,i-2,j) )
1947  call check( __line__, val(k,i+3,j) )
1948 
1949 #endif
1950  vel = ( 0.5_rp * ( mom(k,i,j)+mom(k,i,j+1) ) ) &
1951  / ( 0.25_rp * ( dens(k,i,j)+dens(k,i+1,j)+dens(k,i,j+1)+dens(k,i+1,j+1) ) )
1952  flux(k,i,j) = gsqrt(k,i,j) / mapf(i,j,+2) * vel &
1953  * ( ( f51 * ( val(k,i+3,j)+val(k,i-2,j) ) &
1954  + f52 * ( val(k,i+2,j)+val(k,i-1,j) ) &
1955  + f53 * ( val(k,i+1,j)+val(k,i,j) ) ) &
1956  - ( f51 * ( val(k,i+3,j)-val(k,i-2,j) ) &
1957  + f54 * ( val(k,i+2,j)-val(k,i-1,j) ) &
1958  + f55 * ( val(k,i+1,j)-val(k,i,j) ) ) * sign(1.0_rp,vel) ) &
1959  + gsqrt(k,i,j) * num_diff(k,i,j)
1960  enddo
1961  enddo
1962  enddo
1963 #ifdef DEBUG
1964  k = iundef; i = iundef; j = iundef
1965 #endif
1966 
1967  return
1968  end subroutine atmos_dyn_fvm_fluxx_xvz_ud5
1969 
1970  !-----------------------------------------------------------------------------
1972  subroutine atmos_dyn_fvm_fluxy_xvz_ud5( &
1973  flux, &
1974  mom, val, DENS, &
1975  GSQRT, MAPF, &
1976  num_diff, &
1977  CDZ, &
1978  IIS, IIE, JJS, JJE )
1979  implicit none
1980 
1981  real(RP), intent(inout) :: flux (ka,ia,ja)
1982  real(RP), intent(in) :: mom (ka,ia,ja)
1983  real(RP), intent(in) :: val (ka,ia,ja)
1984  real(RP), intent(in) :: DENS (ka,ia,ja)
1985  real(RP), intent(in) :: GSQRT (ka,ia,ja)
1986  real(RP), intent(in) :: MAPF ( ia,ja,2)
1987  real(RP), intent(in) :: num_diff(ka,ia,ja)
1988  real(RP), intent(in) :: CDZ (ka)
1989  integer, intent(in) :: IIS, IIE, JJS, JJE
1990 
1991  real(RP) :: vel
1992  integer :: k, i, j
1993  !---------------------------------------------------------------------------
1994 
1995  ! note that y-index is added by -1
1996 
1997  !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
1998  !$omp private(vel) &
1999  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,GSQRT,MAPF,num_diff)
2000  do j = jjs, jje+1
2001  do i = iis, iie
2002  do k = ks, ke
2003 #ifdef DEBUG
2004  call check( __line__, mom(k,i ,j) )
2005  call check( __line__, mom(k,i,j-1) )
2006 
2007  call check( __line__, val(k,i,j-1) )
2008  call check( __line__, val(k,i,j) )
2009 
2010  call check( __line__, val(k,i,j-2) )
2011  call check( __line__, val(k,i,j+1) )
2012 
2013  call check( __line__, val(k,i,j-3) )
2014  call check( __line__, val(k,i,j+2) )
2015 
2016 #endif
2017  vel = ( 0.5_rp * ( mom(k,i,j)+mom(k,i,j-1) ) ) &
2018  / ( dens(k,i,j) )
2019  flux(k,i,j-1) = gsqrt(k,i,j) / mapf(i,j,+1) * vel &
2020  * ( ( f51 * ( val(k,i,j+2)+val(k,i,j-3) ) &
2021  + f52 * ( val(k,i,j+1)+val(k,i,j-2) ) &
2022  + f53 * ( val(k,i,j)+val(k,i,j-1) ) ) &
2023  - ( f51 * ( val(k,i,j+2)-val(k,i,j-3) ) &
2024  + f54 * ( val(k,i,j+1)-val(k,i,j-2) ) &
2025  + f55 * ( val(k,i,j)-val(k,i,j-1) ) ) * sign(1.0_rp,vel) ) &
2026  + gsqrt(k,i,j) * num_diff(k,i,j)
2027  enddo
2028  enddo
2029  enddo
2030 #ifdef DEBUG
2031  k = iundef; i = iundef; j = iundef
2032 #endif
2033 
2034  return
2035  end subroutine atmos_dyn_fvm_fluxy_xvz_ud5
2036 
2037 
2038 
2039 
2040 
2041 
2042 
2043 end module scale_atmos_dyn_fvm_flux_ud5
2044 
2045 !--
2046 ! vi:set readonly sw=4 ts=8
2047 !
2048 !Local Variables:
2049 !mode: f90
2050 !buffer-read-only: t
2051 !End:
2052 !
2053 !++
scale_atmos_dyn_fvm_flux_ud5::atmos_dyn_fvm_fluxy_xyw_ud5
subroutine, public atmos_dyn_fvm_fluxy_xyw_ud5(flux, mom, val, DENS, GSQRT, MAPF, num_diff, CDZ, IIS, IIE, JJS, JJE)
calculation Y-flux at XYW
Definition: scale_atmos_dyn_fvm_flux_ud5.F90:841
scale_debug
module DEBUG
Definition: scale_debug.F90:11
scale_atmos_grid_cartesc_index::ia
integer, public ia
of whole cells: x, local, with HALO
Definition: scale_atmos_grid_cartesC_index.F90:48
scale_atmos_dyn_fvm_flux_ud5::atmos_dyn_fvm_fluxj23_uyz_ud5
subroutine, public atmos_dyn_fvm_fluxj23_uyz_ud5(flux, mom, val, DENS, GSQRT, J23G, MAPF, CDZ, IIS, IIE, JJS, JJE)
calculation J23-flux at UYZ
Definition: scale_atmos_dyn_fvm_flux_ud5.F90:1215
scale_atmos_dyn_fvm_flux_ud5::atmos_dyn_fvm_fluxz_xvz_ud5
subroutine, public atmos_dyn_fvm_fluxz_xvz_ud5(flux, mom, val, DENS, GSQRT, J33G, num_diff, CDZ, IIS, IIE, JJS, JJE)
calculation z-flux at XV
Definition: scale_atmos_dyn_fvm_flux_ud5.F90:1487
scale_atmos_grid_cartesc_index::i_xvz
integer, public i_xvz
Definition: scale_atmos_grid_cartesC_index.F90:95
scale_atmos_dyn_fvm_flux_ud5::atmos_dyn_fvm_fluxj23_xvz_ud5
subroutine, public atmos_dyn_fvm_fluxj23_xvz_ud5(flux, mom, val, DENS, GSQRT, J23G, MAPF, CDZ, IIS, IIE, JJS, JJE)
calculation J23-flux at XVZ
Definition: scale_atmos_dyn_fvm_flux_ud5.F90:1776
scale_atmos_dyn_fvm_flux_ud5::atmos_dyn_fvm_fluxj13_xyw_ud5
subroutine, public atmos_dyn_fvm_fluxj13_xyw_ud5(flux, mom, val, DENS, GSQRT, J13G, MAPF, CDZ, IIS, IIE, JJS, JJE)
calculation J13-flux at XYW
Definition: scale_atmos_dyn_fvm_flux_ud5.F90:589
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_debug::check
subroutine, public check(current_line, v)
Undefined value checker.
Definition: scale_debug.F90:56
scale_const::const_undef
real(rp), public const_undef
Definition: scale_const.F90:41
scale_atmos_dyn_fvm_flux_ud5::atmos_dyn_fvm_fluxz_xyw_ud5
subroutine, public atmos_dyn_fvm_fluxz_xyw_ud5(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_ud5.F90:454
scale_atmos_dyn_fvm_flux_ud5
module scale_atmos_dyn_fvm_flux_ud5
Definition: scale_atmos_dyn_fvm_flux_ud5.F90:16
scale_tracer
module TRACER
Definition: scale_tracer.F90:16
scale_index
module Index
Definition: scale_index.F90:11
scale_atmos_dyn_fvm_flux_ud5::atmos_dyn_fvm_fluxj23_xyw_ud5
subroutine, public atmos_dyn_fvm_fluxj23_xyw_ud5(flux, mom, val, DENS, GSQRT, J23G, MAPF, CDZ, IIS, IIE, JJS, JJE)
calculation J23-flux at XYW
Definition: scale_atmos_dyn_fvm_flux_ud5.F90:672
scale_atmos_dyn_fvm_flux_ud5::atmos_dyn_fvm_fluxx_xyw_ud5
subroutine, public atmos_dyn_fvm_fluxx_xyw_ud5(flux, mom, val, DENS, GSQRT, MAPF, num_diff, CDZ, IIS, IIE, JJS, JJE)
calculation X-flux at XYW
Definition: scale_atmos_dyn_fvm_flux_ud5.F90:757
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_prc
module PROCESS
Definition: scale_prc.F90:11
scale_atmos_dyn_fvm_flux_ud5::atmos_dyn_fvm_fluxy_xvz_ud5
subroutine, public atmos_dyn_fvm_fluxy_xvz_ud5(flux, mom, val, DENS, GSQRT, MAPF, num_diff, CDZ, IIS, IIE, JJS, JJE)
calculation Y-flux at XV
Definition: scale_atmos_dyn_fvm_flux_ud5.F90:1979
scale_const::const_undef2
integer, parameter, public const_undef2
undefined value (INT2)
Definition: scale_const.F90:38
scale_atmos_dyn_fvm_flux_ud5::atmos_dyn_fvm_fluxj13_uyz_ud5
subroutine, public atmos_dyn_fvm_fluxj13_uyz_ud5(flux, mom, val, DENS, GSQRT, J13G, MAPF, CDZ, IIS, IIE, JJS, JJE)
calculation J13-flux at UYZ
Definition: scale_atmos_dyn_fvm_flux_ud5.F90:1078
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_ud5::atmos_dyn_fvm_fluxx_uyz_ud5
subroutine, public atmos_dyn_fvm_fluxx_uyz_ud5(flux, mom, val, DENS, GSQRT, MAPF, num_diff, CDZ, IIS, IIE, JJS, JJE)
calculation X-flux at UY
Definition: scale_atmos_dyn_fvm_flux_ud5.F90:1353
scale_const
module CONSTANT
Definition: scale_const.F90:11
scale_atmos_dyn_fvm_flux_ud5::atmos_dyn_fvm_fluxy_uyz_ud5
subroutine, public atmos_dyn_fvm_fluxy_uyz_ud5(flux, mom, val, DENS, GSQRT, MAPF, num_diff, CDZ, IIS, IIE, JJS, JJE)
calculation Y-flux at UY
Definition: scale_atmos_dyn_fvm_flux_ud5.F90:1420
scale_atmos_grid_cartesc_index::i_uyz
integer, public i_uyz
Definition: scale_atmos_grid_cartesC_index.F90:94
scale_atmos_dyn_fvm_flux_ud5::atmos_dyn_fvm_fluxj13_xvz_ud5
subroutine, public atmos_dyn_fvm_fluxj13_xvz_ud5(flux, mom, val, DENS, GSQRT, J13G, MAPF, CDZ, IIS, IIE, JJS, JJE)
calculation J13-flux at XVZ
Definition: scale_atmos_dyn_fvm_flux_ud5.F90:1639
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_precision
module PRECISION
Definition: scale_precision.F90:14
scale_atmos_dyn_fvm_flux_ud5::atmos_dyn_fvm_fluxz_uyz_ud5
subroutine, public atmos_dyn_fvm_fluxz_uyz_ud5(flux, mom, val, DENS, GSQRT, J33G, num_diff, CDZ, IIS, IIE, JJS, JJE)
calculation z-flux at UY
Definition: scale_atmos_dyn_fvm_flux_ud5.F90:926
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_fvm_flux_ud5::atmos_dyn_fvm_fluxz_xyz_ud5
subroutine, public atmos_dyn_fvm_fluxz_xyz_ud5(flux, mflx, val, GSQRT, num_diff, CDZ, IIS, IIE, JJS, JJE)
calculation z-flux at XYZ
Definition: scale_atmos_dyn_fvm_flux_ud5.F90:200
scale_atmos_dyn_fvm_flux_ud5::atmos_dyn_fvm_fluxx_xvz_ud5
subroutine, public atmos_dyn_fvm_fluxx_xvz_ud5(flux, mom, val, DENS, GSQRT, MAPF, num_diff, CDZ, IIS, IIE, JJS, JJE)
calculation X-flux at XV
Definition: scale_atmos_dyn_fvm_flux_ud5.F90:1914
scale_io
module STDIO
Definition: scale_io.F90:10
scale_atmos_dyn_fvm_flux_ud5::atmos_dyn_fvm_fluxx_xyz_ud5
subroutine, public atmos_dyn_fvm_fluxx_xyz_ud5(flux, mflx, val, GSQRT, num_diff, CDZ, IIS, IIE, JJS, JJE)
calculation X-flux at XYZ
Definition: scale_atmos_dyn_fvm_flux_ud5.F90:331
scale_atmos_dyn_fvm_flux_ud5::atmos_dyn_fvm_flux_valuew_z_ud5
subroutine, public atmos_dyn_fvm_flux_valuew_z_ud5(valW, mflx, val, GSQRT, CDZ)
value at XYW
Definition: scale_atmos_dyn_fvm_flux_ud5.F90:108
scale_atmos_dyn_fvm_flux_ud5::atmos_dyn_fvm_fluxy_xyz_ud5
subroutine, public atmos_dyn_fvm_fluxy_xyz_ud5(flux, mflx, val, GSQRT, num_diff, CDZ, IIS, IIE, JJS, JJE)
calculation Y-flux at XYZ
Definition: scale_atmos_dyn_fvm_flux_ud5.F90:391
Generated by   doxygen 1.8.13