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