FV3 Bundle
CRTM_Interpolation.f90
Go to the documentation of this file.
1 !
2 ! CRTM_Interpolation
3 !
4 ! Module containing the generic polynomial interpolation
5 ! routines used in the CRTM
6 !
7 !
8 ! CREATION HISTORY:
9 ! Written by: Paul van Delst, 01-Feb-2007
10 ! paul.vandelst@noaa.gov
11 !
12 MODULE crtm_interpolation
13 
14  ! -----------------
15  ! Environment setup
16  ! -----------------
17  ! Module usage
18  USE type_kinds, ONLY: fp
19  ! Disable implicit typing
20  IMPLICIT NONE
21 
22 
23  ! ------------
24  ! Visibilities
25  ! ------------
26  PRIVATE
27  ! Parameters
28  PUBLIC :: order
29  PUBLIC :: npts
30  ! Derived types and associated procedures
31  PUBLIC :: lpoly_type
32  PUBLIC :: clear_lpoly
33  PUBLIC :: lpoly_init
34  PUBLIC :: lpoly_inspect
35  ! Procedures
36  PUBLIC :: interp_1d
37  PUBLIC :: interp_2d
38  PUBLIC :: interp_3d
39  PUBLIC :: interp_4d
40  PUBLIC :: interp_1d_tl
41  PUBLIC :: interp_2d_tl
42  PUBLIC :: interp_3d_tl
43  PUBLIC :: interp_4d_tl
44  PUBLIC :: interp_1d_ad
45  PUBLIC :: interp_2d_ad
46  PUBLIC :: interp_3d_ad
47  PUBLIC :: interp_4d_ad
48  PUBLIC :: find_index
49  PUBLIC :: lpoly
50  PUBLIC :: lpoly_tl
51  PUBLIC :: lpoly_ad
52 
53 
54  ! -------------------
55  ! Procedure overloads
56  ! -------------------
57  INTERFACE find_index
58  MODULE PROCEDURE find_regular_index
59  MODULE PROCEDURE find_random_index
60  END INTERFACE find_index
61 
62 
63  ! -----------------
64  ! Module parameters
65  ! -----------------
66  CHARACTER(*), PARAMETER :: module_rcs_id=&
67  '$Id: CRTM_Interpolation.f90 60152 2015-08-13 19:19:13Z paul.vandelst@noaa.gov $'
68  REAL(fp), PARAMETER :: zero = 0.0_fp
69  REAL(fp), PARAMETER :: one = 1.0_fp
70  INTEGER, PARAMETER :: order = 2 ! Quadratic
71  INTEGER, PARAMETER :: npoly_pts = order+1 ! No. of points in each polynomial
72  INTEGER, PARAMETER :: npts = npoly_pts+1 ! No. of points total
73 
74 
75  ! -----------------------
76  ! Derived type definition
77  ! -----------------------
78  TYPE :: lpoly_type
79 ! PRIVATE
80  INTEGER :: order=order
81  INTEGER :: npts =npoly_pts
82  ! Left and right side polynomials
83  REAL(fp) :: lp_left(npoly_pts) = zero
84  REAL(fp) :: lp_right(npoly_pts) = zero
85  ! Left and right side weighting factors
86  REAL(fp) :: w_left = zero
87  REAL(fp) :: w_right = zero
88  ! Polynomial numerator differences
89  REAL(fp) :: dxi_left(npoly_pts) = zero
90  REAL(fp) :: dxi_right(npoly_pts) = zero
91  ! Polynomial denominator differences
92  REAL(fp) :: dx_left(npoly_pts) = zero
93  REAL(fp) :: dx_right(npoly_pts) = zero
94  END TYPE lpoly_type
95 
96 
97 CONTAINS
98 
99 
100 !################################################################################
101 !################################################################################
102 !## ##
103 !## ## PUBLIC MODULE ROUTINES ## ##
104 !## ##
105 !################################################################################
106 !################################################################################
107 
108 !--------------------------------------------------------------------------------
109 !
110 ! NAME:
111 ! Clear_LPoly
112 !
113 ! PURPOSE:
114 ! Subroutine to reinitialise the LPoly_type structure
115 !
116 ! CALLING SEQUENCE:
117 ! CALL Clear_LPoly(p)
118 !
119 ! OUTPUT ARGUMENTS:
120 ! p: Reinitialised Lagrangian polynomial structure
121 ! UNITS: N/A
122 ! TYPE: TYPE(LPoly_type)
123 ! DIMENSION: Scalar
124 ! ATTRIBUTES: INTENT(IN OUT)
125 !
126 ! SIDE EFFECTS:
127 ! If the structure contains data on input, it is replaced with the
128 ! reinitialisation value.
129 !
130 !--------------------------------------------------------------------------------
131  SUBROUTINE clear_lpoly(p)
132  TYPE(lpoly_type), INTENT(IN OUT) :: p
133  p%Order = order
134  p%nPts = npoly_pts
135  p%lp_left = zero
136  p%lp_right = zero
137  p%w_left = zero
138  p%w_right = zero
139  p%dxi_left = zero
140  p%dxi_right = zero
141  p%dx_left = zero
142  p%dx_right = zero
143  END SUBROUTINE clear_lpoly
144 
145 
146  SUBROUTINE lpoly_init(self)
147  TYPE(lpoly_type), INTENT(OUT) :: self
148  self%Order = order
149  END SUBROUTINE lpoly_init
150 
151 
152  SUBROUTINE lpoly_inspect( self )
153  TYPE(lpoly_type), INTENT(IN) :: self
154  WRITE(*,'(1x,"LPoly OBJECT")')
155  WRITE(*,'(3x,"Order : ",i0)') self%Order
156  WRITE(*,'(3x,"nPts : ",i0)') self%nPts
157  WRITE(*,'(3x,"Left-side :")')
158  WRITE(*,'(5x,"Weighting factor : ",es13.6)') self%w_left
159  WRITE(*,'(5x,"Polynomial :")')
160  WRITE(*,'(5(1x,es13.6,:))') self%lp_left
161  WRITE(*,'(5x,"Numerator differences :")')
162  WRITE(*,'(5(1x,es13.6,:))') self%dxi_left
163  WRITE(*,'(5x,"Denominator differences :")')
164  WRITE(*,'(5(1x,es13.6,:))') self%dx_left
165  WRITE(*,'(3x,"Right-side :")')
166  WRITE(*,'(5x,"Weighting factor : ",es13.6)') self%w_right
167  WRITE(*,'(5x,"Polynomial :")')
168  WRITE(*,'(5(1x,es13.6,:))') self%lp_right
169  WRITE(*,'(5x,"Numerator differences :")')
170  WRITE(*,'(5(1x,es13.6,:))') self%dxi_right
171  WRITE(*,'(5x,"Denominator differences :")')
172  WRITE(*,'(5(1x,es13.6,:))') self%dx_right
173  END SUBROUTINE lpoly_inspect
174 
175 
176 !--------------------------------------------------------------------------------
177 !
178 ! NAME:
179 ! Interp_1D
180 ! Interp_2D
181 ! Interp_3D
182 ! Interp_4D
183 !
184 ! PURPOSE:
185 ! Subroutines to perform interpolation:
186 ! o 1-D for z=f(u)
187 ! o 2-D for z=f(u,v)
188 ! o 3-D for z=f(u,v,w)
189 ! o 4-D for z=f(u,v,w,x)
190 !
191 ! CALLING SEQUENCE:
192 ! CALL Interp_1D(z, ulp, z_int)
193 ! CALL Interp_2D(z, ulp, vlp, z_int)
194 ! CALL Interp_3D(z, ulp, vlp, wlp, z_int)
195 ! CALL Interp_4D(z, ulp, vlp, wlp, xlp, z_int)
196 !
197 ! INPUT ARGUMENTS:
198 ! z: Data to interpolate
199 ! UNITS: Variable
200 ! TYPE: REAL(fp)
201 ! DIMENSION: Rank-1 (N_PTS), or
202 ! Rank-2 (N_PTS,NPTS), or
203 ! Rank-3 (N_PTS,NPTS,N_PTS), or
204 ! Rank-4 (N_PTS,NPTS,N_PTS,NPTS), or
205 ! ATTRIBUTES: INTENT(IN)
206 !
207 ! ulp, vlp, wlp, xlp: Interpolating polynomial structures for the
208 ! respective dimensions from previous calls to
209 ! the LPoly() subroutine
210 ! UNITS: N/A
211 ! TYPE: LPoly_type
212 ! DIMENSION: Scalar
213 ! ATTRIBUTES: INTENT(IN)
214 !
215 ! OUTPUT ARGUMENTS:
216 ! z_int: Interpolation result
217 ! UNITS: Same as z
218 ! TYPE: REAL(fp)
219 ! DIMENSION: Scalar
220 ! ATTRIBUTES: INTENT(IN OUT)
221 !
222 ! COMMENTS:
223 ! Output z_int argument has INTENT(IN OUT) to prevent default reinitialisation
224 ! purely for computational speed.
225 !
226 !--------------------------------------------------------------------------------
227  ! 1-D routine
228  SUBROUTINE interp_1d(z, ulp, & ! Input
229  z_int ) ! Output
230  ! Arguments
231  REAL(fp), INTENT(IN) :: z(:)
232  TYPE(lpoly_type), INTENT(IN) :: ulp
233  REAL(fp), INTENT(IN OUT) :: z_int ! INTENT(IN OUT) to preclude reinitialisation
234  ! Perform interpolation
235  z_int = ( ulp%w_left * ( ulp%lp_left(1) *z(1) + &
236  ulp%lp_left(2) *z(2) + &
237  ulp%lp_left(3) *z(3) ) ) + &
238  ( ulp%w_right * ( ulp%lp_right(1)*z(2) + &
239  ulp%lp_right(2)*z(3) + &
240  ulp%lp_right(3)*z(4) ) )
241  END SUBROUTINE interp_1d
242 
243  ! 2-D routine
244  SUBROUTINE interp_2d(z, ulp, vlp, & ! Input
245  z_int ) ! Output
246  ! Arguments
247  REAL(fp), INTENT(IN) :: z(:,:)
248  TYPE(lpoly_type), INTENT(IN) :: ulp, vlp
249  REAL(fp), INTENT(IN OUT) :: z_int ! INTENT(IN OUT) to preclude reinitialisation
250  ! Local variables
251  INTEGER :: i
252  REAL(fp) :: a(npts)
253  ! Interpolate z in u dimension for all v
254  DO i = 1, npts
255  CALL interp_1d(z(:,i),ulp,a(i))
256  END DO
257  ! Interpolate z in w dimension
258  CALL interp_1d(a,vlp,z_int)
259  END SUBROUTINE interp_2d
260 
261  ! 3-D routine
262  SUBROUTINE interp_3d(z, ulp, vlp, wlp, & ! Input
263  z_int ) ! Output
264  ! Arguments
265  REAL(fp), INTENT(IN) :: z(:,:,:)
266  TYPE(lpoly_type), INTENT(IN) :: ulp, vlp, wlp
267  REAL(fp), INTENT(IN OUT) :: z_int ! INTENT(IN OUT) to preclude reinitialisation
268  ! Local variables
269  INTEGER :: i
270  REAL(fp) :: a(npts)
271  ! Interpolate z in u,v dimension for all w
272  DO i = 1, npts
273  CALL interp_2d(z(:,:,i),ulp,vlp,a(i))
274  END DO
275  ! Interpolate a in w dimension
276  CALL interp_1d(a,wlp,z_int)
277  END SUBROUTINE interp_3d
278 
279  ! 4-D routine
280  SUBROUTINE interp_4d(z, ulp, vlp, wlp, xlp, & ! Input
281  z_int ) ! Output
282  ! Arguments
283  REAL(fp), INTENT(IN) :: z(:,:,:,:)
284  TYPE(lpoly_type), INTENT(IN) :: ulp, vlp, wlp, xlp
285  REAL(fp), INTENT(IN OUT) :: z_int ! INTENT(IN OUT) to preclude reinitialisation
286  ! Local variables
287  INTEGER :: i
288  REAL(fp) :: a(npts)
289  ! Interpolate z in u,v,w dimension for all x
290  DO i = 1, npts
291  CALL interp_3d(z(:,:,:,i),ulp,vlp,wlp,a(i))
292  END DO
293  ! Interpolate a in x dimension
294  CALL interp_1d(a,xlp,z_int)
295  END SUBROUTINE interp_4d
296 
297 
298 !--------------------------------------------------------------------------------
299 !
300 ! NAME:
301 ! Interp_1D_TL
302 ! Interp_2D_TL
303 ! Interp_3D_TL
304 ! Interp_4D_TL
305 !
306 ! PURPOSE:
307 ! Subroutines to perform tangent-linear interpolation:
308 ! o 1-D for z=f(u)
309 ! o 2-D for z=f(u,v)
310 ! o 3-D for z=f(u,v,w)
311 ! o 4-D for z=f(u,v,w,x)
312 !
313 ! CALLING SEQUENCE:
314 ! CALL Interp_1D_TL(z, ulp, z_TL, ulp_TL, z_int_TL)
315 ! CALL Interp_2D_TL(z, ulp, vlp, z_TL, ulp_TL, vlp_TL, z_int_TL)
316 ! CALL Interp_3D_TL(z, ulp, vlp, wlp, z_TL, ulp_TL, vlp_TL, wlp_TL, z_int_TL)
317 ! CALL Interp_4D_TL(z, ulp, vlp, wlp, xlp, z_TL, ulp_TL, vlp_TL, wlp_TL, xlp_TL, z_int_TL)
318 !
319 ! INPUT ARGUMENTS:
320 ! z: Data to interpolate
321 ! UNITS: Variable
322 ! TYPE: REAL(fp)
323 ! DIMENSION: Rank-1 (N_PTS), or
324 ! Rank-2 (N_PTS,NPTS), or
325 ! Rank-3 (N_PTS,NPTS,N_PTS), or
326 ! Rank-4 (N_PTS,NPTS,N_PTS,NPTS), or
327 ! ATTRIBUTES: INTENT(IN)
328 !
329 ! ulp, vlp, wlp, xlp: Interpolating polynomial structures for the
330 ! respective dimensions from previous calls to
331 ! the LPoly() subroutine
332 ! UNITS: N/A
333 ! TYPE: LPoly_type
334 ! DIMENSION: Scalar
335 ! ATTRIBUTES: INTENT(IN)
336 !
337 ! z_TL: Tangent-linear interpolation data.
338 ! UNITS: Same as z
339 ! TYPE: REAL(fp)
340 ! DIMENSION: Same rank as z input
341 ! ATTRIBUTES: INTENT(IN)
342 !
343 ! ulp_TL, vlp_TL, Tangent-linear interpolating polynomial
344 ! wlp_TL, xlp_TL: structures for the respective dimensions
345 ! from previous calls to the LPoly_TL() subroutine.
346 ! UNITS: N/A
347 ! TYPE: LPoly_type
348 ! DIMENSION: Scalar
349 ! ATTRIBUTES: INTENT(IN)
350 !
351 ! OUTPUT ARGUMENTS:
352 ! z_int_TL: Tangent-linear interpolation result
353 ! UNITS: Same as z
354 ! TYPE: REAL(fp)
355 ! DIMENSION: Scalar
356 ! ATTRIBUTES: INTENT(IN OUT)
357 !
358 ! COMMENTS:
359 ! Output z_int_TL argument has INTENT(IN OUT) to prevent default
360 ! reinitialisation purely for computational speed.
361 !
362 !--------------------------------------------------------------------------------
363  ! 1-D routine
364  SUBROUTINE interp_1d_tl( z , ulp , & ! FWD Input
365  z_TL, ulp_TL, & ! TL Input
366  z_int_TL ) ! TL Output
367  ! Arguments
368  REAL(fp), INTENT(IN) :: z(:)
369  TYPE(lpoly_type), INTENT(IN) :: ulp
370  REAL(fp), INTENT(IN) :: z_tl(:)
371  TYPE(lpoly_type), INTENT(IN) :: ulp_tl
372  REAL(fp), INTENT(IN OUT) :: z_int_tl ! INTENT(IN OUT) to preclude reinitialisation
373  ! Perform TL interpolation
374  z_int_tl = ( ulp%w_left * ( ulp%lp_left(1) * z_tl(1) + &
375  ulp_tl%lp_left(1) * z(1) + &
376  ulp%lp_left(2) * z_tl(2) + &
377  ulp_tl%lp_left(2) * z(2) + &
378  ulp%lp_left(3) * z_tl(3) + &
379  ulp_tl%lp_left(3) * z(3) ) ) + &
380  ( ulp_tl%w_left * ( ulp%lp_left(1) * z(1) + &
381  ulp%lp_left(2) * z(2) + &
382  ulp%lp_left(3) * z(3) ) ) + &
383 
384  ( ulp%w_right * ( ulp%lp_right(1) * z_tl(2) + &
385  ulp_tl%lp_right(1) * z(2) + &
386  ulp%lp_right(2) * z_tl(3) + &
387  ulp_tl%lp_right(2) * z(3) + &
388  ulp%lp_right(3) * z_tl(4) + &
389  ulp_tl%lp_right(3) * z(4) ) ) + &
390  ( ulp_tl%w_right * ( ulp%lp_right(1) * z(2) + &
391  ulp%lp_right(2) * z(3) + &
392  ulp%lp_right(3) * z(4) ) )
393 
394  END SUBROUTINE interp_1d_tl
395 
396  ! 2-D routine
397  SUBROUTINE interp_2d_tl( z, ulp , vlp , & ! FWD Input
398  z_TL, ulp_TL, vlp_TL, & ! TL Input
399  z_int_TL ) ! TL Output
400  REAL(fp), INTENT(IN) :: z(:,:)
401  TYPE(lpoly_type), INTENT(IN) :: ulp, vlp
402  REAL(fp), INTENT(IN) :: z_tl(:,:)
403  TYPE(lpoly_type), INTENT(IN) :: ulp_tl, vlp_tl
404  REAL(fp), INTENT(IN OUT) :: z_int_tl ! INTENT(IN OUT) to preclude reinitialisation
405  ! Local variables
406  INTEGER :: i
407  REAL(fp) :: a(npts), a_tl(npts)
408  ! Interpolate z in v dimension for all w
409  DO i = 1, npts
410  CALL interp_1d(z(:,i),ulp,a(i))
411  CALL interp_1d_tl(z(:,i),ulp,z_tl(:,i),ulp_tl,a_tl(i))
412  END DO
413  ! Interpolate z in w dimension
414  CALL interp_1d_tl(a,vlp,a_tl,vlp_tl,z_int_tl)
415  END SUBROUTINE interp_2d_tl
416 
417  ! 3-D routine
418  SUBROUTINE interp_3d_tl( z , ulp , vlp , wlp , & ! FWD Input
419  z_TL, ulp_TL, vlp_TL, wlp_TL, & ! TL Input
420  z_int_TL ) ! TL Output
421  ! Arguments
422  REAL(fp), INTENT(IN) :: z(:,:,:)
423  TYPE(lpoly_type), INTENT(IN) :: ulp, vlp, wlp
424  REAL(fp), INTENT(IN) :: z_tl(:,:,:)
425  TYPE(lpoly_type), INTENT(IN) :: ulp_tl, vlp_tl, wlp_tl
426  REAL(fp), INTENT(IN OUT) :: z_int_tl ! INTENT(IN OUT) to preclude reinitialisation
427  ! Local variables
428  INTEGER :: i
429  REAL(fp) :: a(npts), a_tl(npts)
430  ! Interpolate z in u,v dimension for all w
431  DO i = 1, npts
432  CALL interp_2d(z(:,:,i),ulp,vlp,a(i))
433  CALL interp_2d_tl(z(:,:,i),ulp,vlp,z_tl(:,:,i),ulp_tl,vlp_tl,a_tl(i))
434  END DO
435  ! Interpolate a in w dimension
436  CALL interp_1d_tl(a,wlp,a_tl,wlp_tl,z_int_tl)
437  END SUBROUTINE interp_3d_tl
438 
439  ! 4-D routine
440  SUBROUTINE interp_4d_tl( z , ulp , vlp , wlp , xlp , & ! FWD Input
441  z_TL, ulp_TL, vlp_TL, wlp_TL, xlp_TL, & ! TL Input
442  z_int_TL ) ! TL Output
443  ! Arguments
444  REAL(fp), INTENT(IN) :: z(:,:,:,:)
445  TYPE(lpoly_type), INTENT(IN) :: ulp, vlp, wlp, xlp
446  REAL(fp), INTENT(IN) :: z_tl(:,:,:,:)
447  TYPE(lpoly_type), INTENT(IN) :: ulp_tl, vlp_tl, wlp_tl, xlp_tl
448  REAL(fp), INTENT(IN OUT) :: z_int_tl ! INTENT(IN OUT) to preclude reinitialisation
449  ! Local variables
450  INTEGER :: i
451  REAL(fp) :: a(npts), a_tl(npts)
452  ! Interpolate z in u,v,w dimension for all x
453  DO i = 1, npts
454  CALL interp_3d(z(:,:,:,i),ulp,vlp,wlp,a(i))
455  CALL interp_3d_tl(z(:,:,:,i),ulp,vlp,wlp,z_tl(:,:,:,i),ulp_tl,vlp_tl,wlp_tl,a_tl(i))
456  END DO
457  ! Interpolate a in x dimension
458  CALL interp_1d_tl(a,xlp,a_tl,xlp_tl,z_int_tl)
459  END SUBROUTINE interp_4d_tl
460 
461 
462 !--------------------------------------------------------------------------------
463 !
464 ! NAME:
465 ! Interp_1D_AD
466 ! Interp_2D_AD
467 ! Interp_3D_AD
468 ! Interp_4D_AD
469 !
470 ! PURPOSE:
471 ! Subroutines to perform adjoint interpolation:
472 ! o 1-D for z=f(u)
473 ! o 2-D for z=f(u,v)
474 ! o 3-D for z=f(u,v,w)
475 ! o 4-D for z=f(u,v,w,x)
476 !
477 ! CALLING SEQUENCE:
478 ! CALL Interp_1D_AD( z , ulp , & ! FWD Input
479 ! z_int_AD , & ! AD Input
480 ! z_AD, ulp_AD ) ! AD Output
481 !
482 ! CALL Interp_2D_AD( z , ulp , vlp , & ! FWD Input
483 ! z_int_AD , & ! AD Input
484 ! z_AD, ulp_AD, vlp_AD ) ! AD Output
485 !
486 ! CALL Interp_3D_AD( z , ulp , vlp , wlp , & ! FWD Input
487 ! z_int_AD , & ! AD Input
488 ! z_AD, ulp_AD, vlp_AD, wlp_AD ) ! AD Output
489 !
490 ! CALL Interp_4D_AD( z , ulp , vlp , wlp , xlp , & ! FWD Input
491 ! z_int_AD , & ! AD Input
492 ! z_AD, ulp_AD, vlp_AD, wlp_AD, xlp_AD ) ! AD Output
493 !
494 ! INPUT ARGUMENTS:
495 ! z: Data to interpolate
496 ! UNITS: Variable
497 ! TYPE: REAL(fp)
498 ! DIMENSION: Rank-1 (N_PTS), or
499 ! Rank-2 (N_PTS,NPTS), or
500 ! Rank-3 (N_PTS,NPTS,N_PTS), or
501 ! Rank-4 (N_PTS,NPTS,N_PTS,NPTS), or
502 ! ATTRIBUTES: INTENT(IN)
503 !
504 ! ulp, vlp, wlp, xlp: Interpolating polynomial structures for the
505 ! respective dimensions from previous calls to
506 ! the LPoly() subroutine
507 ! UNITS: N/A
508 ! TYPE: LPoly_type
509 ! DIMENSION: Scalar
510 ! ATTRIBUTES: INTENT(IN)
511 !
512 ! z_int_AD: Adjoint interpolate.
513 ! UNITS: N/A
514 ! TYPE: REAL(fp)
515 ! DIMENSION: Scalar
516 ! ATTRIBUTES: INTENT(IN OUT)
517 !
518 ! OUTPUT ARGUMENTS
519 ! z_AD: Adjoint interpolation data. Subsequently passed
520 ! into the LPoly_AD() subroutine.
521 ! UNITS: N/A
522 ! TYPE: REAL(fp)
523 ! DIMENSION: Same rank as z input
524 ! ATTRIBUTES: INTENT(IN OUT)
525 !
526 ! ulp_AD, vlp_AD, Adjoint interpolating polynomial
527 ! wlp_AD, xlp_AD: structures for respective dimensions.
528 ! Subsequently passed into the LPoly_AD()
529 ! subroutine.
530 ! UNITS: N/A
531 ! TYPE: LPoly_type
532 ! DIMENSION: Scalar
533 ! ATTRIBUTES: INTENT(IN OUT)
534 !
535 !--------------------------------------------------------------------------------
536  ! 1-D routine
537  SUBROUTINE interp_1d_ad( z , ulp , & ! FWD Input
538  z_int_AD , & ! AD Input
539  z_AD, ulp_AD ) ! AD Output
540  ! Arguments
541  REAL(fp), INTENT(IN) :: z(:)
542  TYPE(lpoly_type), INTENT(IN) :: ulp
543  REAL(fp), INTENT(IN OUT) :: z_int_ad
544  REAL(fp), INTENT(IN OUT) :: z_ad(:)
545  TYPE(lpoly_type), INTENT(IN OUT) :: ulp_ad
546  ! Local variables
547  REAL(fp) :: wl_z_int_ad, wr_z_int_ad
548  ! Perform adjoint interpolation
549  ulp_ad%w_right = ulp_ad%w_right + &
550  ( z_int_ad * ( ( ulp%lp_right(1) * z(2) ) + &
551  ( ulp%lp_right(2) * z(3) ) + &
552  ( ulp%lp_right(3) * z(4) ) ) )
553 
554  ulp_ad%w_left = ulp_ad%w_left + &
555  ( z_int_ad * ( ( ulp%lp_left(1) * z(1) ) + &
556  ( ulp%lp_left(2) * z(2) ) + &
557  ( ulp%lp_left(3) * z(3) ) ) )
558 
559  wr_z_int_ad = ulp%w_right * z_int_ad
560  ulp_ad%lp_right(1) = ulp_ad%lp_right(1) + ( wr_z_int_ad * z(2) )
561  ulp_ad%lp_right(2) = ulp_ad%lp_right(2) + ( wr_z_int_ad * z(3) )
562  ulp_ad%lp_right(3) = ulp_ad%lp_right(3) + ( wr_z_int_ad * z(4) )
563 
564  wl_z_int_ad = ulp%w_left * z_int_ad
565  ulp_ad%lp_left(1) = ulp_ad%lp_left(1) + ( wl_z_int_ad * z(1) )
566  ulp_ad%lp_left(2) = ulp_ad%lp_left(2) + ( wl_z_int_ad * z(2) )
567  ulp_ad%lp_left(3) = ulp_ad%lp_left(3) + ( wl_z_int_ad * z(3) )
568 
569  z_ad(1) = z_ad(1) + ( wl_z_int_ad * ulp%lp_left(1) )
570 
571  z_ad(2) = z_ad(2) + ( wr_z_int_ad * ulp%lp_right(1) ) + &
572  ( wl_z_int_ad * ulp%lp_left(2) )
573 
574  z_ad(3) = z_ad(3) + ( wr_z_int_ad * ulp%lp_right(2) ) + &
575  ( wl_z_int_ad * ulp%lp_left(3) )
576 
577  z_ad(4) = z_ad(4) + ( wr_z_int_ad * ulp%lp_right(3) )
578 
579  END SUBROUTINE interp_1d_ad
580 
581  ! 2-D routine
582  SUBROUTINE interp_2d_ad( z , ulp , vlp , & ! FWD Input
583  z_int_AD , & ! AD Input
584  z_AD, ulp_AD, vlp_AD ) ! AD Output
585  ! Arguments
586  REAL(fp), INTENT(IN) :: z(:,:)
587  TYPE(lpoly_type), INTENT(IN) :: ulp, vlp
588  REAL(fp), INTENT(IN OUT) :: z_int_ad
589  REAL(fp), INTENT(IN OUT) :: z_ad(:,:)
590  TYPE(lpoly_type), INTENT(IN OUT) :: ulp_ad, vlp_ad
591  ! Local variables
592  INTEGER :: i
593  REAL(fp) :: a(npts), a_ad(npts)
594  ! Forward calculations
595  ! Interpolate z in u dimension for all v
596  DO i = 1, npts
597  CALL interp_1d(z(:,i),ulp,a(i))
598  END DO
599  ! Adjoint calculations
600  ! Initialize local AD variables
601  a_ad = zero
602  ! Adjoint of z interpolation in v dimension
603  CALL interp_1d_ad(a,vlp,z_int_ad,a_ad,vlp_ad)
604  ! Adjoint of z interpolation in u dimension for all v
605  DO i = 1, npts
606  CALL interp_1d_ad(z(:,i),ulp,a_ad(i),z_ad(:,i),ulp_ad)
607  END DO
608  END SUBROUTINE interp_2d_ad
609 
610  ! 3-D routine
611  SUBROUTINE interp_3d_ad( z , ulp , vlp , wlp , & ! FWD Input
612  z_int_AD , & ! AD Input
613  z_AD, ulp_AD, vlp_AD, wlp_AD ) ! AD Output
614  ! Arguments
615  REAL(fp), INTENT(IN) :: z(:,:,:)
616  TYPE(lpoly_type), INTENT(IN) :: ulp, vlp, wlp
617  REAL(fp), INTENT(IN OUT) :: z_int_ad
618  REAL(fp), INTENT(IN OUT) :: z_ad(:,:,:)
619  TYPE(lpoly_type), INTENT(IN OUT) :: ulp_ad, vlp_ad, wlp_ad
620  ! Local variables
621  INTEGER :: i
622  REAL(fp) :: a(npts), a_ad(npts)
623 
624  ! Forward calculations
625  ! Interpolate z in u and v dimension for all w
626  DO i = 1, npts
627  CALL interp_2d(z(:,:,i),ulp,vlp,a(i))
628  END DO
629 
630  ! Adjoint calculations
631  ! Initialize local AD variables
632  a_ad = zero
633  ! Adjoint of a interpolation in w dimension
634  CALL interp_1d_ad(a,wlp,z_int_ad,a_ad,wlp_ad)
635  ! Adjoint of z interpolation in u and v dimension for all w
636  DO i = 1, npts
637  CALL interp_2d_ad(z(:,:,i),ulp,vlp,a_ad(i),z_ad(:,:,i),ulp_ad,vlp_ad)
638  END DO
639  END SUBROUTINE interp_3d_ad
640 
641  ! 4-D routine
642  SUBROUTINE interp_4d_ad( z , ulp , vlp , wlp , xlp , & ! FWD Input
643  z_int_AD , & ! AD Input
644  z_AD, ulp_AD, vlp_AD, wlp_AD, xlp_AD ) ! AD Output
645  ! Arguments
646  REAL(fp), INTENT(IN) :: z(:,:,:,:)
647  TYPE(lpoly_type), INTENT(IN) :: ulp, vlp, wlp, xlp
648  REAL(fp), INTENT(IN OUT) :: z_int_ad
649  REAL(fp), INTENT(IN OUT) :: z_ad(:,:,:,:)
650  TYPE(lpoly_type), INTENT(IN OUT) :: ulp_ad, vlp_ad, wlp_ad, xlp_ad
651  ! Local variables
652  INTEGER :: i
653  REAL(fp) :: a(npts), a_ad(npts)
654 
655  ! Forward calculations
656  ! Interpolate z in u,v,w dimension for all x
657  DO i = 1, npts
658  CALL interp_3d(z(:,:,:,i),ulp,vlp,wlp,a(i))
659  END DO
660 
661  ! Adjoint calculations
662  ! Initialize local AD variables
663  a_ad = zero
664  ! Adjoint of a interpolation in x dimension
665  CALL interp_1d_ad(a,xlp,z_int_ad,a_ad,xlp_ad)
666  ! Adjoint of z interpolation in u,v,w dimension for all x
667  DO i = 1, npts
668  CALL interp_3d_ad(z(:,:,:,i),ulp,vlp,wlp,a_ad(i),z_ad(:,:,:,i),ulp_ad,vlp_ad,wlp_ad)
669  END DO
670  END SUBROUTINE interp_4d_ad
671 
672 
673 !--------------------------------------------------------------------------------
674 !
675 ! NAME:
676 ! Find_Index
677 !
678 ! PURPOSE:
679 ! Subroutines to search abscissa data for 4-pt interplation indices.
680 !
681 ! CALLING SEQUENCE:
682 ! For regularly spaced x-data:
683 ! CALL Find_Index( x, dx, & ! Input
684 ! x_int, & ! In/Output
685 ! i1, i2, & ! Output
686 ! out_of_bounds ) ! Output
687 !
688 ! For irregularly spaced x-data:
689 ! CALL Find_Index( x, & ! Input
690 ! x_int, & ! In/Output
691 ! i1, i2, & ! Output
692 ! out_of_bounds ) ! Output
693 !
694 ! INPUTS:
695 ! x: Abscissa data.
696 ! UNITS: Variable
697 ! TYPE: REAL(fp)
698 ! DIMENSION: Rank-1 (N)
699 ! ATTRIBUTES: INTENT(IN)
700 !
701 ! dx: Abscissa data spacing for the regularly spaced case.
702 ! UNITS: Same as x.
703 ! TYPE: REAL(fp)
704 ! DIMENSION: Scalar
705 ! ATTRIBUTES: INTENT(IN)
706 !
707 ! INPUTS/OUTPUTS
708 ! x_int: On input : Abscissa value at which an interpolate
709 ! is desired.
710 ! On output: Valid abscissa value at which an interpolate
711 ! will be computed.
712 ! If x_int < x(1), then x_int = x(1) upon exit.
713 ! x_int > x(N), then x_int = x(N) upon exit
714 ! x(1) <= x_int <= x(N), then x_int is unchanged upon exit.
715 ! Also see the out_of_bounds output argument.
716 ! UNITS: Same as x.
717 ! TYPE: REAL(fp)
718 ! DIMENSION: Scalar
719 ! ATTRIBUTES: INTENT(IN OUT)
720 !
721 ! OUTPUTS
722 ! i1, i2: Begin and end indices in the input x-array to
723 ! use for the 4-pt interpolation at the value x_int.
724 ! Three cases are possible for a x array of length N:
725 ! Normal : x(i1) < x(i1+1) <= x_int <= x(i1+2) < x(i1+3)
726 ! Left edge : x_int < x(2); then i1=1, i2=4
727 ! Right edge: x_int > x(N-1); then i1=N-3, i2=N
728 ! UNITS: N/A
729 ! TYPE: INTEGER
730 ! DIMENSION: Scalar
731 ! ATTRIBUTES: INTENT(IN OUT)
732 !
733 ! out_of_bounds: Logical variable that identifies if the interpolate point,
734 ! x_int, is within the bounds of the search array, x.
735 ! If x(1) <= x_int <= x(n), out_of_bounds == .FALSE.
736 ! Else out_of_bounds == .TRUE.
737 ! UNITS: N/A
738 ! TYPE: LOGICAL
739 ! DIMENSION: Scalar
740 ! ATTRIBUTES: INTENT(IN OUT)
741 !
742 ! COMMENTS:
743 ! Output arguments have INTENT(IN OUT) to prevent default
744 ! reinitialisation purely for computational speed.
745 !
746 !--------------------------------------------------------------------------------
747  ! Find indices for regular spacing
748  SUBROUTINE find_regular_index(x, dx, x_int, i1, i2, out_of_bounds)
749  REAL(fp), INTENT(IN) :: x(:)
750  REAL(fp), INTENT(IN) :: dx
751  REAL(fp), INTENT(IN OUT) :: x_int
752  INTEGER , INTENT(IN OUT) :: i1, i2
753  LOGICAL , INTENT(IN OUT) :: out_of_bounds
754  INTEGER :: n
755  n = SIZE(x)
756  out_of_bounds = .false.
757  IF ( x_int < x(1) .OR. x_int > x(n) ) out_of_bounds = .true.
758  i1 = floor((x_int-x(1))/dx)+1-(npoly_pts/2)
759  i1 = min(max(i1,1),n-npoly_pts)
760  i2 = i1 + npoly_pts
761  IF (out_of_bounds .AND. i1==1) THEN
762  x_int = x(1)
763  ELSE IF (out_of_bounds .AND. i2==n) THEN
764  x_int = x(n)
765  END IF
766  END SUBROUTINE find_regular_index
767 
768  ! Find indices for random spacing.
769  ! Assumption is that x(1) <= xInt <= x(n)
770  ! (despite the MIN/MAX test)
771  SUBROUTINE find_random_index(x, x_int, i1, i2, out_of_bounds)
772  REAL(fp), INTENT(IN) :: x(:)
773  REAL(fp), INTENT(IN OUT) :: x_int
774  INTEGER , INTENT(IN OUT) :: i1, i2
775  LOGICAL , INTENT(IN OUT) :: out_of_bounds
776  INTEGER :: k, n
777  n = SIZE(x)
778  out_of_bounds = .false.
779  IF ( x_int < x(1) .OR. x_int > x(n) ) out_of_bounds = .true.
780  DO k=1,n
781  IF (x_int <= x(k) ) EXIT
782  END DO
783  i1 = min(max(1,k-1-(npoly_pts/2)),n-npoly_pts)
784  i2 = i1 + npoly_pts
785  IF (out_of_bounds .AND. i1==1) THEN
786  x_int = x(1)
787  ELSE IF (out_of_bounds .AND. i2==n) THEN
788  x_int = x(n)
789  END IF
790  END SUBROUTINE find_random_index
791 
792 
793 !--------------------------------------------------------------------------------
794 !
795 ! NAME:
796 ! LPoly
797 !
798 ! PURPOSE:
799 ! Subroutines to compute the Lagrangian polynomial terms for interpolation.
800 !
801 ! CALLING SEQUENCE:
802 ! CALL LPoly( x, x_int, & ! Input
803 ! p ) ! Output
804 !
805 !
806 ! INPUT ARGUMENTS:
807 ! x: 4-pt abscissa data.
808 ! UNITS: Variable
809 ! TYPE: REAL(fp)
810 ! DIMENSION: Rank-1 (4)
811 ! ATTRIBUTES: INTENT(IN)
812 !
813 ! x_int: Abscissa value at which an interpolate is desired.
814 ! Typically, x(2) <= xInt <= x(3), except for the data edges.
815 ! UNITS: Same as x.
816 ! TYPE: REAL(fp)
817 ! DIMENSION: Scalar
818 ! ATTRIBUTES: INTENT(IN)
819 !
820 ! OUTPUT ARGUMENTS:
821 ! p: Lagrangian polynomial structure for subsequent use in the
822 ! various interpolation subroutines.
823 ! UNITS: N/A
824 ! TYPE: TYPE(LPoly_type)
825 ! DIMENSION: Scalar
826 ! ATTRIBUTES: INTENT(IN OUT)
827 !
828 ! COMMENTS:
829 ! Output p argument is INTENT(IN OUT) to prevent default reinitialisation
830 ! purely for computational speed.
831 !
832 !--------------------------------------------------------------------------------
833 
834  SUBROUTINE lpoly(x, x_int, p)
835  REAL(fp), INTENT(IN) :: x(:) ! Input
836  REAL(fp), INTENT(IN) :: x_int ! Input
837  TYPE(lpoly_type), INTENT(IN OUT) :: p ! Output. INTENT(IN OUT) to preclude reinitialisation
838 
839  ! Compute the numerator differences
840  CALL compute_dxi(x(1:3),x_int,p%dxi_left)
841  CALL compute_dxi(x(2:4),x_int,p%dxi_right)
842 
843  ! Compute the denominator differences
844  CALL compute_dx(x(1:3),p%dx_left)
845  CALL compute_dx(x(2:4),p%dx_right)
846 
847  ! Compute the quadratic polynomials
848  CALL compute_qpoly(p%dxi_left , p%dx_left , p%lp_left)
849  CALL compute_qpoly(p%dxi_right, p%dx_right, p%lp_right)
850 
851  ! Polynomial weights
852  IF ( x_int < x(2) ) THEN
853  p%w_right = zero
854  p%w_left = one
855  ELSE IF ( x_int > x(3) ) THEN
856  p%w_right = one
857  p%w_left = zero
858  ELSE
859  p%w_right = p%dxi_left(2) / (-p%dx_left(3))
860  p%w_left = one - p%w_right
861  END IF
862  END SUBROUTINE lpoly
863 
864 
865 !--------------------------------------------------------------------------------
866 !
867 ! NAME:
868 ! LPoly_TL
869 !
870 ! PURPOSE:
871 ! Subroutines to compute the tangent-linear Lagrangian polynomial terms
872 ! for interpolation.
873 !
874 ! CALLING SEQUENCE:
875 ! CALL LPoly_TL( x , x_int , & ! FWD Input
876 ! p , & ! FWD Input
877 ! x_TL, x_int_TL, & ! TL Input
878 ! p_TL ) ! TL Output
879 !
880 !
881 ! INPUT ARGUMENTS:
882 ! x: 4-pt abscissa data.
883 ! UNITS: Variable
884 ! TYPE: REAL(fp)
885 ! DIMENSION: Rank-1 (4)
886 ! ATTRIBUTES: INTENT(IN)
887 !
888 ! x_int: Abscissa value at which an interpolate is desired.
889 ! Typically, x(2) <= xInt <= x(3), except for the data edges.
890 ! UNITS: Same as x.
891 ! TYPE: REAL(fp)
892 ! DIMENSION: Scalar
893 ! ATTRIBUTES: INTENT(IN)
894 !
895 ! p: Lagrangian polynomial structure from previous call to
896 ! the LPoly() subroutine.
897 ! UNITS: N/A
898 ! TYPE: TYPE(LPoly_type)
899 ! DIMENSION: Scalar
900 ! ATTRIBUTES: INTENT(IN)
901 !
902 ! x_TL: Tangent-linear 4-pt abscissa data.
903 ! UNITS: Same as x.
904 ! TYPE: REAL(fp)
905 ! DIMENSION: Rank-1 (4)
906 ! ATTRIBUTES: INTENT(IN)
907 !
908 ! x_int_TL: Tangent-linear abscissa value at which an interpolate
909 ! is desired.
910 ! UNITS: Same as x.
911 ! TYPE: REAL(fp)
912 ! DIMENSION: Scalar
913 ! ATTRIBUTES: INTENT(IN)
914 !
915 ! OUTPUT ARGUMENTS:
916 ! p_TL: Tangent-linear Lagrangian polynomial structure for subsequent
917 ! use in the various tangent-linear interpolation subroutines.
918 ! UNITS: N/A
919 ! TYPE: TYPE(LPoly_type)
920 ! DIMENSION: Scalar
921 ! ATTRIBUTES: INTENT(IN OUT)
922 !
923 ! COMMENTS:
924 ! Output p_TL argument is INTENT(IN OUT) to prevent default reinitialisation
925 ! purely for computational speed.
926 !
927 !--------------------------------------------------------------------------------
928 
929  SUBROUTINE lpoly_tl(x , x_int , p , &
930  x_TL, x_int_TL, p_TL)
931  REAL(fp), INTENT(IN) :: x(:) ! FWD Input
932  REAL(fp), INTENT(IN) :: x_int ! FWD Input
933  TYPE(lpoly_type), INTENT(IN) :: p ! FWD Input
934  REAL(fp), INTENT(IN) :: x_tl(:) ! TL Input
935  REAL(fp), INTENT(IN) :: x_int_tl ! TL Input
936  TYPE(lpoly_type), INTENT(IN OUT) :: p_tl ! TL Output. INTENT(IN OUT) to preclude reinitialisation
937 
938  ! Compute the tangent-linear numerator differences
939  CALL compute_dxi_tl(x_tl(1:3),x_int_tl,p_tl%dxi_left)
940  CALL compute_dxi_tl(x_tl(2:4),x_int_tl,p_tl%dxi_right)
941 
942  ! Compute the tangent-linear denominator differences
943  CALL compute_dx_tl(x_tl(1:3),p_tl%dx_left)
944  CALL compute_dx_tl(x_tl(2:4),p_tl%dx_right)
945 
946  ! Compute the tangent-linear quadratic polynomials
947  CALL compute_qpoly_tl(p%dxi_left , p%dx_left , p%lp_left, &
948  p_tl%dxi_left, p_tl%dx_left , p_tl%lp_left)
949  CALL compute_qpoly_tl(p%dxi_right , p%dx_right , p%lp_right, &
950  p_tl%dxi_right, p_tl%dx_right, p_tl%lp_right)
951 
952  ! Polynomial weights
953  IF ( x_int < x(2) .OR. x_int > x(3) ) THEN
954  p_tl%w_right = zero
955  p_tl%w_left = zero
956  ELSE
957  p_tl%w_right = -( p_tl%dxi_left(2) + (p%w_right*p_tl%dx_left(3)) ) / p%dx_left(3)
958  p_tl%w_left = -p_tl%w_right
959  END IF
960  END SUBROUTINE lpoly_tl
961 
962 
963 !--------------------------------------------------------------------------------
964 !
965 ! NAME:
966 ! LPoly_AD
967 !
968 ! PURPOSE:
969 ! Subroutines to compute the Lagrangian polynomial adjoint terms
970 ! for interpolation.
971 !
972 ! CALLING SEQUENCE:
973 ! CALL LPoly_AD( x , x_int , & ! FWD Input
974 ! p , & ! FWD Input
975 ! p_AD , & ! AD Input
976 ! x_AD, x_int_AD ) ! AD Output
977 !
978 ! INPUT ARGUMENTS:
979 ! x: 4-pt abscissa data.
980 ! UNITS: Variable
981 ! TYPE: REAL(fp)
982 ! DIMENSION: Rank-1 (4)
983 ! ATTRIBUTES: INTENT(IN)
984 !
985 ! x_int: Abscissa value at which an interpolate is desired.
986 ! Typically, x(2) <= xInt <= x(3), except for the data edges.
987 ! UNITS: Same as x.
988 ! TYPE: REAL(fp)
989 ! DIMENSION: Scalar
990 ! ATTRIBUTES: INTENT(IN)
991 !
992 ! p: Lagrangian polynomial structure from previous call to
993 ! the LPoly() subroutine.
994 ! UNITS: N/A
995 ! TYPE: TYPE(LPoly_type)
996 ! DIMENSION: Scalar
997 ! ATTRIBUTES: INTENT(IN)
998 !
999 ! p_AD: Adjoint Lagrangian polynomial structure from previous calls
1000 ! to the various adjoint interpolation subroutines.
1001 ! *Note*: Components are modified (zeroed out) in this routine.
1002 ! UNITS: N/A
1003 ! TYPE: TYPE(LPoly_type)
1004 ! DIMENSION: Scalar
1005 ! ATTRIBUTES: INTENT(IN OUT)
1006 !
1007 ! OUTPUT ARGUMENTS:
1008 ! x_AD: Adjoint 4-pt abscissa data.
1009 ! UNITS: Variable
1010 ! TYPE: REAL(fp)
1011 ! DIMENSION: Rank-1 (4)
1012 ! ATTRIBUTES: INTENT(IN OUT)
1013 !
1014 ! x_int_AD: Adjoint abscissa value at which an interpolate is desired.
1015 ! UNITS: Same as x_AD.
1016 ! TYPE: REAL(fp)
1017 ! DIMENSION: Scalar
1018 ! ATTRIBUTES: INTENT(IN OUT)
1019 !
1020 ! SIDE EFFECTS:
1021 ! The adjoint input argument, p_AD, is modified in this subroutine. Its
1022 ! components are reinitialised (zeroed out).
1023 !
1024 !--------------------------------------------------------------------------------
1025 
1026  SUBROUTINE lpoly_ad(x , x_int , p , &
1027  p_AD, x_AD, x_int_AD)
1028  REAL(fp), INTENT(IN) :: x(:) ! FWD Input
1029  REAL(fp), INTENT(IN) :: x_int ! FWD Input
1030  TYPE(lpoly_type), INTENT(IN) :: p ! FWD Input
1031  TYPE(lpoly_type), INTENT(IN OUT) :: p_ad ! AD Input
1032  REAL(fp), INTENT(IN OUT) :: x_ad(:) ! AD Output
1033  REAL(fp), INTENT(IN OUT) :: x_int_ad ! AD Output
1034 
1035  ! Polynomial weights
1036  IF ( x_int < x(2) .OR. x_int > x(3) ) THEN
1037  p_ad%w_right = zero
1038  p_ad%w_left = zero
1039  ELSE
1040  p_ad%w_right = p_ad%w_right - p_ad%w_left
1041  p_ad%w_left = zero
1042  p_ad%dx_left(3) = p_ad%dx_left(3) - (p%w_right*p_ad%w_right/p%dx_left(3))
1043  p_ad%dxi_left(2) = p_ad%dxi_left(2) - (p_ad%w_right/p%dx_left(3))
1044  p_ad%w_right = zero
1045  END IF
1046 
1047  ! "Right" side quadratic
1048  CALL compute_qpoly_ad(p%dxi_right , p%dx_right, &
1049  p%lp_right , &
1050  p_ad%lp_right , &
1051  p_ad%dxi_right, p_ad%dx_right)
1052 
1053  ! "Left" side quadratic
1054  CALL compute_qpoly_ad(p%dxi_left , p%dx_left, &
1055  p%lp_left , &
1056  p_ad%lp_left , &
1057  p_ad%dxi_left, p_ad%dx_left)
1058 
1059  ! Compute the adjoint denominator differences
1060  CALL compute_dx_ad(p_ad%dx_right,x_ad(2:4))
1061  CALL compute_dx_ad(p_ad%dx_left ,x_ad(1:3))
1062 
1063  ! Compute the adjoint numerator differences
1064  CALL compute_dxi_ad(p_ad%dxi_right,x_ad(2:4),x_int_ad)
1065  CALL compute_dxi_ad(p_ad%dxi_left ,x_ad(1:3),x_int_ad)
1066  END SUBROUTINE lpoly_ad
1067 
1068 
1069 !################################################################################
1070 !################################################################################
1071 !## ##
1072 !## ## PRIVATE MODULE ROUTINES ## ##
1073 !## ##
1074 !################################################################################
1075 !################################################################################
1076 
1077  ! ------------------------------------------------
1078  ! Subroutines to compute the quadratic polynomials
1079  ! ------------------------------------------------
1080  ! Forward model
1081  SUBROUTINE compute_qpoly(dxi, dx, lp)
1082  REAL(fp), INTENT(IN) :: dxi(:) ! Input
1083  REAL(fp), INTENT(IN) :: dx(:) ! Input
1084  REAL(fp), INTENT(IN OUT) :: lp(:) ! Output. INTENT(IN OUT) to preclude reinitialisation
1085  lp(1) = dxi(2)*dxi(3) / ( dx(1)*dx(2))
1086  lp(2) = dxi(1)*dxi(3) / (-dx(1)*dx(3))
1087  lp(3) = dxi(1)*dxi(2) / ( dx(2)*dx(3))
1088  END SUBROUTINE compute_qpoly
1089 
1090  ! Tangent-linear model
1091  SUBROUTINE compute_qpoly_tl(dxi , dx , lp, &
1092  dxi_TL, dx_TL, lp_TL)
1093  REAL(fp), INTENT(IN) :: dxi(:) ! FWD Input
1094  REAL(fp), INTENT(IN) :: dx(:) ! FWD Input
1095  REAL(fp), INTENT(IN) :: lp(:) ! FWD Input
1096  REAL(fp), INTENT(IN) :: dxi_TL(:) ! TL Input
1097  REAL(fp), INTENT(IN) :: dx_TL(:) ! TL Input
1098  REAL(fp), INTENT(IN OUT) :: lp_TL(:) ! TL Output. INTENT(IN OUT) to preclude reinitialisation
1099  lp_tl(1) = ( (dxi(3)*dxi_tl(2) ) + &
1100  (dxi(2)*dxi_tl(3) ) - &
1101  (dx(2) *dx_tl(1) *lp(1)) - &
1102  (dx(1) *dx_tl(2) *lp(1)) ) / (dx(1)*dx(2))
1103  lp_tl(2) = ( (dxi(3)*dxi_tl(1) ) + &
1104  (dxi(1)*dxi_tl(3) ) + &
1105  (dx(3) *dx_tl(1) *lp(2)) + &
1106  (dx(1) *dx_tl(3) *lp(2)) ) / (-dx(1)*dx(3))
1107  lp_tl(3) = ( (dxi(2)*dxi_tl(1) ) + &
1108  (dxi(1)*dxi_tl(2) ) - &
1109  (dx(3) *dx_tl(2) *lp(3)) - &
1110  (dx(2) *dx_tl(3) *lp(3)) ) / (dx(2)*dx(3))
1111  END SUBROUTINE compute_qpoly_tl
1112 
1113  ! Adjoint model
1114  SUBROUTINE compute_qpoly_ad(dxi , dx, &
1115  lp , &
1116  lp_AD , &
1117  dxi_AD, dx_AD)
1118  REAL(fp), INTENT(IN) :: dxi(:) ! FWD Input
1119  REAL(fp), INTENT(IN) :: dx(:) ! FWD Input
1120  REAL(fp), INTENT(IN) :: lp(:) ! FWD Input
1121  REAL(fp), INTENT(IN OUT) :: lp_AD(:) ! AD Input
1122  REAL(fp), INTENT(IN OUT) :: dxi_AD(:) ! AD Output
1123  REAL(fp), INTENT(IN OUT) :: dx_AD(:) ! AD Output
1124  REAL(fp) :: d
1125  ! Adjoint of lp(3)
1126  d = lp_ad(3)/(dx(2)*dx(3))
1127  dxi_ad(1) = dxi_ad(1) + d*dxi(2)
1128  dxi_ad(2) = dxi_ad(2) + d*dxi(1)
1129  dx_ad(2) = dx_ad(2) - d*dx(3)*lp(3)
1130  dx_ad(3) = dx_ad(3) - d*dx(2)*lp(3)
1131  lp_ad(3) = zero
1132  ! Adjoint of lp(2)
1133  d = lp_ad(2)/(-dx(1)*dx(3))
1134  dxi_ad(1) = dxi_ad(1) + d*dxi(3)
1135  dxi_ad(3) = dxi_ad(3) + d*dxi(1)
1136  dx_ad(2) = dx_ad(2) + d*dx(3)*lp(2)
1137  dx_ad(3) = dx_ad(3) + d*dx(2)*lp(2)
1138  lp_ad(2) = zero
1139  ! Adjoint of lp(1)
1140  d = lp_ad(1)/(dx(1)*dx(2))
1141  dxi_ad(2) = dxi_ad(2) + d*dxi(3)
1142  dxi_ad(3) = dxi_ad(3) + d*dxi(2)
1143  dx_ad(1) = dx_ad(1) - d*dx(2)*lp(1)
1144  dx_ad(2) = dx_ad(2) - d*dx(1)*lp(1)
1145  lp_ad(1) = zero
1146  END SUBROUTINE compute_qpoly_ad
1147 
1148 
1149 
1150  ! -------------------------------------------------------------
1151  ! Subroutines to compute the polynomial denominator differences
1152  ! -------------------------------------------------------------
1153  ! Forward model
1154  SUBROUTINE compute_dx(x,dx)
1155  REAL(fp), INTENT(IN) :: x(:) ! Input
1156  REAL(fp), INTENT(IN OUT) :: dx(:) ! Output. INTENT(IN OUT) to preclude reinitialisation
1157  dx(1) = x(1)-x(2)
1158  dx(2) = x(1)-x(3)
1159  dx(3) = x(2)-x(3)
1160  END SUBROUTINE compute_dx
1161 
1162  ! Tangent-linear model
1163  SUBROUTINE compute_dx_tl(x_TL,dx_TL)
1164  REAL(fp), INTENT(IN) :: x_TL(:) ! TL Input
1165  REAL(fp), INTENT(IN OUT) :: dx_TL(:) ! TL Output. INTENT(IN OUT) to preclude reinitialisation
1166  dx_tl(1) = x_tl(1)-x_tl(2)
1167  dx_tl(2) = x_tl(1)-x_tl(3)
1168  dx_tl(3) = x_tl(2)-x_tl(3)
1169  END SUBROUTINE compute_dx_tl
1170 
1171  ! Adjoint model
1172  SUBROUTINE compute_dx_ad(dx_AD,x_AD)
1173  REAL(fp), INTENT(IN OUT) :: dx_AD(:) ! AD Input
1174  REAL(fp), INTENT(IN OUT) :: x_AD(:) ! AD Output
1175  x_ad(3) = x_ad(3) - dx_ad(3)
1176  x_ad(2) = x_ad(2) + dx_ad(3)
1177  dx_ad(3) = zero
1178  x_ad(3) = x_ad(3) - dx_ad(2)
1179  x_ad(1) = x_ad(1) + dx_ad(2)
1180  dx_ad(2) = zero
1181  x_ad(2) = x_ad(2) - dx_ad(1)
1182  x_ad(1) = x_ad(1) + dx_ad(1)
1183  dx_ad(1) = zero
1184  END SUBROUTINE compute_dx_ad
1185 
1186 
1187  ! -----------------------------------------------------------
1188  ! Subroutines to compute the polynomial numerator differences
1189  ! -----------------------------------------------------------
1190  ! Forward model
1191  SUBROUTINE compute_dxi(x,xi,dxi)
1192  REAL(fp), INTENT(IN) :: x(:) ! Input
1193  REAL(fp), INTENT(IN) :: xi ! Input
1194  REAL(fp), INTENT(IN OUT) :: dxi(:) ! Output. INTENT(IN OUT) to preclude reinitialisation
1195  dxi(1) = xi - x(1)
1196  dxi(2) = xi - x(2)
1197  dxi(3) = xi - x(3)
1198  END SUBROUTINE compute_dxi
1199 
1200  ! Tangent-linear model
1201  SUBROUTINE compute_dxi_tl(x_TL,xi_TL,dxi_TL)
1202  REAL(fp), INTENT(IN) :: x_TL(:) ! TL Input
1203  REAL(fp), INTENT(IN) :: xi_TL ! TL Input
1204  REAL(fp), INTENT(IN OUT) :: dxi_TL(:) ! TL Output. INTENT(IN OUT) to preclude reinitialisation
1205  dxi_tl(1) = xi_tl - x_tl(1)
1206  dxi_tl(2) = xi_tl - x_tl(2)
1207  dxi_tl(3) = xi_tl - x_tl(3)
1208  END SUBROUTINE compute_dxi_tl
1209 
1210  ! Adjoint model
1211  SUBROUTINE compute_dxi_ad(dxi_AD,x_AD,xi_AD)
1212  REAL(fp), INTENT(IN OUT) :: dxi_AD(:) ! AD Input
1213  REAL(fp), INTENT(IN OUT) :: x_AD(:) ! AD Output
1214  REAL(fp), INTENT(IN OUT) :: xi_AD ! AD Output
1215  x_ad(1) = x_ad(1) - dxi_ad(1)
1216  xi_ad = xi_ad + dxi_ad(1)
1217  dxi_ad(1) = zero
1218  x_ad(2) = x_ad(2) - dxi_ad(2)
1219  xi_ad = xi_ad + dxi_ad(2)
1220  dxi_ad(2) = zero
1221  x_ad(3) = x_ad(3) - dxi_ad(3)
1222  xi_ad = xi_ad + dxi_ad(3)
1223  dxi_ad(3) = zero
1224  END SUBROUTINE compute_dxi_ad
1225 
1226 END MODULE crtm_interpolation
1227 
crtm_interpolation::interp_3d_ad
subroutine, public interp_3d_ad(z, ulp, vlp, wlp, z_int_AD, z_AD, ulp_AD, vlp_AD, wlp_AD)
Definition: CRTM_Interpolation.f90:614
crtm_interpolation::interp_3d_tl
subroutine, public interp_3d_tl(z, ulp, vlp, wlp, z_TL, ulp_TL, vlp_TL, wlp_TL, z_int_TL)
Definition: CRTM_Interpolation.f90:421
crtm_interpolation::find_regular_index
subroutine find_regular_index(x, dx, x_int, i1, i2, out_of_bounds)
Definition: CRTM_Interpolation.f90:749
crtm_interpolation::module_rcs_id
character(*), parameter module_rcs_id
Definition: CRTM_Interpolation.f90:66
crtm_interpolation::interp_4d
subroutine, public interp_4d(z, ulp, vlp, wlp, xlp, z_int)
Definition: CRTM_Interpolation.f90:282
crtm_interpolation::compute_dxi_ad
subroutine compute_dxi_ad(dxi_AD, x_AD, xi_AD)
Definition: CRTM_Interpolation.f90:1212
type_kinds::fp
integer, parameter, public fp
Definition: Type_Kinds.f90:124
crtm_interpolation::interp_4d_tl
subroutine, public interp_4d_tl(z, ulp, vlp, wlp, xlp, z_TL, ulp_TL, vlp_TL, wlp_TL, xlp_TL, z_int_TL)
Definition: CRTM_Interpolation.f90:443
crtm_interpolation::compute_qpoly_ad
subroutine compute_qpoly_ad(dxi, dx, lp, lp_AD, dxi_AD, dx_AD)
Definition: CRTM_Interpolation.f90:1118
crtm_interpolation::compute_dxi
subroutine compute_dxi(x, xi, dxi)
Definition: CRTM_Interpolation.f90:1192
crtm_interpolation::interp_1d_ad
subroutine, public interp_1d_ad(z, ulp, z_int_AD, z_AD, ulp_AD)
Definition: CRTM_Interpolation.f90:540
crtm_interpolation
Definition: CRTM_Interpolation.f90:12
crtm_interpolation::interp_4d_ad
subroutine, public interp_4d_ad(z, ulp, vlp, wlp, xlp, z_int_AD, z_AD, ulp_AD, vlp_AD, wlp_AD, xlp_AD)
Definition: CRTM_Interpolation.f90:645
crtm_interpolation::clear_lpoly
subroutine, public clear_lpoly(p)
Definition: CRTM_Interpolation.f90:132
crtm_interpolation::lpoly_ad
subroutine, public lpoly_ad(x, x_int, p, p_AD, x_AD, x_int_AD)
Definition: CRTM_Interpolation.f90:1028
crtm_interpolation::compute_qpoly
subroutine compute_qpoly(dxi, dx, lp)
Definition: CRTM_Interpolation.f90:1082
crtm_interpolation::lpoly_type
Definition: CRTM_Interpolation.f90:78
crtm_interpolation::lpoly_inspect
subroutine, public lpoly_inspect(self)
Definition: CRTM_Interpolation.f90:153
crtm_interpolation::interp_1d_tl
subroutine, public interp_1d_tl(z, ulp, z_TL, ulp_TL, z_int_TL)
Definition: CRTM_Interpolation.f90:367
crtm_interpolation::find_index
Definition: CRTM_Interpolation.f90:57
crtm_interpolation::zero
real(fp), parameter zero
Definition: CRTM_Interpolation.f90:68
crtm_interpolation::lpoly_init
subroutine, public lpoly_init(self)
Definition: CRTM_Interpolation.f90:147
crtm_interpolation::interp_2d_ad
subroutine, public interp_2d_ad(z, ulp, vlp, z_int_AD, z_AD, ulp_AD, vlp_AD)
Definition: CRTM_Interpolation.f90:585
crtm_interpolation::lpoly
subroutine, public lpoly(x, x_int, p)
Definition: CRTM_Interpolation.f90:835
crtm_interpolation::npts
integer, parameter, public npts
Definition: CRTM_Interpolation.f90:72
crtm_interpolation::compute_dx
subroutine compute_dx(x, dx)
Definition: CRTM_Interpolation.f90:1155
crtm_interpolation::compute_qpoly_tl
subroutine compute_qpoly_tl(dxi, dx, lp, dxi_TL, dx_TL, lp_TL)
Definition: CRTM_Interpolation.f90:1093
crtm_interpolation::compute_dxi_tl
subroutine compute_dxi_tl(x_TL, xi_TL, dxi_TL)
Definition: CRTM_Interpolation.f90:1202
crtm_interpolation::interp_3d
subroutine, public interp_3d(z, ulp, vlp, wlp, z_int)
Definition: CRTM_Interpolation.f90:264
max
#define max(a, b)
Definition: mosaic_util.h:33
type_kinds
Definition: Type_Kinds.f90:12
crtm_interpolation::order
integer, parameter, public order
Definition: CRTM_Interpolation.f90:70
min
#define min(a, b)
Definition: mosaic_util.h:32
crtm_interpolation::compute_dx_ad
subroutine compute_dx_ad(dx_AD, x_AD)
Definition: CRTM_Interpolation.f90:1173
crtm_interpolation::one
real(fp), parameter one
Definition: CRTM_Interpolation.f90:69
crtm_interpolation::compute_dx_tl
subroutine compute_dx_tl(x_TL, dx_TL)
Definition: CRTM_Interpolation.f90:1164
crtm_interpolation::lpoly_tl
subroutine, public lpoly_tl(x, x_int, p, x_TL, x_int_TL, p_TL)
Definition: CRTM_Interpolation.f90:931
crtm_interpolation::interp_1d
subroutine, public interp_1d(z, ulp, z_int)
Definition: CRTM_Interpolation.f90:230
crtm_interpolation::interp_2d_tl
subroutine, public interp_2d_tl(z, ulp, vlp, z_TL, ulp_TL, vlp_TL, z_int_TL)
Definition: CRTM_Interpolation.f90:400
crtm_interpolation::find_random_index
subroutine find_random_index(x, x_int, i1, i2, out_of_bounds)
Definition: CRTM_Interpolation.f90:772
crtm_interpolation::npoly_pts
integer, parameter npoly_pts
Definition: CRTM_Interpolation.f90:71
crtm_interpolation::interp_2d
subroutine, public interp_2d(z, ulp, vlp, z_int)
Definition: CRTM_Interpolation.f90:246