SCALE-RM
scale_atmos_dyn_fvm_fct.F90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
10 !-------------------------------------------------------------------------------
11 #include "scalelib.h"
12 module scale_atmos_dyn_fvm_fct
13  !-----------------------------------------------------------------------------
14  !
15  !++ used modules
16  !
17  use scale_precision
18  use scale_io
19  use scale_prof
20  use scale_atmos_grid_cartesc_index
21  use scale_index
22  use scale_tracer
23 
24 #ifdef DEBUG
25  use scale_debug, only: &
26  check
27  use scale_const, only: &
28  undef => const_undef, &
29  iundef => const_undef2
30 #endif
31  !-----------------------------------------------------------------------------
32  implicit none
33  private
34  !-----------------------------------------------------------------------------
35  !
36  !++ Public procedure
37  !
38  public :: atmos_dyn_fvm_fct
39 
40  !-----------------------------------------------------------------------------
41  !
42  !++ Public parameters & variables
43  !
44  !-----------------------------------------------------------------------------
45  !
46  !++ Private procedure
47  !
48  private :: get_fact_fct
49 
50  !-----------------------------------------------------------------------------
51  !
52  !++ Private parameters & variables
53  !
54  !-----------------------------------------------------------------------------
55 
56 contains
57 
58  !-----------------------------------------------------------------------------
60 
61  !-----------------------------------------------------------------------------
63  subroutine atmos_dyn_fvm_fct( &
64  qflx_anti, &
65  phi_in, DENS0, DENS, &
66  qflx_hi, qflx_lo, &
67  mflx_hi, &
68  rdz, rdx, rdy, &
69  GSQRT, MAPF, &
70  TwoD, dt, &
71  flag_vect )
72  use scale_const, only: &
73  undef => const_undef, &
74  iundef => const_undef2, &
75  epsilon => const_eps
76  use scale_comm_cartesc, only: &
77  comm_vars8, &
78  comm_wait
79  implicit none
80 
81  real(rp), intent(out) :: qflx_anti(ka,ia,ja,3)
82 
83  real(rp), intent(in) :: phi_in(ka,ia,ja) ! physical quantity
84  real(rp), intent(in) :: dens0(ka,ia,ja)
85  real(rp), intent(in) :: dens (ka,ia,ja)
86 
87  real(rp), intent(in) :: qflx_hi(ka,ia,ja,3)
88  real(rp), intent(in) :: qflx_lo(ka,ia,ja,3)
89  real(rp), intent(in) :: mflx_hi(ka,ia,ja,3)
90 
91  real(rp), intent(in) :: rdz(:)
92  real(rp), intent(in) :: rdx(:)
93  real(rp), intent(in) :: rdy(:)
94 
95  real(rp), intent(in) :: gsqrt(ka,ia,ja)
96  real(rp), intent(in) :: mapf(ia,ja,2)
97 
98  logical, intent(in) :: twod
99  real(rp), intent(in) :: dt
100 
101  logical, intent(in) :: flag_vect
102 
103  ! work for FCT
104  real(rp) :: phi_lo(ka,ia,ja)
105  real(rp) :: pjpls(ka,ia,ja)
106  real(rp) :: pjmns(ka,ia,ja)
107  real(rp) :: qjpls(ka,ia,ja)
108  real(rp) :: qjmns(ka,ia,ja)
109  real(rp) :: rjpls(ka,ia,ja)
110  real(rp) :: rjmns(ka,ia,ja)
111 
112  real(rp) :: qmin, qmax
113  real(rp) :: zerosw, dirsw
114 
115  real(rp) :: fact(0:1,-1:1,-1:1)
116  real(rp) :: rw, ru, rv
117  real(rp) :: qa_in, qb_in
118  real(rp) :: qa_lo, qb_lo
119 
120  integer :: k, i, j, ijs
121  integer :: iis, iie, jjs, jje
122  !---------------------------------------------------------------------------
123 
124 #ifdef DEBUG
125  qflx_anti(:,:,:,:) = undef
126 
127  pjpls(:,:,:) = undef
128  pjmns(:,:,:) = undef
129  qjpls(:,:,:) = undef
130  qjmns(:,:,:) = undef
131  rjpls(:,:,:) = undef
132  rjmns(:,:,:) = undef
133 #endif
134 
135  do jjs = js, je, jblock
136  jje = jjs+jblock-1
137  do iis = is, ie, iblock
138  iie = iis+iblock-1
139 
140  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
141  do j = jjs, jje
142  do i = iis, iie
143  do k = ks, ke-1
144 #ifdef DEBUG
145  call check( __line__, qflx_hi(k,i,j,zdir) )
146  call check( __line__, qflx_lo(k,i,j,zdir) )
147 #endif
148  qflx_anti(k,i,j,zdir) = qflx_hi(k,i,j,zdir) - qflx_lo(k,i,j,zdir)
149  enddo
150  enddo
151  enddo
152 #ifdef DEBUG
153  k = iundef; i = iundef; j = iundef
154 #endif
155  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
156  do j = jjs, jje
157  do i = iis, iie
158  qflx_anti(ks-1,i,j,zdir) = 0.0_rp
159  qflx_anti(ke ,i,j,zdir) = 0.0_rp
160  enddo
161  enddo
162 #ifdef DEBUG
163  k = iundef; i = iundef; j = iundef
164 #endif
165 
166  if ( .not. twod ) then
167  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
168  do j = jjs , jje
169  do i = iis-1, iie
170  do k = ks, ke
171 #ifdef DEBUG
172  call check( __line__, qflx_hi(k,i,j,xdir) )
173  call check( __line__, qflx_lo(k,i,j,xdir) )
174 #endif
175  qflx_anti(k,i,j,xdir) = qflx_hi(k,i,j,xdir) - qflx_lo(k,i,j,xdir)
176  enddo
177  enddo
178  enddo
179 #ifdef DEBUG
180  k = iundef; i = iundef; j = iundef
181 #endif
182  end if
183 
184  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
185  do j = jjs-1, jje
186  do i = iis , iie
187  do k = ks, ke
188 #ifdef DEBUG
189  call check( __line__, qflx_hi(k,i,j,ydir) )
190  call check( __line__, qflx_lo(k,i,j,ydir) )
191 #endif
192  qflx_anti(k,i,j,ydir) = qflx_hi(k,i,j,ydir) - qflx_lo(k,i,j,ydir)
193  enddo
194  enddo
195  enddo
196 #ifdef DEBUG
197  k = iundef; i = iundef; j = iundef
198 #endif
199 
200  !--- update monotone scheme
201  if ( twod ) then
202  !$omp parallel do private(j,k) OMP_SCHEDULE_
203  do j = jjs-1, jje+1
204  do k = ks, ke
205 #ifdef DEBUG
206  call check( __line__, phi_in(k,is,j) )
207  call check( __line__, qflx_lo(k ,is,j ,zdir) )
208  call check( __line__, qflx_lo(k-1,is,j ,zdir) )
209  call check( __line__, qflx_lo(k ,is,j ,ydir) )
210  call check( __line__, qflx_lo(k ,is,j-1,ydir) )
211 #endif
212  phi_lo(k,is,j) = ( phi_in(k,is,j) * dens0(k,is,j) &
213  + dt * ( - ( ( qflx_lo(k,is,j,zdir)-qflx_lo(k-1,is,j ,zdir) ) * rdz(k) &
214  + ( qflx_lo(k,is,j,ydir)-qflx_lo(k ,is,j-1,ydir) ) * rdy(j) &
215  ) * mapf(is,j,2) / gsqrt(k,is,j) ) &
216  ) / dens(k,is,j)
217  enddo
218  enddo
219  else
220  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
221  do j = jjs-1, jje+1
222  do i = iis-1, iie+1
223  do k = ks, ke
224 #ifdef DEBUG
225  call check( __line__, phi_in(k,i,j) )
226  call check( __line__, qflx_lo(k ,i ,j ,zdir) )
227  call check( __line__, qflx_lo(k-1,i ,j ,zdir) )
228  call check( __line__, qflx_lo(k ,i ,j ,xdir) )
229  call check( __line__, qflx_lo(k ,i-1,j ,xdir) )
230  call check( __line__, qflx_lo(k ,i ,j ,ydir) )
231  call check( __line__, qflx_lo(k ,i ,j-1,ydir) )
232 #endif
233  phi_lo(k,i,j) = ( phi_in(k,i,j) * dens0(k,i,j) &
234  + dt * ( - ( ( qflx_lo(k,i,j,zdir)-qflx_lo(k-1,i ,j ,zdir) ) * rdz(k) &
235  + ( qflx_lo(k,i,j,xdir)-qflx_lo(k ,i-1,j ,xdir) ) * rdx(i) &
236  + ( qflx_lo(k,i,j,ydir)-qflx_lo(k ,i ,j-1,ydir) ) * rdy(j) &
237  ) * mapf(i,j,1) * mapf(i,j,2) / gsqrt(k,i,j) ) &
238  ) / dens(k,i,j)
239  enddo
240  enddo
241  enddo
242  end if
243 #ifdef DEBUG
244  k = iundef; i = iundef; j = iundef
245 #endif
246 
247  !--- calc net incoming quantity change by antidiffusive flux
248  if ( twod ) then
249  !$omp parallel do private(j,k) OMP_SCHEDULE_
250  do j = jjs, jje
251  do k = ks, ke
252 #ifdef DEBUG
253  call check( __line__, qflx_anti(k ,is,j ,zdir) )
254  call check( __line__, qflx_anti(k-1,is,j ,zdir) )
255  call check( __line__, qflx_anti(k ,is,j ,ydir) )
256  call check( __line__, qflx_anti(k ,is,j-1,ydir) )
257 #endif
258  pjpls(k,is,j) = dt * ( ( max(0.0_rp,qflx_anti(k-1,is ,j ,zdir)) - min(0.0_rp,qflx_anti(k,is,j,zdir)) ) * rdz(k) &
259  + ( max(0.0_rp,qflx_anti(k ,is,j-1,ydir)) - min(0.0_rp,qflx_anti(k,is,j,ydir)) ) * rdy(j) &
260  ) * mapf(is,j,2) / gsqrt(k,is,j)
261  enddo
262  enddo
263  else
264  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
265  do j = jjs, jje
266  do i = iis, iie
267  do k = ks, ke
268 #ifdef DEBUG
269  call check( __line__, qflx_anti(k ,i ,j ,zdir) )
270  call check( __line__, qflx_anti(k-1,i ,j ,zdir) )
271  call check( __line__, qflx_anti(k ,i ,j ,xdir) )
272  call check( __line__, qflx_anti(k ,i-1,j ,xdir) )
273  call check( __line__, qflx_anti(k ,i ,j ,ydir) )
274  call check( __line__, qflx_anti(k ,i ,j-1,ydir) )
275 #endif
276  pjpls(k,i,j) = dt * ( ( max(0.0_rp,qflx_anti(k-1,i ,j ,zdir)) - min(0.0_rp,qflx_anti(k,i,j,zdir)) ) * rdz(k) &
277  + ( max(0.0_rp,qflx_anti(k ,i-1,j ,xdir)) - min(0.0_rp,qflx_anti(k,i,j,xdir)) ) * rdx(i) &
278  + ( max(0.0_rp,qflx_anti(k ,i ,j-1,ydir)) - min(0.0_rp,qflx_anti(k,i,j,ydir)) ) * rdy(j) &
279  ) * mapf(i,j,1) * mapf(i,j,2) / gsqrt(k,i,j)
280  enddo
281  enddo
282  enddo
283  end if
284 #ifdef DEBUG
285  k = iundef; i = iundef; j = iundef
286 #endif
287 
288  !--- calc net outgoing quantity change by antidiffusive flux
289  if ( twod ) then
290  !$omp parallel do private(j,k) OMP_SCHEDULE_
291  do j = jjs, jje
292  do k = ks, ke
293 #ifdef DEBUG
294  call check( __line__, qflx_anti(k ,is,j ,zdir) )
295  call check( __line__, qflx_anti(k-1,is,j ,zdir) )
296  call check( __line__, qflx_anti(k ,is,j ,ydir) )
297  call check( __line__, qflx_anti(k ,is,j-1,ydir) )
298 #endif
299  pjmns(k,is,j) = dt * ( ( max(0.0_rp,qflx_anti(k,is,j,zdir)) - min(0.0_rp,qflx_anti(k-1,is,j ,zdir)) ) * rdz(k) &
300  + ( max(0.0_rp,qflx_anti(k,is,j,ydir)) - min(0.0_rp,qflx_anti(k ,is,j-1,ydir)) ) * rdy(j) &
301  ) * mapf(is,j,2) / gsqrt(k,is,j)
302  enddo
303  enddo
304  else
305  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
306  do j = jjs, jje
307  do i = iis, iie
308  do k = ks, ke
309 #ifdef DEBUG
310  call check( __line__, qflx_anti(k ,i ,j ,zdir) )
311  call check( __line__, qflx_anti(k-1,i ,j ,zdir) )
312  call check( __line__, qflx_anti(k ,i ,j ,xdir) )
313  call check( __line__, qflx_anti(k ,i-1,j ,xdir) )
314  call check( __line__, qflx_anti(k ,i ,j ,ydir) )
315  call check( __line__, qflx_anti(k ,i ,j-1,ydir) )
316 #endif
317  pjmns(k,i,j) = dt * ( ( max(0.0_rp,qflx_anti(k,i,j,zdir)) - min(0.0_rp,qflx_anti(k-1,i ,j ,zdir)) ) * rdz(k) &
318  + ( max(0.0_rp,qflx_anti(k,i,j,xdir)) - min(0.0_rp,qflx_anti(k ,i-1,j ,xdir)) ) * rdx(i) &
319  + ( max(0.0_rp,qflx_anti(k,i,j,ydir)) - min(0.0_rp,qflx_anti(k ,i ,j-1,ydir)) ) * rdy(j) &
320  ) * mapf(i,j,1) * mapf(i,j,2) / gsqrt(k,i,j)
321  enddo
322  enddo
323  enddo
324  end if
325 #ifdef DEBUG
326  k = iundef; i = iundef; j = iundef
327 #endif
328 
329  !--- calc allowable range or quantity change by antidiffusive flux
330 
331  if (flag_vect) then
332 
333  if ( twod ) then
334  i = is
335  !$omp parallel do private(j,k,rw,ru,rv,fact,qa_in,qb_in,qa_lo,qb_lo,qmax,qmin) OMP_SCHEDULE_
336  do j = jjs, jje
337  do k = ks+1, ke-1
338  rw = (mflx_hi(k,i,j,zdir)+mflx_hi(k-1,i ,j ,zdir)) * rdz(k) ! 2 * rho * w / dz
339  ru = 0.0_rp
340  rv = (mflx_hi(k,i,j,ydir)+mflx_hi(k ,i ,j-1,ydir)) * rdy(j) ! 2 * rho * v / dy
341 
342  call get_fact_fct( fact, & ! (out)
343  rw, ru, rv ) ! (in)
344 
345  qa_in = fact(1, 0, 1) * phi_in(k+1,i ,j+1) &
346  + fact(0, 0, 1) * phi_in(k ,i ,j+1) &
347  + fact(1, 0, 0) * phi_in(k+1,i ,j ) &
348  + fact(1, 0,-1) * phi_in(k+1,i ,j-1) &
349  + fact(0, 0, 0) * phi_in(k ,i ,j )
350  qb_in = fact(1, 0, 1) * phi_in(k-1,i ,j-1) &
351  + fact(0, 0, 1) * phi_in(k ,i ,j-1) &
352  + fact(1, 0, 0) * phi_in(k-1,i ,j ) &
353  + fact(1, 0,-1) * phi_in(k-1,i ,j-1) &
354  + fact(0, 0, 0) * phi_in(k ,i ,j )
355  qa_lo = fact(1, 0, 1) * phi_lo(k+1,i ,j+1) &
356  + fact(0, 0, 1) * phi_lo(k ,i ,j+1) &
357  + fact(1, 0, 0) * phi_lo(k+1,i ,j ) &
358  + fact(1, 0,-1) * phi_lo(k+1,i ,j-1) &
359  + fact(0, 0, 0) * phi_lo(k ,i ,j )
360  qb_lo = fact(1, 0, 1) * phi_lo(k-1,i ,j-1) &
361  + fact(0, 0, 1) * phi_lo(k ,i ,j-1) &
362  + fact(1, 0, 0) * phi_lo(k-1,i ,j ) &
363  + fact(1, 0,-1) * phi_lo(k-1,i ,j-1) &
364  + fact(0, 0, 0) * phi_lo(k ,i ,j )
365 
366  qmax = max( &
367  phi_in(k,i,j), qa_in, qb_in, &
368  phi_lo(k,i,j), qa_lo, qb_lo )
369  qmin = min( &
370  phi_in(k,i,j), qa_in, qb_in, &
371  phi_lo(k,i,j), qa_lo, qb_lo )
372  qjpls(k,i,j) = ( qmax - phi_lo(k,i,j) ) * dens(k,i,j)
373  qjmns(k,i,j) = ( phi_lo(k,i,j) - qmin ) * dens(k,i,j)
374  end do
375  end do
376 
377  ! k = KS
378  !$omp parallel do private(j,rw,ru,rv,fact,qa_in,qb_in,qa_lo,qb_lo,qmax,qmin) OMP_SCHEDULE_
379  do j = jjs, jje
380  rw = (mflx_hi(ks,i,j,zdir) ) * rdz(ks)! 2 * rho * w / dz
381  ru = 0.0_rp
382  rv = (mflx_hi(ks,i,j,ydir)+mflx_hi(ks ,i ,j-1,ydir)) * rdy(j) ! 2 * rho * v / dy
383 
384  call get_fact_fct( fact, & ! (out)
385  rw, ru, rv ) ! (in)
386 
387  qa_in = fact(1, 0, 1) * phi_in(ks+1,i ,j+1) &
388  + fact(0, 0, 1) * phi_in(ks ,i ,j+1) &
389  + fact(1, 0, 0) * phi_in(ks+1,i ,j ) &
390  + fact(1, 0,-1) * phi_in(ks+1,i ,j-1) &
391  + fact(0, 0, 0) * phi_in(ks ,i ,j )
392  qb_in = fact(1, 0, 1) * phi_in(ks ,i ,j-1) &
393  + fact(0, 0, 1) * phi_in(ks ,i ,j-1) &
394  + fact(1, 0, 0) * phi_in(ks ,i ,j ) &
395  + fact(1, 0,-1) * phi_in(ks ,i ,j-1) &
396  + fact(0, 0, 0) * phi_in(ks ,i ,j )
397  qa_lo = fact(1, 0, 1) * phi_lo(ks+1,i ,j+1) &
398  + fact(0, 0, 1) * phi_lo(ks ,i ,j+1) &
399  + fact(1, 0, 0) * phi_lo(ks+1,i ,j ) &
400  + fact(1, 0,-1) * phi_lo(ks+1,i ,j-1) &
401  + fact(0, 0, 0) * phi_lo(ks ,i ,j )
402  qb_lo = fact(1, 0, 1) * phi_lo(ks ,i ,j-1) &
403  + fact(0, 0, 1) * phi_lo(ks ,i ,j-1) &
404  + fact(1, 0, 0) * phi_lo(ks ,i ,j ) &
405  + fact(1, 0,-1) * phi_lo(ks ,i ,j-1) &
406  + fact(0, 0, 0) * phi_lo(ks ,i ,j )
407 
408  qmax = max( &
409  phi_in(ks,i,j), qa_in, qb_in, &
410  phi_lo(ks,i,j), qa_lo, qb_lo )
411  qmin = min( &
412  phi_in(ks,i,j), qa_in, qb_in, &
413  phi_lo(ks,i,j), qa_lo, qb_lo )
414  qjpls(ks,i,j) = ( qmax - phi_lo(ks,i,j) ) * dens(ks,i,j)
415  qjmns(ks,i,j) = ( phi_lo(ks,i,j) - qmin ) * dens(ks,i,j)
416  end do
417 
418  ! k = KE
419  !$omp parallel do private(j,rw,ru,rv,fact,qa_in,qb_in,qa_lo,qb_lo,qmax,qmin) OMP_SCHEDULE_
420  do j = jjs, jje
421  rw = ( mflx_hi(ke-1,i ,j ,zdir)) * rdz(ke)! 2 * rho * w / dz
422  ru = 0.0_rp
423  rv = (mflx_hi(ke,i,j,ydir)+mflx_hi(ke ,i ,j-1,ydir)) * rdy(j) ! 2 * rho * v / dy
424 
425  call get_fact_fct( fact, & ! (out)
426  rw, ru, rv ) ! (in)
427 
428  qa_in = fact(1, 0, 1) * phi_in(ke ,i ,j+1) &
429  + fact(0, 0, 1) * phi_in(ke ,i ,j+1) &
430  + fact(1, 0, 0) * phi_in(ke ,i ,j ) &
431  + fact(1, 0,-1) * phi_in(ke ,i ,j-1) &
432  + fact(0, 0, 0) * phi_in(ke ,i ,j )
433  qb_in = fact(1, 0, 1) * phi_in(ke-1,i ,j-1) &
434  + fact(0, 0, 1) * phi_in(ke ,i ,j-1) &
435  + fact(1, 0, 0) * phi_in(ke-1,i ,j ) &
436  + fact(1, 0,-1) * phi_in(ke-1,i ,j-1) &
437  + fact(0, 0, 0) * phi_in(ke ,i ,j )
438  qa_lo = fact(1, 0, 1) * phi_lo(ke ,i ,j+1) &
439  + fact(0, 0, 1) * phi_lo(ke ,i ,j+1) &
440  + fact(1, 0, 0) * phi_lo(ke ,i ,j ) &
441  + fact(1, 0,-1) * phi_lo(ke ,i ,j-1) &
442  + fact(0, 0, 0) * phi_lo(ke ,i ,j )
443  qb_lo = fact(1, 0, 1) * phi_lo(ke-1,i ,j-1) &
444  + fact(0, 0, 1) * phi_lo(ke ,i ,j-1) &
445  + fact(1, 0, 0) * phi_lo(ke-1,i ,j ) &
446  + fact(1, 0,-1) * phi_lo(ke-1,i ,j-1) &
447  + fact(0, 0, 0) * phi_lo(ke ,i ,j )
448 
449  qmax = max( &
450  phi_in(ke,i,j), qa_in, qb_in, &
451  phi_lo(ke,i,j), qa_lo, qb_lo )
452  qmin = min( &
453  phi_in(ke,i,j), qa_in, qb_in, &
454  phi_lo(ke,i,j), qa_lo, qb_lo )
455  qjpls(ke,i,j) = ( qmax - phi_lo(ke,i,j) ) * dens(ke,i,j)
456  qjmns(ke,i,j) = ( phi_lo(ke,i,j) - qmin ) * dens(ke,i,j)
457  end do
458 
459  else
460  !$omp parallel do private(i,j,k,rw,ru,rv,fact,qa_in,qb_in,qa_lo,qb_lo,qmax,qmin) OMP_SCHEDULE_ collapse(2)
461  do j = jjs, jje
462  do i = iis, iie
463  do k = ks+1, ke-1
464  rw = (mflx_hi(k,i,j,zdir)+mflx_hi(k-1,i ,j ,zdir)) * rdz(k) ! 2 * rho * w / dz
465  ru = (mflx_hi(k,i,j,xdir)+mflx_hi(k ,i-1,j ,xdir)) * rdx(i) ! 2 * rho * u / dx
466  rv = (mflx_hi(k,i,j,ydir)+mflx_hi(k ,i ,j-1,ydir)) * rdy(j) ! 2 * rho * v / dy
467 
468  call get_fact_fct( fact, & ! (out)
469  rw, ru, rv ) ! (in)
470 
471  qa_in = fact(1, 1, 1) * phi_in(k+1,i+1,j+1) &
472  + fact(0, 1, 1) * phi_in(k ,i+1,j+1) &
473  + fact(1, 0, 1) * phi_in(k+1,i ,j+1) &
474  + fact(0, 0, 1) * phi_in(k ,i ,j+1) &
475  + fact(1,-1, 1) * phi_in(k+1,i-1,j+1) &
476  + fact(1, 1, 0) * phi_in(k+1,i+1,j ) &
477  + fact(0, 1, 0) * phi_in(k ,i+1,j ) &
478  + fact(1, 0, 0) * phi_in(k+1,i ,j ) &
479  + fact(1,-1, 0) * phi_in(k+1,i-1,j ) &
480  + fact(1, 1,-1) * phi_in(k+1,i+1,j-1) &
481  + fact(0, 1,-1) * phi_in(k ,i+1,j-1) &
482  + fact(1, 0,-1) * phi_in(k+1,i ,j-1) &
483  + fact(1,-1,-1) * phi_in(k+1,i-1,j-1) &
484  + fact(0, 0, 0) * phi_in(k ,i ,j )
485  qb_in = fact(1, 1, 1) * phi_in(k-1,i-1,j-1) &
486  + fact(0, 1, 1) * phi_in(k ,i-1,j-1) &
487  + fact(1, 0, 1) * phi_in(k-1,i ,j-1) &
488  + fact(0, 0, 1) * phi_in(k ,i ,j-1) &
489  + fact(1,-1, 1) * phi_in(k-1,i+1,j-1) &
490  + fact(1, 1, 0) * phi_in(k-1,i-1,j ) &
491  + fact(0, 1, 0) * phi_in(k ,i-1,j ) &
492  + fact(1, 0, 0) * phi_in(k-1,i ,j ) &
493  + fact(1,-1, 0) * phi_in(k-1,i+1,j ) &
494  + fact(1, 1,-1) * phi_in(k-1,i-1,j+1) &
495  + fact(0, 1,-1) * phi_in(k ,i-1,j-1) &
496  + fact(1, 0,-1) * phi_in(k-1,i ,j-1) &
497  + fact(1,-1,-1) * phi_in(k-1,i+1,j+1) &
498  + fact(0, 0, 0) * phi_in(k ,i ,j )
499  qa_lo = fact(1, 1, 1) * phi_lo(k+1,i+1,j+1) &
500  + fact(0, 1, 1) * phi_lo(k ,i+1,j+1) &
501  + fact(1, 0, 1) * phi_lo(k+1,i ,j+1) &
502  + fact(0, 0, 1) * phi_lo(k ,i ,j+1) &
503  + fact(1,-1, 1) * phi_lo(k+1,i-1,j+1) &
504  + fact(1, 1, 0) * phi_lo(k+1,i+1,j ) &
505  + fact(0, 1, 0) * phi_lo(k ,i+1,j ) &
506  + fact(1, 0, 0) * phi_lo(k+1,i ,j ) &
507  + fact(1,-1, 0) * phi_lo(k+1,i-1,j ) &
508  + fact(1, 1,-1) * phi_lo(k+1,i+1,j-1) &
509  + fact(0, 1,-1) * phi_lo(k ,i+1,j-1) &
510  + fact(1, 0,-1) * phi_lo(k+1,i ,j-1) &
511  + fact(1,-1,-1) * phi_lo(k+1,i-1,j-1) &
512  + fact(0, 0, 0) * phi_lo(k ,i ,j )
513  qb_lo = fact(1, 1, 1) * phi_lo(k-1,i-1,j-1) &
514  + fact(0, 1, 1) * phi_lo(k ,i-1,j-1) &
515  + fact(1, 0, 1) * phi_lo(k-1,i ,j-1) &
516  + fact(0, 0, 1) * phi_lo(k ,i ,j-1) &
517  + fact(1,-1, 1) * phi_lo(k-1,i+1,j-1) &
518  + fact(1, 1, 0) * phi_lo(k-1,i-1,j ) &
519  + fact(0, 1, 0) * phi_lo(k ,i-1,j ) &
520  + fact(1, 0, 0) * phi_lo(k-1,i ,j ) &
521  + fact(1,-1, 0) * phi_lo(k-1,i+1,j ) &
522  + fact(1, 1,-1) * phi_lo(k-1,i-1,j+1) &
523  + fact(0, 1,-1) * phi_lo(k ,i-1,j-1) &
524  + fact(1, 0,-1) * phi_lo(k-1,i ,j-1) &
525  + fact(1,-1,-1) * phi_lo(k-1,i+1,j+1) &
526  + fact(0, 0, 0) * phi_lo(k ,i ,j )
527 
528  qmax = max( &
529  phi_in(k,i,j), qa_in, qb_in, &
530  phi_lo(k,i,j), qa_lo, qb_lo )
531  qmin = min( &
532  phi_in(k,i,j), qa_in, qb_in, &
533  phi_lo(k,i,j), qa_lo, qb_lo )
534  qjpls(k,i,j) = ( qmax - phi_lo(k,i,j) ) * dens(k,i,j)
535  qjmns(k,i,j) = ( phi_lo(k,i,j) - qmin ) * dens(k,i,j)
536  end do
537  end do
538  end do
539 
540  ! k = KS
541  !$omp parallel do private(i,j,rw,ru,rv,fact,qa_in,qb_in,qa_lo,qb_lo,qmax,qmin) OMP_SCHEDULE_ collapse(2)
542  do j = jjs, jje
543  do i = iis, iie
544  rw = (mflx_hi(ks,i,j,zdir) ) * rdz(ks)! 2 * rho * w / dz
545  ru = (mflx_hi(ks,i,j,xdir)+mflx_hi(ks ,i-1,j ,xdir)) * rdx(i) ! 2 * rho * u / dx
546  rv = (mflx_hi(ks,i,j,ydir)+mflx_hi(ks ,i ,j-1,ydir)) * rdy(j) ! 2 * rho * v / dy
547 
548  call get_fact_fct( fact, & ! (out)
549  rw, ru, rv ) ! (in)
550 
551  qa_in = fact(1, 1, 1) * phi_in(ks+1,i+1,j+1) &
552  + fact(0, 1, 1) * phi_in(ks ,i+1,j+1) &
553  + fact(1, 0, 1) * phi_in(ks+1,i ,j+1) &
554  + fact(0, 0, 1) * phi_in(ks ,i ,j+1) &
555  + fact(1,-1, 1) * phi_in(ks+1,i-1,j+1) &
556  + fact(1, 1, 0) * phi_in(ks+1,i+1,j ) &
557  + fact(0, 1, 0) * phi_in(ks ,i+1,j ) &
558  + fact(1, 0, 0) * phi_in(ks+1,i ,j ) &
559  + fact(1,-1, 0) * phi_in(ks+1,i-1,j ) &
560  + fact(1, 1,-1) * phi_in(ks+1,i+1,j-1) &
561  + fact(0, 1,-1) * phi_in(ks ,i+1,j-1) &
562  + fact(1, 0,-1) * phi_in(ks+1,i ,j-1) &
563  + fact(1,-1,-1) * phi_in(ks+1,i-1,j-1) &
564  + fact(0, 0, 0) * phi_in(ks ,i ,j )
565  qb_in = fact(1, 1, 1) * phi_in(ks ,i-1,j-1) &
566  + fact(0, 1, 1) * phi_in(ks ,i-1,j-1) &
567  + fact(1, 0, 1) * phi_in(ks ,i ,j-1) &
568  + fact(0, 0, 1) * phi_in(ks ,i ,j-1) &
569  + fact(1,-1, 1) * phi_in(ks ,i+1,j-1) &
570  + fact(1, 1, 0) * phi_in(ks ,i-1,j ) &
571  + fact(0, 1, 0) * phi_in(ks ,i-1,j ) &
572  + fact(1, 0, 0) * phi_in(ks ,i ,j ) &
573  + fact(1,-1, 0) * phi_in(ks ,i+1,j ) &
574  + fact(1, 1,-1) * phi_in(ks ,i-1,j+1) &
575  + fact(0, 1,-1) * phi_in(ks ,i-1,j-1) &
576  + fact(1, 0,-1) * phi_in(ks ,i ,j-1) &
577  + fact(1,-1,-1) * phi_in(ks ,i+1,j+1) &
578  + fact(0, 0, 0) * phi_in(ks ,i ,j )
579  qa_lo = fact(1, 1, 1) * phi_lo(ks+1,i+1,j+1) &
580  + fact(0, 1, 1) * phi_lo(ks ,i+1,j+1) &
581  + fact(1, 0, 1) * phi_lo(ks+1,i ,j+1) &
582  + fact(0, 0, 1) * phi_lo(ks ,i ,j+1) &
583  + fact(1,-1, 1) * phi_lo(ks+1,i-1,j+1) &
584  + fact(1, 1, 0) * phi_lo(ks+1,i+1,j ) &
585  + fact(0, 1, 0) * phi_lo(ks ,i+1,j ) &
586  + fact(1, 0, 0) * phi_lo(ks+1,i ,j ) &
587  + fact(1,-1, 0) * phi_lo(ks+1,i-1,j ) &
588  + fact(1, 1,-1) * phi_lo(ks+1,i+1,j-1) &
589  + fact(0, 1,-1) * phi_lo(ks ,i+1,j-1) &
590  + fact(1, 0,-1) * phi_lo(ks+1,i ,j-1) &
591  + fact(1,-1,-1) * phi_lo(ks+1,i-1,j-1) &
592  + fact(0, 0, 0) * phi_lo(ks ,i ,j )
593  qb_lo = fact(1, 1, 1) * phi_lo(ks ,i-1,j-1) &
594  + fact(0, 1, 1) * phi_lo(ks ,i-1,j-1) &
595  + fact(1, 0, 1) * phi_lo(ks ,i ,j-1) &
596  + fact(0, 0, 1) * phi_lo(ks ,i ,j-1) &
597  + fact(1,-1, 1) * phi_lo(ks ,i+1,j-1) &
598  + fact(1, 1, 0) * phi_lo(ks ,i-1,j ) &
599  + fact(0, 1, 0) * phi_lo(ks ,i-1,j ) &
600  + fact(1, 0, 0) * phi_lo(ks ,i ,j ) &
601  + fact(1,-1, 0) * phi_lo(ks ,i+1,j ) &
602  + fact(1, 1,-1) * phi_lo(ks ,i-1,j+1) &
603  + fact(0, 1,-1) * phi_lo(ks ,i-1,j-1) &
604  + fact(1, 0,-1) * phi_lo(ks ,i ,j-1) &
605  + fact(1,-1,-1) * phi_lo(ks ,i+1,j+1) &
606  + fact(0, 0, 0) * phi_lo(ks ,i ,j )
607 
608  qmax = max( &
609  phi_in(ks,i,j), qa_in, qb_in, &
610  phi_lo(ks,i,j), qa_lo, qb_lo )
611  qmin = min( &
612  phi_in(ks,i,j), qa_in, qb_in, &
613  phi_lo(ks,i,j), qa_lo, qb_lo )
614  qjpls(ks,i,j) = ( qmax - phi_lo(ks,i,j) ) * dens(ks,i,j)
615  qjmns(ks,i,j) = ( phi_lo(ks,i,j) - qmin ) * dens(ks,i,j)
616  end do
617  end do
618 
619  ! k = KE
620  !$omp parallel do private(i,j,rw,ru,rv,fact,qa_in,qb_in,qa_lo,qb_lo,qmax,qmin) OMP_SCHEDULE_ collapse(2)
621  do j = jjs, jje
622  do i = iis, iie
623  rw = ( mflx_hi(ke-1,i ,j ,zdir)) * rdz(ke)! 2 * rho * w / dz
624  ru = (mflx_hi(ke,i,j,xdir)+mflx_hi(ke ,i-1,j ,xdir)) * rdx(i) ! 2 * rho * u / dx
625  rv = (mflx_hi(ke,i,j,ydir)+mflx_hi(ke ,i ,j-1,ydir)) * rdy(j) ! 2 * rho * v / dy
626 
627  call get_fact_fct( fact, & ! (out)
628  rw, ru, rv ) ! (in)
629 
630  qa_in = fact(1, 1, 1) * phi_in(ke ,i+1,j+1) &
631  + fact(0, 1, 1) * phi_in(ke ,i+1,j+1) &
632  + fact(1, 0, 1) * phi_in(ke ,i ,j+1) &
633  + fact(0, 0, 1) * phi_in(ke ,i ,j+1) &
634  + fact(1,-1, 1) * phi_in(ke ,i-1,j+1) &
635  + fact(1, 1, 0) * phi_in(ke ,i+1,j ) &
636  + fact(0, 1, 0) * phi_in(ke ,i+1,j ) &
637  + fact(1, 0, 0) * phi_in(ke ,i ,j ) &
638  + fact(1,-1, 0) * phi_in(ke ,i-1,j ) &
639  + fact(1, 1,-1) * phi_in(ke ,i+1,j-1) &
640  + fact(0, 1,-1) * phi_in(ke ,i+1,j-1) &
641  + fact(1, 0,-1) * phi_in(ke ,i ,j-1) &
642  + fact(1,-1,-1) * phi_in(ke ,i-1,j-1) &
643  + fact(0, 0, 0) * phi_in(ke ,i ,j )
644  qb_in = fact(1, 1, 1) * phi_in(ke-1,i-1,j-1) &
645  + fact(0, 1, 1) * phi_in(ke ,i-1,j-1) &
646  + fact(1, 0, 1) * phi_in(ke-1,i ,j-1) &
647  + fact(0, 0, 1) * phi_in(ke ,i ,j-1) &
648  + fact(1,-1, 1) * phi_in(ke-1,i+1,j-1) &
649  + fact(1, 1, 0) * phi_in(ke-1,i-1,j ) &
650  + fact(0, 1, 0) * phi_in(ke ,i-1,j ) &
651  + fact(1, 0, 0) * phi_in(ke-1,i ,j ) &
652  + fact(1,-1, 0) * phi_in(ke-1,i+1,j ) &
653  + fact(1, 1,-1) * phi_in(ke-1,i-1,j+1) &
654  + fact(0, 1,-1) * phi_in(ke ,i-1,j-1) &
655  + fact(1, 0,-1) * phi_in(ke-1,i ,j-1) &
656  + fact(1,-1,-1) * phi_in(ke-1,i+1,j+1) &
657  + fact(0, 0, 0) * phi_in(ke ,i ,j )
658  qa_lo = fact(1, 1, 1) * phi_lo(ke ,i+1,j+1) &
659  + fact(0, 1, 1) * phi_lo(ke ,i+1,j+1) &
660  + fact(1, 0, 1) * phi_lo(ke ,i ,j+1) &
661  + fact(0, 0, 1) * phi_lo(ke ,i ,j+1) &
662  + fact(1,-1, 1) * phi_lo(ke ,i-1,j+1) &
663  + fact(1, 1, 0) * phi_lo(ke ,i+1,j ) &
664  + fact(0, 1, 0) * phi_lo(ke ,i+1,j ) &
665  + fact(1, 0, 0) * phi_lo(ke ,i ,j ) &
666  + fact(1,-1, 0) * phi_lo(ke ,i-1,j ) &
667  + fact(1, 1,-1) * phi_lo(ke ,i+1,j-1) &
668  + fact(0, 1,-1) * phi_lo(ke ,i+1,j-1) &
669  + fact(1, 0,-1) * phi_lo(ke ,i ,j-1) &
670  + fact(1,-1,-1) * phi_lo(ke ,i-1,j-1) &
671  + fact(0, 0, 0) * phi_lo(ke ,i ,j )
672  qb_lo = fact(1, 1, 1) * phi_lo(ke-1,i-1,j-1) &
673  + fact(0, 1, 1) * phi_lo(ke ,i-1,j-1) &
674  + fact(1, 0, 1) * phi_lo(ke-1,i ,j-1) &
675  + fact(0, 0, 1) * phi_lo(ke ,i ,j-1) &
676  + fact(1,-1, 1) * phi_lo(ke-1,i+1,j-1) &
677  + fact(1, 1, 0) * phi_lo(ke-1,i-1,j ) &
678  + fact(0, 1, 0) * phi_lo(ke ,i-1,j ) &
679  + fact(1, 0, 0) * phi_lo(ke-1,i ,j ) &
680  + fact(1,-1, 0) * phi_lo(ke-1,i+1,j ) &
681  + fact(1, 1,-1) * phi_lo(ke-1,i-1,j+1) &
682  + fact(0, 1,-1) * phi_lo(ke ,i-1,j-1) &
683  + fact(1, 0,-1) * phi_lo(ke-1,i ,j-1) &
684  + fact(1,-1,-1) * phi_lo(ke-1,i+1,j+1) &
685  + fact(0, 0, 0) * phi_lo(ke ,i ,j )
686 
687  qmax = max( &
688  phi_in(ke,i,j), qa_in, qb_in, &
689  phi_lo(ke,i,j), qa_lo, qb_lo )
690  qmin = min( &
691  phi_in(ke,i,j), qa_in, qb_in, &
692  phi_lo(ke,i,j), qa_lo, qb_lo )
693  qjpls(ke,i,j) = ( qmax - phi_lo(ke,i,j) ) * dens(ke,i,j)
694  qjmns(ke,i,j) = ( phi_lo(ke,i,j) - qmin ) * dens(ke,i,j)
695  end do
696  end do
697  end if
698 
699  else
700 
701  if ( twod ) then
702  i = is
703  !$omp parallel do private(j,k,qmax,qmin) OMP_SCHEDULE_
704  do j = jjs, jje
705  do k = ks+1, ke-1
706 #ifdef DEBUG
707  call check( __line__, phi_in(k ,i ,j ) )
708  call check( __line__, phi_in(k-1,i ,j ) )
709  call check( __line__, phi_in(k+1,i ,j ) )
710  call check( __line__, phi_in(k ,i ,j+1) )
711  call check( __line__, phi_in(k ,i ,j-1) )
712  call check( __line__, phi_lo(k ,i ,j ) )
713  call check( __line__, phi_lo(k-1,i ,j ) )
714  call check( __line__, phi_lo(k+1,i ,j ) )
715  call check( __line__, phi_lo(k ,i ,j+1) )
716  call check( __line__, phi_lo(k ,i ,j-1) )
717 #endif
718  qmax = max( phi_in(k ,i ,j ), &
719  phi_in(k+1,i ,j ), &
720  phi_in(k-1,i ,j ), &
721  phi_in(k ,i ,j+1), &
722  phi_in(k ,i ,j-1), &
723  phi_lo(k ,i ,j ), &
724  phi_lo(k+1,i ,j ), &
725  phi_lo(k-1,i ,j ), &
726  phi_lo(k ,i ,j+1), &
727  phi_lo(k ,i ,j-1) )
728  qmin = min( phi_in(k ,i ,j ), &
729  phi_in(k+1,i ,j ), &
730  phi_in(k-1,i ,j ), &
731  phi_in(k ,i ,j+1), &
732  phi_in(k ,i ,j-1), &
733  phi_lo(k ,i ,j ), &
734  phi_lo(k+1,i ,j ), &
735  phi_lo(k-1,i ,j ), &
736  phi_lo(k ,i ,j+1), &
737  phi_lo(k ,i ,j-1) )
738  qjpls(k,i,j) = ( qmax - phi_lo(k,i,j) ) * dens(k,i,j)
739  qjmns(k,i,j) = ( phi_lo(k,i,j) - qmin ) * dens(k,i,j)
740  enddo
741  enddo
742 #ifdef DEBUG
743  k = iundef; j = iundef
744 #endif
745  !$omp parallel do private(i,j,qmax,qmin) OMP_SCHEDULE_
746  do j = jjs, jje
747 #ifdef DEBUG
748  call check( __line__, phi_in(ks ,i ,j ) )
749  call check( __line__, phi_in(ks+1,i ,j ) )
750  call check( __line__, phi_in(ks ,i ,j+1) )
751  call check( __line__, phi_in(ks ,i ,j-1) )
752  call check( __line__, phi_lo(ks ,i ,j ) )
753  call check( __line__, phi_lo(ks+1,i ,j ) )
754  call check( __line__, phi_lo(ks ,i ,j+1) )
755  call check( __line__, phi_lo(ks ,i ,j-1) )
756  call check( __line__, phi_in(ke ,i ,j ) )
757  call check( __line__, phi_in(ke-1,i ,j ) )
758  call check( __line__, phi_in(ke ,i ,j+1) )
759  call check( __line__, phi_in(ke ,i ,j-1) )
760  call check( __line__, phi_lo(ke ,i ,j ) )
761  call check( __line__, phi_lo(ke-1,i ,j ) )
762  call check( __line__, phi_lo(ke ,i ,j+1) )
763  call check( __line__, phi_lo(ke ,i ,j-1) )
764 #endif
765  qmax = max( phi_in(ks ,i ,j ), &
766  phi_in(ks+1,i ,j ), &
767  phi_in(ks ,i ,j+1), &
768  phi_in(ks ,i ,j-1), &
769  phi_lo(ks ,i ,j ), &
770  phi_lo(ks+1,i ,j ), &
771  phi_lo(ks ,i ,j+1), &
772  phi_lo(ks ,i ,j-1) )
773  qmin = min( phi_in(ks ,i ,j ), &
774  phi_in(ks+1,i ,j ), &
775  phi_in(ks ,i ,j+1), &
776  phi_in(ks ,i ,j-1), &
777  phi_lo(ks ,i ,j ), &
778  phi_lo(ks+1,i ,j ), &
779  phi_lo(ks ,i ,j+1), &
780  phi_lo(ks ,i ,j-1) )
781  qjmns(ks,i,j) = ( phi_lo(ks,i,j) - qmin ) * dens(ks,i,j)
782  qjpls(ks,i,j) = ( qmax - phi_lo(ks,i,j) ) * dens(ks,i,j)
783 
784  qmax = max( phi_in(ke ,i ,j ), &
785  phi_in(ke-1,i ,j ), &
786  phi_in(ke ,i ,j+1), &
787  phi_in(ke ,i ,j-1), &
788  phi_lo(ke ,i ,j ), &
789  phi_lo(ke-1,i ,j ), &
790  phi_lo(ke ,i ,j+1), &
791  phi_lo(ke ,i ,j-1) )
792  qmin = min( phi_in(ke ,i ,j ), &
793  phi_in(ke-1,i ,j ), &
794  phi_in(ke ,i ,j+1), &
795  phi_in(ke ,i ,j-1), &
796  phi_lo(ke ,i ,j ), &
797  phi_lo(ke-1,i ,j ), &
798  phi_lo(ke ,i ,j+1), &
799  phi_lo(ke ,i ,j-1) )
800  qjpls(ke,i,j) = ( qmax - phi_lo(ke,i,j) ) * dens(ke,i,j)
801  qjmns(ke,i,j) = ( phi_lo(ke,i,j) - qmin ) * dens(ke,i,j)
802  enddo
803 #ifdef DEBUG
804  k = iundef; i = iundef; j = iundef
805 #endif
806  else
807  !$omp parallel do private(i,j,k,qmax,qmin) OMP_SCHEDULE_ collapse(2)
808  do j = jjs, jje
809  do i = iis, iie
810  do k = ks+1, ke-1
811 #ifdef DEBUG
812  call check( __line__, phi_in(k ,i ,j ) )
813  call check( __line__, phi_in(k-1,i ,j ) )
814  call check( __line__, phi_in(k+1,i ,j ) )
815  call check( __line__, phi_in(k ,i-1,j ) )
816  call check( __line__, phi_in(k ,i+1,j ) )
817  call check( __line__, phi_in(k ,i ,j+1) )
818  call check( __line__, phi_in(k ,i ,j-1) )
819  call check( __line__, phi_lo(k ,i ,j ) )
820  call check( __line__, phi_lo(k-1,i ,j ) )
821  call check( __line__, phi_lo(k+1,i ,j ) )
822  call check( __line__, phi_lo(k ,i-1,j ) )
823  call check( __line__, phi_lo(k ,i+1,j ) )
824  call check( __line__, phi_lo(k ,i ,j+1) )
825  call check( __line__, phi_lo(k ,i ,j-1) )
826 #endif
827  qmax = max( phi_in(k ,i ,j ), &
828  phi_in(k+1,i ,j ), &
829  phi_in(k-1,i ,j ), &
830  phi_in(k ,i+1,j ), &
831  phi_in(k ,i-1,j ), &
832  phi_in(k ,i ,j+1), &
833  phi_in(k ,i ,j-1), &
834  phi_lo(k ,i ,j ), &
835  phi_lo(k+1,i ,j ), &
836  phi_lo(k-1,i ,j ), &
837  phi_lo(k ,i+1,j ), &
838  phi_lo(k ,i-1,j ), &
839  phi_lo(k ,i ,j+1), &
840  phi_lo(k ,i ,j-1) )
841  qmin = min( phi_in(k ,i ,j ), &
842  phi_in(k+1,i ,j ), &
843  phi_in(k-1,i ,j ), &
844  phi_in(k ,i-1,j ), &
845  phi_in(k ,i+1,j ), &
846  phi_in(k ,i ,j+1), &
847  phi_in(k ,i ,j-1), &
848  phi_lo(k ,i ,j ), &
849  phi_lo(k+1,i ,j ), &
850  phi_lo(k-1,i ,j ), &
851  phi_lo(k ,i-1,j ), &
852  phi_lo(k ,i+1,j ), &
853  phi_lo(k ,i ,j+1), &
854  phi_lo(k ,i ,j-1) )
855  qjpls(k,i,j) = ( qmax - phi_lo(k,i,j) ) * dens(k,i,j)
856  qjmns(k,i,j) = ( phi_lo(k,i,j) - qmin ) * dens(k,i,j)
857  enddo
858  enddo
859  enddo
860 #ifdef DEBUG
861  k = iundef; i = iundef; j = iundef
862 #endif
863  !$omp parallel do private(i,j,qmax,qmin) OMP_SCHEDULE_ collapse(2)
864  do j = jjs, jje
865  do i = iis, iie
866 #ifdef DEBUG
867  call check( __line__, phi_in(ks ,i ,j ) )
868  call check( __line__, phi_in(ks+1,i ,j ) )
869  call check( __line__, phi_in(ks ,i-1,j ) )
870  call check( __line__, phi_in(ks ,i+1,j ) )
871  call check( __line__, phi_in(ks ,i ,j+1) )
872  call check( __line__, phi_in(ks ,i ,j-1) )
873  call check( __line__, phi_lo(ks ,i ,j ) )
874  call check( __line__, phi_lo(ks+1,i ,j ) )
875  call check( __line__, phi_lo(ks ,i-1,j ) )
876  call check( __line__, phi_lo(ks ,i+1,j ) )
877  call check( __line__, phi_lo(ks ,i ,j+1) )
878  call check( __line__, phi_lo(ks ,i ,j-1) )
879  call check( __line__, phi_in(ke ,i ,j ) )
880  call check( __line__, phi_in(ke-1,i ,j ) )
881  call check( __line__, phi_in(ke ,i-1,j ) )
882  call check( __line__, phi_in(ke ,i+1,j ) )
883  call check( __line__, phi_in(ke ,i ,j+1) )
884  call check( __line__, phi_in(ke ,i ,j-1) )
885  call check( __line__, phi_lo(ke ,i ,j ) )
886  call check( __line__, phi_lo(ke-1,i ,j ) )
887  call check( __line__, phi_lo(ke ,i-1,j ) )
888  call check( __line__, phi_lo(ke ,i+1,j ) )
889  call check( __line__, phi_lo(ke ,i ,j+1) )
890  call check( __line__, phi_lo(ke ,i ,j-1) )
891 #endif
892  qmax = max( phi_in(ks ,i ,j ), &
893  phi_in(ks+1,i ,j ), &
894  phi_in(ks ,i+1,j ), &
895  phi_in(ks ,i-1,j ), &
896  phi_in(ks ,i ,j+1), &
897  phi_in(ks ,i ,j-1), &
898  phi_lo(ks ,i ,j ), &
899  phi_lo(ks+1,i ,j ), &
900  phi_lo(ks ,i+1,j ), &
901  phi_lo(ks ,i-1,j ), &
902  phi_lo(ks ,i ,j+1), &
903  phi_lo(ks ,i ,j-1) )
904  qmin = min( phi_in(ks ,i ,j ), &
905  phi_in(ks+1,i ,j ), &
906  phi_in(ks ,i+1,j ), &
907  phi_in(ks ,i-1,j ), &
908  phi_in(ks ,i ,j+1), &
909  phi_in(ks ,i ,j-1), &
910  phi_lo(ks ,i ,j ), &
911  phi_lo(ks+1,i ,j ), &
912  phi_lo(ks ,i+1,j ), &
913  phi_lo(ks ,i-1,j ), &
914  phi_lo(ks ,i ,j+1), &
915  phi_lo(ks ,i ,j-1) )
916  qjmns(ks,i,j) = ( phi_lo(ks,i,j) - qmin ) * dens(ks,i,j)
917  qjpls(ks,i,j) = ( qmax - phi_lo(ks,i,j) ) * dens(ks,i,j)
918 
919  qmax = max( phi_in(ke ,i ,j ), &
920  phi_in(ke-1,i ,j ), &
921  phi_in(ke ,i+1,j ), &
922  phi_in(ke ,i-1,j ), &
923  phi_in(ke ,i ,j+1), &
924  phi_in(ke ,i ,j-1), &
925  phi_lo(ke ,i ,j ), &
926  phi_lo(ke-1,i ,j ), &
927  phi_lo(ke ,i+1,j ), &
928  phi_lo(ke ,i-1,j ), &
929  phi_lo(ke ,i ,j+1), &
930  phi_lo(ke ,i ,j-1) )
931  qmin = min( phi_in(ke ,i ,j ), &
932  phi_in(ke-1,i ,j ), &
933  phi_in(ke ,i-1,j ), &
934  phi_in(ke ,i+1,j ), &
935  phi_in(ke ,i ,j+1), &
936  phi_in(ke ,i ,j-1), &
937  phi_lo(ke ,i ,j ), &
938  phi_lo(ke-1,i ,j ), &
939  phi_lo(ke ,i-1,j ), &
940  phi_lo(ke ,i+1,j ), &
941  phi_lo(ke ,i ,j+1), &
942  phi_lo(ke ,i ,j-1) )
943  qjpls(ke,i,j) = ( qmax - phi_lo(ke,i,j) ) * dens(ke,i,j)
944  qjmns(ke,i,j) = ( phi_lo(ke,i,j) - qmin ) * dens(ke,i,j)
945  enddo
946  enddo
947 #ifdef DEBUG
948  k = iundef; i = iundef; j = iundef
949 #endif
950  end if
951 
952  end if
953 
954  !--- incoming flux limitation factor (0-1)
955  !$omp parallel do private(i,j,k,zerosw) OMP_SCHEDULE_ collapse(2)
956  do j = jjs, jje
957  do i = iis, iie
958  do k = ks, ke
959 #ifdef DEBUG
960  call check( __line__, pjpls(k,i,j) )
961  call check( __line__, qjpls(k,i,j) )
962 #endif
963  ! if pjpls == 0, zerosw = 1 and rjpls = 0
964  zerosw = 0.5_rp - sign( 0.5_rp, pjpls(k,i,j)-epsilon )
965  rjpls(k,i,j) = min( 1.0_rp, qjpls(k,i,j) * ( 1.0_rp-zerosw ) / ( pjpls(k,i,j)-zerosw ) )
966  enddo
967  enddo
968  enddo
969 #ifdef DEBUG
970  k = iundef; i = iundef; j = iundef
971 #endif
972 
973  !--- outgoing flux limitation factor (0-1)
974  !$omp parallel do private(i,j,k,zerosw) OMP_SCHEDULE_ collapse(2)
975  do j = jjs, jje
976  do i = iis, iie
977  do k = ks, ke
978 #ifdef DEBUG
979  call check( __line__, pjmns(k,i,j) )
980  call check( __line__, qjmns(k,i,j) )
981 #endif
982  ! if pjmns == 0, zerosw = 1 and rjmns = 0
983  zerosw = 0.5_rp - sign( 0.5_rp, pjmns(k,i,j)-epsilon )
984  rjmns(k,i,j) = min( 1.0_rp, qjmns(k,i,j) * ( 1.0_rp-zerosw ) / ( pjmns(k,i,j)-zerosw ) )
985  enddo
986  enddo
987  enddo
988 #ifdef DEBUG
989  k = iundef; i = iundef; j = iundef
990 #endif
991 
992  enddo
993  enddo
994 
995  call comm_vars8( rjpls(:,:,:), 1 )
996  call comm_vars8( rjmns(:,:,:), 2 )
997  call comm_wait ( rjpls(:,:,:), 1 )
998  call comm_wait ( rjmns(:,:,:), 2 )
999 
1000  do jjs = js, je, jblock
1001  jje = jjs+jblock-1
1002  do iis = is, ie, iblock
1003  iie = iis+iblock-1
1004 
1005  !--- update high order flux with antidiffusive flux
1006  !$omp parallel do private(i,j,k,dirsw) OMP_SCHEDULE_ collapse(2)
1007  do j = jjs, jje
1008  do i = iis, iie
1009  do k = ks , ke-1
1010 #ifdef DEBUG
1011  call check( __line__, qflx_anti(k,i,j,zdir) )
1012  call check( __line__, rjpls(k ,i,j) )
1013  call check( __line__, rjpls(k+1,i,j) )
1014  call check( __line__, rjmns(k ,i,j) )
1015  call check( __line__, rjmns(k+1,i,j) )
1016 #endif
1017  ! if qflx_anti > 0, dirsw = 1
1018  dirsw = 0.5_rp + sign( 0.5_rp, qflx_anti(k,i,j,zdir) )
1019  qflx_anti(k,i,j,zdir) = qflx_anti(k,i,j,zdir) &
1020  * ( 1.0_rp &
1021  - min( rjpls(k+1,i,j),rjmns(k ,i,j) ) * ( dirsw ) &
1022  - min( rjpls(k ,i,j),rjmns(k+1,i,j) ) * ( 1.0_rp - dirsw ) )
1023  enddo
1024  enddo
1025  enddo
1026 #ifdef DEBUG
1027  k = iundef; i = iundef; j = iundef
1028 #endif
1029  !$omp parallel do private(i,j) OMP_SCHEDULE_ collapse(2)
1030  do j = jjs, jje
1031  do i = iis, iie
1032 #ifdef DEBUG
1033  call check( __line__, qflx_anti(ke,i,j,zdir) )
1034  call check( __line__, rjpls(ke ,i,j) )
1035  call check( __line__, rjmns(ke ,i,j) )
1036 #endif
1037  qflx_anti(ks-1,i,j,zdir) = 0.0_rp ! bottom boundary
1038  qflx_anti(ke ,i,j,zdir) = 0.0_rp ! top boundary
1039  enddo
1040  enddo
1041 #ifdef DEBUG
1042  k = iundef; i = iundef; j = iundef
1043 #endif
1044 
1045  if ( .not. twod ) then
1046  if ( iis == is ) then
1047  ijs = iis-1
1048  else
1049  ijs = iis
1050  end if
1051 
1052  !$omp parallel do private(i,j,k,dirsw) OMP_SCHEDULE_ collapse(2)
1053  do j = jjs, jje
1054  do i = ijs, iie
1055  do k = ks, ke
1056 #ifdef DEBUG
1057  call check( __line__, qflx_anti(k,i,j,xdir) )
1058  call check( __line__, rjpls(k,i ,j) )
1059  call check( __line__, rjpls(k,i+1,j) )
1060  call check( __line__, rjmns(k,i ,j) )
1061  call check( __line__, rjmns(k,i+1,j) )
1062 #endif
1063  ! if qflx_anti > 0, dirsw = 1
1064  dirsw = 0.5_rp + sign( 0.5_rp, qflx_anti(k,i,j,xdir) )
1065  qflx_anti(k,i,j,xdir) = qflx_anti(k,i,j,xdir) &
1066  * ( 1.0_rp &
1067  - min( rjpls(k,i+1,j),rjmns(k,i ,j) ) * ( dirsw ) &
1068  - min( rjpls(k,i ,j),rjmns(k,i+1,j) ) * ( 1.0_rp - dirsw ) )
1069  enddo
1070  enddo
1071  enddo
1072 #ifdef DEBUG
1073  k = iundef; i = iundef; j = iundef
1074 #endif
1075  end if
1076 
1077  if ( jjs == js ) then
1078  ijs = jjs-1
1079  else
1080  ijs = jjs
1081  end if
1082  !$omp parallel do private(i,j,k,dirsw) OMP_SCHEDULE_ collapse(2)
1083  do j = ijs, jje
1084  do i = iis, iie
1085  do k = ks, ke
1086 #ifdef DEBUG
1087  call check( __line__, qflx_anti(k,i,j,ydir) )
1088  call check( __line__, rjpls(k,i,j+1) )
1089  call check( __line__, rjpls(k,i,j ) )
1090  call check( __line__, rjmns(k,i,j ) )
1091  call check( __line__, rjmns(k,i,j+1) )
1092 #endif
1093  ! if qflx_anti > 0, dirsw = 1
1094  dirsw = 0.5_rp + sign( 0.5_rp, qflx_anti(k,i,j,ydir) )
1095  qflx_anti(k,i,j,ydir) = qflx_anti(k,i,j,ydir) &
1096  * ( 1.0_rp &
1097  - min( rjpls(k,i,j+1),rjmns(k,i,j ) ) * ( dirsw ) &
1098  - min( rjpls(k,i,j ),rjmns(k,i,j+1) ) * ( 1.0_rp - dirsw ) )
1099  enddo
1100  enddo
1101  enddo
1102 #ifdef DEBUG
1103  k = iundef; i = iundef; j = iundef
1104 #endif
1105 
1106  enddo
1107  enddo
1108 
1109  return
1110  end subroutine atmos_dyn_fvm_fct
1111 
1112  !-----------------------------------------------------------------------------
1113  ! private procedure
1114  ! get factor for FCT
1115  subroutine get_fact_fct( &
1116  fact, &
1117  rw, ru, rv )
1118  use scale_const, only: &
1119  epsilon => const_eps
1120  implicit none
1121  real(rp), intent(out) :: fact(0:1,-1:1,-1:1)
1122  real(rp), intent(in) :: rw, ru, rv
1123 
1124  real(rp) :: sign_uv, sign_uw, sign_vw ! uv>=0, uw>=0, vw>=0
1125  real(rp) :: ugev, ugew, vgew ! u>=v, u>=w, u>=w
1126  real(rp) :: umax, vmax, wmax
1127  real(rp) :: vu, wu, uv, wv, uw, vw ! |v/u|, |w/u|, ....
1128  real(rp) :: uzero, vzero, wzero
1129  !---------------------------------------------------------------------------
1130 
1131  ugev = sign(0.5_rp, abs(ru)-abs(rv)) + 0.5_rp ! u >= v
1132  ugew = sign(0.5_rp, abs(ru)-abs(rw)) + 0.5_rp ! u >= w
1133  vgew = sign(0.5_rp, abs(rv)-abs(rw)) + 0.5_rp ! v >= w
1134 
1135  uzero = sign(0.5_rp,abs(ru)-epsilon) - 0.5_rp
1136  vzero = sign(0.5_rp,abs(rv)-epsilon) - 0.5_rp
1137  wzero = sign(0.5_rp,abs(rw)-epsilon) - 0.5_rp
1138 
1139  sign_uv = sign(0.5_rp, ru*rv) + 0.5_rp ! uv >= 0
1140  sign_uw = sign(0.5_rp, ru*rw) + 0.5_rp ! uw >= 0
1141  sign_vw = sign(0.5_rp, rv*rw) + 0.5_rp ! vw >= 0
1142 
1143  wu = abs( rw / ( ru+uzero ) * ( 1.0_rp+uzero ) )
1144  vu = abs( rv / ( ru+uzero ) * ( 1.0_rp+uzero ) )
1145  uv = abs( ru / ( rv+vzero ) * ( 1.0_rp+vzero ) )
1146  wv = abs( rw / ( rv+vzero ) * ( 1.0_rp+vzero ) )
1147  uw = abs( ru / ( rw+wzero ) * ( 1.0_rp+wzero ) )
1148  vw = abs( rv / ( rw+wzero ) * ( 1.0_rp+wzero ) )
1149 
1150  umax = ugev * ugew * ( 1.0_rp+uzero ) ! u == max(u,v,w)
1151  vmax = (1.0_rp-ugev) * vgew ! v == max(u,v,w)
1152  wmax = 1.0_rp - ugev * ugew - vmax ! w == max(u,v,w)
1153 
1154  fact(0, 0, 0) = - ugev * ugew * uzero ! 1.0 if max(u,v,w) < epsilon
1155 
1156  fact(1, 0, 0) = wmax * (1.0_rp-uw) * (1.0_rp-vw)
1157  fact(0, 1, 0) = umax * (1.0_rp-vu) * (1.0_rp-wu)
1158  fact(0, 0, 1) = vmax * (1.0_rp-uv) * (1.0_rp-wv)
1159 
1160  fact(1, 1, 1) = sign_uv * sign_uw * ( umax * vu*wu + vmax * uv*wv + wmax * uw*vw )
1161  fact(1,-1, 1) = (1.0_rp-sign_uv) * (1.0_rp-sign_uw) * ( umax * vu*wu + vmax * uv*wv + wmax * uw*vw )
1162  fact(1, 1,-1) = (1.0_rp-sign_uv) * sign_uw * ( umax * vu*wu + vmax * uv*wv + wmax * uw*vw )
1163  fact(1,-1,-1) = sign_uv * (1.0_rp-sign_uw) * ( umax * vu*wu + vmax * uv*wv + wmax * uw*vw )
1164 
1165  fact(1, 1, 0) = sign_uw * (1.0_rp-vmax) * ( ugew * wu * (1.0_rp-vu) + (1.0_rp-ugew) * uw * (1.0_rp-vw) )
1166  fact(1,-1, 0) = (1.0_rp-sign_uw) * (1.0_rp-vmax) * ( ugew * wu * (1.0_rp-vu) + (1.0_rp-ugew) * uw * (1.0_rp-vw) )
1167  fact(1, 0, 1) = sign_vw * (1.0_rp-umax) * ( vgew * wv * (1.0_rp-uv) + (1.0_rp-vgew) * vw * (1.0_rp-uw) )
1168  fact(1, 0,-1) = (1.0_rp-sign_vw) * (1.0_rp-umax) * ( vgew * wv * (1.0_rp-uv) + (1.0_rp-vgew) * vw * (1.0_rp-uw) )
1169  fact(0, 1, 1) = sign_uv * (1.0_rp-wmax) * ( ugev * vu * (1.0_rp-wu) + (1.0_rp-ugev) * uv * (1.0_rp-wv) )
1170  fact(0, 1,-1) = (1.0_rp-sign_uv) * (1.0_rp-wmax) * ( ugev * vu * (1.0_rp-wu) + (1.0_rp-ugev) * uv * (1.0_rp-wv) )
1171 
1172  return
1173  end subroutine get_fact_fct
1174 
1175 end module scale_atmos_dyn_fvm_fct
scale_atmos_grid_cartesc_index::ke
integer, public ke
end point of inner domain: z, local
Definition: scale_atmos_grid_cartesC_index.F90:52
scale_atmos_dyn_fvm_fct::atmos_dyn_fvm_fct
subroutine, public atmos_dyn_fvm_fct(qflx_anti, phi_in, DENS0, DENS, qflx_hi, qflx_lo, mflx_hi, rdz, rdx, rdy, GSQRT, MAPF, TwoD, dt, flag_vect)
Setup.
Definition: scale_atmos_dyn_fvm_fct.F90:72
scale_index
module Index
Definition: scale_index.F90:11
scale_const::const_undef2
integer, parameter, public const_undef2
undefined value (INT2)
Definition: scale_const.F90:40
scale_precision
module PRECISION
Definition: scale_precision.F90:14
scale_atmos_grid_cartesc_index::ka
integer, public ka
Definition: scale_atmos_grid_cartesC_index.F90:47
scale_atmos_grid_cartesc_index::jblock
integer, public jblock
block size for cache blocking: y
Definition: scale_atmos_grid_cartesC_index.F90:41
scale_const::const_eps
real(rp), public const_eps
small number
Definition: scale_const.F90:35
scale_atmos_grid_cartesc_index::iblock
integer, public iblock
block size for cache blocking: x
Definition: scale_atmos_grid_cartesC_index.F90:40
scale_precision::rp
integer, parameter, public rp
Definition: scale_precision.F90:41
scale_atmos_grid_cartesc_index::ie
integer, public ie
end point of inner domain: x, local
Definition: scale_atmos_grid_cartesC_index.F90:54
scale_io
module STDIO
Definition: scale_io.F90:10
scale_tracer::k
real(rp), public k
Definition: scale_tracer.F90:45
scale_atmos_grid_cartesc_index
module atmosphere / grid / cartesC index
Definition: scale_atmos_grid_cartesC_index.F90:12
scale_const
module CONSTANT
Definition: scale_const.F90:11
scale_atmos_grid_cartesc_index::ia
integer, public ia
Definition: scale_atmos_grid_cartesC_index.F90:48
scale_debug::check
subroutine, public check(current_line, v)
Undefined value checker.
Definition: scale_debug.F90:59
scale_atmos_grid_cartesc_index::zdir
integer, parameter, public zdir
Definition: scale_atmos_grid_cartesC_index.F90:32
scale_prof
module profiler
Definition: scale_prof.F90:11
scale_atmos_grid_cartesc_index::is
integer, public is
start point of inner domain: x, local
Definition: scale_atmos_grid_cartesC_index.F90:53
scale_atmos_grid_cartesc_index::ja
integer, public ja
Definition: scale_atmos_grid_cartesC_index.F90:49
scale_tracer
module TRACER
Definition: scale_tracer.F90:12
scale_atmos_grid_cartesc_index::xdir
integer, parameter, public xdir
Definition: scale_atmos_grid_cartesC_index.F90:33
scale_atmos_grid_cartesc_index::ks
integer, public ks
start point of inner domain: z, local
Definition: scale_atmos_grid_cartesC_index.F90:51
scale_debug
module DEBUG
Definition: scale_debug.F90:11
scale_comm_cartesc
module COMMUNICATION
Definition: scale_comm_cartesC.F90:11
scale_atmos_grid_cartesc_index::js
integer, public js
start point of inner domain: y, local
Definition: scale_atmos_grid_cartesC_index.F90:55
scale_atmos_grid_cartesc_index::ydir
integer, parameter, public ydir
Definition: scale_atmos_grid_cartesC_index.F90:34
scale_const::const_undef
real(rp), public const_undef
Definition: scale_const.F90:43
scale_atmos_dyn_fvm_fct
module Atmosphere / Dynamics common
Definition: scale_atmos_dyn_fvm_fct.F90:12
scale_atmos_grid_cartesc_index::je
integer, public je
end point of inner domain: y, local
Definition: scale_atmos_grid_cartesC_index.F90:56
Generated by   doxygen 1.8.17