      PROGRAM COLOCAT                                                                                                               
C     PROGRAM COLOCAT(INF,      DATAD,      BETA,      COF,      OUTPUT,)       D2                                                  
C    1          TAPE1=INF,TAPE2=DATAD,TAPE3=BETA,TAPE4=COF,TAPE6=OUTPUT)                                                            
C     ------------------------------------------------------------------        D2                                                  
C                                                                               D2                                                  
C     SEP 21/2009 - D.YANG (ONLY READ "GRID" KIND RECORDS FROM "INF" INPUT FILE)D2                                                  
C     JUL 03/2003 - A.IRWIN,F.MAJAESS (REVISED FOR F90 AND LINUX PLATFORMS)                                                         
C     APR 22/2002 - B.MIVILLE (REVISED FOR NEW DATA DESCRIPTION FORMAT)                                                             
C     MAR 05/2002 - F.MAJAESS (REVISED FOR "CHAR" KIND RECORDS)                                                                     
C     FEB 18/2002 - B.MIVILLE (CORRECTION FOR READING BETA WHEN COLOCATING W)                                                       
C     JUN 06/2001 - B.MIVILLE - ORIGINAL VERSION                                                                                    
C                                                                               D2                                                  
CCOLOCAT - RECALCULATE U,V,UISO,VISO,W AT TRACERS LOCATION              3  1    D1                                                  
C                                                                               D3                                                  
CAUTHOR  - B. MIVILLE                                                           D3                                                  
C                                                                               D3                                                  
CPURPOSE - INTERPOLATION OF INPUT FIELD TO THE TRACERS LOCATION                 D3                                                  
C                                                                               D3                                                  
CINPUT FILES...                                                                 D3                                                  
C                                                                               D3                                                  
C     INF   = INPUT GRID FIELD                                                  D3                                                  
C     DATAD = DATA DESCRIPTION FILE                                             D3                                                  
C     BETA  = HEAVISIDE FUNCTION (0 LAND, 1 OCEAN, FRACTION FOR SHAVED          D3                                                  
C             BOTTOM CELL)                                                      D3                                                  
C                                                                               D3                                                  
COUTPUT FILE...                                                                 D3                                                  
C                                                                               D3                                                  
C     COF = NEW COLOCATED FIELD                                                 D3                                                  
C                                                                               D3                                                  
C-------------------------------------------------------------------------------                                                    
C                                                                                                                                   
C     * MAXIMUM 2-D GRID SIZE (I+2)*(J+2)=IM*JM --> IJM                                                                             
C     * FOR MAX GRID SIZE OF I*J PLUS ONE MORE ON EACH SIDE.                                                                        
C                                                                                                                                   
C     * MAXIMUM NUMBER OF LEVELS IS SET IN "KM"                                                                                     
C                                                                                                                                   
      IMPLICIT REAL (A-H,O-Z),                                                                                                      
     +INTEGER (I-N)                                                                                                                 
      PARAMETER (IM=771,JM=386,IJM=IM*JM,IJMV=IJM*1)                                                                                
      PARAMETER (KM=200)                                                                                                            
      REAL BETA(IJM),INF(IJM),COF(IJM),WM(IJM)                                                                                      
      REAL LATF(IJM),LATH(IJM),LONF(IJM),LONH(IJM)                                                                                  
      REAL ZF(KM),ZH(KM)                                                                                                            
      CHARACTER*8 DSRC(10),DSOURCE                                                                                                  
C                                                                                                                                   
      LOGICAL OK                                                                                                                    
C                                                                                                                                   
      INTEGER IBUF,IDAT,MAXX,LEVEL,LEVELP,IBUF3,IBUF3P                                                                              
      INTEGER IBUF2,IBUF7,IBUF8                                                                                                     
C                                                                                                                                   
      COMMON/ICOM/IBUF(8),IDAT(IJMV)                                                                                                
      INTEGER NF,NR                                                                                                                 
C                                                                                                                                   
      INTEGER NWDS,L                                                                                                                
      DATA MAXX/IJMV/                                                                                                               
C---------------------------------------------------------------------                                                              
C                                                                                                                                   
      NF=5                                                                                                                          
      CALL JCLPNT(NF,1,2,3,4,6)                                                                                                     
      REWIND 1                                                                                                                      
      REWIND 2                                                                                                                      
      REWIND 3                                                                                                                      
      REWIND 4                                                                                                                      
