c_ ---------------------------------------------------------------------
c_ RCS lines preceded by "c_ "
c_ ---------------------------------------------------------------------
c_
c_ $Source: /home/orr/ocmip/web/OCMIP/phase2/simulations/CFC/boundcond/RCS/cfc_interp.f,v $
c_ $Revision: 1.1 $
c_ $Date: 1998/07/07 15:22:00 $   ;  $State: Exp $
c_ $Author: orr $ ;  $Locker:  $
c_
c_ ---------------------------------------------------------------------
c_ $Log: cfc_interp.f,v $
c_ Revision 1.1  1998/07/07 15:22:00  orr
c_ Initial revision
c_
c_ ---------------------------------------------------------------------
c_
      SUBROUTINE cfc_interp(pnorth, psouth, rlat, imt, jmt, patm) 
c     SUBROUTINE cfc_interp(pnorth, psouth, xintp, rlat, imt, jmt, patm) 

c     ==================================================================
C     ROUTINE to interpolate CFC-11 and CFC-12 DATA for OCMIP-2 model runs

C     Two stations (41S and 45N), each representative of its own hemisphere,
C     except between 10S and 10N WHERE values are interpolated linearly
C     thus there are 3 regions:
C     (1) 90s to 10s, WHERE values take on the value at the station at 41s;
C     (2) 10n to 90n, WHERE values take on the value at the station at 45n; and
C     (3) 10s to 10n, WHERE values are interpolated

c     ---------------------------------------------------------
C     James Orr, LSCE/CEA-CNRS, Saclay, France, 10 June 1998
c     Core algortihm of this ROUTINE by Jean-Claude Dutay, LSCE
c     ---------------------------------------------------------

c     INPUT: 
c     ------
c     Pnorth           REAL   2-member array of CFC-11 and CFC-12 
c                             atmospheric concentrations (mixing ratio) at
c                             the northern station (45N) at time t.
c     Psouth           REAL   2-member array of CFC-11 and CFC-12 
c                             atmospheric concentrations (mixing ratio) at
c                             the southern station (41S) at time t.
c     rlat             REAL   2-D array of latitudes, needed to compute xintp
c     imt           INTEGER   max. dimension in direction i   
c     jmt           INTEGER   max. DIMENSION in direction j

c     OUTPUT:
c     -------
c     Patm             REAL   spatial array for computed atmospheric 
c                             CFC-11 and CFC-12 
c
c     LOCAL:
c     ------
c     xintp            REAL   2-D array of interpolation factors computed on 
c                             1st pass through this routine
c     ientry        INTEGER   Scalar counter to indicate number of calls to
c                             this routine by main PROGRAM
c     ys, yn           REAL   Scaler Latitude limits, within which one does 
c                             a linear interpolation (i.e., 10S to 10N)

c     ==================================================================
      IMPLICIT NONE

      INTEGER nt, ijmax
      PARAMETER (nt=2, ijmax=200*200)

      INTEGER imt, jmt 
      INTEGER ientry
      INTEGER i, j, ij, n     

      REAL ys, yn 
      REAL rlat(imt,jmt), patm(imt,jmt,nt)

      REAL xintp(ijmax) 
      REAL Pnorth(nt), Psouth(nt)

      SAVE ys, yn, ientry, xintp

      DATA ys, yn / -10., 10./
      DATA ientry /0/

      ientry = ientry + 1

c     Test to see if ijmax is large enough
      IF (ijmax .LT. imt*jmt) THEN
          PRINT *," cfc_interp: ERROR -> ijmax must be at least "
     $             ,        imt*jmt
          STOP
      ENDIF

c     IF Block to be executed ONLY on first call to this routine
      IF (ientry .EQ. 1) THEN
          DO j = 1,jmt
            DO i = 1,imt
              ij = (j-1)*imt + i
              IF (rlat(i,j) .GE. yn) THEN
                  xintp(ij) = 1.
              ELSE IF (rlat(i,j) .LE. ys) THEN
                  xintp(ij) = 0.
              ELSE 
                  xintp(ij) =  (rlat(i,j) - ys)/(yn-ys)
              ENDIF 
            END DO
          END DO 
      ENDIF 

c     Block to be executed every pass
      DO n=1,nt
        DO j=1,jmt
          DO i=1,imt
            ij = (j-1)*imt + i
            Patm(i,j,n) =         xintp(ij)  * Pnorth(n) 
     $                   + (1.0 - xintp(ij)) * Psouth(n)
          END DO
        END DO
      END DO
      
      RETURN
      END