C******************* MARCH 2000 ADDITION *******************
      SUBROUTINE THERMIK(TREAD,SREAD,PREAD,RALPHA,RBETA)
C    CALCULATES THE THERMAL EXPANSION COEFFICIENT (RALPHA)
C                      [ (degrees C)**-1 ]
C    CALCULATES THE SALINITY CONTRACTION COEFFICIENT (RBETA)
C                       [ pss**-1 ]
C    Input and Output parameters are REAL
C    ALL COMPUTATIONS in DOUBLE PRECISION
C     **REFERENCE**
C    UNESCO(1987); INTERNATIONAL OCEANOGRAPHIC TABLES 	
C    NO.40
C    JA found ERRATA
C    last term in equation (18) at page 104 must be as followed
C     ... (c1 + c2*t) * S**1.5 (NO POWERS for CONSTANTS c1 and c2!!!)
C
C-----------------------------------------------------------------------        
C     THIS SUBROUTINE CALCULATES THE DENSITY OF SEAWATER USING THE NEW          
C     STANDARD EQUATION OF STATE RECOMMENDED BY UNESCO(1981).                   
C     **UNITS**                                                                 
C       IN-SITU TEMPERATURE(T): DEGREES CENTIGRADE                              
C       SALINITY(S): PERMIL                                                     
C       PRESSURE(P): DECIBARS                                                   
C       DENSITY(RHO): KILOGRAMS PER CUBIC METER                                 
C     **REFERENCE**                                                             
C       UNESCO(1981); TENTH REPORT OF THE JOINT PANEL ON OCEANOGRAPHIC          
C       TABLES AND STANDARDS. UNESCO TECHNICAL PAPERS IN MARINE SCI.            
C       NO.36. UNESCO,PARIS.                                                    
C-----------------------------------------------------------------------        
C                                                                               
      DOUBLE PRECISION RW,RSTO,RHO,KW,KSTO,KSTP                                                         
      DOUBLE PRECISION DAW,DA,DBW,DB,DKW,DKSTO,DKSTP,KMON,ALPHA
      DOUBLE PRECISION BETA,DRSTOS,DAS,DBS,DKSTOS,DKSTPS
C                                                                               
      T=DBLE(TREAD)
      S=DBLE(SREAD)
      P=DBLE(PREAD*0.1)                                                               
      RW =     999.842594D0 + 6.793952D-2*T - 9.095290D-3*T**2               
     1        + 1.001685D-4*T**3 - 1.120083D-6*T**4 + 6.536332D-9*T**5          
C                                                                               
      RSTO = RW + (0.824493D0 - 4.0899D-3*T + 7.6438D-5*T**2                    
     1        - 8.2467D-7*T**3 + 5.3875D-9*T**4) * S                            
     2       + (-5.72466D-3 + 1.0227D-4*T - 1.6546D-6*T**2) * S**1.5            
     3       + 4.8314D-4 * S**2                                                 
C                                                                               
      KW =19652.21D0+148.4206D0*T-2.327105D0*T**2 + 1.360477D-2*T**3         
     1        - 5.155288D-5*T**4                                                
C                                                                               
      KSTO =KW+(54.6746D0-0.603459D0*T+1.09987D-2*T**2                     
     1        - 6.1670D-5*T**3) * S                                             
     2       + (7.944D-2 + 1.6483D-2*T - 5.3009D-4*T**2) * S**1.5               
C                                                                               
      KSTP =KSTO+(3.239908D0+1.43713D-3*T+1.16092D-4*T**2                
     1        - 5.77905D-7*T**3) * P                                            
     2       + (2.2838D-3 - 1.0981D-5*T - 1.6078D-6*T**2) * P * S               
     3       + 1.91075D-4 * P * S**1.5                                          
     4       + (8.50935D-5 - 6.12293D-6*T + 5.2787D-8*T**2) * P**2              
     5       + (-9.9348D-7 + 2.0816D-8*T + 9.1697D-10*T**2) * P**2 * S          
C                                                                               
      RHO =    RSTO/(1.0-P/KSTP)                                                
C
C JA ADDITION 03/17/00, THERMAL EXPANSION COEFFICIENT
C                                                                               
      DRW = 6.793952D-2 - 1.819058D-2*T + 3.005055D-4*T**2
     1     - 4.480332D-6*T**3 + 3.268166D-8*T**4
      DRSTO = DRW + (-4.0899D-3 + 1.52876D-4*T - 2.470401D-6*T**2 
     1       +2.155D-8*T**3 )*S + (1.0227D-4 - 3.3092D-6*T)*S**1.5
C
      DAW = 1.43713D-3 + 2.32184D-4*T - 1.733715D-6*T**2 
      DA  = DAW + (-1.0981D-5 - 3.2156D-6*T) * S
      DBW = -6.12293D-6 + 1.05574D-7*T
      DB  = DBW + (2.0816D-8 + 1.83394D-9*T) * S
      DKW = 148.4206D0 - 4.6542D0*T + 4.081431D-2*T**2
     1     -2.0621152D-4*T**3
      DKSTO = DKW +(-0.603459D0 + 2.19974D-2*T - 1.8501D-4*T**2)*S
     1       + (1.6483D-2 - 1.06018D-3*T) * S**1.5
      DKSTP = DKSTO + DA * P + DB * P**2
