1 2 !----------------------------------------------------------------------- 3 FUNCTION ZINTERP(XP,YP) 4 !----------------------------------------------------------------------- 5 ! 6 ! --- ISCST3 Version dated 97363 (12/29/97) 7 ! R.W. Brode, PES, Inc. 8 ! Programmed from FUNCTION ZTERR in previous versions 9 ! of ISCST3. Modified to avoid array subscript out of 10 ! bounds at edge of grid. 11 ! 12 ! --- ISCST2 Version: 1.0 Level: 931215 ZTERR 13 ! D. Strimaitis, SRC 14 ! 15 ! PURPOSE: Function computes the elevation (m MSL) at the location 16 ! (xp,yp), by interpolating within field of gridded terrain 17 ! elevations. 18 ! 19 ! 20 ! ARGUMENTS: 21 ! PASSED: xp x-coordinate of the point (m) [r] 22 ! yp y-coordinate of the point (m) [r] 23 ! 24 ! RETURNED: zinterp value interpolated at xp,yp (m MSL) [r] 25 ! 26 ! CALLING ROUTINES: F2INT, TGQA 27 ! 28 ! EXTERNAL ROUTINES: none 29 !----------------------------------------------------------------------- 30 USE DEPVAR 31 IMPLICIT NONE 32 33 SAVE 34 INTEGER :: IXLL , IXLLP1 , IYLL , IYLLP1 35 REAL :: XP , YP , ZINTERP , DI , XPOS , YPOS , TT , UU , ONEMT , & 36 & ONEMU 37 38 ! --- Set inverse of the size of a grid-cell 39 DI = 1./SIZEM ! 0 40 41 ! - ll denotes lower left corner of a grid-cell 42 ! - llm denotes lower left corner of grid-cell (1,1) -- this is the 43 ! lower left corner of the master terrain grid 44 ! 45 ! Full development of the algorithm to obtain value at point xp,yp 46 ! -- array index of lower left corner of cell that contains point 47 ! ixll=(xp-xllm)*di+1 48 ! iyll=(yp-yllm)*di+1 49 ! -- position of lower left value 50 ! xll=xllm+sizem*(ixll-1) 51 ! yll=yllm+sizem*(iyll-1) 52 ! -- fractional position of point within cell wrt lower left corner 53 ! tt=(xp-xll)*di 54 ! uu=(yp-yll)*di 55 ! -- interpolated value 56 ! zi=(1.-tt)*(1.-uu)*zarray(ixll,iyll) 57 ! 1 +tt*(1.-uu)*zarray(ixll+1,iyll) 58 ! 2 +tt*uu*zarray(ixll+1,iyll+1) 59 ! 3 +uu*(1.-tt)*zarray(ixll,iyll+1) 60 61 ! --- Compact representation: 62 XPOS = (XP-XLLM)*DI 63 IXLL = INT(XPOS) + 1 64 IF ( IXLL.GE.NTX ) IXLL = IXLL - 1 65 TT = XPOS - (IXLL-1) 66 ONEMT = 1. - TT 67 IXLLP1 = IXLL + 1 68 YPOS = (YP-YLLM)*DI 69 IYLL = INT(YPOS) + 1 70 IF ( IYLL.GE.NTY ) IYLL = IYLL - 1 71 UU = YPOS - (IYLL-1) 72 ONEMU = 1. - UU 73 IYLLP1 = IYLL + 1 74 !JRA ZINTERP = ONEMT*ONEMU*FLOAT(IZARRAY(IXLL,IYLL)) & 75 !JRA & + TT*ONEMU*FLOAT(IZARRAY(IXLLP1,IYLL)) & 76 !JRA & + TT*UU*FLOAT(IZARRAY(IXLLP1,IYLLP1)) & 77 !JRA & + UU*ONEMT*FLOAT(IZARRAY(IXLL,IYLLP1)) 78 !JRA modified 23/09/2005 to avoid compiler failing to recognise FLOAT 79 ZINTERP = ONEMT*ONEMU*IZARRAY(IXLL,IYLL) & 80 & + TT*ONEMU*IZARRAY(IXLLP1,IYLL) & 81 & + TT*UU*IZARRAY(IXLLP1,IYLLP1) & 82 & + UU*ONEMT*IZARRAY(IXLL,IYLLP1) 83 84 85 CONTINUE 86 END
HyperKWIC - Version 1.00DD executed at 20:00 on 1 Mar 2018 | Personal or Academic or Evaluation User | Free for Non-Commercial, Non-Government Use