FV3 Bundle
a2b_edge_tlm.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 !***********************************************************************
20 module a2b_edge_tlm_mod
21 
22  use fv_grid_utils_nlm_mod, only: great_circle_dist
23 #ifdef VAN2
24  use fv_grid_utils_nlm_mod, only: van2
25 #endif
26 
27  use fv_arrays_nlm_mod, only: fv_grid_type, r_grid
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  public :: a2b_ord2_tlm, a2b_ord4_tlm
46 
47 CONTAINS
48  SUBROUTINE a2b_ord4(qin, qout, gridstruct, npx, npy, is, ie, js, je, &
49 & ng, replace)
50  IMPLICIT NONE
51  INTEGER, INTENT(IN) :: npx, npy, is, ie, js, je, ng
52 ! A-grid field
53  REAL, INTENT(INOUT) :: qin(is-ng:ie+ng, js-ng:je+ng)
54 ! Output B-grid field
55  REAL, INTENT(INOUT) :: qout(is-ng:ie+ng, js-ng:je+ng)
56  TYPE(fv_grid_type), INTENT(IN), TARGET :: gridstruct
57  LOGICAL, OPTIONAL, INTENT(IN) :: replace
58 ! local: compact 4-pt cubic
59  REAL, PARAMETER :: c1=2./3.
60  REAL, PARAMETER :: c2=-(1./6.)
61 ! Parabolic spline
62 ! real, parameter:: c1 = 0.75
63 ! real, parameter:: c2 = -0.25
64  REAL :: qx(is:ie+1, js-ng:je+ng)
65  REAL :: qy(is-ng:ie+ng, js:je+1)
66  REAL :: qxx(is-ng:ie+ng, js-ng:je+ng)
67  REAL :: qyy(is-ng:ie+ng, js-ng:je+ng)
68  REAL :: g_in, g_ou
69  REAL :: p0(2)
70  REAL :: q1(is-1:ie+1), q2(js-1:je+1)
71  INTEGER :: i, j, is1, js1, is2, js2, ie1, je1
72  REAL, DIMENSION(:, :, :), POINTER :: grid, agrid
73  REAL, DIMENSION(:, :), POINTER :: dxa, dya
74  REAL(kind=r_grid), DIMENSION(:), POINTER :: edge_w, edge_e, edge_s, &
75 & edge_n
76  INTRINSIC max
77  INTRINSIC min
78  INTRINSIC PRESENT
79  INTEGER :: max1
80  INTEGER :: max2
81  INTEGER :: max3
82  INTEGER :: max4
83  INTEGER :: max5
84  INTEGER :: max6
85  INTEGER :: max7
86  INTEGER :: max8
87  INTEGER :: max9
88  INTEGER :: max10
89  INTEGER :: max11
90  INTEGER :: max12
91  INTEGER :: max13
92  INTEGER :: max14
93  INTEGER :: max15
94  INTEGER :: min1
95  INTEGER :: min2
96  INTEGER :: min3
97  INTEGER :: min4
98  INTEGER :: min5
99  INTEGER :: min6
100  INTEGER :: min7
101  INTEGER :: min8
102  INTEGER :: min9
103  INTEGER :: min10
104  INTEGER :: min11
105  INTEGER :: min12
106  INTEGER :: min13
107  INTEGER :: min14
108  INTEGER :: min15
109  REAL :: result1
110  REAL :: result2
111  REAL :: result3
112  edge_w => gridstruct%edge_w
113  edge_e => gridstruct%edge_e
114  edge_s => gridstruct%edge_s
115  edge_n => gridstruct%edge_n
116  grid => gridstruct%grid
117  agrid => gridstruct%agrid
118  dxa => gridstruct%dxa
119  dya => gridstruct%dya
120  IF (gridstruct%grid_type .LT. 3) THEN
121  IF (1 .LT. is - 1) THEN
122  is1 = is - 1
123  ELSE
124  is1 = 1
125  END IF
126  IF (1 .LT. js - 1) THEN
127  js1 = js - 1
128  ELSE
129  js1 = 1
130  END IF
131  IF (2 .LT. is) THEN
132  is2 = is
133  ELSE
134  is2 = 2
135  END IF
136  IF (2 .LT. js) THEN
137  js2 = js
138  ELSE
139  js2 = 2
140  END IF
141  IF (npx - 1 .GT. ie + 1) THEN
142  ie1 = ie + 1
143  ELSE
144  ie1 = npx - 1
145  END IF
146  IF (npy - 1 .GT. je + 1) THEN
147  je1 = je + 1
148  ELSE
149  je1 = npy - 1
150  END IF
151 ! Corners:
152 ! 3-way extrapolation
153  IF (gridstruct%nested) THEN
154  DO j=js-2,je+2
155  DO i=is,ie+1
156  qx(i, j) = b2*(qin(i-2, j)+qin(i+1, j)) + b1*(qin(i-1, j)+&
157 & qin(i, j))
158  END DO
159  END DO
160  ELSE
161  IF (gridstruct%sw_corner) THEN
162  p0(1:2) = grid(1, 1, 1:2)
163  result1 = extrap_corner(p0, agrid(1, 1, 1:2), agrid(2, 2, 1:2)&
164 & , qin(1, 1), qin(2, 2))
165  result2 = extrap_corner(p0, agrid(0, 1, 1:2), agrid(-1, 2, 1:2&
166 & ), qin(0, 1), qin(-1, 2))
167  result3 = extrap_corner(p0, agrid(1, 0, 1:2), agrid(2, -1, 1:2&
168 & ), qin(1, 0), qin(2, -1))
169  qout(1, 1) = (result1+result2+result3)*r3
170  END IF
171  IF (gridstruct%se_corner) THEN
172  p0(1:2) = grid(npx, 1, 1:2)
173  result1 = extrap_corner(p0, agrid(npx-1, 1, 1:2), agrid(npx-2&
174 & , 2, 1:2), qin(npx-1, 1), qin(npx-2, 2))
175  result2 = extrap_corner(p0, agrid(npx-1, 0, 1:2), agrid(npx-2&
176 & , -1, 1:2), qin(npx-1, 0), qin(npx-2, -1))
177  result3 = extrap_corner(p0, agrid(npx, 1, 1:2), agrid(npx+1, 2&
178 & , 1:2), qin(npx, 1), qin(npx+1, 2))
179  qout(npx, 1) = (result1+result2+result3)*r3
180  END IF
181  IF (gridstruct%ne_corner) THEN
182  p0(1:2) = grid(npx, npy, 1:2)
183  result1 = extrap_corner(p0, agrid(npx-1, npy-1, 1:2), agrid(&
184 & npx-2, npy-2, 1:2), qin(npx-1, npy-1), qin(npx-2, npy-2))
185  result2 = extrap_corner(p0, agrid(npx, npy-1, 1:2), agrid(npx+&
186 & 1, npy-2, 1:2), qin(npx, npy-1), qin(npx+1, npy-2))
187  result3 = extrap_corner(p0, agrid(npx-1, npy, 1:2), agrid(npx-&
188 & 2, npy+1, 1:2), qin(npx-1, npy), qin(npx-2, npy+1))
189  qout(npx, npy) = (result1+result2+result3)*r3
190  END IF
191  IF (gridstruct%nw_corner) THEN
192  p0(1:2) = grid(1, npy, 1:2)
193  result1 = extrap_corner(p0, agrid(1, npy-1, 1:2), agrid(2, npy&
194 & -2, 1:2), qin(1, npy-1), qin(2, npy-2))
195  result2 = extrap_corner(p0, agrid(0, npy-1, 1:2), agrid(-1, &
196 & npy-2, 1:2), qin(0, npy-1), qin(-1, npy-2))
197  result3 = extrap_corner(p0, agrid(1, npy, 1:2), agrid(2, npy+1&
198 & , 1:2), qin(1, npy), qin(2, npy+1))
199  qout(1, npy) = (result1+result2+result3)*r3
200  END IF
201  IF (1 .LT. js - 2) THEN
202  max1 = js - 2
203  ELSE
204  max1 = 1
205  END IF
206  IF (npy - 1 .GT. je + 2) THEN
207  min1 = je + 2
208  ELSE
209  min1 = npy - 1
210  END IF
211 !------------
212 ! X-Interior:
213 !------------
214  DO j=max1,min1
215  IF (3 .LT. is) THEN
216  max2 = is
217  ELSE
218  max2 = 3
219  END IF
220  IF (npx - 2 .GT. ie + 1) THEN
221  min2 = ie + 1
222  ELSE
223  min2 = npx - 2
224  END IF
225  DO i=max2,min2
226  qx(i, j) = b2*(qin(i-2, j)+qin(i+1, j)) + b1*(qin(i-1, j)+&
227 & qin(i, j))
228  END DO
229  END DO
230 ! *** West Edges:
231  IF (is .EQ. 1) THEN
232  DO j=js1,je1
233  q2(j) = (qin(0, j)*dxa(1, j)+qin(1, j)*dxa(0, j))/(dxa(0, j)&
234 & +dxa(1, j))
235  END DO
236  DO j=js2,je1
237  qout(1, j) = edge_w(j)*q2(j-1) + (1.-edge_w(j))*q2(j)
238  END DO
239  IF (1 .LT. js - 2) THEN
240  max3 = js - 2
241  ELSE
242  max3 = 1
243  END IF
244  IF (npy - 1 .GT. je + 2) THEN
245  min3 = je + 2
246  ELSE
247  min3 = npy - 1
248  END IF
249 !
250  DO j=max3,min3
251  g_in = dxa(2, j)/dxa(1, j)
252  g_ou = dxa(-1, j)/dxa(0, j)
253  qx(1, j) = 0.5*(((2.+g_in)*qin(1, j)-qin(2, j))/(1.+g_in)+((&
254 & 2.+g_ou)*qin(0, j)-qin(-1, j))/(1.+g_ou))
255  qx(2, j) = (3.*(g_in*qin(1, j)+qin(2, j))-(g_in*qx(1, j)+qx(&
256 & 3, j)))/(2.+2.*g_in)
257  END DO
258  END IF
259 ! East Edges:
260  IF (ie + 1 .EQ. npx) THEN
261  DO j=js1,je1
262  q2(j) = (qin(npx-1, j)*dxa(npx, j)+qin(npx, j)*dxa(npx-1, j)&
263 & )/(dxa(npx-1, j)+dxa(npx, j))
264  END DO
265  DO j=js2,je1
266  qout(npx, j) = edge_e(j)*q2(j-1) + (1.-edge_e(j))*q2(j)
267  END DO
268  IF (1 .LT. js - 2) THEN
269  max4 = js - 2
270  ELSE
271  max4 = 1
272  END IF
273  IF (npy - 1 .GT. je + 2) THEN
274  min4 = je + 2
275  ELSE
276  min4 = npy - 1
277  END IF
278 !
279  DO j=max4,min4
280  g_in = dxa(npx-2, j)/dxa(npx-1, j)
281  g_ou = dxa(npx+1, j)/dxa(npx, j)
282  qx(npx, j) = 0.5*(((2.+g_in)*qin(npx-1, j)-qin(npx-2, j))/(&
283 & 1.+g_in)+((2.+g_ou)*qin(npx, j)-qin(npx+1, j))/(1.+g_ou))
284  qx(npx-1, j) = (3.*(qin(npx-2, j)+g_in*qin(npx-1, j))-(g_in*&
285 & qx(npx, j)+qx(npx-2, j)))/(2.+2.*g_in)
286  END DO
287  END IF
288  END IF
289 !------------
290 ! Y-Interior:
291 !------------
292  IF (gridstruct%nested) THEN
293  DO j=js,je+1
294  DO i=is-2,ie+2
295  qy(i, j) = b2*(qin(i, j-2)+qin(i, j+1)) + b1*(qin(i, j-1)+&
296 & qin(i, j))
297  END DO
298  END DO
299  ELSE
300  IF (3 .LT. js) THEN
301  max5 = js
302  ELSE
303  max5 = 3
304  END IF
305  IF (npy - 2 .GT. je + 1) THEN
306  min5 = je + 1
307  ELSE
308  min5 = npy - 2
309  END IF
310  DO j=max5,min5
311  IF (1 .LT. is - 2) THEN
312  max6 = is - 2
313  ELSE
314  max6 = 1
315  END IF
316  IF (npx - 1 .GT. ie + 2) THEN
317  min6 = ie + 2
318  ELSE
319  min6 = npx - 1
320  END IF
321  DO i=max6,min6
322  qy(i, j) = b2*(qin(i, j-2)+qin(i, j+1)) + b1*(qin(i, j-1)+&
323 & qin(i, j))
324  END DO
325  END DO
326 ! South Edges:
327  IF (js .EQ. 1) THEN
328  DO i=is1,ie1
329  q1(i) = (qin(i, 0)*dya(i, 1)+qin(i, 1)*dya(i, 0))/(dya(i, 0)&
330 & +dya(i, 1))
331  END DO
332  DO i=is2,ie1
333  qout(i, 1) = edge_s(i)*q1(i-1) + (1.-edge_s(i))*q1(i)
334  END DO
335  IF (1 .LT. is - 2) THEN
336  max7 = is - 2
337  ELSE
338  max7 = 1
339  END IF
340  IF (npx - 1 .GT. ie + 2) THEN
341  min7 = ie + 2
342  ELSE
343  min7 = npx - 1
344  END IF
345 !
346  DO i=max7,min7
347  g_in = dya(i, 2)/dya(i, 1)
348  g_ou = dya(i, -1)/dya(i, 0)
349  qy(i, 1) = 0.5*(((2.+g_in)*qin(i, 1)-qin(i, 2))/(1.+g_in)+((&
350 & 2.+g_ou)*qin(i, 0)-qin(i, -1))/(1.+g_ou))
351  qy(i, 2) = (3.*(g_in*qin(i, 1)+qin(i, 2))-(g_in*qy(i, 1)+qy(&
352 & i, 3)))/(2.+2.*g_in)
353  END DO
354  END IF
355 ! North Edges:
356  IF (je + 1 .EQ. npy) THEN
357  DO i=is1,ie1
358  q1(i) = (qin(i, npy-1)*dya(i, npy)+qin(i, npy)*dya(i, npy-1)&
359 & )/(dya(i, npy-1)+dya(i, npy))
360  END DO
361  DO i=is2,ie1
362  qout(i, npy) = edge_n(i)*q1(i-1) + (1.-edge_n(i))*q1(i)
363  END DO
364  IF (1 .LT. is - 2) THEN
365  max8 = is - 2
366  ELSE
367  max8 = 1
368  END IF
369  IF (npx - 1 .GT. ie + 2) THEN
370  min8 = ie + 2
371  ELSE
372  min8 = npx - 1
373  END IF
374 !
375  DO i=max8,min8
376  g_in = dya(i, npy-2)/dya(i, npy-1)
377  g_ou = dya(i, npy+1)/dya(i, npy)
378  qy(i, npy) = 0.5*(((2.+g_in)*qin(i, npy-1)-qin(i, npy-2))/(&
379 & 1.+g_in)+((2.+g_ou)*qin(i, npy)-qin(i, npy+1))/(1.+g_ou))
380  qy(i, npy-1) = (3.*(qin(i, npy-2)+g_in*qin(i, npy-1))-(g_in*&
381 & qy(i, npy)+qy(i, npy-2)))/(2.+2.*g_in)
382  END DO
383  END IF
384  END IF
385 !--------------------------------------
386  IF (gridstruct%nested) THEN
387  DO j=js,je+1
388  DO i=is,ie+1
389  qxx(i, j) = a2*(qx(i, j-2)+qx(i, j+1)) + a1*(qx(i, j-1)+qx(i&
390 & , j))
391  END DO
392  END DO
393  DO j=js,je+1
394  DO i=is,ie+1
395  qyy(i, j) = a2*(qy(i-2, j)+qy(i+1, j)) + a1*(qy(i-1, j)+qy(i&
396 & , j))
397  END DO
398  DO i=is,ie+1
399 ! averaging
400  qout(i, j) = 0.5*(qxx(i, j)+qyy(i, j))
401  END DO
402  END DO
403  ELSE
404  IF (3 .LT. js) THEN
405  max9 = js
406  ELSE
407  max9 = 3
408  END IF
409  IF (npy - 2 .GT. je + 1) THEN
410  min9 = je + 1
411  ELSE
412  min9 = npy - 2
413  END IF
414  DO j=max9,min9
415  IF (2 .LT. is) THEN
416  max10 = is
417  ELSE
418  max10 = 2
419  END IF
420  IF (npx - 1 .GT. ie + 1) THEN
421  min10 = ie + 1
422  ELSE
423  min10 = npx - 1
424  END IF
425  DO i=max10,min10
426  qxx(i, j) = a2*(qx(i, j-2)+qx(i, j+1)) + a1*(qx(i, j-1)+qx(i&
427 & , j))
428  END DO
429  END DO
430  IF (js .EQ. 1) THEN
431  IF (2 .LT. is) THEN
432  max11 = is
433  ELSE
434  max11 = 2
435  END IF
436  IF (npx - 1 .GT. ie + 1) THEN
437  min11 = ie + 1
438  ELSE
439  min11 = npx - 1
440  END IF
441  DO i=max11,min11
442  qxx(i, 2) = c1*(qx(i, 1)+qx(i, 2)) + c2*(qout(i, 1)+qxx(i, 3&
443 & ))
444  END DO
445  END IF
446  IF (je + 1 .EQ. npy) THEN
447  IF (2 .LT. is) THEN
448  max12 = is
449  ELSE
450  max12 = 2
451  END IF
452  IF (npx - 1 .GT. ie + 1) THEN
453  min12 = ie + 1
454  ELSE
455  min12 = npx - 1
456  END IF
457  DO i=max12,min12
458  qxx(i, npy-1) = c1*(qx(i, npy-2)+qx(i, npy-1)) + c2*(qout(i&
459 & , npy)+qxx(i, npy-2))
460  END DO
461  END IF
462  IF (2 .LT. js) THEN
463  max13 = js
464  ELSE
465  max13 = 2
466  END IF
467  IF (npy - 1 .GT. je + 1) THEN
468  min13 = je + 1
469  ELSE
470  min13 = npy - 1
471  END IF
472  DO j=max13,min13
473  IF (3 .LT. is) THEN
474  max14 = is
475  ELSE
476  max14 = 3
477  END IF
478  IF (npx - 2 .GT. ie + 1) THEN
479  min14 = ie + 1
480  ELSE
481  min14 = npx - 2
482  END IF
483  DO i=max14,min14
484  qyy(i, j) = a2*(qy(i-2, j)+qy(i+1, j)) + a1*(qy(i-1, j)+qy(i&
485 & , j))
486  END DO
487  IF (is .EQ. 1) qyy(2, j) = c1*(qy(1, j)+qy(2, j)) + c2*(qout(1&
488 & , j)+qyy(3, j))
489  IF (ie + 1 .EQ. npx) qyy(npx-1, j) = c1*(qy(npx-2, j)+qy(npx-1&
490 & , j)) + c2*(qout(npx, j)+qyy(npx-2, j))
491  IF (2 .LT. is) THEN
492  max15 = is
493  ELSE
494  max15 = 2
495  END IF
496  IF (npx - 1 .GT. ie + 1) THEN
497  min15 = ie + 1
498  ELSE
499  min15 = npx - 1
500  END IF
501  DO i=max15,min15
502 ! averaging
503  qout(i, j) = 0.5*(qxx(i, j)+qyy(i, j))
504  END DO
505  END DO
506  END IF
507  ELSE
508 ! grid_type>=3
509 !------------------------
510 ! Doubly periodic domain:
511 !------------------------
512 ! X-sweep: PPM
513  DO j=js-2,je+2
514  DO i=is,ie+1
515  qx(i, j) = b1*(qin(i-1, j)+qin(i, j)) + b2*(qin(i-2, j)+qin(i+&
516 & 1, j))
517  END DO
518  END DO
519 ! Y-sweep: PPM
520  DO j=js,je+1
521  DO i=is-2,ie+2
522  qy(i, j) = b1*(qin(i, j-1)+qin(i, j)) + b2*(qin(i, j-2)+qin(i&
523 & , j+1))
524  END DO
525  END DO
526  DO j=js,je+1
527  DO i=is,ie+1
528  qout(i, j) = 0.5*(a1*(qx(i, j-1)+qx(i, j)+qy(i-1, j)+qy(i, j))&
529 & +a2*(qx(i, j-2)+qx(i, j+1)+qy(i-2, j)+qy(i+1, j)))
530  END DO
531  END DO
532  END IF
533  IF (PRESENT(replace)) THEN
534  IF (replace) THEN
535  DO j=js,je+1
536  DO i=is,ie+1
537  qin(i, j) = qout(i, j)
538  END DO
539  END DO
540  END IF
541  END IF
542  END SUBROUTINE a2b_ord4
543 ! Differentiation of a2b_ord4 in forward (tangent) mode:
544 ! variations of useful results: qin qout
545 ! with respect to varying inputs: qin qout
546  SUBROUTINE a2b_ord4_tlm(qin, qin_tl, qout, qout_tl, gridstruct, npx&
547 & , npy, is, ie, js, je, ng, replace)
548  IMPLICIT NONE
549  INTEGER, INTENT(IN) :: npx, npy, is, ie, js, je, ng
550 ! A-grid field
551  REAL, INTENT(INOUT) :: qin(is-ng:ie+ng, js-ng:je+ng)
552  REAL, INTENT(INOUT) :: qin_tl(is-ng:ie+ng, js-ng:je+ng)
553 ! Output B-grid field
554  REAL, INTENT(INOUT) :: qout(is-ng:ie+ng, js-ng:je+ng)
555  REAL, INTENT(INOUT) :: qout_tl(is-ng:ie+ng, js-ng:je+ng)
556  TYPE(fv_grid_type), INTENT(IN), TARGET :: gridstruct
557  LOGICAL, OPTIONAL, INTENT(IN) :: replace
558 ! local: compact 4-pt cubic
559  REAL, PARAMETER :: c1=2./3.
560  REAL, PARAMETER :: c2=-(1./6.)
561 ! Parabolic spline
562 ! real, parameter:: c1 = 0.75
563 ! real, parameter:: c2 = -0.25
564  REAL :: qx(is:ie+1, js-ng:je+ng)
565  REAL :: qx_tl(is:ie+1, js-ng:je+ng)
566  REAL :: qy(is-ng:ie+ng, js:je+1)
567  REAL :: qy_tl(is-ng:ie+ng, js:je+1)
568  REAL :: qxx(is-ng:ie+ng, js-ng:je+ng)
569  REAL :: qxx_tl(is-ng:ie+ng, js-ng:je+ng)
570  REAL :: qyy(is-ng:ie+ng, js-ng:je+ng)
571  REAL :: qyy_tl(is-ng:ie+ng, js-ng:je+ng)
572  REAL :: g_in, g_ou
573  REAL :: p0(2)
574  REAL :: q1(is-1:ie+1), q2(js-1:je+1)
575  REAL :: q1_tl(is-1:ie+1), q2_tl(js-1:je+1)
576  INTEGER :: i, j, is1, js1, is2, js2, ie1, je1
577  REAL, DIMENSION(:, :, :), POINTER :: grid, agrid
578  REAL, DIMENSION(:, :), POINTER :: dxa, dya
579  REAL(kind=r_grid), DIMENSION(:), POINTER :: edge_w, edge_e, edge_s, &
580 & edge_n
581  INTRINSIC max
582  INTRINSIC min
583  INTRINSIC PRESENT
584  INTEGER :: max1
585  INTEGER :: max2
586  INTEGER :: max3
587  INTEGER :: max4
588  INTEGER :: max5
589  INTEGER :: max6
590  INTEGER :: max7
591  INTEGER :: max8
592  INTEGER :: max9
593  INTEGER :: max10
594  INTEGER :: max11
595  INTEGER :: max12
596  INTEGER :: max13
597  INTEGER :: max14
598  INTEGER :: max15
599  INTEGER :: min1
600  INTEGER :: min2
601  INTEGER :: min3
602  INTEGER :: min4
603  INTEGER :: min5
604  INTEGER :: min6
605  INTEGER :: min7
606  INTEGER :: min8
607  INTEGER :: min9
608  INTEGER :: min10
609  INTEGER :: min11
610  INTEGER :: min12
611  INTEGER :: min13
612  INTEGER :: min14
613  INTEGER :: min15
614  REAL :: result1
615  REAL :: result1_tl
616  REAL :: result2
617  REAL :: result2_tl
618  REAL :: result3
619  REAL :: result3_tl
620  edge_w => gridstruct%edge_w
621  edge_e => gridstruct%edge_e
622  edge_s => gridstruct%edge_s
623  edge_n => gridstruct%edge_n
624  grid => gridstruct%grid
625  agrid => gridstruct%agrid
626  dxa => gridstruct%dxa
627  dya => gridstruct%dya
628  IF (gridstruct%grid_type .LT. 3) THEN
629  IF (1 .LT. is - 1) THEN
630  is1 = is - 1
631  ELSE
632  is1 = 1
633  END IF
634  IF (1 .LT. js - 1) THEN
635  js1 = js - 1
636  ELSE
637  js1 = 1
638  END IF
639  IF (2 .LT. is) THEN
640  is2 = is
641  ELSE
642  is2 = 2
643  END IF
644  IF (2 .LT. js) THEN
645  js2 = js
646  ELSE
647  js2 = 2
648  END IF
649  IF (npx - 1 .GT. ie + 1) THEN
650  ie1 = ie + 1
651  ELSE
652  ie1 = npx - 1
653  END IF
654  IF (npy - 1 .GT. je + 1) THEN
655  je1 = je + 1
656  ELSE
657  je1 = npy - 1
658  END IF
659 ! Corners:
660 ! 3-way extrapolation
661  IF (gridstruct%nested) THEN
662  qx_tl = 0.0
663  DO j=js-2,je+2
664  DO i=is,ie+1
665  qx_tl(i, j) = b2*(qin_tl(i-2, j)+qin_tl(i+1, j)) + b1*(&
666 & qin_tl(i-1, j)+qin_tl(i, j))
667  qx(i, j) = b2*(qin(i-2, j)+qin(i+1, j)) + b1*(qin(i-1, j)+&
668 & qin(i, j))
669  END DO
670  END DO
671  ELSE
672  IF (gridstruct%sw_corner) THEN
673  p0(1:2) = grid(1, 1, 1:2)
674  result1_tl = extrap_corner_tlm(p0, agrid(1, 1, 1:2), agrid(&
675 & 2, 2, 1:2), qin(1, 1), qin_tl(1, 1), qin(2, 2), qin_tl(2, 2)&
676 & , result1)
677  result2_tl = extrap_corner_tlm(p0, agrid(0, 1, 1:2), agrid(&
678 & -1, 2, 1:2), qin(0, 1), qin_tl(0, 1), qin(-1, 2), qin_tl(-1&
679 & , 2), result2)
680  result3_tl = extrap_corner_tlm(p0, agrid(1, 0, 1:2), agrid(&
681 & 2, -1, 1:2), qin(1, 0), qin_tl(1, 0), qin(2, -1), qin_tl(2, &
682 & -1), result3)
683  qout_tl(1, 1) = r3*(result1_tl+result2_tl+result3_tl)
684  qout(1, 1) = (result1+result2+result3)*r3
685  END IF
686  IF (gridstruct%se_corner) THEN
687  p0(1:2) = grid(npx, 1, 1:2)
688  result1_tl = extrap_corner_tlm(p0, agrid(npx-1, 1, 1:2), &
689 & agrid(npx-2, 2, 1:2), qin(npx-1, 1), qin_tl(npx-1, 1), qin(&
690 & npx-2, 2), qin_tl(npx-2, 2), result1)
691  result2_tl = extrap_corner_tlm(p0, agrid(npx-1, 0, 1:2), &
692 & agrid(npx-2, -1, 1:2), qin(npx-1, 0), qin_tl(npx-1, 0), qin(&
693 & npx-2, -1), qin_tl(npx-2, -1), result2)
694  result3_tl = extrap_corner_tlm(p0, agrid(npx, 1, 1:2), &
695 & agrid(npx+1, 2, 1:2), qin(npx, 1), qin_tl(npx, 1), qin(npx+1&
696 & , 2), qin_tl(npx+1, 2), result3)
697  qout_tl(npx, 1) = r3*(result1_tl+result2_tl+result3_tl)
698  qout(npx, 1) = (result1+result2+result3)*r3
699  END IF
700  IF (gridstruct%ne_corner) THEN
701  p0(1:2) = grid(npx, npy, 1:2)
702  result1_tl = extrap_corner_tlm(p0, agrid(npx-1, npy-1, 1:2)&
703 & , agrid(npx-2, npy-2, 1:2), qin(npx-1, npy-1), qin_tl(npx-1&
704 & , npy-1), qin(npx-2, npy-2), qin_tl(npx-2, npy-2), result1)
705  result2_tl = extrap_corner_tlm(p0, agrid(npx, npy-1, 1:2), &
706 & agrid(npx+1, npy-2, 1:2), qin(npx, npy-1), qin_tl(npx, npy-1&
707 & ), qin(npx+1, npy-2), qin_tl(npx+1, npy-2), result2)
708  result3_tl = extrap_corner_tlm(p0, agrid(npx-1, npy, 1:2), &
709 & agrid(npx-2, npy+1, 1:2), qin(npx-1, npy), qin_tl(npx-1, npy&
710 & ), qin(npx-2, npy+1), qin_tl(npx-2, npy+1), result3)
711  qout_tl(npx, npy) = r3*(result1_tl+result2_tl+result3_tl)
712  qout(npx, npy) = (result1+result2+result3)*r3
713  END IF
714  IF (gridstruct%nw_corner) THEN
715  p0(1:2) = grid(1, npy, 1:2)
716  result1_tl = extrap_corner_tlm(p0, agrid(1, npy-1, 1:2), &
717 & agrid(2, npy-2, 1:2), qin(1, npy-1), qin_tl(1, npy-1), qin(2&
718 & , npy-2), qin_tl(2, npy-2), result1)
719  result2_tl = extrap_corner_tlm(p0, agrid(0, npy-1, 1:2), &
720 & agrid(-1, npy-2, 1:2), qin(0, npy-1), qin_tl(0, npy-1), qin(&
721 & -1, npy-2), qin_tl(-1, npy-2), result2)
722  result3_tl = extrap_corner_tlm(p0, agrid(1, npy, 1:2), &
723 & agrid(2, npy+1, 1:2), qin(1, npy), qin_tl(1, npy), qin(2, &
724 & npy+1), qin_tl(2, npy+1), result3)
725  qout_tl(1, npy) = r3*(result1_tl+result2_tl+result3_tl)
726  qout(1, npy) = (result1+result2+result3)*r3
727  END IF
728  IF (1 .LT. js - 2) THEN
729  max1 = js - 2
730  ELSE
731  max1 = 1
732  END IF
733  IF (npy - 1 .GT. je + 2) THEN
734  min1 = je + 2
735  qx_tl = 0.0
736  ELSE
737  min1 = npy - 1
738  qx_tl = 0.0
739  END IF
740 !------------
741 ! X-Interior:
742 !------------
743  DO j=max1,min1
744  IF (3 .LT. is) THEN
745  max2 = is
746  ELSE
747  max2 = 3
748  END IF
749  IF (npx - 2 .GT. ie + 1) THEN
750  min2 = ie + 1
751  ELSE
752  min2 = npx - 2
753  END IF
754  DO i=max2,min2
755  qx_tl(i, j) = b2*(qin_tl(i-2, j)+qin_tl(i+1, j)) + b1*(&
756 & qin_tl(i-1, j)+qin_tl(i, j))
757  qx(i, j) = b2*(qin(i-2, j)+qin(i+1, j)) + b1*(qin(i-1, j)+&
758 & qin(i, j))
759  END DO
760  END DO
761 ! *** West Edges:
762  IF (is .EQ. 1) THEN
763  q2_tl = 0.0
764  DO j=js1,je1
765  q2_tl(j) = (dxa(1, j)*qin_tl(0, j)+dxa(0, j)*qin_tl(1, j))/(&
766 & dxa(0, j)+dxa(1, j))
767  q2(j) = (qin(0, j)*dxa(1, j)+qin(1, j)*dxa(0, j))/(dxa(0, j)&
768 & +dxa(1, j))
769  END DO
770  DO j=js2,je1
771  qout_tl(1, j) = edge_w(j)*q2_tl(j-1) + (1.-edge_w(j))*q2_tl(&
772 & j)
773  qout(1, j) = edge_w(j)*q2(j-1) + (1.-edge_w(j))*q2(j)
774  END DO
775  IF (1 .LT. js - 2) THEN
776  max3 = js - 2
777  ELSE
778  max3 = 1
779  END IF
780  IF (npy - 1 .GT. je + 2) THEN
781  min3 = je + 2
782  ELSE
783  min3 = npy - 1
784  END IF
785 !
786  DO j=max3,min3
787  g_in = dxa(2, j)/dxa(1, j)
788  g_ou = dxa(-1, j)/dxa(0, j)
789  qx_tl(1, j) = 0.5*(((2.+g_in)*qin_tl(1, j)-qin_tl(2, j))/(1.&
790 & +g_in)+((2.+g_ou)*qin_tl(0, j)-qin_tl(-1, j))/(1.+g_ou))
791  qx(1, j) = 0.5*(((2.+g_in)*qin(1, j)-qin(2, j))/(1.+g_in)+((&
792 & 2.+g_ou)*qin(0, j)-qin(-1, j))/(1.+g_ou))
793  qx_tl(2, j) = (3.*(g_in*qin_tl(1, j)+qin_tl(2, j))-g_in*&
794 & qx_tl(1, j)-qx_tl(3, j))/(2.+2.*g_in)
795  qx(2, j) = (3.*(g_in*qin(1, j)+qin(2, j))-(g_in*qx(1, j)+qx(&
796 & 3, j)))/(2.+2.*g_in)
797  END DO
798  ELSE
799  q2_tl = 0.0
800  END IF
801 ! East Edges:
802  IF (ie + 1 .EQ. npx) THEN
803  DO j=js1,je1
804  q2_tl(j) = (dxa(npx, j)*qin_tl(npx-1, j)+dxa(npx-1, j)*&
805 & qin_tl(npx, j))/(dxa(npx-1, j)+dxa(npx, j))
806  q2(j) = (qin(npx-1, j)*dxa(npx, j)+qin(npx, j)*dxa(npx-1, j)&
807 & )/(dxa(npx-1, j)+dxa(npx, j))
808  END DO
809  DO j=js2,je1
810  qout_tl(npx, j) = edge_e(j)*q2_tl(j-1) + (1.-edge_e(j))*&
811 & q2_tl(j)
812  qout(npx, j) = edge_e(j)*q2(j-1) + (1.-edge_e(j))*q2(j)
813  END DO
814  IF (1 .LT. js - 2) THEN
815  max4 = js - 2
816  ELSE
817  max4 = 1
818  END IF
819  IF (npy - 1 .GT. je + 2) THEN
820  min4 = je + 2
821  ELSE
822  min4 = npy - 1
823  END IF
824 !
825  DO j=max4,min4
826  g_in = dxa(npx-2, j)/dxa(npx-1, j)
827  g_ou = dxa(npx+1, j)/dxa(npx, j)
828  qx_tl(npx, j) = 0.5*(((2.+g_in)*qin_tl(npx-1, j)-qin_tl(npx-&
829 & 2, j))/(1.+g_in)+((2.+g_ou)*qin_tl(npx, j)-qin_tl(npx+1, j&
830 & ))/(1.+g_ou))
831  qx(npx, j) = 0.5*(((2.+g_in)*qin(npx-1, j)-qin(npx-2, j))/(&
832 & 1.+g_in)+((2.+g_ou)*qin(npx, j)-qin(npx+1, j))/(1.+g_ou))
833  qx_tl(npx-1, j) = (3.*(qin_tl(npx-2, j)+g_in*qin_tl(npx-1, j&
834 & ))-g_in*qx_tl(npx, j)-qx_tl(npx-2, j))/(2.+2.*g_in)
835  qx(npx-1, j) = (3.*(qin(npx-2, j)+g_in*qin(npx-1, j))-(g_in*&
836 & qx(npx, j)+qx(npx-2, j)))/(2.+2.*g_in)
837  END DO
838  END IF
839  END IF
840 !------------
841 ! Y-Interior:
842 !------------
843  IF (gridstruct%nested) THEN
844  qy_tl = 0.0
845  DO j=js,je+1
846  DO i=is-2,ie+2
847  qy_tl(i, j) = b2*(qin_tl(i, j-2)+qin_tl(i, j+1)) + b1*(&
848 & qin_tl(i, j-1)+qin_tl(i, j))
849  qy(i, j) = b2*(qin(i, j-2)+qin(i, j+1)) + b1*(qin(i, j-1)+&
850 & qin(i, j))
851  END DO
852  END DO
853  ELSE
854  IF (3 .LT. js) THEN
855  max5 = js
856  ELSE
857  max5 = 3
858  END IF
859  IF (npy - 2 .GT. je + 1) THEN
860  min5 = je + 1
861  qy_tl = 0.0
862  ELSE
863  min5 = npy - 2
864  qy_tl = 0.0
865  END IF
866  DO j=max5,min5
867  IF (1 .LT. is - 2) THEN
868  max6 = is - 2
869  ELSE
870  max6 = 1
871  END IF
872  IF (npx - 1 .GT. ie + 2) THEN
873  min6 = ie + 2
874  ELSE
875  min6 = npx - 1
876  END IF
877  DO i=max6,min6
878  qy_tl(i, j) = b2*(qin_tl(i, j-2)+qin_tl(i, j+1)) + b1*(&
879 & qin_tl(i, j-1)+qin_tl(i, j))
880  qy(i, j) = b2*(qin(i, j-2)+qin(i, j+1)) + b1*(qin(i, j-1)+&
881 & qin(i, j))
882  END DO
883  END DO
884 ! South Edges:
885  IF (js .EQ. 1) THEN
886  q1_tl = 0.0
887  DO i=is1,ie1
888  q1_tl(i) = (dya(i, 1)*qin_tl(i, 0)+dya(i, 0)*qin_tl(i, 1))/(&
889 & dya(i, 0)+dya(i, 1))
890  q1(i) = (qin(i, 0)*dya(i, 1)+qin(i, 1)*dya(i, 0))/(dya(i, 0)&
891 & +dya(i, 1))
892  END DO
893  DO i=is2,ie1
894  qout_tl(i, 1) = edge_s(i)*q1_tl(i-1) + (1.-edge_s(i))*q1_tl(&
895 & i)
896  qout(i, 1) = edge_s(i)*q1(i-1) + (1.-edge_s(i))*q1(i)
897  END DO
898  IF (1 .LT. is - 2) THEN
899  max7 = is - 2
900  ELSE
901  max7 = 1
902  END IF
903  IF (npx - 1 .GT. ie + 2) THEN
904  min7 = ie + 2
905  ELSE
906  min7 = npx - 1
907  END IF
908 !
909  DO i=max7,min7
910  g_in = dya(i, 2)/dya(i, 1)
911  g_ou = dya(i, -1)/dya(i, 0)
912  qy_tl(i, 1) = 0.5*(((2.+g_in)*qin_tl(i, 1)-qin_tl(i, 2))/(1.&
913 & +g_in)+((2.+g_ou)*qin_tl(i, 0)-qin_tl(i, -1))/(1.+g_ou))
914  qy(i, 1) = 0.5*(((2.+g_in)*qin(i, 1)-qin(i, 2))/(1.+g_in)+((&
915 & 2.+g_ou)*qin(i, 0)-qin(i, -1))/(1.+g_ou))
916  qy_tl(i, 2) = (3.*(g_in*qin_tl(i, 1)+qin_tl(i, 2))-g_in*&
917 & qy_tl(i, 1)-qy_tl(i, 3))/(2.+2.*g_in)
918  qy(i, 2) = (3.*(g_in*qin(i, 1)+qin(i, 2))-(g_in*qy(i, 1)+qy(&
919 & i, 3)))/(2.+2.*g_in)
920  END DO
921  ELSE
922  q1_tl = 0.0
923  END IF
924 ! North Edges:
925  IF (je + 1 .EQ. npy) THEN
926  DO i=is1,ie1
927  q1_tl(i) = (dya(i, npy)*qin_tl(i, npy-1)+dya(i, npy-1)*&
928 & qin_tl(i, npy))/(dya(i, npy-1)+dya(i, npy))
929  q1(i) = (qin(i, npy-1)*dya(i, npy)+qin(i, npy)*dya(i, npy-1)&
930 & )/(dya(i, npy-1)+dya(i, npy))
931  END DO
932  DO i=is2,ie1
933  qout_tl(i, npy) = edge_n(i)*q1_tl(i-1) + (1.-edge_n(i))*&
934 & q1_tl(i)
935  qout(i, npy) = edge_n(i)*q1(i-1) + (1.-edge_n(i))*q1(i)
936  END DO
937  IF (1 .LT. is - 2) THEN
938  max8 = is - 2
939  ELSE
940  max8 = 1
941  END IF
942  IF (npx - 1 .GT. ie + 2) THEN
943  min8 = ie + 2
944  ELSE
945  min8 = npx - 1
946  END IF
947 !
948  DO i=max8,min8
949  g_in = dya(i, npy-2)/dya(i, npy-1)
950  g_ou = dya(i, npy+1)/dya(i, npy)
951  qy_tl(i, npy) = 0.5*(((2.+g_in)*qin_tl(i, npy-1)-qin_tl(i, &
952 & npy-2))/(1.+g_in)+((2.+g_ou)*qin_tl(i, npy)-qin_tl(i, npy+&
953 & 1))/(1.+g_ou))
954  qy(i, npy) = 0.5*(((2.+g_in)*qin(i, npy-1)-qin(i, npy-2))/(&
955 & 1.+g_in)+((2.+g_ou)*qin(i, npy)-qin(i, npy+1))/(1.+g_ou))
956  qy_tl(i, npy-1) = (3.*(qin_tl(i, npy-2)+g_in*qin_tl(i, npy-1&
957 & ))-g_in*qy_tl(i, npy)-qy_tl(i, npy-2))/(2.+2.*g_in)
958  qy(i, npy-1) = (3.*(qin(i, npy-2)+g_in*qin(i, npy-1))-(g_in*&
959 & qy(i, npy)+qy(i, npy-2)))/(2.+2.*g_in)
960  END DO
961  END IF
962  END IF
963 !--------------------------------------
964  IF (gridstruct%nested) THEN
965  qxx_tl = 0.0
966  DO j=js,je+1
967  DO i=is,ie+1
968  qxx_tl(i, j) = a2*(qx_tl(i, j-2)+qx_tl(i, j+1)) + a1*(qx_tl(&
969 & i, j-1)+qx_tl(i, j))
970  qxx(i, j) = a2*(qx(i, j-2)+qx(i, j+1)) + a1*(qx(i, j-1)+qx(i&
971 & , j))
972  END DO
973  END DO
974  qyy_tl = 0.0
975  DO j=js,je+1
976  DO i=is,ie+1
977  qyy_tl(i, j) = a2*(qy_tl(i-2, j)+qy_tl(i+1, j)) + a1*(qy_tl(&
978 & i-1, j)+qy_tl(i, j))
979  qyy(i, j) = a2*(qy(i-2, j)+qy(i+1, j)) + a1*(qy(i-1, j)+qy(i&
980 & , j))
981  END DO
982  DO i=is,ie+1
983 ! averaging
984  qout_tl(i, j) = 0.5*(qxx_tl(i, j)+qyy_tl(i, j))
985  qout(i, j) = 0.5*(qxx(i, j)+qyy(i, j))
986  END DO
987  END DO
988  ELSE
989  IF (3 .LT. js) THEN
990  max9 = js
991  ELSE
992  max9 = 3
993  END IF
994  IF (npy - 2 .GT. je + 1) THEN
995  min9 = je + 1
996  qxx_tl = 0.0
997  ELSE
998  min9 = npy - 2
999  qxx_tl = 0.0
1000  END IF
1001  DO j=max9,min9
1002  IF (2 .LT. is) THEN
1003  max10 = is
1004  ELSE
1005  max10 = 2
1006  END IF
1007  IF (npx - 1 .GT. ie + 1) THEN
1008  min10 = ie + 1
1009  ELSE
1010  min10 = npx - 1
1011  END IF
1012  DO i=max10,min10
1013  qxx_tl(i, j) = a2*(qx_tl(i, j-2)+qx_tl(i, j+1)) + a1*(qx_tl(&
1014 & i, j-1)+qx_tl(i, j))
1015  qxx(i, j) = a2*(qx(i, j-2)+qx(i, j+1)) + a1*(qx(i, j-1)+qx(i&
1016 & , j))
1017  END DO
1018  END DO
1019  IF (js .EQ. 1) THEN
1020  IF (2 .LT. is) THEN
1021  max11 = is
1022  ELSE
1023  max11 = 2
1024  END IF
1025  IF (npx - 1 .GT. ie + 1) THEN
1026  min11 = ie + 1
1027  ELSE
1028  min11 = npx - 1
1029  END IF
1030  DO i=max11,min11
1031  qxx_tl(i, 2) = c1*(qx_tl(i, 1)+qx_tl(i, 2)) + c2*(qout_tl(i&
1032 & , 1)+qxx_tl(i, 3))
1033  qxx(i, 2) = c1*(qx(i, 1)+qx(i, 2)) + c2*(qout(i, 1)+qxx(i, 3&
1034 & ))
1035  END DO
1036  END IF
1037  IF (je + 1 .EQ. npy) THEN
1038  IF (2 .LT. is) THEN
1039  max12 = is
1040  ELSE
1041  max12 = 2
1042  END IF
1043  IF (npx - 1 .GT. ie + 1) THEN
1044  min12 = ie + 1
1045  ELSE
1046  min12 = npx - 1
1047  END IF
1048  DO i=max12,min12
1049  qxx_tl(i, npy-1) = c1*(qx_tl(i, npy-2)+qx_tl(i, npy-1)) + c2&
1050 & *(qout_tl(i, npy)+qxx_tl(i, npy-2))
1051  qxx(i, npy-1) = c1*(qx(i, npy-2)+qx(i, npy-1)) + c2*(qout(i&
1052 & , npy)+qxx(i, npy-2))
1053  END DO
1054  END IF
1055  IF (2 .LT. js) THEN
1056  max13 = js
1057  ELSE
1058  max13 = 2
1059  END IF
1060  IF (npy - 1 .GT. je + 1) THEN
1061  min13 = je + 1
1062  qyy_tl = 0.0
1063  ELSE
1064  min13 = npy - 1
1065  qyy_tl = 0.0
1066  END IF
1067  DO j=max13,min13
1068  IF (3 .LT. is) THEN
1069  max14 = is
1070  ELSE
1071  max14 = 3
1072  END IF
1073  IF (npx - 2 .GT. ie + 1) THEN
1074  min14 = ie + 1
1075  ELSE
1076  min14 = npx - 2
1077  END IF
1078  DO i=max14,min14
1079  qyy_tl(i, j) = a2*(qy_tl(i-2, j)+qy_tl(i+1, j)) + a1*(qy_tl(&
1080 & i-1, j)+qy_tl(i, j))
1081  qyy(i, j) = a2*(qy(i-2, j)+qy(i+1, j)) + a1*(qy(i-1, j)+qy(i&
1082 & , j))
1083  END DO
1084  IF (is .EQ. 1) THEN
1085  qyy_tl(2, j) = c1*(qy_tl(1, j)+qy_tl(2, j)) + c2*(qout_tl(1&
1086 & , j)+qyy_tl(3, j))
1087  qyy(2, j) = c1*(qy(1, j)+qy(2, j)) + c2*(qout(1, j)+qyy(3, j&
1088 & ))
1089  END IF
1090  IF (ie + 1 .EQ. npx) THEN
1091  qyy_tl(npx-1, j) = c1*(qy_tl(npx-2, j)+qy_tl(npx-1, j)) + c2&
1092 & *(qout_tl(npx, j)+qyy_tl(npx-2, j))
1093  qyy(npx-1, j) = c1*(qy(npx-2, j)+qy(npx-1, j)) + c2*(qout(&
1094 & npx, j)+qyy(npx-2, j))
1095  END IF
1096  IF (2 .LT. is) THEN
1097  max15 = is
1098  ELSE
1099  max15 = 2
1100  END IF
1101  IF (npx - 1 .GT. ie + 1) THEN
1102  min15 = ie + 1
1103  ELSE
1104  min15 = npx - 1
1105  END IF
1106  DO i=max15,min15
1107 ! averaging
1108  qout_tl(i, j) = 0.5*(qxx_tl(i, j)+qyy_tl(i, j))
1109  qout(i, j) = 0.5*(qxx(i, j)+qyy(i, j))
1110  END DO
1111  END DO
1112  END IF
1113  ELSE
1114  qx_tl = 0.0
1115 ! grid_type>=3
1116 !------------------------
1117 ! Doubly periodic domain:
1118 !------------------------
1119 ! X-sweep: PPM
1120  DO j=js-2,je+2
1121  DO i=is,ie+1
1122  qx_tl(i, j) = b1*(qin_tl(i-1, j)+qin_tl(i, j)) + b2*(qin_tl(i-&
1123 & 2, j)+qin_tl(i+1, j))
1124  qx(i, j) = b1*(qin(i-1, j)+qin(i, j)) + b2*(qin(i-2, j)+qin(i+&
1125 & 1, j))
1126  END DO
1127  END DO
1128  qy_tl = 0.0
1129 ! Y-sweep: PPM
1130  DO j=js,je+1
1131  DO i=is-2,ie+2
1132  qy_tl(i, j) = b1*(qin_tl(i, j-1)+qin_tl(i, j)) + b2*(qin_tl(i&
1133 & , j-2)+qin_tl(i, j+1))
1134  qy(i, j) = b1*(qin(i, j-1)+qin(i, j)) + b2*(qin(i, j-2)+qin(i&
1135 & , j+1))
1136  END DO
1137  END DO
1138  DO j=js,je+1
1139  DO i=is,ie+1
1140  qout_tl(i, j) = 0.5*(a1*(qx_tl(i, j-1)+qx_tl(i, j)+qy_tl(i-1, &
1141 & j)+qy_tl(i, j))+a2*(qx_tl(i, j-2)+qx_tl(i, j+1)+qy_tl(i-2, j&
1142 & )+qy_tl(i+1, j)))
1143  qout(i, j) = 0.5*(a1*(qx(i, j-1)+qx(i, j)+qy(i-1, j)+qy(i, j))&
1144 & +a2*(qx(i, j-2)+qx(i, j+1)+qy(i-2, j)+qy(i+1, j)))
1145  END DO
1146  END DO
1147  END IF
1148  IF (PRESENT(replace)) THEN
1149  IF (replace) THEN
1150  DO j=js,je+1
1151  DO i=is,ie+1
1152  qin_tl(i, j) = qout_tl(i, j)
1153  qin(i, j) = qout(i, j)
1154  END DO
1155  END DO
1156  END IF
1157  END IF
1158  END SUBROUTINE a2b_ord4_tlm
1159 ! Differentiation of a2b_ord2 in forward (tangent) mode:
1160 ! variations of useful results: qin qout
1161 ! with respect to varying inputs: qin qout
1162  SUBROUTINE a2b_ord2_tlm(qin, qin_tl, qout, qout_tl, gridstruct, npx, &
1163 & npy, is, ie, js, je, ng, replace)
1164  IMPLICIT NONE
1165  INTEGER, INTENT(IN) :: npx, npy, is, ie, js, je, ng
1166 ! A-grid field
1167  REAL, INTENT(INOUT) :: qin(is-ng:ie+ng, js-ng:je+ng)
1168  REAL, INTENT(INOUT) :: qin_tl(is-ng:ie+ng, js-ng:je+ng)
1169 ! Output B-grid field
1170  REAL, INTENT(OUT) :: qout(is-ng:ie+ng, js-ng:je+ng)
1171  REAL, INTENT(OUT) :: qout_tl(is-ng:ie+ng, js-ng:je+ng)
1172  TYPE(fv_grid_type), INTENT(IN), TARGET :: gridstruct
1173  LOGICAL, OPTIONAL, INTENT(IN) :: replace
1174 ! local:
1175  REAL :: q1(npx), q2(npy)
1176  REAL :: q1_tl(npx), q2_tl(npy)
1177  INTEGER :: i, j
1178  INTEGER :: is1, js1, is2, js2, ie1, je1
1179  REAL, DIMENSION(:, :, :), POINTER :: grid, agrid
1180  REAL, DIMENSION(:, :), POINTER :: dxa, dya
1181  REAL(kind=r_grid), DIMENSION(:), POINTER :: edge_w, edge_e, edge_s, &
1182 & edge_n
1183  INTRINSIC max
1184  INTRINSIC min
1185  INTRINSIC PRESENT
1186  edge_w => gridstruct%edge_w
1187  edge_e => gridstruct%edge_e
1188  edge_s => gridstruct%edge_s
1189  edge_n => gridstruct%edge_n
1190  grid => gridstruct%grid
1191  agrid => gridstruct%agrid
1192  dxa => gridstruct%dxa
1193  dya => gridstruct%dya
1194  IF (gridstruct%grid_type .LT. 3) THEN
1195  IF (gridstruct%nested) THEN
1196  DO j=js-2,je+1+2
1197  DO i=is-2,ie+1+2
1198  qout_tl(i, j) = 0.25*(qin_tl(i-1, j-1)+qin_tl(i, j-1)+qin_tl&
1199 & (i-1, j)+qin_tl(i, j))
1200  qout(i, j) = 0.25*(qin(i-1, j-1)+qin(i, j-1)+qin(i-1, j)+qin&
1201 & (i, j))
1202  END DO
1203  END DO
1204  ELSE
1205  IF (1 .LT. is - 1) THEN
1206  is1 = is - 1
1207  ELSE
1208  is1 = 1
1209  END IF
1210  IF (1 .LT. js - 1) THEN
1211  js1 = js - 1
1212  ELSE
1213  js1 = 1
1214  END IF
1215  IF (2 .LT. is) THEN
1216  is2 = is
1217  ELSE
1218  is2 = 2
1219  END IF
1220  IF (2 .LT. js) THEN
1221  js2 = js
1222  ELSE
1223  js2 = 2
1224  END IF
1225  IF (npx - 1 .GT. ie + 1) THEN
1226  ie1 = ie + 1
1227  ELSE
1228  ie1 = npx - 1
1229  END IF
1230  IF (npy - 1 .GT. je + 1) THEN
1231  je1 = je + 1
1232  ELSE
1233  je1 = npy - 1
1234  END IF
1235  DO j=js2,je1
1236  DO i=is2,ie1
1237  qout_tl(i, j) = 0.25*(qin_tl(i-1, j-1)+qin_tl(i, j-1)+qin_tl&
1238 & (i-1, j)+qin_tl(i, j))
1239  qout(i, j) = 0.25*(qin(i-1, j-1)+qin(i, j-1)+qin(i-1, j)+qin&
1240 & (i, j))
1241  END DO
1242  END DO
1243 ! Fix the 4 Corners:
1244  IF (gridstruct%sw_corner) THEN
1245  qout_tl(1, 1) = r3*(qin_tl(1, 1)+qin_tl(1, 0)+qin_tl(0, 1))
1246  qout(1, 1) = r3*(qin(1, 1)+qin(1, 0)+qin(0, 1))
1247  END IF
1248  IF (gridstruct%se_corner) THEN
1249  qout_tl(npx, 1) = r3*(qin_tl(npx-1, 1)+qin_tl(npx-1, 0)+qin_tl&
1250 & (npx, 1))
1251  qout(npx, 1) = r3*(qin(npx-1, 1)+qin(npx-1, 0)+qin(npx, 1))
1252  END IF
1253  IF (gridstruct%ne_corner) THEN
1254  qout_tl(npx, npy) = r3*(qin_tl(npx-1, npy-1)+qin_tl(npx, npy-1&
1255 & )+qin_tl(npx-1, npy))
1256  qout(npx, npy) = r3*(qin(npx-1, npy-1)+qin(npx, npy-1)+qin(npx&
1257 & -1, npy))
1258  END IF
1259  IF (gridstruct%nw_corner) THEN
1260  qout_tl(1, npy) = r3*(qin_tl(1, npy-1)+qin_tl(0, npy-1)+qin_tl&
1261 & (1, npy))
1262  qout(1, npy) = r3*(qin(1, npy-1)+qin(0, npy-1)+qin(1, npy))
1263  END IF
1264 ! *** West Edges:
1265  IF (is .EQ. 1) THEN
1266  q2_tl = 0.0
1267  DO j=js1,je1
1268  q2_tl(j) = 0.5*(qin_tl(0, j)+qin_tl(1, j))
1269  q2(j) = 0.5*(qin(0, j)+qin(1, j))
1270  END DO
1271  DO j=js2,je1
1272  qout_tl(1, j) = edge_w(j)*q2_tl(j-1) + (1.-edge_w(j))*q2_tl(&
1273 & j)
1274  qout(1, j) = edge_w(j)*q2(j-1) + (1.-edge_w(j))*q2(j)
1275  END DO
1276  ELSE
1277  q2_tl = 0.0
1278  END IF
1279 ! East Edges:
1280  IF (ie + 1 .EQ. npx) THEN
1281  DO j=js1,je1
1282  q2_tl(j) = 0.5*(qin_tl(npx-1, j)+qin_tl(npx, j))
1283  q2(j) = 0.5*(qin(npx-1, j)+qin(npx, j))
1284  END DO
1285  DO j=js2,je1
1286  qout_tl(npx, j) = edge_e(j)*q2_tl(j-1) + (1.-edge_e(j))*&
1287 & q2_tl(j)
1288  qout(npx, j) = edge_e(j)*q2(j-1) + (1.-edge_e(j))*q2(j)
1289  END DO
1290  END IF
1291 ! South Edges:
1292  IF (js .EQ. 1) THEN
1293  q1_tl = 0.0
1294  DO i=is1,ie1
1295  q1_tl(i) = 0.5*(qin_tl(i, 0)+qin_tl(i, 1))
1296  q1(i) = 0.5*(qin(i, 0)+qin(i, 1))
1297  END DO
1298  DO i=is2,ie1
1299  qout_tl(i, 1) = edge_s(i)*q1_tl(i-1) + (1.-edge_s(i))*q1_tl(&
1300 & i)
1301  qout(i, 1) = edge_s(i)*q1(i-1) + (1.-edge_s(i))*q1(i)
1302  END DO
1303  ELSE
1304  q1_tl = 0.0
1305  END IF
1306 ! North Edges:
1307  IF (je + 1 .EQ. npy) THEN
1308  DO i=is1,ie1
1309  q1_tl(i) = 0.5*(qin_tl(i, npy-1)+qin_tl(i, npy))
1310  q1(i) = 0.5*(qin(i, npy-1)+qin(i, npy))
1311  END DO
1312  DO i=is2,ie1
1313  qout_tl(i, npy) = edge_n(i)*q1_tl(i-1) + (1.-edge_n(i))*&
1314 & q1_tl(i)
1315  qout(i, npy) = edge_n(i)*q1(i-1) + (1.-edge_n(i))*q1(i)
1316  END DO
1317  END IF
1318  END IF
1319  ELSE
1320  DO j=js,je+1
1321  DO i=is,ie+1
1322  qout_tl(i, j) = 0.25*(qin_tl(i-1, j-1)+qin_tl(i, j-1)+qin_tl(i&
1323 & -1, j)+qin_tl(i, j))
1324  qout(i, j) = 0.25*(qin(i-1, j-1)+qin(i, j-1)+qin(i-1, j)+qin(i&
1325 & , j))
1326  END DO
1327  END DO
1328  END IF
1329  IF (PRESENT(replace)) THEN
1330  IF (replace) THEN
1331  DO j=js,je+1
1332  DO i=is,ie+1
1333  qin_tl(i, j) = qout_tl(i, j)
1334  qin(i, j) = qout(i, j)
1335  END DO
1336  END DO
1337  END IF
1338  END IF
1339  END SUBROUTINE a2b_ord2_tlm
1340  SUBROUTINE a2b_ord2(qin, qout, gridstruct, npx, npy, is, ie, js, je, &
1341 & ng, replace)
1342  IMPLICIT NONE
1343  INTEGER, INTENT(IN) :: npx, npy, is, ie, js, je, ng
1344 ! A-grid field
1345  REAL, INTENT(INOUT) :: qin(is-ng:ie+ng, js-ng:je+ng)
1346 ! Output B-grid field
1347  REAL, INTENT(OUT) :: qout(is-ng:ie+ng, js-ng:je+ng)
1348  TYPE(fv_grid_type), INTENT(IN), TARGET :: gridstruct
1349  LOGICAL, OPTIONAL, INTENT(IN) :: replace
1350 ! local:
1351  REAL :: q1(npx), q2(npy)
1352  INTEGER :: i, j
1353  INTEGER :: is1, js1, is2, js2, ie1, je1
1354  REAL, DIMENSION(:, :, :), POINTER :: grid, agrid
1355  REAL, DIMENSION(:, :), POINTER :: dxa, dya
1356  REAL(kind=r_grid), DIMENSION(:), POINTER :: edge_w, edge_e, edge_s, &
1357 & edge_n
1358  INTRINSIC max
1359  INTRINSIC min
1360  INTRINSIC PRESENT
1361  edge_w => gridstruct%edge_w
1362  edge_e => gridstruct%edge_e
1363  edge_s => gridstruct%edge_s
1364  edge_n => gridstruct%edge_n
1365  grid => gridstruct%grid
1366  agrid => gridstruct%agrid
1367  dxa => gridstruct%dxa
1368  dya => gridstruct%dya
1369  IF (gridstruct%grid_type .LT. 3) THEN
1370  IF (gridstruct%nested) THEN
1371  DO j=js-2,je+1+2
1372  DO i=is-2,ie+1+2
1373  qout(i, j) = 0.25*(qin(i-1, j-1)+qin(i, j-1)+qin(i-1, j)+qin&
1374 & (i, j))
1375  END DO
1376  END DO
1377  ELSE
1378  IF (1 .LT. is - 1) THEN
1379  is1 = is - 1
1380  ELSE
1381  is1 = 1
1382  END IF
1383  IF (1 .LT. js - 1) THEN
1384  js1 = js - 1
1385  ELSE
1386  js1 = 1
1387  END IF
1388  IF (2 .LT. is) THEN
1389  is2 = is
1390  ELSE
1391  is2 = 2
1392  END IF
1393  IF (2 .LT. js) THEN
1394  js2 = js
1395  ELSE
1396  js2 = 2
1397  END IF
1398  IF (npx - 1 .GT. ie + 1) THEN
1399  ie1 = ie + 1
1400  ELSE
1401  ie1 = npx - 1
1402  END IF
1403  IF (npy - 1 .GT. je + 1) THEN
1404  je1 = je + 1
1405  ELSE
1406  je1 = npy - 1
1407  END IF
1408  DO j=js2,je1
1409  DO i=is2,ie1
1410  qout(i, j) = 0.25*(qin(i-1, j-1)+qin(i, j-1)+qin(i-1, j)+qin&
1411 & (i, j))
1412  END DO
1413  END DO
1414 ! Fix the 4 Corners:
1415  IF (gridstruct%sw_corner) qout(1, 1) = r3*(qin(1, 1)+qin(1, 0)+&
1416 & qin(0, 1))
1417  IF (gridstruct%se_corner) qout(npx, 1) = r3*(qin(npx-1, 1)+qin(&
1418 & npx-1, 0)+qin(npx, 1))
1419  IF (gridstruct%ne_corner) qout(npx, npy) = r3*(qin(npx-1, npy-1)&
1420 & +qin(npx, npy-1)+qin(npx-1, npy))
1421  IF (gridstruct%nw_corner) qout(1, npy) = r3*(qin(1, npy-1)+qin(0&
1422 & , npy-1)+qin(1, npy))
1423 ! *** West Edges:
1424  IF (is .EQ. 1) THEN
1425  DO j=js1,je1
1426  q2(j) = 0.5*(qin(0, j)+qin(1, j))
1427  END DO
1428  DO j=js2,je1
1429  qout(1, j) = edge_w(j)*q2(j-1) + (1.-edge_w(j))*q2(j)
1430  END DO
1431  END IF
1432 ! East Edges:
1433  IF (ie + 1 .EQ. npx) THEN
1434  DO j=js1,je1
1435  q2(j) = 0.5*(qin(npx-1, j)+qin(npx, j))
1436  END DO
1437  DO j=js2,je1
1438  qout(npx, j) = edge_e(j)*q2(j-1) + (1.-edge_e(j))*q2(j)
1439  END DO
1440  END IF
1441 ! South Edges:
1442  IF (js .EQ. 1) THEN
1443  DO i=is1,ie1
1444  q1(i) = 0.5*(qin(i, 0)+qin(i, 1))
1445  END DO
1446  DO i=is2,ie1
1447  qout(i, 1) = edge_s(i)*q1(i-1) + (1.-edge_s(i))*q1(i)
1448  END DO
1449  END IF
1450 ! North Edges:
1451  IF (je + 1 .EQ. npy) THEN
1452  DO i=is1,ie1
1453  q1(i) = 0.5*(qin(i, npy-1)+qin(i, npy))
1454  END DO
1455  DO i=is2,ie1
1456  qout(i, npy) = edge_n(i)*q1(i-1) + (1.-edge_n(i))*q1(i)
1457  END DO
1458  END IF
1459  END IF
1460  ELSE
1461  DO j=js,je+1
1462  DO i=is,ie+1
1463  qout(i, j) = 0.25*(qin(i-1, j-1)+qin(i, j-1)+qin(i-1, j)+qin(i&
1464 & , j))
1465  END DO
1466  END DO
1467  END IF
1468  IF (PRESENT(replace)) THEN
1469  IF (replace) THEN
1470  DO j=js,je+1
1471  DO i=is,ie+1
1472  qin(i, j) = qout(i, j)
1473  END DO
1474  END DO
1475  END IF
1476  END IF
1477  END SUBROUTINE a2b_ord2
1478  REAL FUNCTION extrap_corner(p0, p1, p2, q1, q2)
1479  IMPLICIT NONE
1480  REAL, DIMENSION(2), INTENT(IN) :: p0, p1, p2
1481  REAL, INTENT(IN) :: q1, q2
1482  REAL :: x1, x2
1483  INTRINSIC real
1484  x1 = great_circle_dist(REAL(p1, kind=r_grid), REAL(p0, kind=r_grid))
1485  x2 = great_circle_dist(REAL(p2, kind=r_grid), REAL(p0, kind=r_grid))
1486  extrap_corner = q1 + x1/(x2-x1)*(q1-q2)
1487  END FUNCTION extrap_corner
1488 ! Differentiation of extrap_corner in forward (tangent) mode:
1489 ! variations of useful results: extrap_corner
1490 ! with respect to varying inputs: q1 q2
1491  REAL FUNCTION extrap_corner_tlm(p0, p1, p2, q1, q1_tl, q2, q2_tl, &
1492 & extrap_corner)
1493  IMPLICIT NONE
1494  REAL, DIMENSION(2), INTENT(IN) :: p0, p1, p2
1495  REAL, INTENT(IN) :: q1, q2
1496  REAL, INTENT(IN) :: q1_tl, q2_tl
1497  REAL :: x1, x2
1498  INTRINSIC real
1499  REAL :: extrap_corner
1500  x1 = great_circle_dist(REAL(p1, kind=r_grid), REAL(p0, kind=r_grid))
1501  x2 = great_circle_dist(REAL(p2, kind=r_grid), REAL(p0, kind=r_grid))
1502  extrap_corner_tlm = q1_tl + x1*(q1_tl-q2_tl)/(x2-x1)
1503  extrap_corner = q1 + x1/(x2-x1)*(q1-q2)
1504  END FUNCTION extrap_corner_tlm
1505 end module a2b_edge_tlm_mod
a2b_edge_tlm_mod
Definition: a2b_edge_tlm.F90:20
fv_arrays_nlm_mod::fv_grid_type
Definition: fv_arrays_nlm.F90:115
a2b_edge_tlm_mod::b1
real, parameter b1
Definition: a2b_edge_tlm.F90:40
a2b_edge_tlm_mod::extrap_corner
real function extrap_corner(p0, p1, p2, q1, q2)
Definition: a2b_edge_tlm.F90:1479
a2b_edge_tlm_mod::a2b_ord4
subroutine, public a2b_ord4(qin, qout, gridstruct, npx, npy, is, ie, js, je, ng, replace)
Definition: a2b_edge_tlm.F90:50
a2b_edge_tlm_mod::extrap_corner_tlm
real function extrap_corner_tlm(p0, p1, p2, q1, q1_tl, q2, q2_tl, extrap_corner)
Definition: a2b_edge_tlm.F90:1493
fv_grid_utils_nlm_mod
Definition: fv_grid_utils_nlm.F90:20
a2b_ord2
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
fv_arrays_nlm_mod
Definition: fv_arrays_nlm.F90:20
fv_arrays_nlm_mod::r_grid
integer, parameter, public r_grid
Definition: fv_arrays_nlm.F90:32
fv_grid_utils_nlm_mod::great_circle_dist
real function, public great_circle_dist(q1, q2, radius)
Definition: fv_grid_utils_nlm.F90:1968
a2b_edge_tlm_mod::a2b_ord4_tlm
subroutine, public a2b_ord4_tlm(qin, qin_tl, qout, qout_tl, gridstruct, npx, npy, is, ie, js, je, ng, replace)
Definition: a2b_edge_tlm.F90:548
max
#define max(a, b)
Definition: mosaic_util.h:33
a2b_edge_tlm_mod::a1
real, parameter a1
Definition: a2b_edge_tlm.F90:35
min
#define min(a, b)
Definition: mosaic_util.h:32
a2b_edge_tlm_mod::b2
real, parameter b2
Definition: a2b_edge_tlm.F90:41
a2b_edge_tlm_mod::a2
real, parameter a2
Definition: a2b_edge_tlm.F90:36
a2b_edge_tlm_mod::r3
real, parameter r3
Definition: a2b_edge_tlm.F90:31
a2b_edge_tlm_mod::a2b_ord2_tlm
subroutine, public a2b_ord2_tlm(qin, qin_tl, qout, qout_tl, gridstruct, npx, npy, is, ie, js, je, ng, replace)
Definition: a2b_edge_tlm.F90:1164