FV3 Bundle
advect_pv_ad.f90
Go to the documentation of this file.
1 ! (C) Copyright 2009-2016 ECMWF.
2 !
3 ! This software is licensed under the terms of the Apache Licence Version 2.0
4 ! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
5 ! In applying this licence, ECMWF does not waive the privileges and immunities
6 ! granted to it by virtue of its status as an intergovernmental organisation nor
7 ! does it submit to any jurisdiction.
8 
9 !> Advect potential vorticity - Adjoint
10 
11 subroutine advect_pv_ad (qnew,q,q_traj,q_north,q_south, &
12  & u,u_traj,v,v_traj, &
13  & nx,ny,deltax,deltay,dt)
14 
15 use kinds
16 implicit none
17 
18 real(kind=kind_real), intent(inout) :: qnew(nx,ny,2) !< Input adjoint variable for PV
19 real(kind=kind_real), intent(inout) :: q(nx,ny,2) !< Output adjoint variable for PV
20 real(kind=kind_real), intent(in) :: q_traj(nx,ny,2) !< Trajectory potential vorticity
21 real(kind=kind_real), intent(in) :: q_north(nx,2) !< PV on northern wall
22 real(kind=kind_real), intent(in) :: q_south(nx,2) !< PV on southern wall
23 real(kind=kind_real), intent(out) :: u(nx,ny,2) !< Zonal wind adjoint variable
24 real(kind=kind_real), intent(in) :: u_traj(nx,ny,2) !< Trajectory zonal wind
25 real(kind=kind_real), intent(out) :: v(nx,ny,2) !< Meridional wind adjoint variable
26 real(kind=kind_real), intent(in) :: v_traj(nx,ny,2) !< Trajectory meridional wind
27 integer, intent(in) :: nx !< Zonal grid dimension
28 integer, intent(in) :: ny !< Meridional grid dimension
29 real(kind=kind_real), intent(in) :: deltax !< Zonal grid spacing (non-dimensional)
30 real(kind=kind_real), intent(in) :: deltay !< Meridional grid spacing (non-dim)
31 real(kind=kind_real), intent(in) :: dt !< Timestep (non-dimensional)
32 
33 integer :: ii,jj,kk,ixm1,ix,ixp1,ixp2,jym1,jy,jyp1,jyp2
34 real(kind=kind_real) :: ax_traj,ay_traj,qjm1_traj,qj_traj,qjp1_traj,qjp2_traj
35 real(kind=kind_real) :: ax,ay,qjm1,qj,qjp1,qjp2
36 real(kind_real), parameter :: one=1.0_kind_real
37 real(kind_real), parameter :: two=2.0_kind_real
38 real(kind_real), parameter :: half=0.5_kind_real
39 real(kind_real), parameter :: sixth=1.0_kind_real/6.0_kind_real
40 
41 !--- advect
42 
43 do kk=1,2
44  do jj=1,ny
45  do ii=1,nx
46  q(ii,jj,kk) = 0.0_kind_real
47  enddo
48  enddo
49 enddo
50 
51 do kk=1,2
52  do jj=1,ny
53  do ii=1,nx
54 
55 !==== Trajectory calculations
56 
57 !--- find the interpolation point (nonlinear)
58 
59  ax_traj = real(ii,kind_real) - u_traj(ii,jj,kk)*dt/deltax
60  ix = floor(ax_traj)
61  ax_traj = ax_traj-real(ix,kind_real)
62  ixm1 = 1 + modulo(ix-2,nx)
63  ixp1 = 1 + modulo(ix ,nx)
64  ixp2 = 1 + modulo(ix+1,nx)
65  ix = 1 + modulo(ix-1,nx)
66 
67  ay_traj = real(jj,kind_real) - v_traj(ii,jj,kk)*dt/deltay
68  jy = floor(ay_traj)
69  ay_traj = ay_traj-real(jy,kind_real)
70  jym1 = jy-1
71  jyp1 = jy+1
72  jyp2 = jy+2
73 
74  if (jym1 < 1) then
75  qjm1_traj = &
76  & ax_traj *(ax_traj-one)*(ax_traj-two) &
77  & *q_south(ixm1,kk)*(-sixth) + &
78  & (ax_traj+one)*(ax_traj-one)*(ax_traj-two) &
79  & *q_south(ix, kk)*half + &
80  & (ax_traj+one)* ax_traj *(ax_traj-two) &
81  & *q_south(ixp1,kk)*(-half) + &
82  & (ax_traj+one) * ax_traj *(ax_traj-one) &
83  & *q_south(ixp2,kk)*(sixth)
84  else if (jym1 > ny) then
85  qjm1_traj = &
86  & ax_traj *(ax_traj-one)*(ax_traj-two) &
87  & *q_north(ixm1,kk)*(-sixth) + &
88  & (ax_traj+one)*(ax_traj-one)*(ax_traj-two) &
89  & *q_north(ix, kk)*half + &
90  & (ax_traj+one)* ax_traj *(ax_traj-two) &
91  & *q_north(ixp1,kk)*(-half) + &
92  & (ax_traj+one)* ax_traj *(ax_traj-one) &
93  & *q_north(ixp2,kk)*(sixth)
94  else
95  qjm1_traj = &
96  & ax_traj *(ax_traj-one)*(ax_traj-two) &
97  & *q_traj(ixm1,jym1,kk)*(-sixth) + &
98  & (ax_traj+one)*(ax_traj-one)*(ax_traj-two) &
99  & *q_traj(ix, jym1,kk)*half + &
100  & (ax_traj+one)* ax_traj *(ax_traj-two) &
101  & *q_traj(ixp1,jym1,kk)*(-half) + &
102  & (ax_traj+one)* ax_traj *(ax_traj-one) &
103  & *q_traj(ixp2,jym1,kk)*(sixth)
104  endif
105 
106  if (jy < 1) then
107  qj_traj = &
108  & ax_traj *(ax_traj-one)*(ax_traj-two) &
109  & *q_south(ixm1,kk)*(-sixth) + &
110  & (ax_traj+one)*(ax_traj-one)*(ax_traj-two) &
111  & *q_south(ix ,kk)*half + &
112  & (ax_traj+one)* ax_traj *(ax_traj-two) &
113  & *q_south(ixp1,kk)*(-half) + &
114  & (ax_traj+one)* ax_traj *(ax_traj-one) &
115  & *q_south(ixp2,kk)*(sixth)
116  else if (jy > ny) then
117  qj_traj = &
118  & ax_traj *(ax_traj-one)*(ax_traj-two) &
119  & *q_north(ixm1,kk)*(-sixth) + &
120  & (ax_traj+one)*(ax_traj-one)*(ax_traj-two) &
121  & *q_north(ix ,kk)*half + &
122  & (ax_traj+one)* ax_traj *(ax_traj-two) &
123  & *q_north(ixp1,kk)*(-half) + &
124  & (ax_traj+one)* ax_traj *(ax_traj-one) &
125  & *q_north(ixp2,kk)*(sixth)
126  else
127  qj_traj = &
128  & ax_traj *(ax_traj-one)*(ax_traj-two) &
129  & *q_traj(ixm1,jy,kk)*(-sixth) + &
130  & (ax_traj+one)*(ax_traj-one)*(ax_traj-two) &
131  & *q_traj(ix ,jy,kk)*half + &
132  & (ax_traj+one)* ax_traj *(ax_traj-two) &
133  & *q_traj(ixp1,jy,kk)*(-half) + &
134  & (ax_traj+one)* ax_traj *(ax_traj-one) &
135  & *q_traj(ixp2,jy,kk)*(sixth)
136  endif
137 
138  if (jyp1 < 1) then
139  qjp1_traj = &
140  & ax_traj *(ax_traj-one)*(ax_traj-two) &
141  & *q_south(ixm1,kk)*(-sixth) + &
142  & (ax_traj+one)*(ax_traj-one)*(ax_traj-two) &
143  & *q_south(ix ,kk)*half + &
144  & (ax_traj+one)* ax_traj *(ax_traj-two) &
145  & *q_south(ixp1,kk)*(-half) + &
146  & (ax_traj+one)* ax_traj *(ax_traj-one) &
147  & *q_south(ixp2,kk)*(sixth)
148  else if (jyp1 > ny) then
149  qjp1_traj = &
150  & ax_traj *(ax_traj-one)*(ax_traj-two) &
151  & *q_north(ixm1,kk)*(-sixth) + &
152  & (ax_traj+one)*(ax_traj-one)*(ax_traj-two) &
153  & *q_north(ix ,kk)*half + &
154  & (ax_traj+one)* ax_traj *(ax_traj-two) &
155  & *q_north(ixp1,kk)*(-half) + &
156  & (ax_traj+one)* ax_traj *(ax_traj-one) &
157  & *q_north(ixp2,kk)*(sixth)
158  else
159  qjp1_traj = &
160  & ax_traj *(ax_traj-one)*(ax_traj-two) &
161  & *q_traj(ixm1,jyp1,kk)*(-sixth) + &
162  & (ax_traj+one)*(ax_traj-one)*(ax_traj-two) &
163  & *q_traj(ix ,jyp1,kk)*half + &
164  & (ax_traj+one)* ax_traj *(ax_traj-two) &
165  & *q_traj(ixp1,jyp1,kk)*(-half) + &
166  & (ax_traj+one)* ax_traj *(ax_traj-one) &
167  & *q_traj(ixp2,jyp1,kk)*(sixth)
168  endif
169 
170  if (jyp2 < 1) then
171  qjp2_traj = &
172  & ax_traj *(ax_traj-one)*(ax_traj-two) &
173  & *q_south(ixm1,kk)*(-sixth) + &
174  & (ax_traj+one)*(ax_traj-one)*(ax_traj-two) &
175  & *q_south(ix ,kk)*half + &
176  & (ax_traj+one)* ax_traj *(ax_traj-two) &
177  & *q_south(ixp1,kk)*(-half) + &
178  & (ax_traj+one)* ax_traj *(ax_traj-one) &
179  & *q_south(ixp2,kk)*(sixth)
180  else if (jyp2 > ny) then
181  qjp2_traj = &
182  & ax_traj *(ax_traj-one)*(ax_traj-two) &
183  & *q_north(ixm1,kk)*(-sixth) + &
184  & (ax_traj+one)*(ax_traj-one)*(ax_traj-two) &
185  & *q_north(ix ,kk)*half + &
186  & (ax_traj+one)* ax_traj *(ax_traj-two) &
187  & *q_north(ixp1,kk)*(-half) + &
188  & (ax_traj+one)* ax_traj *(ax_traj-one) &
189  & *q_north(ixp2,kk)*(sixth)
190  else
191  qjp2_traj = &
192  & ax_traj *(ax_traj-one)*(ax_traj-two) &
193  & *q_traj(ixm1,jyp2,kk)*(-sixth) + &
194  & (ax_traj+one)*(ax_traj-one)*(ax_traj-two) &
195  & *q_traj(ix ,jyp2,kk)*half + &
196  & (ax_traj+one)* ax_traj *(ax_traj-two) &
197  & *q_traj(ixp1,jyp2,kk)*(-half) + &
198  & (ax_traj+one)* ax_traj *(ax_traj-one) &
199  &*q_traj(ixp2,jyp2,kk)*(sixth)
200  endif
201 
202 !--- Lagrange interpolation in the meridional direction
203 
204 ! qnew_traj(ii,jj,kk) = &
205 ! & ay_traj *(ay_traj-one)*(ay_traj-two)*(-sixth)*qjm1_traj + &
206 ! & (ay_traj+one)*(ay_traj-one)*(ay_traj-two)*half *qj_traj + &
207 ! & (ay_traj+one)* ay_traj *(ay_traj-two)*(-half) *qjp1_traj + &
208 ! & (ay_traj+one)* ay_traj *(ay_traj-one)*(sixth) *qjp2_traj
209 
210 !==== Adjoint calculations
211 
212 !--- Lagrange interpolation in the meridional direction
213 
214  ay = ( &
215  & (ay_traj-one)*(ay_traj-two)*(-sixth)*qjm1_traj + &
216  & ay_traj * (ay_traj-two)*(-sixth)*qjm1_traj + &
217  & ay_traj *(ay_traj-one)* (-sixth)*qjm1_traj + &
218  & (ay_traj-one)*(ay_traj-two)*half *qj_traj + &
219  & (ay_traj+one)* (ay_traj-two)*half *qj_traj + &
220  & (ay_traj+one)*(ay_traj-one)* half *qj_traj + &
221  & ay_traj *(ay_traj-two)*(-half) *qjp1_traj + &
222  & (ay_traj+one)* (ay_traj-two)*(-half) *qjp1_traj + &
223  & (ay_traj+one)* ay_traj * (-half) *qjp1_traj + &
224  & ay_traj *(ay_traj-one)*(sixth) *qjp2_traj + &
225  & (ay_traj+one)* (ay_traj-one)*(sixth) *qjp2_traj + &
226  & (ay_traj+one)* ay_traj * (sixth) *qjp2_traj ) &
227  &*qnew(ii,jj,kk)
228 
229  qjm1 = ( ay_traj *(ay_traj-one)*(ay_traj-two)*(-sixth) ) &
230  & *qnew(ii,jj,kk)
231 
232  qj = ( (ay_traj+one)*(ay_traj-one)*(ay_traj-two)*half ) &
233  & *qnew(ii,jj,kk)
234 
235  qjp1 = ( (ay_traj+one)* ay_traj *(ay_traj-two)*(-half) ) &
236  & *qnew(ii,jj,kk)
237 
238  qjp2 = ( (ay_traj+one)* ay_traj *(ay_traj-one)*(sixth) ) &
239  & *qnew(ii,jj,kk)
240 
241  qnew(ii,jj,kk) = 0.0_kind_real
242 
243 !--- Lagrange interpolation in the zonal direction
244 
245  if (jyp2 < 1) then
246  ax = ( &
247  & (ax_traj-one)*(ax_traj-two) &
248  & *q_south(ixm1,kk)*(-sixth) + &
249  & ax_traj * (ax_traj-two) &
250  & *q_south(ixm1,kk)*(-sixth) + &
251  & ax_traj *(ax_traj-one) &
252  & *q_south(ixm1,kk)*(-sixth) + &
253  & (ax_traj-one)*(ax_traj-two) &
254  & *q_south(ix ,kk)*half + &
255  & (ax_traj+one)* (ax_traj-two) &
256  & *q_south(ix ,kk)*half + &
257  & (ax_traj+one)*(ax_traj-one) &
258  & *q_south(ix ,kk)*half + &
259  & ax_traj *(ax_traj-two) &
260  & *q_south(ixp1,kk)*(-half) + &
261  & (ax_traj+one)* (ax_traj-two) &
262  & *q_south(ixp1,kk)*(-half) + &
263  & (ax_traj+one)* ax_traj &
264  & *q_south(ixp1,kk)*(-half) + &
265  & ax_traj *(ax_traj-one) &
266  & *q_south(ixp2,kk)*(sixth) + &
267  & (ax_traj+one)* (ax_traj-one) &
268  & *q_south(ixp2,kk)*(sixth) + &
269  & (ax_traj+one)* ax_traj &
270  & *q_south(ixp2,kk)*(sixth) &
271  &)*qjp2
272 
273  qjp2 = 0.0_kind_real
274  else if (jyp2 > ny) then
275  ax = ( &
276  & (ax_traj-one)*(ax_traj-two) &
277  & *q_north(ixm1,kk)*(-sixth) + &
278  & ax_traj * (ax_traj-two) &
279  & *q_north(ixm1,kk)*(-sixth) + &
280  & ax_traj *(ax_traj-one) &
281  & *q_north(ixm1,kk)*(-sixth) + &
282  & (ax_traj-one)*(ax_traj-two) &
283  & *q_north(ix ,kk)*half + &
284  & (ax_traj+one)* (ax_traj-two) &
285  & *q_north(ix ,kk)*half + &
286  & (ax_traj+one)*(ax_traj-one) &
287  & *q_north(ix ,kk)*half + &
288  & ax_traj *(ax_traj-two) &
289  & *q_north(ixp1,kk)*(-half) + &
290  & (ax_traj+one)* (ax_traj-two) &
291  & *q_north(ixp1,kk)*(-half) + &
292  & (ax_traj+one)* ax_traj &
293  & *q_north(ixp1,kk)*(-half) + &
294  & ax_traj *(ax_traj-one) &
295  & *q_north(ixp2,kk)*(sixth) + &
296  & (ax_traj+one)* (ax_traj-one) &
297  & *q_north(ixp2,kk)*(sixth) + &
298  & (ax_traj+one)* ax_traj &
299  & *q_north(ixp2,kk)*(sixth) &
300  &)*qjp2
301 
302  qjp2 = 0.0_kind_real
303  else
304  ax = ( &
305  & (ax_traj-one)*(ax_traj-two) &
306  & *q_traj(ixm1,jyp2,kk)*(-sixth) + &
307  & ax_traj * (ax_traj-two) &
308  & *q_traj(ixm1,jyp2,kk)*(-sixth) + &
309  & ax_traj *(ax_traj-one) &
310  & *q_traj(ixm1,jyp2,kk)*(-sixth) + &
311  & (ax_traj-one)*(ax_traj-two) &
312  & *q_traj(ix ,jyp2,kk)*half + &
313  & (ax_traj+one)* (ax_traj-two) &
314  & *q_traj(ix ,jyp2,kk)*half + &
315  & (ax_traj+one)*(ax_traj-one) &
316  & *q_traj(ix ,jyp2,kk)*half + &
317  & ax_traj *(ax_traj-two) &
318  & *q_traj(ixp1,jyp2,kk)*(-half) + &
319  & (ax_traj+one)* (ax_traj-two) &
320  & *q_traj(ixp1,jyp2,kk)*(-half) + &
321  & (ax_traj+one)* ax_traj &
322  & *q_traj(ixp1,jyp2,kk)*(-half) + &
323  & ax_traj *(ax_traj-one) &
324  &*q_traj(ixp2,jyp2,kk)*(sixth) + &
325  & (ax_traj+one)* (ax_traj-one) &
326  &*q_traj(ixp2,jyp2,kk)*(sixth) + &
327  & (ax_traj+one)* ax_traj &
328  &*q_traj(ixp2,jyp2,kk)*(sixth) &
329  &)*qjp2
330 
331  q(ixm1,jyp2,kk) = q(ixm1,jyp2,kk) + ( &
332  & ax_traj *(ax_traj-one)*(ax_traj-two)*(-sixth) &
333  &)*qjp2
334 
335  q(ix ,jyp2,kk) = q(ix ,jyp2,kk) + ( &
336  & (ax_traj+one)*(ax_traj-one)*(ax_traj-two)*half &
337  &)*qjp2
338 
339  q(ixp1,jyp2,kk) = q(ixp1,jyp2,kk) + ( &
340  & (ax_traj+one)* ax_traj *(ax_traj-two)*(-half) &
341  &)*qjp2
342 
343  q(ixp2,jyp2,kk) = q(ixp2,jyp2,kk) + ( &
344  & (ax_traj+one)* ax_traj *(ax_traj-one)*(sixth) &
345  &)*qjp2
346 
347  qjp2 = 0.0_kind_real
348  endif
349 
350  if (jyp1 < 1) then
351  ax = ax + ( &
352  & (ax_traj-one)*(ax_traj-two) &
353  & *q_south(ixm1,kk)*(-sixth) + &
354  & ax_traj * (ax_traj-two) &
355  & *q_south(ixm1,kk)*(-sixth) + &
356  & ax_traj *(ax_traj-one) &
357  & *q_south(ixm1,kk)*(-sixth) + &
358  & (ax_traj-one)*(ax_traj-two) &
359  & *q_south(ix ,kk)*half + &
360  & (ax_traj+one)* (ax_traj-two) &
361  & *q_south(ix ,kk)*half + &
362  & (ax_traj+one)*(ax_traj-one) &
363  & *q_south(ix ,kk)*half + &
364  & ax_traj *(ax_traj-two) &
365  & *q_south(ixp1,kk)*(-half) + &
366  & (ax_traj+one)* (ax_traj-two) &
367  & *q_south(ixp1,kk)*(-half) + &
368  & (ax_traj+one)* ax_traj &
369  & *q_south(ixp1,kk)*(-half) + &
370  & ax_traj *(ax_traj-one) &
371  & *q_south(ixp2,kk)*(sixth) + &
372  & (ax_traj+one)* (ax_traj-one) &
373  & *q_south(ixp2,kk)*(sixth) + &
374  & (ax_traj+one)* ax_traj &
375  & *q_south(ixp2,kk)*(sixth) &
376  &)*qjp1
377 
378  qjp1 = 0.0_kind_real
379  else if (jyp1 > ny) then
380  ax = ax + ( &
381  & (ax_traj-one)*(ax_traj-two) &
382  & *q_north(ixm1,kk)*(-sixth) + &
383  & ax_traj * (ax_traj-two) &
384  & *q_north(ixm1,kk)*(-sixth) + &
385  & ax_traj *(ax_traj-one) &
386  & *q_north(ixm1,kk)*(-sixth) + &
387  & (ax_traj-one)*(ax_traj-two) &
388  & *q_north(ix ,kk)*half + &
389  & (ax_traj+one)* (ax_traj-two) &
390  & *q_north(ix ,kk)*half + &
391  & (ax_traj+one)*(ax_traj-one) &
392  & *q_north(ix ,kk)*half + &
393  & ax_traj *(ax_traj-two) &
394  & *q_north(ixp1,kk)*(-half) + &
395  & (ax_traj+one)* (ax_traj-two) &
396  & *q_north(ixp1,kk)*(-half) + &
397  & (ax_traj+one)* ax_traj &
398  & *q_north(ixp1,kk)*(-half) + &
399  & ax_traj *(ax_traj-one) &
400  & *q_north(ixp2,kk)*(sixth) + &
401  & (ax_traj+one)* (ax_traj-one) &
402  & *q_north(ixp2,kk)*(sixth) + &
403  & (ax_traj+one)* ax_traj &
404  & *q_north(ixp2,kk)*(sixth) &
405  &)*qjp1
406 
407  qjp1 = 0.0_kind_real
408  else
409  ax = ax + ( &
410  & (ax_traj-one)*(ax_traj-two) &
411  & *q_traj(ixm1,jyp1,kk)*(-sixth) + &
412  & ax_traj * (ax_traj-two) &
413  & *q_traj(ixm1,jyp1,kk)*(-sixth) + &
414  & ax_traj *(ax_traj-one) &
415  & *q_traj(ixm1,jyp1,kk)*(-sixth) + &
416  & (ax_traj-one)*(ax_traj-two) &
417  & *q_traj(ix ,jyp1,kk)*half + &
418  & (ax_traj+one)* (ax_traj-two) &
419  & *q_traj(ix ,jyp1,kk)*half + &
420  & (ax_traj+one)*(ax_traj-one) &
421  & *q_traj(ix ,jyp1,kk)*half + &
422  & ax_traj *(ax_traj-two) &
423  & *q_traj(ixp1,jyp1,kk)*(-half) + &
424  & (ax_traj+one)* (ax_traj-two) &
425  & *q_traj(ixp1,jyp1,kk)*(-half) + &
426  & (ax_traj+one)* ax_traj &
427  & *q_traj(ixp1,jyp1,kk)*(-half) + &
428  & ax_traj *(ax_traj-one) &
429  & *q_traj(ixp2,jyp1,kk)*(sixth) + &
430  & (ax_traj+one)* (ax_traj-one) &
431  & *q_traj(ixp2,jyp1,kk)*(sixth) + &
432  & (ax_traj+one)* ax_traj &
433  & *q_traj(ixp2,jyp1,kk)*(sixth) &
434  &)*qjp1
435 
436  q(ixm1,jyp1,kk) = q(ixm1,jyp1,kk) + ( &
437  & ax_traj *(ax_traj-one)*(ax_traj-two)*(-sixth) &
438  &)*qjp1
439 
440 
441  q(ix ,jyp1,kk) = q(ix ,jyp1,kk) + ( &
442  & (ax_traj+one)*(ax_traj-one)*(ax_traj-two)*half &
443  &)*qjp1
444 
445 
446  q(ixp1,jyp1,kk) = q(ixp1,jyp1,kk) + ( &
447  & (ax_traj+one)* ax_traj *(ax_traj-two)*(-half) &
448  &)*qjp1
449 
450  q(ixp2,jyp1,kk) = q(ixp2,jyp1,kk) + ( &
451  & (ax_traj+one)* ax_traj *(ax_traj-one)*(sixth) &
452  &)*qjp1
453  endif
454 
455  if (jy < 1) then
456  ax = ax + ( &
457  & (ax_traj-one)*(ax_traj-two) &
458  & *q_south(ixm1,kk)*(-sixth) + &
459  & ax_traj * (ax_traj-two) &
460  & *q_south(ixm1,kk)*(-sixth) + &
461  & ax_traj *(ax_traj-one) &
462  & *q_south(ixm1,kk)*(-sixth) + &
463  & (ax_traj-one)*(ax_traj-two) &
464  & *q_south(ix ,kk)*half + &
465  & (ax_traj+one)* (ax_traj-two) &
466  & *q_south(ix ,kk)*half + &
467  & (ax_traj+one)*(ax_traj-one) &
468  & *q_south(ix ,kk)*half + &
469  & ax_traj *(ax_traj-two) &
470  & *q_south(ixp1,kk)*(-half) + &
471  & (ax_traj+one)* (ax_traj-two) &
472  & *q_south(ixp1,kk)*(-half) + &
473  & (ax_traj+one)* ax_traj &
474  & *q_south(ixp1,kk)*(-half) + &
475  & ax_traj *(ax_traj-one) &
476  & *q_south(ixp2,kk)*(sixth) + &
477  & (ax_traj+one)* (ax_traj-one) &
478  & *q_south(ixp2,kk)*(sixth) + &
479  & (ax_traj+one)* ax_traj &
480  & *q_south(ixp2,kk)*(sixth) &
481  &)*qj
482 
483  qj = 0.0_kind_real
484  else if (jy > ny) then
485  ax = ax + ( &
486  & (ax_traj-one)*(ax_traj-two) &
487  & *q_north(ixm1,kk)*(-sixth) + &
488  & ax_traj * (ax_traj-two) &
489  & *q_north(ixm1,kk)*(-sixth) + &
490  & ax_traj *(ax_traj-one) &
491  & *q_north(ixm1,kk)*(-sixth) + &
492  & (ax_traj-one)*(ax_traj-two) &
493  & *q_north(ix ,kk)*half + &
494  & (ax_traj+one)* (ax_traj-two) &
495  & *q_north(ix ,kk)*half + &
496  & (ax_traj+one)*(ax_traj-one) &
497  & *q_north(ix ,kk)*half + &
498  & ax_traj *(ax_traj-two) &
499  & *q_north(ixp1,kk)*(-half) + &
500  & (ax_traj+one)* (ax_traj-two) &
501  & *q_north(ixp1,kk)*(-half) + &
502  & (ax_traj+one)* ax_traj &
503  & *q_north(ixp1,kk)*(-half) + &
504  & ax_traj *(ax_traj-one) &
505  & *q_north(ixp2,kk)*(sixth) + &
506  & (ax_traj+one)* (ax_traj-one) &
507  & *q_north(ixp2,kk)*(sixth) + &
508  & (ax_traj+one)* ax_traj &
509  & *q_north(ixp2,kk)*(sixth) &
510  &)*qj
511 
512  qj = 0.0_kind_real
513  else
514  ax = ax + ( &
515  & (ax_traj-one)*(ax_traj-two) &
516  & *q_traj(ixm1,jy,kk)*(-sixth) + &
517  & ax_traj * (ax_traj-two) &
518  & *q_traj(ixm1,jy,kk)*(-sixth) + &
519  & ax_traj *(ax_traj-one) &
520  & *q_traj(ixm1,jy,kk)*(-sixth) + &
521  & (ax_traj-one)*(ax_traj-two) &
522  & *q_traj(ix ,jy,kk)*half + &
523  & (ax_traj+one)* (ax_traj-two) &
524  & *q_traj(ix ,jy,kk)*half + &
525  & (ax_traj+one)*(ax_traj-one) &
526  & *q_traj(ix ,jy,kk)*half + &
527  & ax_traj *(ax_traj-two) &
528  & *q_traj(ixp1,jy,kk)*(-half) + &
529  & (ax_traj+one)* (ax_traj-two) &
530  & *q_traj(ixp1,jy,kk)*(-half) + &
531  & (ax_traj+one)* ax_traj &
532  & *q_traj(ixp1,jy,kk)*(-half) + &
533  & ax_traj *(ax_traj-one) &
534  & *q_traj(ixp2,jy,kk)*(sixth) + &
535  & (ax_traj+one)* (ax_traj-one) &
536  & *q_traj(ixp2,jy,kk)*(sixth) + &
537  & (ax_traj+one)* ax_traj &
538  & *q_traj(ixp2,jy,kk)*(sixth) &
539  &)*qj
540 
541  q(ixm1,jy,kk) = q(ixm1,jy,kk) + ( &
542  & ax_traj *(ax_traj-one)*(ax_traj-two)*(-sixth) &
543  &)*qj
544 
545  q(ix ,jy,kk) = q(ix ,jy,kk) + ( &
546  & (ax_traj+one)*(ax_traj-one)*(ax_traj-two)*half &
547  &)*qj
548 
549  q(ixp1,jy,kk) = q(ixp1,jy,kk) + ( &
550  & (ax_traj+one)* ax_traj *(ax_traj-two)*(-half) &
551  &)*qj
552 
553  q(ixp2,jy,kk) = q(ixp2,jy,kk) + ( &
554  & (ax_traj+one)* ax_traj *(ax_traj-one)*(sixth) &
555  &)*qj
556 
557  qj = 0.0_kind_real
558  endif
559 
560  if (jym1 < 1) then
561  ax = ax + ( &
562  & (ax_traj-one)*(ax_traj-two) &
563  & *q_south(ixm1,kk)*(-sixth) + &
564  & ax_traj * (ax_traj-two) &
565  & *q_south(ixm1,kk)*(-sixth) + &
566  & ax_traj *(ax_traj-one) &
567  & *q_south(ixm1,kk)*(-sixth) + &
568  & (ax_traj-one)*(ax_traj-two) &
569  & *q_south(ix, kk)*half + &
570  & (ax_traj+one)* (ax_traj-two) &
571  & *q_south(ix, kk)*half + &
572  & (ax_traj+one)*(ax_traj-one) &
573  & *q_south(ix, kk)*half + &
574  & ax_traj *(ax_traj-two) &
575  & *q_south(ixp1,kk)*(-half) + &
576  & (ax_traj+one)* (ax_traj-two) &
577  & *q_south(ixp1,kk)*(-half) + &
578  & (ax_traj+one)* ax_traj &
579  & *q_south(ixp1,kk)*(-half) + &
580  & ax_traj *(ax_traj-one) &
581  & *q_south(ixp2,kk)*(sixth) + &
582  & (ax_traj+one) * (ax_traj-one) &
583  & *q_south(ixp2,kk)*(sixth) + &
584  & (ax_traj+one) * ax_traj &
585  & *q_south(ixp2,kk)*(sixth) &
586  &)*qjm1
587 
588  qjm1 = 0.0_kind_real
589  else if (jym1 > ny) then
590  ax = ax + ( &
591  & (ax_traj-one)*(ax_traj-two) &
592  & *q_north(ixm1,kk)*(-sixth) + &
593  & ax_traj * (ax_traj-two) &
594  & *q_north(ixm1,kk)*(-sixth) + &
595  & ax_traj *(ax_traj-one) &
596  & *q_north(ixm1,kk)*(-sixth) + &
597  & (ax_traj-one)*(ax_traj-two) &
598  & *q_north(ix, kk)*half + &
599  & (ax_traj+one)* (ax_traj-two) &
600  & *q_north(ix, kk)*half + &
601  & (ax_traj+one)*(ax_traj-one) &
602  & *q_north(ix, kk)*half + &
603  & ax_traj *(ax_traj-two) &
604  & *q_north(ixp1,kk)*(-half) + &
605  & (ax_traj+one)* (ax_traj-two) &
606  & *q_north(ixp1,kk)*(-half) + &
607  & (ax_traj+one)* ax_traj &
608  & *q_north(ixp1,kk)*(-half) + &
609  & ax_traj *(ax_traj-one) &
610  & *q_north(ixp2,kk)*(sixth) + &
611  & (ax_traj+one)* (ax_traj-one) &
612  & *q_north(ixp2,kk)*(sixth) + &
613  & (ax_traj+one)* ax_traj &
614  & *q_north(ixp2,kk)*(sixth) &
615  &)*qjm1
616 
617  qjm1 = 0.0_kind_real
618  else
619  ax = ax + ( &
620  & (ax_traj-one)*(ax_traj-two) &
621  & *q_traj(ixm1,jym1,kk)*(-sixth) + &
622  & ax_traj * (ax_traj-two) &
623  & *q_traj(ixm1,jym1,kk)*(-sixth) + &
624  & ax_traj *(ax_traj-one) &
625  & *q_traj(ixm1,jym1,kk)*(-sixth) + &
626  & (ax_traj-one)*(ax_traj-two) &
627  & *q_traj(ix, jym1,kk)*half + &
628  & (ax_traj+one)* (ax_traj-two) &
629  & *q_traj(ix, jym1,kk)*half + &
630  & (ax_traj+one)*(ax_traj-one) &
631  & *q_traj(ix, jym1,kk)*half + &
632  & ax_traj *(ax_traj-two) &
633  & *q_traj(ixp1,jym1,kk)*(-half) + &
634  & (ax_traj+one)* (ax_traj-two) &
635  & *q_traj(ixp1,jym1,kk)*(-half) + &
636  & (ax_traj+one)* ax_traj &
637  & *q_traj(ixp1,jym1,kk)*(-half) + &
638  & ax_traj *(ax_traj-one) &
639  & *q_traj(ixp2,jym1,kk)*(sixth) + &
640  & (ax_traj+one)* (ax_traj-one) &
641  & *q_traj(ixp2,jym1,kk)*(sixth) + &
642  & (ax_traj+one)* ax_traj &
643  & *q_traj(ixp2,jym1,kk)*(sixth) &
644  &)*qjm1
645 
646  q(ixm1,jym1,kk) = q(ixm1,jym1,kk) + ( &
647  & ax_traj *(ax_traj-one)*(ax_traj-two)*(-sixth) &
648  &)*qjm1
649 
650  q(ix, jym1,kk) = q(ix, jym1,kk) + ( &
651  & (ax_traj+one)*(ax_traj-one)*(ax_traj-two)*half &
652  &)*qjm1
653 
654  q(ixp1,jym1,kk) = q(ixp1,jym1,kk) + ( &
655  & (ax_traj+one)* ax_traj *(ax_traj-two)*(-half) &
656  &)*qjm1
657 
658  q(ixp2,jym1,kk) = q(ixp2,jym1,kk) + ( &
659  & (ax_traj+one)* ax_traj *(ax_traj-one)*(sixth) &
660  &)*qjm1
661 
662  qjm1 = 0.0_kind_real
663  endif
664 
665 !--- find the interpolation point (adjoint)
666 
667  v(ii,jj,kk) = ay * (-dt/deltay)
668  ay = 0.0_kind_real
669 
670  u(ii,jj,kk) = ax * (-dt/deltax)
671  ax = 0.0_kind_real
672 
673  enddo
674  enddo
675 enddo
676 
677 end subroutine advect_pv_ad
advect_pv_ad
subroutine advect_pv_ad(qnew, q, q_traj, q_north, q_south, u, u_traj, v, v_traj, nx, ny, deltax, deltay, dt)
Advect potential vorticity - Adjoint.
Definition: advect_pv_ad.f90:14