dataintp.f

Go to the documentation of this file.

      SUBROUTINE dataintp(LL,LAT,LON,F1,DT1,F2,DT2,IPT,NBAND,RNG,DEF
     .  ,INTORDER,DUMMY,FOUT, func, int_bad,IR,JC)
C**********************************************************************
C  Data interpolation (spatial and temporal) subroutine.       
C  Spatial interpolation is weighted by distance and temporal interpolation
C  is a simple linear approximation.       
C  [Only valid surrounding data will be used.  If no valid surrounding value 
C  found, the user specified values will be assigned as output values.]
C------------------------------------------------------------
C  Input parameters:
C  -----------------
C  LL        R*4(2)   Actual pixel location, latitude (-90.0 to 90.0) and
C                     longitude (-180. to 180.)
C  LAT       R*4(x)   Surrounding latitudes (x=IPT=number of data points 
C                     surrounding LL), -90.0 to 90.0.
C  LON       R*4(x)   Surrounding longitudes (x=IPT=number of data points
C                     surrounding LL), -180.0 to 180.0.
C  F1        R*4()    Time step 1 surrounding data values.  For the case of 
C                     each data point has only one value (such as P,SST,...),
C                     the callig program declares F1 as a one dimensional array
C                     , F1(IPT).
C                     For multiple values data (e.g. ozone optical thicknesses),
C                     F1 will be declared as a two dimensional array,
C                     F1(IPT,NBAND). 
C  DT1       R*4      Time difference bt. actual time and time step 1. (DT1>=0.)
C                     [Note: actual time is a time bt. time steps 1 & 2].
C  F2        R*4()    Same as F1 except for time step 2.
C  DT2       R*4      Time difference bt. actual time and time step 2. (DT2>=0.)
C  IPT       I*4      Number of data points surrounding position LL at any
C                     single time step.  (1<=IPT<=25)
C  NBAND     I*4      Band number of each data point.  If each data point has
C                     only one corresponding value, set NBAND=1. (1<=NBAND<=15)
C  RNG       R*4(2)   Data thresholds.  Any surrounding data outside the range
C                     (i.e. data less than or data greater than specified
C                     range) will be excluded for data interpolation.
C  DEF       R*4()    User specified default value(s) which will only be used
C                     if no data is valid for interpolation.  If NBAND=1 the
C                     calling program declares DEF as a scalar.  Otherwise,
C                     DEF has the size of NBAND, DEF(NBAND).
C  INTORDER  I*4      Interpolation order.  Valid sequences are:
C                     112: space then time;    
C                     121: time then space;       
C                     110: spatial interpolation only;
C                     101: temporal interpolation only.
C                     [Note: If all data are good, results from 112 and 121 
C                     are identical] 
C  DUMMY     R*4()    A dummy array which has the size as that of F1.
C  IR,JC     I*4      Sizes, Row and column, of output parameter FOUT.
C                     1)If you expect the output as a single value, set
C                       IR=1, JC=1 and FOUT as a scalar;
C                     2)Output data contain NBAND values, set
C                       IR=1, JC=NBAND and FOUT(NBAND);
C                     3)Output data contain IPT values (e.g. INTORDER=101), set
C                       IR=IPT, JC=1 and FOUT(IPT);  
C                     4)Results contain IPT & NBAND values (e.g. INTORDER=101), 
C                       IR=IPT, JC=NBAND and FOUT(IPT,NBAND).
C  Output parameters:
C  -----------------
C  FOUT      R*4()    Interpolated results for actual location and/or time from
C                     surrounding points. (Reference: parameters IR,JC)
C  func      R*4()    Uncertainty estimate - for 112 intorder, the difference 
C                     between the 2 time values
C  int_bad   I*4      flag to indicate that the space interpolation used
C                     the default or one or both of the inputs to the time 
C                     interpolation was a default - cause for ATMWARN
C                     0 if good, 1 if bad
C------------------------------------------------------------
C
C  Created by Eueng-nan Yeh, 11/17/92, GSC
C  W. Robinson, SAIC 6 Dec 2013  add the int_bad argument
C  W. Robinson, SAIC 4 Aug 2016  add an uncertainty estimate and pass back
C
C**********************************************************************
C%%
 
      INTEGER*4 MBAND
      parameter(mband=15)
      INTEGER*4 NBAND,IPT,INTORDER,IR,JC, int_bad, int_bad1, int_bad2
      real*4    ll(1),lat(1),lon(1),f1(ipt,1),f2(ipt,1),rng(*),def(1)
     .         ,fout(1),fout1(mband),fout2(mband),dummy(ipt,1)
     .         ,mimx(2)
      real*4    func(1)
      real*8    dt1,dt2
 
      int_bad = 0
      IF(rng(2).GT.rng(1))THEN
        mimx(2)=rng(2)       ! maximum      
        mimx(1)=rng(1)       ! minimum
      ELSE
        mimx(1)=rng(2)
        mimx(2)=rng(1)
      END IF 
   
      IF(intorder .EQ. 112)THEN            ! 112
        CALL spaceint(ll,lat,lon,f1,ipt,nband,mimx,def,fout1, int_bad1 )
        CALL spaceint(ll,lat,lon,f2,ipt,nband,mimx,def,fout2, int_bad2 )
        CALL timeint(fout1,fout2,dt1,dt2,1,nband,mimx,def,fout,func,ir,jc)
        IF( ( int_bad1 .EQ. 1 ) .OR. ( int_bad2 .EQ. 1 ) ) int_bad = 1
      ELSE
        ntot = nband * ir
        DO 10 i = 1, ntot
          func(i) = 0.
   10     CONTINUE
        IF(intorder.EQ.110)THEN            ! 110
          CALL spaceint(ll,lat,lon,f1,ipt,nband,mimx,def,fout, int_bad )
        ELSE
          if(intorder.EQ.101)then          ! 101
            CALL timeint(f1,f2,dt1,dt2,ipt,nband,mimx,def,fout,func,ir,jc)
          else
            IF(intorder.EQ.121)THEN        ! 121
              CALL timeint(f1,f2,dt1,dt2,ipt,nband,mimx,def,dummy,func
     .          ,ipt,nband)
              CALL spaceint(ll,lat,lon,dummy,ipt,nband,mimx,def,fout, int_bad )
            END IF
          end if
        END IF
      END IF
      RETURN
      END
 
      SUBROUTINE timeint(F1,F2,DT1,DT2,IPT,NBAND,MIMX,DEF,FOUT,func
     .   ,IR,JC)
 
      INTEGER*4 NBAND,IPT,IR,JC,N,I,K,IJ
      REAL*4    F1(IPT,1),F2(IPT,1),MIMX(*)
     .         ,def(1),fout(1),func(1)
      real*8    dt1,dt2,t1,t2,w1,w2
