c c------------------------------------------------------------------ c Begin file /home/sting/foster/proclib/tgcmint.f c------------------------------------------------------------------ c SUBROUTINE TGCMINT(GCMIN,GCMHT,RLAT,NLAT,RLON,NLON,HTS,NHTS, + LOGHT,IXFLAG,GCMOUT,IDIM1,IDIM2,IDIM3, + NPRES,PRES,IER,IGRID,SPVAL,IPRNT) C C Ben Foster (NCAR High Altitude Observatory) C foster@ncar.ucar.edu C 303-497-1595 C C PURPOSE: Interpolate tgcm grid data to new lat(s)/lon(s) (linear) C and/or to heights (km) (linear or log) C C ON INPUT: C C GCMIN(73,25,36) = Tgcm parameter data at full grid to interpolate C GCMHT(73,25,36) = Tgcm heights at full grid C RLAT(NLAT) = Latitudes (deg) at which to do linear interp C (-87.5 <= RLAT(I) <= +87.5) C NLAT = Number of latitudes in RLAT (>=1) C RLON(NLON) = Longitudes (deg) at which to do linear interp C (-180 <= RLON(I) <= +180) C NLON = Number of longitudes in RLON (>=1) C HTS(NHTS) = Heights at which to do log or linear interpolation C NHTS = Number of heights in HTS C (>=0, if = 0, see NPRES and PRES) C LOGHT = 0 or 1: if=0 -> do linear interpolation in height C if=1 -> do log interpolation in height C (is usually 1 if doing densities) C IXFLAG = 1 -> 6 depending on index scheme for GCMOUT: C 1 -> GCMOUT(lon,ht,lat) C 2 -> GCMOUT(lon,lat,ht) C 3 -> GCMOUT(lat,lon,ht) C 4 -> GCMOUT(lat,ht,lon) C 5 -> GCMOUT(ht,lat,lon) C 6 -> GCMOUT(ht,lon,lat) C NPRES = Meaningful only if NHTS=0. If NHTS=0 and NPRES>0 C then user wants interpolation in location only. C In this case, IXFLAG must be either 2 or 3, and C pressures given in PRES(NPRES) are returned in C the last dimension of GCMOUT. C (as of 11/88, it is assumed that PRES pressures C correspond to TGCM pressure grid points, i.e., C this routine not yet interpolating to constant C pressure levels. C PRES(NPRES) = Used only when NHTS <=0, and PRES > 0 (see above) C C SPVAL = special value, used when desired height is above maximum C gcmht column at current grid point. Used as follows: C if spval != 0 -> return spval when ht > max(gcmht(i,j)) C if spval = 0 -> extrapolate above max(gcmht(i,j)) to find value C C ON OUPUT: C C GCMOUT = Contains interpolated data as described in IXFLAG C and NPRES above. C IER = 0 if no error has occurred, >0 if an error has C occurred, in which case routine has written C an explanation to FT06. C PARAMETER(IMX=73,KMX=29,JMX=36) DIMENSION GCMIN(IMX,KMX,JMX),GCMHT(IMX,KMX,JMX),RLAT(NLAT), + RLON(NLON),HTS(NHTS),COLINT(KMX),COLHTS(KMX), + GLAT(JMX),GLON(IMX),P(KMX),PRES(NPRES) DIMENSION GCMOUT(IDIM1,IDIM2,IDIM3) DATA NCALLS/0/, DLON/5./, DLAT/5./, DPRES/0.5/ SAVE GLAT,GLON,P C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC TEMP CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC if (iprnt.gt.0) then WRITE(6,"(' ENTER TGCMINT: NLAT NLON=',2I3,' NHTS=',I3,' HTS=', + /(9F8.2))") NLAT,NLON,NHTS,(HTS(K),K=1,NHTS) WRITE(6,"(' IDIM1 2 3=',3I3,' IXFLAG=',I3)") IDIM1,IDIM2,IDIM3, + IXFLAG WRITE(6,"(' LOGHT=',I3,' NPRES=',I3)") LOGHT,NPRES DO 1 I=1,IMX,20 DO 1 J=1,JMX,6 1 WRITE(6,"(' I J=',2I3,' GCMIN=',/(7E10.3))") I,J, + (GCMIN(I,K,J),K=1,KMX) DO 2 I=1,IMX,20 DO 2 J=1,JMX,6 2 WRITE(6,"(' I J=',2I3,' GCMHT=',/(7E10.3))") I,J, + (GCMHT(I,K,J),K=1,KMX) endif CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C NCALLS=NCALLS+1 IER=0 LOGH = LOGHT IF (IXFLAG.LT.1.OR.IXFLAG.GT.6) THEN WRITE(6,"(' >>> TGCMINT: BAD IXFLAG=',I6)") IXFLAG IER=1 RETURN ENDIF IF (LOGH.NE.0.AND.LOGH.NE.1) THEN WRITE(6,"(' >>> TGCMINT: BAD LOGHT=',I6,' WILL DEFAULT TO ', + 'LOGHT=0 (will do linear height interpolations)')") LOGH LOGH=0 ENDIF IF (NLON.LE.0) THEN WRITE(6,"(' >>> TGCMINT: BAD NLON=',I6)") NLON IER=2 RETURN ENDIF IF (NLAT.LE.0) THEN WRITE(6,"(' >>> TGCMINT: BAD NLAT=',I6)") NLAT IER=3 RETURN ENDIF DO 100 I=1,NLON IF (RLON(I).LT.-180..OR.RLON(I).GT.+180.) THEN WRITE(6,"(' >>> TGCMINT: BAD RLON(',I2,')=',F10.2)") I,RLON(I) IER=4 RETURN ENDIF 100 CONTINUE DO 105 J=1,NLAT IF (RLAT(J).LT.-90..OR.RLAT(J).GT.+90.) THEN WRITE(6,"(' >>> TGCMINT: BAD RLAT(',I2,')=',F10.2)") J,RLAT(J) IER=5 RETURN ENDIF 105 CONTINUE IF (NHTS.LE.0) THEN IF (NPRES.LE.0) THEN WRITE(6,"(' >>> TGCMINT: BAD NPRES=',I6)") NPRES IER=6 RETURN ENDIF IF (IXFLAG.LT.2.OR.IXFLAG.GT.3) THEN WRITE(6,"(' >>> TGCMINT: NHTS=',I2,' BAD IXFLAG=',I6)") + NHTS,IXFLAG IER=7 RETURN ENDIF ENDIF C C Set up tgcm grid: C IF (NCALLS.EQ.1) THEN DO 110 I=1,IMX 110 GLON(I) = -185.+I*DLON DO 115 J=1,JMX 115 GLAT(J) = -92.5+J*DLAT DO 120 K=1,KMX 120 P(K) = -7.5+K*DPRES ENDIF C C If IGRID > 0, then user wants to use the tgcm lat/lon grid C -- therefore skip finding indices and skip interp to lat/lon points C C Find longitude indices: DO 200 I=1,NLON IF (IGRID.EQ.1) GOTO 206 DO 205 II=1,IMX-1 IF (RLON(I).GE.GLON(II).AND.RLON(I).LE.GLON(II+1)) THEN LON1 = II LON2 = II+1 GOTO 206 ENDIF 205 CONTINUE WRITE(6,"(' >>> TGCMINT: COULD NOT FIND LON INDEX FOR I=',I3, + ' RLON(I)=',F8.2)") I,RLON(I) 206 CONTINUE C C Find latitude indices: DO 250 J=1,NLAT IF (IGRID.EQ.1) GOTO 266 DO 255 JJ=1,JMX IF (RLAT(J).GE.GLAT(JJ).AND.RLAT(J).LE.GLAT(JJ+1)) THEN LAT1 = JJ LAT2 = JJ+1 GOTO 256 ENDIF 255 CONTINUE WRITE(6,"(' >>> TGCMINT: COULD NOT FIND LAT INDEX FOR J=',I3, + ' RLAT(J)=',F8.2)") J,RLAT(J) 256 CONTINUE C C Do bilinear interp to current RLAT,RLON to get interp press column C of the tgcm input param: DO 260 K=1,KMX FLON = (RLON(I)-GLON(LON1)) / DLON FLAT = (RLAT(J)-GLAT(LAT1)) / DLAT F1 = GCMIN(LON1,K,LAT1) F2 = GCMIN(LON2,K,LAT1) F3 = GCMIN(LON1,K,LAT2) F4 = GCMIN(LON2,K,LAT2) COLINT(K) = FLON*(FLAT*F4+(1.-FLAT)*F2)+ + (1.-FLON)*(FLAT*F3+(1.-FLAT)*F1) 260 CONTINUE C C Do bilinear interp to current RLAT,RLON to get interp press column C of the tgcm input heights: DO 265 K=1,KMX FLON = (RLON(I)-GLON(LON1)) / DLON FLAT = (RLAT(J)-GLAT(LAT1)) / DLAT F1 = GCMHT(LON1,K,LAT1) F2 = GCMHT(LON2,K,LAT1) F3 = GCMHT(LON1,K,LAT2) F4 = GCMHT(LON2,K,LAT2) COLHTS(K) = FLON*(FLAT*F4+(1.-FLAT)*F2)+ + (1.-FLON)*(FLAT*F3+(1.-FLAT)*F1) 265 CONTINUE 266 IF (IGRID.EQ.1) THEN DO 270 K=1,KMX COLHTS(K) = GCMHT(I,K,J) COLINT(K) = GCMIN(I,K,J) 270 CONTINUE ENDIF C C Height loop: IF (NHTS.GT.0) THEN DO 300 K=1,NHTS HTMAX = -1.E36 DO 302 KK=1,KMX 302 HTMAX = AMAX1(HTMAX,COLHTS(KK)) c c If desired height is above max height in column at this grid point c and spval != 0 (i.e, its a real special value), then return spval: c IF (HTS(K).GT.HTMAX) THEN IER=8 if (spval.ne.0.) then temp = spval goto 290 c c If desired height is above max height in column at this grid point c and spval = 0, then spval is a flag that means extrapolate: c (i.e., interpolate from between top two heights) c else ialt1 = kmx-1 ialt2 = kmx goto 306 endif ENDIF c c Bracket desired height: DO 305 KK=1,KMX-1 IF (HTS(K).GE.COLHTS(KK).AND.HTS(K).LE.COLHTS(KK+1)) THEN IALT1 = KK IALT2 = KK+1 GOTO 306 ENDIF 305 CONTINUE c WRITE(6,"(' >>> TGCMINT: COULD NOT FIND INDEX FOR HTS(',I2, c + ')=',F8.2)") K,HTS(K) temp = spval goto 290 306 CONTINUE c c Do the interpolation: IF (LOGH.EQ.1) THEN EXPARG = (ALOG(COLINT(IALT2) / COLINT(IALT1)) / + (COLHTS(IALT2)-COLHTS(IALT1))) * + (HTS(K)-COLHTS(IALT1)) TEMP = COLINT(IALT1)*EXP(EXPARG) ELSE F1 = (HTS(K)-COLHTS(IALT1)) / (COLHTS(IALT2)-COLHTS(IALT1)) TEMP = F1*COLINT(IALT2)+(1.-F1)*COLINT(IALT1) ENDIF C C Transfer to output array: 290 continue IF (IXFLAG.EQ.1) GCMOUT(I,K,J) = TEMP IF (IXFLAG.EQ.2) GCMOUT(I,J,K) = TEMP IF (IXFLAG.EQ.3) GCMOUT(J,I,K) = TEMP IF (IXFLAG.EQ.4) GCMOUT(J,K,I) = TEMP IF (IXFLAG.EQ.5) GCMOUT(K,J,I) = TEMP IF (IXFLAG.EQ.6) GCMOUT(K,I,J) = TEMP 300 CONTINUE C C NHTS.LE.0 -> doing interp in lat/lon only C (PRES(NPRES) in last dim of GCMOUT) C ELSE DO 350 K=1,NPRES DO 355 KK=1,KMX IF (PRES(K).GE.P(KK)-.001.AND.PRES(K).LE.P(KK)+.001) THEN IPRES=KK GOTO 356 ENDIF 355 CONTINUE WRITE(6,"(' >>> TGCMINT: COULD NOT FIND INDEX FOR PRES(', + I2,')=',F8.2)") K,PRES(K) 356 CONTINUE IF (IXFLAG.EQ.2) GCMOUT(I,J,K) = COLINT(IPRES) IF (IXFLAG.EQ.3) GCMOUT(J,I,K) = COLINT(IPRES) 350 CONTINUE ENDIF 250 CONTINUE 200 CONTINUE RETURN END