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