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 
101 contains
102  !-----------------------------------------------------------------------------
104 !OCL SERIAL
105  subroutine atmos_dyn_fvm_flux_valuew_z_ud5( &
106  valW, &
107  mflx, val, GSQRT, &
108  CDZ )
109  !$acc routine vector
110  implicit none
111 
112  real(rp), intent(out) :: valw (ka)
113  real(rp), intent(in) :: mflx (ka)
114  real(rp), intent(in) :: val (ka)
115  real(rp), intent(in) :: gsqrt(ka)
116  real(rp), intent(in) :: cdz (ka)
117 
118  integer :: k
119  !---------------------------------------------------------------------------
120 
121  do k = ks+2, ke-3
122 #ifdef DEBUG
123  call check( __line__, mflx(k) )
124 
125  call check( __line__, val(k) )
126  call check( __line__, val(k+1) )
127 
128  call check( __line__, val(k-1) )
129  call check( __line__, val(k+2) )
130 
131  call check( __line__, val(k-2) )
132  call check( __line__, val(k+3) )
133 
134 #endif
135  valw(k) = ( f51 * ( val(k+3)+val(k-2) ) &
136  + f52 * ( val(k+2)+val(k-1) ) &
137  + f53 * ( val(k+1)+val(k) ) ) &
138  - ( f51 * ( val(k+3)-val(k-2) ) &
139  + f54 * ( val(k+2)-val(k-1) ) &
140  + f55 * ( val(k+1)-val(k) ) ) * sign(1.0_rp,mflx(k))
141  enddo
142 #ifdef DEBUG
143  k = iundef
144 #endif
145 
146 #ifdef DEBUG
147 
148  call check( __line__, mflx(ks) )
149  call check( __line__, val(ks ) )
150  call check( __line__, val(ks+1) )
151  call check( __line__, mflx(ke-1) )
152  call check( __line__, val(ke ) )
153  call check( __line__, val(ke-1) )
154 
155  call check( __line__, mflx(ks+1) )
156  call check( __line__, val(ks+2 ) )
157  call check( __line__, val(ks+3) )
158  call check( __line__, mflx(ke-2) )
159  call check( __line__, val(ke-2 ) )
160  call check( __line__, val(ke-3) )
161 
162 #endif
163 
164  valw(ks) = f2 * ( val(ks+1)+val(ks) ) &
165  * ( 0.5_rp + sign(0.5_rp,mflx(ks)) ) &
166  + ( 2.0_rp * val(ks) + 5.0_rp * val(ks+1) - val(ks+2) ) / 6.0_rp &
167  * ( 0.5_rp - sign(0.5_rp,mflx(ks)) )
168  valw(ke-1) = ( 2.0_rp * val(ke) + 5.0_rp * val(ke-1) - val(ke-2) ) / 6.0_rp &
169  * ( 0.5_rp + sign(0.5_rp,mflx(ke-1)) ) &
170  + f2 * ( val(ke)+val(ke-1) ) &
171  * ( 0.5_rp - sign(0.5_rp,mflx(ke-1)) )
172 
173  valw(ks+1) = ( 2.0_rp * val(ks+2) + 5.0_rp * val(ks+1) - val(ks) ) / 6.0_rp &
174  * ( 0.5_rp + sign(0.5_rp,mflx(ks+1)) ) &
175  + ( - 3.0_rp * val(ks) &
176  + 27.0_rp * val(ks+1) &
177  + 47.0_rp * val(ks+2) &
178  - 13.0_rp * val(ks+3) &
179  + 2.0_rp * val(ks+4) ) / 60.0_rp &
180  * ( 0.5_rp - sign(0.5_rp,mflx(ks+1)) )
181  valw(ke-2) = ( - 3.0_rp * val(ke) &
182  + 27.0_rp * val(ke-1) &
183  + 47.0_rp * val(ke-2) &
184  - 13.0_rp * val(ke-3) &
185  + 2.0_rp * val(ke-4) ) / 60.0_rp &
186  * ( 0.5_rp + sign(0.5_rp,mflx(ke-2)) ) &
187  + ( 2.0_rp * val(ke-2) + 5.0_rp * val(ke-1) - val(ke) ) / 6.0_rp &
188  * ( 0.5_rp - sign(0.5_rp,mflx(ke-2)) )
189 
190 
191  return
192  end subroutine atmos_dyn_fvm_flux_valuew_z_ud5
193 
194  !-----------------------------------------------------------------------------
196  subroutine atmos_dyn_fvm_fluxz_xyz_ud5( &
197  flux, &
198  mflx, val, GSQRT, &
199  num_diff, &
200  CDZ, &
201  IIS, IIE, JJS, JJE )
202  use scale_const, only: &
203  eps => const_eps
204  implicit none
205 
206  real(rp), intent(inout) :: flux (ka,ia,ja)
207  real(rp), intent(in) :: mflx (ka,ia,ja)
208  real(rp), intent(in) :: val (ka,ia,ja)
209  real(rp), intent(in) :: gsqrt (ka,ia,ja)
210  real(rp), intent(in) :: num_diff(ka,ia,ja)
211  real(rp), intent(in) :: cdz (ka)
212  integer, intent(in) :: iis, iie, jjs, jje
213 
214  real(rp) :: vel
215  integer :: k, i, j
216  !---------------------------------------------------------------------------
217 
218  !$omp parallel default(none) private(i,j,k, vel) &
219  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mflx,val,flux,GSQRT,num_diff,EPS)
220 
221  !$acc data copy(flux) copyin(mflx, val, GSQRT, num_diff, CDZ)
222 
223  !$omp do OMP_SCHEDULE_ collapse(2)
224  !$acc kernels
225  do j = jjs, jje
226  do i = iis, iie
227  do k = ks+2, ke-3
228 #ifdef DEBUG
229  call check( __line__, mflx(k,i,j) )
230 
231  call check( __line__, val(k,i,j) )
232  call check( __line__, val(k+1,i,j) )
233 
234  call check( __line__, val(k-1,i,j) )
235  call check( __line__, val(k+2,i,j) )
236 
237  call check( __line__, val(k-2,i,j) )
238  call check( __line__, val(k+3,i,j) )
239 
240 #endif
241  vel = mflx(k,i,j)
242  flux(k,i,j) = vel &
243  * ( ( f51 * ( val(k+3,i,j)+val(k-2,i,j) ) &
244  + f52 * ( val(k+2,i,j)+val(k-1,i,j) ) &
245  + f53 * ( val(k+1,i,j)+val(k,i,j) ) ) &
246  - ( f51 * ( val(k+3,i,j)-val(k-2,i,j) ) &
247  + f54 * ( val(k+2,i,j)-val(k-1,i,j) ) &
248  + f55 * ( val(k+1,i,j)-val(k,i,j) ) ) * sign(1.0_rp,vel) ) &
249  + gsqrt(k,i,j) * num_diff(k,i,j)
250  enddo
251  enddo
252  enddo
253  !$acc end kernels
254  !$omp end do nowait
255 #ifdef DEBUG
256  k = iundef; i = iundef; j = iundef
257 #endif
258 
259  !$omp do OMP_SCHEDULE_ collapse(2)
260  !$acc kernels
261  do j = jjs, jje
262  do i = iis, iie
263 #ifdef DEBUG
264 
265  call check( __line__, mflx(ks,i,j) )
266  call check( __line__, val(ks ,i,j) )
267  call check( __line__, val(ks+1,i,j) )
268  call check( __line__, mflx(ke-1,i,j) )
269  call check( __line__, val(ke ,i,j) )
270  call check( __line__, val(ke-1,i,j) )
271 
272  call check( __line__, mflx(ks+1,i,j) )
273  call check( __line__, val(ks+2 ,i,j) )
274  call check( __line__, val(ks+3,i,j) )
275  call check( __line__, mflx(ke-2,i,j) )
276  call check( __line__, val(ke-2 ,i,j) )
277  call check( __line__, val(ke-3,i,j) )
278 
279 #endif
280  flux(ks-1,i,j) = 0.0_rp
281 
282  vel = mflx(ks,i,j)
283  flux(ks,i,j) = vel &
284  * ( f2 * ( val(ks+1,i,j)+val(ks,i,j) ) &
285  * ( 0.5_rp + sign(0.5_rp,vel) ) &
286  + ( 2.0_rp * val(ks,i,j) + 5.0_rp * val(ks+1,i,j) - val(ks+2,i,j) ) / 6.0_rp &
287  * ( 0.5_rp - sign(0.5_rp,vel) ) ) &
288  + gsqrt(ks,i,j) * num_diff(ks,i,j)
289  vel = mflx(ke-1,i,j)
290  flux(ke-1,i,j) = vel &
291  * ( ( 2.0_rp * val(ke,i,j) + 5.0_rp * val(ke-1,i,j) - val(ke-2,i,j) ) / 6.0_rp &
292  * ( 0.5_rp + sign(0.5_rp,vel) ) &
293  + f2 * ( val(ke,i,j)+val(ke-1,i,j) ) &
294  * ( 0.5_rp - sign(0.5_rp,vel) ) ) &
295  + gsqrt(ke-1,i,j) * num_diff(ke-1,i,j)
296 
297  vel = mflx(ks+1,i,j)
298  flux(ks+1,i,j) = vel &
299  * ( ( 2.0_rp * val(ks+2,i,j) + 5.0_rp * val(ks+1,i,j) - val(ks,i,j) ) / 6.0_rp &
300  * ( 0.5_rp + sign(0.5_rp,vel) ) &
301  + ( - 3.0_rp * val(ks,i,j) &
302  + 27.0_rp * val(ks+1,i,j) &
303  + 47.0_rp * val(ks+2,i,j) &
304  - 13.0_rp * val(ks+3,i,j) &
305  + 2.0_rp * val(ks+4,i,j) ) / 60.0_rp &
306  * ( 0.5_rp - sign(0.5_rp,vel) ) ) &
307  + gsqrt(ks+1,i,j) * num_diff(ks+1,i,j)
308  vel = mflx(ke-2,i,j)
309  flux(ke-2,i,j) = vel &
310  * ( ( - 3.0_rp * val(ke,i,j) &
311  + 27.0_rp * val(ke-1,i,j) &
312  + 47.0_rp * val(ke-2,i,j) &
313  - 13.0_rp * val(ke-3,i,j) &
314  + 2.0_rp * val(ke-4,i,j) ) / 60.0_rp &
315  * ( 0.5_rp + sign(0.5_rp,vel) ) &
316  + ( 2.0_rp * val(ke-2,i,j) + 5.0_rp * val(ke-1,i,j) - val(ke,i,j) ) / 6.0_rp &
317  * ( 0.5_rp - sign(0.5_rp,vel) ) ) &
318  + gsqrt(ke-2,i,j) * num_diff(ke-2,i,j)
319 
320  flux(ke ,i,j) = 0.0_rp
321  enddo
322  enddo
323  !$acc end kernels
324  !$omp end do nowait
325 
326  !$acc end data
327 
328  !$omp end parallel
329 #ifdef DEBUG
330  k = iundef; i = iundef; j = iundef
331 #endif
332 
333  return
334  end subroutine atmos_dyn_fvm_fluxz_xyz_ud5
335 
336  !-----------------------------------------------------------------------------
338  subroutine atmos_dyn_fvm_fluxx_xyz_ud5( &
339  flux, &
340  mflx, val, GSQRT, &
341  num_diff, &
342  CDZ, &
343  IIS, IIE, JJS, JJE )
344  implicit none
345 
346  real(rp), intent(inout) :: flux (ka,ia,ja)
347  real(rp), intent(in) :: mflx (ka,ia,ja)
348  real(rp), intent(in) :: val (ka,ia,ja)
349  real(rp), intent(in) :: gsqrt (ka,ia,ja)
350  real(rp), intent(in) :: num_diff(ka,ia,ja)
351  real(rp), intent(in) :: cdz(ka)
352  integer, intent(in) :: iis, iie, jjs, jje
353 
354  real(rp) :: vel
355  integer :: k, i, j
356  !---------------------------------------------------------------------------
357 
358  !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
359  !$omp private(vel) &
360  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mflx,val,flux,GSQRT,num_diff)
361  !$acc kernels
362  do j = jjs, jje
363  do i = iis-1, iie
364  do k = ks, ke
365 #ifdef DEBUG
366  call check( __line__, mflx(k,i,j) )
367 
368  call check( __line__, val(k,i,j) )
369  call check( __line__, val(k,i+1,j) )
370 
371  call check( __line__, val(k,i-1,j) )
372  call check( __line__, val(k,i+2,j) )
373 
374  call check( __line__, val(k,i-2,j) )
375  call check( __line__, val(k,i+3,j) )
376 
377 #endif
378  vel = mflx(k,i,j)
379  flux(k,i,j) = vel &
380  * ( ( f51 * ( val(k,i+3,j)+val(k,i-2,j) ) &
381  + f52 * ( val(k,i+2,j)+val(k,i-1,j) ) &
382  + f53 * ( val(k,i+1,j)+val(k,i,j) ) ) &
383  - ( f51 * ( val(k,i+3,j)-val(k,i-2,j) ) &
384  + f54 * ( val(k,i+2,j)-val(k,i-1,j) ) &
385  + f55 * ( val(k,i+1,j)-val(k,i,j) ) ) * sign(1.0_rp,vel) ) &
386  + gsqrt(k,i,j) * num_diff(k,i,j)
387  enddo
388  enddo
389  enddo
390  !$acc end kernels
391 #ifdef DEBUG
392  k = iundef; i = iundef; j = iundef
393 #endif
394 
395  return
396  end subroutine atmos_dyn_fvm_fluxx_xyz_ud5
397 
398  !-----------------------------------------------------------------------------
400  subroutine atmos_dyn_fvm_fluxy_xyz_ud5( &
401  flux, &
402  mflx, val, GSQRT, &
403  num_diff, &
404  CDZ, &
405  IIS, IIE, JJS, JJE )
406  implicit none
407 
408  real(rp), intent(inout) :: flux (ka,ia,ja)
409  real(rp), intent(in) :: mflx (ka,ia,ja)
410  real(rp), intent(in) :: val (ka,ia,ja)
411  real(rp), intent(in) :: gsqrt (ka,ia,ja)
412  real(rp), intent(in) :: num_diff(ka,ia,ja)
413  real(rp), intent(in) :: cdz(ka)
414  integer, intent(in) :: iis, iie, jjs, jje
415 
416  real(rp) :: vel
417  integer :: k, i, j
418  !---------------------------------------------------------------------------
419 
420  !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
421  !$omp private(vel) &
422  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mflx,val,flux,GSQRT,num_diff)
423  !$acc kernels
424  do j = jjs-1, jje
425  do i = iis, iie
426  do k = ks, ke
427 #ifdef DEBUG
428  call check( __line__, mflx(k,i,j) )
429 
430  call check( __line__, val(k,i,j) )
431  call check( __line__, val(k,i,j+1) )
432 
433  call check( __line__, val(k,i,j-1) )
434  call check( __line__, val(k,i,j+2) )
435 
436  call check( __line__, val(k,i,j-2) )
437  call check( __line__, val(k,i,j+3) )
438 
439 #endif
440  vel = mflx(k,i,j)
441  flux(k,i,j) = vel &
442  * ( ( f51 * ( val(k,i,j+3)+val(k,i,j-2) ) &
443  + f52 * ( val(k,i,j+2)+val(k,i,j-1) ) &
444  + f53 * ( val(k,i,j+1)+val(k,i,j) ) ) &
445  - ( f51 * ( val(k,i,j+3)-val(k,i,j-2) ) &
446  + f54 * ( val(k,i,j+2)-val(k,i,j-1) ) &
447  + f55 * ( val(k,i,j+1)-val(k,i,j) ) ) * sign(1.0_rp,vel) ) &
448  + gsqrt(k,i,j) * num_diff(k,i,j)
449  enddo
450  enddo
451  enddo
452  !$acc end kernels
453 #ifdef DEBUG
454  k = iundef; i = iundef; j = iundef
455 #endif
456 
457  return
458  end subroutine atmos_dyn_fvm_fluxy_xyz_ud5
459 
460 
461  !-----------------------------------------------------------------------------
463  subroutine atmos_dyn_fvm_fluxz_xyw_ud5( &
464  flux, &
465  mom, val, DENS, &
466  GSQRT, J33G, &
467  num_diff, &
468  CDZ, FDZ, &
469  dtrk, &
470  IIS, IIE, JJS, JJE )
471  implicit none
472 
473  real(rp), intent(inout) :: flux (ka,ia,ja)
474  real(rp), intent(in) :: mom (ka,ia,ja)
475  real(rp), intent(in) :: val (ka,ia,ja)
476  real(rp), intent(in) :: dens (ka,ia,ja)
477  real(rp), intent(in) :: gsqrt (ka,ia,ja)
478  real(rp), intent(in) :: j33g
479  real(rp), intent(in) :: num_diff(ka,ia,ja)
480  real(rp), intent(in) :: cdz (ka)
481  real(rp), intent(in) :: fdz (ka-1)
482  real(rp), intent(in) :: dtrk
483  integer, intent(in) :: iis, iie, jjs, jje
484 
485  real(rp) :: vel
486  integer :: k, i, j
487  !---------------------------------------------------------------------------
488 
489  ! note than z-index is added by -1
490 
491  !$omp parallel default(none) private(i,j,k,vel) &
492  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,flux,J33G,GSQRT,num_diff,DENS,FDZ,dtrk)
493 
494  !$acc data copy(flux) copyin(mom, val, DENS, GSQRT, num_diff, CDZ, FDZ)
495 
496  !$omp do OMP_SCHEDULE_ collapse(2)
497  !$acc kernels
498  do j = jjs, jje
499  do i = iis, iie
500  do k = ks+3, ke-2
501 #ifdef DEBUG
502  call check( __line__, mom(k-1,i,j) )
503  call check( __line__, mom(k ,i,j) )
504 
505  call check( __line__, val(k-1,i,j) )
506  call check( __line__, val(k,i,j) )
507 
508  call check( __line__, val(k-2,i,j) )
509  call check( __line__, val(k+1,i,j) )
510 
511  call check( __line__, val(k-3,i,j) )
512  call check( __line__, val(k+2,i,j) )
513 
514 #endif
515  vel = ( 0.5_rp * ( mom(k-1,i,j) &
516  + mom(k,i,j) ) ) &
517  / dens(k,i,j)
518  flux(k-1,i,j) = j33g * vel &
519  * ( ( f51 * ( val(k+2,i,j)+val(k-3,i,j) ) &
520  + f52 * ( val(k+1,i,j)+val(k-2,i,j) ) &
521  + f53 * ( val(k,i,j)+val(k-1,i,j) ) ) &
522  - ( f51 * ( val(k+2,i,j)-val(k-3,i,j) ) &
523  + f54 * ( val(k+1,i,j)-val(k-2,i,j) ) &
524  + f55 * ( val(k,i,j)-val(k-1,i,j) ) ) * sign(1.0_rp,vel) ) &
525  + gsqrt(k,i,j) * num_diff(k,i,j)
526  enddo
527  enddo
528  enddo
529  !$acc end kernels
530  !$omp end do nowait
531 #ifdef DEBUG
532  k = iundef; i = iundef; j = iundef
533 #endif
534 
535  !$omp do OMP_SCHEDULE_ collapse(2)
536  !$acc kernels
537  do j = jjs, jje
538  do i = iis, iie
539 #ifdef DEBUG
540 
541  call check( __line__, val(ks,i,j) )
542  call check( __line__, val(ks+1,i,j) )
543  call check( __line__, val(ks+2,i,j) )
544  call check( __line__, val(ks+3,i,j) )
545  call check( __line__, val(ks+4,i,j) )
546 
547 
548  call check( __line__, val(ke-4,i,j) )
549  call check( __line__, val(ke-3,i,j) )
550  call check( __line__, val(ke-2,i,j) )
551  call check( __line__, val(ke-1,i,j) )
552 
553 #endif
554  ! The boundary condition is qflx_hi + qflxJ13 + qfluxJ23 = 0 at KS.
555  ! The flux at KS can be non-zero.
556  ! To reduce calculations, all the fluxes are set to zero.
557  flux(ks-1,i,j) = 0.0_rp ! k = KS
558 
559  vel = ( 0.5_rp * ( mom(ks,i,j) &
560  + mom(ks+1,i,j) ) ) &
561  / dens(ks+1,i,j)
562  flux(ks,i,j) = j33g * vel &
563  * ( f2 * ( val(ks+1,i,j)+val(ks,i,j) ) &
564  * ( 0.5_rp + sign(0.5_rp,vel) ) &
565  + ( 2.0_rp * val(ks,i,j) + 5.0_rp * val(ks+1,i,j) - val(ks+2,i,j) ) / 6.0_rp &
566  * ( 0.5_rp - sign(0.5_rp,vel) ) ) &
567  + gsqrt(ks+1,i,j) * num_diff(ks+1,i,j) ! k = KS+1
568 
569  vel = ( 0.5_rp * ( mom(ks+1,i,j) &
570  + mom(ks+2,i,j) ) ) &
571  / dens(ks+2,i,j)
572  flux(ks+1,i,j) = j33g * vel &
573  * ( ( 2.0_rp * val(ks+2,i,j) + 5.0_rp * val(ks+1,i,j) - val(ks,i,j) ) / 6.0_rp &
574  * ( 0.5_rp + sign(0.5_rp,vel) ) &
575  + ( - 3.0_rp * val(ks,i,j) &
576  + 27.0_rp * val(ks+1,i,j) &
577  + 47.0_rp * val(ks+2,i,j) &
578  - 13.0_rp * val(ks+3,i,j) &
579  + 2.0_rp * val(ks+4,i,j) ) / 60.0_rp &
580  * ( 0.5_rp - sign(0.5_rp,vel) ) ) &
581  + gsqrt(ks+2,i,j) * num_diff(ks+2,i,j) ! k = KS+2
582 
583 
584 
585  vel = ( 0.5_rp * ( mom(ke-2,i,j) &
586  + mom(ke-1,i,j) ) ) &
587  / dens(ke-1,i,j)
588  flux(ke-2,i,j) = j33g * vel &
589  * ( ( - 3.0_rp * val(ke,i,j) &
590  + 27.0_rp * val(ke-1,i,j) &
591  + 47.0_rp * val(ke-2,i,j) &
592  - 13.0_rp * val(ke-3,i,j) &
593  + 2.0_rp * val(ke-4,i,j) ) / 60.0_rp &
594  * ( 0.5_rp + sign(0.5_rp,vel) ) &
595  + ( 2.0_rp * val(ke-2,i,j) + 5.0_rp * val(ke-1,i,j) - val(ke,i,j) ) / 6.0_rp &
596  * ( 0.5_rp - sign(0.5_rp,vel) ) ) &
597  + gsqrt(ke-1,i,j) * num_diff(ke-1,i,j) ! k = KE-1
598 
599  flux(ke-1,i,j) = 0.0_rp ! k = KE
600  flux(ke ,i,j) = 0.0_rp ! k = KE+1
601  enddo
602  enddo
603  !$acc end kernels
604  !$omp end do nowait
605 
606  !$acc end data
607 
608  !$omp end parallel
609 
610  return
611  end subroutine atmos_dyn_fvm_fluxz_xyw_ud5
612 
613 
614  !-----------------------------------------------------------------------------
616  subroutine atmos_dyn_fvm_fluxj13_xyw_ud5( &
617  flux, &
618  mom, val, DENS, &
619  GSQRT, J13G, MAPF, &
620  CDZ, TwoD, &
621  IIS, IIE, JJS, JJE )
622  implicit none
623 
624  real(rp), intent(inout) :: flux (ka,ia,ja)
625  real(rp), intent(in) :: mom (ka,ia,ja)
626  real(rp), intent(in) :: val (ka,ia,ja)
627  real(rp), intent(in) :: dens (ka,ia,ja)
628  real(rp), intent(in) :: gsqrt (ka,ia,ja)
629  real(rp), intent(in) :: j13g (ka,ia,ja)
630  real(rp), intent(in) :: mapf ( ia,ja,2)
631  real(rp), intent(in) :: cdz (ka)
632  logical, intent(in) :: twod
633  integer, intent(in) :: iis, iie, jjs, jje
634 
635  real(rp) :: vel
636  integer :: k, i, j
637  !---------------------------------------------------------------------------
638 
639  !$omp parallel default(none) private(i,j,k,vel) &
640  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,J13G,MAPF)
641 
642  !$acc data copy(flux) copyin(mom, val, DENS, GSQRT, J13G, MAPF, CDZ)
643 
644  !$omp do OMP_SCHEDULE_ collapse(2)
645  !$acc kernels
646  do j = jjs, jje
647  do i = iis, iie
648  do k = ks+2, ke-2
649  vel = ( 0.5_rp * ( mom(k,i,j)+mom(k,i-1,j) ) ) &
650  / dens(k,i,j)
651  vel = vel * j13g(k,i,j)
652  flux(k-1,i,j) = vel / mapf(i,j,+2) &
653  * ( ( f51 * ( val(k+2,i,j)+val(k-3,i,j) ) &
654  + f52 * ( val(k+1,i,j)+val(k-2,i,j) ) &
655  + f53 * ( val(k,i,j)+val(k-1,i,j) ) ) &
656  - ( f51 * ( val(k+2,i,j)-val(k-3,i,j) ) &
657  + f54 * ( val(k+1,i,j)-val(k-2,i,j) ) &
658  + f55 * ( val(k,i,j)-val(k-1,i,j) ) ) * sign(1.0_rp,vel) )
659  enddo
660  enddo
661  enddo
662  !$acc end kernels
663  !$omp end do nowait
664 
665  !$omp do OMP_SCHEDULE_ collapse(2)
666  !$acc kernels
667  do j = jjs, jje
668  do i = iis, iie
669  ! The boundary condition is qflx_hi + qflxJ13 + qfluxJ23 = 0 at KS.
670  ! The flux at KS can be non-zero.
671  ! To reduce calculations, all the fluxes are set to zero.
672  flux(ks-1,i,j) = 0.0_rp ! k = KS
673 
674  ! physically incorrect but for numerical stability
675  vel = ( ( 0.5_rp * ( mom(ks+1,i,j)+mom(ks+1,i-1,j) ) ) / dens(ks+1,i,j) &
676  + ( 0.5_rp * ( mom(ks,i,j)+mom(ks,i-1,j) ) ) / dens(ks ,i,j) ) * 0.5_rp
677 ! vel = ( 0.5_RP * ( mom(KS+1,i,j)+mom(KS+1,i-1,j) ) ) &
678 ! / DENS(KS+1,i,j)
679  vel = vel * j13g(ks+1,i,j)
680  flux(ks,i,j) = vel / mapf(i,j,+2) &
681  * ( f2 * ( val(ks+1,i,j)+val(ks,i,j) ) &
682  * ( 0.5_rp + sign(0.5_rp,vel) ) &
683  + ( 2.0_rp * val(ks,i,j) + 5.0_rp * val(ks+1,i,j) - val(ks+2,i,j) ) / 6.0_rp &
684  * ( 0.5_rp - sign(0.5_rp,vel) ) ) ! k = KS+1
685 
686 
687  vel = ( 0.5_rp * ( mom(ke-1,i,j)+mom(ke-1,i-1,j) ) ) &
688  / dens(ke-1,i,j)
689  vel = vel * j13g(ke-1,i,j)
690  flux(ke-2,i,j) = vel / mapf(i,j,+2) &
691  * ( ( 2.0_rp * val(ke-1,i,j) + 5.0_rp * val(ke-2,i,j) - val(ke-3,i,j) ) / 6.0_rp &
692  * ( 0.5_rp + sign(0.5_rp,vel) ) &
693  + f2 * ( val(ke-1,i,j)+val(ke-2,i,j) ) &
694  * ( 0.5_rp - sign(0.5_rp,vel) ) ) ! k = KE-3
695 
696  flux(ke-1,i,j) = 0.0_rp
697  enddo
698  enddo
699  !$acc end kernels
700  !$omp end do nowait
701 
702  !$acc end data
703 
704  !$omp end parallel
705 
706  return
707  end subroutine atmos_dyn_fvm_fluxj13_xyw_ud5
708 
709  !-----------------------------------------------------------------------------
711  subroutine atmos_dyn_fvm_fluxj23_xyw_ud5( &
712  flux, &
713  mom, val, DENS, &
714  GSQRT, J23G, MAPF, &
715  CDZ, TwoD, &
716  IIS, IIE, JJS, JJE )
717  implicit none
718 
719  real(rp), intent(inout) :: flux (ka,ia,ja)
720  real(rp), intent(in) :: mom (ka,ia,ja)
721  real(rp), intent(in) :: val (ka,ia,ja)
722  real(rp), intent(in) :: dens (ka,ia,ja)
723  real(rp), intent(in) :: gsqrt (ka,ia,ja)
724  real(rp), intent(in) :: j23g (ka,ia,ja)
725  real(rp), intent(in) :: mapf ( ia,ja,2)
726  real(rp), intent(in) :: cdz (ka)
727  logical, intent(in) :: twod
728  integer, intent(in) :: iis, iie, jjs, jje
729 
730  real(rp) :: vel
731  integer :: k, i, j
732  !---------------------------------------------------------------------------
733 
734  !$omp parallel default(none) private(i,j,k,vel) &
735  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,J23G,MAPF)
736 
737  !$acc data copy(flux) copyin(mom, val, DENS, GSQRT, J23G, MAPF, CDZ)
738 
739  !$omp do OMP_SCHEDULE_ collapse(2)
740  !$acc kernels
741  do j = jjs, jje
742  do i = iis, iie
743  do k = ks+2, ke-2
744  vel = ( 0.5_rp * ( mom(k,i,j)+mom(k,i,j-1) ) ) &
745  / dens(k,i,j)
746  vel = vel * j23g(k,i,j)
747  flux(k-1,i,j) = vel / mapf(i,j,+1) &
748  * ( ( f51 * ( val(k+2,i,j)+val(k-3,i,j) ) &
749  + f52 * ( val(k+1,i,j)+val(k-2,i,j) ) &
750  + f53 * ( val(k,i,j)+val(k-1,i,j) ) ) &
751  - ( f51 * ( val(k+2,i,j)-val(k-3,i,j) ) &
752  + f54 * ( val(k+1,i,j)-val(k-2,i,j) ) &
753  + f55 * ( val(k,i,j)-val(k-1,i,j) ) ) * sign(1.0_rp,vel) )
754  enddo
755  enddo
756  enddo
757  !$acc end kernels
758  !$omp end do nowait
759 
760  !$omp do OMP_SCHEDULE_ collapse(2)
761  !$acc kernels
762  do j = jjs, jje
763  do i = iis, iie
764  ! The boundary condition is qflx_hi + qflxJ13 + qfluxJ23 = 0 at KS.
765  ! The flux at KS can be non-zero.
766  ! To reduce calculations, all the fluxes are set to zero.
767  flux(ks-1,i,j) = 0.0_rp ! k = KS
768 
769  ! physically incorrect but for numerical stability
770  vel = ( ( 0.5_rp * ( mom(ks+1,i,j)+mom(ks+1,i,j-1) ) ) / dens(ks+1,i,j) &
771  + ( 0.5_rp * ( mom(ks,i,j)+mom(ks,i,j-1) ) ) / dens(ks ,i,j) ) * 0.5_rp
772 ! vel = ( 0.5_RP * ( mom(KS+1,i,j)+mom(KS+1,i,j-1) ) ) &
773 ! / DENS(KS+1,i,j)
774  vel = vel * j23g(ks+1,i,j)
775  flux(ks,i,j) = vel / mapf(i,j,+1) &
776  * ( f2 * ( val(ks+1,i,j)+val(ks,i,j) ) &
777  * ( 0.5_rp + sign(0.5_rp,vel) ) &
778  + ( 2.0_rp * val(ks,i,j) + 5.0_rp * val(ks+1,i,j) - val(ks+2,i,j) ) / 6.0_rp &
779  * ( 0.5_rp - sign(0.5_rp,vel) ) ) ! k = KS+1
780 
781 
782  vel = ( 0.5_rp * ( mom(ke-1,i,j)+mom(ke-1,i,j-1) ) ) &
783  / dens(ke-1,i,j)
784  vel = vel * j23g(ke-1,i,j)
785  flux(ke-2,i,j) = vel / mapf(i,j,+1) &
786  * ( ( 2.0_rp * val(ke-1,i,j) + 5.0_rp * val(ke-2,i,j) - val(ke-3,i,j) ) / 6.0_rp &
787  * ( 0.5_rp + sign(0.5_rp,vel) ) &
788  + f2 * ( val(ke-1,i,j)+val(ke-2,i,j) ) &
789  * ( 0.5_rp - sign(0.5_rp,vel) ) ) ! k = KE-3
790 
791  flux(ke-1,i,j) = 0.0_rp
792  enddo
793  enddo
794  !$acc end kernels
795  !$omp end do nowait
796 
797  !$acc end data
798 
799  !$omp end parallel
800 
801  return
802  end subroutine atmos_dyn_fvm_fluxj23_xyw_ud5
803 
804 
805  !-----------------------------------------------------------------------------
807  subroutine atmos_dyn_fvm_fluxx_xyw_ud5( &
808  flux, &
809  mom, val, DENS, &
810  GSQRT, MAPF, &
811  num_diff, &
812  CDZ, TwoD, &
813  IIS, IIE, JJS, JJE )
814  implicit none
815 
816  real(rp), intent(inout) :: flux (ka,ia,ja)
817  real(rp), intent(in) :: mom (ka,ia,ja)
818  real(rp), intent(in) :: val (ka,ia,ja)
819  real(rp), intent(in) :: dens (ka,ia,ja)
820  real(rp), intent(in) :: gsqrt (ka,ia,ja)
821  real(rp), intent(in) :: mapf ( ia,ja,2)
822  real(rp), intent(in) :: num_diff(ka,ia,ja)
823  real(rp), intent(in) :: cdz (ka)
824  logical, intent(in) :: twod
825  integer, intent(in) :: iis, iie, jjs, jje
826 
827  real(rp) :: vel
828  integer :: k, i, j
829  !---------------------------------------------------------------------------
830 
831  !$omp parallel default(none) private(i,j,k,vel) &
832  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,GSQRT,MAPF,num_diff) &
833  !$omp shared(CDZ)
834 
835  !$acc data copy(flux) copyin(mom, val, DENS, GSQRT, MAPF, num_diff, CDZ)
836 
837  !$omp do OMP_SCHEDULE_ collapse(2)
838  !$acc kernels
839  do j = jjs, jje
840  do i = iis-1, iie
841  do k = ks, ke-1
842 #ifdef DEBUG
843  call check( __line__, mom(k ,i,j) )
844  call check( __line__, mom(k+1,i,j) )
845 
846  call check( __line__, val(k,i,j) )
847  call check( __line__, val(k,i+1,j) )
848 
849  call check( __line__, val(k,i-1,j) )
850  call check( __line__, val(k,i+2,j) )
851 
852  call check( __line__, val(k,i-2,j) )
853  call check( __line__, val(k,i+3,j) )
854 
855 #endif
856  vel = ( f2h(k,1,i_uyz) &
857  * mom(k+1,i,j) &
858  + f2h(k,2,i_uyz) &
859  * mom(k,i,j) ) &
860  / ( f2h(k,1,i_uyz) &
861  * 0.5_rp * ( dens(k+1,i,j)+dens(k+1,i+1,j) ) &
862  + f2h(k,2,i_uyz) &
863  * 0.5_rp * ( dens(k,i,j)+dens(k,i+1,j) ) )
864  flux(k,i,j) = gsqrt(k,i,j) / mapf(i,j,+2) * vel &
865  * ( ( f51 * ( val(k,i+3,j)+val(k,i-2,j) ) &
866  + f52 * ( val(k,i+2,j)+val(k,i-1,j) ) &
867  + f53 * ( val(k,i+1,j)+val(k,i,j) ) ) &
868  - ( f51 * ( val(k,i+3,j)-val(k,i-2,j) ) &
869  + f54 * ( val(k,i+2,j)-val(k,i-1,j) ) &
870  + f55 * ( val(k,i+1,j)-val(k,i,j) ) ) * sign(1.0_rp,vel) ) &
871  + gsqrt(k,i,j) * num_diff(k,i,j)
872  enddo
873  enddo
874  enddo
875  !$acc end kernels
876  !$omp end do nowait
877 #ifdef DEBUG
878  k = iundef; i = iundef; j = iundef
879 #endif
880 
881  !$omp do OMP_SCHEDULE_ collapse(2)
882  !$acc kernels
883  do j = jjs, jje
884  do i = iis-1, iie
885  flux(ke,i,j) = 0.0_rp
886  enddo
887  enddo
888  !$acc end kernels
889  !$omp end do nowait
890 
891  !$acc end data
892 
893  !$omp end parallel
894 #ifdef DEBUG
895  k = iundef; i = iundef; j = iundef
896 #endif
897 
898  return
899  end subroutine atmos_dyn_fvm_fluxx_xyw_ud5
900 
901  !-----------------------------------------------------------------------------
903  subroutine atmos_dyn_fvm_fluxy_xyw_ud5( &
904  flux, &
905  mom, val, DENS, &
906  GSQRT, MAPF, &
907  num_diff, &
908  CDZ, TwoD, &
909  IIS, IIE, JJS, JJE )
910  implicit none
911 
912  real(rp), intent(inout) :: flux (ka,ia,ja)
913  real(rp), intent(in) :: mom (ka,ia,ja)
914  real(rp), intent(in) :: val (ka,ia,ja)
915  real(rp), intent(in) :: dens (ka,ia,ja)
916  real(rp), intent(in) :: gsqrt (ka,ia,ja)
917  real(rp), intent(in) :: mapf ( ia,ja,2)
918  real(rp), intent(in) :: num_diff(ka,ia,ja)
919  real(rp), intent(in) :: cdz (ka)
920  logical, intent(in) :: twod
921  integer, intent(in) :: iis, iie, jjs, jje
922 
923  real(rp) :: vel
924  integer :: k, i, j
925  !---------------------------------------------------------------------------
926 
927  !$omp parallel default(none) private(i,j,k,vel) &
928  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,GSQRT,MAPF,num_diff) &
929  !$omp shared(CDZ)
930 
931  !$acc data copy(flux) copyin(mom, val, DENS, GSQRT, MAPF, num_diff, CDZ)
932 
933  !$omp do OMP_SCHEDULE_ collapse(2)
934  !$acc kernels
935  do j = jjs-1, jje
936  do i = iis, iie
937  do k = ks, ke-1
938 #ifdef DEBUG
939  call check( __line__, mom(k ,i,j) )
940  call check( __line__, mom(k+1,i,j) )
941 
942  call check( __line__, val(k,i,j) )
943  call check( __line__, val(k,i,j+1) )
944 
945  call check( __line__, val(k,i,j-1) )
946  call check( __line__, val(k,i,j+2) )
947 
948  call check( __line__, val(k,i,j-2) )
949  call check( __line__, val(k,i,j+3) )
950 
951 #endif
952  vel = ( f2h(k,1,i_xvz) &
953  * mom(k+1,i,j) &
954  + f2h(k,2,i_xvz) &
955  * mom(k,i,j) ) &
956  / ( f2h(k,1,i_xvz) &
957  * 0.5_rp * ( dens(k+1,i,j)+dens(k+1,i,j+1) ) &
958  + f2h(k,2,i_xvz) &
959  * 0.5_rp * ( dens(k,i,j)+dens(k,i,j+1) ) )
960  flux(k,i,j) = gsqrt(k,i,j) / mapf(i,j,+1) * vel &
961  * ( ( f51 * ( val(k,i,j+3)+val(k,i,j-2) ) &
962  + f52 * ( val(k,i,j+2)+val(k,i,j-1) ) &
963  + f53 * ( val(k,i,j+1)+val(k,i,j) ) ) &
964  - ( f51 * ( val(k,i,j+3)-val(k,i,j-2) ) &
965  + f54 * ( val(k,i,j+2)-val(k,i,j-1) ) &
966  + f55 * ( val(k,i,j+1)-val(k,i,j) ) ) * sign(1.0_rp,vel) ) &
967  + gsqrt(k,i,j) * num_diff(k,i,j)
968  enddo
969  enddo
970  enddo
971  !$acc end kernels
972  !$omp end do nowait
973 #ifdef DEBUG
974  k = iundef; i = iundef; j = iundef
975 #endif
976 
977  !$omp do OMP_SCHEDULE_ collapse(2)
978  !$acc kernels
979  do j = jjs-1, jje
980  do i = iis, iie
981  flux(ke,i,j) = 0.0_rp
982  enddo
983  enddo
984  !$acc end kernels
985  !$omp end do nowait
986 
987  !$acc end data
988 
989  !$omp end parallel
990 #ifdef DEBUG
991  k = iundef; i = iundef; j = iundef
992 #endif
993 
994  return
995  end subroutine atmos_dyn_fvm_fluxy_xyw_ud5
996 
997 
998  !-----------------------------------------------------------------------------
1000  subroutine atmos_dyn_fvm_fluxz_uyz_ud5( &
1001  flux, &
1002  mom, val, DENS, &
1003  GSQRT, J33G, &
1004  num_diff, &
1005  CDZ, TwoD, &
1006  IIS, IIE, JJS, JJE )
1007  implicit none
1008 
1009  real(rp), intent(inout) :: flux (ka,ia,ja)
1010  real(rp), intent(in) :: mom (ka,ia,ja)
1011  real(rp), intent(in) :: val (ka,ia,ja)
1012  real(rp), intent(in) :: dens (ka,ia,ja)
1013  real(rp), intent(in) :: gsqrt (ka,ia,ja)
1014  real(rp), intent(in) :: j33g
1015  real(rp), intent(in) :: num_diff(ka,ia,ja)
1016  real(rp), intent(in) :: cdz (ka)
1017  logical, intent(in) :: twod
1018  integer, intent(in) :: iis, iie, jjs, jje
1019 
1020  real(rp) :: vel
1021  integer :: k, i, j
1022  !---------------------------------------------------------------------------
1023 
1024  !$omp parallel default(none) private(i,j,k,vel) &
1025  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,J33G,GSQRT,num_diff) &
1026  !$omp shared(CDZ,TwoD)
1027 
1028  !$acc data copy(flux) copyin(mom, val, DENS, GSQRT, num_diff, CDZ)
1029 
1030 
1031  if ( twod ) then
1032 
1033  !$omp do OMP_SCHEDULE_
1034  !$acc kernels
1035  do j = jjs, jje
1036  do k = ks+2, ke-3
1037  i = iis
1038 #ifdef DEBUG
1039  call check( __line__, mom(k,i,j) )
1040 
1041  call check( __line__, val(k,i,j) )
1042  call check( __line__, val(k+1,i,j) )
1043 
1044  call check( __line__, val(k-1,i,j) )
1045  call check( __line__, val(k+2,i,j) )
1046 
1047  call check( __line__, val(k-2,i,j) )
1048  call check( __line__, val(k+3,i,j) )
1049 
1050 #endif
1051  vel = ( mom(k,i,j) ) &
1052  / ( f2h(k,1,i_xyz) &
1053  * dens(k+1,i,j) &
1054  + f2h(k,2,i_xyz) &
1055  * dens(k,i,j) )
1056  flux(k,i,j) = j33g * vel &
1057  * ( ( f51 * ( val(k+3,i,j)+val(k-2,i,j) ) &
1058  + f52 * ( val(k+2,i,j)+val(k-1,i,j) ) &
1059  + f53 * ( val(k+1,i,j)+val(k,i,j) ) ) &
1060  - ( f51 * ( val(k+3,i,j)-val(k-2,i,j) ) &
1061  + f54 * ( val(k+2,i,j)-val(k-1,i,j) ) &
1062  + f55 * ( val(k+1,i,j)-val(k,i,j) ) ) * sign(1.0_rp,vel) ) &
1063  + gsqrt(k,i,j) * num_diff(k,i,j)
1064  enddo
1065  enddo
1066  !$acc end kernels
1067  !$omp end do nowait
1068 #ifdef DEBUG
1069  k = iundef; i = iundef; j = iundef
1070 #endif
1071 
1072  !$omp do OMP_SCHEDULE_
1073  !$acc kernels
1074  do j = jjs, jje
1075  i = iis
1076 #ifdef DEBUG
1077 
1078  call check( __line__, mom(ks,i ,j) )
1079  call check( __line__, val(ks+1,i,j) )
1080  call check( __line__, val(ks,i,j) )
1081 
1082  call check( __line__, mom(ks+1,i ,j) )
1083  call check( __line__, val(ks+3,i,j) )
1084  call check( __line__, val(ks+2,i,j) )
1085 
1086 #endif
1087  ! The boundary condition is qflx_hi + qflxJ13 + qfluxJ23 = 0 at KS-1.
1088  ! The flux at KS-1 can be non-zero.
1089  ! To reduce calculations, all the fluxes are set to zero.
1090  flux(ks-1,i,j) = 0.0_rp
1091 
1092  vel = ( mom(ks,i,j) ) &
1093  / ( f2h(ks,1,i_xyz) &
1094  * dens(ks+1,i,j) &
1095  + f2h(ks,2,i_xyz) &
1096  * dens(ks,i,j) )
1097  flux(ks,i,j) = j33g * vel &
1098  * ( f2 * ( val(ks+1,i,j)+val(ks,i,j) ) &
1099  * ( 0.5_rp + sign(0.5_rp,vel) ) &
1100  + ( 2.0_rp * val(ks,i,j) + 5.0_rp * val(ks+1,i,j) - val(ks+2,i,j) ) / 6.0_rp &
1101  * ( 0.5_rp - sign(0.5_rp,vel) ) ) &
1102  + gsqrt(ks,i,j) * num_diff(ks,i,j)
1103  vel = ( mom(ke-1,i,j) ) &
1104  / ( f2h(ke-1,1,i_xyz) &
1105  * dens(ke,i,j) &
1106  + f2h(ke-1,2,i_xyz) &
1107  * dens(ke-1,i,j) )
1108  flux(ke-1,i,j) = j33g * vel &
1109  * ( ( 2.0_rp * val(ke,i,j) + 5.0_rp * val(ke-1,i,j) - val(ke-2,i,j) ) / 6.0_rp &
1110  * ( 0.5_rp + sign(0.5_rp,vel) ) &
1111  + f2 * ( val(ke,i,j)+val(ke-1,i,j) ) &
1112  * ( 0.5_rp - sign(0.5_rp,vel) ) ) &
1113  + gsqrt(ke-1,i,j) * num_diff(ke-1,i,j)
1114 
1115  vel = ( mom(ks+1,i,j) ) &
1116  / ( f2h(ks+1,1,i_xyz) &
1117  * dens(ks+2,i,j) &
1118  + f2h(ks+1,2,i_xyz) &
1119  * dens(ks+1,i,j) )
1120  flux(ks+1,i,j) = j33g * vel &
1121  * ( ( 2.0_rp * val(ks+2,i,j) + 5.0_rp * val(ks+1,i,j) - val(ks,i,j) ) / 6.0_rp &
1122  * ( 0.5_rp + sign(0.5_rp,vel) ) &
1123  + ( - 3.0_rp * val(ks,i,j) &
1124  + 27.0_rp * val(ks+1,i,j) &
1125  + 47.0_rp * val(ks+2,i,j) &
1126  - 13.0_rp * val(ks+3,i,j) &
1127  + 2.0_rp * val(ks+4,i,j) ) / 60.0_rp &
1128  * ( 0.5_rp - sign(0.5_rp,vel) ) ) &
1129  + gsqrt(ks+1,i,j) * num_diff(ks+1,i,j)
1130  vel = ( mom(ke-2,i,j) ) &
1131  / ( f2h(ke-2,1,i_xyz) &
1132  * dens(ke-1,i,j) &
1133  + f2h(ke-2,2,i_xyz) &
1134  * dens(ke-2,i,j) )
1135  flux(ke-2,i,j) = j33g * vel &
1136  * ( ( - 3.0_rp * val(ke,i,j) &
1137  + 27.0_rp * val(ke-1,i,j) &
1138  + 47.0_rp * val(ke-2,i,j) &
1139  - 13.0_rp * val(ke-3,i,j) &
1140  + 2.0_rp * val(ke-4,i,j) ) / 60.0_rp &
1141  * ( 0.5_rp + sign(0.5_rp,vel) ) &
1142  + ( 2.0_rp * val(ke-2,i,j) + 5.0_rp * val(ke-1,i,j) - val(ke,i,j) ) / 6.0_rp &
1143  * ( 0.5_rp - sign(0.5_rp,vel) ) ) &
1144  + gsqrt(ke-2,i,j) * num_diff(ke-2,i,j)
1145 
1146  flux(ke,i,j) = 0.0_rp
1147  enddo
1148  !$acc end kernels
1149  !$omp end do nowait
1150 
1151  else
1152 
1153 
1154  !$omp do OMP_SCHEDULE_ collapse(2)
1155  !$acc kernels
1156  do j = jjs, jje
1157  do i = iis, iie
1158  do k = ks+2, ke-3
1159 #ifdef DEBUG
1160  call check( __line__, mom(k,i,j) )
1161  call check( __line__, mom(k,i+1,j) )
1162 
1163  call check( __line__, val(k,i,j) )
1164  call check( __line__, val(k+1,i,j) )
1165 
1166  call check( __line__, val(k-1,i,j) )
1167  call check( __line__, val(k+2,i,j) )
1168 
1169  call check( __line__, val(k-2,i,j) )
1170  call check( __line__, val(k+3,i,j) )
1171 
1172 #endif
1173  vel = ( 0.5_rp * ( mom(k,i,j)+mom(k,i+1,j) ) ) &
1174  / ( f2h(k,1,i_uyz) &
1175  * 0.5_rp * ( dens(k+1,i,j)+dens(k+1,i+1,j) ) &
1176  + f2h(k,2,i_uyz) &
1177  * 0.5_rp * ( dens(k,i,j)+dens(k,i+1,j) ) )
1178  flux(k,i,j) = j33g * vel &
1179  * ( ( f51 * ( val(k+3,i,j)+val(k-2,i,j) ) &
1180  + f52 * ( val(k+2,i,j)+val(k-1,i,j) ) &
1181  + f53 * ( val(k+1,i,j)+val(k,i,j) ) ) &
1182  - ( f51 * ( val(k+3,i,j)-val(k-2,i,j) ) &
1183  + f54 * ( val(k+2,i,j)-val(k-1,i,j) ) &
1184  + f55 * ( val(k+1,i,j)-val(k,i,j) ) ) * sign(1.0_rp,vel) ) &
1185  + gsqrt(k,i,j) * num_diff(k,i,j)
1186  enddo
1187  enddo
1188  enddo
1189  !$acc end kernels
1190  !$omp end do nowait
1191 #ifdef DEBUG
1192  k = iundef; i = iundef; j = iundef
1193 #endif
1194 
1195  !$omp do OMP_SCHEDULE_ collapse(2)
1196  !$acc kernels
1197  do j = jjs, jje
1198  do i = iis, iie
1199 #ifdef DEBUG
1200 
1201  call check( __line__, mom(ks,i ,j) )
1202  call check( __line__, mom(ks,i+1,j) )
1203  call check( __line__, val(ks+1,i,j) )
1204  call check( __line__, val(ks,i,j) )
1205 
1206  call check( __line__, mom(ks+1,i ,j) )
1207  call check( __line__, mom(ks+1,i+1,j) )
1208  call check( __line__, val(ks+3,i,j) )
1209  call check( __line__, val(ks+2,i,j) )
1210 
1211 #endif
1212  ! The boundary condition is qflx_hi + qflxJ13 + qfluxJ23 = 0 at KS-1.
1213  ! The flux at KS-1 can be non-zero.
1214  ! To reduce calculations, all the fluxes are set to zero.
1215  flux(ks-1,i,j) = 0.0_rp
1216 
1217  vel = ( 0.5_rp * ( mom(ks,i,j)+mom(ks,i+1,j) ) ) &
1218  / ( f2h(ks,1,i_uyz) &
1219  * 0.5_rp * ( dens(ks+1,i,j)+dens(ks+1,i+1,j) ) &
1220  + f2h(ks,2,i_uyz) &
1221  * 0.5_rp * ( dens(ks,i,j)+dens(ks,i+1,j) ) )
1222  flux(ks,i,j) = j33g * vel &
1223  * ( f2 * ( val(ks+1,i,j)+val(ks,i,j) ) &
1224  * ( 0.5_rp + sign(0.5_rp,vel) ) &
1225  + ( 2.0_rp * val(ks,i,j) + 5.0_rp * val(ks+1,i,j) - val(ks+2,i,j) ) / 6.0_rp &
1226  * ( 0.5_rp - sign(0.5_rp,vel) ) ) &
1227  + gsqrt(ks,i,j) * num_diff(ks,i,j)
1228  vel = ( 0.5_rp * ( mom(ke-1,i,j)+mom(ke-1,i+1,j) ) ) &
1229  / ( f2h(ke-1,1,i_uyz) &
1230  * 0.5_rp * ( dens(ke,i,j)+dens(ke,i+1,j) ) &
1231  + f2h(ke-1,2,i_uyz) &
1232  * 0.5_rp * ( dens(ke-1,i,j)+dens(ke-1,i+1,j) ) )
1233  flux(ke-1,i,j) = j33g * vel &
1234  * ( ( 2.0_rp * val(ke,i,j) + 5.0_rp * val(ke-1,i,j) - val(ke-2,i,j) ) / 6.0_rp &
1235  * ( 0.5_rp + sign(0.5_rp,vel) ) &
1236  + f2 * ( val(ke,i,j)+val(ke-1,i,j) ) &
1237  * ( 0.5_rp - sign(0.5_rp,vel) ) ) &
1238  + gsqrt(ke-1,i,j) * num_diff(ke-1,i,j)
1239 
1240  vel = ( 0.5_rp * ( mom(ks+1,i,j)+mom(ks+1,i+1,j) ) ) &
1241  / ( f2h(ks+1,1,i_uyz) &
1242  * 0.5_rp * ( dens(ks+2,i,j)+dens(ks+2,i+1,j) ) &
1243  + f2h(ks+1,2,i_uyz) &
1244  * 0.5_rp * ( dens(ks+1,i,j)+dens(ks+1,i+1,j) ) )
1245  flux(ks+1,i,j) = j33g * vel &
1246  * ( ( 2.0_rp * val(ks+2,i,j) + 5.0_rp * val(ks+1,i,j) - val(ks,i,j) ) / 6.0_rp &
1247  * ( 0.5_rp + sign(0.5_rp,vel) ) &
1248  + ( - 3.0_rp * val(ks,i,j) &
1249  + 27.0_rp * val(ks+1,i,j) &
1250  + 47.0_rp * val(ks+2,i,j) &
1251  - 13.0_rp * val(ks+3,i,j) &
1252  + 2.0_rp * val(ks+4,i,j) ) / 60.0_rp &
1253  * ( 0.5_rp - sign(0.5_rp,vel) ) ) &
1254  + gsqrt(ks+1,i,j) * num_diff(ks+1,i,j)
1255  vel = ( 0.5_rp * ( mom(ke-2,i,j)+mom(ke-2,i+1,j) ) ) &
1256  / ( f2h(ke-2,1,i_uyz) &
1257  * 0.5_rp * ( dens(ke-1,i,j)+dens(ke-1,i+1,j) ) &
1258  + f2h(ke-2,2,i_uyz) &
1259  * 0.5_rp * ( dens(ke-2,i,j)+dens(ke-2,i+1,j) ) )
1260  flux(ke-2,i,j) = j33g * vel &
1261  * ( ( - 3.0_rp * val(ke,i,j) &
1262  + 27.0_rp * val(ke-1,i,j) &
1263  + 47.0_rp * val(ke-2,i,j) &
1264  - 13.0_rp * val(ke-3,i,j) &
1265  + 2.0_rp * val(ke-4,i,j) ) / 60.0_rp &
1266  * ( 0.5_rp + sign(0.5_rp,vel) ) &
1267  + ( 2.0_rp * val(ke-2,i,j) + 5.0_rp * val(ke-1,i,j) - val(ke,i,j) ) / 6.0_rp &
1268  * ( 0.5_rp - sign(0.5_rp,vel) ) ) &
1269  + gsqrt(ke-2,i,j) * num_diff(ke-2,i,j)
1270 
1271  flux(ke,i,j) = 0.0_rp
1272  enddo
1273  enddo
1274  !$acc end kernels
1275  !$omp end do nowait
1276 
1277  end if
1278 
1279 
1280  !$acc end data
1281 
1282  !$omp end parallel
1283 #ifdef DEBUG
1284  k = iundef; i = iundef; j = iundef
1285 #endif
1286 
1287  return
1288  end subroutine atmos_dyn_fvm_fluxz_uyz_ud5
1289 
1290  !-----------------------------------------------------------------------------
1292  subroutine atmos_dyn_fvm_fluxj13_uyz_ud5( &
1293  flux, &
1294  mom, val, DENS, &
1295  GSQRT, J13G, MAPF, &
1296  CDZ, TwoD, &
1297  IIS, IIE, JJS, JJE )
1298  implicit none
1299 
1300  real(rp), intent(inout) :: flux (ka,ia,ja)
1301  real(rp), intent(in) :: mom (ka,ia,ja)
1302  real(rp), intent(in) :: val (ka,ia,ja)
1303  real(rp), intent(in) :: dens (ka,ia,ja)
1304  real(rp), intent(in) :: gsqrt (ka,ia,ja)
1305  real(rp), intent(in) :: j13g (ka,ia,ja)
1306  real(rp), intent(in) :: mapf ( ia,ja,2)
1307  real(rp), intent(in) :: cdz (ka)
1308  logical, intent(in) :: twod
1309  integer, intent(in) :: iis, iie, jjs, jje
1310 
1311  real(rp) :: vel
1312  integer :: k, i, j
1313  !---------------------------------------------------------------------------
1314 
1315  !$omp parallel default(none) private(i,j,k,vel) &
1316  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,J13G,MAPF) &
1317  !$omp shared(GSQRT,CDZ,TwoD)
1318 
1319  !$acc data copy(flux) copyin(mom, val, DENS, GSQRT, J13G, MAPF, CDZ)
1320 
1321 
1322 
1323  !$omp do OMP_SCHEDULE_ collapse(2)
1324  !$acc kernels
1325  do j = jjs, jje
1326  do i = iis, iie
1327  do k = ks+2, ke-3
1328  vel = ( f2h(k,1,i_uyz) &
1329  * mom(k+1,i,j) &
1330  + f2h(k,2,i_uyz) &
1331  * mom(k,i,j) ) &
1332  / ( f2h(k,1,i_uyz) &
1333  * 0.5_rp * ( dens(k+1,i,j)+dens(k+1,i+1,j) ) &
1334  + f2h(k,2,i_uyz) &
1335  * 0.5_rp * ( dens(k,i,j)+dens(k,i+1,j) ) )
1336  vel = vel * j13g(k,i,j)
1337  flux(k,i,j) = vel / mapf(i,j,+2) &
1338  * ( ( f51 * ( val(k+3,i,j)+val(k-2,i,j) ) &
1339  + f52 * ( val(k+2,i,j)+val(k-1,i,j) ) &
1340  + f53 * ( val(k+1,i,j)+val(k,i,j) ) ) &
1341  - ( f51 * ( val(k+3,i,j)-val(k-2,i,j) ) &
1342  + f54 * ( val(k+2,i,j)-val(k-1,i,j) ) &
1343  + f55 * ( val(k+1,i,j)-val(k,i,j) ) ) * sign(1.0_rp,vel) )
1344  enddo
1345  enddo
1346  enddo
1347  !$acc end kernels
1348  !$omp end do nowait
1349 
1350  !$omp do OMP_SCHEDULE_ collapse(2)
1351  !$acc kernels
1352  do j = jjs, jje
1353  do i = iis, iie
1354  ! The boundary condition is qflx_hi + qflxJ13 + qfluxJ23 = 0 at KS-1.
1355  ! The flux at KS-1 can be non-zero.
1356  ! To reduce calculations, all the fluxes are set to zero.
1357  flux(ks-1,i,j) = 0.0_rp
1358 
1359  vel = ( f2h(ks,1,i_uyz) &
1360  * mom(ks+1,i,j) &
1361  + f2h(ks,2,i_uyz) &
1362  * mom(ks,i,j) ) &
1363  / ( f2h(ks,1,i_uyz) &
1364  * 0.5_rp * ( dens(ks+1,i,j)+dens(ks+1,i+1,j) ) &
1365  + f2h(ks,2,i_uyz) &
1366  * 0.5_rp * ( dens(ks,i,j)+dens(ks,i+1,j) ) )
1367  vel = vel * j13g(ks,i,j)
1368  flux(ks,i,j) = vel / mapf(i,j,+2) &
1369  * ( f2 * ( val(ks+1,i,j)+val(ks,i,j) ) &
1370  * ( 0.5_rp + sign(0.5_rp,vel) ) &
1371  + ( 2.0_rp * val(ks,i,j) + 5.0_rp * val(ks+1,i,j) - val(ks+2,i,j) ) / 6.0_rp &
1372  * ( 0.5_rp - sign(0.5_rp,vel) ) )
1373 
1374  vel = ( f2h(ke-1,1,i_uyz) &
1375  * mom(ke,i,j) &
1376  + f2h(ke-1,2,i_uyz) &
1377  * mom(ke-1,i,j) ) &
1378  / ( f2h(ke-1,1,i_uyz) &
1379  * 0.5_rp * ( dens(ke,i,j)+dens(ke,i+1,j) ) &
1380  + f2h(ke-1,2,i_uyz) &
1381  * 0.5_rp * ( dens(ke-1,i,j)+dens(ke-1,i+1,j) ) )
1382  vel = vel * j13g(ke-1,i,j)
1383  flux(ke-1,i,j) = vel / mapf(i,j,+2) &
1384  * ( ( 2.0_rp * val(ke,i,j) + 5.0_rp * val(ke-1,i,j) - val(ke-2,i,j) ) / 6.0_rp &
1385  * ( 0.5_rp + sign(0.5_rp,vel) ) &
1386  + f2 * ( val(ke,i,j)+val(ke-1,i,j) ) &
1387  * ( 0.5_rp - sign(0.5_rp,vel) ) )
1388 
1389  vel = ( f2h(ks+1,1,i_uyz) &
1390  * mom(ks+2,i,j) &
1391  + f2h(ks+1,2,i_uyz) &
1392  * mom(ks+1,i,j) ) &
1393  / ( f2h(ks+1,1,i_uyz) &
1394  * 0.5_rp * ( dens(ks+2,i,j)+dens(ks+2,i+1,j) ) &
1395  + f2h(ks+1,2,i_uyz) &
1396  * 0.5_rp * ( dens(ks+1,i,j)+dens(ks+1,i+1,j) ) )
1397  vel = vel * j13g(ks+1,i,j)
1398  flux(ks+1,i,j) = vel / mapf(i,j,+2) &
1399  * ( ( 2.0_rp * val(ks+2,i,j) + 5.0_rp * val(ks+1,i,j) - val(ks,i,j) ) / 6.0_rp &
1400  * ( 0.5_rp + sign(0.5_rp,vel) ) &
1401  + ( - 3.0_rp * val(ks,i,j) &
1402  + 27.0_rp * val(ks+1,i,j) &
1403  + 47.0_rp * val(ks+2,i,j) &
1404  - 13.0_rp * val(ks+3,i,j) &
1405  + 2.0_rp * val(ks+4,i,j) ) / 60.0_rp &
1406  * ( 0.5_rp - sign(0.5_rp,vel) ) )
1407 
1408  vel = ( f2h(ke-2,1,i_uyz) &
1409  * mom(ke-1,i,j) &
1410  + f2h(ke-2,2,i_uyz) &
1411  * mom(ke-2,i,j) ) &
1412  / ( f2h(ke-2,1,i_uyz) &
1413  * 0.5_rp * ( dens(ke-1,i,j)+dens(ke-1,i+1,j) ) &
1414  + f2h(ke-2,2,i_uyz) &
1415  * 0.5_rp * ( dens(ke-2,i,j)+dens(ke-2,i+1,j) ) )
1416  vel = vel * j13g(ke-2,i,j)
1417  flux(ke-2,i,j) = vel / mapf(i,j,+2) &
1418  * ( ( - 3.0_rp * val(ke,i,j) &
1419  + 27.0_rp * val(ke-1,i,j) &
1420  + 47.0_rp * val(ke-2,i,j) &
1421  - 13.0_rp * val(ke-3,i,j) &
1422  + 2.0_rp * val(ke-4,i,j) ) / 60.0_rp &
1423  * ( 0.5_rp + sign(0.5_rp,vel) ) &
1424  + ( 2.0_rp * val(ke-2,i,j) + 5.0_rp * val(ke-1,i,j) - val(ke,i,j) ) / 6.0_rp &
1425  * ( 0.5_rp - sign(0.5_rp,vel) ) )
1426 
1427  flux(ke ,i,j) = 0.0_rp
1428  enddo
1429  enddo
1430  !$acc end kernels
1431  !$omp end do nowait
1432 
1433 
1434 
1435  !$acc end data
1436 
1437  !$omp end parallel
1438  return
1439  end subroutine atmos_dyn_fvm_fluxj13_uyz_ud5
1440 
1441  !-----------------------------------------------------------------------------
1443  subroutine atmos_dyn_fvm_fluxj23_uyz_ud5( &
1444  flux, &
1445  mom, val, DENS, &
1446  GSQRT, J23G, MAPF, &
1447  CDZ, TwoD, &
1448  IIS, IIE, JJS, JJE )
1449  implicit none
1450 
1451  real(rp), intent(inout) :: flux (ka,ia,ja)
1452  real(rp), intent(in) :: mom (ka,ia,ja)
1453  real(rp), intent(in) :: val (ka,ia,ja)
1454  real(rp), intent(in) :: dens (ka,ia,ja)
1455  real(rp), intent(in) :: gsqrt (ka,ia,ja)
1456  real(rp), intent(in) :: j23g (ka,ia,ja)
1457  real(rp), intent(in) :: mapf ( ia,ja,2)
1458  real(rp), intent(in) :: cdz (ka)
1459  logical, intent(in) :: twod
1460  integer, intent(in) :: iis, iie, jjs, jje
1461 
1462  real(rp) :: vel
1463  integer :: k, i, j
1464  !---------------------------------------------------------------------------
1465 
1466  !$omp parallel default(none) private(i,j,k,vel) &
1467  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,J23G,MAPF) &
1468  !$omp shared(GSQRT,CDZ,TwoD)
1469 
1470  !$acc data copy(flux) copyin(mom, val, DENS, GSQRT, J23G, MAPF, CDZ)
1471 
1472 
1473  if ( twod ) then
1474 
1475  !$omp do OMP_SCHEDULE_
1476  !$acc kernels
1477  do j = jjs, jje
1478  do k = ks+2, ke-3
1479  i = iis
1480  vel = ( f2h(k,1,i_xyz) &
1481  * 0.5_rp * ( mom(k+1,i,j)+mom(k+1,i,j-1) ) &
1482  + f2h(k,2,i_xyz) &
1483  * 0.5_rp * ( mom(k,i,j)+mom(k,i,j-1) ) ) &
1484  / ( f2h(k,1,i_xyz) &
1485  * dens(k+1,i,j) &
1486  + f2h(k,2,i_xyz) &
1487  * dens(k,i,j) )
1488  vel = vel * j23g(k,i,j)
1489  flux(k,i,j) = vel * ( ( f51 * ( val(k+3,i,j)+val(k-2,i,j) ) &
1490  + f52 * ( val(k+2,i,j)+val(k-1,i,j) ) &
1491  + f53 * ( val(k+1,i,j)+val(k,i,j) ) ) &
1492  - ( f51 * ( val(k+3,i,j)-val(k-2,i,j) ) &
1493  + f54 * ( val(k+2,i,j)-val(k-1,i,j) ) &
1494  + f55 * ( val(k+1,i,j)-val(k,i,j) ) ) * sign(1.0_rp,vel) )
1495  enddo
1496  enddo
1497  !$acc end kernels
1498  !$omp end do nowait
1499 
1500  !$omp do OMP_SCHEDULE_
1501  !$acc kernels
1502  do j = jjs, jje
1503  i = iis
1504  ! The boundary condition is qflx_hi + qflxJ13 + qfluxJ23 = 0 at KS-1.
1505  ! The flux at KS-1 can be non-zero.
1506  ! To reduce calculations, all the fluxes are set to zero.
1507  flux(ks-1,i,j) = 0.0_rp
1508 
1509  vel = ( f2h(ks,1,i_xyz) &
1510  * 0.5_rp * ( mom(ks+1,i,j)+mom(ks+1,i,j-1) ) &
1511  + f2h(ks,2,i_xyz) &
1512  * 0.5_rp * ( mom(ks,i,j)+mom(ks,i,j-1) ) ) &
1513  / ( f2h(ks,1,i_xyz) &
1514  * dens(ks+1,i,j) &
1515  + f2h(ks,2,i_xyz) &
1516  * dens(ks,i,j) )
1517  vel = vel * j23g(ks,i,j)
1518  flux(ks,i,j) = vel / mapf(i,j,+1) &
1519  * ( f2 * ( val(ks+1,i,j)+val(ks,i,j) ) &
1520  * ( 0.5_rp + sign(0.5_rp,vel) ) &
1521  + ( 2.0_rp * val(ks,i,j) + 5.0_rp * val(ks+1,i,j) - val(ks+2,i,j) ) / 6.0_rp &
1522  * ( 0.5_rp - sign(0.5_rp,vel) ) )
1523 
1524  vel = ( f2h(ke-1,1,i_xyz) &
1525  * 0.5_rp * ( mom(ke,i,j)+mom(ke,i,j-1) ) &
1526  + f2h(ke-1,2,i_xyz) &
1527  * 0.5_rp * ( mom(ke-1,i,j)+mom(ke-1,i,j-1) ) ) &
1528  / ( f2h(ke-1,1,i_xyz) &
1529  * dens(ke,i,j) &
1530  + f2h(ke-1,2,i_xyz) &
1531  * dens(ke-1,i,j) )
1532  vel = vel * j23g(ke-1,i,j)
1533  flux(ke-1,i,j) = vel / mapf(i,j,+1) &
1534  * ( ( 2.0_rp * val(ke,i,j) + 5.0_rp * val(ke-1,i,j) - val(ke-2,i,j) ) / 6.0_rp &
1535  * ( 0.5_rp + sign(0.5_rp,vel) ) &
1536  + f2 * ( val(ke,i,j)+val(ke-1,i,j) ) &
1537  * ( 0.5_rp - sign(0.5_rp,vel) ) )
1538 
1539  vel = ( f2h(ks+1,1,i_xyz) &
1540  * 0.5_rp * ( mom(ks+2,i,j)+mom(ks+2,i,j-1) ) &
1541  + f2h(ks+1,2,i_xyz) &
1542  * 0.5_rp * ( mom(ks+1,i,j)+mom(ks+1,i,j-1) ) ) &
1543  / ( f2h(ks+1,1,i_xyz) &
1544  * dens(ks+2,i,j) &
1545  + f2h(ks+1,2,i_xyz) &
1546  * dens(ks+1,i,j) )
1547  vel = vel * j23g(ks+1,i,j)
1548  flux(ks+1,i,j) = vel / mapf(i,j,+1) &
1549  * ( ( 2.0_rp * val(ks+2,i,j) + 5.0_rp * val(ks+1,i,j) - val(ks,i,j) ) / 6.0_rp &
1550  * ( 0.5_rp + sign(0.5_rp,vel) ) &
1551  + ( - 3.0_rp * val(ks,i,j) &
1552  + 27.0_rp * val(ks+1,i,j) &
1553  + 47.0_rp * val(ks+2,i,j) &
1554  - 13.0_rp * val(ks+3,i,j) &
1555  + 2.0_rp * val(ks+4,i,j) ) / 60.0_rp &
1556  * ( 0.5_rp - sign(0.5_rp,vel) ) )
1557 
1558  vel = ( f2h(ke-2,1,i_xyz) &
1559  * 0.5_rp * ( mom(ke-1,i,j)+mom(ke-1,i,j-1) ) &
1560  + f2h(ke-2,2,i_xyz) &
1561  * 0.5_rp * ( mom(ke-2,i,j)+mom(ke-2,i,j-1) ) ) &
1562  / ( f2h(ke-2,1,i_xyz) &
1563  * dens(ke-1,i,j) &
1564  + f2h(ke-2,2,i_xyz) &
1565  * dens(ke-2,i,j) )
1566  vel = vel * j23g(ke-2,i,j)
1567  flux(ke-2,i,j) = vel / mapf(i,j,+1) &
1568  * ( ( - 3.0_rp * val(ke,i,j) &
1569  + 27.0_rp * val(ke-1,i,j) &
1570  + 47.0_rp * val(ke-2,i,j) &
1571  - 13.0_rp * val(ke-3,i,j) &
1572  + 2.0_rp * val(ke-4,i,j) ) / 60.0_rp &
1573  * ( 0.5_rp + sign(0.5_rp,vel) ) &
1574  + ( 2.0_rp * val(ke-2,i,j) + 5.0_rp * val(ke-1,i,j) - val(ke,i,j) ) / 6.0_rp &
1575  * ( 0.5_rp - sign(0.5_rp,vel) ) )
1576 
1577  flux(ke ,i,j) = 0.0_rp
1578  enddo
1579  !$acc end kernels
1580  !$omp end do nowait
1581 
1582  else
1583 
1584 
1585  !$omp do OMP_SCHEDULE_ collapse(2)
1586  !$acc kernels
1587  do j = jjs, jje
1588  do i = iis, iie
1589  do k = ks+2, ke-3
1590  vel = ( f2h(k,1,i_uyz) &
1591  * 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) ) &
1592  + f2h(k,2,i_uyz) &
1593  * 0.25_rp * ( mom(k,i,j)+mom(k,i+1,j)+mom(k,i,j-1)+mom(k,i+1,j-1) ) ) &
1594  / ( f2h(k,1,i_uyz) &
1595  * 0.5_rp * ( dens(k+1,i,j)+dens(k+1,i+1,j) ) &
1596  + f2h(k,2,i_uyz) &
1597  * 0.5_rp * ( dens(k,i,j)+dens(k,i+1,j) ) )
1598  vel = vel * j23g(k,i,j)
1599  flux(k,i,j) = vel / mapf(i,j,+1) &
1600  * ( ( f51 * ( val(k+3,i,j)+val(k-2,i,j) ) &
1601  + f52 * ( val(k+2,i,j)+val(k-1,i,j) ) &
1602  + f53 * ( val(k+1,i,j)+val(k,i,j) ) ) &
1603  - ( f51 * ( val(k+3,i,j)-val(k-2,i,j) ) &
1604  + f54 * ( val(k+2,i,j)-val(k-1,i,j) ) &
1605  + f55 * ( val(k+1,i,j)-val(k,i,j) ) ) * sign(1.0_rp,vel) )
1606  enddo
1607  enddo
1608  enddo
1609  !$acc end kernels
1610  !$omp end do nowait
1611 
1612  !$omp do OMP_SCHEDULE_ collapse(2)
1613  !$acc kernels
1614  do j = jjs, jje
1615  do i = iis, iie
1616  ! The boundary condition is qflx_hi + qflxJ13 + qfluxJ23 = 0 at KS-1.
1617  ! The flux at KS-1 can be non-zero.
1618  ! To reduce calculations, all the fluxes are set to zero.
1619  flux(ks-1,i,j) = 0.0_rp
1620 
1621  vel = ( f2h(ks,1,i_uyz) &
1622  * 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) ) &
1623  + f2h(ks,2,i_uyz) &
1624  * 0.25_rp * ( mom(ks,i,j)+mom(ks,i+1,j)+mom(ks,i,j-1)+mom(ks,i+1,j-1) ) ) &
1625  / ( f2h(ks,1,i_uyz) &
1626  * 0.5_rp * ( dens(ks+1,i,j)+dens(ks+1,i+1,j) ) &
1627  + f2h(ks,2,i_uyz) &
1628  * 0.5_rp * ( dens(ks,i,j)+dens(ks,i+1,j) ) )
1629  vel = vel * j23g(ks,i,j)
1630  flux(ks,i,j) = vel / mapf(i,j,+1) &
1631  * ( f2 * ( val(ks+1,i,j)+val(ks,i,j) ) &
1632  * ( 0.5_rp + sign(0.5_rp,vel) ) &
1633  + ( 2.0_rp * val(ks,i,j) + 5.0_rp * val(ks+1,i,j) - val(ks+2,i,j) ) / 6.0_rp &
1634  * ( 0.5_rp - sign(0.5_rp,vel) ) )
1635 
1636  vel = ( f2h(ke-1,1,i_uyz) &
1637  * 0.25_rp * ( mom(ke,i,j)+mom(ke,i+1,j)+mom(ke,i,j-1)+mom(ke,i+1,j-1) ) &
1638  + f2h(ke-1,2,i_uyz) &
1639  * 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) ) ) &
1640  / ( f2h(ke-1,1,i_uyz) &
1641  * 0.5_rp * ( dens(ke,i,j)+dens(ke,i+1,j) ) &
1642  + f2h(ke-1,2,i_uyz) &
1643  * 0.5_rp * ( dens(ke-1,i,j)+dens(ke-1,i+1,j) ) )
1644  vel = vel * j23g(ke-1,i,j)
1645  flux(ke-1,i,j) = vel / mapf(i,j,+1) &
1646  * ( ( 2.0_rp * val(ke,i,j) + 5.0_rp * val(ke-1,i,j) - val(ke-2,i,j) ) / 6.0_rp &
1647  * ( 0.5_rp + sign(0.5_rp,vel) ) &
1648  + f2 * ( val(ke,i,j)+val(ke-1,i,j) ) &
1649  * ( 0.5_rp - sign(0.5_rp,vel) ) )
1650 
1651  vel = ( f2h(ks+1,1,i_uyz) &
1652  * 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) ) &
1653  + f2h(ks+1,2,i_uyz) &
1654  * 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) ) ) &
1655  / ( f2h(ks+1,1,i_uyz) &
1656  * 0.5_rp * ( dens(ks+2,i,j)+dens(ks+2,i+1,j) ) &
1657  + f2h(ks+1,2,i_uyz) &
1658  * 0.5_rp * ( dens(ks+1,i,j)+dens(ks+1,i+1,j) ) )
1659  vel = vel * j23g(ks+1,i,j)
1660  flux(ks+1,i,j) = vel / mapf(i,j,+1) &
1661  * ( ( 2.0_rp * val(ks+2,i,j) + 5.0_rp * val(ks+1,i,j) - val(ks,i,j) ) / 6.0_rp &
1662  * ( 0.5_rp + sign(0.5_rp,vel) ) &
1663  + ( - 3.0_rp * val(ks,i,j) &
1664  + 27.0_rp * val(ks+1,i,j) &
1665  + 47.0_rp * val(ks+2,i,j) &
1666  - 13.0_rp * val(ks+3,i,j) &
1667  + 2.0_rp * val(ks+4,i,j) ) / 60.0_rp &
1668  * ( 0.5_rp - sign(0.5_rp,vel) ) )
1669 
1670  vel = ( f2h(ke-2,1,i_uyz) &
1671  * 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) ) &
1672  + f2h(ke-2,2,i_uyz) &
1673  * 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) ) ) &
1674  / ( f2h(ke-2,1,i_uyz) &
1675  * 0.5_rp * ( dens(ke-1,i,j)+dens(ke-1,i+1,j) ) &
1676  + f2h(ke-2,2,i_uyz) &
1677  * 0.5_rp * ( dens(ke-2,i,j)+dens(ke-2,i+1,j) ) )
1678  vel = vel * j23g(ke-2,i,j)
1679  flux(ke-2,i,j) = vel / mapf(i,j,+1) &
1680  * ( ( - 3.0_rp * val(ke,i,j) &
1681  + 27.0_rp * val(ke-1,i,j) &
1682  + 47.0_rp * val(ke-2,i,j) &
1683  - 13.0_rp * val(ke-3,i,j) &
1684  + 2.0_rp * val(ke-4,i,j) ) / 60.0_rp &
1685  * ( 0.5_rp + sign(0.5_rp,vel) ) &
1686  + ( 2.0_rp * val(ke-2,i,j) + 5.0_rp * val(ke-1,i,j) - val(ke,i,j) ) / 6.0_rp &
1687  * ( 0.5_rp - sign(0.5_rp,vel) ) )
1688 
1689  flux(ke ,i,j) = 0.0_rp
1690  enddo
1691  enddo
1692  !$acc end kernels
1693  !$omp end do nowait
1694 
1695 
1696  end if
1697 
1698 
1699  !$acc end data
1700 
1701  !$omp end parallel
1702  return
1703  end subroutine atmos_dyn_fvm_fluxj23_uyz_ud5
1704 
1705  !-----------------------------------------------------------------------------
1707  subroutine atmos_dyn_fvm_fluxx_uyz_ud5( &
1708  flux, &
1709  mom, val, DENS, &
1710  GSQRT, MAPF, &
1711  num_diff, &
1712  CDZ, TwoD, &
1713  IIS, IIE, JJS, JJE )
1714  implicit none
1715 
1716  real(rp), intent(inout) :: flux (ka,ia,ja)
1717  real(rp), intent(in) :: mom (ka,ia,ja)
1718  real(rp), intent(in) :: val (ka,ia,ja)
1719  real(rp), intent(in) :: dens (ka,ia,ja)
1720  real(rp), intent(in) :: gsqrt (ka,ia,ja)
1721  real(rp), intent(in) :: mapf ( ia,ja,2)
1722  real(rp), intent(in) :: num_diff(ka,ia,ja)
1723  real(rp), intent(in) :: cdz (ka)
1724  logical, intent(in) :: twod
1725  integer, intent(in) :: iis, iie, jjs, jje
1726 
1727  real(rp) :: vel
1728  integer :: k, i, j
1729  !---------------------------------------------------------------------------
1730 
1731  ! note that x-index is added by -1
1732 
1733 
1734 
1735  !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
1736  !$omp private(vel) &
1737  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,GSQRT,MAPF,num_diff)
1738  !$acc kernels
1739  do j = jjs, jje
1740  do i = iis, iie+1
1741  do k = ks, ke
1742 #ifdef DEBUG
1743  call check( __line__, mom(k,i ,j) )
1744  call check( __line__, mom(k,i-1,j) )
1745 
1746  call check( __line__, val(k,i-1,j) )
1747  call check( __line__, val(k,i,j) )
1748 
1749  call check( __line__, val(k,i-2,j) )
1750  call check( __line__, val(k,i+1,j) )
1751 
1752  call check( __line__, val(k,i-3,j) )
1753  call check( __line__, val(k,i+2,j) )
1754 
1755 #endif
1756  vel = ( 0.5_rp * ( mom(k,i,j)+mom(k,i-1,j) ) ) &
1757  / ( dens(k,i,j) )
1758  flux(k,i-1,j) = gsqrt(k,i,j) / mapf(i,j,+2) * vel &
1759  * ( ( f51 * ( val(k,i+2,j)+val(k,i-3,j) ) &
1760  + f52 * ( val(k,i+1,j)+val(k,i-2,j) ) &
1761  + f53 * ( val(k,i,j)+val(k,i-1,j) ) ) &
1762  - ( f51 * ( val(k,i+2,j)-val(k,i-3,j) ) &
1763  + f54 * ( val(k,i+1,j)-val(k,i-2,j) ) &
1764  + f55 * ( val(k,i,j)-val(k,i-1,j) ) ) * sign(1.0_rp,vel) ) &
1765  + gsqrt(k,i,j) * num_diff(k,i,j)
1766  enddo
1767  enddo
1768  enddo
1769  !$acc end kernels
1770 #ifdef DEBUG
1771  k = iundef; i = iundef; j = iundef
1772 #endif
1773 
1774 
1775 
1776  return
1777  end subroutine atmos_dyn_fvm_fluxx_uyz_ud5
1778 
1779  !-----------------------------------------------------------------------------
1781  subroutine atmos_dyn_fvm_fluxy_uyz_ud5( &
1782  flux, &
1783  mom, val, DENS, &
1784  GSQRT, MAPF, &
1785  num_diff, &
1786  CDZ, TwoD, &
1787  IIS, IIE, JJS, JJE )
1788  implicit none
1789 
1790  real(rp), intent(inout) :: flux (ka,ia,ja)
1791  real(rp), intent(in) :: mom (ka,ia,ja)
1792  real(rp), intent(in) :: val (ka,ia,ja)
1793  real(rp), intent(in) :: dens (ka,ia,ja)
1794  real(rp), intent(in) :: gsqrt (ka,ia,ja)
1795  real(rp), intent(in) :: mapf ( ia,ja,2)
1796  real(rp), intent(in) :: num_diff(ka,ia,ja)
1797  real(rp), intent(in) :: cdz (ka)
1798  logical, intent(in) :: twod
1799  integer, intent(in) :: iis, iie, jjs, jje
1800 
1801  real(rp) :: vel
1802  integer :: k, i, j
1803  !---------------------------------------------------------------------------
1804 
1805 
1806 
1807  if ( twod ) then
1808 
1809  !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ &
1810  !$omp private(vel) &
1811  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,GSQRT,MAPF,num_diff,TwoD)
1812  !$acc kernels
1813  do j = jjs-1, jje
1814  do k = ks, ke
1815  i = iis
1816 #ifdef DEBUG
1817  call check( __line__, mom(k,i ,j) )
1818 
1819  call check( __line__, val(k,i,j) )
1820  call check( __line__, val(k,i,j+1) )
1821 
1822  call check( __line__, val(k,i,j-1) )
1823  call check( __line__, val(k,i,j+2) )
1824 
1825  call check( __line__, val(k,i,j-2) )
1826  call check( __line__, val(k,i,j+3) )
1827 
1828 #endif
1829  vel = ( mom(k,i,j) ) &
1830  / ( 0.5_rp * ( dens(k,i,j)+dens(k,i,j+1) ) )
1831  flux(k,i,j) = gsqrt(k,i,j) / mapf(i,j,+1) * vel &
1832  * ( ( f51 * ( val(k,i,j+3)+val(k,i,j-2) ) &
1833  + f52 * ( val(k,i,j+2)+val(k,i,j-1) ) &
1834  + f53 * ( val(k,i,j+1)+val(k,i,j) ) ) &
1835  - ( f51 * ( val(k,i,j+3)-val(k,i,j-2) ) &
1836  + f54 * ( val(k,i,j+2)-val(k,i,j-1) ) &
1837  + f55 * ( val(k,i,j+1)-val(k,i,j) ) ) * sign(1.0_rp,vel) ) &
1838  + gsqrt(k,i,j) * num_diff(k,i,j)
1839  enddo
1840  enddo
1841  !$acc end kernels
1842 #ifdef DEBUG
1843  k = iundef; i = iundef; j = iundef
1844 #endif
1845 
1846  else
1847 
1848 
1849  !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
1850  !$omp private(vel) &
1851  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,GSQRT,MAPF,num_diff)
1852  !$acc kernels
1853  do j = jjs-1, jje
1854  do i = iis, iie
1855  do k = ks, ke
1856 #ifdef DEBUG
1857  call check( __line__, mom(k,i ,j) )
1858  call check( __line__, mom(k,i-1,j) )
1859 
1860  call check( __line__, val(k,i,j) )
1861  call check( __line__, val(k,i,j+1) )
1862 
1863  call check( __line__, val(k,i,j-1) )
1864  call check( __line__, val(k,i,j+2) )
1865 
1866  call check( __line__, val(k,i,j-2) )
1867  call check( __line__, val(k,i,j+3) )
1868 
1869 #endif
1870  vel = ( 0.5_rp * ( mom(k,i,j)+mom(k,i+1,j) ) ) &
1871  / ( 0.25_rp * ( dens(k,i,j)+dens(k,i+1,j)+dens(k,i,j+1)+dens(k,i+1,j+1) ) )
1872  flux(k,i,j) = gsqrt(k,i,j) / mapf(i,j,+1) * vel &
1873  * ( ( f51 * ( val(k,i,j+3)+val(k,i,j-2) ) &
1874  + f52 * ( val(k,i,j+2)+val(k,i,j-1) ) &
1875  + f53 * ( val(k,i,j+1)+val(k,i,j) ) ) &
1876  - ( f51 * ( val(k,i,j+3)-val(k,i,j-2) ) &
1877  + f54 * ( val(k,i,j+2)-val(k,i,j-1) ) &
1878  + f55 * ( val(k,i,j+1)-val(k,i,j) ) ) * sign(1.0_rp,vel) ) &
1879  + gsqrt(k,i,j) * num_diff(k,i,j)
1880  enddo
1881  enddo
1882  enddo
1883  !$acc end kernels
1884 #ifdef DEBUG
1885  k = iundef; i = iundef; j = iundef
1886 #endif
1887 
1888 
1889  end if
1890 
1891 
1892  return
1893  end subroutine atmos_dyn_fvm_fluxy_uyz_ud5
1894 
1895 
1896 
1897  !-----------------------------------------------------------------------------
1899  subroutine atmos_dyn_fvm_fluxz_xvz_ud5( &
1900  flux, &
1901  mom, val, DENS, &
1902  GSQRT, J33G, &
1903  num_diff, &
1904  CDZ, TwoD, &
1905  IIS, IIE, JJS, JJE )
1906  implicit none
1907 
1908  real(rp), intent(inout) :: flux (ka,ia,ja)
1909  real(rp), intent(in) :: mom (ka,ia,ja)
1910  real(rp), intent(in) :: val (ka,ia,ja)
1911  real(rp), intent(in) :: dens (ka,ia,ja)
1912  real(rp), intent(in) :: gsqrt (ka,ia,ja)
1913  real(rp), intent(in) :: j33g
1914  real(rp), intent(in) :: num_diff(ka,ia,ja)
1915  real(rp), intent(in) :: cdz (ka)
1916  logical, intent(in) :: twod
1917  integer, intent(in) :: iis, iie, jjs, jje
1918 
1919  real(rp) :: vel
1920  integer :: k, i, j
1921  !---------------------------------------------------------------------------
1922 
1923  !$omp parallel default(none) private(i,j,k,vel) &
1924  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,J33G,GSQRT,num_diff) &
1925  !$omp shared(CDZ,TwoD)
1926 
1927  !$acc data copy(flux) copyin(mom, val, DENS, GSQRT, num_diff, CDZ)
1928 
1929 
1930  !$omp do OMP_SCHEDULE_ collapse(2)
1931  !$acc kernels
1932  do j = jjs, jje
1933  do i = iis, iie
1934  do k = ks+2, ke-3
1935 #ifdef DEBUG
1936  call check( __line__, mom(k,i,j) )
1937  call check( __line__, mom(k,i,j+1) )
1938 
1939  call check( __line__, val(k,i,j) )
1940  call check( __line__, val(k+1,i,j) )
1941 
1942  call check( __line__, val(k-1,i,j) )
1943  call check( __line__, val(k+2,i,j) )
1944 
1945  call check( __line__, val(k-2,i,j) )
1946  call check( __line__, val(k+3,i,j) )
1947 
1948 #endif
1949  vel = ( 0.5_rp * ( mom(k,i,j)+mom(k,i,j+1) ) ) &
1950  / ( f2h(k,1,i_xvz) &
1951  * 0.5_rp * ( dens(k+1,i,j)+dens(k+1,i,j+1) ) &
1952  + f2h(k,2,i_xvz) &
1953  * 0.5_rp * ( dens(k,i,j)+dens(k,i,j+1) ) )
1954  flux(k,i,j) = j33g * vel &
1955  * ( ( f51 * ( val(k+3,i,j)+val(k-2,i,j) ) &
1956  + f52 * ( val(k+2,i,j)+val(k-1,i,j) ) &
1957  + f53 * ( val(k+1,i,j)+val(k,i,j) ) ) &
1958  - ( f51 * ( val(k+3,i,j)-val(k-2,i,j) ) &
1959  + f54 * ( val(k+2,i,j)-val(k-1,i,j) ) &
1960  + f55 * ( val(k+1,i,j)-val(k,i,j) ) ) * sign(1.0_rp,vel) ) &
1961  + gsqrt(k,i,j) * num_diff(k,i,j)
1962  enddo
1963  enddo
1964  enddo
1965  !$acc end kernels
1966  !$omp end do nowait
1967 #ifdef DEBUG
1968  k = iundef; i = iundef; j = iundef
1969 #endif
1970 
1971  !$omp do OMP_SCHEDULE_ collapse(2)
1972  !$acc kernels
1973  do j = jjs, jje
1974  do i = iis, iie
1975 #ifdef DEBUG
1976 
1977  call check( __line__, mom(ks,i ,j) )
1978  call check( __line__, mom(ks,i,j+1) )
1979  call check( __line__, val(ks+1,i,j) )
1980  call check( __line__, val(ks,i,j) )
1981 
1982  call check( __line__, mom(ks+1,i ,j) )
1983  call check( __line__, mom(ks+1,i,j+1) )
1984  call check( __line__, val(ks+3,i,j) )
1985  call check( __line__, val(ks+2,i,j) )
1986 
1987 #endif
1988  ! The boundary condition is qflx_hi + qflxJ13 + qfluxJ23 = 0 at KS-1.
1989  ! The flux at KS-1 can be non-zero.
1990  ! To reduce calculations, all the fluxes are set to zero.
1991  flux(ks-1,i,j) = 0.0_rp
1992 
1993  vel = ( 0.5_rp * ( mom(ks,i,j)+mom(ks,i,j+1) ) ) &
1994  / ( f2h(ks,1,i_xvz) &
1995  * 0.5_rp * ( dens(ks+1,i,j)+dens(ks+1,i,j+1) ) &
1996  + f2h(ks,2,i_xvz) &
1997  * 0.5_rp * ( dens(ks,i,j)+dens(ks,i,j+1) ) )
1998  flux(ks,i,j) = j33g * vel &
1999  * ( f2 * ( val(ks+1,i,j)+val(ks,i,j) ) &
2000  * ( 0.5_rp + sign(0.5_rp,vel) ) &
2001  + ( 2.0_rp * val(ks,i,j) + 5.0_rp * val(ks+1,i,j) - val(ks+2,i,j) ) / 6.0_rp &
2002  * ( 0.5_rp - sign(0.5_rp,vel) ) ) &
2003  + gsqrt(ks,i,j) * num_diff(ks,i,j)
2004  vel = ( 0.5_rp * ( mom(ke-1,i,j)+mom(ke-1,i,j+1) ) ) &
2005  / ( f2h(ke-1,1,i_xvz) &
2006  * 0.5_rp * ( dens(ke,i,j)+dens(ke,i,j+1) ) &
2007  + f2h(ke-1,2,i_xvz) &
2008  * 0.5_rp * ( dens(ke-1,i,j)+dens(ke-1,i,j+1) ) )
2009  flux(ke-1,i,j) = j33g * vel &
2010  * ( ( 2.0_rp * val(ke,i,j) + 5.0_rp * val(ke-1,i,j) - val(ke-2,i,j) ) / 6.0_rp &
2011  * ( 0.5_rp + sign(0.5_rp,vel) ) &
2012  + f2 * ( val(ke,i,j)+val(ke-1,i,j) ) &
2013  * ( 0.5_rp - sign(0.5_rp,vel) ) ) &
2014  + gsqrt(ke-1,i,j) * num_diff(ke-1,i,j)
2015 
2016  vel = ( 0.5_rp * ( mom(ks+1,i,j)+mom(ks+1,i,j+1) ) ) &
2017  / ( f2h(ks+1,1,i_xvz) &
2018  * 0.5_rp * ( dens(ks+2,i,j)+dens(ks+2,i,j+1) ) &
2019  + f2h(ks+1,2,i_xvz) &
2020  * 0.5_rp * ( dens(ks+1,i,j)+dens(ks+1,i,j+1) ) )
2021  flux(ks+1,i,j) = j33g * vel &
2022  * ( ( 2.0_rp * val(ks+2,i,j) + 5.0_rp * val(ks+1,i,j) - val(ks,i,j) ) / 6.0_rp &
2023  * ( 0.5_rp + sign(0.5_rp,vel) ) &
2024  + ( - 3.0_rp * val(ks,i,j) &
2025  + 27.0_rp * val(ks+1,i,j) &
2026  + 47.0_rp * val(ks+2,i,j) &
2027  - 13.0_rp * val(ks+3,i,j) &
2028  + 2.0_rp * val(ks+4,i,j) ) / 60.0_rp &
2029  * ( 0.5_rp - sign(0.5_rp,vel) ) ) &
2030  + gsqrt(ks+1,i,j) * num_diff(ks+1,i,j)
2031  vel = ( 0.5_rp * ( mom(ke-2,i,j)+mom(ke-2,i,j+1) ) ) &
2032  / ( f2h(ke-2,1,i_xvz) &
2033  * 0.5_rp * ( dens(ke-1,i,j)+dens(ke-1,i,j+1) ) &
2034  + f2h(ke-2,2,i_xvz) &
2035  * 0.5_rp * ( dens(ke-2,i,j)+dens(ke-2,i,j+1) ) )
2036  flux(ke-2,i,j) = j33g * vel &
2037  * ( ( - 3.0_rp * val(ke,i,j) &
2038  + 27.0_rp * val(ke-1,i,j) &
2039  + 47.0_rp * val(ke-2,i,j) &
2040  - 13.0_rp * val(ke-3,i,j) &
2041  + 2.0_rp * val(ke-4,i,j) ) / 60.0_rp &
2042  * ( 0.5_rp + sign(0.5_rp,vel) ) &
2043  + ( 2.0_rp * val(ke-2,i,j) + 5.0_rp * val(ke-1,i,j) - val(ke,i,j) ) / 6.0_rp &
2044  * ( 0.5_rp - sign(0.5_rp,vel) ) ) &
2045  + gsqrt(ke-2,i,j) * num_diff(ke-2,i,j)
2046 
2047  flux(ke,i,j) = 0.0_rp
2048  enddo
2049  enddo
2050  !$acc end kernels
2051  !$omp end do nowait
2052 
2053 
2054  !$acc end data
2055 
2056  !$omp end parallel
2057 #ifdef DEBUG
2058  k = iundef; i = iundef; j = iundef
2059 #endif
2060 
2061  return
2062  end subroutine atmos_dyn_fvm_fluxz_xvz_ud5
2063 
2064  !-----------------------------------------------------------------------------
2066  subroutine atmos_dyn_fvm_fluxj13_xvz_ud5( &
2067  flux, &
2068  mom, val, DENS, &
2069  GSQRT, J13G, MAPF, &
2070  CDZ, TwoD, &
2071  IIS, IIE, JJS, JJE )
2072  implicit none
2073 
2074  real(rp), intent(inout) :: flux (ka,ia,ja)
2075  real(rp), intent(in) :: mom (ka,ia,ja)
2076  real(rp), intent(in) :: val (ka,ia,ja)
2077  real(rp), intent(in) :: dens (ka,ia,ja)
2078  real(rp), intent(in) :: gsqrt (ka,ia,ja)
2079  real(rp), intent(in) :: j13g (ka,ia,ja)
2080  real(rp), intent(in) :: mapf ( ia,ja,2)
2081  real(rp), intent(in) :: cdz (ka)
2082  logical, intent(in) :: twod
2083  integer, intent(in) :: iis, iie, jjs, jje
2084 
2085  real(rp) :: vel
2086  integer :: k, i, j
2087  !---------------------------------------------------------------------------
2088 
2089  !$omp parallel default(none) private(i,j,k,vel) &
2090  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,J13G,MAPF) &
2091  !$omp shared(GSQRT,CDZ,TwoD)
2092 
2093  !$acc data copy(flux) copyin(mom, val, DENS, GSQRT, J13G, MAPF, CDZ)
2094 
2095 
2096 
2097  !$omp do OMP_SCHEDULE_ collapse(2)
2098  !$acc kernels
2099  do j = jjs, jje
2100  do i = iis, iie
2101  do k = ks+2, ke-3
2102  vel = ( f2h(k,1,i_xvz) &
2103  * 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) ) &
2104  + f2h(k,2,i_xvz) &
2105  * 0.25_rp * ( mom(k,i,j)+mom(k,i-1,j)+mom(k,i,j+1)+mom(k,i-1,j+1) ) ) &
2106  / ( f2h(k,1,i_xvz) &
2107  * 0.5_rp * ( dens(k+1,i,j)+dens(k+1,i,j+1) ) &
2108  + f2h(k,2,i_xvz) &
2109  * 0.5_rp * ( dens(k,i,j)+dens(k,i,j+1) ) )
2110  vel = vel * j13g(k,i,j)
2111  flux(k,i,j) = vel / mapf(i,j,+2) &
2112  * ( ( f51 * ( val(k+3,i,j)+val(k-2,i,j) ) &
2113  + f52 * ( val(k+2,i,j)+val(k-1,i,j) ) &
2114  + f53 * ( val(k+1,i,j)+val(k,i,j) ) ) &
2115  - ( f51 * ( val(k+3,i,j)-val(k-2,i,j) ) &
2116  + f54 * ( val(k+2,i,j)-val(k-1,i,j) ) &
2117  + f55 * ( val(k+1,i,j)-val(k,i,j) ) ) * sign(1.0_rp,vel) )
2118  enddo
2119  enddo
2120  enddo
2121  !$acc end kernels
2122  !$omp end do nowait
2123 
2124  !$omp do OMP_SCHEDULE_ collapse(2)
2125  !$acc kernels
2126  do j = jjs, jje
2127  do i = iis, iie
2128  ! The boundary condition is qflx_hi + qflxJ13 + qfluxJ23 = 0 at KS-1.
2129  ! The flux at KS-1 can be non-zero.
2130  ! To reduce calculations, all the fluxes are set to zero.
2131  flux(ks-1,i,j) = 0.0_rp
2132 
2133  vel = ( f2h(ks,1,i_xvz) &
2134  * 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) ) &
2135  + f2h(ks,2,i_xvz) &
2136  * 0.25_rp * ( mom(ks,i,j)+mom(ks,i-1,j)+mom(ks,i,j+1)+mom(ks,i-1,j+1) ) ) &
2137  / ( f2h(ks,1,i_xvz) &
2138  * 0.5_rp * ( dens(ks+1,i,j)+dens(ks+1,i,j+1) ) &
2139  + f2h(ks,2,i_xvz) &
2140  * 0.5_rp * ( dens(ks,i,j)+dens(ks,i,j+1) ) )
2141  vel = vel * j13g(ks,i,j)
2142  flux(ks,i,j) = vel / mapf(i,j,+2) &
2143  * ( f2 * ( val(ks+1,i,j)+val(ks,i,j) ) &
2144  * ( 0.5_rp + sign(0.5_rp,vel) ) &
2145  + ( 2.0_rp * val(ks,i,j) + 5.0_rp * val(ks+1,i,j) - val(ks+2,i,j) ) / 6.0_rp &
2146  * ( 0.5_rp - sign(0.5_rp,vel) ) )
2147 
2148  vel = ( f2h(ke-1,1,i_xvz) &
2149  * 0.25_rp * ( mom(ke,i,j)+mom(ke,i-1,j)+mom(ke,i,j+1)+mom(ke,i-1,j+1) ) &
2150  + f2h(ke-1,2,i_xvz) &
2151  * 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) ) ) &
2152  / ( f2h(ke-1,1,i_xvz) &
2153  * 0.5_rp * ( dens(ke,i,j)+dens(ke,i,j+1) ) &
2154  + f2h(ke-1,2,i_xvz) &
2155  * 0.5_rp * ( dens(ke-1,i,j)+dens(ke-1,i,j+1) ) )
2156  vel = vel * j13g(ke-1,i,j)
2157  flux(ke-1,i,j) = vel / mapf(i,j,+2) &
2158  * ( ( 2.0_rp * val(ke,i,j) + 5.0_rp * val(ke-1,i,j) - val(ke-2,i,j) ) / 6.0_rp &
2159  * ( 0.5_rp + sign(0.5_rp,vel) ) &
2160  + f2 * ( val(ke,i,j)+val(ke-1,i,j) ) &
2161  * ( 0.5_rp - sign(0.5_rp,vel) ) )
2162 
2163  vel = ( f2h(ks+1,1,i_xvz) &
2164  * 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) ) &
2165  + f2h(ks+1,2,i_xvz) &
2166  * 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) ) ) &
2167  / ( f2h(ks+1,1,i_xvz) &
2168  * 0.5_rp * ( dens(ks+2,i,j)+dens(ks+2,i,j+1) ) &
2169  + f2h(ks+1,2,i_xvz) &
2170  * 0.5_rp * ( dens(ks+1,i,j)+dens(ks+1,i,j+1) ) )
2171  vel = vel * j13g(ks+1,i,j)
2172  flux(ks+1,i,j) = vel / mapf(i,j,+2) &
2173  * ( ( 2.0_rp * val(ks+2,i,j) + 5.0_rp * val(ks+1,i,j) - val(ks,i,j) ) / 6.0_rp &
2174  * ( 0.5_rp + sign(0.5_rp,vel) ) &
2175  + ( - 3.0_rp * val(ks,i,j) &
2176  + 27.0_rp * val(ks+1,i,j) &
2177  + 47.0_rp * val(ks+2,i,j) &
2178  - 13.0_rp * val(ks+3,i,j) &
2179  + 2.0_rp * val(ks+4,i,j) ) / 60.0_rp &
2180  * ( 0.5_rp - sign(0.5_rp,vel) ) )
2181 
2182  vel = ( f2h(ke-2,1,i_xvz) &
2183  * 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) ) &
2184  + f2h(ke-2,2,i_xvz) &
2185  * 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) ) ) &
2186  / ( f2h(ke-2,1,i_xvz) &
2187  * 0.5_rp * ( dens(ke-1,i,j)+dens(ke-1,i,j+1) ) &
2188  + f2h(ke-2,2,i_xvz) &
2189  * 0.5_rp * ( dens(ke-2,i,j)+dens(ke-2,i,j+1) ) )
2190  vel = vel * j13g(ke-2,i,j)
2191  flux(ke-2,i,j) = vel / mapf(i,j,+2) &
2192  * ( ( - 3.0_rp * val(ke,i,j) &
2193  + 27.0_rp * val(ke-1,i,j) &
2194  + 47.0_rp * val(ke-2,i,j) &
2195  - 13.0_rp * val(ke-3,i,j) &
2196  + 2.0_rp * val(ke-4,i,j) ) / 60.0_rp &
2197  * ( 0.5_rp + sign(0.5_rp,vel) ) &
2198  + ( 2.0_rp * val(ke-2,i,j) + 5.0_rp * val(ke-1,i,j) - val(ke,i,j) ) / 6.0_rp &
2199  * ( 0.5_rp - sign(0.5_rp,vel) ) )
2200 
2201  flux(ke ,i,j) = 0.0_rp
2202  enddo
2203  enddo
2204  !$acc end kernels
2205  !$omp end do nowait
2206 
2207 
2208 
2209  !$acc end data
2210 
2211  !$omp end parallel
2212  return
2213  end subroutine atmos_dyn_fvm_fluxj13_xvz_ud5
2214 
2215  !-----------------------------------------------------------------------------
2217  subroutine atmos_dyn_fvm_fluxj23_xvz_ud5( &
2218  flux, &
2219  mom, val, DENS, &
2220  GSQRT, J23G, MAPF, &
2221  CDZ, TwoD, &
2222  IIS, IIE, JJS, JJE )
2223  implicit none
2224 
2225  real(rp), intent(inout) :: flux (ka,ia,ja)
2226  real(rp), intent(in) :: mom (ka,ia,ja)
2227  real(rp), intent(in) :: val (ka,ia,ja)
2228  real(rp), intent(in) :: dens (ka,ia,ja)
2229  real(rp), intent(in) :: gsqrt (ka,ia,ja)
2230  real(rp), intent(in) :: j23g (ka,ia,ja)
2231  real(rp), intent(in) :: mapf ( ia,ja,2)
2232  real(rp), intent(in) :: cdz (ka)
2233  logical, intent(in) :: twod
2234  integer, intent(in) :: iis, iie, jjs, jje
2235 
2236  real(rp) :: vel
2237  integer :: k, i, j
2238  !---------------------------------------------------------------------------
2239 
2240  !$omp parallel default(none) private(i,j,k,vel) &
2241  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,J23G,MAPF) &
2242  !$omp shared(GSQRT,CDZ,TwoD)
2243 
2244  !$acc data copy(flux) copyin(mom, val, DENS, GSQRT, J23G, MAPF, CDZ)
2245 
2246 
2247 
2248  !$omp do OMP_SCHEDULE_ collapse(2)
2249  !$acc kernels
2250  do j = jjs, jje
2251  do i = iis, iie
2252  do k = ks+2, ke-3
2253  vel = ( f2h(k,1,i_xvz) &
2254  * mom(k+1,i,j) &
2255  + f2h(k,2,i_xvz) &
2256  * mom(k,i,j) ) &
2257  / ( f2h(k,1,i_xvz) &
2258  * 0.5_rp * ( dens(k+1,i,j)+dens(k+1,i,j+1) ) &
2259  + f2h(k,2,i_xvz) &
2260  * 0.5_rp * ( dens(k,i,j)+dens(k,i,j+1) ) )
2261  vel = vel * j23g(k,i,j)
2262  flux(k,i,j) = vel / mapf(i,j,+1) &
2263  * ( ( f51 * ( val(k+3,i,j)+val(k-2,i,j) ) &
2264  + f52 * ( val(k+2,i,j)+val(k-1,i,j) ) &
2265  + f53 * ( val(k+1,i,j)+val(k,i,j) ) ) &
2266  - ( f51 * ( val(k+3,i,j)-val(k-2,i,j) ) &
2267  + f54 * ( val(k+2,i,j)-val(k-1,i,j) ) &
2268  + f55 * ( val(k+1,i,j)-val(k,i,j) ) ) * sign(1.0_rp,vel) )
2269  enddo
2270  enddo
2271  enddo
2272  !$acc end kernels
2273  !$omp end do nowait
2274 
2275  !$omp do OMP_SCHEDULE_ collapse(2)
2276  !$acc kernels
2277  do j = jjs, jje
2278  do i = iis, iie
2279  ! The boundary condition is qflx_hi + qflxJ13 + qfluxJ23 = 0 at KS-1.
2280  ! The flux at KS-1 can be non-zero.
2281  ! To reduce calculations, all the fluxes are set to zero.
2282  flux(ks-1,i,j) = 0.0_rp
2283 
2284  vel = ( f2h(ks,1,i_xvz) &
2285  * mom(ks+1,i,j) &
2286  + f2h(ks,2,i_xvz) &
2287  * mom(ks,i,j) ) &
2288  / ( f2h(ks,1,i_xvz) &
2289  * 0.5_rp * ( dens(ks+1,i,j)+dens(ks+1,i,j+1) ) &
2290  + f2h(ks,2,i_xvz) &
2291  * 0.5_rp * ( dens(ks,i,j)+dens(ks,i,j+1) ) )
2292  vel = vel * j23g(ks,i,j)
2293  flux(ks,i,j) = vel / mapf(i,j,+1) &
2294  * ( f2 * ( val(ks+1,i,j)+val(ks,i,j) ) &
2295  * ( 0.5_rp + sign(0.5_rp,vel) ) &
2296  + ( 2.0_rp * val(ks,i,j) + 5.0_rp * val(ks+1,i,j) - val(ks+2,i,j) ) / 6.0_rp &
2297  * ( 0.5_rp - sign(0.5_rp,vel) ) )
2298 
2299  vel = ( f2h(ke-1,1,i_xvz) &
2300  * mom(ke,i,j) &
2301  + f2h(ke-1,2,i_xvz) &
2302  * mom(ke-1,i,j) ) &
2303  / ( f2h(ke-1,1,i_xvz) &
2304  * 0.5_rp * ( dens(ke,i,j)+dens(ke,i,j+1) ) &
2305  + f2h(ke-1,2,i_xvz) &
2306  * 0.5_rp * ( dens(ke-1,i,j)+dens(ke-1,i,j+1) ) )
2307  vel = vel * j23g(ke-1,i,j)
2308  flux(ke-1,i,j) = vel / mapf(i,j,+1) &
2309  * ( ( 2.0_rp * val(ke,i,j) + 5.0_rp * val(ke-1,i,j) - val(ke-2,i,j) ) / 6.0_rp &
2310  * ( 0.5_rp + sign(0.5_rp,vel) ) &
2311  + f2 * ( val(ke,i,j)+val(ke-1,i,j) ) &
2312  * ( 0.5_rp - sign(0.5_rp,vel) ) )
2313 
2314  vel = ( f2h(ks+1,1,i_xvz) &
2315  * mom(ks+2,i,j) &
2316  + f2h(ks+1,2,i_xvz) &
2317  * mom(ks+1,i,j) ) &
2318  / ( f2h(ks+1,1,i_xvz) &
2319  * 0.5_rp * ( dens(ks+2,i,j)+dens(ks+2,i,j+1) ) &
2320  + f2h(ks+1,2,i_xvz) &
2321  * 0.5_rp * ( dens(ks+1,i,j)+dens(ks+1,i,j+1) ) )
2322  vel = vel * j23g(ks+1,i,j)
2323  flux(ks+1,i,j) = vel / mapf(i,j,+1) &
2324  * ( ( 2.0_rp * val(ks+2,i,j) + 5.0_rp * val(ks+1,i,j) - val(ks,i,j) ) / 6.0_rp &
2325  * ( 0.5_rp + sign(0.5_rp,vel) ) &
2326  + ( - 3.0_rp * val(ks,i,j) &
2327  + 27.0_rp * val(ks+1,i,j) &
2328  + 47.0_rp * val(ks+2,i,j) &
2329  - 13.0_rp * val(ks+3,i,j) &
2330  + 2.0_rp * val(ks+4,i,j) ) / 60.0_rp &
2331  * ( 0.5_rp - sign(0.5_rp,vel) ) )
2332 
2333  vel = ( f2h(ke-2,1,i_xvz) &
2334  * mom(ke-1,i,j) &
2335  + f2h(ke-2,2,i_xvz) &
2336  * mom(ke-2,i,j) ) &
2337  / ( f2h(ke-2,1,i_xvz) &
2338  * 0.5_rp * ( dens(ke-1,i,j)+dens(ke-1,i,j+1) ) &
2339  + f2h(ke-2,2,i_xvz) &
2340  * 0.5_rp * ( dens(ke-2,i,j)+dens(ke-2,i,j+1) ) )
2341  vel = vel * j23g(ke-2,i,j)
2342  flux(ke-2,i,j) = vel / mapf(i,j,+1) &
2343  * ( ( - 3.0_rp * val(ke,i,j) &
2344  + 27.0_rp * val(ke-1,i,j) &
2345  + 47.0_rp * val(ke-2,i,j) &
2346  - 13.0_rp * val(ke-3,i,j) &
2347  + 2.0_rp * val(ke-4,i,j) ) / 60.0_rp &
2348  * ( 0.5_rp + sign(0.5_rp,vel) ) &
2349  + ( 2.0_rp * val(ke-2,i,j) + 5.0_rp * val(ke-1,i,j) - val(ke,i,j) ) / 6.0_rp &
2350  * ( 0.5_rp - sign(0.5_rp,vel) ) )
2351 
2352  flux(ke ,i,j) = 0.0_rp
2353  enddo
2354  enddo
2355  !$acc end kernels
2356  !$omp end do nowait
2357 
2358 
2359 
2360  !$acc end data
2361 
2362  !$omp end parallel
2363  return
2364  end subroutine atmos_dyn_fvm_fluxj23_xvz_ud5
2365 
2366  !-----------------------------------------------------------------------------
2368  subroutine atmos_dyn_fvm_fluxx_xvz_ud5( &
2369  flux, &
2370  mom, val, DENS, &
2371  GSQRT, MAPF, &
2372  num_diff, &
2373  CDZ, TwoD, &
2374  IIS, IIE, JJS, JJE )
2375  implicit none
2376 
2377  real(rp), intent(inout) :: flux (ka,ia,ja)
2378  real(rp), intent(in) :: mom (ka,ia,ja)
2379  real(rp), intent(in) :: val (ka,ia,ja)
2380  real(rp), intent(in) :: dens (ka,ia,ja)
2381  real(rp), intent(in) :: gsqrt (ka,ia,ja)
2382  real(rp), intent(in) :: mapf ( ia,ja,2)
2383  real(rp), intent(in) :: num_diff(ka,ia,ja)
2384  real(rp), intent(in) :: cdz (ka)
2385  logical, intent(in) :: twod
2386  integer, intent(in) :: iis, iie, jjs, jje
2387 
2388  real(rp) :: vel
2389  integer :: k, i, j
2390  !---------------------------------------------------------------------------
2391 
2392 
2393 
2394  !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
2395  !$omp private(vel) &
2396  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,GSQRT,MAPF,num_diff)
2397  !$acc kernels
2398  do j = jjs, jje
2399  do i = iis-1, iie
2400  do k = ks, ke
2401 #ifdef DEBUG
2402  call check( __line__, mom(k,i ,j) )
2403  call check( __line__, mom(k,i,j-1) )
2404 
2405  call check( __line__, val(k,i,j) )
2406  call check( __line__, val(k,i+1,j) )
2407 
2408  call check( __line__, val(k,i-1,j) )
2409  call check( __line__, val(k,i+2,j) )
2410 
2411  call check( __line__, val(k,i-2,j) )
2412  call check( __line__, val(k,i+3,j) )
2413 
2414 #endif
2415  vel = ( 0.5_rp * ( mom(k,i,j)+mom(k,i,j+1) ) ) &
2416  / ( 0.25_rp * ( dens(k,i,j)+dens(k,i+1,j)+dens(k,i,j+1)+dens(k,i+1,j+1) ) )
2417  flux(k,i,j) = gsqrt(k,i,j) / mapf(i,j,+2) * vel &
2418  * ( ( f51 * ( val(k,i+3,j)+val(k,i-2,j) ) &
2419  + f52 * ( val(k,i+2,j)+val(k,i-1,j) ) &
2420  + f53 * ( val(k,i+1,j)+val(k,i,j) ) ) &
2421  - ( f51 * ( val(k,i+3,j)-val(k,i-2,j) ) &
2422  + f54 * ( val(k,i+2,j)-val(k,i-1,j) ) &
2423  + f55 * ( val(k,i+1,j)-val(k,i,j) ) ) * sign(1.0_rp,vel) ) &
2424  + gsqrt(k,i,j) * num_diff(k,i,j)
2425  enddo
2426  enddo
2427  enddo
2428  !$acc end kernels
2429 #ifdef DEBUG
2430  k = iundef; i = iundef; j = iundef
2431 #endif
2432 
2433 
2434 
2435  return
2436  end subroutine atmos_dyn_fvm_fluxx_xvz_ud5
2437 
2438  !-----------------------------------------------------------------------------
2440  subroutine atmos_dyn_fvm_fluxy_xvz_ud5( &
2441  flux, &
2442  mom, val, DENS, &
2443  GSQRT, MAPF, &
2444  num_diff, &
2445  CDZ, TwoD, &
2446  IIS, IIE, JJS, JJE )
2447  implicit none
2448 
2449  real(rp), intent(inout) :: flux (ka,ia,ja)
2450  real(rp), intent(in) :: mom (ka,ia,ja)
2451  real(rp), intent(in) :: val (ka,ia,ja)
2452  real(rp), intent(in) :: dens (ka,ia,ja)
2453  real(rp), intent(in) :: gsqrt (ka,ia,ja)
2454  real(rp), intent(in) :: mapf ( ia,ja,2)
2455  real(rp), intent(in) :: num_diff(ka,ia,ja)
2456  real(rp), intent(in) :: cdz (ka)
2457  logical, intent(in) :: twod
2458  integer, intent(in) :: iis, iie, jjs, jje
2459 
2460  real(rp) :: vel
2461  integer :: k, i, j
2462  !---------------------------------------------------------------------------
2463 
2464  ! note that y-index is added by -1
2465 
2466 
2467 
2468  !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
2469  !$omp private(vel) &
2470  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,GSQRT,MAPF,num_diff)
2471  !$acc kernels
2472  do j = jjs, jje+1
2473  do i = iis, iie
2474  do k = ks, ke
2475 #ifdef DEBUG
2476  call check( __line__, mom(k,i ,j) )
2477  call check( __line__, mom(k,i,j-1) )
2478 
2479  call check( __line__, val(k,i,j-1) )
2480  call check( __line__, val(k,i,j) )
2481 
2482  call check( __line__, val(k,i,j-2) )
2483  call check( __line__, val(k,i,j+1) )
2484 
2485  call check( __line__, val(k,i,j-3) )
2486  call check( __line__, val(k,i,j+2) )
2487 
2488 #endif
2489  vel = ( 0.5_rp * ( mom(k,i,j)+mom(k,i,j-1) ) ) &
2490  / ( dens(k,i,j) )
2491  flux(k,i,j-1) = gsqrt(k,i,j) / mapf(i,j,+1) * vel &
2492  * ( ( f51 * ( val(k,i,j+2)+val(k,i,j-3) ) &
2493  + f52 * ( val(k,i,j+1)+val(k,i,j-2) ) &
2494  + f53 * ( val(k,i,j)+val(k,i,j-1) ) ) &
2495  - ( f51 * ( val(k,i,j+2)-val(k,i,j-3) ) &
2496  + f54 * ( val(k,i,j+1)-val(k,i,j-2) ) &
2497  + f55 * ( val(k,i,j)-val(k,i,j-1) ) ) * sign(1.0_rp,vel) ) &
2498  + gsqrt(k,i,j) * num_diff(k,i,j)
2499  enddo
2500  enddo
2501  enddo
2502  !$acc end kernels
2503 #ifdef DEBUG
2504  k = iundef; i = iundef; j = iundef
2505 #endif
2506 
2507 
2508 
2509  return
2510  end subroutine atmos_dyn_fvm_fluxy_xvz_ud5
2511 
2512 
2513 
2514 
2515 
2516 
2517 
2518 end module scale_atmos_dyn_fvm_flux_ud5
2519 
2520 !--
2521 ! vi:set readonly sw=4 ts=8
2522 !
2523 !Local Variables:
2524 !mode: f90
2525 !buffer-read-only: t
2526 !End:
2527 !
2528 !++
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_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, TwoD, IIS, IIE, JJS, JJE)
calculation X-flux at XYW
Definition: scale_atmos_dyn_fvm_flux_ud5.F90:814
scale_index
module Index
Definition: scale_index.F90:11
scale_const::const_undef2
integer, parameter, public const_undef2
undefined value (INT2)
Definition: scale_const.F90:40
scale_precision
module PRECISION
Definition: scale_precision.F90:14
scale_atmos_grid_cartesc_index::ka
integer, public ka
Definition: scale_atmos_grid_cartesC_index.F90:47
scale_const::const_eps
real(rp), public const_eps
small number
Definition: scale_const.F90:35
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:109
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:471
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, TwoD, IIS, IIE, JJS, JJE)
calculation X-flux at UY
Definition: scale_atmos_dyn_fvm_flux_ud5.F90:1714
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:406
scale_atmos_grid_cartesc_index::i_xyz
integer, public i_xyz
Definition: scale_atmos_grid_cartesC_index.F90:91
scale_prc
module PROCESS
Definition: scale_prc.F90:11
scale_atmos_grid_cartesc_index::i_uyz
integer, public i_uyz
Definition: scale_atmos_grid_cartesC_index.F90:95
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:344
scale_tracer::k
real(rp), public k
Definition: scale_tracer.F90:45
scale_atmos_grid_cartesc_index
module atmosphere / grid / cartesC index
Definition: scale_atmos_grid_cartesC_index.F90:12
scale_const
module CONSTANT
Definition: scale_const.F90:11
scale_atmos_grid_cartesc_index::ia
integer, public ia
Definition: scale_atmos_grid_cartesC_index.F90:48
scale_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, TwoD, IIS, IIE, JJS, JJE)
calculation J13-flux at UYZ
Definition: scale_atmos_dyn_fvm_flux_ud5.F90:1298
scale_debug::check
subroutine, public check(current_line, v)
Undefined value checker.
Definition: scale_debug.F90:59
scale_atmos_grid_cartesc_index::i_xvz
integer, public i_xvz
Definition: scale_atmos_grid_cartesC_index.F90:96
scale_prof
module profiler
Definition: scale_prof.F90:11
scale_atmos_grid_cartesc_index::ja
integer, public ja
Definition: scale_atmos_grid_cartesC_index.F90:49
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, TwoD, IIS, IIE, JJS, JJE)
calculation z-flux at UY
Definition: scale_atmos_dyn_fvm_flux_ud5.F90:1007
scale_tracer
module TRACER
Definition: scale_tracer.F90:12
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, TwoD, IIS, IIE, JJS, JJE)
calculation J13-flux at XYW
Definition: scale_atmos_dyn_fvm_flux_ud5.F90:622
scale_atmos_grid_cartesc_index::ks
integer, public ks
start point of inner domain: z, local
Definition: scale_atmos_grid_cartesC_index.F90:51
scale_debug
module DEBUG
Definition: scale_debug.F90:11
scale_atmos_dyn_fvm_flux_ud5::atmos_dyn_fvm_fluxj23_xvz_ud5
subroutine, public atmos_dyn_fvm_fluxj23_xvz_ud5(flux, mom, val, DENS, GSQRT, J23G, MAPF, CDZ, TwoD, IIS, IIE, JJS, JJE)
calculation J23-flux at XVZ
Definition: scale_atmos_dyn_fvm_flux_ud5.F90:2223
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, TwoD, IIS, IIE, JJS, JJE)
calculation Y-flux at XV
Definition: scale_atmos_dyn_fvm_flux_ud5.F90:2447
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, TwoD, IIS, IIE, JJS, JJE)
calculation z-flux at XV
Definition: scale_atmos_dyn_fvm_flux_ud5.F90:1906
scale_const::const_undef
real(rp), public const_undef
Definition: scale_const.F90:43
scale_atmos_dyn_fvm_flux_ud5
module scale_atmos_dyn_fvm_flux_ud5
Definition: scale_atmos_dyn_fvm_flux_ud5.F90:16
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, TwoD, IIS, IIE, JJS, JJE)
calculation J23-flux at UYZ
Definition: scale_atmos_dyn_fvm_flux_ud5.F90:1449
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:202
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, TwoD, IIS, IIE, JJS, JJE)
calculation J23-flux at XYW
Definition: scale_atmos_dyn_fvm_flux_ud5.F90:717
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, TwoD, IIS, IIE, JJS, JJE)
calculation Y-flux at UY
Definition: scale_atmos_dyn_fvm_flux_ud5.F90:1788
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, TwoD, IIS, IIE, JJS, JJE)
calculation J13-flux at XVZ
Definition: scale_atmos_dyn_fvm_flux_ud5.F90:2072
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, TwoD, IIS, IIE, JJS, JJE)
calculation Y-flux at XYW
Definition: scale_atmos_dyn_fvm_flux_ud5.F90:910
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, TwoD, IIS, IIE, JJS, JJE)
calculation X-flux at XV
Definition: scale_atmos_dyn_fvm_flux_ud5.F90:2375
Generated by   doxygen 1.8.17