SCALE-RM
scale_atmos_dyn_common.F90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
13 !-------------------------------------------------------------------------------
14 #include "inc_openmp.h"
15 module scale_atmos_dyn_common
16  !-----------------------------------------------------------------------------
17  !
18  !++ used modules
19  !
20  use scale_precision
21  use scale_stdio
22  use scale_prof
23  use scale_grid_index
24  use scale_index
25  use scale_tracer
26 
27 #ifdef DEBUG
28  use scale_debug, only: &
29  check
30  use scale_const, only: &
31  undef => const_undef, &
32  iundef => const_undef2
33 #endif
34  !-----------------------------------------------------------------------------
35  implicit none
36  private
37  !-----------------------------------------------------------------------------
38  !
39  !++ Public procedure
40  !
41  public :: atmos_dyn_filter_setup
42  public :: atmos_dyn_wdamp_setup
43  public :: atmos_dyn_numfilter_coef
44  public :: atmos_dyn_numfilter_coef_q
45  public :: atmos_dyn_filter_tend
46  public :: atmos_dyn_fct
47  public :: atmos_dyn_copy_boundary
48  public :: atmos_dyn_copy_boundary_tracer
49  public :: atmos_dyn_divergence
50 
51  !-----------------------------------------------------------------------------
52  !
53  !++ Public parameters & variables
54  !
55  !-----------------------------------------------------------------------------
56  !
57  !++ Private procedure
58  !
59  private :: get_fact_fct
60 
61  !-----------------------------------------------------------------------------
62  !
63  !++ Private parameters & variables
64  !
65  !-----------------------------------------------------------------------------
66  real(RP), allocatable :: CNZ3(:,:,:)
67  real(RP), allocatable :: CNX3(:,:,:)
68  real(RP), allocatable :: CNY3(:,:,:)
69  real(RP), allocatable :: CNZ4(:,:,:)
70  real(RP), allocatable :: CNX4(:,:,:)
71  real(RP), allocatable :: CNY4(:,:,:)
72 
73  integer :: I_COMM_DENS_Z = 1
74  integer :: I_COMM_DENS_X = 2
75  integer :: I_COMM_DENS_Y = 3
76  integer :: I_COMM_MOMZ_Z = 4
77  integer :: I_COMM_MOMZ_X = 5
78  integer :: I_COMM_MOMZ_Y = 6
79  integer :: I_COMM_MOMX_Z = 7
80  integer :: I_COMM_MOMX_X = 8
81  integer :: I_COMM_MOMX_Y = 9
82  integer :: I_COMM_MOMY_Z = 10
83  integer :: I_COMM_MOMY_X = 11
84  integer :: I_COMM_MOMY_Y = 12
85  integer :: I_COMM_RHOT_Z = 13
86  integer :: I_COMM_RHOT_X = 14
87  integer :: I_COMM_RHOT_Y = 15
88  integer :: I_COMM_QTRC_Z = 1
89  integer :: I_COMM_QTRC_X = 2
90  integer :: I_COMM_QTRC_Y = 3
91 
92 contains
93  !-----------------------------------------------------------------------------
95  subroutine atmos_dyn_filter_setup( &
96  num_diff, num_diff_q, &
97  CDZ, CDX, CDY, FDZ, FDX, FDY )
98  use scale_process, only: &
99  prc_mpistop
100  use scale_comm, only: &
101  comm_vars8_init
102  implicit none
103  real(RP), intent(inout) :: num_diff(ka,ia,ja,5,3)
104  real(RP), intent(inout) :: num_diff_q(ka,ia,ja,3)
105  real(RP), intent(in) :: cdz(ka)
106  real(RP), intent(in) :: cdx(ia)
107  real(RP), intent(in) :: cdy(ja)
108  real(RP), intent(in) :: fdz(ka-1)
109  real(RP), intent(in) :: fdx(ia-1)
110  real(RP), intent(in) :: fdy(ja-1)
111 
112  integer :: k, i, j
113  !---------------------------------------------------------------------------
114 
115  if ( ihalo < 2 .or. jhalo < 2 .or. khalo < 2 ) then
116  write(*,*) 'xxx number of HALO must be at least 2 for numrical filter'
117  call prc_mpistop
118  end if
119 
120  ! allocation
121  allocate( cnz3(3,ka,2) )
122  allocate( cnx3(3,ia,2) )
123  allocate( cny3(3,ja,2) )
124  allocate( cnz4(5,ka,2) )
125  allocate( cnx4(5,ia,2) )
126  allocate( cny4(5,ja,2) )
127 
128 
129  call comm_vars8_init( 'num_diff_DENS_Z', num_diff(:,:,:,i_dens,zdir), i_comm_dens_z )
130  call comm_vars8_init( 'num_diff_DENS_X', num_diff(:,:,:,i_dens,xdir), i_comm_dens_x )
131  call comm_vars8_init( 'num_diff_DENS_Y', num_diff(:,:,:,i_dens,ydir), i_comm_dens_y )
132  call comm_vars8_init( 'num_diff_MOMZ_Z', num_diff(:,:,:,i_momz,zdir), i_comm_momz_z )
133  call comm_vars8_init( 'num_diff_MOMZ_X', num_diff(:,:,:,i_momz,xdir), i_comm_momz_x )
134  call comm_vars8_init( 'num_diff_MOMZ_Y', num_diff(:,:,:,i_momz,ydir), i_comm_momz_y )
135  call comm_vars8_init( 'num_diff_MOMX_Z', num_diff(:,:,:,i_momx,zdir), i_comm_momx_z )
136  call comm_vars8_init( 'num_diff_MOMX_X', num_diff(:,:,:,i_momx,xdir), i_comm_momx_x )
137  call comm_vars8_init( 'num_diff_MOMX_Y', num_diff(:,:,:,i_momx,ydir), i_comm_momx_y )
138  call comm_vars8_init( 'num_diff_MOMY_Z', num_diff(:,:,:,i_momy,zdir), i_comm_momy_z )
139  call comm_vars8_init( 'num_diff_MOMY_X', num_diff(:,:,:,i_momy,xdir), i_comm_momy_x )
140  call comm_vars8_init( 'num_diff_MOMY_Y', num_diff(:,:,:,i_momy,ydir), i_comm_momy_y )
141  call comm_vars8_init( 'num_diff_RHOT_Z', num_diff(:,:,:,i_rhot,zdir), i_comm_rhot_z )
142  call comm_vars8_init( 'num_diff_RHOT_X', num_diff(:,:,:,i_rhot,xdir), i_comm_rhot_x )
143  call comm_vars8_init( 'num_diff_RHOT_Y', num_diff(:,:,:,i_rhot,ydir), i_comm_rhot_y )
144 
145  call comm_vars8_init( 'num_diff_QTRC_Z', num_diff_q(:,:,:,zdir), i_comm_qtrc_z )
146  call comm_vars8_init( 'num_diff_QTRC_X', num_diff_q(:,:,:,xdir), i_comm_qtrc_x )
147  call comm_vars8_init( 'num_diff_QTRC_Y', num_diff_q(:,:,:,ydir), i_comm_qtrc_y )
148 
149 #ifdef DEBUG
150  cnz3(:,:,:) = undef
151  cnx3(:,:,:) = undef
152  cny3(:,:,:) = undef
153  cnz4(:,:,:) = undef
154  cnx4(:,:,:) = undef
155  cny4(:,:,:) = undef
156 #endif
157 
158  ! z direction
159  do k = ks+1, ke-1
160  cnz3(1,k,1) = 1.0_rp / ( fdz(k ) * cdz(k ) * fdz(k-1) )
161  cnz3(2,k,1) = 1.0_rp / ( fdz(k ) * cdz(k ) * fdz(k-1) ) &
162  + 1.0_rp / ( fdz(k-1) * cdz(k ) * fdz(k-1) ) &
163  + 1.0_rp / ( fdz(k-1) * cdz(k-1) * fdz(k-1) )
164  enddo
165  do k = ks+2, ke
166  cnz3(3,k,1) = 1.0_rp / ( fdz(k-1) * cdz(k ) * fdz(k-1) ) &
167  + 1.0_rp / ( fdz(k-1) * cdz(k-1) * fdz(k-1) ) &
168  + 1.0_rp / ( fdz(k-1) * cdz(k-1) * fdz(k-2) )
169  enddo
170  cnz3(1,ks-1,1) = 1.0_rp / ( fdz(ks ) * cdz(ks+1) * fdz(ks ) )
171  cnz3(1,ks ,1) = 1.0_rp / ( fdz(ks ) * cdz(ks ) * fdz(ks ) )
172  cnz3(2,ks ,1) = 1.0_rp / ( fdz(ks ) * cdz(ks ) * fdz(ks ) ) &
173  + 1.0_rp / ( fdz(ks ) * cdz(ks ) * fdz(ks ) ) &
174  + 1.0_rp / ( fdz(ks ) * cdz(ks+1) * fdz(ks ) )
175  cnz3(3,ks ,1) = 1.0_rp / ( fdz(ks ) * cdz(ks ) * fdz(ks ) ) &
176  + 1.0_rp / ( fdz(ks ) * cdz(ks+1) * fdz(ks ) ) &
177  + 1.0_rp / ( fdz(ks ) * cdz(ks+1) * fdz(ks+1) )
178  cnz3(3,ks+1,1) = 1.0_rp / ( fdz(ks ) * cdz(ks+1) * fdz(ks ) ) &
179  + 1.0_rp / ( fdz(ks ) * cdz(ks ) * fdz(ks ) ) &
180  + 1.0_rp / ( fdz(ks ) * cdz(ks ) * fdz(ks ) )
181  cnz3(1,ke ,1) = 1.0_rp / ( fdz(ke-1) * cdz(ke ) * fdz(ke-1) )
182  cnz3(2,ke ,1) = 1.0_rp / ( fdz(ke-1) * cdz(ke ) * fdz(ke-1) ) &
183  + 1.0_rp / ( fdz(ke-1) * cdz(ke ) * fdz(ke-1) ) &
184  + 1.0_rp / ( fdz(ke-1) * cdz(ke-1) * fdz(ke-1) )
185  cnz3(1,ke+1,1) = 1.0_rp / ( fdz(ke-2) * cdz(ke-1) * fdz(ke-1) )
186  cnz3(2,ke+1,1) = 1.0_rp / ( fdz(ke-2) * cdz(ke-1) * fdz(ke-1) ) &
187  + 1.0_rp / ( fdz(ke-1) * cdz(ke-1) * fdz(ke-1) ) &
188  + 1.0_rp / ( fdz(ke-1) * cdz(ke ) * fdz(ke-1) )
189  cnz3(3,ke+1,1) = 1.0_rp / ( fdz(ke-1) * cdz(ke+1) * fdz(ke-1) ) &
190  + 1.0_rp / ( fdz(ke-1) * cdz(ke ) * fdz(ke-1) ) &
191  + 1.0_rp / ( fdz(ke-1) * cdz(ke ) * fdz(ke-1) )
192 
193  do k = ks, ke-1
194  cnz4(1,k,1) = ( cnz3(1,k+1,1) ) / cdz(k)
195  cnz4(2,k,1) = ( cnz3(2,k+1,1) + cnz3(1,k,1) ) / cdz(k)
196  cnz4(3,k,1) = ( cnz3(3,k+1,1) + cnz3(2,k,1) ) / cdz(k)
197  cnz4(4,k,1) = ( cnz3(1,k ,1) + cnz3(3,k,1) ) / cdz(k)
198  cnz4(5,k,1) = ( cnz3(1,k-1,1) ) / cdz(k)
199  enddo
200 
201  do k = ks+1, ke-1
202  cnz3(1,k,2) = 1.0_rp / ( cdz(k+1) * fdz(k ) * cdz(k ) )
203  cnz3(2,k,2) = 1.0_rp / ( cdz(k+1) * fdz(k ) * cdz(k ) ) &
204  + 1.0_rp / ( cdz(k ) * fdz(k ) * cdz(k ) ) &
205  + 1.0_rp / ( cdz(k ) * fdz(k-1) * cdz(k ) )
206  cnz3(3,k,2) = 1.0_rp / ( cdz(k ) * fdz(k ) * cdz(k ) ) &
207  + 1.0_rp / ( cdz(k ) * fdz(k-1) * cdz(k ) ) &
208  + 1.0_rp / ( cdz(k ) * fdz(k-1) * cdz(k-1) )
209  enddo
210  cnz3(1,ks-1,2) = 1.0_rp / ( cdz(ks ) * fdz(ks ) * cdz(ks ) )
211  cnz3(1,ks ,2) = 1.0_rp / ( cdz(ks+1) * fdz(ks ) * cdz(ks ) )
212  cnz3(2,ks ,2) = 1.0_rp / ( cdz(ks+1) * fdz(ks ) * cdz(ks ) ) &
213  + 1.0_rp / ( cdz(ks ) * fdz(ks ) * cdz(ks ) ) &
214  + 1.0_rp / ( cdz(ks ) * fdz(ks ) * cdz(ks ) )
215  cnz3(3,ks ,2) = 1.0_rp / ( cdz(ks ) * fdz(ks ) * cdz(ks ) ) &
216  + 1.0_rp / ( cdz(ks ) * fdz(ks ) * cdz(ks ) ) &
217  + 1.0_rp / ( cdz(ks ) * fdz(ks ) * cdz(ks+1) )
218  cnz3(1,ke ,2) = 1.0_rp / ( cdz(ke-1) * fdz(ke-1) * cdz(ke ) )
219  cnz3(2,ke ,2) = 1.0_rp / ( cdz(ke-1) * fdz(ke-1) * cdz(ke ) ) &
220  + 1.0_rp / ( cdz(ke ) * fdz(ke-1) * cdz(ke ) ) &
221  + 1.0_rp / ( cdz(ke ) * fdz(ke-1) * cdz(ke ) )
222  cnz3(3,ke ,2) = 1.0_rp / ( cdz(ke ) * fdz(ke-1) * cdz(ke ) ) &
223  + 1.0_rp / ( cdz(ke ) * fdz(ke-1) * cdz(ke ) ) &
224  + 1.0_rp / ( cdz(ke ) * fdz(ke-1) * cdz(ke-1) )
225  cnz3(1,ke+1,2) = 1.0_rp / ( cdz(ke-2) * fdz(ke-2) * cdz(ke-1) )
226  cnz3(2,ke+1,2) = 1.0_rp / ( cdz(ke-2) * fdz(ke-2) * cdz(ke-1) ) &
227  + 1.0_rp / ( cdz(ke-1) * fdz(ke-2) * cdz(ke-1) ) &
228  + 1.0_rp / ( cdz(ke-1) * fdz(ke-1) * cdz(ke-1) )
229  cnz3(3,ke+1,2) = 1.0_rp / ( cdz(ke-1) * fdz(ke-2) * cdz(ke-1) ) &
230  + 1.0_rp / ( cdz(ke-1) * fdz(ke-1) * cdz(ke-1) ) &
231  + 1.0_rp / ( cdz(ke-1) * fdz(ke-1) * cdz(ke ) )
232 
233  do k = ks, ke-1
234  cnz4(1,k,2) = ( cnz3(1,k+1,2) ) / fdz(k)
235  cnz4(2,k,2) = ( cnz3(2,k+1,2) + cnz3(1,k,2) ) / fdz(k)
236  cnz4(3,k,2) = ( cnz3(3,k+1,2) + cnz3(2,k,2) ) / fdz(k)
237  cnz4(4,k,2) = ( cnz3(1,k ,2) + cnz3(3,k,2) ) / fdz(k)
238  cnz4(5,k,2) = ( cnz3(1,k-1,2) ) / fdz(k)
239  enddo
240 ! CNZ4(1,KE,2) = ( CNZ3(1,KE+1,2) ) / FDZ(KE-1)
241  cnz4(2,ke,2) = ( cnz3(2,ke+1,2) + cnz3(1,ke,2) ) / fdz(ke-1)
242  cnz4(3,ke,2) = ( cnz3(3,ke+1,2) + cnz3(2,ke,2) ) / fdz(ke-1)
243  cnz4(4,ke,2) = ( cnz3(1,ke ,2) + cnz3(3,ke,2) ) / fdz(ke-1)
244 
245  ! x direction
246  cnx3(1,is-1,1) = 1.0_rp / ( fdx(is-1) * cdx(is-1) * fdx(is-2) )
247  do i = is, ie+1
248  cnx3(1,i,1) = 1.0_rp / ( fdx(i ) * cdx(i ) * fdx(i-1) )
249  cnx3(2,i,1) = 1.0_rp / ( fdx(i ) * cdx(i ) * fdx(i-1) ) &
250  + 1.0_rp / ( fdx(i-1) * cdx(i ) * fdx(i-1) ) &
251  + 1.0_rp / ( fdx(i-1) * cdx(i-1) * fdx(i-1) )
252  cnx3(3,i,1) = 1.0_rp / ( fdx(i-1) * cdx(i ) * fdx(i-1) ) &
253  + 1.0_rp / ( fdx(i-1) * cdx(i-1) * fdx(i-1) ) &
254  + 1.0_rp / ( fdx(i-1) * cdx(i-1) * fdx(i-2) )
255  enddo
256 
257  do i = is, ie
258  cnx4(1,i,1) = ( cnx3(1,i+1,1) ) / cdx(i)
259  cnx4(2,i,1) = ( cnx3(2,i+1,1) + cnx3(1,i,1) ) / cdx(i)
260  cnx4(3,i,1) = ( cnx3(3,i+1,1) + cnx3(2,i,1) ) / cdx(i)
261  cnx4(4,i,1) = ( cnx3(1,i ,1) + cnx3(3,i,1) ) / cdx(i)
262  cnx4(5,i,1) = ( cnx3(1,i-1,1) ) / cdx(i)
263  enddo
264 
265  do i = is-1, ie+1
266  cnx3(1,i,2) = 1.0_rp / ( cdx(i+1) * fdx(i ) * cdx(i ) )
267  cnx3(2,i,2) = 1.0_rp / ( cdx(i+1) * fdx(i ) * cdx(i ) ) &
268  + 1.0_rp / ( cdx(i ) * fdx(i ) * cdx(i ) ) &
269  + 1.0_rp / ( cdx(i ) * fdx(i-1) * cdx(i ) )
270  cnx3(3,i,2) = 1.0_rp / ( cdx(i ) * fdx(i ) * cdx(i ) ) &
271  + 1.0_rp / ( cdx(i ) * fdx(i-1) * cdx(i ) ) &
272  + 1.0_rp / ( cdx(i ) * fdx(i-1) * cdx(i-1) )
273  enddo
274 
275  do i = is, ie
276  cnx4(1,i,2) = ( cnx3(1,i+1,2) ) / fdx(i)
277  cnx4(2,i,2) = ( cnx3(2,i+1,2) + cnx3(1,i,2) ) / fdx(i)
278  cnx4(3,i,2) = ( cnx3(3,i+1,2) + cnx3(2,i,2) ) / fdx(i)
279  cnx4(4,i,2) = ( cnx3(1,i ,2) + cnx3(3,i,2) ) / fdx(i)
280  cnx4(5,i,2) = ( cnx3(1,i-1,2) ) / fdx(i)
281  enddo
282 
283  ! y direction
284  cny3(1,js-1,1) = 1.0_rp / ( fdy(js-1) * cdy(js-1) * fdy(js-2) )
285  do j = js, je+1
286  cny3(1,j,1) = 1.0_rp / ( fdy(j ) * cdy(j ) * fdy(j-1) )
287  cny3(2,j,1) = 1.0_rp / ( fdy(j ) * cdy(j ) * fdy(j-1) ) &
288  + 1.0_rp / ( fdy(j-1) * cdy(j ) * fdy(j-1) ) &
289  + 1.0_rp / ( fdy(j-1) * cdy(j-1) * fdy(j-1) )
290  cny3(3,j,1) = 1.0_rp / ( fdy(j-1) * cdy(j ) * fdy(j-1) ) &
291  + 1.0_rp / ( fdy(j-1) * cdy(j-1) * fdy(j-1) ) &
292  + 1.0_rp / ( fdy(j-1) * cdy(j-1) * fdy(j-2) )
293  enddo
294 
295  do j = js, je
296  cny4(1,j,1) = ( cny3(1,j+1,1) ) / cdy(j)
297  cny4(2,j,1) = ( cny3(2,j+1,1) + cny3(1,j,1) ) / cdy(j)
298  cny4(3,j,1) = ( cny3(3,j+1,1) + cny3(2,j,1) ) / cdy(j)
299  cny4(4,j,1) = ( cny3(1,j ,1) + cny3(3,j,1) ) / cdy(j)
300  cny4(5,j,1) = ( cny3(1,j-1,1) ) / cdy(j)
301  enddo
302 
303  do j = js-1, je+1
304  cny3(1,j,2) = 1.0_rp / ( cdy(j+1) * fdy(j ) * cdy(j ) )
305  cny3(2,j,2) = 1.0_rp / ( cdy(j+1) * fdy(j ) * cdy(j ) ) &
306  + 1.0_rp / ( cdy(j ) * fdy(j ) * cdy(j ) ) &
307  + 1.0_rp / ( cdy(j ) * fdy(j-1) * cdy(j ) )
308  cny3(3,j,2) = 1.0_rp / ( cdy(j ) * fdy(j ) * cdy(j ) ) &
309  + 1.0_rp / ( cdy(j ) * fdy(j-1) * cdy(j ) ) &
310  + 1.0_rp / ( cdy(j ) * fdy(j-1) * cdy(j-1) )
311  enddo
312 
313  do j = js, je
314  cny4(1,j,2) = ( cny3(1,j+1,2) ) / fdy(j)
315  cny4(2,j,2) = ( cny3(2,j+1,2) + cny3(1,j,2) ) / fdy(j)
316  cny4(3,j,2) = ( cny3(3,j+1,2) + cny3(2,j,2) ) / fdy(j)
317  cny4(4,j,2) = ( cny3(1,j ,2) + cny3(3,j,2) ) / fdy(j)
318  cny4(5,j,2) = ( cny3(1,j-1,2) ) / fdy(j)
319  enddo
320 
321  return
322  end subroutine atmos_dyn_filter_setup
323 
324  !-----------------------------------------------------------------------------
326  subroutine atmos_dyn_wdamp_setup( &
327  wdamp_coef, &
328  wdamp_tau, wdamp_height, &
329  FZ )
330  use scale_const, only: &
331  pi => const_pi, &
332  eps => const_eps
333  implicit none
334 
335  real(RP), intent(inout) :: wdamp_coef(ka)
336  real(RP), intent(in) :: wdamp_tau
337  real(RP), intent(in) :: wdamp_height
338  real(RP), intent(in) :: fz(0:ka)
339 
340  real(RP) :: alpha, sw
341 
342  integer :: k
343  !---------------------------------------------------------------------------
344 
345  if ( wdamp_height < 0.0_rp ) then
346  wdamp_coef(:) = 0.0_rp
347  elseif( fz(ke)-wdamp_height < eps ) then
348  wdamp_coef(:) = 0.0_rp
349  else
350  alpha = 1.0_rp / wdamp_tau
351 
352  do k = ks, ke
353  sw = 0.5_rp + sign( 0.5_rp, fz(k)-wdamp_height )
354 
355  wdamp_coef(k) = alpha * sw &
356  * 0.5_rp * ( 1.0_rp - cos( pi * (fz(k)-wdamp_height) / (fz(ke)-wdamp_height)) )
357  enddo
358  wdamp_coef( 1:ks-1) = wdamp_coef(ks)
359  wdamp_coef(ke+1:ka ) = wdamp_coef(ke)
360 
361  if( io_l ) write(io_fid_log,*)
362  if( io_l ) write(io_fid_log,*) ' *** Setup Rayleigh damping coefficient ***'
363  if( io_l ) write(io_fid_log,'(1x,A)') '|=== Rayleigh Damping Coef ===|'
364  if( io_l ) write(io_fid_log,'(1x,A)') '| k zh[m] coef[/s] |'
365  do k = ka, ke+1, -1
366  if( io_l ) write(io_fid_log,'(1x,A,I5,F10.2,ES12.4,A)') '| ',k, fz(k), wdamp_coef(k),' |'
367  enddo
368  k = ke
369  if( io_l ) write(io_fid_log,'(1x,A,I5,F10.2,ES12.4,A)') '| ',k, fz(k), wdamp_coef(k),' | KE = TOA'
370  do k = ke-1, ks, -1
371  if( io_l ) write(io_fid_log,'(1x,A,I5,F10.2,ES12.4,A)') '| ',k, fz(k), wdamp_coef(k),' |'
372  enddo
373  k = ks-1
374  if( io_l ) write(io_fid_log,'(1x,A,I5,F10.2,ES12.4,A)') '| ',k, fz(k), wdamp_coef(k),' | KS-1 = surface'
375  do k = ks-2, 1, -1
376  if( io_l ) write(io_fid_log,'(1x,A,I5,F10.2,ES12.4,A)') '| ',k, fz(k), wdamp_coef(k),' |'
377  enddo
378  k = 0
379  if( io_l ) write(io_fid_log,'(1x,A,I5,F10.2,12x,A)') '| ',k, fz(k), ' |'
380  if( io_l ) write(io_fid_log,'(1x,A)') '|=============================|'
381  endif
382 
383  return
384  end subroutine atmos_dyn_wdamp_setup
385 
386  !-----------------------------------------------------------------------------
388  subroutine atmos_dyn_numfilter_coef( &
389  num_diff, &
390  DENS, MOMZ, MOMX, MOMY, RHOT, &
391  CDZ, CDX, CDY, FDZ, FDX, FDY, DT, &
392  REF_dens, REF_pott, &
393  ND_COEF, ND_ORDER, ND_SFC_FACT, ND_USE_RS )
394  use scale_comm, only: &
395  comm_vars8, &
396  comm_wait
397  implicit none
398 
399  real(RP), intent(out) :: num_diff(ka,ia,ja,5,3)
400 
401  real(RP), intent(in) :: dens(ka,ia,ja)
402  real(RP), intent(in) :: momz(ka,ia,ja)
403  real(RP), intent(in) :: momx(ka,ia,ja)
404  real(RP), intent(in) :: momy(ka,ia,ja)
405  real(RP), intent(in) :: rhot(ka,ia,ja)
406 
407  real(RP), intent(in) :: cdz(ka)
408  real(RP), intent(in) :: cdx(ia)
409  real(RP), intent(in) :: cdy(ja)
410  real(RP), intent(in) :: fdz(ka-1)
411  real(RP), intent(in) :: fdx(ia-1)
412  real(RP), intent(in) :: fdy(ja-1)
413 
414  real(RP), intent(in) :: dt
415 
416  real(RP), intent(in) :: ref_dens(ka,ia,ja)
417  real(RP), intent(in) :: ref_pott(ka,ia,ja)
418 
419  real(RP), intent(in) :: nd_coef
420  integer, intent(in) :: nd_order
421  real(RP), intent(in) :: nd_sfc_fact
422  logical, intent(in) :: nd_use_rs
423 
424  ! diagnostic variables
425  real(RP) :: velz (ka,ia,ja) ! velocity w [m/s]
426  real(RP) :: velx (ka,ia,ja) ! velocity u [m/s]
427  real(RP) :: vely (ka,ia,ja) ! velocity v [m/s]
428  real(RP) :: pott (ka,ia,ja) ! potential temperature [K]
429 
430  real(RP) :: dens_diff(ka,ia,ja) ! anomary of density
431  real(RP) :: pott_diff(ka,ia,ja) ! anomary of rho * pott
432 
433  real(RP) :: work(ka,ia,ja,3,2)
434  integer :: iwork
435 
436  real(RP) :: diff4
437  integer :: nd_order4
438  real(RP) :: nd_coef_cdz(ka)
439  real(RP) :: nd_coef_cdx(ia)
440  real(RP) :: nd_coef_cdy(ja)
441  real(RP) :: nd_coef_fdz(ka-1)
442  real(RP) :: nd_coef_fdx(ia-1)
443  real(RP) :: nd_coef_fdy(ja-1)
444 
445  integer :: k, i, j
446  !---------------------------------------------------------------------------
447 
448  ! numerical diffusion
449  nd_order4 = nd_order * 4
450  diff4 = nd_coef / ( 2**(nd_order4) * dt )
451  do k = ks-1, ke
452  nd_coef_cdz(k) = diff4 * cdz(k)**nd_order4
453  end do
454  do k = ks+1, ke-1
455  nd_coef_fdz(k) = diff4 * fdz(k)**nd_order4
456  end do
457  do i = is, ie
458  nd_coef_cdx(i) = diff4 * cdx(i)**nd_order4
459  nd_coef_fdx(i) = diff4 * fdx(i)**nd_order4
460  end do
461  do j = js, je
462  nd_coef_cdy(j) = diff4 * cdy(j)**nd_order4
463  nd_coef_fdy(j) = diff4 * fdy(j)**nd_order4
464  end do
465 
466 
467  !###########################################################################
468  ! 1st order coefficients
469  !###########################################################################
470 
471  if ( .NOT. nd_use_rs ) then
472 
473  call prof_rapstart("NumFilter_Main", 3)
474 
475  do j = js-2, je+2
476  do i = is-2, ie+2
477  do k = ks, ke
478  pott(k,i,j) = rhot(k,i,j) / dens(k,i,j)
479  end do
480  end do
481  end do
482 
483  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
484  do j = js, je
485  do i = is, ie
486  do k = ks+1, ke-1
487  dens_diff(k,i,j) = ( ( dens(k,i,j) ) * 3.0_rp &
488  + ( dens(k,i+1,j)+dens(k,i-1,j)+dens(k,i,j+1)+dens(k,i,j-1) ) * 2.0_rp &
489  + ( dens(k,i+2,j)+dens(k,i-2,j)+dens(k,i,j+2)+dens(k,i,j-2) ) &
490  + ( dens(k+1,i,j)+dens(k-1,i,j) ) * 2.0_rp &
491  ) / 19.0_rp
492 
493  pott_diff(k,i,j) = ( ( pott(k,i,j) ) * 3.0_rp &
494  + ( pott(k,i+1,j)+pott(k,i-1,j)+pott(k,i,j+1)+pott(k,i,j-1) ) * 2.0_rp &
495  + ( pott(k,i+2,j)+pott(k,i-2,j)+pott(k,i,j+2)+pott(k,i,j-2) ) &
496  + ( pott(k+1,i,j)+pott(k-1,i,j) ) * 2.0_rp &
497  ) / 19.0_rp
498  enddo
499  enddo
500  enddo
501 
502  do j = js, je
503  do i = is, ie
504  dens_diff(ks,i,j) = ( ( dens(ks,i,j) ) * 3.0_rp &
505  + ( dens(ks,i+1,j)+dens(ks,i-1,j)+dens(ks,i,j+1)+dens(ks,i,j-1) ) * 2.0_rp &
506  + ( dens(ks,i+2,j)+dens(ks,i-2,j)+dens(ks,i,j+2)+dens(ks,i,j-2) ) &
507  + ( dens(ks+1,i,j) ) * 2.0_rp &
508  ) / 17.0_rp
509  dens_diff(ke,i,j) = ( ( dens(ke,i,j) ) * 3.0_rp &
510  + ( dens(ke,i+1,j)+dens(ke,i-1,j)+dens(ke,i,j+1)+dens(ke,i,j-1) ) * 2.0_rp &
511  + ( dens(ke,i+2,j)+dens(ke,i-2,j)+dens(ke,i,j+2)+dens(ke,i,j-2) ) &
512  + ( dens(ke-1,i,j) ) * 2.0_rp &
513  ) / 17.0_rp
514 
515  pott_diff(ks,i,j) = ( ( pott(ks,i,j) ) * 3.0_rp &
516  + ( pott(ks,i+1,j)+pott(ks,i-1,j)+pott(ks,i,j+1)+pott(ks,i,j-1) ) * 2.0_rp &
517  + ( pott(ks,i+2,j)+pott(ks,i-2,j)+pott(ks,i,j+2)+pott(ks,i,j-2) ) &
518  + ( pott(ks+1,i,j) ) * 2.0_rp &
519  ) / 17.0_rp
520  pott_diff(ke,i,j) = ( ( pott(ke,i,j) ) * 3.0_rp &
521  + ( pott(ke,i+1,j)+pott(ke,i-1,j)+pott(ke,i,j+1)+pott(ke,i,j-1) ) * 2.0_rp &
522  + ( pott(ke,i+2,j)+pott(ke,i-2,j)+pott(ke,i,j+2)+pott(ke,i,j-2) ) &
523  + ( pott(ke-1,i,j) ) * 2.0_rp &
524  ) / 17.0_rp
525  end do
526  end do
527 
528  call prof_rapend ("NumFilter_Main", 3)
529 
530  call prof_rapstart("NumFilter_Comm", 3)
531 
532  call comm_vars8( dens_diff, 1 )
533  call comm_vars8( pott_diff, 2 )
534 
535  call comm_wait ( dens_diff, 1 )
536  call comm_wait ( pott_diff, 2 )
537 
538  call prof_rapend ("NumFilter_Comm", 3)
539 
540  end if
541 
542 
543  !-----< density >-----
544 
545  if ( nd_use_rs ) then
546 
547  call prof_rapstart("NumFilter_Main", 3)
548 
549  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
550  do j = js-1, je+2
551  do i = is-1, ie+2
552  do k = ks, ke
553  dens_diff(k,i,j) = dens(k,i,j) - ref_dens(k,i,j)
554  enddo
555  enddo
556  enddo
557 
558  call prof_rapend("NumFilter_Main", 3)
559 
560  endif
561 
562  call calc_numdiff( work, iwork, & ! (out)
563  dens_diff, & ! (in)
564  nd_order, & ! (in)
565  0, 0, 0, ke )
566 
567  call prof_rapstart("NumFilter_Main", 3)
568 
569  !-----< density >-----
570 
571  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
572  do j = js, je
573  do i = is, ie
574  do k = ks-1, ke
575  num_diff(k,i,j,i_dens,zdir) = work(k,i,j,zdir,iwork) * nd_coef_cdz(k)
576  enddo
577  enddo
578  enddo
579  do j = js, je
580  do i = is, ie
581  num_diff( 1:ks-2,i,j,i_dens,zdir) = 0.0_rp
582  num_diff(ke+1:ka ,i,j,i_dens,zdir) = 0.0_rp
583  enddo
584  enddo
585 
586  do j = js, je
587  do i = is, ie
588  do k = ks, ke
589  num_diff(k,i,j,i_dens,xdir) = work(k,i,j,xdir,iwork) * nd_coef_cdx(i)
590  enddo
591  enddo
592  enddo
593  do j = js, je
594  do i = is, ie
595  num_diff( 1:ks-1,i,j,i_dens,xdir) = 0.0_rp
596  num_diff(ks ,i,j,i_dens,xdir) = num_diff(ks ,i,j,i_dens,xdir) * nd_sfc_fact
597  num_diff(ks+1,i,j,i_dens,xdir) = num_diff(ks+1,i,j,i_dens,xdir) * (1.0_rp + nd_sfc_fact) * 0.5_rp
598  num_diff(ke+1:ka ,i,j,i_dens,xdir) = 0.0_rp
599  enddo
600  enddo
601 
602  do j = js, je
603  do i = is, ie
604  do k = ks, ke
605  num_diff(k,i,j,i_dens,ydir) = work(k,i,j,ydir,iwork) * nd_coef_cdy(j)
606  enddo
607  enddo
608  enddo
609  do j = js, je
610  do i = is, ie
611  num_diff( 1:ks-1,i,j,i_dens,ydir) = 0.0_rp
612  num_diff(ks ,i,j,i_dens,ydir) = num_diff(ks ,i,j,i_dens,ydir) * nd_sfc_fact
613  num_diff(ks+1,i,j,i_dens,ydir) = num_diff(ks+1,i,j,i_dens,ydir) * (1.0_rp + nd_sfc_fact) * 0.5_rp
614  num_diff(ke+1:ka ,i,j,i_dens,ydir) = 0.0_rp
615  enddo
616  enddo
617 
618  call prof_rapend ("NumFilter_Main", 3)
619 
620  call prof_rapstart("NumFilter_Comm", 3)
621 
622  call comm_vars8( num_diff(:,:,:,i_dens,zdir), i_comm_dens_z )
623  call comm_vars8( num_diff(:,:,:,i_dens,xdir), i_comm_dens_x )
624  call comm_vars8( num_diff(:,:,:,i_dens,ydir), i_comm_dens_y )
625 
626  call prof_rapend ("NumFilter_Comm", 3)
627 
628 
629  !-----< z-momentum >-----
630 
631  call prof_rapstart("NumFilter_Main", 3)
632 
633  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
634  do j = js-2, je+2
635  do i = is-2, ie+2
636  do k = ks, ke-1
637  velz(k,i,j) = 2.0_rp * momz(k,i,j) / ( dens(k+1,i,j)+dens(k,i,j) )
638  enddo
639  enddo
640  enddo
641 
642  call prof_rapend ("NumFilter_Main", 3)
643 
644  call calc_numdiff( work, iwork, & ! (out)
645  velz, & ! (in)
646  nd_order, & ! (in)
647  1, 0, 0, ke-1 )
648 
649  call prof_rapstart("NumFilter_Main", 3)
650 
651  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
652  do j = js, je
653  do i = is, ie
654  do k = ks+1, ke-1
655  num_diff(k,i,j,i_momz,zdir) = work(k,i,j,zdir,iwork) * nd_coef_fdz(k) &
656  * dens(k,i,j)
657  enddo
658  enddo
659  enddo
660  do j = js, je
661  do i = is, ie
662  num_diff( 1:ks,i,j,i_momz,zdir) = 0.0_rp
663  num_diff(ke:ka,i,j,i_momz,zdir) = 0.0_rp
664  enddo
665  enddo
666 
667  do j = js, je
668  do i = is, ie
669  do k = ks, ke-1
670  num_diff(k,i,j,i_momz,xdir) = work(k,i,j,xdir,iwork) * nd_coef_cdx(i) &
671  * 0.25_rp * ( dens(k+1,i+1,j)+dens(k+1,i,j)+dens(k,i+1,j)+dens(k,i,j) )
672  enddo
673  enddo
674  enddo
675  do j = js, je
676  do i = is, ie
677  num_diff( 1:ks-1,i,j,i_momz,xdir) = 0.0_rp
678  num_diff(ks ,i,j,i_momz,xdir) = num_diff(ks ,i,j,i_momz,xdir) * nd_sfc_fact
679  num_diff(ks+1,i,j,i_momz,xdir) = num_diff(ks+1,i,j,i_momz,xdir) * (1.0_rp + nd_sfc_fact) * 0.5_rp
680  num_diff(ke:ka ,i,j,i_momz,xdir) = 0.0_rp
681  enddo
682  enddo
683 
684  do j = js, je
685  do i = is, ie
686  do k = ks, ke-1
687  num_diff(k,i,j,i_momz,ydir) = work(k,i,j,ydir,iwork) * nd_coef_cdy(j) &
688  * 0.25_rp * ( dens(k+1,i,j+1)+dens(k+1,i,j)+dens(k,i,j+1)+dens(k,i,j) )
689  enddo
690  enddo
691  enddo
692  do j = js, je
693  do i = is, ie
694  num_diff( 1:ks-1,i,j,i_momz,ydir) = 0.0_rp
695  num_diff(ks ,i,j,i_momz,ydir) = num_diff(ks ,i,j,i_momz,ydir) * nd_sfc_fact
696  num_diff(ks+1,i,j,i_momz,ydir) = num_diff(ks+1,i,j,i_momz,ydir) * (1.0_rp + nd_sfc_fact) * 0.5_rp
697 
698  num_diff(ke:ka ,i,j,i_momz,ydir) = 0.0_rp
699  enddo
700  enddo
701 
702  call prof_rapend ("NumFilter_Main", 3)
703 
704  call prof_rapstart("NumFilter_Comm", 3)
705 
706  call comm_vars8( num_diff(:,:,:,i_momz,zdir), i_comm_momz_z )
707  call comm_vars8( num_diff(:,:,:,i_momz,xdir), i_comm_momz_x )
708  call comm_vars8( num_diff(:,:,:,i_momz,ydir), i_comm_momz_y )
709 
710  call prof_rapend ("NumFilter_Comm", 3)
711 
712 
713  !-----< x-momentum >-----
714 
715  call prof_rapstart("NumFilter_Main", 3)
716 
717  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
718  do j = js-2, je+2
719  do i = is-2, ie+1
720  do k = ks, ke
721  velx(k,i,j) = 2.0_rp * momx(k,i,j) / ( dens(k,i+1,j)+dens(k,i,j) )
722  enddo
723  enddo
724  enddo
725 
726  call prof_rapend ("NumFilter_Main", 3)
727 
728  call calc_numdiff( work, iwork, & ! (out)
729  velx, & ! (in)
730  nd_order, & ! (in)
731  0, 1, 0, ke )
732 
733 
734  call prof_rapstart("NumFilter_Main", 3)
735 
736  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
737  do j = js, je
738  do i = is, ie
739  do k = ks, ke-1
740  num_diff(k,i,j,i_momx,zdir) = work(k,i,j,zdir,iwork) * nd_coef_cdz(k) &
741  * 0.25_rp * ( dens(k+1,i+1,j)+dens(k+1,i,j)+dens(k,i+1,j)+dens(k,i,j) )
742  enddo
743  enddo
744  enddo
745  do j = js, je
746  do i = is, ie
747  num_diff( 1:ks-1,i,j,i_momx,zdir) = 0.0_rp
748  num_diff(ke:ka ,i,j,i_momx,zdir) = 0.0_rp
749  enddo
750  enddo
751 
752  do j = js, je
753  do i = is, ie
754  do k = ks, ke
755  num_diff(k,i,j,i_momx,xdir) = work(k,i,j,xdir,iwork) * nd_coef_fdx(i) &
756  * dens(k,i,j)
757  enddo
758  enddo
759  enddo
760  do j = js, je
761  do i = is, ie
762  num_diff( 1:ks-1,i,j,i_momx,xdir) = 0.0_rp
763  num_diff(ks ,i,j,i_momx,xdir) = num_diff(ks ,i,j,i_momx,xdir) * nd_sfc_fact
764  num_diff(ks+1,i,j,i_momx,xdir) = num_diff(ks+1,i,j,i_momx,xdir) * (1.0_rp + nd_sfc_fact) * 0.5_rp
765 
766  num_diff(ke+1:ka ,i,j,i_momx,xdir) = 0.0_rp
767  enddo
768  enddo
769 
770  do j = js, je
771  do i = is, ie
772  do k = ks, ke
773  num_diff(k,i,j,i_momx,ydir) = work(k,i,j,ydir,iwork) * nd_coef_cdy(j) &
774  * 0.25_rp * ( dens(k,i+1,j+1)+dens(k,i+1,j)+dens(k,i,j+1)+dens(k,i,j) )
775  enddo
776  enddo
777  enddo
778  do j = js, je
779  do i = is, ie
780  num_diff( 1:ks-1,i,j,i_momx,ydir) = 0.0_rp
781  num_diff(ks ,i,j,i_momx,ydir) = num_diff(ks ,i,j,i_momx,ydir) * nd_sfc_fact
782  num_diff(ks+1,i,j,i_momx,ydir) = num_diff(ks+1,i,j,i_momx,ydir) * (1.0_rp + nd_sfc_fact) * 0.5_rp
783  num_diff(ke+1:ka ,i,j,i_momx,ydir) = 0.0_rp
784  enddo
785  enddo
786 
787  call prof_rapend ("NumFilter_Main", 3)
788 
789  call prof_rapstart("NumFilter_Comm", 3)
790 
791  call comm_vars8( num_diff(:,:,:,i_momx,zdir), i_comm_momx_z )
792  call comm_vars8( num_diff(:,:,:,i_momx,xdir), i_comm_momx_x )
793  call comm_vars8( num_diff(:,:,:,i_momx,ydir), i_comm_momx_y )
794 
795  call prof_rapend ("NumFilter_Comm", 3)
796 
797 
798  !-----< y-momentum >-----
799 
800  call prof_rapstart("NumFilter_Main", 3)
801 
802  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
803  do j = js-2, je+1
804  do i = is-2, ie+2
805  do k = ks, ke
806  vely(k,i,j) = 2.0_rp * momy(k,i,j) / ( dens(k,i,j+1)+dens(k,i,j) )
807  enddo
808  enddo
809  enddo
810 
811  call prof_rapend ("NumFilter_Main", 3)
812 
813  call calc_numdiff( work, iwork, & ! (out)
814  vely, & ! (in)
815  nd_order, & ! (in)
816  0, 0, 1, ke )
817 
818  call prof_rapstart("NumFilter_Main", 3)
819 
820  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
821  do j = js, je
822  do i = is, ie
823  do k = ks, ke-1
824  num_diff(k,i,j,i_momy,zdir) = work(k,i,j,zdir,iwork) * nd_coef_cdz(k) &
825  * 0.25_rp * ( dens(k+1,i,j+1)+dens(k+1,i,j)+dens(k,i,j+1)+dens(k,i,j) )
826  enddo
827  enddo
828  enddo
829  do j = js, je
830  do i = is, ie
831  num_diff( 1:ks-1,i,j,i_momy,zdir) = 0.0_rp
832  num_diff(ke:ka ,i,j,i_momy,zdir) = 0.0_rp
833  end do
834  end do
835 
836  do j = js, je
837  do i = is, ie
838  do k = ks, ke
839  num_diff(k,i,j,i_momy,xdir) = work(k,i,j,xdir,iwork) * nd_coef_cdx(i) &
840  * 0.25_rp * ( dens(k,i+1,j+1)+dens(k,i,j+1)+dens(k,i+1,j)+dens(k,i,j) )
841  enddo
842  enddo
843  enddo
844  do j = js, je
845  do i = is, ie
846  num_diff( 1:ks-1,i,j,i_momy,xdir) = 0.0_rp
847  num_diff(ks ,i,j,i_momy,xdir) = num_diff(ks ,i,j,i_momy,xdir) * nd_sfc_fact
848  num_diff(ks+1,i,j,i_momy,xdir) = num_diff(ks+1,i,j,i_momy,xdir) * (1.0_rp + nd_sfc_fact) * 0.5_rp
849  num_diff(ke+1:ka ,i,j,i_momy,xdir) = 0.0_rp
850  enddo
851  enddo
852 
853  do j = js, je
854  do i = is, ie
855  do k = ks, ke
856  num_diff(k,i,j,i_momy,ydir) = work(k,i,j,ydir,iwork) * nd_coef_fdy(j) &
857  * dens(k,i,j)
858  enddo
859  enddo
860  enddo
861  do j = js, je
862  do i = is, ie
863  num_diff( 1:ks-1,i,j,i_momy,ydir) = 0.0_rp
864  num_diff(ks ,i,j,i_momy,ydir) = num_diff(ks ,i,j,i_momy,ydir) * nd_sfc_fact
865  num_diff(ks+1,i,j,i_momy,ydir) = num_diff(ks+1,i,j,i_momy,ydir) * (1.0_rp + nd_sfc_fact) * 0.5_rp
866  num_diff(ke+1:ka ,i,j,i_momy,ydir) = 0.0_rp
867  enddo
868  enddo
869 
870  call prof_rapend ("NumFilter_Main", 3)
871 
872  call prof_rapstart("NumFilter_Comm", 3)
873 
874  call comm_vars8( num_diff(:,:,:,i_momy,zdir), i_comm_momy_z )
875  call comm_vars8( num_diff(:,:,:,i_momy,xdir), i_comm_momy_x )
876  call comm_vars8( num_diff(:,:,:,i_momy,ydir), i_comm_momy_y )
877 
878  call prof_rapend ("NumFilter_Comm", 3)
879 
880  !-----< rho * theta >-----
881 
882  call prof_rapstart("NumFilter_Main", 3)
883 
884  if ( nd_use_rs ) then
885  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
886  do j = js-1, je+2
887  do i = is-1, ie+2
888  do k = ks, ke
889  pott_diff(k,i,j) = rhot(k,i,j) / dens(k,i,j) - ref_pott(k,i,j)
890  enddo
891  enddo
892  enddo
893  endif
894 
895  call prof_rapend ("NumFilter_Main", 3)
896 
897  call calc_numdiff( work, iwork, & ! (out)
898  pott_diff, & ! (in)
899  nd_order, & ! (in)
900  0, 0, 0, ke )
901 
902  call prof_rapstart("NumFilter_Main", 3)
903 
904  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
905  do j = js, je
906  do i = is, ie
907  do k = ks, ke-1
908  num_diff(k,i,j,i_rhot,zdir) = work(k,i,j,zdir,iwork) * nd_coef_cdz(k) &
909  * 0.5_rp * ( dens(k+1,i,j)+dens(k,i,j) )
910  enddo
911  enddo
912  enddo
913  !$omp parallel do default(none) &
914  !$omp shared(JS,JE,IS,IE,KS,KE,KA,num_diff,work,DENS,iwork,nd_coef_cdz) &
915  !$omp private(i,j) OMP_SCHEDULE_ collapse(2)
916  do j = js, je
917  do i = is, ie
918  num_diff( 1:ks-2,i,j,i_rhot,zdir) = 0.0_rp
919  num_diff(ks-1,i,j,i_rhot,zdir) = work(ks-1,i,j,zdir,iwork) * nd_coef_cdz(ks-1) &
920  * dens(ks,i,j)
921  num_diff(ke ,i,j,i_rhot,zdir) = work(ke ,i,j,zdir,iwork) * nd_coef_cdz(ke ) &
922  * dens(ke,i,j)
923  num_diff(ke+1:ka ,i,j,i_rhot,zdir) = 0.0_rp
924  enddo
925  enddo
926 
927  do j = js, je
928  do i = is, ie
929  do k = ks, ke
930  num_diff(k,i,j,i_rhot,xdir) = work(k,i,j,xdir,iwork) * nd_coef_cdx(i) &
931  * 0.5_rp * ( dens(k,i+1,j)+dens(k,i,j) )
932  enddo
933  enddo
934  enddo
935  do j = js, je
936  do i = is, ie
937  num_diff( 1:ks-1,i,j,i_rhot,xdir) = 0.0_rp
938  num_diff(ks ,i,j,i_rhot,xdir) = num_diff(ks ,i,j,i_rhot,xdir) * nd_sfc_fact
939  num_diff(ks+1,i,j,i_rhot,xdir) = num_diff(ks+1,i,j,i_rhot,xdir) * (1.0_rp + nd_sfc_fact) * 0.5_rp
940  num_diff(ke+1:ka ,i,j,i_rhot,xdir) = 0.0_rp
941  enddo
942  enddo
943 
944  do j = js, je
945  do i = is, ie
946  do k = ks, ke
947  num_diff(k,i,j,i_rhot,ydir) = work(k,i,j,ydir,iwork) * nd_coef_cdy(j) &
948  * 0.5_rp * ( dens(k,i,j+1)+dens(k,i,j) )
949  enddo
950  enddo
951  enddo
952  do j = js, je
953  do i = is, ie
954  num_diff( 1:ks-1,i,j,i_rhot,ydir) = 0.0_rp
955  num_diff(ks ,i,j,i_rhot,ydir) = num_diff(ks ,i,j,i_rhot,ydir) * nd_sfc_fact
956  num_diff(ks+1,i,j,i_rhot,ydir) = num_diff(ks+1,i,j,i_rhot,ydir) * (1.0_rp + nd_sfc_fact) * 0.5_rp
957  num_diff(ke+1:ka ,i,j,i_rhot,ydir) = 0.0_rp
958  enddo
959  enddo
960 
961  call prof_rapend ("NumFilter_Main", 3)
962 
963  call prof_rapstart("NumFilter_Comm", 3)
964 
965  call comm_vars8( num_diff(:,:,:,i_rhot,zdir), i_comm_rhot_z )
966  call comm_vars8( num_diff(:,:,:,i_rhot,xdir), i_comm_rhot_x )
967  call comm_vars8( num_diff(:,:,:,i_rhot,ydir), i_comm_rhot_y )
968 
969  call comm_wait ( num_diff(:,:,:,i_dens,zdir), i_comm_dens_z )
970  call comm_wait ( num_diff(:,:,:,i_dens,xdir), i_comm_dens_x )
971  call comm_wait ( num_diff(:,:,:,i_dens,ydir), i_comm_dens_y )
972  call comm_wait ( num_diff(:,:,:,i_momz,zdir), i_comm_momz_z )
973  call comm_wait ( num_diff(:,:,:,i_momz,xdir), i_comm_momz_x )
974  call comm_wait ( num_diff(:,:,:,i_momz,ydir), i_comm_momz_y )
975  call comm_wait ( num_diff(:,:,:,i_momx,zdir), i_comm_momx_z )
976  call comm_wait ( num_diff(:,:,:,i_momx,xdir), i_comm_momx_x )
977  call comm_wait ( num_diff(:,:,:,i_momx,ydir), i_comm_momx_y )
978  call comm_wait ( num_diff(:,:,:,i_momy,zdir), i_comm_momy_z )
979  call comm_wait ( num_diff(:,:,:,i_momy,xdir), i_comm_momy_x )
980  call comm_wait ( num_diff(:,:,:,i_momy,ydir), i_comm_momy_y )
981  call comm_wait ( num_diff(:,:,:,i_rhot,zdir), i_comm_rhot_z )
982  call comm_wait ( num_diff(:,:,:,i_rhot,xdir), i_comm_rhot_x )
983  call comm_wait ( num_diff(:,:,:,i_rhot,ydir), i_comm_rhot_y )
984 
985  call prof_rapend ("NumFilter_Comm", 3)
986 
987  return
988  end subroutine atmos_dyn_numfilter_coef
989 
990  !-----------------------------------------------------------------------------
992  subroutine atmos_dyn_numfilter_coef_q( &
993  num_diff_q, &
994  DENS, QTRC, &
995  CDZ, CDX, CDY, dt, &
996  REF_qv, iq, &
997  ND_COEF, ND_ORDER, ND_SFC_FACT, ND_USE_RS )
998  use scale_comm, only: &
999  comm_vars8, &
1000  comm_wait
1001  use scale_atmos_hydrometeor, only: &
1002  i_qv
1003  implicit none
1004 
1005  real(RP), intent(out) :: num_diff_q(ka,ia,ja,3)
1006 
1007  real(RP), intent(in) :: dens(ka,ia,ja)
1008  real(RP), intent(in) :: qtrc(ka,ia,ja)
1009 
1010  real(RP), intent(in) :: cdz(ka)
1011  real(RP), intent(in) :: cdx(ia)
1012  real(RP), intent(in) :: cdy(ja)
1013 
1014  real(RP), intent(in) :: dt
1015 
1016  real(RP), intent(in) :: ref_qv(ka,ia,ja)
1017  integer, intent(in) :: iq
1018 
1019  real(RP), intent(in) :: nd_coef
1020  integer, intent(in) :: nd_order
1021  real(RP), intent(in) :: nd_sfc_fact
1022  logical, intent(in) :: nd_use_rs
1023 
1024  real(RP) :: qv_diff(ka,ia,ja) ! anomary of water vapor
1025 
1026  real(RP) :: work(ka,ia,ja,3,2)
1027  integer :: iwork
1028 
1029  real(RP) :: diff4
1030  integer :: nd_order4
1031  real(RP) :: nd_coef_cdz(ka)
1032  real(RP) :: nd_coef_cdx(ia)
1033  real(RP) :: nd_coef_cdy(ja)
1034 
1035  integer :: k, i, j
1036  !---------------------------------------------------------------------------
1037 
1038  !###########################################################################
1039  ! 1st order coefficients
1040  !###########################################################################
1041 
1042  nd_order4 = nd_order * 4
1043  diff4 = nd_coef / ( 2**(nd_order4) * dt )
1044  do k = ks-1, ke
1045  nd_coef_cdz(k) = diff4 * cdz(k)**nd_order4
1046  end do
1047  do i = is, ie
1048  nd_coef_cdx(i) = diff4 * cdx(i)**nd_order4
1049  end do
1050  do j = js, je
1051  nd_coef_cdy(j) = diff4 * cdy(j)**nd_order4
1052  end do
1053 
1054  if ( iq == i_qv .AND. (.NOT. nd_use_rs) ) then
1055 
1056  call prof_rapstart("NumFilter_Main", 3)
1057 
1058  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
1059  do j = js-1, je+2
1060  do i = is-1, ie+2
1061  do k = ks+1, ke-1
1062  qv_diff(k,i,j) = ( ( qtrc(k,i,j) ) * 3.0_rp &
1063  + ( qtrc(k,i+1,j)+qtrc(k,i-1,j)+qtrc(k,i,j+1)+qtrc(k,i,j-1) ) * 2.0_rp &
1064  + ( qtrc(k,i+2,j)+qtrc(k,i-2,j)+qtrc(k,i,j+2)+qtrc(k,i,j-2) ) &
1065  + ( qtrc(k+1,i,j)+qtrc(k-1,i,j) ) * 2.0_rp &
1066  ) / 19.0_rp
1067  enddo
1068  enddo
1069  enddo
1070 
1071  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
1072  do j = js-1, je+2
1073  do i = is-1, ie+2
1074  qv_diff(ks,i,j) = ( ( qtrc(ks,i,j) ) * 3.0_rp &
1075  + ( qtrc(ks,i+1,j)+qtrc(ks,i-1,j)+qtrc(ks,i,j+1)+qtrc(ks,i,j-1) ) * 2.0_rp &
1076  + ( qtrc(ks,i+2,j)+qtrc(ks,i-2,j)+qtrc(ks,i,j+2)+qtrc(ks,i,j-2) ) &
1077  + ( qtrc(ks+1,i,j) ) * 2.0_rp &
1078  ) / 17.0_rp
1079  qv_diff(ke,i,j) = ( ( qtrc(ke,i,j) ) * 3.0_rp &
1080  + ( qtrc(ke,i+1,j)+qtrc(ke,i-1,j)+qtrc(ke,i,j+1)+qtrc(ke,i,j-1) ) * 2.0_rp &
1081  + ( qtrc(ke,i+2,j)+qtrc(ke,i-2,j)+qtrc(ke,i,j+2)+qtrc(ke,i,j-2) ) &
1082  + ( qtrc(ke-1,i,j) ) * 2.0_rp &
1083  ) / 17.0_rp
1084  end do
1085  end do
1086 
1087  call prof_rapend ("NumFilter_Main", 3)
1088 
1089  call prof_rapstart("NumFilter_Comm", 3)
1090 
1091  call comm_vars8(qv_diff, 1)
1092  call comm_wait (qv_diff, 1)
1093 
1094  call prof_rapend ("NumFilter_Comm", 3)
1095 
1096  end if
1097 
1098  if ( iq == i_qv ) then
1099 
1100  if ( nd_use_rs ) then
1101 
1102  call prof_rapstart("NumFilter_Main", 3)
1103 
1104  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
1105  do j = js-1, je+2
1106  do i = is-1, ie+2
1107  do k = ks, ke
1108  qv_diff(k,i,j) = qtrc(k,i,j) - ref_qv(k,i,j)
1109  enddo
1110  enddo
1111  enddo
1112 
1113  call prof_rapend ("NumFilter_Main", 3)
1114 
1115  endif
1116 
1117  call calc_numdiff( work, iwork, & ! (out)
1118  qv_diff, & ! (in)
1119  nd_order, & ! (in)
1120  0, 0, 0, ke )
1121 
1122  else ! iq /= I_QV
1123 
1124  call calc_numdiff( work, iwork, & ! (out)
1125  qtrc, & ! (in)
1126  nd_order, & ! (in)
1127  0, 0, 0, ke )
1128 
1129  endif ! QV or not?
1130 
1131 
1132  call prof_rapstart("NumFilter_Main", 3)
1133 
1134 
1135  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
1136  do j = js, je
1137  do i = is, ie
1138  do k = ks, ke-1
1139  num_diff_q(k,i,j,zdir) = work(k,i,j,zdir,iwork) * nd_coef_cdz(k) &
1140  * 0.5_rp * ( dens(k+1,i,j)+dens(k,i,j) )
1141  enddo
1142  enddo
1143  enddo
1144  do j = js, je
1145  do i = is, ie
1146  num_diff_q(1:ks-2,i,j,zdir) = 0.0_rp
1147  num_diff_q(ks-1,i,j,zdir) = work(ks-1,i,j,zdir,iwork) * nd_coef_cdz(ks-1) &
1148  * dens(ks,i,j)
1149  num_diff_q(ke ,i,j,zdir) = work(ke ,i,j,zdir,iwork) * nd_coef_cdz(ke ) &
1150  * dens(ke,i,j)
1151  num_diff_q(ke+1:ka,i,j,zdir) = 0.0_rp
1152  enddo
1153  enddo
1154 
1155  do j = js, je
1156  do i = is, ie
1157  do k = ks, ke
1158  num_diff_q(k,i,j,xdir) = work(k,i,j,xdir,iwork) * nd_coef_cdx(i) &
1159  * 0.5_rp * ( dens(k,i+1,j)+dens(k,i,j) )
1160  enddo
1161  enddo
1162  enddo
1163  do j = js, je
1164  do i = is, ie
1165  num_diff_q(1:ks-1,i,j,xdir) = 0.0_rp
1166  num_diff_q(ks ,i,j,xdir) = num_diff_q(ks ,i,j,xdir) * nd_sfc_fact
1167  num_diff_q(ks+1,i,j,xdir) = num_diff_q(ks+1,i,j,xdir) * (1.0_rp + nd_sfc_fact) * 0.5_rp
1168  num_diff_q(ke+1:ka,i,j,xdir) = 0.0_rp
1169  enddo
1170  enddo
1171 
1172  do j = js, je
1173  do i = is, ie
1174  do k = ks, ke
1175  num_diff_q(k,i,j,ydir) = work(k,i,j,ydir,iwork) * nd_coef_cdy(j) &
1176  * 0.5_rp * ( dens(k,i,j+1)+dens(k,i,j) )
1177  enddo
1178  enddo
1179  enddo
1180  do j = js, je
1181  do i = is, ie
1182  num_diff_q(1:ks-1,i,j,ydir) = 0.0_rp
1183  num_diff_q(ks ,i,j,ydir) = num_diff_q(ks ,i,j,ydir) * nd_sfc_fact
1184  num_diff_q(ks+1,i,j,ydir) = num_diff_q(ks+1,i,j,ydir) * (1.0_rp + nd_sfc_fact) * 0.5_rp
1185  num_diff_q(ke+1:ka,i,j,ydir) = 0.0_rp
1186  enddo
1187  enddo
1188 
1189  call prof_rapend ("NumFilter_Main", 3)
1190 
1191  call prof_rapstart("NumFilter_Comm", 3)
1192 
1193  call comm_vars8( num_diff_q(:,:,:,zdir), i_comm_qtrc_z )
1194  call comm_vars8( num_diff_q(:,:,:,xdir), i_comm_qtrc_x )
1195  call comm_vars8( num_diff_q(:,:,:,ydir), i_comm_qtrc_y )
1196 
1197  call comm_wait ( num_diff_q(:,:,:,zdir), i_comm_qtrc_z )
1198  call comm_wait ( num_diff_q(:,:,:,xdir), i_comm_qtrc_x )
1199  call comm_wait ( num_diff_q(:,:,:,ydir), i_comm_qtrc_y )
1200 
1201  call prof_rapend ("NumFilter_Comm", 3)
1202 
1203  return
1204  end subroutine atmos_dyn_numfilter_coef_q
1205 
1206  !-----------------------------------------------------------------------------
1207  subroutine atmos_dyn_filter_tend( &
1208  phi_t, &
1209  phi, &
1210  rdz, rdx, rdy, &
1211  KO, IO, JO )
1212  use scale_comm, only: &
1213  comm_vars8, &
1214  comm_wait
1215  implicit none
1216  real(RP), intent(out) :: phi_t(ka,ia,ja)
1217  real(RP), intent(in ) :: phi (ka,ia,ja)
1218  real(RP), intent(in ) :: rdz(:)
1219  real(RP), intent(in ) :: rdx(:)
1220  real(RP), intent(in ) :: rdy(:)
1221  integer , intent(in ) :: ko
1222  integer , intent(in ) :: io
1223  integer , intent(in ) :: jo
1224 
1225  real(RP) :: flux(ka,ia,ja,3)
1226 
1227  integer :: k, i, j
1228 
1229  call calc_diff3( flux, & ! (out)
1230  phi, & ! (in)
1231  ko, io, jo )
1232 
1233  call comm_vars8( flux(:,:,:,xdir), 1 )
1234  call comm_vars8( flux(:,:,:,ydir), 2 )
1235  call comm_wait ( flux(:,:,:,xdir), 1 )
1236  call comm_wait ( flux(:,:,:,ydir), 2 )
1237 
1238  do j = js, je
1239  do i = is, ie
1240  do k = ks, ke
1241  phi_t(k,i,j) = ( flux(k+ko,i,j,zdir) - flux(k-1+ko,i,j,zdir) ) * rdz(k) &
1242  + ( flux(k,i+io,j,xdir) - flux(k,i-1+io,j,xdir) ) * rdx(i) &
1243  + ( flux(k,i,j+jo,ydir) - flux(k,i,j-1+jo,ydir) ) * rdy(j)
1244  end do
1245  end do
1246  end do
1247 
1248  return
1249  end subroutine atmos_dyn_filter_tend
1250 
1251  !-----------------------------------------------------------------------------
1252  subroutine atmos_dyn_copy_boundary( &
1253  DENS, MOMZ, MOMX, MOMY, RHOT, PROG, &
1254  DENS0, MOMZ0, MOMX0, MOMY0, RHOT0, PROG0, &
1255  BND_W, BND_E, BND_S, BND_N )
1256  implicit none
1257  real(RP), intent(inout) :: dens (ka,ia,ja)
1258  real(RP), intent(inout) :: momz (ka,ia,ja)
1259  real(RP), intent(inout) :: momx (ka,ia,ja)
1260  real(RP), intent(inout) :: momy (ka,ia,ja)
1261  real(RP), intent(inout) :: rhot (ka,ia,ja)
1262  real(RP), intent(inout) :: prog (ka,ia,ja,va)
1263  real(RP), intent(in) :: dens0(ka,ia,ja)
1264  real(RP), intent(in) :: momz0(ka,ia,ja)
1265  real(RP), intent(in) :: momx0(ka,ia,ja)
1266  real(RP), intent(in) :: momy0(ka,ia,ja)
1267  real(RP), intent(in) :: rhot0(ka,ia,ja)
1268  real(RP), intent(in) :: prog0(ka,ia,ja,va)
1269  logical, intent(in) :: bnd_w
1270  logical, intent(in) :: bnd_e
1271  logical, intent(in) :: bnd_s
1272  logical, intent(in) :: bnd_n
1273 
1274  integer :: k, i, j, iv
1275 
1276  if ( bnd_w ) then
1277  !$omp parallel do default(none) private(j,k) OMP_SCHEDULE_ collapse(2) &
1278  !$omp private(i,iv) &
1279  !$omp shared(JA,IS,KS,KE,DENS,DENS0,MOMZ,MOMZ0,MOMX,MOMX0,MOMY,MOMY0,RHOT,RHOT0,VA,PROG,PROG0)
1280 !OCL XFILL
1281  do j = 1, ja
1282  do i = 1, is-1
1283  do k = ks, ke
1284  dens(k,i,j) = dens0(k,i,j)
1285  momz(k,i,j) = momz0(k,i,j)
1286  momx(k,i,j) = momx0(k,i,j)
1287  momy(k,i,j) = momy0(k,i,j)
1288  rhot(k,i,j) = rhot0(k,i,j)
1289  do iv = 1, va
1290  prog(k,i,j,iv) = prog0(k,i,j,iv)
1291  end do
1292  enddo
1293  enddo
1294  enddo
1295  end if
1296  if ( bnd_e ) then
1297  !$omp parallel do default(none) private(j,k) OMP_SCHEDULE_ collapse(2) &
1298  !$omp private(i,iv) &
1299  !$omp shared(JA,IE,IA,KS,KE,DENS,DENS0,MOMZ,MOMZ0,MOMX,MOMX0,MOMY,MOMY0,RHOT,RHOT0,VA,PROG,PROG0)
1300 !OCL XFILL
1301  do j = 1, ja
1302  do i = ie+1, ia
1303  do k = ks, ke
1304  dens(k,i,j) = dens0(k,i,j)
1305  momz(k,i,j) = momz0(k,i,j)
1306  momx(k,i,j) = momx0(k,i,j)
1307  momy(k,i,j) = momy0(k,i,j)
1308  rhot(k,i,j) = rhot0(k,i,j)
1309  do iv = 1, va
1310  prog(k,i,j,iv) = prog0(k,i,j,iv)
1311  end do
1312  enddo
1313  enddo
1314  enddo
1315  !$omp parallel do private(j,k) OMP_SCHEDULE_ collapse(2)
1316 !OCL XFILL
1317  do j = 1, ja
1318  do k = ks, ke
1319  momx(k,ie,j) = momx0(k,ie,j)
1320  enddo
1321  enddo
1322  end if
1323  if ( bnd_s ) then
1324  !$omp parallel do default(none) private(j,k) OMP_SCHEDULE_ collapse(2) &
1325  !$omp private(i,iv) &
1326  !$omp shared(JS,IA,KS,KE,DENS,DENS0,MOMZ,MOMZ0,MOMX,MOMX0,MOMY,MOMY0,RHOT,RHOT0,VA,PROG,PROG0)
1327 !OCL XFILL
1328  do j = 1, js-1
1329  do i = 1, ia
1330  do k = ks, ke
1331  dens(k,i,j) = dens0(k,i,j)
1332  momz(k,i,j) = momz0(k,i,j)
1333  momx(k,i,j) = momx0(k,i,j)
1334  momy(k,i,j) = momy0(k,i,j)
1335  rhot(k,i,j) = rhot0(k,i,j)
1336  do iv = 1, va
1337  prog(k,i,j,iv) = prog0(k,i,j,iv)
1338  end do
1339  enddo
1340  enddo
1341  enddo
1342  end if
1343  if ( bnd_n ) then
1344  !$omp parallel do default(none) private(j,k) OMP_SCHEDULE_ collapse(2) &
1345  !$omp private(i,iv) &
1346  !$omp shared(JA,JE,IA,KS,KE,DENS,DENS0,MOMZ,MOMZ0,MOMX,MOMX0,MOMY,MOMY0,RHOT,RHOT0,VA,PROG,PROG0)
1347 !OCL XFILL
1348  do j = je+1, ja
1349  do i = 1, ia
1350  do k = ks, ke
1351  dens(k,i,j) = dens0(k,i,j)
1352  momz(k,i,j) = momz0(k,i,j)
1353  momx(k,i,j) = momx0(k,i,j)
1354  momy(k,i,j) = momy0(k,i,j)
1355  rhot(k,i,j) = rhot0(k,i,j)
1356  do iv = 1, va
1357  prog(k,i,j,iv) = prog0(k,i,j,iv)
1358  end do
1359  enddo
1360  enddo
1361  enddo
1362  !$omp parallel do private(i,k) OMP_SCHEDULE_ collapse(2)
1363 !OCL XFILL
1364  do i = 1, ia
1365  do k = ks, ke
1366  momy(k,i,je) = momy0(k,i,je)
1367  enddo
1368  enddo
1369  end if
1370 
1371  return
1372  end subroutine atmos_dyn_copy_boundary
1373 
1374  !-----------------------------------------------------------------------------
1375  subroutine atmos_dyn_copy_boundary_tracer( &
1376  QTRC, QTRC0, &
1377  BND_W, BND_E, BND_S, BND_N )
1378  implicit none
1379  real(RP), intent(inout) :: qtrc (ka,ia,ja)
1380  real(RP), intent(in) :: qtrc0(ka,ia,ja)
1381  logical, intent(in) :: bnd_w
1382  logical, intent(in) :: bnd_e
1383  logical, intent(in) :: bnd_s
1384  logical, intent(in) :: bnd_n
1385 
1386  integer :: k, i, j
1387 
1388  if ( bnd_w ) then
1389  !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
1390  !$omp shared(JA,IS,KS,KE,QTRC,QTRC0)
1391 !OCL XFILL
1392  do j = 1, ja
1393  do i = 1, is-1
1394  do k = ks, ke
1395  qtrc(k,i,j) = qtrc0(k,i,j)
1396  enddo
1397  enddo
1398  enddo
1399  end if
1400  if ( bnd_e ) then
1401  !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
1402  !$omp shared(JA,IE,IA,KS,KE,QTRC,QTRC0)
1403 !OCL XFILL
1404  do j = 1, ja
1405  do i = ie+1, ia
1406  do k = ks, ke
1407  qtrc(k,i,j) = qtrc0(k,i,j)
1408  enddo
1409  enddo
1410  enddo
1411  end if
1412  if ( bnd_s ) then
1413  !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
1414  !$omp shared(JS,IA,KS,KE,QTRC,QTRC0)
1415 !OCL XFILL
1416  do j = 1, js-1
1417  do i = 1, ia
1418  do k = ks, ke
1419  qtrc(k,i,j) = qtrc0(k,i,j)
1420  enddo
1421  enddo
1422  enddo
1423  end if
1424  if ( bnd_n ) then
1425  !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
1426  !$omp shared(JA,JE,IA,KS,KE,QTRC,QTRC0)
1427 !OCL XFILL
1428  do j = je+1, ja
1429  do i = 1, ia
1430  do k = ks, ke
1431  qtrc(k,i,j) = qtrc0(k,i,j)
1432  enddo
1433  enddo
1434  enddo
1435  end if
1436 
1437  return
1438  end subroutine atmos_dyn_copy_boundary_tracer
1439 
1440  !-----------------------------------------------------------------------------
1441  subroutine atmos_dyn_divergence( &
1442  DDIV, &
1443  MOMZ, MOMX, MOMY, &
1444  GSQRT, J13G, J23G, J33G, MAPF, &
1445  RCDZ, RCDX, RCDY, RFDZ, FDZ )
1446  use scale_gridtrans, only: &
1447  i_xyz, &
1448  i_xyw, &
1449  i_uyz, &
1450  i_xvz, &
1451  i_xy, &
1452  i_uy, &
1453  i_xv
1454  implicit none
1455  real(RP), intent(out) :: ddiv(ka,ia,ja)
1456  real(RP), intent(in) :: momz(ka,ia,ja)
1457  real(RP), intent(in) :: momx(ka,ia,ja)
1458  real(RP), intent(in) :: momy(ka,ia,ja)
1459  real(RP), intent(in) :: gsqrt(ka,ia,ja,7)
1460  real(RP), intent(in) :: j13g(ka,ia,ja,7)
1461  real(RP), intent(in) :: j23g(ka,ia,ja,7)
1462  real(RP), intent(in) :: j33g
1463  real(RP), intent(in) :: mapf(ia,ja,2,7)
1464  real(RP), intent(in) :: rcdz(ka)
1465  real(RP), intent(in) :: rcdx(ia)
1466  real(RP), intent(in) :: rcdy(ja)
1467  real(RP), intent(in) :: rfdz(ka-1)
1468  real(RP), intent(in) :: fdz(ka-1)
1469 
1470  integer :: k, i, j
1471 
1472  call prof_rapstart("DYN_divercence", 2)
1473 
1474  ! 3D divergence
1475 
1476  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
1477  do j = js, je+1
1478  do i = is, ie+1
1479  do k = ks-1, ke+1
1480  ddiv(k,i,j) = j33g * ( momz(k,i,j) - momz(k-1,i ,j ) ) * rcdz(k) &
1481  + ( ( momx(k+1,i,j) + momx(k+1,i-1,j ) ) * j13g(k+1,i,j,i_xyw) &
1482  - ( momx(k-1,i,j) + momx(k-1,i-1,j ) ) * j13g(k-1,i,j,i_xyw) &
1483  + ( momy(k+1,i,j) + momy(k+1,i ,j-1) ) * j23g(k+1,i,j,i_xyw) &
1484  - ( momy(k-1,i,j) + momy(k-1,i ,j-1) ) * j23g(k-1,i,j,i_xyw) ) / ( fdz(k)+fdz(k-1) ) &
1485  + mapf(i,j,1,i_xy) * mapf(i,j,2,i_xy) &
1486  * ( ( momx(k,i ,j ) * gsqrt(k,i ,j ,i_uyz) / mapf(i ,j ,2,i_uy) &
1487  - momx(k,i-1,j ) * gsqrt(k,i-1,j ,i_uyz) / mapf(i-1,j ,2,i_uy) ) * rcdx(i) &
1488  + ( momy(k,i ,j ) * gsqrt(k,i ,j ,i_xvz) / mapf(i ,j ,1,i_xv) &
1489  - momy(k,i, j-1) * gsqrt(k,i ,j-1,i_xvz) / mapf(i ,j-1,1,i_xv) ) * rcdy(j) )
1490  enddo
1491  enddo
1492  enddo
1493 #ifdef DEBUG
1494  k = iundef; i = iundef; j = iundef
1495 #endif
1496  !$omp parallel do private(i,j) OMP_SCHEDULE_ collapse(2)
1497  do j = js, je+1
1498  do i = is, ie+1
1499  ddiv(ks,i,j) = j33g * ( momz(ks,i,j) ) * rcdz(ks) &
1500  + ( ( momx(ks+1,i,j) + momx(ks+1,i-1,j ) ) * j13g(ks+1,i,j,i_xyw) &
1501  - ( momx(ks-1,i,j) + momx(ks ,i-1,j ) ) * j13g(ks ,i,j,i_xyw) &
1502  + ( momy(ks+1,i,j) + momy(ks+1,i ,j-1) ) * j23g(ks+1,i,j,i_xyw) &
1503  - ( momy(ks ,i,j) + momy(ks ,i ,j-1) ) * j23g(ks ,i,j,i_xyw) ) * rfdz(ks) &
1504  + mapf(i,j,1,i_xy) * mapf(i,j,2,i_xy) &
1505  * ( ( momx(ks,i ,j ) * gsqrt(ks,i ,j ,i_uyz) / mapf(i ,j ,2,i_uy) &
1506  - momx(ks,i-1,j ) * gsqrt(ks,i-1,j ,i_uyz) / mapf(i-1,j ,2,i_uy) ) * rcdx(i) &
1507  + ( momy(ks,i ,j ) * gsqrt(ks,i ,j ,i_xvz) / mapf(i ,j ,1,i_xv) &
1508  - momy(ks,i, j-1) * gsqrt(ks,i ,j-1,i_xvz) / mapf(i ,j-1,1,i_xv) ) * rcdy(j) )
1509  ddiv(ke,i,j) = j33g * ( - momz(ke-1,i ,j ) ) * rcdz(ke) &
1510  + ( ( momx(ke ,i,j) + momx(ke ,i-1,j ) ) * j13g(ke ,i,j,i_xyw) &
1511  - ( momx(ke-1,i,j) + momx(ke-1,i-1,j ) ) * j13g(ke-1,i,j,i_xyw) &
1512  + ( momy(ke ,i,j) + momy(ke ,i ,j-1) ) * j23g(ke ,i,j,i_xyw) &
1513  - ( momy(ke-1,i,j) + momy(ke-1,i ,j-1) ) * j23g(ke-1,i,j,i_xyw) ) * rfdz(ke) &
1514  + mapf(i,j,1,i_xy) * mapf(i,j,2,i_xy) &
1515  * ( ( momx(ke,i ,j ) * gsqrt(ke,i ,j ,i_uyz) / mapf(i ,j ,2,i_uy) &
1516  - momx(ke,i-1,j ) * gsqrt(ke,i-1,j ,i_uyz) / mapf(i-1,j ,2,i_uy) ) * rcdx(i) &
1517  + ( momy(ke,i ,j ) * gsqrt(ke,i ,j ,i_xvz) / mapf(i ,j ,1,i_xv) &
1518  - momy(ke,i, j-1) * gsqrt(ke,i ,j-1,i_xvz) / mapf(i ,j-1,1,i_xv) ) * rcdy(j) )
1519  enddo
1520  enddo
1521 #ifdef DEBUG
1522  k = iundef; i = iundef; j = iundef
1523 #endif
1524  call prof_rapend ("DYN_divercence", 2)
1525 
1526  return
1527  end subroutine atmos_dyn_divergence
1528 
1529  !-----------------------------------------------------------------------------
1530  subroutine calc_numdiff(&
1531  work, iwork, &
1532  data, &
1533  nd_order, &
1534  KO, IO, JO, KEE )
1535  use scale_comm, only: &
1536  comm_vars8, &
1537  comm_wait
1538  implicit none
1539  real(RP), intent(out) :: work(KA,IA,JA,3,2)
1540  integer, intent(out) :: iwork
1541  real(RP), intent(in) :: data(KA,IA,JA)
1542  integer, intent(in) :: nd_order
1543  integer, intent(in) :: KO
1544  integer, intent(in) :: IO
1545  integer, intent(in) :: JO
1546  integer, intent(in) :: KEE
1547 
1548  integer :: i_in, i_out, i_tmp
1549 
1550  integer :: no
1551 
1552  call prof_rapstart("NumFilter_Main", 3)
1553 
1554  call calc_diff3( work(:,:,:,:,1), & ! (out)
1555  data, & ! (in)
1556  ko, io, jo ) ! (in)
1557 
1558  call prof_rapend ("NumFilter_Main", 3)
1559 
1560  !###########################################################################
1561  ! High order coefficients
1562  !###########################################################################
1563 
1564  i_in = 1
1565  i_out = 2
1566 
1567  do no = 2, nd_order
1568 
1569  call prof_rapstart("NumFilter_Comm", 3)
1570 
1571  call comm_vars8( work(:,:,:,zdir,i_in), 16 )
1572  call comm_vars8( work(:,:,:,xdir,i_in), 17 )
1573  call comm_vars8( work(:,:,:,ydir,i_in), 18 )
1574 
1575  call comm_wait ( work(:,:,:,zdir,i_in), 16 )
1576  call comm_wait ( work(:,:,:,xdir,i_in), 17 )
1577  call comm_wait ( work(:,:,:,ydir,i_in), 18 )
1578 
1579  call prof_rapend ("NumFilter_Comm", 3)
1580 
1581  call prof_rapstart("NumFilter_Main", 3)
1582 
1583  call calc_diff4( work(:,:,:,:,i_out), & ! (out)
1584  work(:,:,:,:,i_in), & ! (in)
1585  cnz4(:,:,1+ko), & ! (in)
1586  cnx4(:,:,1+io), & ! (in)
1587  cny4(:,:,1+jo), & ! (in)
1588  kee ) ! (in)
1589 
1590  call prof_rapend ("NumFilter_Main", 3)
1591 
1592  ! swap pointer target
1593  i_tmp = i_in
1594  i_in = i_out
1595  i_out = i_tmp
1596  enddo
1597 
1598  iwork = i_in
1599 
1600  return
1601  end subroutine calc_numdiff
1602 
1603  !-----------------------------------------------------------------------------
1604  subroutine calc_diff3( &
1605  diff, &
1606  phi, &
1607  KO, IO, JO )
1608  implicit none
1609  real(RP), intent(out) :: diff(KA,IA,JA,3)
1610  real(RP), intent(in ) :: phi(KA,IA,JA)
1611  integer , intent(in ) :: KO
1612  integer , intent(in ) :: IO
1613  integer , intent(in ) :: JO
1614 
1615  integer :: kee
1616  integer :: k, i, j
1617 
1618  kee = ke-ko
1619 
1620  if ( ko == 0 ) then
1621 
1622  !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
1623  !$omp shared(JS,JE,IS,IE,KS,KE,phi,diff,CNZ3)
1624  do j = js, je
1625  do i = is, ie
1626  do k = ks+1, ke-2
1627 #ifdef DEBUG
1628  call check( __line__, phi(k+2,i,j) )
1629  call check( __line__, phi(k+1,i,j) )
1630  call check( __line__, phi(k ,i,j) )
1631  call check( __line__, phi(k-1,i,j) )
1632 #endif
1633  diff(k,i,j,zdir) = ( + cnz3(1,k+1,1) * phi(k+2,i,j) &
1634  - cnz3(2,k+1,1) * phi(k+1,i,j) &
1635  + cnz3(3,k+1,1) * phi(k ,i,j) &
1636  - cnz3(1,k ,1) * phi(k-1,i,j) )
1637  enddo
1638  enddo
1639  enddo
1640 
1641  !$omp parallel do private(i,j) OMP_SCHEDULE_ collapse(2)
1642  do j = js, je
1643  do i = is, ie
1644 #ifdef DEBUG
1645  call check( __line__, phi(ks+2,i,j) )
1646  call check( __line__, phi(ks+1,i,j) )
1647  call check( __line__, phi(ks,i,j) )
1648  call check( __line__, phi(ke,i,j) )
1649  call check( __line__, phi(ke-1,i,j) )
1650  call check( __line__, phi(ke-2,i,j) )
1651 #endif
1652  diff(ks,i,j,zdir) = ( + cnz3(1,ks+1,1) * phi(ks+2,i,j) &
1653  - cnz3(2,ks+1,1) * phi(ks+1,i,j) &
1654  + cnz3(3,ks+1,1) * phi(ks ,i,j) &
1655  - cnz3(1,ks ,1) * phi(ks+1,i,j) )
1656  diff(ks-1,i,j,zdir) = - diff(ks ,i,j,zdir)
1657  diff(ks-2,i,j,zdir) = - diff(ks+1,i,j,zdir)
1658  diff(ke-1,i,j,zdir) = ( + cnz3(1,ke ,1) * phi(ke-1,i,j) &
1659  - cnz3(2,ke ,1) * phi(ke ,i,j) &
1660  + cnz3(3,ke ,1) * phi(ke-1,i,j) &
1661  - cnz3(1,ke-1,1) * phi(ke-2,i,j) )
1662  diff(ke ,i,j,zdir) = - diff(ke-1,i,j,zdir)
1663  diff(ke+1,i,j,zdir) = - diff(ke-2,i,j,zdir)
1664  diff(ke+2,i,j,zdir) = 0.0_rp
1665  end do
1666  end do
1667 
1668  else ! K0=1
1669 
1670  !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
1671  !$omp shared(JS,JE,IS,IE,KS,KE,phi,diff,CNZ3)
1672  do j = js, je
1673  do i = is, ie
1674  do k = ks+2, ke-2
1675 #ifdef DEBUG
1676  call check( __line__, phi(k+1,i,j) )
1677  call check( __line__, phi(k ,i,j) )
1678  call check( __line__, phi(k-1,i,j) )
1679  call check( __line__, phi(k-2,i,j) )
1680 #endif
1681  diff(k,i,j,zdir) = ( + cnz3(1,k ,2) * phi(k+1,i,j) &
1682  - cnz3(2,k ,2) * phi(k ,i,j) &
1683  + cnz3(3,k ,2) * phi(k-1,i,j) &
1684  - cnz3(1,k-1,2) * phi(k-2,i,j) )
1685  enddo
1686  enddo
1687  enddo
1688 
1689  !$omp parallel do private(i,j) OMP_SCHEDULE_ collapse(2)
1690  do j = js, je
1691  do i = is, ie
1692 #ifdef DEBUG
1693  call check( __line__, phi(ks+2,i,j) )
1694  call check( __line__, phi(ks+1,i,j) )
1695  call check( __line__, phi(ks,i,j) )
1696  call check( __line__, phi(ks+1,i,j) )
1697  call check( __line__, phi(ks ,i,j) )
1698  call check( __line__, phi(ke-1,i,j) )
1699  call check( __line__, phi(ke-2,i,j) )
1700  call check( __line__, phi(ke-3,i,j) )
1701 #endif
1702  diff(ks+1,i,j,zdir) = ( + cnz3(1,ks+1,2) * phi(ks+2,i,j) &
1703  - cnz3(2,ks+1,2) * phi(ks+1,i,j) &
1704  + cnz3(3,ks+1,2) * phi(ks ,i,j) &
1705  - cnz3(1,ks ,2) * phi(ks+1,i,j) )
1706  diff(ks ,i,j,zdir) = - diff(ks+1,i,j,zdir)
1707  diff(ks-1,i,j,zdir) = - diff(ks+2,i,j,zdir)
1708  diff(ks-2,i,j,zdir) = - diff(ks+3,i,j,zdir)
1709  diff(ke-1,i,j,zdir) = ( - cnz3(2,ke-1,2) * phi(ke-1,i,j) &
1710  + cnz3(3,ke-1,2) * phi(ke-2,i,j) &
1711  - cnz3(1,ke-2,2) * phi(ke-3,i,j) )
1712  diff(ke ,i,j,zdir) = ( + cnz3(1,ke ,2) * phi(ke-1,i,j) &
1713  + cnz3(3,ke ,2) * phi(ke-1,i,j) &
1714  - cnz3(1,ke-1,2) * phi(ke-2,i,j) )
1715  diff(ke+1,i,j,zdir) = - diff(ke ,i,j,zdir)
1716  diff(ke+2,i,j,zdir) = - diff(ke-1,i,j,zdir)
1717  end do
1718  end do
1719 
1720  end if
1721 
1722  if ( io == 0 ) then
1723  !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
1724  !$omp shared(JS,JE,IS,IE,KS,KEE,phi,diff,CNX3)
1725  do j = js, je
1726  do i = is, ie
1727  do k = ks, kee
1728 #ifdef DEBUG
1729  call check( __line__, phi(k,i+2,j) )
1730  call check( __line__, phi(k,i+1,j) )
1731  call check( __line__, phi(k,i ,j) )
1732  call check( __line__, phi(k,i-1,j) )
1733 #endif
1734  diff(k,i,j,xdir) = ( + cnx3(1,i+1,1) * phi(k,i+2,j) &
1735  - cnx3(2,i+1,1) * phi(k,i+1,j) &
1736  + cnx3(3,i+1,1) * phi(k,i ,j) &
1737  - cnx3(1,i ,1) * phi(k,i-1,j) )
1738  enddo
1739  enddo
1740  enddo
1741  else
1742  !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ collapse(2) &
1743  !$omp shared(JS,JE,IS,IE,KS,KEE,phi,diff,CNX3)
1744  do j = js, je
1745  do i = is, ie
1746  do k = ks, kee
1747 #ifdef DEBUG
1748  call check( __line__, phi(k,i+1,j) )
1749  call check( __line__, phi(k,i ,j) )
1750  call check( __line__, phi(k,i-1,j) )
1751  call check( __line__, phi(k,i-2,j) )
1752 #endif
1753  diff(k,i,j,xdir) = ( + cnx3(1,i ,2) * phi(k,i+1,j) &
1754  - cnx3(2,i ,2) * phi(k,i ,j) &
1755  + cnx3(3,i ,2) * phi(k,i-1,j) &
1756  - cnx3(1,i-1,2) * phi(k,i-2,j) )
1757  enddo
1758  enddo
1759  enddo
1760  end if
1761 
1762  !$omp parallel do private(i,j) OMP_SCHEDULE_ collapse(2)
1763  do j = js, je
1764  do i = is, ie
1765  diff( 1:ks-1,i,j,xdir) = 0.0_rp
1766  diff(ke+1:ka ,i,j,xdir) = 0.0_rp
1767  enddo
1768  enddo
1769 
1770  if ( jo == 0 ) then
1771  !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ &
1772  !$omp shared(JS,JE,IS,IE,KS,KEE,phi,diff,CNY3)
1773  do j = js, je
1774  do i = is, ie
1775  do k = ks, kee
1776 #ifdef DEBUG
1777  call check( __line__, phi(k,i,j+2) )
1778  call check( __line__, phi(k,i,j+1) )
1779  call check( __line__, phi(k,i,j ) )
1780  call check( __line__, phi(k,i,j-1) )
1781 #endif
1782  diff(k,i,j,ydir) = ( + cny3(1,j+1,1) * phi(k,i,j+2) &
1783  - cny3(2,j+1,1) * phi(k,i,j+1) &
1784  + cny3(3,j+1,1) * phi(k,i,j ) &
1785  - cny3(1,j ,1) * phi(k,i,j-1) )
1786  enddo
1787  enddo
1788  enddo
1789  else
1790  !$omp parallel do default(none) private(i,j,k) OMP_SCHEDULE_ &
1791  !$omp shared(JS,JE,IS,IE,KS,KEE,phi,diff,CNY3)
1792  do j = js, je
1793  do i = is, ie
1794  do k = ks, kee
1795 #ifdef DEBUG
1796  call check( __line__, phi(k,i,j+1) )
1797  call check( __line__, phi(k,i,j ) )
1798  call check( __line__, phi(k,i,j-1) )
1799  call check( __line__, phi(k,i,j-2) )
1800 #endif
1801  diff(k,i,j,ydir) = ( + cny3(1,j ,2) * phi(k,i,j+1) &
1802  - cny3(2,j ,2) * phi(k,i,j ) &
1803  + cny3(3,j ,2) * phi(k,i,j-1) &
1804  - cny3(1,j-1,2) * phi(k,i,j-2) )
1805  enddo
1806  enddo
1807  enddo
1808  end if
1809 
1810  !$omp parallel do private(i,j) OMP_SCHEDULE_ collapse(2)
1811  do j = js, je
1812  do i = is, ie
1813  diff( 1:ks-1,i,j,ydir) = 0.0_rp
1814  diff(ke+1:ka ,i,j,ydir) = 0.0_rp
1815  enddo
1816  enddo
1817 
1818  return
1819  end subroutine calc_diff3
1820 
1821  !-----------------------------------------------------------------------------
1822  subroutine calc_diff4( &
1823  num_diff_pt1, &
1824  num_diff_pt0, &
1825  CNZ4, &
1826  CNX4, &
1827  CNY4, &
1828  k1 )
1829  implicit none
1830 
1831  real(RP), intent(out) :: num_diff_pt1(KA,IA,JA,3)
1832  real(RP), intent(in) :: num_diff_pt0(KA,IA,JA,3)
1833  real(RP), intent(in) :: CNZ4(5,KA)
1834  real(RP), intent(in) :: CNX4(5,IA)
1835  real(RP), intent(in) :: CNY4(5,JA)
1836  integer, intent(in) :: k1
1837 
1838  integer :: i, j, k
1839  !---------------------------------------------------------------------------
1840 
1841  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
1842  do j = js, je
1843  do i = is, ie
1844  do k = ks, ke-1
1845 #ifdef DEBUG
1846  call check( __line__, cnz4(1,k) )
1847  call check( __line__, cnz4(2,k) )
1848  call check( __line__, cnz4(3,k) )
1849  call check( __line__, cnz4(4,k) )
1850  call check( __line__, cnz4(5,k) )
1851  call check( __line__, num_diff_pt0(k+2,i,j,zdir) )
1852  call check( __line__, num_diff_pt0(k+1,i,j,zdir) )
1853  call check( __line__, num_diff_pt0(k ,i,j,zdir) )
1854  call check( __line__, num_diff_pt0(k-1,i,j,zdir) )
1855  call check( __line__, num_diff_pt0(k-2,i,j,zdir) )
1856 #endif
1857  num_diff_pt1(k,i,j,zdir) = &
1858  ( cnz4(1,k) * num_diff_pt0(k+2,i,j,zdir) &
1859  - cnz4(2,k) * num_diff_pt0(k+1,i,j,zdir) &
1860  + cnz4(3,k) * num_diff_pt0(k ,i,j,zdir) &
1861  - cnz4(4,k) * num_diff_pt0(k-1,i,j,zdir) &
1862  + cnz4(5,k) * num_diff_pt0(k-2,i,j,zdir) )
1863  enddo
1864  enddo
1865  enddo
1866 
1867  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
1868  do j = js, je
1869  do i = is, ie
1870  num_diff_pt1(ks-1,i,j,zdir) = - num_diff_pt1(ks ,i,j,zdir)
1871  num_diff_pt1(ks-2,i,j,zdir) = - num_diff_pt1(ks+1,i,j,zdir)
1872  num_diff_pt1(ke ,i,j,zdir) = - num_diff_pt1(ke-1,i,j,zdir)
1873  num_diff_pt1(ke+1,i,j,zdir) = - num_diff_pt1(ke-2,i,j,zdir)
1874  enddo
1875  enddo
1876 
1877  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
1878  do j = js, je
1879  do i = is, ie
1880  do k = ks, k1
1881 #ifdef DEBUG
1882  call check( __line__, cnx4(1,i) )
1883  call check( __line__, cnx4(2,i) )
1884  call check( __line__, cnx4(3,i) )
1885  call check( __line__, cnx4(4,i) )
1886  call check( __line__, cnx4(5,i) )
1887  call check( __line__, num_diff_pt0(k,i-2,j,xdir) )
1888  call check( __line__, num_diff_pt0(k,i+1,j,xdir) )
1889  call check( __line__, num_diff_pt0(k,i ,j,xdir) )
1890  call check( __line__, num_diff_pt0(k,i-1,j,xdir) )
1891  call check( __line__, num_diff_pt0(k,i-2,j,xdir) )
1892 #endif
1893  num_diff_pt1(k,i,j,xdir) = &
1894  ( cnx4(1,i) * num_diff_pt0(k,i+2,j,xdir) &
1895  - cnx4(2,i) * num_diff_pt0(k,i+1,j,xdir) &
1896  + cnx4(3,i) * num_diff_pt0(k,i ,j,xdir) &
1897  - cnx4(4,i) * num_diff_pt0(k,i-1,j,xdir) &
1898  + cnx4(5,i) * num_diff_pt0(k,i-2,j,xdir) )
1899  enddo
1900  enddo
1901  enddo
1902 
1903  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
1904  do j = js, je
1905  do i = is, ie
1906  do k = ks, k1
1907 #ifdef DEBUG
1908  call check( __line__, cny4(1,j) )
1909  call check( __line__, cny4(2,j) )
1910  call check( __line__, cny4(3,j) )
1911  call check( __line__, cny4(4,j) )
1912  call check( __line__, cny4(5,j) )
1913  call check( __line__, num_diff_pt0(k,i,j-2,ydir) )
1914  call check( __line__, num_diff_pt0(k,i,j+1,ydir) )
1915  call check( __line__, num_diff_pt0(k,i,j ,ydir) )
1916  call check( __line__, num_diff_pt0(k,i,j-1,ydir) )
1917  call check( __line__, num_diff_pt0(k,i,j-2,ydir) )
1918 #endif
1919  num_diff_pt1(k,i,j,ydir) = &
1920  ( cny4(1,j) * num_diff_pt0(k,i,j+2,ydir) &
1921  - cny4(2,j) * num_diff_pt0(k,i,j+1,ydir) &
1922  + cny4(3,j) * num_diff_pt0(k,i,j ,ydir) &
1923  - cny4(4,j) * num_diff_pt0(k,i,j-1,ydir) &
1924  + cny4(5,j) * num_diff_pt0(k,i,j-2,ydir) )
1925  enddo
1926  enddo
1927  enddo
1928 
1929  return
1930  end subroutine calc_diff4
1931 
1932  !-----------------------------------------------------------------------------
1934  subroutine atmos_dyn_fct( &
1935  qflx_anti, &
1936  phi_in, DENS0, DENS, &
1937  qflx_hi, qflx_lo, &
1938  mflx_hi, &
1939  rdz, rdx, rdy, &
1940  GSQRT, MAPF, dt, &
1941  flag_vect )
1942  use scale_const, only: &
1943  undef => const_undef, &
1944  iundef => const_undef2, &
1945  epsilon => const_eps
1946  use scale_comm, only: &
1947  comm_vars8, &
1948  comm_wait
1949  implicit none
1950 
1951  real(RP), intent(out) :: qflx_anti(ka,ia,ja,3)
1952 
1953  real(RP), intent(in) :: phi_in(ka,ia,ja) ! physical quantity
1954  real(RP), intent(in) :: dens0(ka,ia,ja)
1955  real(RP), intent(in) :: dens (ka,ia,ja)
1956 
1957  real(RP), intent(in) :: qflx_hi(ka,ia,ja,3)
1958  real(RP), intent(in) :: qflx_lo(ka,ia,ja,3)
1959  real(RP), intent(in) :: mflx_hi(ka,ia,ja,3)
1960 
1961  real(RP), intent(in) :: rdz(:)
1962  real(RP), intent(in) :: rdx(:)
1963  real(RP), intent(in) :: rdy(:)
1964 
1965  real(RP), intent(in) :: gsqrt(ka,ia,ja)
1966  real(RP), intent(in) :: mapf(ia,ja,2)
1967 
1968  real(RP), intent(in) :: dt
1969 
1970  logical, intent(in) :: flag_vect
1971 
1972  ! work for FCT
1973  real(RP) :: phi_lo(ka,ia,ja)
1974  real(RP) :: pjpls(ka,ia,ja)
1975  real(RP) :: pjmns(ka,ia,ja)
1976  real(RP) :: qjpls(ka,ia,ja)
1977  real(RP) :: qjmns(ka,ia,ja)
1978  real(RP) :: rjpls(ka,ia,ja)
1979  real(RP) :: rjmns(ka,ia,ja)
1980 
1981  real(RP) :: qmin, qmax
1982  real(RP) :: zerosw, dirsw
1983 
1984  real(RP) :: fact(0:1,-1:1,-1:1)
1985  real(RP) :: rw, ru, rv
1986  real(RP) :: qa_in, qb_in
1987  real(RP) :: qa_lo, qb_lo
1988 
1989  integer :: k, i, j, ijs
1990  integer :: iis, iie, jjs, jje
1991  !---------------------------------------------------------------------------
1992 
1993 #ifdef DEBUG
1994  qflx_anti(:,:,:,:) = undef
1995 
1996  pjpls(:,:,:) = undef
1997  pjmns(:,:,:) = undef
1998  qjpls(:,:,:) = undef
1999  qjmns(:,:,:) = undef
2000  rjpls(:,:,:) = undef
2001  rjmns(:,:,:) = undef
2002 #endif
2003 
2004  do jjs = js, je, jblock
2005  jje = jjs+jblock-1
2006  do iis = is, ie, iblock
2007  iie = iis+iblock-1
2008 
2009  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
2010  do j = jjs, jje
2011  do i = iis, iie
2012  do k = ks, ke-1
2013 #ifdef DEBUG
2014  call check( __line__, qflx_hi(k,i,j,zdir) )
2015  call check( __line__, qflx_lo(k,i,j,zdir) )
2016 #endif
2017  qflx_anti(k,i,j,zdir) = qflx_hi(k,i,j,zdir) - qflx_lo(k,i,j,zdir)
2018  enddo
2019  enddo
2020  enddo
2021 #ifdef DEBUG
2022  k = iundef; i = iundef; j = iundef
2023 #endif
2024  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
2025  do j = jjs, jje
2026  do i = iis, iie
2027  qflx_anti(ks-1,i,j,zdir) = 0.0_rp
2028  qflx_anti(ke ,i,j,zdir) = 0.0_rp
2029  enddo
2030  enddo
2031 #ifdef DEBUG
2032  k = iundef; i = iundef; j = iundef
2033 #endif
2034  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
2035  do j = jjs , jje
2036  do i = iis-1, iie
2037  do k = ks, ke
2038 #ifdef DEBUG
2039  call check( __line__, qflx_hi(k,i,j,xdir) )
2040  call check( __line__, qflx_lo(k,i,j,xdir) )
2041 #endif
2042  qflx_anti(k,i,j,xdir) = qflx_hi(k,i,j,xdir) - qflx_lo(k,i,j,xdir)
2043  enddo
2044  enddo
2045  enddo
2046 #ifdef DEBUG
2047  k = iundef; i = iundef; j = iundef
2048 #endif
2049  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
2050  do j = jjs-1, jje
2051  do i = iis , iie
2052  do k = ks, ke
2053 #ifdef DEBUG
2054  call check( __line__, qflx_hi(k,i,j,ydir) )
2055  call check( __line__, qflx_lo(k,i,j,ydir) )
2056 #endif
2057  qflx_anti(k,i,j,ydir) = qflx_hi(k,i,j,ydir) - qflx_lo(k,i,j,ydir)
2058  enddo
2059  enddo
2060  enddo
2061 #ifdef DEBUG
2062  k = iundef; i = iundef; j = iundef
2063 #endif
2064 
2065  !--- update monotone scheme
2066  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
2067  do j = jjs-1, jje+1
2068  do i = iis-1, iie+1
2069  do k = ks, ke
2070 #ifdef DEBUG
2071  call check( __line__, phi_in(k,i,j) )
2072  call check( __line__, qflx_lo(k ,i ,j ,zdir) )
2073  call check( __line__, qflx_lo(k-1,i ,j ,zdir) )
2074  call check( __line__, qflx_lo(k ,i ,j ,xdir) )
2075  call check( __line__, qflx_lo(k ,i-1,j ,xdir) )
2076  call check( __line__, qflx_lo(k ,i ,j ,ydir) )
2077  call check( __line__, qflx_lo(k ,i ,j-1,ydir) )
2078 #endif
2079  phi_lo(k,i,j) = ( phi_in(k,i,j) * dens0(k,i,j) &
2080  + dt * ( - ( ( qflx_lo(k,i,j,zdir)-qflx_lo(k-1,i ,j ,zdir) ) * rdz(k) &
2081  + ( qflx_lo(k,i,j,xdir)-qflx_lo(k ,i-1,j ,xdir) ) * rdx(i) &
2082  + ( qflx_lo(k,i,j,ydir)-qflx_lo(k ,i ,j-1,ydir) ) * rdy(j) &
2083  ) * mapf(i,j,1) * mapf(i,j,2) / gsqrt(k,i,j) ) &
2084  ) / dens(k,i,j)
2085  enddo
2086  enddo
2087  enddo
2088 #ifdef DEBUG
2089  k = iundef; i = iundef; j = iundef
2090 #endif
2091 
2092  !--- calc net incoming quantity change by antidiffusive flux
2093  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
2094  do j = jjs, jje
2095  do i = iis, iie
2096  do k = ks, ke
2097 #ifdef DEBUG
2098  call check( __line__, qflx_anti(k ,i ,j ,zdir) )
2099  call check( __line__, qflx_anti(k-1,i ,j ,zdir) )
2100  call check( __line__, qflx_anti(k ,i ,j ,xdir) )
2101  call check( __line__, qflx_anti(k ,i-1,j ,xdir) )
2102  call check( __line__, qflx_anti(k ,i ,j ,ydir) )
2103  call check( __line__, qflx_anti(k ,i ,j-1,ydir) )
2104 #endif
2105  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) &
2106  + ( max(0.0_rp,qflx_anti(k ,i-1,j ,xdir)) - min(0.0_rp,qflx_anti(k,i,j,xdir)) ) * rdx(i) &
2107  + ( max(0.0_rp,qflx_anti(k ,i ,j-1,ydir)) - min(0.0_rp,qflx_anti(k,i,j,ydir)) ) * rdy(j) &
2108  ) * mapf(i,j,1) * mapf(i,j,2) / gsqrt(k,i,j)
2109  enddo
2110  enddo
2111  enddo
2112 #ifdef DEBUG
2113  k = iundef; i = iundef; j = iundef
2114 #endif
2115 
2116  !--- calc net outgoing quantity change by antidiffusive flux
2117  !$omp parallel do private(i,j,k) OMP_SCHEDULE_ collapse(2)
2118  do j = jjs, jje
2119  do i = iis, iie
2120  do k = ks, ke
2121 #ifdef DEBUG
2122  call check( __line__, qflx_anti(k ,i ,j ,zdir) )
2123  call check( __line__, qflx_anti(k-1,i ,j ,zdir) )
2124  call check( __line__, qflx_anti(k ,i ,j ,xdir) )
2125  call check( __line__, qflx_anti(k ,i-1,j ,xdir) )
2126  call check( __line__, qflx_anti(k ,i ,j ,ydir) )
2127  call check( __line__, qflx_anti(k ,i ,j-1,ydir) )
2128 #endif
2129  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) &
2130  + ( max(0.0_rp,qflx_anti(k,i,j,xdir)) - min(0.0_rp,qflx_anti(k ,i-1,j ,xdir)) ) * rdx(i) &
2131  + ( max(0.0_rp,qflx_anti(k,i,j,ydir)) - min(0.0_rp,qflx_anti(k ,i ,j-1,ydir)) ) * rdy(j) &
2132  ) * mapf(i,j,1) * mapf(i,j,2) / gsqrt(k,i,j)
2133  enddo
2134  enddo
2135  enddo
2136 #ifdef DEBUG
2137  k = iundef; i = iundef; j = iundef
2138 #endif
2139 
2140  !--- calc allowable range or quantity change by antidiffusive flux
2141 
2142  if (flag_vect) then
2143 
2144  !$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)
2145  do j = jjs, jje
2146  do i = iis, iie
2147  do k = ks+1, ke-1
2148  rw = (mflx_hi(k,i,j,zdir)+mflx_hi(k-1,i ,j ,zdir)) * rdz(k) ! 2 * rho * w / dz
2149  ru = (mflx_hi(k,i,j,xdir)+mflx_hi(k ,i-1,j ,xdir)) * rdx(i) ! 2 * rho * u / dx
2150  rv = (mflx_hi(k,i,j,ydir)+mflx_hi(k ,i ,j-1,ydir)) * rdy(j) ! 2 * rho * v / dy
2151 
2152  call get_fact_fct( fact, & ! (out)
2153  rw, ru, rv ) ! (in)
2154 
2155  qa_in = fact(1, 1, 1) * phi_in(k+1,i+1,j+1) &
2156  + fact(0, 1, 1) * phi_in(k ,i+1,j+1) &
2157  + fact(1, 0, 1) * phi_in(k+1,i ,j+1) &
2158  + fact(0, 0, 1) * phi_in(k ,i ,j+1) &
2159  + fact(1,-1, 1) * phi_in(k+1,i-1,j+1) &
2160  + fact(1, 1, 0) * phi_in(k+1,i+1,j ) &
2161  + fact(0, 1, 0) * phi_in(k ,i+1,j ) &
2162  + fact(1, 0, 0) * phi_in(k+1,i ,j ) &
2163  + fact(1,-1, 0) * phi_in(k+1,i-1,j ) &
2164  + fact(1, 1,-1) * phi_in(k+1,i+1,j-1) &
2165  + fact(0, 1,-1) * phi_in(k ,i+1,j-1) &
2166  + fact(1, 0,-1) * phi_in(k+1,i ,j-1) &
2167  + fact(1,-1,-1) * phi_in(k+1,i-1,j-1) &
2168  + fact(0, 0, 0) * phi_in(k ,i ,j )
2169  qb_in = fact(1, 1, 1) * phi_in(k-1,i-1,j-1) &
2170  + fact(0, 1, 1) * phi_in(k ,i-1,j-1) &
2171  + fact(1, 0, 1) * phi_in(k-1,i ,j-1) &
2172  + fact(0, 0, 1) * phi_in(k ,i ,j-1) &
2173  + fact(1,-1, 1) * phi_in(k-1,i+1,j-1) &
2174  + fact(1, 1, 0) * phi_in(k-1,i-1,j ) &
2175  + fact(0, 1, 0) * phi_in(k ,i-1,j ) &
2176  + fact(1, 0, 0) * phi_in(k-1,i ,j ) &
2177  + fact(1,-1, 0) * phi_in(k-1,i+1,j ) &
2178  + fact(1, 1,-1) * phi_in(k-1,i-1,j+1) &
2179  + fact(0, 1,-1) * phi_in(k ,i-1,j-1) &
2180  + fact(1, 0,-1) * phi_in(k-1,i ,j-1) &
2181  + fact(1,-1,-1) * phi_in(k-1,i+1,j+1) &
2182  + fact(0, 0, 0) * phi_in(k ,i ,j )
2183  qa_lo = fact(1, 1, 1) * phi_lo(k+1,i+1,j+1) &
2184  + fact(0, 1, 1) * phi_lo(k ,i+1,j+1) &
2185  + fact(1, 0, 1) * phi_lo(k+1,i ,j+1) &
2186  + fact(0, 0, 1) * phi_lo(k ,i ,j+1) &
2187  + fact(1,-1, 1) * phi_lo(k+1,i-1,j+1) &
2188  + fact(1, 1, 0) * phi_lo(k+1,i+1,j ) &
2189  + fact(0, 1, 0) * phi_lo(k ,i+1,j ) &
2190  + fact(1, 0, 0) * phi_lo(k+1,i ,j ) &
2191  + fact(1,-1, 0) * phi_lo(k+1,i-1,j ) &
2192  + fact(1, 1,-1) * phi_lo(k+1,i+1,j-1) &
2193  + fact(0, 1,-1) * phi_lo(k ,i+1,j-1) &
2194  + fact(1, 0,-1) * phi_lo(k+1,i ,j-1) &
2195  + fact(1,-1,-1) * phi_lo(k+1,i-1,j-1) &
2196  + fact(0, 0, 0) * phi_lo(k ,i ,j )
2197  qb_lo = fact(1, 1, 1) * phi_lo(k-1,i-1,j-1) &
2198  + fact(0, 1, 1) * phi_lo(k ,i-1,j-1) &
2199  + fact(1, 0, 1) * phi_lo(k-1,i ,j-1) &
2200  + fact(0, 0, 1) * phi_lo(k ,i ,j-1) &
2201  + fact(1,-1, 1) * phi_lo(k-1,i+1,j-1) &
2202  + fact(1, 1, 0) * phi_lo(k-1,i-1,j ) &
2203  + fact(0, 1, 0) * phi_lo(k ,i-1,j ) &
2204  + fact(1, 0, 0) * phi_lo(k-1,i ,j ) &
2205  + fact(1,-1, 0) * phi_lo(k-1,i+1,j ) &
2206  + fact(1, 1,-1) * phi_lo(k-1,i-1,j+1) &
2207  + fact(0, 1,-1) * phi_lo(k ,i-1,j-1) &
2208  + fact(1, 0,-1) * phi_lo(k-1,i ,j-1) &
2209  + fact(1,-1,-1) * phi_lo(k-1,i+1,j+1) &
2210  + fact(0, 0, 0) * phi_lo(k ,i ,j )
2211 
2212  qmax = max( &
2213  phi_in(k,i,j), qa_in, qb_in, &
2214  phi_lo(k,i,j), qa_lo, qb_lo )
2215  qmin = min( &
2216  phi_in(k,i,j), qa_in, qb_in, &
2217  phi_lo(k,i,j), qa_lo, qb_lo )
2218  qjpls(k,i,j) = ( qmax - phi_lo(k,i,j) ) * dens(k,i,j)
2219  qjmns(k,i,j) = ( phi_lo(k,i,j) - qmin ) * dens(k,i,j)
2220  end do
2221  end do
2222  end do
2223 
2224  ! k = KS
2225  !$omp parallel do private(i,j,rw,ru,rv,fact,qa_in,qb_in,qa_lo,qb_lo,qmax,qmin) OMP_SCHEDULE_ collapse(2)
2226  do j = jjs, jje
2227  do i = iis, iie
2228  rw = (mflx_hi(ks,i,j,zdir) ) * rdz(ks)! 2 * rho * w / dz
2229  ru = (mflx_hi(ks,i,j,xdir)+mflx_hi(ks ,i-1,j ,xdir)) * rdx(i) ! 2 * rho * u / dx
2230  rv = (mflx_hi(ks,i,j,ydir)+mflx_hi(ks ,i ,j-1,ydir)) * rdy(j) ! 2 * rho * v / dy
2231 
2232  call get_fact_fct( fact, & ! (out)
2233  rw, ru, rv ) ! (in)
2234 
2235  qa_in = fact(1, 1, 1) * phi_in(ks+1,i+1,j+1) &
2236  + fact(0, 1, 1) * phi_in(ks ,i+1,j+1) &
2237  + fact(1, 0, 1) * phi_in(ks+1,i ,j+1) &
2238  + fact(0, 0, 1) * phi_in(ks ,i ,j+1) &
2239  + fact(1,-1, 1) * phi_in(ks+1,i-1,j+1) &
2240  + fact(1, 1, 0) * phi_in(ks+1,i+1,j ) &
2241  + fact(0, 1, 0) * phi_in(ks ,i+1,j ) &
2242  + fact(1, 0, 0) * phi_in(ks+1,i ,j ) &
2243  + fact(1,-1, 0) * phi_in(ks+1,i-1,j ) &
2244  + fact(1, 1,-1) * phi_in(ks+1,i+1,j-1) &
2245  + fact(0, 1,-1) * phi_in(ks ,i+1,j-1) &
2246  + fact(1, 0,-1) * phi_in(ks+1,i ,j-1) &
2247  + fact(1,-1,-1) * phi_in(ks+1,i-1,j-1) &
2248  + fact(0, 0, 0) * phi_in(ks ,i ,j )
2249  qb_in = fact(1, 1, 1) * phi_in(ks ,i-1,j-1) &
2250  + fact(0, 1, 1) * phi_in(ks ,i-1,j-1) &
2251  + fact(1, 0, 1) * phi_in(ks ,i ,j-1) &
2252  + fact(0, 0, 1) * phi_in(ks ,i ,j-1) &
2253  + fact(1,-1, 1) * phi_in(ks ,i+1,j-1) &
2254  + fact(1, 1, 0) * phi_in(ks ,i-1,j ) &
2255  + fact(0, 1, 0) * phi_in(ks ,i-1,j ) &
2256  + fact(1, 0, 0) * phi_in(ks ,i ,j ) &
2257  + fact(1,-1, 0) * phi_in(ks ,i+1,j ) &
2258  + fact(1, 1,-1) * phi_in(ks ,i-1,j+1) &
2259  + fact(0, 1,-1) * phi_in(ks ,i-1,j-1) &
2260  + fact(1, 0,-1) * phi_in(ks ,i ,j-1) &
2261  + fact(1,-1,-1) * phi_in(ks ,i+1,j+1) &
2262  + fact(0, 0, 0) * phi_in(ks ,i ,j )
2263  qa_lo = fact(1, 1, 1) * phi_lo(ks+1,i+1,j+1) &
2264  + fact(0, 1, 1) * phi_lo(ks ,i+1,j+1) &
2265  + fact(1, 0, 1) * phi_lo(ks+1,i ,j+1) &
2266  + fact(0, 0, 1) * phi_lo(ks ,i ,j+1) &
2267  + fact(1,-1, 1) * phi_lo(ks+1,i-1,j+1) &
2268  + fact(1, 1, 0) * phi_lo(ks+1,i+1,j ) &
2269  + fact(0, 1, 0) * phi_lo(ks ,i+1,j ) &
2270  + fact(1, 0, 0) * phi_lo(ks+1,i ,j ) &
2271  + fact(1,-1, 0) * phi_lo(ks+1,i-1,j ) &
2272  + fact(1, 1,-1) * phi_lo(ks+1,i+1,j-1) &
2273  + fact(0, 1,-1) * phi_lo(ks ,i+1,j-1) &
2274  + fact(1, 0,-1) * phi_lo(ks+1,i ,j-1) &
2275  + fact(1,-1,-1) * phi_lo(ks+1,i-1,j-1) &
2276  + fact(0, 0, 0) * phi_lo(ks ,i ,j )
2277  qb_lo = fact(1, 1, 1) * phi_lo(ks ,i-1,j-1) &
2278  + fact(0, 1, 1) * phi_lo(ks ,i-1,j-1) &
2279  + fact(1, 0, 1) * phi_lo(ks ,i ,j-1) &
2280  + fact(0, 0, 1) * phi_lo(ks ,i ,j-1) &
2281  + fact(1,-1, 1) * phi_lo(ks ,i+1,j-1) &
2282  + fact(1, 1, 0) * phi_lo(ks ,i-1,j ) &
2283  + fact(0, 1, 0) * phi_lo(ks ,i-1,j ) &
2284  + fact(1, 0, 0) * phi_lo(ks ,i ,j ) &
2285  + fact(1,-1, 0) * phi_lo(ks ,i+1,j ) &
2286  + fact(1, 1,-1) * phi_lo(ks ,i-1,j+1) &
2287  + fact(0, 1,-1) * phi_lo(ks ,i-1,j-1) &
2288  + fact(1, 0,-1) * phi_lo(ks ,i ,j-1) &
2289  + fact(1,-1,-1) * phi_lo(ks ,i+1,j+1) &
2290  + fact(0, 0, 0) * phi_lo(ks ,i ,j )
2291 
2292  qmax = max( &
2293  phi_in(ks,i,j), qa_in, qb_in, &
2294  phi_lo(ks,i,j), qa_lo, qb_lo )
2295  qmin = min( &
2296  phi_in(ks,i,j), qa_in, qb_in, &
2297  phi_lo(ks,i,j), qa_lo, qb_lo )
2298  qjpls(ks,i,j) = ( qmax - phi_lo(ks,i,j) ) * dens(ks,i,j)
2299  qjmns(ks,i,j) = ( phi_lo(ks,i,j) - qmin ) * dens(ks,i,j)
2300  end do
2301  end do
2302 
2303  ! k = KE
2304  !$omp parallel do private(i,j,rw,ru,rv,fact,qa_in,qb_in,qa_lo,qb_lo,qmax,qmin) OMP_SCHEDULE_ collapse(2)
2305  do j = jjs, jje
2306  do i = iis, iie
2307  rw = ( mflx_hi(ke-1,i ,j ,zdir)) * rdz(ke)! 2 * rho * w / dz
2308  ru = (mflx_hi(ke,i,j,xdir)+mflx_hi(ke ,i-1,j ,xdir)) * rdx(i) ! 2 * rho * u / dx
2309  rv = (mflx_hi(ke,i,j,ydir)+mflx_hi(ke ,i ,j-1,ydir)) * rdy(j) ! 2 * rho * v / dy
2310 
2311  call get_fact_fct( fact, & ! (out)
2312  rw, ru, rv ) ! (in)
2313 
2314  qa_in = fact(1, 1, 1) * phi_in(ke ,i+1,j+1) &
2315  + fact(0, 1, 1) * phi_in(ke ,i+1,j+1) &
2316  + fact(1, 0, 1) * phi_in(ke ,i ,j+1) &
2317  + fact(0, 0, 1) * phi_in(ke ,i ,j+1) &
2318  + fact(1,-1, 1) * phi_in(ke ,i-1,j+1) &
2319  + fact(1, 1, 0) * phi_in(ke ,i+1,j ) &
2320  + fact(0, 1, 0) * phi_in(ke ,i+1,j ) &
2321  + fact(1, 0, 0) * phi_in(ke ,i ,j ) &
2322  + fact(1,-1, 0) * phi_in(ke ,i-1,j ) &
2323  + fact(1, 1,-1) * phi_in(ke ,i+1,j-1) &
2324  + fact(0, 1,-1) * phi_in(ke ,i+1,j-1) &
2325  + fact(1, 0,-1) * phi_in(ke ,i ,j-1) &
2326  + fact(1,-1,-1) * phi_in(ke ,i-1,j-1) &
2327  + fact(0, 0, 0) * phi_in(ke ,i ,j )
2328  qb_in = fact(1, 1, 1) * phi_in(ke-1,i-1,j-1) &
2329  + fact(0, 1, 1) * phi_in(ke ,i-1,j-1) &
2330  + fact(1, 0, 1) * phi_in(ke-1,i ,j-1) &
2331  + fact(0, 0, 1) * phi_in(ke ,i ,j-1) &
2332  + fact(1,-1, 1) * phi_in(ke-1,i+1,j-1) &
2333  + fact(1, 1, 0) * phi_in(ke-1,i-1,j ) &
2334  + fact(0, 1, 0) * phi_in(ke ,i-1,j ) &
2335  + fact(1, 0, 0) * phi_in(ke-1,i ,j ) &
2336  + fact(1,-1, 0) * phi_in(ke-1,i+1,j ) &
2337  + fact(1, 1,-1) * phi_in(ke-1,i-1,j+1) &
2338  + fact(0, 1,-1) * phi_in(ke ,i-1,j-1) &
2339  + fact(1, 0,-1) * phi_in(ke-1,i ,j-1) &
2340  + fact(1,-1,-1) * phi_in(ke-1,i+1,j+1) &
2341  + fact(0, 0, 0) * phi_in(ke ,i ,j )
2342  qa_lo = fact(1, 1, 1) * phi_lo(ke ,i+1,j+1) &
2343  + fact(0, 1, 1) * phi_lo(ke ,i+1,j+1) &
2344  + fact(1, 0, 1) * phi_lo(ke ,i ,j+1) &
2345  + fact(0, 0, 1) * phi_lo(ke ,i ,j+1) &
2346  + fact(1,-1, 1) * phi_lo(ke ,i-1,j+1) &
2347  + fact(1, 1, 0) * phi_lo(ke ,i+1,j ) &
2348  + fact(0, 1, 0) * phi_lo(ke ,i+1,j ) &
2349  + fact(1, 0, 0) * phi_lo(ke ,i ,j ) &
2350  + fact(1,-1, 0) * phi_lo(ke ,i-1,j ) &
2351  + fact(1, 1,-1) * phi_lo(ke ,i+1,j-1) &
2352  + fact(0, 1,-1) * phi_lo(ke ,i+1,j-1) &
2353  + fact(1, 0,-1) * phi_lo(ke ,i ,j-1) &
2354  + fact(1,-1,-1) * phi_lo(ke ,i-1,j-1) &
2355  + fact(0, 0, 0) * phi_lo(ke ,i ,j )
2356  qb_lo = fact(1, 1, 1) * phi_lo(ke-1,i-1,j-1) &
2357  + fact(0, 1, 1) * phi_lo(ke ,i-1,j-1) &
2358  + fact(1, 0, 1) * phi_lo(ke-1,i ,j-1) &
2359  + fact(0, 0, 1) * phi_lo(ke ,i ,j-1) &
2360  + fact(1,-1, 1) * phi_lo(ke-1,i+1,j-1) &
2361  + fact(1, 1, 0) * phi_lo(ke-1,i-1,j ) &
2362  + fact(0, 1, 0) * phi_lo(ke ,i-1,j ) &
2363  + fact(1, 0, 0) * phi_lo(ke-1,i ,j ) &
2364  + fact(1,-1, 0) * phi_lo(ke-1,i+1,j ) &
2365  + fact(1, 1,-1) * phi_lo(ke-1,i-1,j+1) &
2366  + fact(0, 1,-1) * phi_lo(ke ,i-1,j-1) &
2367  + fact(1, 0,-1) * phi_lo(ke-1,i ,j-1) &
2368  + fact(1,-1,-1) * phi_lo(ke-1,i+1,j+1) &
2369  + fact(0, 0, 0) * phi_lo(ke ,i ,j )
2370 
2371  qmax = max( &
2372  phi_in(ke,i,j), qa_in, qb_in, &
2373  phi_lo(ke,i,j), qa_lo, qb_lo )
2374  qmin = min( &
2375  phi_in(ke,i,j), qa_in, qb_in, &
2376  phi_lo(ke,i,j), qa_lo, qb_lo )
2377  qjpls(ke,i,j) = ( qmax - phi_lo(ke,i,j) ) * dens(ke,i,j)
2378  qjmns(ke,i,j) = ( phi_lo(ke,i,j) - qmin ) * dens(ke,i,j)
2379  end do
2380  end do
2381 
2382  else
2383 
2384  !$omp parallel do private(i,j,k,qmax,qmin) OMP_SCHEDULE_ collapse(2)
2385  do j = jjs, jje
2386  do i = iis, iie
2387  do k = ks+1, ke-1
2388 #ifdef DEBUG
2389  call check( __line__, phi_in(k ,i ,j ) )
2390  call check( __line__, phi_in(k-1,i ,j ) )
2391  call check( __line__, phi_in(k+1,i ,j ) )
2392  call check( __line__, phi_in(k ,i-1,j ) )
2393  call check( __line__, phi_in(k ,i+1,j ) )
2394  call check( __line__, phi_in(k ,i ,j+1) )
2395  call check( __line__, phi_in(k ,i ,j-1) )
2396  call check( __line__, phi_lo(k ,i ,j ) )
2397  call check( __line__, phi_lo(k-1,i ,j ) )
2398  call check( __line__, phi_lo(k+1,i ,j ) )
2399  call check( __line__, phi_lo(k ,i-1,j ) )
2400  call check( __line__, phi_lo(k ,i+1,j ) )
2401  call check( __line__, phi_lo(k ,i ,j+1) )
2402  call check( __line__, phi_lo(k ,i ,j-1) )
2403 #endif
2404  qmax = max( phi_in(k ,i ,j ), &
2405  phi_in(k+1,i ,j ), &
2406  phi_in(k-1,i ,j ), &
2407  phi_in(k ,i+1,j ), &
2408  phi_in(k ,i-1,j ), &
2409  phi_in(k ,i ,j+1), &
2410  phi_in(k ,i ,j-1), &
2411  phi_lo(k ,i ,j ), &
2412  phi_lo(k+1,i ,j ), &
2413  phi_lo(k-1,i ,j ), &
2414  phi_lo(k ,i+1,j ), &
2415  phi_lo(k ,i-1,j ), &
2416  phi_lo(k ,i ,j+1), &
2417  phi_lo(k ,i ,j-1) )
2418  qmin = min( phi_in(k ,i ,j ), &
2419  phi_in(k+1,i ,j ), &
2420  phi_in(k-1,i ,j ), &
2421  phi_in(k ,i-1,j ), &
2422  phi_in(k ,i+1,j ), &
2423  phi_in(k ,i ,j+1), &
2424  phi_in(k ,i ,j-1), &
2425  phi_lo(k ,i ,j ), &
2426  phi_lo(k+1,i ,j ), &
2427  phi_lo(k-1,i ,j ), &
2428  phi_lo(k ,i-1,j ), &
2429  phi_lo(k ,i+1,j ), &
2430  phi_lo(k ,i ,j+1), &
2431  phi_lo(k ,i ,j-1) )
2432  qjpls(k,i,j) = ( qmax - phi_lo(k,i,j) ) * dens(k,i,j)
2433  qjmns(k,i,j) = ( phi_lo(k,i,j) - qmin ) * dens(k,i,j)
2434  enddo
2435  enddo
2436  enddo
2437 #ifdef DEBUG
2438  k = iundef; i = iundef; j = iundef
2439 #endif
2440  !$omp parallel do private(i,j,qmax,qmin) OMP_SCHEDULE_ collapse(2)
2441  do j = jjs, jje
2442  do i = iis, iie
2443 #ifdef DEBUG
2444  call check( __line__, phi_in(ks ,i ,j ) )
2445  call check( __line__, phi_in(ks+1,i ,j ) )
2446  call check( __line__, phi_in(ks ,i-1,j ) )
2447  call check( __line__, phi_in(ks ,i+1,j ) )
2448  call check( __line__, phi_in(ks ,i ,j+1) )
2449  call check( __line__, phi_in(ks ,i ,j-1) )
2450  call check( __line__, phi_lo(ks ,i ,j ) )
2451  call check( __line__, phi_lo(ks+1,i ,j ) )
2452  call check( __line__, phi_lo(ks ,i-1,j ) )
2453  call check( __line__, phi_lo(ks ,i+1,j ) )
2454  call check( __line__, phi_lo(ks ,i ,j+1) )
2455  call check( __line__, phi_lo(ks ,i ,j-1) )
2456  call check( __line__, phi_in(ke ,i ,j ) )
2457  call check( __line__, phi_in(ke-1,i ,j ) )
2458  call check( __line__, phi_in(ke ,i-1,j ) )
2459  call check( __line__, phi_in(ke ,i+1,j ) )
2460  call check( __line__, phi_in(ke ,i ,j+1) )
2461  call check( __line__, phi_in(ke ,i ,j-1) )
2462  call check( __line__, phi_lo(ke ,i ,j ) )
2463  call check( __line__, phi_lo(ke-1,i ,j ) )
2464  call check( __line__, phi_lo(ke ,i-1,j ) )
2465  call check( __line__, phi_lo(ke ,i+1,j ) )
2466  call check( __line__, phi_lo(ke ,i ,j+1) )
2467  call check( __line__, phi_lo(ke ,i ,j-1) )
2468 #endif
2469  qmax = max( phi_in(ks ,i ,j ), &
2470  phi_in(ks+1,i ,j ), &
2471  phi_in(ks ,i+1,j ), &
2472  phi_in(ks ,i-1,j ), &
2473  phi_in(ks ,i ,j+1), &
2474  phi_in(ks ,i ,j-1), &
2475  phi_lo(ks ,i ,j ), &
2476  phi_lo(ks+1,i ,j ), &
2477  phi_lo(ks ,i+1,j ), &
2478  phi_lo(ks ,i-1,j ), &
2479  phi_lo(ks ,i ,j+1), &
2480  phi_lo(ks ,i ,j-1) )
2481  qmin = min( phi_in(ks ,i ,j ), &
2482  phi_in(ks+1,i ,j ), &
2483  phi_in(ks ,i+1,j ), &
2484  phi_in(ks ,i-1,j ), &
2485  phi_in(ks ,i ,j+1), &
2486  phi_in(ks ,i ,j-1), &
2487  phi_lo(ks ,i ,j ), &
2488  phi_lo(ks+1,i ,j ), &
2489  phi_lo(ks ,i+1,j ), &
2490  phi_lo(ks ,i-1,j ), &
2491  phi_lo(ks ,i ,j+1), &
2492  phi_lo(ks ,i ,j-1) )
2493  qjmns(ks,i,j) = ( phi_lo(ks,i,j) - qmin ) * dens(ks,i,j)
2494  qjpls(ks,i,j) = ( qmax - phi_lo(ks,i,j) ) * dens(ks,i,j)
2495 
2496  qmax = max( phi_in(ke ,i ,j ), &
2497  phi_in(ke-1,i ,j ), &
2498  phi_in(ke ,i+1,j ), &
2499  phi_in(ke ,i-1,j ), &
2500  phi_in(ke ,i ,j+1), &
2501  phi_in(ke ,i ,j-1), &
2502  phi_lo(ke ,i ,j ), &
2503  phi_lo(ke-1,i ,j ), &
2504  phi_lo(ke ,i+1,j ), &
2505  phi_lo(ke ,i-1,j ), &
2506  phi_lo(ke ,i ,j+1), &
2507  phi_lo(ke ,i ,j-1) )
2508  qmin = min( phi_in(ke ,i ,j ), &
2509  phi_in(ke-1,i ,j ), &
2510  phi_in(ke ,i-1,j ), &
2511  phi_in(ke ,i+1,j ), &
2512  phi_in(ke ,i ,j+1), &
2513  phi_in(ke ,i ,j-1), &
2514  phi_lo(ke ,i ,j ), &
2515  phi_lo(ke-1,i ,j ), &
2516  phi_lo(ke ,i-1,j ), &
2517  phi_lo(ke ,i+1,j ), &
2518  phi_lo(ke ,i ,j+1), &
2519  phi_lo(ke ,i ,j-1) )
2520  qjpls(ke,i,j) = ( qmax - phi_lo(ke,i,j) ) * dens(ke,i,j)
2521  qjmns(ke,i,j) = ( phi_lo(ke,i,j) - qmin ) * dens(ke,i,j)
2522  enddo
2523  enddo
2524 #ifdef DEBUG
2525  k = iundef; i = iundef; j = iundef
2526 #endif
2527  end if
2528 
2529  !--- incoming flux limitation factor (0-1)
2530  !$omp parallel do private(i,j,k,zerosw) OMP_SCHEDULE_ collapse(2)
2531  do j = jjs, jje
2532  do i = iis, iie
2533  do k = ks, ke
2534 #ifdef DEBUG
2535  call check( __line__, pjpls(k,i,j) )
2536  call check( __line__, qjpls(k,i,j) )
2537 #endif
2538  ! if pjpls == 0, zerosw = 1 and rjpls = 0
2539  zerosw = 0.5_rp - sign( 0.5_rp, pjpls(k,i,j)-epsilon )
2540  rjpls(k,i,j) = min( 1.0_rp, qjpls(k,i,j) * ( 1.0_rp-zerosw ) / ( pjpls(k,i,j)-zerosw ) )
2541  enddo
2542  enddo
2543  enddo
2544 #ifdef DEBUG
2545  k = iundef; i = iundef; j = iundef
2546 #endif
2547 
2548  !--- outgoing flux limitation factor (0-1)
2549  !$omp parallel do private(i,j,k,zerosw) OMP_SCHEDULE_ collapse(2)
2550  do j = jjs, jje
2551  do i = iis, iie
2552  do k = ks, ke
2553 #ifdef DEBUG
2554  call check( __line__, pjmns(k,i,j) )
2555  call check( __line__, qjmns(k,i,j) )
2556 #endif
2557  ! if pjmns == 0, zerosw = 1 and rjmns = 0
2558  zerosw = 0.5_rp - sign( 0.5_rp, pjmns(k,i,j)-epsilon )
2559  rjmns(k,i,j) = min( 1.0_rp, qjmns(k,i,j) * ( 1.0_rp-zerosw ) / ( pjmns(k,i,j)-zerosw ) )
2560  enddo
2561  enddo
2562  enddo
2563 #ifdef DEBUG
2564  k = iundef; i = iundef; j = iundef
2565 #endif
2566 
2567  enddo
2568  enddo
2569 
2570  call comm_vars8( rjpls(:,:,:), 1 )
2571  call comm_vars8( rjmns(:,:,:), 2 )
2572  call comm_wait ( rjpls(:,:,:), 1 )
2573  call comm_wait ( rjmns(:,:,:), 2 )
2574 
2575  do jjs = js, je, jblock
2576  jje = jjs+jblock-1
2577  do iis = is, ie, iblock
2578  iie = iis+iblock-1
2579 
2580  !--- update high order flux with antidiffusive flux
2581  !$omp parallel do private(i,j,k,dirsw) OMP_SCHEDULE_ collapse(2)
2582  do j = jjs, jje
2583  do i = iis, iie
2584  do k = ks , ke-1
2585 #ifdef DEBUG
2586  call check( __line__, qflx_anti(k,i,j,zdir) )
2587  call check( __line__, rjpls(k ,i,j) )
2588  call check( __line__, rjpls(k+1,i,j) )
2589  call check( __line__, rjmns(k ,i,j) )
2590  call check( __line__, rjmns(k+1,i,j) )
2591 #endif
2592  ! if qflx_anti > 0, dirsw = 1
2593  dirsw = 0.5_rp + sign( 0.5_rp, qflx_anti(k,i,j,zdir) )
2594  qflx_anti(k,i,j,zdir) = qflx_anti(k,i,j,zdir) &
2595  * ( 1.0_rp &
2596  - min( rjpls(k+1,i,j),rjmns(k ,i,j) ) * ( dirsw ) &
2597  - min( rjpls(k ,i,j),rjmns(k+1,i,j) ) * ( 1.0_rp - dirsw ) )
2598  enddo
2599  enddo
2600  enddo
2601 #ifdef DEBUG
2602  k = iundef; i = iundef; j = iundef
2603 #endif
2604  !$omp parallel do private(i,j) OMP_SCHEDULE_ collapse(2)
2605  do j = jjs, jje
2606  do i = iis, iie
2607 #ifdef DEBUG
2608  call check( __line__, qflx_anti(ke,i,j,zdir) )
2609  call check( __line__, rjpls(ke ,i,j) )
2610  call check( __line__, rjmns(ke ,i,j) )
2611 #endif
2612  qflx_anti(ks-1,i,j,zdir) = 0.0_rp ! bottom boundary
2613  qflx_anti(ke ,i,j,zdir) = 0.0_rp ! top boundary
2614  enddo
2615  enddo
2616 #ifdef DEBUG
2617  k = iundef; i = iundef; j = iundef
2618 #endif
2619 
2620  if ( iis == is ) then
2621  ijs = iis-1
2622  else
2623  ijs = iis
2624  end if
2625 
2626  !$omp parallel do private(i,j,k,dirsw) OMP_SCHEDULE_ collapse(2)
2627  do j = jjs, jje
2628  do i = ijs, iie
2629  do k = ks, ke
2630 #ifdef DEBUG
2631  call check( __line__, qflx_anti(k,i,j,xdir) )
2632  call check( __line__, rjpls(k,i ,j) )
2633  call check( __line__, rjpls(k,i+1,j) )
2634  call check( __line__, rjmns(k,i ,j) )
2635  call check( __line__, rjmns(k,i+1,j) )
2636 #endif
2637  ! if qflx_anti > 0, dirsw = 1
2638  dirsw = 0.5_rp + sign( 0.5_rp, qflx_anti(k,i,j,xdir) )
2639  qflx_anti(k,i,j,xdir) = qflx_anti(k,i,j,xdir) &
2640  * ( 1.0_rp &
2641  - min( rjpls(k,i+1,j),rjmns(k,i ,j) ) * ( dirsw ) &
2642  - min( rjpls(k,i ,j),rjmns(k,i+1,j) ) * ( 1.0_rp - dirsw ) )
2643  enddo
2644  enddo
2645  enddo
2646 #ifdef DEBUG
2647  k = iundef; i = iundef; j = iundef
2648 #endif
2649 
2650  if ( jjs == js ) then
2651  ijs = jjs-1
2652  else
2653  ijs = jjs
2654  end if
2655  !$omp parallel do private(i,j,k,dirsw) OMP_SCHEDULE_ collapse(2)
2656  do j = ijs, jje
2657  do i = iis, iie
2658  do k = ks, ke
2659 #ifdef DEBUG
2660  call check( __line__, qflx_anti(k,i,j,ydir) )
2661  call check( __line__, rjpls(k,i,j+1) )
2662  call check( __line__, rjpls(k,i,j ) )
2663  call check( __line__, rjmns(k,i,j ) )
2664  call check( __line__, rjmns(k,i,j+1) )
2665 #endif
2666  ! if qflx_anti > 0, dirsw = 1
2667  dirsw = 0.5_rp + sign( 0.5_rp, qflx_anti(k,i,j,ydir) )
2668  qflx_anti(k,i,j,ydir) = qflx_anti(k,i,j,ydir) &
2669  * ( 1.0_rp &
2670  - min( rjpls(k,i,j+1),rjmns(k,i,j ) ) * ( dirsw ) &
2671  - min( rjpls(k,i,j ),rjmns(k,i,j+1) ) * ( 1.0_rp - dirsw ) )
2672  enddo
2673  enddo
2674  enddo
2675 #ifdef DEBUG
2676  k = iundef; i = iundef; j = iundef
2677 #endif
2678 
2679  enddo
2680  enddo
2681 
2682  return
2683  end subroutine atmos_dyn_fct
2684 
2685  !-----------------------------------------------------------------------------
2686  ! private procedure
2687  ! get factor for FCT
2688  subroutine get_fact_fct( &
2689  fact, &
2690  rw, ru, rv )
2691  use scale_const, only: &
2692  epsilon => const_eps
2693  implicit none
2694  real(RP), intent(out) :: fact(0:1,-1:1,-1:1)
2695  real(RP), intent(in) :: rw, ru, rv
2696 
2697  real(RP) :: sign_uv, sign_uw, sign_vw ! uv>=0, uw>=0, vw>=0
2698  real(RP) :: ugev, ugew, vgew ! u>=v, u>=w, u>=w
2699  real(RP) :: umax, vmax, wmax
2700  real(RP) :: vu, wu, uv, wv, uw, vw ! |v/u|, |w/u|, ....
2701  real(RP) :: uzero, vzero, wzero
2702  !---------------------------------------------------------------------------
2703 
2704  ugev = sign(0.5_rp, abs(ru)-abs(rv)) + 0.5_rp ! u >= v
2705  ugew = sign(0.5_rp, abs(ru)-abs(rw)) + 0.5_rp ! u >= w
2706  vgew = sign(0.5_rp, abs(rv)-abs(rw)) + 0.5_rp ! v >= w
2707 
2708  uzero = sign(0.5_rp,abs(ru)-epsilon) - 0.5_rp
2709  vzero = sign(0.5_rp,abs(rv)-epsilon) - 0.5_rp
2710  wzero = sign(0.5_rp,abs(rw)-epsilon) - 0.5_rp
2711 
2712  sign_uv = sign(0.5_rp, ru*rv) + 0.5_rp ! uv >= 0
2713  sign_uw = sign(0.5_rp, ru*rw) + 0.5_rp ! uw >= 0
2714  sign_vw = sign(0.5_rp, rv*rw) + 0.5_rp ! vw >= 0
2715 
2716  wu = abs( rw / ( ru+uzero ) * ( 1.0_rp+uzero ) )
2717  vu = abs( rv / ( ru+uzero ) * ( 1.0_rp+uzero ) )
2718  uv = abs( ru / ( rv+vzero ) * ( 1.0_rp+vzero ) )
2719  wv = abs( rw / ( rv+vzero ) * ( 1.0_rp+vzero ) )
2720  uw = abs( ru / ( rw+wzero ) * ( 1.0_rp+wzero ) )
2721  vw = abs( rv / ( rw+wzero ) * ( 1.0_rp+wzero ) )
2722 
2723  umax = ugev * ugew * ( 1.0_rp+uzero ) ! u == max(u,v,w)
2724  vmax = (1.0_rp-ugev) * vgew ! v == max(u,v,w)
2725  wmax = 1.0_rp - ugev * ugew - vmax ! w == max(u,v,w)
2726 
2727  fact(0, 0, 0) = - ugev * ugew * uzero ! 1.0 if max(u,v,w) < epsilon
2728 
2729  fact(1, 0, 0) = wmax * (1.0_rp-uw) * (1.0_rp-vw)
2730  fact(0, 1, 0) = umax * (1.0_rp-vu) * (1.0_rp-wu)
2731  fact(0, 0, 1) = vmax * (1.0_rp-uv) * (1.0_rp-wv)
2732 
2733  fact(1, 1, 1) = sign_uv * sign_uw * ( umax * vu*wu + vmax * uv*wv + wmax * uw*vw )
2734  fact(1,-1, 1) = (1.0_rp-sign_uv) * (1.0_rp-sign_uw) * ( umax * vu*wu + vmax * uv*wv + wmax * uw*vw )
2735  fact(1, 1,-1) = (1.0_rp-sign_uv) * sign_uw * ( umax * vu*wu + vmax * uv*wv + wmax * uw*vw )
2736  fact(1,-1,-1) = sign_uv * (1.0_rp-sign_uw) * ( umax * vu*wu + vmax * uv*wv + wmax * uw*vw )
2737 
2738  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) )
2739  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) )
2740  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) )
2741  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) )
2742  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) )
2743  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) )
2744 
2745  return
2746  end subroutine get_fact_fct
2747 
2748 end module scale_atmos_dyn_common
scale_index::i_rhot
integer, parameter, public i_rhot
Definition: scale_index.F90:35
scale_grid_index::is
integer, public is
start point of inner domain: x, local
Definition: scale_grid_index.F90:60
scale_debug
module DEBUG
Definition: scale_debug.F90:13
scale_gridtrans::i_xvz
integer, public i_xvz
Definition: scale_gridtrans.F90:54
scale_grid_index::je
integer, public je
end point of inner domain: y, local
Definition: scale_grid_index.F90:63
scale_process::prc_mpistop
subroutine, public prc_mpistop
Abort MPI.
Definition: scale_process.F90:267
scale_atmos_dyn_common::atmos_dyn_numfilter_coef_q
subroutine, public atmos_dyn_numfilter_coef_q(num_diff_q, DENS, QTRC, CDZ, CDX, CDY, dt, REF_qv, iq, ND_COEF, ND_ORDER, ND_SFC_FACT, ND_USE_RS)
Calc coefficient of numerical filter.
Definition: scale_atmos_dyn_common.F90:998
scale_index::va
integer, public va
Definition: scale_index.F90:38
scale_comm::comm_vars8_init
subroutine, public comm_vars8_init(varname, var, vid)
Register variables.
Definition: scale_comm.F90:313
scale_atmos_dyn_common::atmos_dyn_numfilter_coef
subroutine, public atmos_dyn_numfilter_coef(num_diff, DENS, MOMZ, MOMX, MOMY, RHOT, CDZ, CDX, CDY, FDZ, FDX, FDY, DT, REF_dens, REF_pott, ND_COEF, ND_ORDER, ND_SFC_FACT, ND_USE_RS)
Calc coefficient of numerical filter.
Definition: scale_atmos_dyn_common.F90:394
scale_grid_index::iblock
integer, public iblock
block size for cache blocking: x
Definition: scale_grid_index.F90:47
scale_index::i_momx
integer, parameter, public i_momx
Definition: scale_index.F90:33
scale_stdio::io_l
logical, public io_l
output log or not? (this process)
Definition: scale_stdio.F90:61
scale_grid_index::zdir
integer, parameter, public zdir
Definition: scale_grid_index.F90:35
scale_index::i_momz
integer, parameter, public i_momz
Definition: scale_index.F90:32
scale_stdio
module STDIO
Definition: scale_stdio.F90:12
scale_grid_index::ydir
integer, parameter, public ydir
Definition: scale_grid_index.F90:37
scale_grid_index::ke
integer, public ke
end point of inner domain: z, local
Definition: scale_grid_index.F90:59
scale_grid_index::xdir
integer, parameter, public xdir
Definition: scale_grid_index.F90:36
scale_atmos_dyn_common::atmos_dyn_filter_tend
subroutine, public atmos_dyn_filter_tend(phi_t, phi, rdz, rdx, rdy, KO, IO, JO)
Definition: scale_atmos_dyn_common.F90:1212
scale_gridtrans::i_xy
integer, public i_xy
Definition: scale_gridtrans.F90:57
scale_debug::check
subroutine, public check(current_line, v)
Undefined value checker.
Definition: scale_debug.F90:58
scale_index::i_dens
integer, parameter, public i_dens
Definition: scale_index.F90:31
scale_index::i_momy
integer, parameter, public i_momy
Definition: scale_index.F90:34
scale_atmos_dyn_common::atmos_dyn_wdamp_setup
subroutine, public atmos_dyn_wdamp_setup(wdamp_coef, wdamp_tau, wdamp_height, FZ)
Setup.
Definition: scale_atmos_dyn_common.F90:330
scale_const::const_undef
real(rp), public const_undef
Definition: scale_const.F90:43
scale_grid_index
module grid index
Definition: scale_grid_index.F90:14
scale_tracer
module TRACER
Definition: scale_tracer.F90:15
scale_index
module Index
Definition: scale_index.F90:14
scale_grid_index::ia
integer, public ia
of whole cells: x, local, with HALO
Definition: scale_grid_index.F90:55
scale_atmos_hydrometeor
module Hydrometeor
Definition: scale_atmos_sub_hydrometeor.F90:14
scale_gridtrans
module GRIDTRANS
Definition: scale_gridtrans.F90:15
scale_grid_index::ka
integer, public ka
of whole cells: z, local, with HALO
Definition: scale_grid_index.F90:54
scale_gridtrans::i_uy
integer, public i_uy
Definition: scale_gridtrans.F90:58
scale_atmos_dyn_common::calc_diff3
subroutine calc_diff3(diff, phi, KO, IO, JO)
Definition: scale_atmos_dyn_common.F90:1608
scale_grid_index::jblock
integer, public jblock
block size for cache blocking: y
Definition: scale_grid_index.F90:48
scale_gridtrans::i_xyw
integer, public i_xyw
Definition: scale_gridtrans.F90:50
scale_atmos_dyn_common::atmos_dyn_filter_setup
subroutine, public atmos_dyn_filter_setup(num_diff, num_diff_q, CDZ, CDX, CDY, FDZ, FDX, FDY)
Setup.
Definition: scale_atmos_dyn_common.F90:98
scale_grid_index::jhalo
integer, public jhalo
of halo cells: y
Definition: scale_grid_index.F90:52
scale_comm
module COMMUNICATION
Definition: scale_comm.F90:23
scale_grid_index::js
integer, public js
start point of inner domain: y, local
Definition: scale_grid_index.F90:62
scale_const::const_undef2
integer, parameter, public const_undef2
undefined value (INT2)
Definition: scale_const.F90:40
scale_atmos_dyn_common
module Atmosphere / Dynamics common
Definition: scale_atmos_dyn_common.F90:15
scale_process
module PROCESS
Definition: scale_process.F90:14
scale_const
module CONSTANT
Definition: scale_const.F90:14
scale_grid_index::ks
integer, public ks
start point of inner domain: z, local
Definition: scale_grid_index.F90:58
scale_grid_index::khalo
integer, parameter, public khalo
of halo cells: z
Definition: scale_grid_index.F90:50
scale_atmos_dyn_common::calc_numdiff
subroutine calc_numdiff(work, iwork, data, nd_order, KO, IO, JO, KEE)
Definition: scale_atmos_dyn_common.F90:1535
scale_prof::prof_rapstart
subroutine, public prof_rapstart(rapname_base, level)
Start raptime.
Definition: scale_prof.F90:156
scale_atmos_dyn_common::atmos_dyn_copy_boundary
subroutine, public atmos_dyn_copy_boundary(DENS, MOMZ, MOMX, MOMY, RHOT, PROG, DENS0, MOMZ0, MOMX0, MOMY0, RHOT0, PROG0, BND_W, BND_E, BND_S, BND_N)
Definition: scale_atmos_dyn_common.F90:1256
scale_gridtrans::i_xv
integer, public i_xv
Definition: scale_gridtrans.F90:59
scale_prof
module profiler
Definition: scale_prof.F90:10
scale_grid_index::ie
integer, public ie
end point of inner domain: x, local
Definition: scale_grid_index.F90:61
scale_const::const_eps
real(rp), public const_eps
small number
Definition: scale_const.F90:36
scale_atmos_hydrometeor::i_qv
integer, public i_qv
Definition: scale_atmos_sub_hydrometeor.F90:68
scale_gridtrans::i_uyz
integer, public i_uyz
Definition: scale_gridtrans.F90:53
scale_precision
module PRECISION
Definition: scale_precision.F90:13
scale_atmos_dyn_common::atmos_dyn_fct
subroutine, public atmos_dyn_fct(qflx_anti, phi_in, DENS0, DENS, qflx_hi, qflx_lo, mflx_hi, rdz, rdx, rdy, GSQRT, MAPF, dt, flag_vect)
Flux Correction Transport Limiter.
Definition: scale_atmos_dyn_common.F90:1942
scale_atmos_dyn_common::atmos_dyn_divergence
subroutine, public atmos_dyn_divergence(DDIV, MOMZ, MOMX, MOMY, GSQRT, J13G, J23G, J33G, MAPF, RCDZ, RCDX, RCDY, RFDZ, FDZ)
Definition: scale_atmos_dyn_common.F90:1446
scale_const::const_pi
real(rp), public const_pi
pi
Definition: scale_const.F90:34
scale_gridtrans::i_xyz
integer, public i_xyz
Definition: scale_gridtrans.F90:49
scale_stdio::io_fid_log
integer, public io_fid_log
Log file ID.
Definition: scale_stdio.F90:56
scale_prof::prof_rapend
subroutine, public prof_rapend(rapname_base, level)
Save raptime.
Definition: scale_prof.F90:204
scale_atmos_dyn_common::atmos_dyn_copy_boundary_tracer
subroutine, public atmos_dyn_copy_boundary_tracer(QTRC, QTRC0, BND_W, BND_E, BND_S, BND_N)
Definition: scale_atmos_dyn_common.F90:1378
scale_grid_index::ihalo
integer, public ihalo
of halo cells: x
Definition: scale_grid_index.F90:51
scale_grid_index::ja
integer, public ja
of whole cells: y, local, with HALO
Definition: scale_grid_index.F90:56
Generated by   doxygen 1.8.14