MODULE MODULE_INTP !------------------------------------------------------------------------------! ! Horizontal (2D) bilinear and vertical (1D) linear interpolation ! ! F. VANDENBERGHE, March 2001 !------------------------------------------------------------------------------! CONTAINS !==============================================================================! SUBROUTINE INTPLIN_HORIZ (PVAL, PI, PJ, PFIELD2, KIX, KJX, KCRS) !------------------------------------------------------------------------------! ! ! ROUTINE INTPLIN_HORIZ ! *********************** ! ! PURPOSE: ! ------- ! FAST 2D INTERPOLATION OF MM5 2D FIELD AT OBSERVATIONS MM5 INDECES ! ! METHOD: ! ------- ! PERFORM A BILINEAR INTERPOLATION OF MM5 GRIDDED 2D FIELD AT OBSERVATION ! LOCATION DEFINED IN MM5 INDECES ! ASSUME THAT OBSERVATION LATITUDE AND LONGITUDE HAVE ALREADY BEEN ! CONVERTED INTO MM5 INDECES ! ! BE CAREFULL WITH THE MM5 INDEXATION: ! ------------------------------------ ! ! i (latitude) ! /|\ /\ k (vertical) ! | / ! | / ! -----> j (longitude) ! ! ! INPUT: ! ----- ! CONST: PI, PJ INTERPOLATION LOCATION IN MM5 INDECES ! ACTIVE: PFIELD2 2D MM5 ARRAY DEFINED ON LATITUDE AND LONGITUDE ! CONST: KIX, KJX DIMENSION OF THE MM5 2D ARRAY TO INTERPOLATE ! CONST: KCRS FLAG = 1 WHEN PFIELD2 IS DEFINED AT MM5 CROSS PT ! FLAG = 0 WHEN PFIELD2 IS DEFINED AT MM5 DOT PT ! (NOT USED FOR PLUG COMPATIBILITY WITH BINT) ! OUTPUT: ! ------ ! ACTIVE: PVAL INTERPOLATED VALUE OF PFIELD2 AT (PI,PJ) ! ! ! COMMON: NO ! ------- ! EXTERNAL: NO ! -------- ! REFERENCES: NO ! ---------- ! ! MODIFICATIONS: ! -------------- ! ORIGINAL : 98-07 (F. VANDENBERGHE) ! ADDITIONS : 98-11 Norm DOCTOR (F. VANDENBERGHE) !------------------------------------------------------------------------------! IMPLICIT NONE REAL :: PVAL REAL :: PI, PJ INTEGER :: KIX, KJX REAL, DIMENSION (KIX,KJX) :: PFIELD2 INTEGER :: KCRS INTEGER :: II, IJ REAL :: ZDI, ZDJ, ZDIM, ZDJM !------------------------------------------------------------------------------C ! 1. MODEL HORIZONTAL INDECES ! ============================ II = IFIX (PI) IJ = IFIX (PJ) IF ((1 .LE. II) .AND. (II .LE. KIX-1) .AND. & (1 .LE. IJ) .AND. (IJ .LE. KJX-1)) THEN ! 2. HORIZONTAL INTERPOLATION COEFICIENTS ! ======================================== ! ! I+1,J | | I+1,J+1 ! --+----------+--- ! | | /\ ! | | DIM ! | | \/ ! +PI * + -- ! | | /\ ! | | DI ! | PJ | \/ ! --+----+-----+--- ! I,J ||| I,J+1 ZDI = PI - FLOAT (II) ZDJ = PJ - FLOAT (IJ) ZDIM = 1. - ZDI ZDJM = 1. - ZDJ ! 3. 2D FIELD HORIZONTAL BILINEAIRE INTERPOLATION ! =============================================== PVAL = ZDJM*( ZDIM * PFIELD2 (II, IJ ) & +ZDI * PFIELD2 (II+1,IJ )) & + ZDJ *( ZDIM * PFIELD2 (II, IJ+1) & +ZDI * PFIELD2 (II+1,IJ+1)) ELSE ! 4. ERROR PROCESSING ! ==================== WRITE (0,'(/,A)') ' Warning in sub. intplin_horiz:' WRITE (0,'(A,F7.3,A,I3,A)') ' yi = ',pi,' outside [1 ',kix,'] or' WRITE (0,'(A,F7.3,A,I3,A)') ' xj = ',pj,' outside [1 ',kjx,']' WRITE (0,'(A,/)') ' No interpolation was performed' ! STOP ' in intplin_horiz.F90' PVAL = 0. ENDIF ! 5. INTERPOLATION IS DONE EXIT ! ============================= RETURN END SUBROUTINE INTPLIN_HORIZ SUBROUTINE INTPLIN_VERTI (PVAL, PZK, PROFILE, KKX) !------------------------------------------------------------------------------C IMPLICIT NONE REAL, INTENT (OUT) :: PVAL REAL, INTENT (IN ) :: PZK INTEGER, INTENT (IN) :: KKX REAL, DIMENSION (KKX), INTENT (IN) :: PROFILE ! REAL :: ZF_ABOVE, ZF_BELOW, ZDZ_ABOVE INTEGER IK !------------------------------------------------------------------------------C ! ! 1. PRE-PROCESSING ! ================= ! 1.1 VERTICAL INDEX ! -------------- ! IK = IFIX (PZK) ! ! 1.2 IF POINT IS ABOVE THE MODEL LID, TAKE THE VALUE AT LID AND EXIT ! ------------------------------------------------------------- ! IF (IK .LT. 1) THEN PVAL = PROFILE (1) WRITE (0,'(A,F6.3)') ' Point above model lid z = ',PZK RETURN ENDIF ! 1.3 IF POINT IS ON THE MODEL SURFACE, TAKE SURFACE VALUE AND EXIT ! ------------------------------------------------------------- ! IF (IK .GE. KKX) THEN PVAL = PROFILE (KKX) WRITE (0,'(A,F6.3)') ' Point below model surface z = ',PZK RETURN ENDIF ! ! 1.4 VALUE AT LEVEL JUST ABOVE AND BELOW INTERPOLATION POINT ! ------------------------------------------------------- ! ZF_ABOVE = PROFILE (IK) ZF_BELOW = PROFILE (IK+1) ! ! ! 2. VERTICAL INTERPOLATION ! ========================== ! ! 2.1 DISTANCE BETWEEN POINT AND LEVEL ABOVE ! -------------------------------------- ! ZDZ_ABOVE = PZK - REAL (IK) ! ! ! 2.2 3D FIELD VERTICAL LINEAR INTERPOLATION ! -------------------------------------- ! ! (F (IK+1) - F (IK)) ! VAL = ------------------- * DZ + F (IK) ! IK+1 - IK ! PVAL = (ZF_BELOW - ZF_ABOVE) * ZDZ_ABOVE + ZF_ABOVE ! ! ! ! 3. END ! ======= RETURN END SUBROUTINE INTPLIN_VERTI !==============================================================================! FUNCTION intplin (x,xx,yy) RESULT (val) !------------------------------------------------------------------------------! IMPLICIT NONE REAL, DIMENSION (:) :: xx, yy REAL :: x REAL :: val INTEGER :: n,m,jl !------------------------------------------------------------------------------! n = size (xx) m = size (yy) IF (n .NE. m) THEN WRITE (UNIT = 0, FMT = '(A)' ) ' ERROR arrays must have same size' STOP 'in intplin.F' ENDIF jl = locate (x,xx) IF (jl .LE. 0) THEN val = yy (1) ELSE IF (jl .GE. n) THEN val = yy (n) ELSE val = (xx (jl+1) - x) * yy (jl) + (x - xx (jl)) * yy (jl+1) val = val / (xx (jl+1) - xx (jl)) ENDIF END FUNCTION intplin FUNCTION intplog (x,xx,yy) RESULT (val) !------------------------------------------------------------------------------! IMPLICIT NONE REAL, DIMENSION (:) :: xx, yy REAL :: x REAL :: val INTEGER :: n,m,jl !------------------------------------------------------------------------------! n = size (xx) m = size (yy) IF (n .NE. m) THEN WRITE (UNIT = 0, FMT = '(A)' ) ' ERROR arrays must have same size' STOP 'in intplin.F' ENDIF jl = locate (x,xx) IF (jl .LE. 0) THEN val = yy (1) ELSE IF (jl .GE. n) THEN val = yy (n) ELSE val = log (xx (jl+1) / x) * yy (jl) + log (x / xx (jl)) * yy (jl+1) val = val / log (xx (jl+1) / xx (jl)) ENDIF END FUNCTION intplog FUNCTION locate (x,xx) RESULT (index) !------------------------------------------------------------------------------! IMPLICIT NONE REAL, DIMENSION (:) :: xx REAL :: x INTEGER :: index INTEGER :: n,jl,jm,ju LOGICAL :: ascnd !------------------------------------------------------------------------------! n = size (xx) ascnd = (xx (n) >= xx (1)) ! True if ascending order, false otherwise jl = 0 ! Initialize lower limit ju = n+1 ! Initialize upper limit DO IF (ju-jl <= 1) EXIT ! Repeat until this condition is satisfied jm = (ju+jl) / 2. ! Compute a mid point IF (ascnd .EQV. (x >= xx (jm))) THEN jl = jm ! Replace mid point by lower limit ELSE ju = jm ! Replace mid point by upper limit ENDIF ENDDO IF (x .EQ. xx (1)) THEN ! Set the output, being carefull with endpoints index = 1 ELSE IF (x .EQ. xx (n)) THEN index = n-1 ELSE index = jl ENDIF END FUNCTION LOCATE END MODULE MODULE_INTP