REAL FUNCTION FPRESS(DEPTH,LAT,DLAST,PLAST,N)
*********************************************************************
c     Function fpress computes pressure from depth using
c     the algorithm outlined in Processing of Oceanographic
c     Station Data (JPOTS)
c
c     DEPTH - depth 
c     LAT   - latitude
c     DLAST  - last recorded depth
c     PLAST - last recorded pressure
c     N - number of iterations of trapezoidal integration
c
c     Check Value: DEPTH=5000, LAT=15, DLAST=4500, PLAST=4575.04
c                  P = 5089.22
c
*********************************************************************

      REAL LAT
      FPRESS=0.0
      IT=0
      X=SIN(LAT/57.29578)
      X=X*X
c RHOLAST - density of seawater at 35 PSU and 0 degC
c           for last recorded pressure
      RHOLAST = FDENSE(0.,35.,PLAST)+1000.
c GRLAST - Gravity at last recorded pressure
      GRLAST=9.780318*(1.0+(5.2788E-3+2.36E-5*X)*X)+1.092E-6*PLAST
c P=Pressure estimate calculated using last recorded pressure
      P=PLAST + ((RHOLAST*GRLAST)*(DEPTH-DLAST))/1.E4

      DO 50 NN=1,N

      IF ( NN .EQ. 1 ) THEN

      GR0=9.780318*(1.0+(5.2788E-3+2.36E-5*X)*X)
      RHO0 = FDENSE(0.,35.,0.) + 1000.

* GR=Gravity variation with Latitude: ANON (1970) Bulletin Geodesique
      GR=9.780318*(1.0+(5.2788E-3+2.36E-5*X)*X)+1.092E-6*P
* RHO=Density of seawater at 35 PSU and 0 degC, pressure P
      RHO = FDENSE(0.,35.,P) + 1000.
      FPRESS= .5 * (DEPTH/1.E4)*((RHO*GR)+(RHO0*GR0))
      IT=1

      ELSE

      TNM=IT
      DEL=P/TNM
      Y=0.5*DEL
      SUM=0.

      DO 100 J=1,IT
      GR=9.780318*(1.0+(5.2788E-3+2.36E-5*X)*X)+1.092E-6*Y
      RHO = FDENSE(0.,35.,Y) + 1000.
      SUM = SUM + (RHO*GR)
      Y=Y+DEL
100   CONTINUE

      FPRESS= 0.5*(FPRESS+(DEPTH/1.E4)*SUM/TNM)
      IT=2*IT
      
      ENDIF

50    CONTINUE

c***********************************************************
c
c   Modification as detailed by R. Millard (personal communication)
c    Iterative process to invert depth to pressure, using the
c    calculated pressure (fpress) as a starting point.
c
c************************************************************

      if ( depth .eq. 0. ) return

      p= fpress
      pp= p
      dp= .01

      do 60 m=1,20

       z1= fdepth(p,lat)
       dz= abs(z1 - depth)/depth
       if ( dz .ne. 0 ) p= p + p*dz

       if ( abs( p - pp ) .le. dp ) then
        fpress= p
        return
       endif

       pp= p

60    continue

      FPRESS= p
    
      RETURN
      END