C                                                                                                                                   
      NR=0                                                                                                                          
      LEVELP=99999                                                                                                                  
      LEVELB=0                                                                                                                      
      NU=-2                                                                                                                         
C                                                                                                                                   
C     * READ THE DATA DESCRIPTION FILE                                                                                              
C                                                                                                                                   
C     * DATA SOURCE                                                                                                                 
C                                                                                                                                   
      CALL GETFLD2 (NU,DSRC,NC4TO8("CHAR"),-1,NC4TO8("DSRC"),-1,                                                                    
     +                                             IBUF,MAXX,OK)                                                                    
      IF (.NOT.OK) CALL                            XIT('COLOCAT',-1)                                                                
C                                                                                                                                   
C     * VERIFY THAT DSRC IS ONLY OGCM1,2 OR 3 SINCE IT ASSUMES A CERTAIN                                                            
C     * TYPE OF STAGERED GRID                                                                                                       
C                                                                                                                                   
      WRITE(DSOURCE,7000) DSRC(1)                                                                                                   
      IF(DSOURCE.NE.'OGCM1   '.AND.DSOURCE.NE.'OGCM2   '.AND.                                                                       
     1   DSOURCE.NE.'OGCM3   ') THEN                                                                                                
         WRITE(6,6002) DSOURCE                                                                                                      
         CALL                                      XIT('COLOCAT',-2)                                                                
      ENDIF                                                                                                                         
C                                                                                                                                   
C     * LATITUDE OF FULL CELL (TRACER)                                                                                              
C                                                                                                                                   
      CALL GETFLD2 (NU,LATF,NC4TO8("GRID"),-1,NC4TO8(" LAT"),-1,                                                                    
     +                                             IBUF,MAXX,OK)                                                                    
      IF (.NOT.OK) CALL                            XIT('COLOCAT',-3)                                                                
C                                                                                                                                   
      NLON = IBUF(5)                                                                                                                
      NLAT = IBUF(6)                                                                                                                
C                                                                                                                                   
C     * LATITUDE OF HALF CELL (U,V)                                                                                                 
C                                                                                                                                   
      CALL GETFLD2 (NU,LATH,NC4TO8("GRID"),-1,NC4TO8("LATH"),-1,                                                                    
     +                                             IBUF,MAXX,OK)                                                                    
      IF (.NOT.OK) CALL                            XIT('COLOCAT',-4)                                                                
C                                                                                                                                   
      IF((IBUF(5).NE.NLON).OR.(IBUF(6).NE.NLAT))                                                                                    
     1   CALL                                      XIT('COLOCAT',-5)                                                                
C                                                                                                                                   
C     * LONGITUDE OF FULL CELL (TRACER)                                                                                             
C                                                                                                                                   
      CALL GETFLD2 (NU,LONF,NC4TO8("GRID"),-1,NC4TO8(" LON"),-1,                                                                    
     +                                             IBUF,MAXX,OK)                                                                    
      IF (.NOT.OK) CALL                            XIT('COLOCAT',-6)                                                                
C                                                                                                                                   
C     * LONGITUDE OF HALF CELL (U,V)                                                                                                
C                                                                                                                                   
      CALL GETFLD2 (NU,LONH,NC4TO8("GRID"),-1,NC4TO8("LONH"),-1,                                                                    
     +                                             IBUF,MAXX,OK)                                                                    
      IF (.NOT.OK) CALL                            XIT('COLOCAT',-7)                                                                
C                                                                                                                                   
      IF((IBUF(5).NE.NLON).OR.(IBUF(6).NE.NLAT))                                                                                    
     1   CALL                                      XIT('COLOCAT',-8)                                                                
C                                                                                                                                   
C     * DEPTH LEVEL OF TRACER                                                                                                       
C                                                                                                                                   
      CALL GETFLD2 (NU,ZF,NC4TO8("GRID"),-1,NC4TO8("   Z"),-1,                                                                      
     +                                           IBUF,MAXX,OK)                                                                      
      IF (.NOT.OK) CALL                            XIT('COLOCAT',-9)                                                                
C                                                                                                                                   
      NLEV = IBUF(5)                                                                                                                
C                                                                                                                                   
C     * DEPTH LEVEL OF OMEGA (W)                                                                                                    
C                                                                                                                                   
      CALL GETFLD2 (NU,ZH,NC4TO8("GRID"),-1,NC4TO8("  ZH"),-1,                                                                      
     +                                           IBUF,MAXX,OK)                                                                      
      IF (.NOT.OK) CALL                            XIT('COLOCAT',-10)                                                               
