      PROGRAM DENSITY                                                                                                               
C     PROGRAM DENSITY(T,       S,       ODENSE,       INPUT,        OUTPUT )    C2                                                  
C    1          TAPE1=T, TAPE2=S, TAPE3=ODENSE, TAPE5=INPUT,  TAPE6=OUTPUT)                                                         
C     ---------------------------------------------------------------------     C2                                                  
C                                                                               C2                                                  
C     JUL 03/2003 - A.IRWIN,F.MAJAESS (REVISED FOR F90 AND LINUX PLATFORMS)     C2                                                  
C     JUN 07/2001 - B. MIVILLE - ORIGINAL VERSION BASED ON OLDER VERSION FROM   C2                                                  
C                                RAMSDEN AND PAL                                C2                                                  
C                                                                               C2                                                  
CDENSITY - CALCULATES OCEAN DENSITY                                     2  1    C1                                                  
C                                                                               C3                                                  
CAUTHOR  - B. MIVILLE                                                           C3                                                  
C                                                                               C3                                                  
CPURPOSE - CALCULATES DENSITY FROM POTENTIAL TEMPERATURE, SALINITY, AND DEPTH   C3                                                  
C                                                                               C3                                                  
C          NOTE: MODEL TEMPERATURE OUTPUT ARE ALREADY POTENTIAL TEMPERATURE.    C3                                                  
C                NO NEED TO CONVERT THEM.                                       C3                                                  
C                                                                               C3                                                  
CINPUT FILES...                                                                 C3                                                  
C                                                                               C3                                                  
C      T = FILE WITH POTENTIAL TEMPERATURE RECORDS (KELVIN ONLY)                C3                                                  
C      S = FILE WITH SALINTY RECORDS (KG/KG)                                    C3                                                  
C                                                                               C3                                                  
COUTPUT FILE...                                                                 C3                                                  
C                                                                               C3                                                  
C      ODENSE = DENSITY (KG/M^3)  (E.G. 1026.1)                                 C3                                                  
C-------------------------------------------------------------------------------                                                    
C                                                                                                                                   
C     * MAXIMUM 2D GRID SIZE IM*JM --> IJM                                                                                          
C                                                                                                                                   
      IMPLICIT REAL (A-H,O-Z),                                                                                                      
     +INTEGER (I-N)                                                                                                                 
      PARAMETER (IM=771,JM=386,IJM=IM*JM,IJMV=IJM*1)                                                                                
                                                                                                                                    
C     * G    = GRAVITY = 9.806 M/S**2                                                                                               
C     * DC   = CONSTANT DENSITY = 1025 KG/M**3                                                                                      
C     *        IT IS AN AVERAGE DENSITY AS USED IN THE OCEAN MODEL                                                                  
C                                                                                                                                   
      PARAMETER (G=9.806E0,DC=1025.0E0)                                                                                             
C                                                                                                                                   
      LOGICAL OK,SPEC                                                                                                               
      REAL T(IJM),S(IJM),ODENSE(IJM)                                                                                                
      REAL UNODEN                                                                                                                   
      INTEGER IBUF,IDAT,MAXX                                                                                                        
      INTEGER TBUF(8)                                                                                                               
C                                                                                                                                   
      COMMON/ICOM/IBUF(8),IDAT(IJMV)                                                                                                
      DATA MAXX/IJMV/                                                                                                               
C---------------------------------------------------------------------                                                              
C                                                                                                                                   
      NFF=4                                                                                                                         
      CALL JCLPNT(NFF,1,2,3,6)                                                                                                      
      REWIND 1                                                                                                                      
      REWIND 2                                                                                                                      
      REWIND 3                                                                                                                      
C                                                                                                                                   
      NR=0                                                                                                                          
