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