| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real(kind=r64) | :: | XX | ||||
| real(kind=r64) | :: | YY | ||||
| real(kind=r64), | DIMENSION(:) | :: | X | |||
| real(kind=r64), | DIMENSION(:) | :: | Y | |||
| real(kind=r64), | DIMENSION(:,:) | :: | Z | |||
| integer | :: | NX | ||||
| integer | :: | NY | ||||
| integer | :: | M | ||||
| integer | :: | IEXTX | ||||
| integer | :: | IEXTY |
Nodes of different colours represent the following:
Solid arrows point from a procedure to one which it calls. Dashed arrows point from an interface to procedures which implement that interface. This could include the module procedures in a generic interface or the implementation in a submodule of an interface in a parent module. Where possible, edges connecting nodes are given different colours to make them easier to distinguish in large graphs.
Nodes of different colours represent the following:
Solid arrows point from a procedure to one which it calls. Dashed arrows point from an interface to procedures which implement that interface. This could include the module procedures in a generic interface or the implementation in a submodule of an interface in a parent module. Where possible, edges connecting nodes are given different colours to make them easier to distinguish in large graphs.
FUNCTION DLAG(XX,YY,X,Y,Z,NX,NY,M,IEXTX,IEXTY)
! SUBROUTINE INFORMATION:
! AUTHOR AUTHOR: F.D.HAMMERLING, COMPUTING TECHNOLOGY CENTER, ORNL
! DATE WRITTEN 2010
! MODIFIED Richard Raustad, FSEC
! RE-ENGINEERED na
! PURPOSE OF THIS SUBROUTINE:
!
!
! DLAG IS A DOUBLE-PRECISION, TWO-DIMENSIONAL,
! LAGRANGIAN INTERPOLATION
!
! INPUT VARIABLES:
! XX X-COORDINATE OF THE DESIRED INTERPOLATED POINT
! YY Y-COORDINATE OF THE DESIRED INTERPOLATED POINT
! X SINGLY DIMENSIONED ARRAY OF X-COORDINATES
! Y SINGLY DIMENSIONED ARRAY OF Y-COORDINATES
! Z DOUBLY DIMENSIONED ARRAY OF FUNCTION VALUES,
! I.E. Z(I,J) = F( X(I), Y(J) )
! NX NUMBER OF ELEMENTS IN THE X-ARRAY
! NY NUMBER OF ELEMENTS IN THE Y-ARRAY
! M THE SQUARE ROOT OF THE NUMBER OF POINTS TO BE
! CONSIDERED IN THE INTERPOLATION - NUMBER OF
! POINTS IN EACH DIRECTION
! ID THE FIRST DIMENSION OF THE Z-ARRAY (AT LEAST NX)
!
! OUTPUT VARIABLES:
! IEXTX =1, IF EXTRAPOLATION OCCURED ABOVE THE X-ARRAY
! =0, IF INTERPOLATION OCCURED
! =-1, IF EXTRAPOLATION OCCURED BELOW THE X-ARRAY
! IEXTY SAME FOR THE Y-ARRAY AS IEXTX IS FOR THE X-ARRAY
!
! THIS PROGRAM WAS MODIFIED AUGUST 1984 BY CJ EMERSON TO INSURE
! THAT ITERATIVE CALLS TO DLAG USE THE SAME POINTS FOR INTERPOLATION.
! ISXPT(ISYPT) ARE CHOSEN FROM THE UPPER END OF THE INTERVAL.
! ALSO THE PROGRAM WAS MODIFIED SO THAT EXTRAPOLATION ALWAYS USES
! AT MOST TWO POINTS.
!
!
!
! METHODOLOGY EMPLOYED:
! na
! REFERENCES:
! na
! USE STATEMENTS:
! na
IMPLICIT NONE ! Enforce explicit typing of all variables in this routine
! SUBROUTINE ARGUMENT DEFINITIONS:
REAL(r64) :: DLAG
REAL(r64) :: XX, YY
REAL(r64), DIMENSION(:) :: X, Y
REAL(r64), DIMENSION(:,:) :: Z
!DIMENSION Z(ID,NY),X(NX),Y(NY),XLAG(100)
INTEGER :: NX, NY, IEXTX, IEXTY, M
! SUBROUTINE PARAMETER DEFINITIONS:
! na
! INTERFACE BLOCK SPECIFICATIONS
! na
! DERIVED TYPE DEFINITIONS
! na
! SUBROUTINE LOCAL VARIABLE DECLARATIONS:!
LOGICAL :: QUITX, QUITY
INTEGER :: I, ISXPT, IEXPT, J, ISYPT, IEYPT, K, L, M1
REAL(r64) :: MIDX, MIDY
REAL(r64), ALLOCATABLE, DIMENSION(:) :: XLAG, YLAG
! INITIALIZE
!
L = SIZE(Z(:,1))
L = MAX(L,SIZE(Z(1,:)))+1
ALLOCATE(XLAG(L))
ALLOCATE(YLAG(L))
QUITX=.FALSE.
QUITY=.FALSE.
IEXTX=0
IEXTY=0
L=0
!
! The following code has been upgraded to current Fortran standards
! See Starteam Revision 33, August 17, 2010 for legacy code if comparison is needed
!
! FIND THE RANGE OF INTERPOLATION ALONG X
!
M1=M ! number of points to be interpolated (X direction is first)
IF (M1 .GT. NX) M1=NX ! limit to number of X points if necessary
! loop through X data and find the first x-coordinate less than the interpolated point
! if the interpolation point is less than or greater than the X data then linearly extrapolate
! linear extrapolation uses only 2 points (M1=2)
DO I = 1, NX
IF(XX-X(I) .LT. 0.0D0)THEN
MIDX=I ! found X point just greater than interpolation point
IF (MIDX .EQ. 1) THEN
IEXTX=-1 ! extrapolating at the lower bound of x
IF(M1 .GT. 2)M1=2 ! limit to linear extrapolation
END IF
ISXPT=MIDX-((M1+1)/2) ! calculate starting point in X array
IF (ISXPT .LE. 0)ISXPT=1 ! limit to first element in X array
IEXPT=ISXPT+M1-1 ! calculate ending point in X array
IF (IEXPT .GT. NX) THEN
ISXPT=NX-M1+1 ! if upper X array boundary exceeded, recalculate starting point
IEXPT=NX ! limit ending point to upper boundary of X array
END IF
EXIT
ELSE IF(XX-X(I) .EQ. 0.0D0)THEN ! interpolation point is equal to element in X array
QUITX=.TRUE. ! exact interpolation point found in X array, do not interpolate
EXIT
ELSE IF(I .EQ. NX)THEN ! interpolation point is greater than max X value
IEXTX=1 ! extrapolating at the upper bound of X
IF(M1 .GT. 2)M1=2 ! limit to linear extrapolation
ISXPT=NX-M1+1 ! calculate starting point in X array
IEXPT=NX ! ending point equals upper bound of X array
EXIT
END IF
END DO
M1=M ! number of points to be interpolated (Y direction is second)
IF (M1 .GT. NY) M1=NY ! limit to number of Y points if necessary
DO J=1,NY
IF(YY-Y(J) .LT. 0.0D0)THEN
MIDY=J ! found Y point just greater than interpolation point
IF (MIDY .LE. 1) THEN
IEXTY=-1 ! extrapolating at the lower bound of y
IF(M1 .GT. 2)M1=2 ! limit to linear extrapolation
END IF
ISYPT=MIDY-((M1+1)/2) ! calculate starting point in Y array
IF (ISYPT .LE. 0) ISYPT=1 ! limit to first element in array
IEYPT=ISYPT+M1-1 ! calculate ending point in X array
IF (IEYPT .GT. NY) THEN
ISYPT=NY-M1+1 ! if upper Y array boundary exceeded, recalculate starting point
IEYPT=NY ! limit ending point to upper boundary of Y array
END IF
EXIT
ELSE IF(YY-Y(J) .EQ. 0.0D0)THEN ! interpolation point is equal to element in Y array
QUITY=.TRUE. ! exact interpolation point found in Y array, do not interpolate
EXIT
ELSE IF(J .EQ. NY)THEN ! interpolation point is greater than max Y value
IEXTY=1 ! extrapolating at the upper bound of Y
IF(M1 .GT. 2)M1=2 ! limit to linear extrapolation
ISYPT=NY-M1+1 ! calculate starting point in Y array
IEYPT=NY ! ending point equals upper bound of Y array
EXIT
END IF
END DO
IF (QUITX .AND. QUITY) THEN
DLAG=Z(I,J) ! found exact X and Y point in Z array
ELSE IF (QUITX .AND. .NOT.QUITY) THEN ! only interpolate in Y direction
DO L=ISYPT,IEYPT
XLAG(L)=Z(I,L) ! store X's at each Y (I = midpoint of array from above)
END DO
CALL Interpolate_Lagrange(YY,XLAG,Y,ISYPT,IEYPT,DLAG) ! now interpolate these X's
ELSE IF (.NOT.QUITX .AND. QUITY) THEN ! only interpolate in X direction
CALL Interpolate_Lagrange(XX,Z(:,J),X,ISXPT,IEXPT,DLAG) ! (:,J) interpolate X array at fixed Y (J here)
ELSE ! else interpolate in X and Y directions
DO K=ISYPT,IEYPT
CALL Interpolate_Lagrange(XX,Z(:,K),X,ISXPT,IEXPT,XLAG(K)) ! (:,K) interpolate X array at all Y's (K here)
END DO
CALL Interpolate_Lagrange(YY,XLAG,Y,ISYPT,IEYPT,DLAG) ! final interpolation of X array
END IF
DEALLOCATE(XLAG)
DEALLOCATE(YLAG)
RETURN
END FUNCTION DLAG