C                                                                                                                                   
      IF(IBUF(5).NE.NLEV)CALL                      XIT('COLOCAT',-11)                                                               
C                                                                                                                                   
C     * DO LOOP THROUGH THE INPUT FILE                                                                                              
C                                                                                                                                   
      K=1                                                                                                                           
C                                                                                                                                   
 200  CONTINUE                                                                                                                      
C                                                                                                                                   
C     * READ INPUT GRID FIELD                                                                                                       
C                                                                                                                                   
      CALL GETFLD2 (1,INF,NC4TO8("GRID"),-1,-1,-1,IBUF,MAXX,OK)                                                                     
      IF (.NOT.OK)THEN                                                                                                              
         IF (NR.EQ.0) CALL                         XIT('COLOCAT',-12)                                                               
         CALL                                      XIT('COLOCAT',0)                                                                 
      ENDIF                                                                                                                         
C                                                                                                                                   
      IBUF3=IBUF(3)                                                                                                                 
      IBUF8=IBUF(8)                                                                                                                 
C                                                                                                                                   
C     * MAKE SURE INFIELD HAS SAME DIMENSION AS OTHER VARIABLES                                                                     
C                                                                                                                                   
      IF((IBUF(5).NE.NLON).OR.(IBUF(6).NE.NLAT))                                                                                    
     1  CALL                                       XIT('COLOCAT',-13)                                                               
C                                                                                                                                   
C                                                                                                                                   
C     * SAVE SOME IBUF SO WE CAN WRITE TO FILE USING THE SAME VALUES                                                                
C                                                                                                                                   
      IBUF2=IBUF(2)                                                                                                                 
      IBUF3=IBUF(3)                                                                                                                 
      IBUF7=IBUF(7)                                                                                                                 
      LEVEL=IBUF(4)                                                                                                                 
C                                                                                                                                   
      IF(NR.EQ.0)THEN                                                                                                               
C                                                                                                                                   
C        * ORIGINAL INPUT FIELD LABEL                                                                                               
C                                                                                                                                   
         IBUF3P=IBUF3                                                                                                               
         NLON=IBUF(5)                                                                                                               
         NLAT=IBUF(6)                                                                                                               
         NWDS=NLON*NLAT                                                                                                             
C                                                                                                                                   
C        * CHECK IF FIELD HAS REPEATED VALUE AT REPEATED LONGITUDE                                                                  
C                                                                                                                                   
         IF(INF(1).NE.INF(NLON))THEN                                                                                                
            WRITE(6,6020)'FIELD DOES NOT REPEAT AT REPEATED LONGITUDE'                                                              
            WRITE(6,6010)'FIELD(1): ',INF(1),' FIELD(NLON): ',                                                                      
     1           INF(NLON)                                                                                                          
            CALL                                   XIT('COLOCAT',-14)                                                               
         ENDIF                                                                                                                      
C                                                                                                                                   
C        * CHECK IF NLON OF FIELD IS ODD MEANING IT HAS REPEATED LONGITUDE                                                          
C                                                                                                                                   
         IODD=MOD(NLON,2)                                                                                                           
         IF(IODD.EQ.0)THEN                                                                                                          
            WRITE(6,6020)'LAST LONGITUDE IS NOT A REPEAT OF FIRST ONE'                                                              
            WRITE(6,6010)'NLON: ',NLON                                                                                              
            CALL                                   XIT('COLOCAT',-15)                                                               
         ENDIF                                                                                                                      
      ELSE                                                                                                                          
C                                                                                                                                   
C        * CHECK IF FIELD CHANGED LABEL                                                                                             
C                                                                                                                                   
         IF(IBUF3.NE.IBUF3P) THEN                                                                                                   
            WRITE(6,6020)' * FIELD NAME HAS CHANGED *'                                                                              
            WRITE(6,6010)'RECORD NUMBER: ',NR+1                                                                                     
            CALL                                   XIT('COLOCAT',-16)                                                               
         ENDIF                                                                                                                      
