OOPS
qg_interp_mod.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 module qg_interp_mod
10 
11 use fckit_log_module,only: fckit_log
12 use kinds
13 use qg_constants_mod
14 use qg_geom_mod
15 use qg_projection_mod
16 use qg_tools_mod
17 
18 implicit none
19 
20 private
21 public :: qg_interp_trilinear,qg_interp_trilinear_ad, &
22  & qg_interp_bicubic,qg_interp_bicubic_tl,qg_interp_bicubic_ad, &
23  & find_x_indices,find_y_indices
24 ! ------------------------------------------------------------------------------
25 contains
26 ! ------------------------------------------------------------------------------
27 !> Trilinear interpolation
28 subroutine qg_interp_trilinear(geom,lon,lat,z,field,val)
29 
30 ! Passed variables
31 type(qg_geom),intent(in) :: geom !< Geometry
32 real(kind_real),intent(in) :: lon !< Longitude
33 real(kind_real),intent(in) :: lat !< Latitude
34 real(kind_real),intent(in) :: z !< Altitude
35 real(kind_real),intent(in) :: field(geom%nx,geom%ny,geom%nz) !< Field
36 real(kind_real),intent(out) :: val !< Value
37 
38 ! Local variables
39 integer :: jxm1,jxo,jxp1,jxp2
40 integer :: jym1,jyo,jyp1,jyp2
41 integer :: jzo,jzp1
42 real(kind_real) :: x,y
43 real(kind_real) :: ax,ay,az
44 real(kind_real) :: oo,op1,p1o,p1p1,o,p1
45 
46 ! Convert lon/lat to x/y
47 call lonlat_to_xy(lon,lat,x,y)
48 
49 ! Find indices
50 call find_x_indices(geom,x,jxm1,jxo,jxp1,jxp2,ax)
51 call find_y_indices(geom,y,jym1,jyo,jyp1,jyp2,ay)
52 call find_z_indices(geom,z,jzo,jzp1,az)
53 
54 ! Extrapolate along y if needed
55 if (jyo<1) then
56  ay = ay+real(jyo-1,kind_real)
57  jyo = 1
58  jyp1 = 2
59 endif
60 if (jyp1>geom%ny) then
61  ay = ay+real(jyp1-geom%ny,kind_real)
62  jyo = geom%ny-1
63  jyp1 = geom%ny
64 endif
65 
66 ! Interpolate along x
67 oo = (1.0-ax)*field(jxo,jyo,jzo)+ax*field(jxp1,jyo,jzo)
68 op1 = (1.0-ax)*field(jxo,jyo,jzp1)+ax*field(jxp1,jyo,jzp1)
69 p1o = (1.0-ax)*field(jxo,jyp1,jzo)+ax*field(jxp1,jyp1,jzo)
70 p1p1 = (1.0-ax)*field(jxo,jyp1,jzp1)+ax*field(jxp1,jyp1,jzp1)
71 
72 ! Interpolate along y
73 o = (1.0-ay)*oo+ay*p1o
74 p1 = (1.0-ay)*op1+ay*p1p1
75 
76 ! Interpolate along z
77 val = (1.0-az)*o+az*p1
78 
79 end subroutine qg_interp_trilinear
80 ! ------------------------------------------------------------------------------
81 !> Interpolation - adjoint
82 subroutine qg_interp_trilinear_ad(geom,lon,lat,z,val,field)
83 
84 ! Passed variables
85 type(qg_geom),intent(in) :: geom !< Geometry
86 real(kind_real),intent(in) :: lon !< Longitude
87 real(kind_real),intent(in) :: lat !< Latitude
88 real(kind_real),intent(in) :: z !< Altitude
89 real(kind_real),intent(in) :: val !< Value
90 real(kind_real),intent(inout) :: field(geom%nx,geom%ny,geom%nz) !< Field
91 
92 ! Local variables
93 integer :: jxm1,jxo,jxp1,jxp2
94 integer :: jym1,jyo,jyp1,jyp2
95 integer :: jzo,jzp1
96 real(kind_real) :: x,y
97 real(kind_real) :: ax,ay,az
98 real(kind_real) :: oo,op1,p1o,p1p1,o,p1
99 
100 ! Convert lon/lat to x/y
101 call lonlat_to_xy(lon,lat,x,y)
102 
103 ! Find indices
104 call find_x_indices(geom,x,jxm1,jxo,jxp1,jxp2,ax)
105 call find_y_indices(geom,y,jym1,jyo,jyp1,jyp2,ay)
106 call find_z_indices(geom,z,jzo,jzp1,az)
107 
108 ! Extrapolate along y if needed
109 if (jyo<1) then
110  ay = ay+real(jyo-1,kind_real)
111  jyo = 1
112  jyp1 = 2
113 endif
114 if (jyp1>geom%ny) then
115  ay = ay+real(jyp1-geom%ny,kind_real)
116  jyo = geom%ny-1
117  jyp1 = geom%ny
118 endif
119 
120 ! Interpolate along z
121 o = (1.0-az)*val
122 p1 = az*val
123 
124 ! Interpolate along y
125 oo = (1.0-ay)*o
126 p1o = ay*o
127 op1 = (1.0-ay)*p1
128 p1p1 = ay*p1
129 
130 ! Interpolate along x
131 field(jxo,jyo,jzo) = field(jxo,jyo,jzo)+(1.0-ax)*oo
132 field(jxp1,jyo,jzo) = field(jxp1,jyo,jzo)+ax*oo
133 field(jxo,jyo,jzp1) = field(jxo,jyo,jzp1)+(1.0-ax)*op1
134 field(jxp1,jyo,jzp1) = field(jxp1,jyo,jzp1)+ax*op1
135 field(jxo,jyp1,jzo) = field(jxo,jyp1,jzo)+(1.0-ax)*p1o
136 field(jxp1,jyp1,jzo) = field(jxp1,jyp1,jzo)+ax*p1o
137 field(jxo,jyp1,jzp1) = field(jxo,jyp1,jzp1)+(1.0-ax)*p1p1
138 field(jxp1,jyp1,jzp1) = field(jxp1,jyp1,jzp1)+ax*p1p1
139 
140 end subroutine qg_interp_trilinear_ad
141 ! ------------------------------------------------------------------------------
142 !> Bicubic horizontal interpolation
143 subroutine qg_interp_bicubic(geom,x,y,gfld2dext,val)
144 
145 ! Passed variables
146 type(qg_geom),intent(in) :: geom !< Geometry
147 real(kind_real),intent(in) :: x !< X value
148 real(kind_real),intent(in) :: y !< Y value
149 real(kind_real),intent(in) :: gfld2dext(geom%nx,-1:geom%ny+2) !< Extended 2D field
150 real(kind_real),intent(out) :: val !< Value
151 
152 ! Local variables
153 integer :: jxm1,jxo,jxp1,jxp2
154 integer :: jym1,jyo,jyp1,jyp2
155 real(kind_real) :: ax,ay,m1,o,p1,p2
156 
157 ! Find indices
158 call find_x_indices(geom,x,jxm1,jxo,jxp1,jxp2,ax)
159 call find_y_indices(geom,y,jym1,jyo,jyp1,jyp2,ay)
160 
161 ! Interpolation along x
162 call cubic(ax,gfld2dext(jxm1,jym1),gfld2dext(jxo,jym1),gfld2dext(jxp1,jym1),gfld2dext(jxp2,jym1),m1)
163 call cubic(ax,gfld2dext(jxm1,jyo),gfld2dext(jxo,jyo),gfld2dext(jxp1,jyo),gfld2dext(jxp2,jyo),o)
164 call cubic(ax,gfld2dext(jxm1,jyp1),gfld2dext(jxo,jyp1),gfld2dext(jxp1,jyp1),gfld2dext(jxp2,jyp1),p1)
165 call cubic(ax,gfld2dext(jxm1,jyp2),gfld2dext(jxo,jyp2),gfld2dext(jxp1,jyp2),gfld2dext(jxp2,jyp2),p2)
166 
167 ! Interpolation along y
168 call cubic(ay,m1,o,p1,p2,val)
169 
170 end subroutine qg_interp_bicubic
171 ! ------------------------------------------------------------------------------
172 !> Bicubic horizontal interpolation - tangent linear
173 subroutine qg_interp_bicubic_tl(geom,x,y,gfld2dext_traj,dx,dy,gfld2dext,val)
174 
175 ! Passed variables
176 type(qg_geom),intent(in) :: geom !< Geometry
177 real(kind_real),intent(in) :: x !< X value
178 real(kind_real),intent(in) :: y !< Y value
179 real(kind_real),intent(in) :: gfld2dext_traj(geom%nx,-1:geom%ny+2) !< Extended 2D trajectory
180 real(kind_real),intent(in) :: dx !< X perturbation
181 real(kind_real),intent(in) :: dy !< Y perturbation
182 real(kind_real),intent(in) :: gfld2dext(geom%nx,-1:geom%ny+2) !< Extended 2D perturbation
183 real(kind_real),intent(out) :: val !< Value
184 
185 ! Local variables
186 integer :: jxm1,jxo,jxp1,jxp2
187 integer :: jym1,jyo,jyp1,jyp2
188 real(kind_real) :: ax_traj,ay_traj,m1_traj,o_traj,p1_traj,p2_traj,ax,ay,m1,o,p1,p2
189 
190 ! Find indices
191 call find_x_indices(geom,x,jxm1,jxo,jxp1,jxp2,ax_traj)
192 call find_y_indices(geom,y,jym1,jyo,jyp1,jyp2,ay_traj)
193 
194 ! Interpolation along x (trajectory)
195 call cubic(ax_traj,gfld2dext_traj(jxm1,jym1),gfld2dext_traj(jxo,jym1),gfld2dext_traj(jxp1,jym1),gfld2dext_traj(jxp2,jym1),m1_traj)
196 call cubic(ax_traj,gfld2dext_traj(jxm1,jyo),gfld2dext_traj(jxo,jyo),gfld2dext_traj(jxp1,jyo),gfld2dext_traj(jxp2,jyo),o_traj)
197 call cubic(ax_traj,gfld2dext_traj(jxm1,jyp1),gfld2dext_traj(jxo,jyp1),gfld2dext_traj(jxp1,jyp1),gfld2dext_traj(jxp2,jyp1),p1_traj)
198 call cubic(ax_traj,gfld2dext_traj(jxm1,jyp2),gfld2dext_traj(jxo,jyp2),gfld2dext_traj(jxp1,jyp2),gfld2dext_traj(jxp2,jyp2),p2_traj)
199 
200 ! Compute linearization coefficients
201 ax = dx/geom%deltax
202 ay = dy/geom%deltay
203 
204 ! Interpolation along x (perturbation)
205 call cubic_tl(ax_traj,gfld2dext_traj(jxm1,jym1),gfld2dext_traj(jxo,jym1),gfld2dext_traj(jxp1,jym1),gfld2dext_traj(jxp2,jym1), &
206  & ax,gfld2dext(jxm1,jym1),gfld2dext(jxo,jym1),gfld2dext(jxp1,jym1),gfld2dext(jxp2,jym1),m1)
207 call cubic_tl(ax_traj,gfld2dext_traj(jxm1,jyo),gfld2dext_traj(jxo,jyo),gfld2dext_traj(jxp1,jyo),gfld2dext_traj(jxp2,jyo), &
208  & ax,gfld2dext(jxm1,jyo),gfld2dext(jxo,jyo),gfld2dext(jxp1,jyo),gfld2dext(jxp2,jyo),o)
209 call cubic_tl(ax_traj,gfld2dext_traj(jxm1,jyp1),gfld2dext_traj(jxo,jyp1),gfld2dext_traj(jxp1,jyp1),gfld2dext_traj(jxp2,jyp1), &
210  & ax,gfld2dext(jxm1,jyp1),gfld2dext(jxo,jyp1),gfld2dext(jxp1,jyp1),gfld2dext(jxp2,jyp1),p1)
211 call cubic_tl(ax_traj,gfld2dext_traj(jxm1,jyp2),gfld2dext_traj(jxo,jyp2),gfld2dext_traj(jxp1,jyp2),gfld2dext_traj(jxp2,jyp2), &
212  & ax,gfld2dext(jxm1,jyp2),gfld2dext(jxo,jyp2),gfld2dext(jxp1,jyp2),gfld2dext(jxp2,jyp2),p2)
213 
214 ! Interpolation along y (perturbation)
215 call cubic_tl(ay_traj,m1_traj,o_traj,p1_traj,p2_traj,ay,m1,o,p1,p2,val)
216 
217 end subroutine qg_interp_bicubic_tl
218 ! ------------------------------------------------------------------------------
219 !> Bicubic horizontal interpolation - adjoint
220 subroutine qg_interp_bicubic_ad(geom,x,y,gfld2dext_traj,val,dx,dy,gfld2dext)
221 
222 ! Passed variables
223 type(qg_geom),intent(in) :: geom !< Geometry
224 real(kind_real),intent(in) :: x !< X value
225 real(kind_real),intent(in) :: y !< Y value
226 real(kind_real),intent(in) :: gfld2dext_traj(geom%nx,-1:geom%ny+2) !< Extended 2D trajectory
227 real(kind_real),intent(in) :: val !< Value
228 real(kind_real),intent(inout) :: dx !< X perturbation
229 real(kind_real),intent(inout) :: dy !< Y perturbation
230 real(kind_real),intent(inout) :: gfld2dext(geom%nx,-1:geom%ny+2) !< Extended 2D perturbation
231 
232 ! Local variables
233 integer :: jxm1,jxo,jxp1,jxp2
234 integer :: jym1,jyo,jyp1,jyp2
235 real(kind_real) :: ax_traj,ay_traj,m1_traj,o_traj,p1_traj,p2_traj,ax,ay,m1,o,p1,p2
236 
237 ! Find indices
238 call find_x_indices(geom,x,jxm1,jxo,jxp1,jxp2,ax_traj)
239 call find_y_indices(geom,y,jym1,jyo,jyp1,jyp2,ay_traj)
240 
241 ! Interpolation along x (trajectory)
242 call cubic(ax_traj,gfld2dext_traj(jxm1,jym1),gfld2dext_traj(jxo,jym1),gfld2dext_traj(jxp1,jym1),gfld2dext_traj(jxp2,jym1),m1_traj)
243 call cubic(ax_traj,gfld2dext_traj(jxm1,jyo),gfld2dext_traj(jxo,jyo),gfld2dext_traj(jxp1,jyo),gfld2dext_traj(jxp2,jyo),o_traj)
244 call cubic(ax_traj,gfld2dext_traj(jxm1,jyp1),gfld2dext_traj(jxo,jyp1),gfld2dext_traj(jxp1,jyp1),gfld2dext_traj(jxp2,jyp1),p1_traj)
245 call cubic(ax_traj,gfld2dext_traj(jxm1,jyp2),gfld2dext_traj(jxo,jyp2),gfld2dext_traj(jxp1,jyp2),gfld2dext_traj(jxp2,jyp2),p2_traj)
246 
247 ! Initialization
248 m1 = 0.0
249 o = 0.0
250 p1 = 0.0
251 p2 = 0.0
252 ax = 0.0
253 ay = 0.0
254 
255 ! Interpolation along y (perturbation)
256 call cubic_ad(ay_traj,m1_traj,o_traj,p1_traj,p2_traj,val,ay,m1,o,p1,p2)
257 
258 ! Interpolation along x (perturbation)
259 call cubic_ad(ax_traj,gfld2dext_traj(jxm1,jym1),gfld2dext_traj(jxo,jym1),gfld2dext_traj(jxp1,jym1),gfld2dext_traj(jxp2,jym1), &
260  & m1,ax,gfld2dext(jxm1,jym1),gfld2dext(jxo,jym1),gfld2dext(jxp1,jym1),gfld2dext(jxp2,jym1))
261 call cubic_ad(ax_traj,gfld2dext_traj(jxm1,jyo),gfld2dext_traj(jxo,jyo),gfld2dext_traj(jxp1,jyo),gfld2dext_traj(jxp2,jyo), &
262  & o,ax,gfld2dext(jxm1,jyo),gfld2dext(jxo,jyo),gfld2dext(jxp1,jyo),gfld2dext(jxp2,jyo))
263 call cubic_ad(ax_traj,gfld2dext_traj(jxm1,jyp1),gfld2dext_traj(jxo,jyp1),gfld2dext_traj(jxp1,jyp1),gfld2dext_traj(jxp2,jyp1), &
264  & p1,ax,gfld2dext(jxm1,jyp1),gfld2dext(jxo,jyp1),gfld2dext(jxp1,jyp1),gfld2dext(jxp2,jyp1))
265 call cubic_ad(ax_traj,gfld2dext_traj(jxm1,jyp2),gfld2dext_traj(jxo,jyp2),gfld2dext_traj(jxp1,jyp2),gfld2dext_traj(jxp2,jyp2), &
266  & p2,ax,gfld2dext(jxm1,jyp2),gfld2dext(jxo,jyp2),gfld2dext(jxp1,jyp2),gfld2dext(jxp2,jyp2))
267 
268 ! Compute linearization coefficients
269 dx = dx+ax/geom%deltax
270 dy = dy+ay/geom%deltay
271 
272 end subroutine qg_interp_bicubic_ad
273 ! ------------------------------------------------------------------------------
274 !> Find x indices
275 subroutine find_x_indices(geom,x,jxm1,jxo,jxp1,jxp2,ax)
276 
277 ! Passed variables
278 type(qg_geom),intent(in) :: geom !< Geometry
279 real(kind_real),intent(in) :: x !< X value
280 integer,intent(out) :: jxm1 !< Minus 1 index
281 integer,intent(out) :: jxo !< Origin index
282 integer,intent(out) :: jxp1 !< Plus 1 index
283 integer,intent(out) :: jxp2 !< Plus 2 index
284 real(kind_real),intent(out) :: ax !< Coefficient
285 
286 ! Local variables
287 integer :: ix
288 real(kind_real) :: xx
289 
290 ! Initialization
291 if (x<0.0) then
292  xx = x+domain_zonal
293 elseif (x>domain_zonal) then
294  xx = x-domain_zonal
295 else
296  xx = x
297 endif
298 
299 ! Find base index and weight
300 if ((0.0<=xx).and.(xx<geom%x(1))) then
301  jxo = geom%nx
302  ax = (domain_zonal-geom%x(jxo)+xx)/geom%deltax
303 elseif ((geom%x(geom%nx)<=xx).and.(xx<=domain_zonal)) then
304  jxo = geom%nx
305  ax = (xx-geom%x(jxo))/geom%deltax
306 else
307  jxo = -1
308  do ix=1,geom%nx-1
309  if ((geom%x(ix)<=xx).and.(xx<geom%x(ix+1))) then
310  jxo = ix
311  exit
312  end if
313  end do
314  if (jxo==-1) call abor1_ftn('find_x_indices: cannot find jxo')
315  ax = (xx-geom%x(jxo))/geom%deltax
316 end if
317 
318 ! Define other indices
319 jxm1 = jxo-1
320 if (jxm1<1) jxm1 = jxm1+geom%nx
321 jxp1 = jxo+1
322 if (jxp1>geom%nx) jxp1 = jxp1-geom%nx
323 jxp2 = jxo+2
324 if (jxp2>geom%nx) jxp2 = jxp2-geom%nx
325 
326 end subroutine find_x_indices
327 ! ------------------------------------------------------------------------------
328 !> Find y indices
329 subroutine find_y_indices(geom,y,jym1,jyo,jyp1,jyp2,ay)
330 
331 ! Passed variables
332 type(qg_geom),intent(in) :: geom !< Geometry
333 real(kind_real),intent(in) :: y !< Y value
334 integer,intent(out) :: jym1 !< Minus 1 index
335 integer,intent(out) :: jyo !< Origin index
336 integer,intent(out) :: jyp1 !< Plus 1 index
337 integer,intent(out) :: jyp2 !< Plus 2 index
338 real(kind_real),intent(out) :: ay !< Coefficient
339 
340 ! Local variables
341 integer :: iy
342 real(kind_real) :: yy
343 
344 ! Initialization
345 yy = max(0.0_kind_real,min(y,domain_meridional))
346 
347 ! Find base index and weight
348 if ((0.0<=yy).and.(yy<geom%y(1))) then
349  jyo = 0
350  ay = yy/geom%deltay
351 elseif ((geom%y(geom%ny)<=yy).and.(yy<=domain_meridional)) then
352  jyo = geom%ny
353  ay = (yy-geom%y(jyo))/geom%deltay
354 else
355  jyo = -1
356  do iy=1,geom%ny-1
357  if ((geom%y(iy)<=yy).and.(yy<geom%y(iy+1))) then
358  jyo = iy
359  exit
360  end if
361  end do
362  if (jyo==-1) call abor1_ftn('find_y_indices: cannot find jyo')
363  ay = (yy-geom%y(jyo))/geom%deltay
364 end if
365 
366 ! Define other indices
367 jym1 = jyo-1
368 jyp1 = jyo+1
369 jyp2 = jyo+2
370 
371 end subroutine find_y_indices
372 ! ------------------------------------------------------------------------------
373 !> Find z indices
374 subroutine find_z_indices(geom,z,jzo,jzp1,az)
375 
376 ! Passed variables
377 type(qg_geom),intent(in) :: geom !< Geometry
378 real(kind_real),intent(in) :: z !< Z value
379 integer,intent(out) :: jzo !< Origin index
380 integer,intent(out) :: jzp1 !< Plus 1 index
381 real(kind_real),intent(out) :: az !< Coefficient
382 
383 ! Local variables
384 integer :: iz
385 real(kind_real) :: zz
386 
387 ! Initialization
388 zz = max(0.0_kind_real,min(z,domain_depth))
389 
390 ! Find base index and weight
391 if ((0.0<=zz).and.(zz<geom%z(1))) then
392  jzo = 1
393 elseif ((geom%z(geom%nz)<=zz).and.(zz<=domain_depth)) then
394  jzo = geom%nz-1
395 else
396  jzo = -1
397  do iz=1,geom%nz-1
398  if ((geom%z(iz)<=zz).and.(zz<geom%z(iz+1))) then
399  jzo = iz
400  exit
401  end if
402  end do
403  if (jzo==-1) call abor1_ftn('find_z_indices: cannot find jzo')
404 end if
405 az = (zz-geom%z(jzo))/(geom%z(jzo+1)-geom%z(jzo))
406 
407 ! Define other indices
408 jzp1 = jzo+1
409 
410 end subroutine find_z_indices
411 ! ------------------------------------------------------------------------------
412 !> Cubic interpolation
413 subroutine cubic(a,m1,o,p1,p2,res)
414 
415 ! Passed variables
416 real(kind_real),intent(in) :: a !< Coefficient
417 real(kind_real),intent(in) :: m1 !< Minus 1 value
418 real(kind_real),intent(in) :: o !< Origin value
419 real(kind_real),intent(in) :: p1 !< Plus 1 value
420 real(kind_real),intent(in) :: p2 !< Plus 2 value
421 real(kind_real),intent(out) :: res !< Result
422 
423 res = -a*(a-1.0)*(a-2.0)/6.0*m1 &
424  & +(a+1.0)*(a-1.0)*(a-2.0)/2.0*o &
425  & -(a+1.0)*a*(a-2.0)/2.0*p1 &
426  & +(a+1.0)*a*(a-1.0)/6.0*p2
427 
428 end subroutine cubic
429 ! ------------------------------------------------------------------------------
430 !> Cubic interpolation - tangent linear
431 subroutine cubic_tl(a_traj,m1_traj,o_traj,p1_traj,p2_traj,a,m1,o,p1,p2,res)
432 
433 ! Passed variables
434 real(kind_real),intent(in) :: a_traj !< Coefficient trajectory
435 real(kind_real),intent(in) :: m1_traj !< Minus 1 trajectory
436 real(kind_real),intent(in) :: o_traj !< Origin value
437 real(kind_real),intent(in) :: p1_traj !< Plus 1 value
438 real(kind_real),intent(in) :: p2_traj !< Plus 2 value
439 real(kind_real),intent(in) :: a !< Coefficient perturbation
440 real(kind_real),intent(in) :: m1 !< Minus perturbation
441 real(kind_real),intent(in) :: o !< Origin perturbation
442 real(kind_real),intent(in) :: p1 !< Plus 1 perturbation
443 real(kind_real),intent(in) :: p2 !< Plus 2 perturbation
444 real(kind_real),intent(out) :: res !< Result
445 
446 res = -a*(a_traj-1.0)*(a_traj-2.0)/6.0*m1_traj &
447  & -a_traj*a*(a_traj-2.0)/6.0*m1_traj &
448  & -a_traj*(a_traj-1.0)*a/6.0*m1_traj &
449  & -a_traj*(a_traj-1.0)*(a_traj-2.0)/6.0*m1 &
450  & +a*(a_traj-1.0)*(a_traj-2.0)/2.0*o_traj &
451  & +(a_traj+1.0)*a*(a_traj-2.0)/2.0*o_traj &
452  & +(a_traj+1.0)*(a_traj-1.0)*a/2.0*o_traj &
453  & +(a_traj+1.0)*(a_traj-1.0)*(a_traj-2.0)/2.0*o &
454  & -a*a_traj*(a_traj-2.0)/2.0*p1_traj &
455  & -(a_traj+1.0)*a*(a_traj-2.0)/2.0*p1_traj &
456  & -(a_traj+1.0)*a_traj*a/2.0*p1_traj &
457  & -(a_traj+1.0)*a_traj*(a_traj-2.0)/2.0*p1 &
458  & +a*a_traj*(a_traj-1.0)/6.0*p2_traj &
459  & +(a_traj+1.0)*a*(a_traj-1.0)/6.0*p2_traj &
460  & +(a_traj+1.0)*a_traj*a/6.0*p2_traj &
461  & +(a_traj+1.0)*a_traj*(a_traj-1.0)/6.0*p2
462 
463 end subroutine cubic_tl
464 ! ------------------------------------------------------------------------------
465 !> Cubic interpolation - adjoint
466 subroutine cubic_ad(a_traj,m1_traj,o_traj,p1_traj,p2_traj,res,a,m1,o,p1,p2)
467 
468 ! Passed variables
469 real(kind_real),intent(in) :: a_traj !< Coefficient trajectory
470 real(kind_real),intent(in) :: m1_traj !< Minus 1 trajectory
471 real(kind_real),intent(in) :: o_traj !< Origin value
472 real(kind_real),intent(in) :: p1_traj !< Plus 1 value
473 real(kind_real),intent(in) :: p2_traj !< Plus 2 value
474 real(kind_real),intent(in) :: res !< Result
475 real(kind_real),intent(inout) :: a !< Coefficient perturbation
476 real(kind_real),intent(inout) :: m1 !< Minus perturbation
477 real(kind_real),intent(inout) :: o !< Origin perturbation
478 real(kind_real),intent(inout) :: p1 !< Plus 1 perturbation
479 real(kind_real),intent(inout) :: p2 !< Plus 2 perturbation
480 
481 a = a+(-(a_traj-1.0)*(a_traj-2.0)/6.0*m1_traj &
482  & -a_traj*(a_traj-2.0)/6.0*m1_traj &
483  & -a_traj*(a_traj-1.0)/6.0*m1_traj &
484  & +(a_traj-1.0)*(a_traj-2.0)/2.0*o_traj &
485  & +(a_traj+1.0)*(a_traj-2.0)/2.0*o_traj &
486  & +(a_traj+1.0)*(a_traj-1.0)/2.0*o_traj &
487  & -a_traj*(a_traj-2.0)/2.0*p1_traj &
488  & -(a_traj+1.0)*(a_traj-2.0)/2.0*p1_traj &
489  & -(a_traj+1.0)*a_traj/2.0*p1_traj &
490  & +a_traj*(a_traj-1.0)/6.0*p2_traj &
491  & +(a_traj+1.0)*(a_traj-1.0)/6.0*p2_traj &
492  & +(a_traj+1.0)*a_traj/6.0*p2_traj) &
493  & *res
494 m1 = m1-a_traj*(a_traj-1.0)*(a_traj-2.0)/6.0*res
495 o = o+(a_traj+1.0)*(a_traj-1.0)*(a_traj-2.0)/2.0*res
496 p1 = p1-(a_traj+1.0)*a_traj*(a_traj-2.0)/2.0*res
497 p2 = p2+(a_traj+1.0)*a_traj*(a_traj-1.0)/6.0*res
498 
499 end subroutine cubic_ad
500 ! ------------------------------------------------------------------------------
501 end module qg_interp_mod
qg_constants_mod
Definition: qg_constants_mod.F90:9
qg_constants_mod::domain_depth
real(kind_real), parameter, public domain_depth
Model depth (m)
Definition: qg_constants_mod.F90:34
qg_constants_mod::domain_meridional
real(kind_real), parameter, public domain_meridional
Model domain in meridional direction (m)
Definition: qg_constants_mod.F90:33
qg_constants_mod::domain_zonal
real(kind_real), parameter, public domain_zonal
Model domain in zonal direction (m)
Definition: qg_constants_mod.F90:32
qg_geom_mod
Definition: qg_geom_mod.F90:9
qg_interp_mod
Definition: qg_interp_mod.F90:9
qg_interp_mod::find_x_indices
subroutine, public find_x_indices(geom, x, jxm1, jxo, jxp1, jxp2, ax)
Find x indices.
Definition: qg_interp_mod.F90:276
qg_interp_mod::qg_interp_bicubic_ad
subroutine, public qg_interp_bicubic_ad(geom, x, y, gfld2dext_traj, val, dx, dy, gfld2dext)
Bicubic horizontal interpolation - adjoint.
Definition: qg_interp_mod.F90:221
qg_interp_mod::find_z_indices
subroutine find_z_indices(geom, z, jzo, jzp1, az)
Find z indices.
Definition: qg_interp_mod.F90:375
qg_interp_mod::cubic
subroutine cubic(a, m1, o, p1, p2, res)
Cubic interpolation.
Definition: qg_interp_mod.F90:414
qg_interp_mod::qg_interp_bicubic_tl
subroutine, public qg_interp_bicubic_tl(geom, x, y, gfld2dext_traj, dx, dy, gfld2dext, val)
Bicubic horizontal interpolation - tangent linear.
Definition: qg_interp_mod.F90:174
qg_interp_mod::cubic_ad
subroutine cubic_ad(a_traj, m1_traj, o_traj, p1_traj, p2_traj, res, a, m1, o, p1, p2)
Cubic interpolation - adjoint.
Definition: qg_interp_mod.F90:467
qg_interp_mod::qg_interp_trilinear_ad
subroutine, public qg_interp_trilinear_ad(geom, lon, lat, z, val, field)
Interpolation - adjoint.
Definition: qg_interp_mod.F90:83
qg_interp_mod::find_y_indices
subroutine, public find_y_indices(geom, y, jym1, jyo, jyp1, jyp2, ay)
Find y indices.
Definition: qg_interp_mod.F90:330
qg_interp_mod::qg_interp_trilinear
subroutine, public qg_interp_trilinear(geom, lon, lat, z, field, val)
Trilinear interpolation.
Definition: qg_interp_mod.F90:29
qg_interp_mod::qg_interp_bicubic
subroutine, public qg_interp_bicubic(geom, x, y, gfld2dext, val)
Bicubic horizontal interpolation.
Definition: qg_interp_mod.F90:144
qg_interp_mod::cubic_tl
subroutine cubic_tl(a_traj, m1_traj, o_traj, p1_traj, p2_traj, a, m1, o, p1, p2, res)
Cubic interpolation - tangent linear.
Definition: qg_interp_mod.F90:432
qg_projection_mod
Definition: qg_projection_mod.F90:9
qg_projection_mod::lonlat_to_xy
subroutine, public lonlat_to_xy(lon, lat, x, y)
Convert lon/lat to x/y.
Definition: qg_projection_mod.F90:57
qg_tools_mod
Definition: qg_tools_mod.F90:9
qg_geom_mod::qg_geom
Definition: qg_geom_mod.F90:26