C:  FOUT(1)=FOUT(# data point, # of band)   
      t1=abs(dt1)
      t2=abs(dt2)
      IF(t1.EQ.0. .AND. t2.EQ.0.) t2=1
C: Interpolating data between T1 and T2
      w1=t2/(t1+t2)
      w2=t1/(t1+t2)
      IF (w1.GE.1.) w2=0.
      IF (w2.GE.1.) w1=0.
      DO 200 n=1,nband
        ij=(n-1)*ir
        DO 201 i=1,ipt
          k=ij+i
           IF(f1(i,n).GE.mimx(1).AND.f1(i,n).LE.mimx(2))THEN  ! F1 is OK
            if(f2(i,n).GE.mimx(1).AND.f2(i,n).LE.mimx(2))then ! F1 OK & F2 OK
              fout(k)=w1*f1(i,n)+w2*f2(i,n)
              func(k)=abs(f1(i,n)-f2(i,n))
            else
              fout(k)=f1(i,n)           ! F1 OK & F2 NG
              func(k)=0
            end if
          ELSE                     ! F1 is NG
            if(f2(i,n).GE.mimx(1).AND.f2(i,n).LE.mimx(2))then
              fout(k)=f2(i,n)           ! F1 NG & F2 OK
              func(k)=0
            else
              fout(k)=def(n)            ! F1 NG & F2 NG
              func(k)=0
            end if
          END IF
201     CONTINUE
200   CONTINUE    
 
9000  CONTINUE
      RETURN
      END
c
c
c
      SUBROUTINE spaceint(LL,LAT,LON,F,IPT,NBAND,MIMX,DEF,FOUT, int_bad )
C  Rectangular bi-linear interpolation (IPT==4).
C  by Eueng-nan Yeh           GSC/SAIC                 11/17/95
c
      INTEGER*4 NBAND,IPT,I,N,NNG, int_bad
      INTEGER*4 NCC,NC(4),N1(3,4),N2(5,6),N3(2,4)
      real*4    ll(*),lat(*),lon(*),f(ipt,1),mimx(*),def(1)
     $         ,fout(1),llft(2),urht(2),dx,dy,dd,xp,yp,a(4),g(4)
      DATA nc/1,10,100,1000/
      DATA n1/1,2,4,10,1,3,100,2,4,1000,1,3/
      DATA n2/11,1,2,4,3,101,1,3,4,2,1001,1,4,2,3,110,2,3,1,4
     $     ,1010,2,4,3,1,1100,3,4,2,1/
      DATA n3/111,4,1110,1,1101,2,1011,3/
      llft(1)=lat(1)
      llft(2)=lon(1)
      urht(1)=lat(3)
      urht(2)=lon(3)
C: Rectangular bi-linear interpolation.
      xp=ll(2)-llft(2)
      IF(abs(xp) .gt. 180.) xp=sign(360-abs(xp), xp)
      yp=ll(1)-llft(1)
      dx=urht(2)-llft(2)
      IF(abs(dx) .gt. 180.) dx=sign(360-abs(dx), dx)
      dy=urht(1)-llft(1)
      dd=dx*dy
      DO 110 n=1,nband
        ncc=0
        nng=0
        int_bad = 0
        DO 115 i=1,ipt   !IPT==4
          g(i)=f(i,n)    !assign F to local variable G
          IF(g(i).LT.mimx(1) .OR. g(i).GT.mimx(2))THEN  !F is ng
            ncc=ncc+nc(i)
            nng=nng+1
          ENDIF
115     CONTINUE
        IF(nng .eq. 0)GOTO 9000
 
        IF(nng .eq. 2)THEN    !Two bad points
          DO i=1,6   ! WDR replace IPT with 6 to correct
            IF(ncc .eq. n2(1,i))THEN
              g(n2(2,i))=g(n2(4,i))
              g(n2(3,i))=g(n2(5,i))
              GOTO 9000
            ENDIF
          ENDDO
        ELSE
          IF(nng .eq.1)THEN   !One bad point
            DO i=1,ipt
              IF(ncc .eq. n1(1,i))THEN
                g(i)=(g(n1(2,i))+g(n1(3,i)))/2.
                GOTO 9000
              ENDIF
            ENDDO
          ENDIF
        ENDIF
        IF(nng .eq. 4)THEN    !Four bad points
          fout(n)=def(n)
          int_bad = 1
          GOTO 110
        ELSE
          IF(nng .eq. 3) THEN !Three bad points
            DO i=1,ipt
              IF(ncc .eq. n3(1,i))THEN
                fout(n)=g(n3(2,i))
                GOTO 110
              ENDIF
            ENDDO
          ENDIF
        ENDIF
9000    CONTINUE
        a(1)=g(1)
        IF (abs(dx) .le. 1.e-9) THEN
           a(2) = 0
        ELSE
           a(2)=(g(4)-a(1))/dx
        ENDIF 
        IF (abs(dy) .le. 1.e-9) THEN
           a(3) = 0
        ELSE
           a(3)=(g(2)-a(1))/dy
        ENDIF
        IF (abs(dd) .le. 1.e-9) THEN
            a(4) = 0
        ELSE    
            a(4)=(a(1)-g(2)+g(3)-g(4))/dd
        ENDIF
        fout(n)=a(1)+a(2)*xp+yp*(a(3)+a(4)*xp)
110   CONTINUE
      RETURN
      END