C                                                                                                                                   
C        * CHECK IF FIELD HAS SAME DIMENSION THROUGHT THE INPUT FILE                                                                
C                                                                                                                                   
         IF((IBUF(5).NE.NLON).OR.(IBUF(6).NE.NLAT))THEN                                                                             
            WRITE(6,6020) ' * DIMENSION OF FIELD HAS CHANGED *'                                                                     
            WRITE(6,6010) 'RECORD NUMBER: ',NR+1                                                                                    
            CALL                                   XIT('COLOCAT',-17)                                                               
         ENDIF                                                                                                                      
C                                                                                                                                   
C        * CHECK IF LEVEL IS LESS DEEP THAN PREVIOUS LEVEL                                                                          
C        * WE WILL THEN REWIND THE BETA FILE                                                                                        
C                                                                                                                                   
         IF(LEVEL.LT.LEVELP) THEN                                                                                                   
            K=1                                                                                                                     
            REWIND 3                                                                                                                
            LEVELB=0                                                                                                                
         ENDIF                                                                                                                      
C                                                                                                                                   
      ENDIF                                                                                                                         
C                                                                                                                                   
C     * INPUT BETA                                                                                                                  
C                                                                                                                                   
C     * NOTE: W IS LOCATED AT THE BOTTOM OF A VERTICAL GRID BOX WHILE                                                               
C     *       BETA IS LOCATED AT THE CENTER OF A VERTICAL GRID BOX.                                                                 
C                                                                                                                                   
      IF(IBUF3.EQ.NC4TO8("   W"))THEN                                                                                               
         CALL GETFLD2 (3,BETA,NC4TO8("GRID"),-1,NC4TO8("BETA"),-1,                                                                  
     +                                               IBUF,MAXX,OK)                                                                  
         IF(IBUF(4).GT.LEVEL) THEN                                                                                                  
            WRITE(6,6020) ' * BETA LEVEL SHOULD BE LESS THAN FIELD *'                                                               
            WRITE(6,6010) ' * FIELD LEVEL: ',LEVEL                                                                                  
            WRITE(6,6010) ' * BETA  LEVEL: ',IBUF(4)                                                                                
         ENDIF                                                                                                                      
         IF (.NOT.OK) CALL                         XIT('COLOCAT',-18)                                                               
         LEVELB=IBUF(4)                                                                                                             
      ELSE                                                                                                                          
         IF (LEVEL.NE.LEVELP)THEN                                                                                                   
            CALL GETFLD2 (3,BETA,NC4TO8("GRID"),-1,NC4TO8("BETA"),                                                                  
     +                                         LEVEL,IBUF,MAXX,OK)                                                                  
            IF (.NOT.OK)THEN                                                                                                        
               WRITE(6,6020) ' * COULD BE RELATED TO NON MATCH LEVELS *'                                                            
               WRITE(6,6010) ' * FIELD LEVEL: ',LEVEL                                                                               
               WRITE(6,6010) ' * BETA  LEVEL: ',IBUF(4)                                                                             
C              * ERROR MESSAGE IF NO FIELD WAS READ                                                                                 
               CALL                                XIT('COLOCAT',-19)                                                               
            ENDIF                                                                                                                   
         ENDIF                                                                                                                      
      ENDIF                                                                                                                         
C                                                                                                                                   
C     * MAKE SURE IT HAS SAME DIMENSION AS ALL OTHER VARIABLES                                                                      
C                                                                                                                                   
      IF(IBUF(5).NE.NLON)CALL                      XIT('COLOCAT',-20)                                                               
      IF(IBUF(6).NE.NLAT)CALL                      XIT('COLOCAT',-21)                                                               
C                                                                                                                                   
C     * MAKE SURE THAT SHAVED CELL BETA ARE SET TO 1                                                                                
C     * SET W(1/2) = 0 (AT K=1)                                                                                                     
C                                                                                                                                   
      DO 230 L=1,NWDS                                                                                                               
         IF(BETA(L).GT.0.00001E0) THEN                                                                                              
            BETA(L)=1.E0                                                                                                            
         ELSE                                                                                                                       
            BETA(L)=0.E0                                                                                                            
         ENDIF                                                                                                                      
         IF (K.EQ.1) WM(L)=0.E0                                                                                                     
 230  CONTINUE                                                                                                                      
C                                                                                                                                   
C     * KEEP PREVIOUS LEVEL                                                                                                         
C                                                                                                                                   
      LEVELP=LEVEL                                                                                                                  