C                                                                                                                                   
C     * READ IN POTENTIAL TEMPERATURE (MODEL OUTPUT IS ALREADY IN POTENTIAL T)                                                      
C                                                                                                                                   
 100  CONTINUE                                                                                                                      
      CALL GETFLD2(1,T,NC4TO8("GRID") ,-1,NC4TO8("TEMP"),-1,                                                                        
     +                                         IBUF,MAXX,OK)                                                                        
      IF(.NOT.OK)THEN                                                                                                               
         IF(NR.EQ.0) CALL                          XIT('DENSITY',-1)                                                                
         WRITE(6,6010) NR                                                                                                           
         CALL                                      XIT('DENSITY',0)                                                                 
      ENDIF                                                                                                                         
      IF (NR.EQ.0) CALL PRTLAB(IBUF)                                                                                                
C                                                                                                                                   
C     * SAVE IBUF IN TBUF FOR CHECKING TIMES, LEVEL AND DIMENSIONS                                                                  
C                                                                                                                                   
      DO 110 I=1,8                                                                                                                  
         TBUF(I)=IBUF(I)                                                                                                            
 110  CONTINUE                                                                                                                      
C                                                                                                                                   
C     * READ IN SALT                                                                                                                
C                                                                                                                                   
      CALL GETFLD2(2,S,NC4TO8("GRID"),-1,NC4TO8("SALT"),-1,IBUF,MAXX,OK)                                                            
      IF(.NOT.OK) CALL                             XIT('DENSITY',-2)                                                                
      IF(NR.EQ.0) CALL PRTLAB(IBUF)                                                                                                 
C                                                                                                                                   
C     * MAKE SURE T AND S MATCH IN TIME, SIZE AND DEPTH                                                                             
C                                                                                                                                   
      IF((IBUF(2).NE.TBUF(2)).OR.(IBUF(4).NE.TBUF(4)).OR.                                                                           
     1   (IBUF(5).NE.TBUF(5)).OR.(IBUF(6).NE.TBUF(6)))                                                                              
     2    CALL                                     XIT('DENSITY',-3)                                                                
C                                                                                                                                   
      NWDS=IBUF(5)*IBUF(6)                                                                                                          
C                                                                                                                                   
      IF(NWDS.GT.IJM) CALL                         XIT('DENSITY',-4)                                                                
C                                                                                                                                   
      NR=NR+1                                                                                                                       
C                                                                                                                                   
C     * CALCULATE APPROXIMATE PRESSURE IN DECIBARS FOR THIS LEVEL                                                                   
C                                                                                                                                   
C     * P    = PRESSURE IN DECIBARS = P_PA * 1.0E-4                                                                                 
C     * P_PA = PRESSURE IN PASCALS (M/M^2) = Z * DC * G                                                                             
C     * Z    = DEPTH IN METERS = IBUF(4)/10.0                                                                                       
C                                                                                                                                   
      Z    = FLOAT(IBUF(4))/10.0E0                                                                                                  
      P_PA = Z*DC*G                                                                                                                 
      P    = P_PA * 1.0E-4                                                                                                          
C                                                                                                                                   
C     * CONVERT T,S                                                                                                                 
C                                                                                                                                   
C     * CALCULATE DENSITY                                                                                                           
C                                                                                                                                   
      DO 130 I=1,NWDS                                                                                                               
C                                                                                                                                   
         T(I)=T(I)-273.16E0                                                                                                         
         S(I)=1000.E0*S(I)                                                                                                          
C                                                                                                                                   
         IF(T(I).GT.-200.E0)THEN                                                                                                    
            ODENSE(I)=UNODEN(T(I), S(I), P)                                                                                         
         ELSE                                                                                                                       
            ODENSE(I)=0.E0                                                                                                          
         ENDIF                                                                                                                      
C                                                                                                                                   
 130  CONTINUE                                                                                                                      
C                                                                                                                                   
      IBUF(3)=NC4TO8("DENS")                                                                                                        
      CALL PUTFLD2(3,ODENSE,IBUF,MAXX)                                                                                              
C                                                                                                                                   
      GO TO 100                                                                                                                     
C                                                                                                                                   
C------------------------------------------------------------------------                                                           
 5001 FORMAT(10X,I5)                                                                                                                
 6001 FORMAT('  DENSITY PARAMETER: KIND=',I5)                                                                                       
 6010 FORMAT('0',I6,'  RECORDS PROCESSED')                                                                                          
      END                                                                                                                           