C
      KMON=1.D0-P/KSTP
      ALPHA= DRSTO/KMON - ((RSTO*P*DKSTP)/KSTP/KMON)/KSTP/KMON
      ALPHA = -1.D0*ALPHA/RHO
      RALPHA=REAL(ALPHA)
C
C SALINITY CONTRACTION COEFFICIENT
C
      DRSTOS = 8.24493D-1 - 4.0899D-3 * T + 7.6438D-5 * T**2  
     1        - 8.2467D-7 * T**3 + 5.3875D-9 * T**4
     2        + 1.5D0 * (-5.72466D-3 + 1.0227D-4 * T
     3        - 1.6546D-6 * T**2) * S**0.5
     4        + 9.6628D-4 * S

      DAS = 2.2838D-3 - 1.0981D-5*T - 1.6078D-6*T**2
     1     +2.866125D-4*S**0.5
      DBS = -9.9348D-7 + 2.0816D-8*T + 9.1697D-10*T**2
      DKSTOS = 54.6746D0 - 0.603459D0*T + 1.09987D-2*T**2
     1        - 6.1670D-5*T**3 + 1.5D0*(7.944D-2
     2        + 1.6483D-2*T - 5.3009D-4*T**2)*S**0.5
      DKSTPS = DKSTOS + DAS * P + DBS * P**2
C
      BETA = DRSTOS/KMON - ((RSTO*P*DKSTPS)/KSTP/KMON)/KSTP/KMON
      BETA = BETA/RHO
      RBETA = REAL(BETA)

      
      RETURN                                                                    
      END                                                                       
C******************* JULY 1997 ADDITION ********************
c     real function atg(s,t,p)
C adiabatic temperature gradient DEG C per DECIBAR
C Bryden, H., 1973, Deep-sea Res., 20, 401-408
C   PRESSURE     P    DECIBARS
C   TEMPERATURE  T    DEG C (IPTS-68)
C   SALINITY     S    (PSS-78)
C   ADIABATIC    ATG   DEG C/DECIBAR
Code after UNESCO #44 (1983), page 40
C***********************************************************
c     ds = s - 35.0
c     atg = (((-2.1687e-16*t + 1.8676e-14)*t-4.6206e-13)*p
c    1      ((2.7759e-12*t-1.1351e-10)*ds+((-5.4481e-14*t
c    2      +8.733e-12)*t-6.7795e-10)*t+1.8741e-8))*p
c    3      +(-4.2393e-8*t+1.8932e-6)*ds
c    4      +((6.6228e-10*t-6.836e-8)*t+8.5258e-6)*t+3.5803e-5
c     return
c     end
C******************** JULY 1997 ADDITION *********************
c***********************************************************
      function deplat(x,xlat)
c pressure in DECIBARS from depth in METERS using
c LEROY's method J.Acoust.Soc.Amer.,1969,46,N 1(part 2),216-226
c  depth     X      in  meters
c  latitude  XLAT   in  degrees
c  pressure  DEPLAT in  decibars ( 100 db =  MPa)
      y= sin(xlat)
      y=y*y
      y=1.02E-1 + (1 + 5.28E-3*y)* 1.00524E-2* x
      deplat=1.E2*(y + 2.475E-8* x* x)
      return
      end
c***********************************************************
      function depkal(x,xlat)
c pressure in DECIBARS from depth in METERS using
c KALASHNIKOV(1985) formula
c
      sqlat=xlat*xlat
      y=(1.00574E-2 - 6.5279E-8*xlat + 2.1889E-8*sqlat - 
     1   1.5709E-10*xlat*sqlat)*x + 2.3486E-8*x*x
      depkal=1.E2*y
      return
      end
      SUBROUTINE DIS(T,S,PREAD,RHO)                                             
C                                                                               
C-----------------------------------------------------------------------        
C     THIS SUBROUTINE CALCULATES THE DENSITY OF SEAWATER USING THE NEW          
C     STANDARD EQUATION OF STATE RECOMMENDED BY UNESCO(1981).                   
C     **UNITS**                                                                 
C       IN-SITU TEMPERATURE(T): DEGREES CENTIGRADE                              
C       SALINITY(S): PERMIL                                                     
C       PRESSURE(P): DECIBARS                                                   
C       DENSITY(RHO): KILOGRAMS PER CUBIC METER                                 
C     **REFERENCE**                                                             
C       UNESCO(1981); TENTH REPORT OF THE JOINT PANEL ON OCEANOGRAPHIC          
C       TABLES AND STANDARDS. UNESCO TECHNICAL PAPERS IN MARINE SCI.            
C       NO.36. UNESCO,PARIS.                                                    
C-----------------------------------------------------------------------        
C                                                                               
      REAL KW,KSTO,KSTP                                                         
