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