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