C                                                                               
      P=PREAD*0.1                                                               
C                                                                               
      RW =     999.842594 + 6.793952E-2*T - 9.095290E-3*T**2                    
     1        + 1.001685E-4*T**3 - 1.120083E-6*T**4 + 6.536332E-9*T**5          
C                                                                               
      RSTO =   RW + (0.824493 - 4.0899E-3*T + 7.6438E-5*T**2                    
     1        - 8.2467E-7*T**3 + 5.3875E-9*T**4) * S                            
     2       + (-5.72466E-3 + 1.0227E-4*T - 1.6546E-6*T**2) * S**1.5            
     3       + 4.8314E-4 * S**2                                                 
C                                                                               
      KW =     19652.21 + 148.4206*T - 2.327105*T**2 + 1.360477E-2*T**3         
     1        - 5.155288E-5*T**4                                                
C                                                                               
      KSTO =   KW + (54.6746 - 0.603459*T + 1.09987E-2*T**2                     
     1        - 6.1670E-5*T**3) * S                                             
     2       + (7.944E-2 + 1.6483E-2*T - 5.3009E-4*T**2) * S**1.5               
C                                                                               
      KSTP =   KSTO + (3.239908 + 1.43713E-3*T + 1.16092E-4*T**2                
     1        - 5.77905E-7*T**3) * P                                            
     2       + (2.2838E-3 - 1.0981E-5*T - 1.6078E-6*T**2) * P * S               
     3       + 1.91075E-4 * P * S**1.5                                          
     4       + (8.50935E-5 - 6.12293E-6*T + 5.2787E-8*T**2) * P**2              
     5       + (-9.9348E-7 + 2.0816E-8*T + 9.1697E-10*T**2) * P**2 * S          
C                                                                               
      RHO =    RSTO/(1.0-P/KSTP)                                                
C                                                                               
      RETURN                                                                    
      END                                                                       
C                                                                               
      SUBROUTINE PTS(T,S,P,PT)                                                  
C                                                                               
C-----------------------------------------------------------------------        
C     THIS SUBROUTINE CALCULATES THE POTENTIAL TEMPERATURE REFERRED TO          
C     THE SURFACE. IT USES BRYDEN'S(1973) POLYNOMIAL FORMULA.                   
C     **UNITS**                                                                 
C       IN-SITU TEMPERATURE(T): DEGREES CENTIGRADE                              
C       SALINITY(S): PERMIL                                                     
C       PRESSURE(P): DECIBARS                                                   
C       POTENTIAL TEMPERATURE(PT): DEGREES CENTIGRADE                           
C     ** REFERENCE**                                                            
C       H.L.BRYDEN(1973); NEW POLYNOMIALS FOR THERMAL EXPANSION, ADIA-          
C       BATIC TEMPERATURE GRADIENT AND POTENTIAL TEMPERATURE OF SEA             
C       WATER. DEEP-SEA RESEARCH, VOL.20, PP.401-408.                           
C-----------------------------------------------------------------------        
C                                                                               
      PT = T - (0.36504E-4+0.83198E-5*T-0.54065E-7*T**2+0.40274E-9*T**3)        
     1        *P                                                                
     2       - (0.17439E-5-0.29778E-7*T)*P*(S-35.0)                             
     3       + 0.41057E-10*P**2*(S-35.0)                                        
     4       - (0.89309E-8-0.31628E-9*T+0.21987E-11*T**2)*P**2                  
     5       - (-0.16056E-12+0.50484E-14*T)*P**3                                
C                                                                               
      RETURN                                                                    
      END                                                                       
      SUBROUTINE SIGMA0(T,S,P,PD)                                               
C                                                                               
C-----------------------------------------------------------------------        
C     THIS SUBROUTINE CALCULATES THE POTENTIAL DENSITY OF A WATER PARCEL        
C     REFERRED TO THE SURFACE. IT USES BRYDEN'S(1973) FORMULA FOR SUR-          
C     FACE POTENTIAL TEMPERATURE AND UNESCO(1981) EQUATION OF STATE FOR         
C     SEA WATER.                                                                
C     **UNITS**                                                                 
C       IN-SITU TEMPERATURE(T): DEGREES CENTIGRADE                              
C       SALINITY(S): PERMIL                                                     
C       PRESSURE(P): DECIBARS                                                   
C       POTENTIAL DENSITY(PD): GIVEN AS DENSITY ANOMALY (POT.DENSITY            
C         - 1000.0 KG. PER CUBIC METER)                                         
C-----------------------------------------------------------------------        
C                                                                               
      CALL PTS(T,S,P,PT)                                                        
      CALL DIS(PT,S,0.0,RHO)                                                    
C                                                                               
      PD=RHO-1000.0                                                             
C                                                                               
      RETURN                                                                    
      END
