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