C                                                                                                                                   
C     * DO COLOCATION BASED ON NAME (LABEL) OF VARIABLES TO BE COLOCATED                                                            
C                                                                                                                                   
      IF(IBUF3.EQ.NC4TO8("   U").OR.IBUF3.EQ.NC4TO8("   V")) THEN                                                                   
         CALL COUV(INF,BETA,NLON,NLAT,LATF,LONF,LATH,LONH,COF)                                                                      
      ELSEIF(IBUF3.EQ.NC4TO8("UISO")) THEN                                                                                          
         CALL COUISO(INF,BETA,NLON,NLAT,LATF,LONF,LATH,LONH,COF)                                                                    
      ELSEIF(IBUF3.EQ.NC4TO8("VISO")) THEN                                                                                          
         CALL COVISO(INF,BETA,NLON,NLAT,LATF,LONF,LATH,LONH,COF)                                                                    
      ELSEIF(IBUF3.EQ.NC4TO8("   W")) THEN                                                                                          
C                                                                                                                                   
C     * COLOCATION OF OMEGA (W)                                                                                                     
C                                                                                                                                   
C     *                                                                                                                             
C     *     -------------W----------- Z(K-1/2)=ZH(K-1)                                                                              
C     *                                                                                                                             
C     *                                                                                                                             
C     *     -------------T----------- Z(K)=ZF(K)                                                                                    
C     *                                                                                                                             
C     *                                                                                                                             
C     *     -------------W----------- Z(K+1/2)=ZH(K)                                                                                
C     *                                                                                                                             
C     *     WE DO AN INTERPOLATION USING                                                                                            
C     *                                                                                                                             
C     *     COW=(W(K-1/2)* DZ(K-1/4) + W(K+1/2)* DZ(K+1/4))                                                                         
C     *          ------------------------------------------                                                                         
C     *                         (DZ(K))                                                                                             
C     *                                                                                                                             
C     *     DZ(K+1/4) = Z(K+1/2) - Z(K)                                                                                             
C     *     DZ(K-1/4) = Z(K) - Z(K-1/2)                                                                                             
C     *                                                                                                                             
C     *     Z(1/2) = 0   SURFACE                                                                                                    
C     *                                                                                                                             
         IF(K.EQ.1) THEN                                                                                                            
            DZM = ZF(K)                                                                                                             
            DZ = ZH(K)                                                                                                              
         ELSE                                                                                                                       
            DZM = ZF(K) - ZH(K-1)                                                                                                   
            DZ = ZH(K) - ZH(K-1)                                                                                                    
         ENDIF                                                                                                                      
         DZP = ZH(K) - ZF(K)                                                                                                        
         DO 250 L=1,NWDS                                                                                                            
            COF(L) = BETA(L)*(INF(L)*DZP + WM(L)*DZM)/DZ                                                                            
            WM(L)=INF(L)                                                                                                            
 250     CONTINUE                                                                                                                   
      ELSE                                                                                                                          
         WRITE(6,6020) 'FIELD LABEL IS NOT VALID '                                                                                  
         CALL                                      XIT('COLOCAT',-22)                                                               
      ENDIF                                                                                                                         
C                                                                                                                                   
C     * WRITE COLOCATED FIELD TO FILE                                                                                               
C                                                                                                                                   
      IF(IBUF3.EQ.NC4TO8("   W")) THEN                                                                                              
         CALL SETLAB(IBUF,NC4TO8("GRID"),IBUF2,IBUF3,LEVELB,NLON,NLAT,                                                              
     1               IBUF7,IBUF8)                                                                                                   
      ELSE                                                                                                                          
         CALL SETLAB(IBUF,NC4TO8("GRID"),IBUF2,IBUF3,LEVEL,NLON,NLAT,                                                               
     1               IBUF7,IBUF8)                                                                                                   
      ENDIF                                                                                                                         
C                                                                                                                                   
      CALL PUTFLD2(4,COF,IBUF,MAXX)                                                                                                 
C                                                                                                                                   
      NR=NR+1                                                                                                                       
C                                                                                                                                   
      IF(K.GT.NLEV) CALL                           XIT('COLOCAT',-23)                                                               
      K=K+1                                                                                                                         
      GOTO 200                                                                                                                      
C                                                                                                                                   
C---------------------------------------------------------------------                                                              
 6002 FORMAT('DATA SOURCE WRONG: ',A8)                                                                                              
 6010 FORMAT(A,I8)                                                                                                                  
 6020 FORMAT(A)                                                                                                                     
 7000 FORMAT(A8)                                                                                                                    
      END                                                                                                                           
