FV3 Bundle
a2b_edge_nlm.F90
Go to the documentation of this file.
1 !***********************************************************************
2 !* GNU General Public License *
3 !* This file is a part of fvGFS. *
4 !* *
5 !* fvGFS is free software; you can redistribute it and/or modify it *
6 !* and are expected to follow the terms of the GNU General Public *
7 !* License as published by the Free Software Foundation; either *
8 !* version 2 of the License, or (at your option) any later version. *
9 !* *
10 !* fvGFS is distributed in the hope that it will be useful, but *
11 !* WITHOUT ANY WARRANTY; without even the implied warranty of *
12 !* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
13 !* General Public License for more details. *
14 !* *
15 !* For the full text of the GNU General Public License, *
16 !* write to: Free Software Foundation, Inc., *
17 !* 675 Mass Ave, Cambridge, MA 02139, USA. *
18 !* or see: http://www.gnu.org/licenses/gpl.html *
19 !***********************************************************************
21 
23 #ifdef VAN2
24  use fv_grid_utils_nlm_mod, only: van2
25 #endif
26 
28 
29  implicit none
30 
31  real, parameter:: r3 = 1./3.
32 !----------------------------
33 ! 4-pt Lagrange interpolation
34 !----------------------------
35  real, parameter:: a1 = 0.5625 ! 9/16
36  real, parameter:: a2 = -0.0625 ! -1/16
37 !----------------------
38 ! PPM volume mean form:
39 !----------------------
40  real, parameter:: b1 = 7./12. ! 0.58333333
41  real, parameter:: b2 = -1./12.
42 
43  private
44  public :: a2b_ord2, a2b_ord4
45 
46 contains
47 
48 #ifndef USE_OLD_ALGORITHM
49  subroutine a2b_ord4(qin, qout, gridstruct, npx, npy, is, ie, js, je, ng, replace)
50  integer, intent(IN):: npx, npy, is, ie, js, je, ng
51  real, intent(INOUT):: qin(is-ng:ie+ng,js-ng:je+ng) ! A-grid field
52  real, intent(INOUT):: qout(is-ng:ie+ng,js-ng:je+ng) ! Output B-grid field
53  type(fv_grid_type), intent(IN), target :: gridstruct
54  logical, optional, intent(IN):: replace
55 ! local: compact 4-pt cubic
56  real, parameter:: c1 = 2./3.
57  real, parameter:: c2 = -1./6.
58 ! Parabolic spline
59 ! real, parameter:: c1 = 0.75
60 ! real, parameter:: c2 = -0.25
61 
62  real qx(is:ie+1,js-ng:je+ng)
63  real qy(is-ng:ie+ng,js:je+1)
64  real qxx(is-ng:ie+ng,js-ng:je+ng)
65  real qyy(is-ng:ie+ng,js-ng:je+ng)
66  real g_in, g_ou
67  real:: p0(2)
68  real:: q1(is-1:ie+1), q2(js-1:je+1)
69  integer:: i, j, is1, js1, is2, js2, ie1, je1
70 
71  real, pointer, dimension(:,:,:) :: grid, agrid
72  real, pointer, dimension(:,:) :: dxa, dya
73  real(kind=R_GRID), pointer, dimension(:) :: edge_w, edge_e, edge_s, edge_n
74 
75  edge_w => gridstruct%edge_w
76  edge_e => gridstruct%edge_e
77  edge_s => gridstruct%edge_s
78  edge_n => gridstruct%edge_n
79 
80  grid => gridstruct%grid
81  agrid => gridstruct%agrid
82  dxa => gridstruct%dxa
83  dya => gridstruct%dya
84 
85  if (gridstruct%grid_type < 3) then
86 
87  is1 = max(1,is-1)
88  js1 = max(1,js-1)
89  is2 = max(2,is)
90  js2 = max(2,js)
91 
92  ie1 = min(npx-1,ie+1)
93  je1 = min(npy-1,je+1)
94 
95 ! Corners:
96 ! 3-way extrapolation
97  if (gridstruct%nested) then
98 
99  do j=js-2,je+2
100  do i=is,ie+1
101  qx(i,j) = b2*(qin(i-2,j)+qin(i+1,j)) + b1*(qin(i-1,j)+qin(i,j))
102  enddo
103  enddo
104 
105 
106  else
107 
108  if ( gridstruct%sw_corner ) then
109  p0(1:2) = grid(1,1,1:2)
110  qout(1,1) = (extrap_corner(p0, agrid(1,1,1:2), agrid( 2, 2,1:2), qin(1,1), qin( 2, 2)) + &
111  extrap_corner(p0, agrid(0,1,1:2), agrid(-1, 2,1:2), qin(0,1), qin(-1, 2)) + &
112  extrap_corner(p0, agrid(1,0,1:2), agrid( 2,-1,1:2), qin(1,0), qin( 2,-1)))*r3
113 
114  endif
115  if ( gridstruct%se_corner ) then
116  p0(1:2) = grid(npx,1,1:2)
117  qout(npx,1) = (extrap_corner(p0, agrid(npx-1,1,1:2), agrid(npx-2, 2,1:2), qin(npx-1,1), qin(npx-2, 2)) + &
118  extrap_corner(p0, agrid(npx-1,0,1:2), agrid(npx-2,-1,1:2), qin(npx-1,0), qin(npx-2,-1)) + &
119  extrap_corner(p0, agrid(npx ,1,1:2), agrid(npx+1, 2,1:2), qin(npx ,1), qin(npx+1, 2)))*r3
120  endif
121  if ( gridstruct%ne_corner ) then
122  p0(1:2) = grid(npx,npy,1:2)
123  qout(npx,npy) = (extrap_corner(p0, agrid(npx-1,npy-1,1:2), agrid(npx-2,npy-2,1:2), qin(npx-1,npy-1), qin(npx-2,npy-2)) + &
124  extrap_corner(p0, agrid(npx ,npy-1,1:2), agrid(npx+1,npy-2,1:2), qin(npx ,npy-1), qin(npx+1,npy-2)) + &
125  extrap_corner(p0, agrid(npx-1,npy ,1:2), agrid(npx-2,npy+1,1:2), qin(npx-1,npy ), qin(npx-2,npy+1)))*r3
126  endif
127  if ( gridstruct%nw_corner ) then
128  p0(1:2) = grid(1,npy,1:2)
129  qout(1,npy) = (extrap_corner(p0, agrid(1,npy-1,1:2), agrid( 2,npy-2,1:2), qin(1,npy-1), qin( 2,npy-2)) + &
130  extrap_corner(p0, agrid(0,npy-1,1:2), agrid(-1,npy-2,1:2), qin(0,npy-1), qin(-1,npy-2)) + &
131  extrap_corner(p0, agrid(1,npy, 1:2), agrid( 2,npy+1,1:2), qin(1,npy ), qin( 2,npy+1)))*r3
132  endif
133 
134 !------------
135 ! X-Interior:
136 !------------
137  do j=max(1,js-2),min(npy-1,je+2)
138  do i=max(3,is), min(npx-2,ie+1)
139  qx(i,j) = b2*(qin(i-2,j)+qin(i+1,j)) + b1*(qin(i-1,j)+qin(i,j))
140  enddo
141  enddo
142 
143  ! *** West Edges:
144  if ( is==1 ) then
145  do j=js1, je1
146  q2(j) = (qin(0,j)*dxa(1,j) + qin(1,j)*dxa(0,j))/(dxa(0,j) + dxa(1,j))
147  enddo
148  do j=js2, je1
149  qout(1,j) = edge_w(j)*q2(j-1) + (1.-edge_w(j))*q2(j)
150  enddo
151 !
152  do j=max(1,js-2),min(npy-1,je+2)
153  g_in = dxa(2,j) / dxa(1,j)
154  g_ou = dxa(-1,j) / dxa(0,j)
155  qx(1,j) = 0.5*( ((2.+g_in)*qin(1,j)-qin( 2,j))/(1.+g_in) + &
156  ((2.+g_ou)*qin(0,j)-qin(-1,j))/(1.+g_ou) )
157  qx(2,j) = ( 3.*(g_in*qin(1,j)+qin(2,j))-(g_in*qx(1,j)+qx(3,j)) ) / (2.+2.*g_in)
158  enddo
159  endif
160 
161  ! East Edges:
162  if ( (ie+1)==npx ) then
163  do j=js1, je1
164  q2(j) = (qin(npx-1,j)*dxa(npx,j) + qin(npx,j)*dxa(npx-1,j))/(dxa(npx-1,j) + dxa(npx,j))
165  enddo
166  do j=js2, je1
167  qout(npx,j) = edge_e(j)*q2(j-1) + (1.-edge_e(j))*q2(j)
168  enddo
169 !
170  do j=max(1,js-2),min(npy-1,je+2)
171  g_in = dxa(npx-2,j) / dxa(npx-1,j)
172  g_ou = dxa(npx+1,j) / dxa(npx,j)
173  qx(npx,j) = 0.5*( ((2.+g_in)*qin(npx-1,j)-qin(npx-2,j))/(1.+g_in) + &
174  ((2.+g_ou)*qin(npx, j)-qin(npx+1,j))/(1.+g_ou) )
175  qx(npx-1,j) = (3.*(qin(npx-2,j)+g_in*qin(npx-1,j)) - (g_in*qx(npx,j)+qx(npx-2,j)))/(2.+2.*g_in)
176  enddo
177  endif
178 
179  end if
180 !------------
181 ! Y-Interior:
182 !------------
183 
184  if (gridstruct%nested) then
185 
186 
187  do j=js,je+1
188  do i=is-2,ie+2
189  qy(i,j) = b2*(qin(i,j-2)+qin(i,j+1)) + b1*(qin(i,j-1) + qin(i,j))
190  enddo
191  enddo
192 
193  else
194 
195  do j=max(3,js),min(npy-2,je+1)
196  do i=max(1,is-2), min(npx-1,ie+2)
197  qy(i,j) = b2*(qin(i,j-2)+qin(i,j+1)) + b1*(qin(i,j-1) + qin(i,j))
198  enddo
199  enddo
200 
201  ! South Edges:
202  if ( js==1 ) then
203  do i=is1, ie1
204  q1(i) = (qin(i,0)*dya(i,1) + qin(i,1)*dya(i,0))/(dya(i,0) + dya(i,1))
205  enddo
206  do i=is2, ie1
207  qout(i,1) = edge_s(i)*q1(i-1) + (1.-edge_s(i))*q1(i)
208  enddo
209 !
210  do i=max(1,is-2),min(npx-1,ie+2)
211  g_in = dya(i,2) / dya(i,1)
212  g_ou = dya(i,-1) / dya(i,0)
213  qy(i,1) = 0.5*( ((2.+g_in)*qin(i,1)-qin(i,2))/(1.+g_in) + &
214  ((2.+g_ou)*qin(i,0)-qin(i,-1))/(1.+g_ou) )
215  qy(i,2) = (3.*(g_in*qin(i,1)+qin(i,2)) - (g_in*qy(i,1)+qy(i,3)))/(2.+2.*g_in)
216  enddo
217  endif
218 
219  ! North Edges:
220  if ( (je+1)==npy ) then
221  do i=is1, ie1
222  q1(i) = (qin(i,npy-1)*dya(i,npy) + qin(i,npy)*dya(i,npy-1))/(dya(i,npy-1)+dya(i,npy))
223  enddo
224  do i=is2, ie1
225  qout(i,npy) = edge_n(i)*q1(i-1) + (1.-edge_n(i))*q1(i)
226  enddo
227 !
228  do i=max(1,is-2),min(npx-1,ie+2)
229  g_in = dya(i,npy-2) / dya(i,npy-1)
230  g_ou = dya(i,npy+1) / dya(i,npy)
231  qy(i,npy) = 0.5*( ((2.+g_in)*qin(i,npy-1)-qin(i,npy-2))/(1.+g_in) + &
232  ((2.+g_ou)*qin(i,npy )-qin(i,npy+1))/(1.+g_ou) )
233  qy(i,npy-1) = (3.*(qin(i,npy-2)+g_in*qin(i,npy-1)) - (g_in*qy(i,npy)+qy(i,npy-2)))/(2.+2.*g_in)
234  enddo
235  endif
236 
237  end if
238 !--------------------------------------
239 
240  if (gridstruct%nested) then
241 
242  do j=js, je+1
243  do i=is,ie+1
244  qxx(i,j) = a2*(qx(i,j-2)+qx(i,j+1)) + a1*(qx(i,j-1)+qx(i,j))
245  enddo
246  enddo
247 
248  do j=js,je+1
249  do i=is,ie+1
250  qyy(i,j) = a2*(qy(i-2,j)+qy(i+1,j)) + a1*(qy(i-1,j)+qy(i,j))
251  enddo
252 
253  do i=is,ie+1
254  qout(i,j) = 0.5*(qxx(i,j) + qyy(i,j)) ! averaging
255  enddo
256  enddo
257 
258 
259 
260  else
261 
262  do j=max(3,js),min(npy-2,je+1)
263  do i=max(2,is),min(npx-1,ie+1)
264  qxx(i,j) = a2*(qx(i,j-2)+qx(i,j+1)) + a1*(qx(i,j-1)+qx(i,j))
265  enddo
266  enddo
267 
268  if ( js==1 ) then
269  do i=max(2,is),min(npx-1,ie+1)
270  qxx(i,2) = c1*(qx(i,1)+qx(i,2))+c2*(qout(i,1)+qxx(i,3))
271  enddo
272  endif
273  if ( (je+1)==npy ) then
274  do i=max(2,is),min(npx-1,ie+1)
275  qxx(i,npy-1) = c1*(qx(i,npy-2)+qx(i,npy-1))+c2*(qout(i,npy)+qxx(i,npy-2))
276  enddo
277  endif
278 
279 
280  do j=max(2,js),min(npy-1,je+1)
281  do i=max(3,is),min(npx-2,ie+1)
282  qyy(i,j) = a2*(qy(i-2,j)+qy(i+1,j)) + a1*(qy(i-1,j)+qy(i,j))
283  enddo
284  if ( is==1 ) qyy(2,j) = c1*(qy(1,j)+qy(2,j))+c2*(qout(1,j)+qyy(3,j))
285  if((ie+1)==npx) qyy(npx-1,j) = c1*(qy(npx-2,j)+qy(npx-1,j))+c2*(qout(npx,j)+qyy(npx-2,j))
286 
287  do i=max(2,is),min(npx-1,ie+1)
288  qout(i,j) = 0.5*(qxx(i,j) + qyy(i,j)) ! averaging
289  enddo
290  enddo
291 
292  end if
293 
294  else ! grid_type>=3
295 !------------------------
296 ! Doubly periodic domain:
297 !------------------------
298 ! X-sweep: PPM
299  do j=js-2,je+2
300  do i=is,ie+1
301  qx(i,j) = b1*(qin(i-1,j)+qin(i,j)) + b2*(qin(i-2,j)+qin(i+1,j))
302  enddo
303  enddo
304 ! Y-sweep: PPM
305  do j=js,je+1
306  do i=is-2,ie+2
307  qy(i,j) = b1*(qin(i,j-1)+qin(i,j)) + b2*(qin(i,j-2)+qin(i,j+1))
308  enddo
309  enddo
310 
311  do j=js,je+1
312  do i=is,ie+1
313  qout(i,j) = 0.5*( a1*(qx(i,j-1)+qx(i,j ) + qy(i-1,j)+qy(i, j)) + &
314  a2*(qx(i,j-2)+qx(i,j+1) + qy(i-2,j)+qy(i+1,j)) )
315  enddo
316  enddo
317  endif
318 
319  if ( present(replace) ) then
320  if ( replace ) then
321  do j=js,je+1
322  do i=is,ie+1
323  qin(i,j) = qout(i,j)
324  enddo
325  enddo
326  endif
327  endif
328 
329  end subroutine a2b_ord4
330 
331 #else
332 
333 ! Working version:
334  subroutine a2b_ord4(qin, qout, gridstruct, npx, npy, is, ie, js, je, ng, replace)
335  integer, intent(IN):: npx, npy, is, ie, js, je, ng
336  real, intent(INOUT):: qin(is-ng:ie+ng,js-ng:je+ng) ! A-grid field
337  real, intent(INOUT):: qout(is-ng:ie+ng,js-ng:je+ng) ! Output B-grid field
338  type(fv_grid_type), intent(IN), target :: gridstruct
339  logical, optional, intent(IN):: replace
340 ! local: compact 4-pt cubic
341  real, parameter:: c1 = 2./3.
342  real, parameter:: c2 = -1./6.
343 ! Parabolic spline
344 ! real, parameter:: c1 = 0.75
345 ! real, parameter:: c2 = -0.25
346 !-----------------------------
347 ! 6-pt corner interpolation:
348 !-----------------------------
349  real, parameter:: d1 = 0.375 ! 0.5
350  real, parameter:: d2 = -1./24. ! -1./6.
351 
352  real qx(is:ie+1,js-ng:je+ng)
353  real qy(is-ng:ie+ng,js:je+1)
354  real qxx(is-ng:ie+ng,js-ng:je+ng)
355  real qyy(is-ng:ie+ng,js-ng:je+ng)
356  real gratio
357  real g_in, g_ou
358  real:: p0(2)
359  real q1(npx), q2(npy)
360  integer:: i, j, is1, js1, is2, js2, ie1, je1
361 
362  real, pointer, dimension(:,:,:) :: grid, agrid
363  real, pointer, dimension(:,:) :: dxa, dya
364  real, pointer, dimension(:) :: edge_w, edge_e, edge_s, edge_n
365 
366  edge_w => gridstruct%edge_w
367  edge_e => gridstruct%edge_e
368  edge_s => gridstruct%edge_s
369  edge_n => gridstruct%edge_n
370 
371 
372  grid => gridstruct%grid
373  agrid => gridstruct%agrid
374  dxa => gridstruct%dxa
375  dya => gridstruct%dya
376 
377  if (gridstruct%grid_type < 3) then
378 
379  is1 = max(1,is-1)
380  js1 = max(1,js-1)
381  is2 = max(2,is)
382  js2 = max(2,js)
383 
384  ie1 = min(npx-1,ie+1)
385  je1 = min(npy-1,je+1)
386 
387 ! Corners:
388 #ifdef USE_3PT
389  if ( gridstruct%sw_corner ) qout(1, 1) = r3*(qin(1, 1)+qin(1, 0)+qin(0, 1))
390  if ( gridstruct%se_corner ) qout(npx, 1) = r3*(qin(npx-1, 1)+qin(npx-1, 0)+qin(npx, 1))
391  if ( gridstruct%ne_corner ) qout(npx,npy) = r3*(qin(npx-1,npy-1)+qin(npx,npy-1)+qin(npx-1,npy))
392  if ( gridstruct%nw_corner ) qout(1, npy) = r3*(qin(1, npy-1)+qin(0, npy-1)+qin(1, npy))
393 #else
394 
395 #ifdef SYMM_GRID
396 ! Symmetrical 6-point formular:
397  if ( gridstruct%sw_corner ) then
398  qout(1,1) = d1*(qin(1, 0) + qin( 0,1) + qin(1,1)) + &
399  d2*(qin(2,-1) + qin(-1,2) + qin(2,2))
400  endif
401  if ( gridstruct%se_corner ) then
402  qout(npx,1) = d1*(qin(npx-1, 0) + qin(npx-1,1) + qin(npx, 1)) + &
403  d2*(qin(npx-2,-1) + qin(npx-2,2) + qin(npx+1,2))
404  endif
405  if ( gridstruct%ne_corner ) then
406  qout(npx,npy) = d1*(qin(npx-1,npy-1) + qin(npx, npy-1) + qin(npx-1,npy)) + &
407  d2*(qin(npx-2,npy-2) + qin(npx+1,npy-2) + qin(npx-2,npy+1))
408  endif
409  if ( gridstruct%nw_corner ) then
410  qout(1,npy) = d1*(qin( 0,npy-1) + qin(1,npy-1) + qin(1,npy)) + &
411  d2*(qin(-1,npy-2) + qin(2,npy-2) + qin(2,npy+1))
412  endif
413 #else
414 ! 3-way extrapolation
415  if ( gridstruct%sw_corner ) then
416  p0(1:2) = grid(1,1,1:2)
417  qout(1,1) = (extrap_corner(p0, agrid(1,1,1:2), agrid( 2, 2,1:2), qin(1,1), qin( 2, 2)) + &
418  extrap_corner(p0, agrid(0,1,1:2), agrid(-1, 2,1:2), qin(0,1), qin(-1, 2)) + &
419  extrap_corner(p0, agrid(1,0,1:2), agrid( 2,-1,1:2), qin(1,0), qin( 2,-1)))*r3
420 
421  endif
422  if ( gridstruct%se_corner ) then
423  p0(1:2) = grid(npx,1,1:2)
424  qout(npx,1) = (extrap_corner(p0, agrid(npx-1,1,1:2), agrid(npx-2, 2,1:2), qin(npx-1,1), qin(npx-2, 2)) + &
425  extrap_corner(p0, agrid(npx-1,0,1:2), agrid(npx-2,-1,1:2), qin(npx-1,0), qin(npx-2,-1)) + &
426  extrap_corner(p0, agrid(npx ,1,1:2), agrid(npx+1, 2,1:2), qin(npx ,1), qin(npx+1, 2)))*r3
427  endif
428  if ( gridstruct%ne_corner ) then
429  p0(1:2) = grid(npx,npy,1:2)
430  qout(npx,npy) = (extrap_corner(p0, agrid(npx-1,npy-1,1:2), agrid(npx-2,npy-2,1:2), qin(npx-1,npy-1), qin(npx-2,npy-2)) + &
431  extrap_corner(p0, agrid(npx ,npy-1,1:2), agrid(npx+1,npy-2,1:2), qin(npx ,npy-1), qin(npx+1,npy-2)) + &
432  extrap_corner(p0, agrid(npx-1,npy ,1:2), agrid(npx-2,npy+1,1:2), qin(npx-1,npy ), qin(npx-2,npy+1)))*r3
433  endif
434  if ( gridstruct%nw_corner ) then
435  p0(1:2) = grid(1,npy,1:2)
436  qout(1,npy) = (extrap_corner(p0, agrid(1,npy-1,1:2), agrid( 2,npy-2,1:2), qin(1,npy-1), qin( 2,npy-2)) + &
437  extrap_corner(p0, agrid(0,npy-1,1:2), agrid(-1,npy-2,1:2), qin(0,npy-1), qin(-1,npy-2)) + &
438  extrap_corner(p0, agrid(1,npy, 1:2), agrid( 2,npy+1,1:2), qin(1,npy ), qin( 2,npy+1)))*r3
439  endif
440 #endif
441 #endif
442 
443 !------------
444 ! X-Interior:
445 !------------
446  if (gridstruct%nested) then
447 
448  do j=js-2,je+2
449  do i=is, ie+1
450  qx(i,j) = b2*(qin(i-2,j)+qin(i+1,j)) + b1*(qin(i-1,j)+qin(i,j))
451  enddo
452  enddo
453 
454 
455  else
456 
457  do j=max(1,js-2),min(npy-1,je+2)
458  do i=max(3,is), min(npx-2,ie+1)
459  qx(i,j) = b2*(qin(i-2,j)+qin(i+1,j)) + b1*(qin(i-1,j)+qin(i,j))
460  enddo
461  enddo
462 
463 ! West Edges:
464  if ( is==1 ) then
465 
466  do j=max(1,js-2),min(npy-1,je+2)
467  gratio = dxa(2,j) / dxa(1,j)
468 #ifdef SYMM_GRID
469  qx(1,j) = 0.5*((2.+gratio)*(qin(0,j)+qin(1,j))-(qin(-1,j)+qin(2,j))) / (1.+gratio)
470 #else
471  g_in = gratio
472  g_ou = dxa(-1,j) / dxa(0,j)
473  qx(1,j) = 0.5*( ((2.+g_in)*qin(1,j)-qin( 2,j))/(1.+g_in) + &
474  ((2.+g_ou)*qin(0,j)-qin(-1,j))/(1.+g_ou) )
475 #endif
476  qx(2,j) = ( 3.*(gratio*qin(1,j)+qin(2,j))-(gratio*qx(1,j)+qx(3,j)) ) / (2.+2.*gratio)
477  enddo
478 
479  do j=max(3,js),min(npy-2,je+1)
480  qout(1,j) = a2*(qx(1,j-2)+qx(1,j+1)) + a1*(qx(1,j-1)+qx(1,j))
481  enddo
482 
483  if( js==1 ) qout(1, 2) = c1*(qx(1,1)+qx(1,2)) + c2*(qout(1,1)+qout(1,3))
484  if((je+1)==npy) qout(1,npy-1) = c1*(qx(1,npy-2)+qx(1,npy-1)) + c2*(qout(1,npy-2)+qout(1,npy))
485  endif
486 
487 ! East Edges:
488  if ( (ie+1)==npx ) then
489 
490  do j=max(1,js-2),min(npy-1,je+2)
491  gratio = dxa(npx-2,j) / dxa(npx-1,j)
492 #ifdef SYMM_GRID
493  qx(npx,j) = 0.5*((2.+gratio)*(qin(npx-1,j)+qin(npx,j)) &
494  - (qin(npx-2,j)+qin(npx+1,j))) / (1.+gratio )
495 #else
496  g_in = gratio
497  g_ou = dxa(npx+1,j) / dxa(npx,j)
498  qx(npx,j) = 0.5*( ((2.+g_in)*qin(npx-1,j)-qin(npx-2,j))/(1.+g_in) + &
499  ((2.+g_ou)*qin(npx, j)-qin(npx+1,j))/(1.+g_ou) )
500 #endif
501  qx(npx-1,j) = (3.*(qin(npx-2,j)+gratio*qin(npx-1,j)) - (gratio*qx(npx,j)+qx(npx-2,j)))/(2.+2.*gratio)
502  enddo
503 
504  do j=max(3,js),min(npy-2,je+1)
505  qout(npx,j) = a2*(qx(npx,j-2)+qx(npx,j+1)) + a1*(qx(npx,j-1)+qx(npx,j))
506  enddo
507 
508  if(js==1) qout(npx,2) = c1*(qx(npx,1)+qx(npx,2))+c2*(qout(npx,1)+qout(npx,3))
509  if((je+1)==npy) qout(npx,npy-1) = &
510  c1*(qx(npx,npy-2)+qx(npx,npy-1))+c2*(qout(npx,npy-2)+qout(npx,npy))
511  endif
512 
513  endif
514 !------------
515 ! Y-Interior:
516 !------------
517  if (gridstruct%nested) then
518 
519  do j=js,je+1
520  do i=is-2, ie+2
521  qy(i,j) = b2*(qin(i,j-2)+qin(i,j+1)) + b1*(qin(i,j-1)+qin(i,j))
522  enddo
523  enddo
524 
525  else
526 
527  do j=max(3,js),min(npy-2,je+1)
528  do i=max(1,is-2), min(npx-1,ie+2)
529  qy(i,j) = b2*(qin(i,j-2)+qin(i,j+1)) + b1*(qin(i,j-1)+qin(i,j))
530  enddo
531  enddo
532 
533 ! South Edges:
534  if ( js==1 ) then
535 
536  do i=max(1,is-2),min(npx-1,ie+2)
537  gratio = dya(i,2) / dya(i,1)
538 #ifdef SYMM_GRID
539  qy(i,1) = 0.5*((2.+gratio)*(qin(i,0)+qin(i,1)) &
540  - (qin(i,-1)+qin(i,2))) / (1.+gratio )
541 #else
542  g_in = gratio
543  g_ou = dya(i,-1) / dya(i,0)
544  qy(i,1) = 0.5*( ((2.+g_in)*qin(i,1)-qin(i,2))/(1.+g_in) + &
545  ((2.+g_ou)*qin(i,0)-qin(i,-1))/(1.+g_ou) )
546 #endif
547  qy(i,2) = (3.*(gratio*qin(i,1)+qin(i,2)) - (gratio*qy(i,1)+qy(i,3)))/(2.+2.*gratio)
548  enddo
549 
550  do i=max(3,is),min(npx-2,ie+1)
551  qout(i,1) = a2*(qy(i-2,1)+qy(i+1,1)) + a1*(qy(i-1,1)+qy(i,1))
552  enddo
553 
554  if( is==1 ) qout(2,1) = c1*(qy(1,1)+qy(2,1))+c2*(qout(1,1)+qout(3,1))
555  if((ie+1)==npx) qout(npx-1,1) = c1*(qy(npx-2,1)+qy(npx-1,1))+c2*(qout(npx-2,1)+qout(npx,1))
556  endif
557 
558 
559 ! North Edges:
560  if ( (je+1)==npy ) then
561  do i=max(1,is-2),min(npx-1,ie+2)
562  gratio = dya(i,npy-2) / dya(i,npy-1)
563 #ifdef SYMM_GRID
564  qy(i,npy ) = 0.5*((2.+gratio)*(qin(i,npy-1)+qin(i,npy)) &
565  - (qin(i,npy-2)+qin(i,npy+1))) / (1.+gratio)
566 #else
567  g_in = gratio
568  g_ou = dya(i,npy+1) / dya(i,npy)
569  qy(i,npy) = 0.5*( ((2.+g_in)*qin(i,npy-1)-qin(i,npy-2))/(1.+g_in) + &
570  ((2.+g_ou)*qin(i,npy )-qin(i,npy+1))/(1.+g_ou) )
571 #endif
572  qy(i,npy-1) = (3.*(qin(i,npy-2)+gratio*qin(i,npy-1)) - (gratio*qy(i,npy)+qy(i,npy-2)))/(2.+2.*gratio)
573  enddo
574 
575  do i=max(3,is),min(npx-2,ie+1)
576  qout(i,npy) = a2*(qy(i-2,npy)+qy(i+1,npy)) + a1*(qy(i-1,npy)+qy(i,npy))
577  enddo
578 
579  if( is==1 ) qout(2,npy) = c1*(qy(1,npy)+qy(2,npy))+c2*(qout(1,npy)+qout(3,npy))
580  if((ie+1)==npx) qout(npx-1,npy) = c1*(qy(npx-2,npy)+qy(npx-1,npy))+c2*(qout(npx-2,npy)+qout(npx,npy))
581  endif
582 
583  end if
584 
585  if (gridstruct%nested) then
586 
587  do j=js,je+1
588  do i=is,ie+1
589  qxx(i,j) = a2*(qx(i,j-2)+qx(i,j+1)) + a1*(qx(i,j-1)+qx(i,j))
590  enddo
591  enddo
592 
593  do j=js,je+1
594  do i=is,ie+1
595  qyy(i,j) = a2*(qy(i-2,j)+qy(i+1,j)) + a1*(qy(i-1,j)+qy(i,j))
596  enddo
597 
598  do i=is,ie+1
599  qout(i,j) = 0.5*(qxx(i,j) + qyy(i,j)) ! averaging
600  enddo
601  enddo
602 
603 
604  else
605 
606  do j=max(3,js),min(npy-2,je+1)
607  do i=max(2,is),min(npx-1,ie+1)
608  qxx(i,j) = a2*(qx(i,j-2)+qx(i,j+1)) + a1*(qx(i,j-1)+qx(i,j))
609  enddo
610  enddo
611 
612  if ( js==1 ) then
613  do i=max(2,is),min(npx-1,ie+1)
614  qxx(i,2) = c1*(qx(i,1)+qx(i,2))+c2*(qout(i,1)+qxx(i,3))
615  enddo
616  endif
617  if ( (je+1)==npy ) then
618  do i=max(2,is),min(npx-1,ie+1)
619  qxx(i,npy-1) = c1*(qx(i,npy-2)+qx(i,npy-1))+c2*(qout(i,npy)+qxx(i,npy-2))
620  enddo
621  endif
622 
623 
624  do j=max(2,js),min(npy-1,je+1)
625  do i=max(3,is),min(npx-2,ie+1)
626  qyy(i,j) = a2*(qy(i-2,j)+qy(i+1,j)) + a1*(qy(i-1,j)+qy(i,j))
627  enddo
628  if ( is==1 ) qyy(2,j) = c1*(qy(1,j)+qy(2,j))+c2*(qout(1,j)+qyy(3,j))
629  if((ie+1)==npx) qyy(npx-1,j) = c1*(qy(npx-2,j)+qy(npx-1,j))+c2*(qout(npx,j)+qyy(npx-2,j))
630 
631  do i=max(2,is),min(npx-1,ie+1)
632  qout(i,j) = 0.5*(qxx(i,j) + qyy(i,j)) ! averaging
633  enddo
634  enddo
635 
636  endif
637 
638  else ! grid_type>=3
639 !------------------------
640 ! Doubly periodic domain:
641 !------------------------
642 ! X-sweep: PPM
643  do j=js-2,je+2
644  do i=is,ie+1
645  qx(i,j) = b1*(qin(i-1,j)+qin(i,j)) + b2*(qin(i-2,j)+qin(i+1,j))
646  enddo
647  enddo
648 ! Y-sweep: PPM
649  do j=js,je+1
650  do i=is-2,ie+2
651  qy(i,j) = b1*(qin(i,j-1)+qin(i,j)) + b2*(qin(i,j-2)+qin(i,j+1))
652  enddo
653  enddo
654 
655  do j=js,je+1
656  do i=is,ie+1
657  qout(i,j) = 0.5*( a1*(qx(i,j-1)+qx(i,j ) + qy(i-1,j)+qy(i, j)) + &
658  a2*(qx(i,j-2)+qx(i,j+1) + qy(i-2,j)+qy(i+1,j)) )
659  enddo
660  enddo
661  endif
662 
663  if ( present(replace) ) then
664  if ( replace ) then
665  do j=js,je+1
666  do i=is,ie+1
667  qin(i,j) = qout(i,j)
668  enddo
669  enddo
670  endif
671  endif
672 
673  end subroutine a2b_ord4
674 #endif
675 
676 
677  subroutine a2b_ord2(qin, qout, gridstruct, npx, npy, is, ie, js, je, ng, replace)
678  integer, intent(IN ) :: npx, npy, is, ie, js, je, ng
679  real , intent(INOUT) :: qin(is-ng:ie+ng,js-ng:je+ng) ! A-grid field
680  real , intent( OUT) :: qout(is-ng:ie+ng,js-ng:je+ng) ! Output B-grid field
681  type(fv_grid_type), intent(IN), target :: gridstruct
682  logical, optional, intent(IN) :: replace
683  ! local:
684  real q1(npx), q2(npy)
685  integer :: i,j
686  integer :: is1, js1, is2, js2, ie1, je1
687 
688  real, pointer, dimension(:,:,:) :: grid, agrid
689  real, pointer, dimension(:,:) :: dxa, dya
690 
691  real(kind=R_GRID), pointer, dimension(:) :: edge_w, edge_e, edge_s, edge_n
692 
693  edge_w => gridstruct%edge_w
694  edge_e => gridstruct%edge_e
695  edge_s => gridstruct%edge_s
696  edge_n => gridstruct%edge_n
697 
698  grid => gridstruct%grid
699  agrid => gridstruct%agrid
700  dxa => gridstruct%dxa
701  dya => gridstruct%dya
702 
703  if (gridstruct%grid_type < 3) then
704 
705  if (gridstruct%nested) then
706 
707  do j=js-2,je+1+2
708  do i=is-2,ie+1+2
709  qout(i,j) = 0.25*(qin(i-1,j-1)+qin(i,j-1)+qin(i-1,j)+qin(i,j))
710  enddo
711  enddo
712 
713  else
714 
715  is1 = max(1,is-1)
716  js1 = max(1,js-1)
717  is2 = max(2,is)
718  js2 = max(2,js)
719 
720  ie1 = min(npx-1,ie+1)
721  je1 = min(npy-1,je+1)
722 
723  do j=js2,je1
724  do i=is2,ie1
725  qout(i,j) = 0.25*(qin(i-1,j-1)+qin(i,j-1)+qin(i-1,j)+qin(i,j))
726  enddo
727  enddo
728 
729 ! Fix the 4 Corners:
730  if ( gridstruct%sw_corner ) qout(1, 1) = r3*(qin(1, 1)+qin(1, 0)+qin(0, 1))
731  if ( gridstruct%se_corner ) qout(npx, 1) = r3*(qin(npx-1, 1)+qin(npx-1, 0)+qin(npx, 1))
732  if ( gridstruct%ne_corner ) qout(npx,npy) = r3*(qin(npx-1,npy-1)+qin(npx,npy-1)+qin(npx-1,npy))
733  if ( gridstruct%nw_corner ) qout(1, npy) = r3*(qin(1, npy-1)+qin(0, npy-1)+qin(1, npy))
734 
735  ! *** West Edges:
736  if ( is==1 ) then
737  do j=js1, je1
738  q2(j) = 0.5*(qin(0,j) + qin(1,j))
739  enddo
740  do j=js2, je1
741  qout(1,j) = edge_w(j)*q2(j-1) + (1.-edge_w(j))*q2(j)
742  enddo
743  endif
744 
745  ! East Edges:
746  if ( (ie+1)==npx ) then
747  do j=js1, je1
748  q2(j) = 0.5*(qin(npx-1,j) + qin(npx,j))
749  enddo
750  do j=js2, je1
751  qout(npx,j) = edge_e(j)*q2(j-1) + (1.-edge_e(j))*q2(j)
752  enddo
753  endif
754 
755  ! South Edges:
756  if ( js==1 ) then
757  do i=is1, ie1
758  q1(i) = 0.5*(qin(i,0) + qin(i,1))
759  enddo
760  do i=is2, ie1
761  qout(i,1) = edge_s(i)*q1(i-1) + (1.-edge_s(i))*q1(i)
762  enddo
763  endif
764 
765  ! North Edges:
766  if ( (je+1)==npy ) then
767  do i=is1, ie1
768  q1(i) = 0.5*(qin(i,npy-1) + qin(i,npy))
769  enddo
770  do i=is2, ie1
771  qout(i,npy) = edge_n(i)*q1(i-1) + (1.-edge_n(i))*q1(i)
772  enddo
773  endif
774 
775  end if
776 
777  else
778 
779  do j=js,je+1
780  do i=is,ie+1
781  qout(i,j) = 0.25*(qin(i-1,j-1)+qin(i,j-1)+qin(i-1,j)+qin(i,j))
782  enddo
783  enddo
784 
785  endif
786 
787 
788  if ( present(replace) ) then
789  if ( replace ) then
790  do j=js,je+1
791  do i=is,ie+1
792  qin(i,j) = qout(i,j)
793  enddo
794  enddo
795  endif
796  endif
797 
798  end subroutine a2b_ord2
799 
800  real function extrap_corner ( p0, p1, p2, q1, q2 )
801  real, intent(in ), dimension(2):: p0, p1, p2
802  real, intent(in ):: q1, q2
803  real:: x1, x2
804 
805  x1 = great_circle_dist( real(p1,kind=R_GRID), real(p0,kind=R_GRID) )
806  x2 = great_circle_dist( real(p2,kind=R_GRID), real(p0,kind=R_GRID) )
807 
808  extrap_corner = q1 + x1/(x2-x1) * (q1-q2)
809 
810  end function extrap_corner
811 
812 #ifdef TEST_VAND2
813  subroutine a2b_ord4(qin, qout, grid, agrid, npx, npy, is, ie, js, je, ng, replace)
814 ! use tp_core_nlm_mod, only: copy_corners
815  integer, intent(IN):: npx, npy, is, ie, js, je, ng
816  real, intent(INOUT):: qin(is-ng:ie+ng,js-ng:je+ng) ! A-grid field
817  real, intent(INOUT):: qout(is-ng:ie+ng,js-ng:je+ng) ! Output B-grid field
818  real, intent(in) :: grid(is-ng:ie+ng+1,js-ng:je+ng+1,2)
819  real, intent(in) :: agrid(is-ng:ie+ng,js-ng:je+ng,2)
820  logical, optional, intent(IN):: replace
821  real qx(is:ie+1,js-ng:je+ng)
822  real qy(is-ng:ie+ng,js:je+1)
823  real:: p0(2)
824  integer :: i, j
825 
826  real, pointer, dimension(:,:,:) :: grid, agrid
827  real, pointer, dimension(:,:) :: dxa, dya
828 
829  real, pointer, dimension(:) :: edge_w, edge_e, edge_s, edge_n
830 
831  edge_w => gridstruct%edge_w
832  edge_e => gridstruct%edge_e
833  edge_s => gridstruct%edge_s
834  edge_n => gridstruct%edge_n
835 
836  grid => gridstruct%grid
837  agrid => gridstruct%agrid
838  dxa => gridstruct%dxa
839  dya => gridstruct%dya
840 
841 
842  if (gridstruct%grid_type < 3) then
843 
844 !------------------------------------------
845 ! Copy fields to the phantom corner region:
846 !------------------------------------------
847 ! call copy_corners(qin, npx, npy, 1)
848 
849  do j=js,je+1
850  do i=is,ie+1
851 !SW:
852  if ( i==1 .and. j==1 ) goto 123
853  if ( i==2 .and. j==1 ) then
854  qin(0,-1) = qin(-1,2)
855  qin(0, 0) = qin(-1,1)
856  endif
857  if ( i==1 .and. j==2 ) then
858  qin(-1,0) = qin(2,-1)
859  qin( 0,0) = qin(1,-1)
860  endif
861  if ( i==2 .and. j==2 ) then
862  qin( 0,0) = qin(4,4)
863  endif
864 !SE:
865  if ( i==npx .and. j==1 ) goto 123
866  if ( i==npx-1 .and. j==1 ) then
867  qin(npx,-1) = qin(npx+1,2)
868  qin(npx, 0) = qin(npx+1,1)
869  endif
870  if ( i==npx-1 .and. j==2 ) then
871  qin(npx,0) = qin(npx-4,4)
872  endif
873  if ( i==npx .and. j==2 ) then
874  qin(npx+1,0) = qin(npx-2,-1)
875  qin(npx, 0) = qin(npx-1,-1)
876  endif
877 !NE:
878  if ( i==npx .and. j==npy ) goto 123
879  if ( i==npx-1 .and. j==npy-1 ) then
880  qin(npx,npy) = qin(npx-4,npy-4)
881  endif
882  if ( i==npx .and. j==npy-1 ) then
883  qin(npx+1,npy) = qin(npx-2,npy+1)
884  qin(npx, npy) = qin(npx-1,npy+1)
885  endif
886  if ( i==npx-1 .and. j==npy ) then
887  qin(npx,npy+1) = qin(npx+1,npy-2)
888  qin(npx,npy ) = qin(npx+1,npy-1)
889  endif
890 !NW:
891  if ( i==1 .and. j==npy ) goto 123
892  if ( i==1 .and. j==npy-1 ) then
893  qin(-1,npy) = qin(2,npy+1)
894  qin( 0,npy) = qin(1,npy+1)
895  endif
896  if ( i==2 .and. j==npy-1 ) then
897  qin(0,npy) = qin(4,npy-4)
898  endif
899  if ( i==2 .and. j==npy ) then
900  qin(0,npy+1) = qin(-1,npy-2)
901  qin(0,npy ) = qin(-1,npy-1)
902  endif
903 
904  qout(i,j) = van2(1, i,j)*qin(i-2,j-2) + van2(2, i,j)*qin(i-1,j-2) + &
905  van2(3, i,j)*qin(i ,j-2) + van2(4, i,j)*qin(i+1,j-2) + &
906  van2(5, i,j)*qin(i-2,j-1) + van2(6, i,j)*qin(i-1,j-1) + &
907  van2(7, i,j)*qin(i ,j-1) + van2(8, i,j)*qin(i+1,j-1) + &
908  van2(9, i,j)*qin(i-2,j ) + van2(10,i,j)*qin(i-1,j ) + &
909  van2(11,i,j)*qin(i ,j ) + van2(12,i,j)*qin(i+1,j ) + &
910  van2(13,i,j)*qin(i-2,j+1) + van2(14,i,j)*qin(i-1,j+1) + &
911  van2(15,i,j)*qin(i ,j+1) + van2(16,i,j)*qin(i+1,j+1)
912 123 continue
913  enddo
914  enddo
915 
916 ! 3-way extrapolation
917  if ( gridstruct%sw_corner ) then
918  p0(1:2) = grid(1,1,1:2)
919  qout(1,1) = (extrap_corner(p0, agrid(1,1,1:2), agrid( 2, 2,1:2), qin(1,1), qin( 2, 2)) + &
920  extrap_corner(p0, agrid(0,1,1:2), agrid(-1, 2,1:2), qin(0,1), qin(-1, 2)) + &
921  extrap_corner(p0, agrid(1,0,1:2), agrid( 2,-1,1:2), qin(1,0), qin( 2,-1)))*r3
922 
923  endif
924  if ( gridstruct%se_corner ) then
925  p0(1:2) = grid(npx,1,1:2)
926  qout(npx,1) = (extrap_corner(p0, agrid(npx-1,1,1:2), agrid(npx-2, 2,1:2), qin(npx-1,1), qin(npx-2, 2)) + &
927  extrap_corner(p0, agrid(npx-1,0,1:2), agrid(npx-2,-1,1:2), qin(npx-1,0), qin(npx-2,-1)) + &
928  extrap_corner(p0, agrid(npx ,1,1:2), agrid(npx+1, 2,1:2), qin(npx ,1), qin(npx+1, 2)))*r3
929  endif
930  if ( gridstruct%ne_corner ) then
931  p0(1:2) = grid(npx,npy,1:2)
932  qout(npx,npy) = (extrap_corner(p0, agrid(npx-1,npy-1,1:2), agrid(npx-2,npy-2,1:2), qin(npx-1,npy-1), qin(npx-2,npy-2)) + &
933  extrap_corner(p0, agrid(npx ,npy-1,1:2), agrid(npx+1,npy-2,1:2), qin(npx ,npy-1), qin(npx+1,npy-2)) + &
934  extrap_corner(p0, agrid(npx-1,npy ,1:2), agrid(npx-2,npy+1,1:2), qin(npx-1,npy ), qin(npx-2,npy+1)))*r3
935  endif
936  if ( gridstruct%nw_corner ) then
937  p0(1:2) = grid(1,npy,1:2)
938  qout(1,npy) = (extrap_corner(p0, agrid(1,npy-1,1:2), agrid( 2,npy-2,1:2), qin(1,npy-1), qin( 2,npy-2)) + &
939  extrap_corner(p0, agrid(0,npy-1,1:2), agrid(-1,npy-2,1:2), qin(0,npy-1), qin(-1,npy-2)) + &
940  extrap_corner(p0, agrid(1,npy, 1:2), agrid( 2,npy+1,1:2), qin(1,npy ), qin( 2,npy+1)))*r3
941  endif
942 
943  else ! grid_type>=3
944 
945 !------------------------
946 ! Doubly periodic domain:
947 !------------------------
948 ! X-sweep: PPM
949  do j=js-2,je+2
950  do i=is,ie+1
951  qx(i,j) = b1*(qin(i-1,j)+qin(i,j)) + b2*(qin(i-2,j)+qin(i+1,j))
952  enddo
953  enddo
954 ! Y-sweep: PPM
955  do j=js,je+1
956  do i=is-2,ie+2
957  qy(i,j) = b1*(qin(i,j-1)+qin(i,j)) + b2*(qin(i,j-2)+qin(i,j+1))
958  enddo
959  enddo
960 
961  do j=js,je+1
962  do i=is,ie+1
963  qout(i,j) = 0.5*( a1*(qx(i,j-1)+qx(i,j ) + qy(i-1,j)+qy(i, j)) + &
964  a2*(qx(i,j-2)+qx(i,j+1) + qy(i-2,j)+qy(i+1,j)) )
965  enddo
966  enddo
967 
968  endif
969 
970  if ( present(replace) ) then
971  if ( replace ) then
972  do j=js,je+1
973  do i=is,ie+1
974  qin(i,j) = qout(i,j)
975  enddo
976  enddo
977  endif
978  endif
979 
980  end subroutine a2b_ord4
981 #endif
982 
983 end module a2b_edge_nlm_mod
real, parameter a1
real, parameter r3
real, parameter b1
void a2b_ord2(int nx, int ny, const double *qin, const double *edge_w, const double *edge_e, const double *edge_s, const double *edge_n, double *qout, int on_west_edge, int on_east_edge, int on_south_edge, int on_north_edge)
Definition: gradient_c2l.c:119
subroutine, public a2b_ord4(qin, qout, gridstruct, npx, npy, is, ie, js, je, ng, replace)
real, parameter a2
integer, parameter, public r_grid
real function, public great_circle_dist(q1, q2, radius)
#define max(a, b)
Definition: mosaic_util.h:33
#define min(a, b)
Definition: mosaic_util.h:32
real function extrap_corner(p0, p1, p2, q1, q2)
real, parameter b2