SCALE-RM
scale_atmos_dyn_fvm_flux_ud3.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_ud3
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_ud3
43 
44  public :: atmos_dyn_fvm_fluxz_xyz_ud3
45  public :: atmos_dyn_fvm_fluxx_xyz_ud3
46  public :: atmos_dyn_fvm_fluxy_xyz_ud3
47 
48  public :: atmos_dyn_fvm_fluxz_xyw_ud3
49  public :: atmos_dyn_fvm_fluxj13_xyw_ud3
50  public :: atmos_dyn_fvm_fluxj23_xyw_ud3
51  public :: atmos_dyn_fvm_fluxx_xyw_ud3
52  public :: atmos_dyn_fvm_fluxy_xyw_ud3
53 
54  public :: atmos_dyn_fvm_fluxz_uyz_ud3
55  public :: atmos_dyn_fvm_fluxj13_uyz_ud3
56  public :: atmos_dyn_fvm_fluxj23_uyz_ud3
57  public :: atmos_dyn_fvm_fluxx_uyz_ud3
58  public :: atmos_dyn_fvm_fluxy_uyz_ud3
59 
60  public :: atmos_dyn_fvm_fluxz_xvz_ud3
61  public :: atmos_dyn_fvm_fluxj13_xvz_ud3
62  public :: atmos_dyn_fvm_fluxj23_xvz_ud3
63  public :: atmos_dyn_fvm_fluxx_xvz_ud3
64  public :: atmos_dyn_fvm_fluxy_xvz_ud3
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 
93 
94 contains
95  !-----------------------------------------------------------------------------
97 !OCL SERIAL
98  subroutine atmos_dyn_fvm_flux_valuew_z_ud3( &
99  valW, &
100  mflx, val, GSQRT, &
101  CDZ )
102  implicit none
103 
104  real(RP), intent(out) :: valW (ka)
105  real(RP), intent(in) :: mflx (ka)
106  real(RP), intent(in) :: val (ka)
107  real(RP), intent(in) :: GSQRT(ka)
108  real(RP), intent(in) :: CDZ (ka)
109 
110  integer :: k
111  !---------------------------------------------------------------------------
112 
113  do k = ks+1, ke-2
114 #ifdef DEBUG
115  call check( __line__, mflx(k) )
116 
117  call check( __line__, val(k) )
118  call check( __line__, val(k+1) )
119 
120  call check( __line__, val(k-1) )
121  call check( __line__, val(k+2) )
122 
123 #endif
124  valw(k) = ( f31 * ( val(k+2)+val(k-1) ) + f32 * ( val(k+1)+val(k) ) ) &
125  - ( f31 * ( val(k+2)-val(k-1) ) + f33 * ( val(k+1)-val(k) ) ) * sign(1.0_rp,mflx(k))
126  enddo
127 #ifdef DEBUG
128  k = iundef
129 #endif
130 
131 #ifdef DEBUG
132 
133  call check( __line__, mflx(ks) )
134  call check( __line__, val(ks ) )
135  call check( __line__, val(ks+1) )
136  call check( __line__, mflx(ke-1) )
137  call check( __line__, val(ke ) )
138  call check( __line__, val(ke-1) )
139 
140 #endif
141 
142  valw(ks) = f2 * ( val(ks+1)+val(ks) ) &
143  * ( 0.5_rp + sign(0.5_rp,mflx(ks)) ) &
144  + ( 2.0_rp * val(ks) + 5.0_rp * val(ks+1) - val(ks+2) ) / 6.0_rp &
145  * ( 0.5_rp - sign(0.5_rp,mflx(ks)) )
146  valw(ke-1) = ( 2.0_rp * val(ke) + 5.0_rp * val(ke-1) - val(ke-2) ) / 6.0_rp &
147  * ( 0.5_rp + sign(0.5_rp,mflx(ke-1)) ) &
148  + f2 * ( val(ke)+val(ke-1) ) &
149  * ( 0.5_rp - sign(0.5_rp,mflx(ke-1)) )
150 
151 
152  return
153  end subroutine atmos_dyn_fvm_flux_valuew_z_ud3
154 
155  !-----------------------------------------------------------------------------
157  subroutine atmos_dyn_fvm_fluxz_xyz_ud3( &
158  flux, &
159  mflx, val, GSQRT, &
160  num_diff, &
161  CDZ, &
162  IIS, IIE, JJS, JJE )
163  use scale_const, only: &
164  eps => const_eps
165  implicit none
166 
167  real(RP), intent(inout) :: flux (ka,ia,ja)
168  real(RP), intent(in) :: mflx (ka,ia,ja)
169  real(RP), intent(in) :: val (ka,ia,ja)
170  real(RP), intent(in) :: GSQRT (ka,ia,ja)
171  real(RP), intent(in) :: num_diff(ka,ia,ja)
172  real(RP), intent(in) :: CDZ (ka)
173  integer, intent(in) :: IIS, IIE, JJS, JJE
174 
175  real(RP) :: vel
176  integer :: k, i, j
177  !---------------------------------------------------------------------------
178 
179  !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
180  !$omp private(vel) &
181  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mflx,val,flux,GSQRT,num_diff)
182  do j = jjs, jje
183  do i = iis, iie
184  do k = ks+1, ke-2
185 #ifdef DEBUG
186  call check( __line__, mflx(k,i,j) )
187 
188  call check( __line__, val(k,i,j) )
189  call check( __line__, val(k+1,i,j) )
190 
191  call check( __line__, val(k-1,i,j) )
192  call check( __line__, val(k+2,i,j) )
193 
194 #endif
195  vel = mflx(k,i,j)
196  flux(k,i,j) = vel &
197  * ( ( f31 * ( val(k+2,i,j)+val(k-1,i,j) ) + f32 * ( val(k+1,i,j)+val(k,i,j) ) ) &
198  - ( f31 * ( val(k+2,i,j)-val(k-1,i,j) ) + f33 * ( val(k+1,i,j)-val(k,i,j) ) ) * sign(1.0_rp,vel) ) &
199  + gsqrt(k,i,j) * num_diff(k,i,j)
200  enddo
201  enddo
202  enddo
203 #ifdef DEBUG
204  k = iundef; i = iundef; j = iundef
205 #endif
206 
207  !$omp parallel do default(none) private(i,j) OMP_SCHEDULE_ collapse(2) &
208  !$omp private(vel) &
209  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mflx,val,flux,GSQRT,num_diff,EPS)
210  do j = jjs, jje
211  do i = iis, iie
212 #ifdef DEBUG
213 
214  call check( __line__, mflx(ks,i,j) )
215  call check( __line__, val(ks ,i,j) )
216  call check( __line__, val(ks+1,i,j) )
217  call check( __line__, mflx(ke-1,i,j) )
218  call check( __line__, val(ke ,i,j) )
219  call check( __line__, val(ke-1,i,j) )
220 
221 #endif
222  flux(ks-1,i,j) = 0.0_rp
223 
224  vel = mflx(ks,i,j)
225  flux(ks,i,j) = vel &
226  * ( f2 * ( val(ks+1,i,j)+val(ks,i,j) ) &
227  * ( 0.5_rp + sign(0.5_rp,vel) ) &
228  + ( 2.0_rp * val(ks,i,j) + 5.0_rp * val(ks+1,i,j) - val(ks+2,i,j) ) / 6.0_rp &
229  * ( 0.5_rp - sign(0.5_rp,vel) ) ) &
230  + gsqrt(ks,i,j) * num_diff(ks,i,j)
231  vel = mflx(ke-1,i,j)
232  flux(ke-1,i,j) = vel &
233  * ( ( 2.0_rp * val(ke,i,j) + 5.0_rp * val(ke-1,i,j) - val(ke-2,i,j) ) / 6.0_rp &
234  * ( 0.5_rp + sign(0.5_rp,vel) ) &
235  + f2 * ( val(ke,i,j)+val(ke-1,i,j) ) &
236  * ( 0.5_rp - sign(0.5_rp,vel) ) ) &
237  + gsqrt(ke-1,i,j) * num_diff(ke-1,i,j)
238 
239  flux(ke ,i,j) = 0.0_rp
240  enddo
241  enddo
242 #ifdef DEBUG
243  k = iundef; i = iundef; j = iundef
244 #endif
245 
246  return
247  end subroutine atmos_dyn_fvm_fluxz_xyz_ud3
248 
249  !-----------------------------------------------------------------------------
251  subroutine atmos_dyn_fvm_fluxx_xyz_ud3( &
252  flux, &
253  mflx, val, GSQRT, &
254  num_diff, &
255  CDZ, &
256  IIS, IIE, JJS, JJE )
257  implicit none
258 
259  real(RP), intent(inout) :: flux (ka,ia,ja)
260  real(RP), intent(in) :: mflx (ka,ia,ja)
261  real(RP), intent(in) :: val (ka,ia,ja)
262  real(RP), intent(in) :: GSQRT (ka,ia,ja)
263  real(RP), intent(in) :: num_diff(ka,ia,ja)
264  real(RP), intent(in) :: CDZ(ka)
265  integer, intent(in) :: IIS, IIE, JJS, JJE
266 
267  real(RP) :: vel
268  integer :: k, i, j
269  !---------------------------------------------------------------------------
270 
271  !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
272  !$omp private(vel) &
273  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mflx,val,flux,GSQRT,num_diff)
274  do j = jjs, jje
275  do i = iis-1, iie
276  do k = ks, ke
277 #ifdef DEBUG
278  call check( __line__, mflx(k,i,j) )
279 
280  call check( __line__, val(k,i,j) )
281  call check( __line__, val(k,i+1,j) )
282 
283  call check( __line__, val(k,i-1,j) )
284  call check( __line__, val(k,i+2,j) )
285 
286 #endif
287  vel = mflx(k,i,j)
288  flux(k,i,j) = vel &
289  * ( ( f31 * ( val(k,i+2,j)+val(k,i-1,j) ) + f32 * ( val(k,i+1,j)+val(k,i,j) ) ) &
290  - ( f31 * ( val(k,i+2,j)-val(k,i-1,j) ) + f33 * ( val(k,i+1,j)-val(k,i,j) ) ) * sign(1.0_rp,vel) ) &
291  + gsqrt(k,i,j) * num_diff(k,i,j)
292  enddo
293  enddo
294  enddo
295 #ifdef DEBUG
296  k = iundef; i = iundef; j = iundef
297 #endif
298 
299  return
300  end subroutine atmos_dyn_fvm_fluxx_xyz_ud3
301 
302  !-----------------------------------------------------------------------------
304  subroutine atmos_dyn_fvm_fluxy_xyz_ud3( &
305  flux, &
306  mflx, val, GSQRT, &
307  num_diff, &
308  CDZ, &
309  IIS, IIE, JJS, JJE )
310  implicit none
311 
312  real(RP), intent(inout) :: flux (ka,ia,ja)
313  real(RP), intent(in) :: mflx (ka,ia,ja)
314  real(RP), intent(in) :: val (ka,ia,ja)
315  real(RP), intent(in) :: GSQRT (ka,ia,ja)
316  real(RP), intent(in) :: num_diff(ka,ia,ja)
317  real(RP), intent(in) :: CDZ(ka)
318  integer, intent(in) :: IIS, IIE, JJS, JJE
319 
320  real(RP) :: vel
321  integer :: k, i, j
322  !---------------------------------------------------------------------------
323 
324  !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
325  !$omp private(vel) &
326  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mflx,val,flux,GSQRT,num_diff)
327  do j = jjs-1, jje
328  do i = iis, iie
329  do k = ks, ke
330 #ifdef DEBUG
331  call check( __line__, mflx(k,i,j) )
332 
333  call check( __line__, val(k,i,j) )
334  call check( __line__, val(k,i,j+1) )
335 
336  call check( __line__, val(k,i,j-1) )
337  call check( __line__, val(k,i,j+2) )
338 
339 #endif
340  vel = mflx(k,i,j)
341  flux(k,i,j) = vel &
342  * ( ( f31 * ( val(k,i,j+2)+val(k,i,j-1) ) + f32 * ( val(k,i,j+1)+val(k,i,j) ) ) &
343  - ( f31 * ( val(k,i,j+2)-val(k,i,j-1) ) + f33 * ( val(k,i,j+1)-val(k,i,j) ) ) * sign(1.0_rp,vel) ) &
344  + gsqrt(k,i,j) * num_diff(k,i,j)
345  enddo
346  enddo
347  enddo
348 #ifdef DEBUG
349  k = iundef; i = iundef; j = iundef
350 #endif
351 
352  return
353  end subroutine atmos_dyn_fvm_fluxy_xyz_ud3
354 
355 
356  !-----------------------------------------------------------------------------
358  subroutine atmos_dyn_fvm_fluxz_xyw_ud3( &
359  flux, &
360  mom, val, DENS, &
361  GSQRT, J33G, &
362  num_diff, &
363  CDZ, FDZ, &
364  dtrk, &
365  IIS, IIE, JJS, JJE )
366  implicit none
367 
368  real(RP), intent(inout) :: flux (ka,ia,ja)
369  real(RP), intent(in) :: mom (ka,ia,ja)
370  real(RP), intent(in) :: val (ka,ia,ja)
371  real(RP), intent(in) :: DENS (ka,ia,ja)
372  real(RP), intent(in) :: GSQRT (ka,ia,ja)
373  real(RP), intent(in) :: J33G
374  real(RP), intent(in) :: num_diff(ka,ia,ja)
375  real(RP), intent(in) :: CDZ (ka)
376  real(RP), intent(in) :: FDZ (ka-1)
377  real(RP), intent(in) :: dtrk
378  integer, intent(in) :: IIS, IIE, JJS, JJE
379 
380  real(RP) :: vel
381  integer :: k, i, j
382  !---------------------------------------------------------------------------
383 
384  ! note than z-index is added by -1
385 
386  !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
387  !$omp private(vel) &
388  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,flux,J33G,GSQRT,num_diff,DENS)
389  do j = jjs, jje
390  do i = iis, iie
391  do k = ks+2, ke-1
392 #ifdef DEBUG
393  call check( __line__, mom(k-1,i,j) )
394  call check( __line__, mom(k ,i,j) )
395 
396  call check( __line__, val(k-1,i,j) )
397  call check( __line__, val(k,i,j) )
398 
399  call check( __line__, val(k-2,i,j) )
400  call check( __line__, val(k+1,i,j) )
401 
402 #endif
403  vel = ( 0.5_rp * ( mom(k-1,i,j) &
404  + mom(k,i,j) ) ) &
405  / dens(k,i,j)
406  flux(k-1,i,j) = j33g * vel &
407  * ( ( f31 * ( val(k+1,i,j)+val(k-2,i,j) ) + f32 * ( val(k,i,j)+val(k-1,i,j) ) ) &
408  - ( f31 * ( val(k+1,i,j)-val(k-2,i,j) ) + f33 * ( val(k,i,j)-val(k-1,i,j) ) ) * sign(1.0_rp,vel) ) &
409  + gsqrt(k,i,j) * num_diff(k,i,j)
410  enddo
411  enddo
412  enddo
413 #ifdef DEBUG
414  k = iundef; i = iundef; j = iundef
415 #endif
416 
417  !$omp parallel do default(none) private(i,j) OMP_SCHEDULE_ collapse(2) &
418  !$omp private(vel) &
419  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,J33G,GSQRT,num_diff,FDZ,dtrk)
420  do j = jjs, jje
421  do i = iis, iie
422 #ifdef DEBUG
423 
424  call check( __line__, val(ks ,i,j) )
425  call check( __line__, val(ks+1,i,j) )
426 
427 
428 #endif
429  ! The boundary condition is qflx_hi + qflxJ13 + qfluxJ23 = 0 at KS.
430  ! The flux at KS can be non-zero.
431  ! To reduce calculations, all the fluxes are set to zero.
432  flux(ks-1,i,j) = 0.0_rp ! k = KS
433 
434  vel = ( 0.5_rp * ( mom(ks,i,j) &
435  + mom(ks+1,i,j) ) ) &
436  / dens(ks+1,i,j)
437  flux(ks,i,j) = j33g * vel &
438  * ( f2 * ( val(ks+1,i,j)+val(ks,i,j) ) &
439  * ( 0.5_rp + sign(0.5_rp,vel) ) &
440  + ( 2.0_rp * val(ks,i,j) + 5.0_rp * val(ks+1,i,j) - val(ks+2,i,j) ) / 6.0_rp &
441  * ( 0.5_rp - sign(0.5_rp,vel) ) ) &
442  + gsqrt(ks+1,i,j) * num_diff(ks+1,i,j) ! k = KS+1
443 
444 
445 
446  flux(ke-1,i,j) = 0.0_rp ! k = KE
447  flux(ke ,i,j) = 0.0_rp ! k = KE+1
448  enddo
449  enddo
450 
451  return
452  end subroutine atmos_dyn_fvm_fluxz_xyw_ud3
453 
454 
455  !-----------------------------------------------------------------------------
457  subroutine atmos_dyn_fvm_fluxj13_xyw_ud3( &
458  flux, &
459  mom, val, DENS, &
460  GSQRT, J13G, MAPF, &
461  CDZ, &
462  IIS, IIE, JJS, JJE )
463  implicit none
464 
465  real(RP), intent(inout) :: flux (ka,ia,ja)
466  real(RP), intent(in) :: mom (ka,ia,ja)
467  real(RP), intent(in) :: val (ka,ia,ja)
468  real(RP), intent(in) :: DENS (ka,ia,ja)
469  real(RP), intent(in) :: GSQRT (ka,ia,ja)
470  real(RP), intent(in) :: J13G (ka,ia,ja)
471  real(RP), intent(in) :: MAPF ( ia,ja,2)
472  real(RP), intent(in) :: CDZ (ka)
473  integer, intent(in) :: IIS, IIE, JJS, JJE
474 
475  real(RP) :: vel
476  integer :: k, i, j
477  !---------------------------------------------------------------------------
478 
479  !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
480  !$omp private(vel) &
481  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,J13G,MAPF)
482  do j = jjs, jje
483  do i = iis, iie
484  do k = ks+2, ke-1
485  vel = ( 0.5_rp * ( mom(k,i,j)+mom(k,i-1,j) ) ) &
486  / dens(k,i,j)
487  vel = vel * j13g(k,i,j)
488  flux(k-1,i,j) = vel / mapf(i,j,+2) &
489  * ( ( f31 * ( val(k+1,i,j)+val(k-2,i,j) ) + f32 * ( val(k,i,j)+val(k-1,i,j) ) ) &
490  - ( f31 * ( val(k+1,i,j)-val(k-2,i,j) ) + f33 * ( val(k,i,j)-val(k-1,i,j) ) ) * sign(1.0_rp,vel) )
491  enddo
492  enddo
493  enddo
494 
495  !$omp parallel do default(none) private(i,j) OMP_SCHEDULE_ collapse(2) &
496  !$omp private(vel) &
497  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,J13G,MAPF)
498  do j = jjs, jje
499  do i = iis, iie
500  ! The boundary condition is qflx_hi + qflxJ13 + qfluxJ23 = 0 at KS.
501  ! The flux at KS can be non-zero.
502  ! To reduce calculations, all the fluxes are set to zero.
503  flux(ks-1,i,j) = 0.0_rp ! k = KS
504 
505  ! physically incorrect but for numerical stability
506  vel = ( ( 0.5_rp * ( mom(ks+1,i,j)+mom(ks+1,i-1,j) ) ) / dens(ks+1,i,j) &
507  + ( 0.5_rp * ( mom(ks,i,j)+mom(ks,i-1,j) ) ) / dens(ks ,i,j) ) * 0.5_rp
508 ! vel = ( 0.5_RP * ( mom(KS+1,i,j)+mom(KS+1,i-1,j) ) ) &
509 ! / DENS(KS+1,i,j)
510  vel = vel * j13g(ks+1,i,j)
511  flux(ks,i,j) = vel / mapf(i,j,+2) &
512  * ( f2 * ( val(ks+1,i,j)+val(ks,i,j) ) &
513  * ( 0.5_rp + sign(0.5_rp,vel) ) &
514  + ( 2.0_rp * val(ks,i,j) + 5.0_rp * val(ks+1,i,j) - val(ks+2,i,j) ) / 6.0_rp &
515  * ( 0.5_rp - sign(0.5_rp,vel) ) ) ! k = KS+1
516 
517 
518  flux(ke-1,i,j) = 0.0_rp
519  enddo
520  enddo
521 
522  return
523  end subroutine atmos_dyn_fvm_fluxj13_xyw_ud3
524 
525  !-----------------------------------------------------------------------------
527  subroutine atmos_dyn_fvm_fluxj23_xyw_ud3( &
528  flux, &
529  mom, val, DENS, &
530  GSQRT, J23G, MAPF, &
531  CDZ, &
532  IIS, IIE, JJS, JJE )
533  implicit none
534 
535  real(RP), intent(inout) :: flux (ka,ia,ja)
536  real(RP), intent(in) :: mom (ka,ia,ja)
537  real(RP), intent(in) :: val (ka,ia,ja)
538  real(RP), intent(in) :: DENS (ka,ia,ja)
539  real(RP), intent(in) :: GSQRT (ka,ia,ja)
540  real(RP), intent(in) :: J23G (ka,ia,ja)
541  real(RP), intent(in) :: MAPF ( ia,ja,2)
542  real(RP), intent(in) :: CDZ (ka)
543  integer, intent(in) :: IIS, IIE, JJS, JJE
544 
545  real(RP) :: vel
546  integer :: k, i, j
547  !---------------------------------------------------------------------------
548 
549  !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
550  !$omp private(vel) &
551  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,J23G,MAPF)
552  do j = jjs, jje
553  do i = iis, iie
554  do k = ks+2, ke-1
555  vel = ( 0.5_rp * ( mom(k,i,j)+mom(k,i,j-1) ) ) &
556  / dens(k,i,j)
557  vel = vel * j23g(k,i,j)
558  flux(k-1,i,j) = vel / mapf(i,j,+1) &
559  * ( ( f31 * ( val(k+1,i,j)+val(k-2,i,j) ) + f32 * ( val(k,i,j)+val(k-1,i,j) ) ) &
560  - ( f31 * ( val(k+1,i,j)-val(k-2,i,j) ) + f33 * ( val(k,i,j)-val(k-1,i,j) ) ) * sign(1.0_rp,vel) )
561  enddo
562  enddo
563  enddo
564 
565  !$omp parallel do default(none) private(i,j) OMP_SCHEDULE_ collapse(2) &
566  !$omp private(vel) &
567  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,J23G,MAPF)
568  do j = jjs, jje
569  do i = iis, iie
570  ! The boundary condition is qflx_hi + qflxJ13 + qfluxJ23 = 0 at KS.
571  ! The flux at KS can be non-zero.
572  ! To reduce calculations, all the fluxes are set to zero.
573  flux(ks-1,i,j) = 0.0_rp ! k = KS
574 
575  ! physically incorrect but for numerical stability
576  vel = ( ( 0.5_rp * ( mom(ks+1,i,j)+mom(ks+1,i,j-1) ) ) / dens(ks+1,i,j) &
577  + ( 0.5_rp * ( mom(ks,i,j)+mom(ks,i,j-1) ) ) / dens(ks ,i,j) ) * 0.5_rp
578 ! vel = ( 0.5_RP * ( mom(KS+1,i,j)+mom(KS+1,i,j-1) ) ) &
579 ! / DENS(KS+1,i,j)
580  vel = vel * j23g(ks+1,i,j)
581  flux(ks,i,j) = vel / mapf(i,j,+1) &
582  * ( f2 * ( val(ks+1,i,j)+val(ks,i,j) ) &
583  * ( 0.5_rp + sign(0.5_rp,vel) ) &
584  + ( 2.0_rp * val(ks,i,j) + 5.0_rp * val(ks+1,i,j) - val(ks+2,i,j) ) / 6.0_rp &
585  * ( 0.5_rp - sign(0.5_rp,vel) ) ) ! k = KS+1
586 
587 
588  flux(ke-1,i,j) = 0.0_rp
589  enddo
590  enddo
591 
592  return
593  end subroutine atmos_dyn_fvm_fluxj23_xyw_ud3
594 
595 
596  !-----------------------------------------------------------------------------
598  subroutine atmos_dyn_fvm_fluxx_xyw_ud3( &
599  flux, &
600  mom, val, DENS, &
601  GSQRT, MAPF, &
602  num_diff, &
603  CDZ, &
604  IIS, IIE, JJS, JJE )
605  implicit none
606 
607  real(RP), intent(inout) :: flux (ka,ia,ja)
608  real(RP), intent(in) :: mom (ka,ia,ja)
609  real(RP), intent(in) :: val (ka,ia,ja)
610  real(RP), intent(in) :: DENS (ka,ia,ja)
611  real(RP), intent(in) :: GSQRT (ka,ia,ja)
612  real(RP), intent(in) :: MAPF ( ia,ja,2)
613  real(RP), intent(in) :: num_diff(ka,ia,ja)
614  real(RP), intent(in) :: CDZ (ka)
615  integer, intent(in) :: IIS, IIE, JJS, JJE
616 
617  real(RP) :: vel
618  integer :: k, i, j
619  !---------------------------------------------------------------------------
620 
621  !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
622  !$omp private(vel) &
623  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,GSQRT,MAPF,num_diff) &
624  !$omp shared(CDZ)
625  do j = jjs, jje
626  do i = iis-1, iie
627  do k = ks, ke-1
628 #ifdef DEBUG
629  call check( __line__, mom(k ,i,j) )
630  call check( __line__, mom(k+1,i,j) )
631 
632  call check( __line__, val(k,i,j) )
633  call check( __line__, val(k,i+1,j) )
634 
635  call check( __line__, val(k,i-1,j) )
636  call check( __line__, val(k,i+2,j) )
637 
638 #endif
639  vel = ( f2h(k,1,i_uyz) &
640  * mom(k+1,i,j) &
641  + f2h(k,2,i_uyz) &
642  * mom(k,i,j) ) &
643  / ( f2h(k,1,i_uyz) &
644  * 0.5_rp * ( dens(k+1,i,j)+dens(k+1,i+1,j) ) &
645  + f2h(k,2,i_uyz) &
646  * 0.5_rp * ( dens(k,i,j)+dens(k,i+1,j) ) )
647  flux(k,i,j) = gsqrt(k,i,j) / mapf(i,j,+2) * vel &
648  * ( ( f31 * ( val(k,i+2,j)+val(k,i-1,j) ) + f32 * ( val(k,i+1,j)+val(k,i,j) ) ) &
649  - ( f31 * ( val(k,i+2,j)-val(k,i-1,j) ) + f33 * ( val(k,i+1,j)-val(k,i,j) ) ) * sign(1.0_rp,vel) ) &
650  + gsqrt(k,i,j) * num_diff(k,i,j)
651  enddo
652  enddo
653  enddo
654 #ifdef DEBUG
655  k = iundef; i = iundef; j = iundef
656 #endif
657 
658  !$omp parallel do default(none) private(i,j) OMP_SCHEDULE_ collapse(2) &
659  !$omp private(vel) &
660  !$omp shared(JJS,JJE,IIS,IIE,KE,flux)
661  do j = jjs, jje
662  do i = iis-1, iie
663  flux(ke,i,j) = 0.0_rp
664  enddo
665  enddo
666 #ifdef DEBUG
667  k = iundef; i = iundef; j = iundef
668 #endif
669 
670  return
671  end subroutine atmos_dyn_fvm_fluxx_xyw_ud3
672 
673  !-----------------------------------------------------------------------------
675  subroutine atmos_dyn_fvm_fluxy_xyw_ud3( &
676  flux, &
677  mom, val, DENS, &
678  GSQRT, MAPF, &
679  num_diff, &
680  CDZ, &
681  IIS, IIE, JJS, JJE )
682  implicit none
683 
684  real(RP), intent(inout) :: flux (ka,ia,ja)
685  real(RP), intent(in) :: mom (ka,ia,ja)
686  real(RP), intent(in) :: val (ka,ia,ja)
687  real(RP), intent(in) :: DENS (ka,ia,ja)
688  real(RP), intent(in) :: GSQRT (ka,ia,ja)
689  real(RP), intent(in) :: MAPF ( ia,ja,2)
690  real(RP), intent(in) :: num_diff(ka,ia,ja)
691  real(RP), intent(in) :: CDZ (ka)
692  integer, intent(in) :: IIS, IIE, JJS, JJE
693 
694  real(RP) :: vel
695  integer :: k, i, j
696  !---------------------------------------------------------------------------
697 
698  !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
699  !$omp private(vel) &
700  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,GSQRT,MAPF,num_diff) &
701  !$omp shared(CDZ)
702  do j = jjs-1, jje
703  do i = iis, iie
704  do k = ks, ke-1
705 #ifdef DEBUG
706  call check( __line__, mom(k ,i,j) )
707  call check( __line__, mom(k+1,i,j) )
708 
709  call check( __line__, val(k,i,j) )
710  call check( __line__, val(k,i,j+1) )
711 
712  call check( __line__, val(k,i,j-1) )
713  call check( __line__, val(k,i,j+2) )
714 
715 #endif
716  vel = ( f2h(k,1,i_xvz) &
717  * mom(k+1,i,j) &
718  + f2h(k,2,i_xvz) &
719  * mom(k,i,j) ) &
720  / ( f2h(k,1,i_xvz) &
721  * 0.5_rp * ( dens(k+1,i,j)+dens(k+1,i,j+1) ) &
722  + f2h(k,2,i_xvz) &
723  * 0.5_rp * ( dens(k,i,j)+dens(k,i,j+1) ) )
724  flux(k,i,j) = gsqrt(k,i,j) / mapf(i,j,+1) * vel &
725  * ( ( f31 * ( val(k,i,j+2)+val(k,i,j-1) ) + f32 * ( val(k,i,j+1)+val(k,i,j) ) ) &
726  - ( f31 * ( val(k,i,j+2)-val(k,i,j-1) ) + f33 * ( val(k,i,j+1)-val(k,i,j) ) ) * sign(1.0_rp,vel) ) &
727  + gsqrt(k,i,j) * num_diff(k,i,j)
728  enddo
729  enddo
730  enddo
731 #ifdef DEBUG
732  k = iundef; i = iundef; j = iundef
733 #endif
734 
735  !$omp parallel do default(none) private(i,j) OMP_SCHEDULE_ collapse(2) &
736  !$omp private(vel) &
737  !$omp shared(JJS,JJE,IIS,IIE,KE,flux)
738  do j = jjs-1, jje
739  do i = iis, iie
740  flux(ke,i,j) = 0.0_rp
741  enddo
742  enddo
743 #ifdef DEBUG
744  k = iundef; i = iundef; j = iundef
745 #endif
746 
747  return
748  end subroutine atmos_dyn_fvm_fluxy_xyw_ud3
749 
750 
751  !-----------------------------------------------------------------------------
753  subroutine atmos_dyn_fvm_fluxz_uyz_ud3( &
754  flux, &
755  mom, val, DENS, &
756  GSQRT, J33G, &
757  num_diff, &
758  CDZ, &
759  IIS, IIE, JJS, JJE )
760  implicit none
761 
762  real(RP), intent(inout) :: flux (ka,ia,ja)
763  real(RP), intent(in) :: mom (ka,ia,ja)
764  real(RP), intent(in) :: val (ka,ia,ja)
765  real(RP), intent(in) :: DENS (ka,ia,ja)
766  real(RP), intent(in) :: GSQRT (ka,ia,ja)
767  real(RP), intent(in) :: J33G
768  real(RP), intent(in) :: num_diff(ka,ia,ja)
769  real(RP), intent(in) :: CDZ (ka)
770  integer, intent(in) :: IIS, IIE, JJS, JJE
771 
772  real(RP) :: vel
773  integer :: k, i, j
774  !---------------------------------------------------------------------------
775 
776  !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
777  !$omp private(vel) &
778  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,J33G,GSQRT,num_diff) &
779  !$omp shared(CDZ)
780  do j = jjs, jje
781  do i = iis, iie
782  do k = ks+1, ke-2
783 #ifdef DEBUG
784  call check( __line__, mom(k,i,j) )
785  call check( __line__, mom(k,i+1,j) )
786 
787  call check( __line__, val(k,i,j) )
788  call check( __line__, val(k+1,i,j) )
789 
790  call check( __line__, val(k-1,i,j) )
791  call check( __line__, val(k+2,i,j) )
792 
793 #endif
794  vel = ( 0.5_rp * ( mom(k,i,j)+mom(k,i+1,j) ) ) &
795  / ( f2h(k,1,i_uyz) &
796  * 0.5_rp * ( dens(k+1,i,j)+dens(k+1,i+1,j) ) &
797  + f2h(k,2,i_uyz) &
798  * 0.5_rp * ( dens(k,i,j)+dens(k,i+1,j) ) )
799  flux(k,i,j) = j33g * vel &
800  * ( ( f31 * ( val(k+2,i,j)+val(k-1,i,j) ) + f32 * ( val(k+1,i,j)+val(k,i,j) ) ) &
801  - ( f31 * ( val(k+2,i,j)-val(k-1,i,j) ) + f33 * ( val(k+1,i,j)-val(k,i,j) ) ) * sign(1.0_rp,vel) ) &
802  + gsqrt(k,i,j) * num_diff(k,i,j)
803  enddo
804  enddo
805  enddo
806 #ifdef DEBUG
807  k = iundef; i = iundef; j = iundef
808 #endif
809 
810  !$omp parallel do default(none) private(i,j) OMP_SCHEDULE_ collapse(2) &
811  !$omp private(vel) &
812  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,flux,J33G,GSQRT,num_diff,DENS,CDZ)
813  do j = jjs, jje
814  do i = iis, iie
815 #ifdef DEBUG
816 
817  call check( __line__, mom(ks,i ,j) )
818  call check( __line__, mom(ks,i+1,j) )
819  call check( __line__, val(ks+1,i,j) )
820  call check( __line__, val(ks,i,j) )
821 
822 #endif
823  ! The boundary condition is qflx_hi + qflxJ13 + qfluxJ23 = 0 at KS-1.
824  ! The flux at KS-1 can be non-zero.
825  ! To reduce calculations, all the fluxes are set to zero.
826  flux(ks-1,i,j) = 0.0_rp
827 
828  vel = ( 0.5_rp * ( mom(ks,i,j)+mom(ks,i+1,j) ) ) &
829  / ( f2h(ks,1,i_uyz) &
830  * 0.5_rp * ( dens(ks+1,i,j)+dens(ks+1,i+1,j) ) &
831  + f2h(ks,2,i_uyz) &
832  * 0.5_rp * ( dens(ks,i,j)+dens(ks,i+1,j) ) )
833  flux(ks,i,j) = j33g * vel &
834  * ( f2 * ( val(ks+1,i,j)+val(ks,i,j) ) &
835  * ( 0.5_rp + sign(0.5_rp,vel) ) &
836  + ( 2.0_rp * val(ks,i,j) + 5.0_rp * val(ks+1,i,j) - val(ks+2,i,j) ) / 6.0_rp &
837  * ( 0.5_rp - sign(0.5_rp,vel) ) ) &
838  + gsqrt(ks,i,j) * num_diff(ks,i,j)
839  vel = ( 0.5_rp * ( mom(ke-1,i,j)+mom(ke-1,i+1,j) ) ) &
840  / ( f2h(ke-1,1,i_uyz) &
841  * 0.5_rp * ( dens(ke,i,j)+dens(ke,i+1,j) ) &
842  + f2h(ke-1,2,i_uyz) &
843  * 0.5_rp * ( dens(ke-1,i,j)+dens(ke-1,i+1,j) ) )
844  flux(ke-1,i,j) = j33g * vel &
845  * ( ( 2.0_rp * val(ke,i,j) + 5.0_rp * val(ke-1,i,j) - val(ke-2,i,j) ) / 6.0_rp &
846  * ( 0.5_rp + sign(0.5_rp,vel) ) &
847  + f2 * ( val(ke,i,j)+val(ke-1,i,j) ) &
848  * ( 0.5_rp - sign(0.5_rp,vel) ) ) &
849  + gsqrt(ke-1,i,j) * num_diff(ke-1,i,j)
850 
851  flux(ke,i,j) = 0.0_rp
852  enddo
853  enddo
854 #ifdef DEBUG
855  k = iundef; i = iundef; j = iundef
856 #endif
857 
858  return
859  end subroutine atmos_dyn_fvm_fluxz_uyz_ud3
860 
861  !-----------------------------------------------------------------------------
863  subroutine atmos_dyn_fvm_fluxj13_uyz_ud3( &
864  flux, &
865  mom, val, DENS, &
866  GSQRT, J13G, MAPF, &
867  CDZ, &
868  IIS, IIE, JJS, JJE )
869  implicit none
870 
871  real(RP), intent(inout) :: flux (ka,ia,ja)
872  real(RP), intent(in) :: mom (ka,ia,ja)
873  real(RP), intent(in) :: val (ka,ia,ja)
874  real(RP), intent(in) :: DENS (ka,ia,ja)
875  real(RP), intent(in) :: GSQRT (ka,ia,ja)
876  real(RP), intent(in) :: J13G (ka,ia,ja)
877  real(RP), intent(in) :: MAPF ( ia,ja,2)
878  real(RP), intent(in) :: CDZ (ka)
879  integer, intent(in) :: IIS, IIE, JJS, JJE
880 
881  real(RP) :: vel
882  integer :: k, i, j
883  !---------------------------------------------------------------------------
884 
885  !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
886  !$omp private(vel) &
887  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,J13G,MAPF) &
888  !$omp shared(GSQRT,CDZ)
889  do j = jjs, jje
890  do i = iis, iie
891  do k = ks+1, ke-2
892  vel = ( f2h(k,1,i_uyz) &
893  * mom(k+1,i,j) &
894  + f2h(k,2,i_uyz) &
895  * mom(k,i,j) ) &
896  / ( f2h(k,1,i_uyz) &
897  * 0.5_rp * ( dens(k+1,i,j)+dens(k+1,i+1,j) ) &
898  + f2h(k,2,i_uyz) &
899  * 0.5_rp * ( dens(k,i,j)+dens(k,i+1,j) ) )
900  vel = vel * j13g(k,i,j)
901  flux(k,i,j) = vel / mapf(i,j,+2) &
902  * ( ( f31 * ( val(k+2,i,j)+val(k-1,i,j) ) + f32 * ( val(k+1,i,j)+val(k,i,j) ) ) &
903  - ( f31 * ( val(k+2,i,j)-val(k-1,i,j) ) + f33 * ( val(k+1,i,j)-val(k,i,j) ) ) * sign(1.0_rp,vel) )
904  enddo
905  enddo
906  enddo
907 
908  !$omp parallel do default(none) private(i,j) OMP_SCHEDULE_ collapse(2) &
909  !$omp private(vel) &
910  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,J13G,MAPF) &
911  !$omp shared(GSQRT,CDZ)
912  do j = jjs, jje
913  do i = iis, iie
914  ! The boundary condition is qflx_hi + qflxJ13 + qfluxJ23 = 0 at KS-1.
915  ! The flux at KS-1 can be non-zero.
916  ! To reduce calculations, all the fluxes are set to zero.
917  flux(ks-1,i,j) = 0.0_rp
918 
919  vel = ( f2h(ks,1,i_uyz) &
920  * mom(ks+1,i,j) &
921  + f2h(ks,2,i_uyz) &
922  * mom(ks,i,j) ) &
923  / ( f2h(ks,1,i_uyz) &
924  * 0.5_rp * ( dens(ks+1,i,j)+dens(ks+1,i+1,j) ) &
925  + f2h(ks,2,i_uyz) &
926  * 0.5_rp * ( dens(ks,i,j)+dens(ks,i+1,j) ) )
927  vel = vel * j13g(ks,i,j)
928  flux(ks,i,j) = vel / mapf(i,j,+2) &
929  * ( f2 * ( val(ks+1,i,j)+val(ks,i,j) ) &
930  * ( 0.5_rp + sign(0.5_rp,vel) ) &
931  + ( 2.0_rp * val(ks,i,j) + 5.0_rp * val(ks+1,i,j) - val(ks+2,i,j) ) / 6.0_rp &
932  * ( 0.5_rp - sign(0.5_rp,vel) ) )
933 
934  vel = ( f2h(ke-1,1,i_uyz) &
935  * mom(ke,i,j) &
936  + f2h(ke-1,2,i_uyz) &
937  * mom(ke-1,i,j) ) &
938  / ( f2h(ke-1,1,i_uyz) &
939  * 0.5_rp * ( dens(ke,i,j)+dens(ke,i+1,j) ) &
940  + f2h(ke-1,2,i_uyz) &
941  * 0.5_rp * ( dens(ke-1,i,j)+dens(ke-1,i+1,j) ) )
942  vel = vel * j13g(ke-1,i,j)
943  flux(ke-1,i,j) = vel / mapf(i,j,+2) &
944  * ( ( 2.0_rp * val(ke,i,j) + 5.0_rp * val(ke-1,i,j) - val(ke-2,i,j) ) / 6.0_rp &
945  * ( 0.5_rp + sign(0.5_rp,vel) ) &
946  + f2 * ( val(ke,i,j)+val(ke-1,i,j) ) &
947  * ( 0.5_rp - sign(0.5_rp,vel) ) )
948 
949  flux(ke ,i,j) = 0.0_rp
950  enddo
951  enddo
952 
953  return
954  end subroutine atmos_dyn_fvm_fluxj13_uyz_ud3
955 
956  !-----------------------------------------------------------------------------
958  subroutine atmos_dyn_fvm_fluxj23_uyz_ud3( &
959  flux, &
960  mom, val, DENS, &
961  GSQRT, J23G, MAPF, &
962  CDZ, &
963  IIS, IIE, JJS, JJE )
964  implicit none
965 
966  real(RP), intent(inout) :: flux (ka,ia,ja)
967  real(RP), intent(in) :: mom (ka,ia,ja)
968  real(RP), intent(in) :: val (ka,ia,ja)
969  real(RP), intent(in) :: DENS (ka,ia,ja)
970  real(RP), intent(in) :: GSQRT (ka,ia,ja)
971  real(RP), intent(in) :: J23G (ka,ia,ja)
972  real(RP), intent(in) :: MAPF ( ia,ja,2)
973  real(RP), intent(in) :: CDZ (ka)
974  integer, intent(in) :: IIS, IIE, JJS, JJE
975 
976  real(RP) :: vel
977  integer :: k, i, j
978  !---------------------------------------------------------------------------
979 
980  !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
981  !$omp private(vel) &
982  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,J23G,MAPF) &
983  !$omp shared(GSQRT,CDZ)
984  do j = jjs, jje
985  do i = iis, iie
986  do k = ks+1, ke-2
987  vel = ( f2h(k,1,i_uyz) &
988  * 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) ) &
989  + f2h(k,2,i_uyz) &
990  * 0.25_rp * ( mom(k,i,j)+mom(k,i+1,j)+mom(k,i,j-1)+mom(k,i+1,j-1) ) ) &
991  / ( f2h(k,1,i_uyz) &
992  * 0.5_rp * ( dens(k+1,i,j)+dens(k+1,i+1,j) ) &
993  + f2h(k,2,i_uyz) &
994  * 0.5_rp * ( dens(k,i,j)+dens(k,i+1,j) ) )
995  vel = vel * j23g(k,i,j)
996  flux(k,i,j) = vel / mapf(i,j,+1) &
997  * ( ( f31 * ( val(k+2,i,j)+val(k-1,i,j) ) + f32 * ( val(k+1,i,j)+val(k,i,j) ) ) &
998  - ( f31 * ( val(k+2,i,j)-val(k-1,i,j) ) + f33 * ( val(k+1,i,j)-val(k,i,j) ) ) * sign(1.0_rp,vel) )
999  enddo
1000  enddo
1001  enddo
1002 
1003  !$omp parallel do default(none) private(i,j) OMP_SCHEDULE_ collapse(2) &
1004  !$omp private(vel) &
1005  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,J23G,MAPF) &
1006  !$omp shared(GSQRT,CDZ)
1007  do j = jjs, jje
1008  do i = iis, iie
1009  ! The boundary condition is qflx_hi + qflxJ13 + qfluxJ23 = 0 at KS-1.
1010  ! The flux at KS-1 can be non-zero.
1011  ! To reduce calculations, all the fluxes are set to zero.
1012  flux(ks-1,i,j) = 0.0_rp
1013 
1014  vel = ( f2h(ks,1,i_uyz) &
1015  * 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) ) &
1016  + f2h(ks,2,i_uyz) &
1017  * 0.25_rp * ( mom(ks,i,j)+mom(ks,i+1,j)+mom(ks,i,j-1)+mom(ks,i+1,j-1) ) ) &
1018  / ( f2h(ks,1,i_uyz) &
1019  * 0.5_rp * ( dens(ks+1,i,j)+dens(ks+1,i+1,j) ) &
1020  + f2h(ks,2,i_uyz) &
1021  * 0.5_rp * ( dens(ks,i,j)+dens(ks,i+1,j) ) )
1022  vel = vel * j23g(ks,i,j)
1023  flux(ks,i,j) = vel / mapf(i,j,+1) &
1024  * ( f2 * ( val(ks+1,i,j)+val(ks,i,j) ) &
1025  * ( 0.5_rp + sign(0.5_rp,vel) ) &
1026  + ( 2.0_rp * val(ks,i,j) + 5.0_rp * val(ks+1,i,j) - val(ks+2,i,j) ) / 6.0_rp &
1027  * ( 0.5_rp - sign(0.5_rp,vel) ) )
1028 
1029  vel = ( f2h(ke-1,1,i_uyz) &
1030  * 0.25_rp * ( mom(ke,i,j)+mom(ke,i+1,j)+mom(ke,i,j-1)+mom(ke,i+1,j-1) ) &
1031  + f2h(ke-1,2,i_uyz) &
1032  * 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) ) ) &
1033  / ( f2h(ke-1,1,i_uyz) &
1034  * 0.5_rp * ( dens(ke,i,j)+dens(ke,i+1,j) ) &
1035  + f2h(ke-1,2,i_uyz) &
1036  * 0.5_rp * ( dens(ke-1,i,j)+dens(ke-1,i+1,j) ) )
1037  vel = vel * j23g(ke-1,i,j)
1038  flux(ke-1,i,j) = vel / mapf(i,j,+1) &
1039  * ( ( 2.0_rp * val(ke,i,j) + 5.0_rp * val(ke-1,i,j) - val(ke-2,i,j) ) / 6.0_rp &
1040  * ( 0.5_rp + sign(0.5_rp,vel) ) &
1041  + f2 * ( val(ke,i,j)+val(ke-1,i,j) ) &
1042  * ( 0.5_rp - sign(0.5_rp,vel) ) )
1043 
1044  flux(ke ,i,j) = 0.0_rp
1045  enddo
1046  enddo
1047 
1048  return
1049  end subroutine atmos_dyn_fvm_fluxj23_uyz_ud3
1050 
1051  !-----------------------------------------------------------------------------
1053  subroutine atmos_dyn_fvm_fluxx_uyz_ud3( &
1054  flux, &
1055  mom, val, DENS, &
1056  GSQRT, MAPF, &
1057  num_diff, &
1058  CDZ, &
1059  IIS, IIE, JJS, JJE )
1060  implicit none
1061 
1062  real(RP), intent(inout) :: flux (ka,ia,ja)
1063  real(RP), intent(in) :: mom (ka,ia,ja)
1064  real(RP), intent(in) :: val (ka,ia,ja)
1065  real(RP), intent(in) :: DENS (ka,ia,ja)
1066  real(RP), intent(in) :: GSQRT (ka,ia,ja)
1067  real(RP), intent(in) :: MAPF ( ia,ja,2)
1068  real(RP), intent(in) :: num_diff(ka,ia,ja)
1069  real(RP), intent(in) :: CDZ (ka)
1070  integer, intent(in) :: IIS, IIE, JJS, JJE
1071 
1072  real(RP) :: vel
1073  integer :: k, i, j
1074  !---------------------------------------------------------------------------
1075 
1076  ! note that x-index is added by -1
1077 
1078  !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
1079  !$omp private(vel) &
1080  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,GSQRT,MAPF,num_diff)
1081  do j = jjs, jje
1082  do i = iis, iie+1
1083  do k = ks, ke
1084 #ifdef DEBUG
1085  call check( __line__, mom(k,i ,j) )
1086  call check( __line__, mom(k,i-1,j) )
1087 
1088  call check( __line__, val(k,i-1,j) )
1089  call check( __line__, val(k,i,j) )
1090 
1091  call check( __line__, val(k,i-2,j) )
1092  call check( __line__, val(k,i+1,j) )
1093 
1094 #endif
1095  vel = ( 0.5_rp * ( mom(k,i,j)+mom(k,i-1,j) ) ) &
1096  / ( dens(k,i,j) )
1097  flux(k,i-1,j) = gsqrt(k,i,j) / mapf(i,j,+2) * vel &
1098  * ( ( f31 * ( val(k,i+1,j)+val(k,i-2,j) ) + f32 * ( val(k,i,j)+val(k,i-1,j) ) ) &
1099  - ( f31 * ( val(k,i+1,j)-val(k,i-2,j) ) + f33 * ( val(k,i,j)-val(k,i-1,j) ) ) * sign(1.0_rp,vel) ) &
1100  + gsqrt(k,i,j) * num_diff(k,i,j)
1101  enddo
1102  enddo
1103  enddo
1104 #ifdef DEBUG
1105  k = iundef; i = iundef; j = iundef
1106 #endif
1107 
1108  return
1109  end subroutine atmos_dyn_fvm_fluxx_uyz_ud3
1110 
1111  !-----------------------------------------------------------------------------
1113  subroutine atmos_dyn_fvm_fluxy_uyz_ud3( &
1114  flux, &
1115  mom, val, DENS, &
1116  GSQRT, MAPF, &
1117  num_diff, &
1118  CDZ, &
1119  IIS, IIE, JJS, JJE )
1120  implicit none
1121 
1122  real(RP), intent(inout) :: flux (ka,ia,ja)
1123  real(RP), intent(in) :: mom (ka,ia,ja)
1124  real(RP), intent(in) :: val (ka,ia,ja)
1125  real(RP), intent(in) :: DENS (ka,ia,ja)
1126  real(RP), intent(in) :: GSQRT (ka,ia,ja)
1127  real(RP), intent(in) :: MAPF ( ia,ja,2)
1128  real(RP), intent(in) :: num_diff(ka,ia,ja)
1129  real(RP), intent(in) :: CDZ (ka)
1130  integer, intent(in) :: IIS, IIE, JJS, JJE
1131 
1132  real(RP) :: vel
1133  integer :: k, i, j
1134  !---------------------------------------------------------------------------
1135 
1136  !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
1137  !$omp private(vel) &
1138  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,GSQRT,MAPF,num_diff)
1139  do j = jjs-1, jje
1140  do i = iis, iie
1141  do k = ks, ke
1142 #ifdef DEBUG
1143  call check( __line__, mom(k,i ,j) )
1144  call check( __line__, mom(k,i-1,j) )
1145 
1146  call check( __line__, val(k,i,j) )
1147  call check( __line__, val(k,i,j+1) )
1148 
1149  call check( __line__, val(k,i,j-1) )
1150  call check( __line__, val(k,i,j+2) )
1151 
1152 #endif
1153  vel = ( 0.5_rp * ( mom(k,i,j)+mom(k,i+1,j) ) ) &
1154  / ( 0.25_rp * ( dens(k,i,j)+dens(k,i+1,j)+dens(k,i,j+1)+dens(k,i+1,j+1) ) )
1155  flux(k,i,j) = gsqrt(k,i,j) / mapf(i,j,+1) * vel &
1156  * ( ( f31 * ( val(k,i,j+2)+val(k,i,j-1) ) + f32 * ( val(k,i,j+1)+val(k,i,j) ) ) &
1157  - ( f31 * ( val(k,i,j+2)-val(k,i,j-1) ) + f33 * ( val(k,i,j+1)-val(k,i,j) ) ) * sign(1.0_rp,vel) ) &
1158  + gsqrt(k,i,j) * num_diff(k,i,j)
1159  enddo
1160  enddo
1161  enddo
1162 #ifdef DEBUG
1163  k = iundef; i = iundef; j = iundef
1164 #endif
1165 
1166  return
1167  end subroutine atmos_dyn_fvm_fluxy_uyz_ud3
1168 
1169 
1170 
1171  !-----------------------------------------------------------------------------
1173  subroutine atmos_dyn_fvm_fluxz_xvz_ud3( &
1174  flux, &
1175  mom, val, DENS, &
1176  GSQRT, J33G, &
1177  num_diff, &
1178  CDZ, &
1179  IIS, IIE, JJS, JJE )
1180  implicit none
1181 
1182  real(RP), intent(inout) :: flux (ka,ia,ja)
1183  real(RP), intent(in) :: mom (ka,ia,ja)
1184  real(RP), intent(in) :: val (ka,ia,ja)
1185  real(RP), intent(in) :: DENS (ka,ia,ja)
1186  real(RP), intent(in) :: GSQRT (ka,ia,ja)
1187  real(RP), intent(in) :: J33G
1188  real(RP), intent(in) :: num_diff(ka,ia,ja)
1189  real(RP), intent(in) :: CDZ (ka)
1190  integer, intent(in) :: IIS, IIE, JJS, JJE
1191 
1192  real(RP) :: vel
1193  integer :: k, i, j
1194  !---------------------------------------------------------------------------
1195 
1196  !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
1197  !$omp private(vel) &
1198  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,J33G,GSQRT,num_diff) &
1199  !$omp shared(CDZ)
1200  do j = jjs, jje
1201  do i = iis, iie
1202  do k = ks+1, ke-2
1203 #ifdef DEBUG
1204  call check( __line__, mom(k,i,j) )
1205  call check( __line__, mom(k,i,j+1) )
1206 
1207  call check( __line__, val(k,i,j) )
1208  call check( __line__, val(k+1,i,j) )
1209 
1210  call check( __line__, val(k-1,i,j) )
1211  call check( __line__, val(k+2,i,j) )
1212 
1213 #endif
1214  vel = ( 0.5_rp * ( mom(k,i,j)+mom(k,i,j+1) ) ) &
1215  / ( f2h(k,1,i_xvz) &
1216  * 0.5_rp * ( dens(k+1,i,j)+dens(k+1,i,j+1) ) &
1217  + f2h(k,2,i_xvz) &
1218  * 0.5_rp * ( dens(k,i,j)+dens(k,i,j+1) ) )
1219  flux(k,i,j) = j33g * vel &
1220  * ( ( f31 * ( val(k+2,i,j)+val(k-1,i,j) ) + f32 * ( val(k+1,i,j)+val(k,i,j) ) ) &
1221  - ( f31 * ( val(k+2,i,j)-val(k-1,i,j) ) + f33 * ( val(k+1,i,j)-val(k,i,j) ) ) * sign(1.0_rp,vel) ) &
1222  + gsqrt(k,i,j) * num_diff(k,i,j)
1223  enddo
1224  enddo
1225  enddo
1226 #ifdef DEBUG
1227  k = iundef; i = iundef; j = iundef
1228 #endif
1229 
1230  !$omp parallel do default(none) private(i,j) OMP_SCHEDULE_ collapse(2) &
1231  !$omp private(vel) &
1232  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,flux,J33G,GSQRT,num_diff,DENS,CDZ)
1233  do j = jjs, jje
1234  do i = iis, iie
1235 #ifdef DEBUG
1236 
1237  call check( __line__, mom(ks,i ,j) )
1238  call check( __line__, mom(ks,i,j+1) )
1239  call check( __line__, val(ks+1,i,j) )
1240  call check( __line__, val(ks,i,j) )
1241 
1242 #endif
1243  ! The boundary condition is qflx_hi + qflxJ13 + qfluxJ23 = 0 at KS-1.
1244  ! The flux at KS-1 can be non-zero.
1245  ! To reduce calculations, all the fluxes are set to zero.
1246  flux(ks-1,i,j) = 0.0_rp
1247 
1248  vel = ( 0.5_rp * ( mom(ks,i,j)+mom(ks,i,j+1) ) ) &
1249  / ( f2h(ks,1,i_xvz) &
1250  * 0.5_rp * ( dens(ks+1,i,j)+dens(ks+1,i,j+1) ) &
1251  + f2h(ks,2,i_xvz) &
1252  * 0.5_rp * ( dens(ks,i,j)+dens(ks,i,j+1) ) )
1253  flux(ks,i,j) = j33g * vel &
1254  * ( f2 * ( val(ks+1,i,j)+val(ks,i,j) ) &
1255  * ( 0.5_rp + sign(0.5_rp,vel) ) &
1256  + ( 2.0_rp * val(ks,i,j) + 5.0_rp * val(ks+1,i,j) - val(ks+2,i,j) ) / 6.0_rp &
1257  * ( 0.5_rp - sign(0.5_rp,vel) ) ) &
1258  + gsqrt(ks,i,j) * num_diff(ks,i,j)
1259  vel = ( 0.5_rp * ( mom(ke-1,i,j)+mom(ke-1,i,j+1) ) ) &
1260  / ( f2h(ke-1,1,i_xvz) &
1261  * 0.5_rp * ( dens(ke,i,j)+dens(ke,i,j+1) ) &
1262  + f2h(ke-1,2,i_xvz) &
1263  * 0.5_rp * ( dens(ke-1,i,j)+dens(ke-1,i,j+1) ) )
1264  flux(ke-1,i,j) = j33g * vel &
1265  * ( ( 2.0_rp * val(ke,i,j) + 5.0_rp * val(ke-1,i,j) - val(ke-2,i,j) ) / 6.0_rp &
1266  * ( 0.5_rp + sign(0.5_rp,vel) ) &
1267  + f2 * ( val(ke,i,j)+val(ke-1,i,j) ) &
1268  * ( 0.5_rp - sign(0.5_rp,vel) ) ) &
1269  + gsqrt(ke-1,i,j) * num_diff(ke-1,i,j)
1270 
1271  flux(ke,i,j) = 0.0_rp
1272  enddo
1273  enddo
1274 #ifdef DEBUG
1275  k = iundef; i = iundef; j = iundef
1276 #endif
1277 
1278  return
1279  end subroutine atmos_dyn_fvm_fluxz_xvz_ud3
1280 
1281  !-----------------------------------------------------------------------------
1283  subroutine atmos_dyn_fvm_fluxj13_xvz_ud3( &
1284  flux, &
1285  mom, val, DENS, &
1286  GSQRT, J13G, MAPF, &
1287  CDZ, &
1288  IIS, IIE, JJS, JJE )
1289  implicit none
1290 
1291  real(RP), intent(inout) :: flux (ka,ia,ja)
1292  real(RP), intent(in) :: mom (ka,ia,ja)
1293  real(RP), intent(in) :: val (ka,ia,ja)
1294  real(RP), intent(in) :: DENS (ka,ia,ja)
1295  real(RP), intent(in) :: GSQRT (ka,ia,ja)
1296  real(RP), intent(in) :: J13G (ka,ia,ja)
1297  real(RP), intent(in) :: MAPF ( ia,ja,2)
1298  real(RP), intent(in) :: CDZ (ka)
1299  integer, intent(in) :: IIS, IIE, JJS, JJE
1300 
1301  real(RP) :: vel
1302  integer :: k, i, j
1303  !---------------------------------------------------------------------------
1304 
1305  !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
1306  !$omp private(vel) &
1307  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,J13G,MAPF) &
1308  !$omp shared(GSQRT,CDZ)
1309  do j = jjs, jje
1310  do i = iis, iie
1311  do k = ks+1, ke-2
1312  vel = ( f2h(k,1,i_xvz) &
1313  * 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) ) &
1314  + f2h(k,2,i_xvz) &
1315  * 0.25_rp * ( mom(k,i,j)+mom(k,i-1,j)+mom(k,i,j+1)+mom(k,i-1,j+1) ) ) &
1316  / ( f2h(k,1,i_xvz) &
1317  * 0.5_rp * ( dens(k+1,i,j)+dens(k+1,i,j+1) ) &
1318  + f2h(k,2,i_xvz) &
1319  * 0.5_rp * ( dens(k,i,j)+dens(k,i,j+1) ) )
1320  vel = vel * j13g(k,i,j)
1321  flux(k,i,j) = vel / mapf(i,j,+2) &
1322  * ( ( f31 * ( val(k+2,i,j)+val(k-1,i,j) ) + f32 * ( val(k+1,i,j)+val(k,i,j) ) ) &
1323  - ( f31 * ( val(k+2,i,j)-val(k-1,i,j) ) + f33 * ( val(k+1,i,j)-val(k,i,j) ) ) * sign(1.0_rp,vel) )
1324  enddo
1325  enddo
1326  enddo
1327 
1328  !$omp parallel do default(none) private(i,j) OMP_SCHEDULE_ collapse(2) &
1329  !$omp private(vel) &
1330  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,J13G,MAPF) &
1331  !$omp shared(GSQRT,CDZ)
1332  do j = jjs, jje
1333  do i = iis, iie
1334  ! The boundary condition is qflx_hi + qflxJ13 + qfluxJ23 = 0 at KS-1.
1335  ! The flux at KS-1 can be non-zero.
1336  ! To reduce calculations, all the fluxes are set to zero.
1337  flux(ks-1,i,j) = 0.0_rp
1338 
1339  vel = ( f2h(ks,1,i_xvz) &
1340  * 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) ) &
1341  + f2h(ks,2,i_xvz) &
1342  * 0.25_rp * ( mom(ks,i,j)+mom(ks,i-1,j)+mom(ks,i,j+1)+mom(ks,i-1,j+1) ) ) &
1343  / ( f2h(ks,1,i_xvz) &
1344  * 0.5_rp * ( dens(ks+1,i,j)+dens(ks+1,i,j+1) ) &
1345  + f2h(ks,2,i_xvz) &
1346  * 0.5_rp * ( dens(ks,i,j)+dens(ks,i,j+1) ) )
1347  vel = vel * j13g(ks,i,j)
1348  flux(ks,i,j) = vel / mapf(i,j,+2) &
1349  * ( f2 * ( val(ks+1,i,j)+val(ks,i,j) ) &
1350  * ( 0.5_rp + sign(0.5_rp,vel) ) &
1351  + ( 2.0_rp * val(ks,i,j) + 5.0_rp * val(ks+1,i,j) - val(ks+2,i,j) ) / 6.0_rp &
1352  * ( 0.5_rp - sign(0.5_rp,vel) ) )
1353 
1354  vel = ( f2h(ke-1,1,i_xvz) &
1355  * 0.25_rp * ( mom(ke,i,j)+mom(ke,i-1,j)+mom(ke,i,j+1)+mom(ke,i-1,j+1) ) &
1356  + f2h(ke-1,2,i_xvz) &
1357  * 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) ) ) &
1358  / ( f2h(ke-1,1,i_xvz) &
1359  * 0.5_rp * ( dens(ke,i,j)+dens(ke,i,j+1) ) &
1360  + f2h(ke-1,2,i_xvz) &
1361  * 0.5_rp * ( dens(ke-1,i,j)+dens(ke-1,i,j+1) ) )
1362  vel = vel * j13g(ke-1,i,j)
1363  flux(ke-1,i,j) = vel / mapf(i,j,+2) &
1364  * ( ( 2.0_rp * val(ke,i,j) + 5.0_rp * val(ke-1,i,j) - val(ke-2,i,j) ) / 6.0_rp &
1365  * ( 0.5_rp + sign(0.5_rp,vel) ) &
1366  + f2 * ( val(ke,i,j)+val(ke-1,i,j) ) &
1367  * ( 0.5_rp - sign(0.5_rp,vel) ) )
1368 
1369  flux(ke ,i,j) = 0.0_rp
1370  enddo
1371  enddo
1372 
1373  return
1374  end subroutine atmos_dyn_fvm_fluxj13_xvz_ud3
1375 
1376  !-----------------------------------------------------------------------------
1378  subroutine atmos_dyn_fvm_fluxj23_xvz_ud3( &
1379  flux, &
1380  mom, val, DENS, &
1381  GSQRT, J23G, MAPF, &
1382  CDZ, &
1383  IIS, IIE, JJS, JJE )
1384  implicit none
1385 
1386  real(RP), intent(inout) :: flux (ka,ia,ja)
1387  real(RP), intent(in) :: mom (ka,ia,ja)
1388  real(RP), intent(in) :: val (ka,ia,ja)
1389  real(RP), intent(in) :: DENS (ka,ia,ja)
1390  real(RP), intent(in) :: GSQRT (ka,ia,ja)
1391  real(RP), intent(in) :: J23G (ka,ia,ja)
1392  real(RP), intent(in) :: MAPF ( ia,ja,2)
1393  real(RP), intent(in) :: CDZ (ka)
1394  integer, intent(in) :: IIS, IIE, JJS, JJE
1395 
1396  real(RP) :: vel
1397  integer :: k, i, j
1398  !---------------------------------------------------------------------------
1399 
1400  !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
1401  !$omp private(vel) &
1402  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,J23G,MAPF) &
1403  !$omp shared(GSQRT,CDZ)
1404  do j = jjs, jje
1405  do i = iis, iie
1406  do k = ks+1, ke-2
1407  vel = ( f2h(k,1,i_xvz) &
1408  * mom(k+1,i,j) &
1409  + f2h(k,2,i_xvz) &
1410  * mom(k,i,j) ) &
1411  / ( f2h(k,1,i_xvz) &
1412  * 0.5_rp * ( dens(k+1,i,j)+dens(k+1,i,j+1) ) &
1413  + f2h(k,2,i_xvz) &
1414  * 0.5_rp * ( dens(k,i,j)+dens(k,i,j+1) ) )
1415  vel = vel * j23g(k,i,j)
1416  flux(k,i,j) = vel / mapf(i,j,+1) &
1417  * ( ( f31 * ( val(k+2,i,j)+val(k-1,i,j) ) + f32 * ( val(k+1,i,j)+val(k,i,j) ) ) &
1418  - ( f31 * ( val(k+2,i,j)-val(k-1,i,j) ) + f33 * ( val(k+1,i,j)-val(k,i,j) ) ) * sign(1.0_rp,vel) )
1419  enddo
1420  enddo
1421  enddo
1422 
1423  !$omp parallel do default(none) private(i,j) OMP_SCHEDULE_ collapse(2) &
1424  !$omp private(vel) &
1425  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,J23G,MAPF) &
1426  !$omp shared(GSQRT,CDZ)
1427  do j = jjs, jje
1428  do i = iis, iie
1429  ! The boundary condition is qflx_hi + qflxJ13 + qfluxJ23 = 0 at KS-1.
1430  ! The flux at KS-1 can be non-zero.
1431  ! To reduce calculations, all the fluxes are set to zero.
1432  flux(ks-1,i,j) = 0.0_rp
1433 
1434  vel = ( f2h(ks,1,i_xvz) &
1435  * mom(ks+1,i,j) &
1436  + f2h(ks,2,i_xvz) &
1437  * mom(ks,i,j) ) &
1438  / ( f2h(ks,1,i_xvz) &
1439  * 0.5_rp * ( dens(ks+1,i,j)+dens(ks+1,i,j+1) ) &
1440  + f2h(ks,2,i_xvz) &
1441  * 0.5_rp * ( dens(ks,i,j)+dens(ks,i,j+1) ) )
1442  vel = vel * j23g(ks,i,j)
1443  flux(ks,i,j) = vel / mapf(i,j,+1) &
1444  * ( f2 * ( val(ks+1,i,j)+val(ks,i,j) ) &
1445  * ( 0.5_rp + sign(0.5_rp,vel) ) &
1446  + ( 2.0_rp * val(ks,i,j) + 5.0_rp * val(ks+1,i,j) - val(ks+2,i,j) ) / 6.0_rp &
1447  * ( 0.5_rp - sign(0.5_rp,vel) ) )
1448 
1449  vel = ( f2h(ke-1,1,i_xvz) &
1450  * mom(ke,i,j) &
1451  + f2h(ke-1,2,i_xvz) &
1452  * mom(ke-1,i,j) ) &
1453  / ( f2h(ke-1,1,i_xvz) &
1454  * 0.5_rp * ( dens(ke,i,j)+dens(ke,i,j+1) ) &
1455  + f2h(ke-1,2,i_xvz) &
1456  * 0.5_rp * ( dens(ke-1,i,j)+dens(ke-1,i,j+1) ) )
1457  vel = vel * j23g(ke-1,i,j)
1458  flux(ke-1,i,j) = vel / mapf(i,j,+1) &
1459  * ( ( 2.0_rp * val(ke,i,j) + 5.0_rp * val(ke-1,i,j) - val(ke-2,i,j) ) / 6.0_rp &
1460  * ( 0.5_rp + sign(0.5_rp,vel) ) &
1461  + f2 * ( val(ke,i,j)+val(ke-1,i,j) ) &
1462  * ( 0.5_rp - sign(0.5_rp,vel) ) )
1463 
1464  flux(ke ,i,j) = 0.0_rp
1465  enddo
1466  enddo
1467 
1468  return
1469  end subroutine atmos_dyn_fvm_fluxj23_xvz_ud3
1470 
1471  !-----------------------------------------------------------------------------
1473  subroutine atmos_dyn_fvm_fluxx_xvz_ud3( &
1474  flux, &
1475  mom, val, DENS, &
1476  GSQRT, MAPF, &
1477  num_diff, &
1478  CDZ, &
1479  IIS, IIE, JJS, JJE )
1480  implicit none
1481 
1482  real(RP), intent(inout) :: flux (ka,ia,ja)
1483  real(RP), intent(in) :: mom (ka,ia,ja)
1484  real(RP), intent(in) :: val (ka,ia,ja)
1485  real(RP), intent(in) :: DENS (ka,ia,ja)
1486  real(RP), intent(in) :: GSQRT (ka,ia,ja)
1487  real(RP), intent(in) :: MAPF ( ia,ja,2)
1488  real(RP), intent(in) :: num_diff(ka,ia,ja)
1489  real(RP), intent(in) :: CDZ (ka)
1490  integer, intent(in) :: IIS, IIE, JJS, JJE
1491 
1492  real(RP) :: vel
1493  integer :: k, i, j
1494  !---------------------------------------------------------------------------
1495 
1496  !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
1497  !$omp private(vel) &
1498  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,GSQRT,MAPF,num_diff)
1499  do j = jjs, jje
1500  do i = iis-1, iie
1501  do k = ks, ke
1502 #ifdef DEBUG
1503  call check( __line__, mom(k,i ,j) )
1504  call check( __line__, mom(k,i,j-1) )
1505 
1506  call check( __line__, val(k,i,j) )
1507  call check( __line__, val(k,i+1,j) )
1508 
1509  call check( __line__, val(k,i-1,j) )
1510  call check( __line__, val(k,i+2,j) )
1511 
1512 #endif
1513  vel = ( 0.5_rp * ( mom(k,i,j)+mom(k,i,j+1) ) ) &
1514  / ( 0.25_rp * ( dens(k,i,j)+dens(k,i+1,j)+dens(k,i,j+1)+dens(k,i+1,j+1) ) )
1515  flux(k,i,j) = gsqrt(k,i,j) / mapf(i,j,+2) * vel &
1516  * ( ( f31 * ( val(k,i+2,j)+val(k,i-1,j) ) + f32 * ( val(k,i+1,j)+val(k,i,j) ) ) &
1517  - ( f31 * ( val(k,i+2,j)-val(k,i-1,j) ) + f33 * ( val(k,i+1,j)-val(k,i,j) ) ) * sign(1.0_rp,vel) ) &
1518  + gsqrt(k,i,j) * num_diff(k,i,j)
1519  enddo
1520  enddo
1521  enddo
1522 #ifdef DEBUG
1523  k = iundef; i = iundef; j = iundef
1524 #endif
1525 
1526  return
1527  end subroutine atmos_dyn_fvm_fluxx_xvz_ud3
1528 
1529  !-----------------------------------------------------------------------------
1531  subroutine atmos_dyn_fvm_fluxy_xvz_ud3( &
1532  flux, &
1533  mom, val, DENS, &
1534  GSQRT, MAPF, &
1535  num_diff, &
1536  CDZ, &
1537  IIS, IIE, JJS, JJE )
1538  implicit none
1539 
1540  real(RP), intent(inout) :: flux (ka,ia,ja)
1541  real(RP), intent(in) :: mom (ka,ia,ja)
1542  real(RP), intent(in) :: val (ka,ia,ja)
1543  real(RP), intent(in) :: DENS (ka,ia,ja)
1544  real(RP), intent(in) :: GSQRT (ka,ia,ja)
1545  real(RP), intent(in) :: MAPF ( ia,ja,2)
1546  real(RP), intent(in) :: num_diff(ka,ia,ja)
1547  real(RP), intent(in) :: CDZ (ka)
1548  integer, intent(in) :: IIS, IIE, JJS, JJE
1549 
1550  real(RP) :: vel
1551  integer :: k, i, j
1552  !---------------------------------------------------------------------------
1553 
1554  ! note that y-index is added by -1
1555 
1556  !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
1557  !$omp private(vel) &
1558  !$omp shared(JJS,JJE,IIS,IIE,KS,KE,mom,val,DENS,flux,GSQRT,MAPF,num_diff)
1559  do j = jjs, jje+1
1560  do i = iis, iie
1561  do k = ks, ke
1562 #ifdef DEBUG
1563  call check( __line__, mom(k,i ,j) )
1564  call check( __line__, mom(k,i,j-1) )
1565 
1566  call check( __line__, val(k,i,j-1) )
1567  call check( __line__, val(k,i,j) )
1568 
1569  call check( __line__, val(k,i,j-2) )
1570  call check( __line__, val(k,i,j+1) )
1571 
1572 #endif
1573  vel = ( 0.5_rp * ( mom(k,i,j)+mom(k,i,j-1) ) ) &
1574  / ( dens(k,i,j) )
1575  flux(k,i,j-1) = gsqrt(k,i,j) / mapf(i,j,+1) * vel &
1576  * ( ( f31 * ( val(k,i,j+1)+val(k,i,j-2) ) + f32 * ( val(k,i,j)+val(k,i,j-1) ) ) &
1577  - ( f31 * ( val(k,i,j+1)-val(k,i,j-2) ) + f33 * ( val(k,i,j)-val(k,i,j-1) ) ) * sign(1.0_rp,vel) ) &
1578  + gsqrt(k,i,j) * num_diff(k,i,j)
1579  enddo
1580  enddo
1581  enddo
1582 #ifdef DEBUG
1583  k = iundef; i = iundef; j = iundef
1584 #endif
1585 
1586  return
1587  end subroutine atmos_dyn_fvm_fluxy_xvz_ud3
1588 
1589 
1590 
1591 
1592 
1593 
1594 
1595 end module scale_atmos_dyn_fvm_flux_ud3
1596 
1597 !--
1598 ! vi:set readonly sw=4 ts=8
1599 !
1600 !Local Variables:
1601 !mode: f90
1602 !buffer-read-only: t
1603 !End:
1604 !
1605 !++
scale_debug
module DEBUG
Definition: scale_debug.F90:11
scale_atmos_dyn_fvm_flux_ud3::atmos_dyn_fvm_fluxj13_xvz_ud3
subroutine, public atmos_dyn_fvm_fluxj13_xvz_ud3(flux, mom, val, DENS, GSQRT, J13G, MAPF, CDZ, IIS, IIE, JJS, JJE)
calculation J13-flux at XVZ
Definition: scale_atmos_dyn_fvm_flux_ud3.F90:1289
scale_atmos_grid_cartesc_index::ia
integer, public ia
of whole cells: x, local, with HALO
Definition: scale_atmos_grid_cartesC_index.F90:48
scale_atmos_dyn_fvm_flux_ud3::atmos_dyn_fvm_fluxj23_uyz_ud3
subroutine, public atmos_dyn_fvm_fluxj23_uyz_ud3(flux, mom, val, DENS, GSQRT, J23G, MAPF, CDZ, IIS, IIE, JJS, JJE)
calculation J23-flux at UYZ
Definition: scale_atmos_dyn_fvm_flux_ud3.F90:964
scale_atmos_grid_cartesc_index::i_xvz
integer, public i_xvz
Definition: scale_atmos_grid_cartesC_index.F90:95
scale_atmos_dyn_fvm_flux_ud3::atmos_dyn_fvm_fluxz_xyz_ud3
subroutine, public atmos_dyn_fvm_fluxz_xyz_ud3(flux, mflx, val, GSQRT, num_diff, CDZ, IIS, IIE, JJS, JJE)
calculation z-flux at XYZ
Definition: scale_atmos_dyn_fvm_flux_ud3.F90:163
scale_atmos_dyn_fvm_flux_ud3::atmos_dyn_fvm_flux_valuew_z_ud3
subroutine, public atmos_dyn_fvm_flux_valuew_z_ud3(valW, mflx, val, GSQRT, CDZ)
value at XYW
Definition: scale_atmos_dyn_fvm_flux_ud3.F90:102
scale_atmos_grid_cartesc_index::ja
integer, public ja
of whole cells: y, local, with HALO
Definition: scale_atmos_grid_cartesC_index.F90:49
scale_debug::check
subroutine, public check(current_line, v)
Undefined value checker.
Definition: scale_debug.F90:56
scale_atmos_dyn_fvm_flux_ud3::atmos_dyn_fvm_fluxx_xyw_ud3
subroutine, public atmos_dyn_fvm_fluxx_xyw_ud3(flux, mom, val, DENS, GSQRT, MAPF, num_diff, CDZ, IIS, IIE, JJS, JJE)
calculation X-flux at XYW
Definition: scale_atmos_dyn_fvm_flux_ud3.F90:605
scale_const::const_undef
real(rp), public const_undef
Definition: scale_const.F90:41
scale_tracer
module TRACER
Definition: scale_tracer.F90:16
scale_index
module Index
Definition: scale_index.F90:11
scale_atmos_dyn_fvm_flux_ud3::atmos_dyn_fvm_fluxj23_xvz_ud3
subroutine, public atmos_dyn_fvm_fluxj23_xvz_ud3(flux, mom, val, DENS, GSQRT, J23G, MAPF, CDZ, IIS, IIE, JJS, JJE)
calculation J23-flux at XVZ
Definition: scale_atmos_dyn_fvm_flux_ud3.F90:1384
scale_atmos_dyn_fvm_flux_ud3::atmos_dyn_fvm_fluxx_xvz_ud3
subroutine, public atmos_dyn_fvm_fluxx_xvz_ud3(flux, mom, val, DENS, GSQRT, MAPF, num_diff, CDZ, IIS, IIE, JJS, JJE)
calculation X-flux at XV
Definition: scale_atmos_dyn_fvm_flux_ud3.F90:1480
scale_atmos_grid_cartesc_index
module atmosphere / grid / cartesC index
Definition: scale_atmos_grid_cartesC_index.F90:12
scale_atmos_grid_cartesc_index::ke
integer, public ke
end point of inner domain: z, local
Definition: scale_atmos_grid_cartesC_index.F90:52
scale_atmos_dyn_fvm_flux_ud3::atmos_dyn_fvm_fluxz_xvz_ud3
subroutine, public atmos_dyn_fvm_fluxz_xvz_ud3(flux, mom, val, DENS, GSQRT, J33G, num_diff, CDZ, IIS, IIE, JJS, JJE)
calculation z-flux at XV
Definition: scale_atmos_dyn_fvm_flux_ud3.F90:1180
scale_atmos_dyn_fvm_flux_ud3::atmos_dyn_fvm_fluxx_xyz_ud3
subroutine, public atmos_dyn_fvm_fluxx_xyz_ud3(flux, mflx, val, GSQRT, num_diff, CDZ, IIS, IIE, JJS, JJE)
calculation X-flux at XYZ
Definition: scale_atmos_dyn_fvm_flux_ud3.F90:257
scale_prc
module PROCESS
Definition: scale_prc.F90:11
scale_atmos_dyn_fvm_flux_ud3::atmos_dyn_fvm_fluxy_xyw_ud3
subroutine, public atmos_dyn_fvm_fluxy_xyw_ud3(flux, mom, val, DENS, GSQRT, MAPF, num_diff, CDZ, IIS, IIE, JJS, JJE)
calculation Y-flux at XYW
Definition: scale_atmos_dyn_fvm_flux_ud3.F90:682
scale_const::const_undef2
integer, parameter, public const_undef2
undefined value (INT2)
Definition: scale_const.F90:38
scale_atmos_grid_cartesc_index::ks
integer, public ks
start point of inner domain: z, local
Definition: scale_atmos_grid_cartesC_index.F90:51
scale_atmos_dyn_fvm_flux_ud3::atmos_dyn_fvm_fluxy_xyz_ud3
subroutine, public atmos_dyn_fvm_fluxy_xyz_ud3(flux, mflx, val, GSQRT, num_diff, CDZ, IIS, IIE, JJS, JJE)
calculation Y-flux at XYZ
Definition: scale_atmos_dyn_fvm_flux_ud3.F90:310
scale_atmos_dyn_fvm_flux_ud3::atmos_dyn_fvm_fluxj23_xyw_ud3
subroutine, public atmos_dyn_fvm_fluxj23_xyw_ud3(flux, mom, val, DENS, GSQRT, J23G, MAPF, CDZ, IIS, IIE, JJS, JJE)
calculation J23-flux at XYW
Definition: scale_atmos_dyn_fvm_flux_ud3.F90:533
scale_const
module CONSTANT
Definition: scale_const.F90:11
scale_atmos_dyn_fvm_flux_ud3::atmos_dyn_fvm_fluxy_xvz_ud3
subroutine, public atmos_dyn_fvm_fluxy_xvz_ud3(flux, mom, val, DENS, GSQRT, MAPF, num_diff, CDZ, IIS, IIE, JJS, JJE)
calculation Y-flux at XV
Definition: scale_atmos_dyn_fvm_flux_ud3.F90:1538
scale_atmos_dyn_fvm_flux_ud3::atmos_dyn_fvm_fluxy_uyz_ud3
subroutine, public atmos_dyn_fvm_fluxy_uyz_ud3(flux, mom, val, DENS, GSQRT, MAPF, num_diff, CDZ, IIS, IIE, JJS, JJE)
calculation Y-flux at UY
Definition: scale_atmos_dyn_fvm_flux_ud3.F90:1120
scale_atmos_grid_cartesc_index::i_uyz
integer, public i_uyz
Definition: scale_atmos_grid_cartesC_index.F90:94
scale_prof
module profiler
Definition: scale_prof.F90:11
scale_atmos_dyn_fvm_flux_ud3::atmos_dyn_fvm_fluxz_uyz_ud3
subroutine, public atmos_dyn_fvm_fluxz_uyz_ud3(flux, mom, val, DENS, GSQRT, J33G, num_diff, CDZ, IIS, IIE, JJS, JJE)
calculation z-flux at UY
Definition: scale_atmos_dyn_fvm_flux_ud3.F90:760
scale_const::const_eps
real(rp), public const_eps
small number
Definition: scale_const.F90:33
scale_atmos_dyn_fvm_flux_ud3::atmos_dyn_fvm_fluxj13_xyw_ud3
subroutine, public atmos_dyn_fvm_fluxj13_xyw_ud3(flux, mom, val, DENS, GSQRT, J13G, MAPF, CDZ, IIS, IIE, JJS, JJE)
calculation J13-flux at XYW
Definition: scale_atmos_dyn_fvm_flux_ud3.F90:463
scale_atmos_dyn_fvm_flux_ud3::atmos_dyn_fvm_fluxx_uyz_ud3
subroutine, public atmos_dyn_fvm_fluxx_uyz_ud3(flux, mom, val, DENS, GSQRT, MAPF, num_diff, CDZ, IIS, IIE, JJS, JJE)
calculation X-flux at UY
Definition: scale_atmos_dyn_fvm_flux_ud3.F90:1060
scale_precision
module PRECISION
Definition: scale_precision.F90:14
scale_atmos_grid_cartesc_index::ka
integer, public ka
of whole cells: z, local, with HALO
Definition: scale_atmos_grid_cartesC_index.F90:47
scale_atmos_dyn_fvm_flux_ud3::atmos_dyn_fvm_fluxj13_uyz_ud3
subroutine, public atmos_dyn_fvm_fluxj13_uyz_ud3(flux, mom, val, DENS, GSQRT, J13G, MAPF, CDZ, IIS, IIE, JJS, JJE)
calculation J13-flux at UYZ
Definition: scale_atmos_dyn_fvm_flux_ud3.F90:869
scale_io
module STDIO
Definition: scale_io.F90:10
scale_atmos_dyn_fvm_flux_ud3::atmos_dyn_fvm_fluxz_xyw_ud3
subroutine, public atmos_dyn_fvm_fluxz_xyw_ud3(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_ud3.F90:366
scale_atmos_dyn_fvm_flux_ud3
module scale_atmos_dyn_fvm_flux_ud3
Definition: scale_atmos_dyn_fvm_flux_ud3.F90:16
Generated by   doxygen 1.8.13