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