c
      program msmain
c   
c
c  This main program is used to exercise the marginal sea model,
c  mseabc, which computes the inflows and outflows to a marginal sea.  
c  No need to link with anything; needs a data file goczts.dat to
c  define hydrographic conditions in the Gulf of Cadiz.
c
c  Written by Jim Price, WHOI, jprice@whoi.edu, 508-289-2526.
c
c  Version 2.1 in September, 1996, is for the  Mediterranean only. 
c  Last modified in April 1998; cleaned up, width handled slightly
c  differently.  This version 2.3 gives answers very slightly different 
c  from those of Price and Yang (1998).
c
c  The subroutine mseabc computes inflow and outflow boundary conditions for 
c  an ocean GCM as if a marginal sea were present on that boundary.  Air-sea
c  exchange over the marginal sea and the exchange with the ocean
c  are computed by the Bryden and Kinder (1991, DSR) model of maximal exchange
c  by a marginal sea.  Entrainment and descent of the resulting outflow
c  are treated by the Price and Baringer (1994, PinO) model.  
c
c  To run this little gem you have to provide an awful lot of data; 
c  the air-sea fluxes over the marginal sea and the surface area of 
c  the marginal sea, the depth and width of the connecting strait, and
c  the temperature and salinity profile of the ocean just outside the 
c  marginal sea.  The latter would presumably come from the ocean GCM in
c  a case where this subroutine was used in a GCM.  As a means to test the
c  model, a data set,  goczts.dat, can be used to specify 
c  conditions found in the Gulf of Cadiz (Mediterranean outflow).
c
c
c  This common is used to pass data back and forth to the marginal sea
c  subroutine, mseabc.
c
      dimension zocn(100), tocn(100), socn(100)
      dimension tms(4),sms(4),sigims(4),sigpms(4),dms(4),hms(4),
     c wms(4),vms(4),tranms(4), wkm(4), transv(4)
c
      common/marsea/
     & ams,hflux,emp,dsill,width,f,dssb,slope,xssb,
     & zocn,tocn,socn,nocn,
     & tms,sms,sigims,sigpms,dms,hms,wms,vms,tranms,
     & msflag
c
      character*80 aaa
c
c  Variables and their units in this main program are as follows:
c
c  ams is the surface area of the marginal sea [millions of square km],
c  hflux is the heat flux from the sea to the atmosphere [W/m**2],
c  emp is the evaporation minus precipitation rate [m/year],
c  dsill is the depth of the sill [m],
c  width is the width of the strait at the surface [km],
c  lat is the latitude of the shelf-slope break [deg],
c  f is the usual Coriolis parameter [1/sec],
c  dssb is the depth of the shelf-slope break [m],
c  slope is the maximum slope just over the shelf-slope break [non-d],
c  xssb is the distance downstream of the shelf-slope break [km],
c  zocn, tocn, socn are the oceanic z, t, s at the shelf-slope break,
c  nocn is the number of actual data values in the zocn, etc. arrays.
c
c
      new = 0
      write (6,60) 
   60 format (1x,//,1x,' Welcome to the water type factory.',/,
     c'  Version 2.2  (Mediterranean only)',//,
     c 1x,' Are you a new user ? (enter 1 if yes)')
      read (5,*) new
c
c
      if(new.eq.1) then
      write (6,61) 
   61 format (1x,/,1x,
     c' This program exercises the marginal sea BC program. You',/,
     c'  first have to enter a file name of data that define',/,
     c'  the oceanic T-S profiles.  Then you have to specify',/,
     c'  1) the area of the marginal sea, the air-sea fluxes',/, 
     c'  over the marginal sea and its latitude, and,',/,
     c'  2) the geometry of the strait, the depth of the',/,   
     c'  shelf-slope break, and the continental slope.',/,
     c'  Use the mixed units noted below.',/) 
c
      end if
c
c  This main program reads in the oceanic z,t,s profile from an
c  ASCII data file.  If this were a GCM, then this data would come
c  from the ocean just outside the marginal sea strait.
c
      open(unit=1,file='goczts.dat',status='old')
c
c  Note that this file has one line of character data, then free 
c  format triples (z,s,t) with z assumed positive down, i.e., as
c  if depth were a pressure (> 0). Temperature is assumed to be
c  potential temperature.
c
      read (1,33) aaa
   33 format (a80)
      write (6,34) aaa
   34 format (1x,a80)
      do 30 j=1,100
      read (1,*,end=31) zocn(j),tocn(j),socn(j)
c
c  Compute potential and insitu densities.
c
      dpot = sigma(socn(j), tocn(j), 0.)
      tinsit = theta(socn(j),tocn(j),0., zocn(j))
      dinsit = sigma(socn(j), tinsit, zocn(j))
c
      write (6,32) zocn(j),tocn(j),socn(j), dpot, dinsit
      write (2,32) zocn(j),tocn(j),socn(j), dpot, dinsit
   32 format (1x,' oceanic z, t, s, d, di are', 5f8.2)
   30 continue
   31 continue
      nocn =  j - 1
      write (6,35)
   35 format (1x,//)
c
c  Note that some of these are mixed units for now:
c
c  ams is m*m/1.e12,
c  emp is m/year,
c  width is km.
c
  100 continue 
c
c
c_______________________  the Mediterranean outflow  _______________________
c
c  The values set below are for the Mediterranean Sea:
c
      ams = 2.5                ! 10**12 m*m
      hflux = 0.               ! W/m*m, at any rate it is pretty small.
      emp = 0.7                ! m/year, uncertain to about 0.1 
      rlat = 36.
c
      dsill = 280.             ! the shallowest sill
      width = 20.              ! the width of the strait at the surface
      dssb = 400.              ! dept of the shelf-slope break
      slope = 12.e-3           ! bottom slope just over the ssb
      xssb = 20.e3             ! distance from the strait to the ssb 
c 
c
      if(new.eq.1) then
      write (6,62) 
   62 format (1x, 
     c ' ams is the area of the marginal sea [M km**2], ',/,  
     c '  hflux is the annual average heat flux [W/m**2], ',/,  
     c '  emp is the evaporation minus precipitation [m/yr], ',/,  
     c '  lat is the latitude [deg]. ') 
      end if
c   
      write (6,10) ams, hflux, emp, rlat
   10 format (1x,/,1x,' ams, hflux, emp, rlat are now', 4f9.2)
      write (6,19) 
   19 format (1x,' enter new values, or default with ,,,,,')   
      read (5,*) ams, hflux, emp, rlat
c
      write (2,10) ams, hflux, emp, rlat
c 
      if(new.eq.1) then
      write (6,63) 
   63 format (1x,/,1x,
     c' dsill is the depth of the sill [m],',/,
     c'  width is the width of the strait [km],',/, 
     c'  dssb is the depth of the shelf-slope break [m], ',/,   
     c'  slope is the slope of the continental slope [non-d].') 
      end if
c
      write (6,11) dsill, width, dssb, slope
   11 format (1x,/,1x,' dsill, width, dssb, slope are now', 
     c  3f9.2,f9.5)
      write (6,19)
      read (5,*) dsill, width, dssb, slope
c
  211 format (1x,/,1x,' dsill, width, dssb, slope are now', 
     c  3f9.2,f9.5,/)
      write (2,211) dsill, width, dssb, slope
c 
c  Convert to MKS units where required.
c
      ams = ams*1.e12
      emp = emp/3.154e7
      width = width*1.e3
c
      f = 2.*7.292e-5*sin(rlat*3.1415/180.)
      if(dssb.lt.dsill) dssb = dsill
c
c  Everything has now been specified; call the marginal sea routine. 
c
      call mseabc
c
c  The output comes back in these arrays:
c
c  tms is the temperature [C],
c  sms is the salinity [non-d],
c  sigms is the insitu density [kg/m*3],
c  hms is the thickness [m],
c  dms is the depth (assumes that z is positive down) [m],
c  wms is the width [m]
c  vms is the current, (negative means flowing out of the ocean) [m/sec],
c  tranms is the transport, sign as above [m^3/sec]
c
c  Array elements are as follows:
c
c  tms(1) is the temperature of the oceanic water that inflows 
c    to the marginal sea.
c  tms(2) is the temperature of the water that outflows from
c    the marginal sea into the ocean. This is the initial condition
c    for the outflow model. 
c  tms(3) is the temperature of the oceanic water that is entrained 
c    by the descending outflow. This is the oceanic water at the depth
c    of the shelf-slope break. 
c  tms(4) is the final temperature of the outflow, i.e., the temperature
c    of the product water that will be on the bottom or neutrally buoyant
c    in the open ocean.
c
c  If a variable doesn't have a defined value (i.e., the speed and the
c  thickness of the final outflow depends upon the bottom slope somewhere 
c  beyond the shelf-slope break) then it will be set to a default, 999.
c
c  Store the output data to a file called mseabc.dat (unit 2).
c
      open(unit=2,file='mseabc.dat',form='formatted',status='unknown')
c
      do 40 j=1,4
      write (2,41) dms(j),tms(j),sms(j),sigpms(j),sigims(j),
     c hms(j),vms(j),wms(j),tranms(j)
   41 format (1x, 7f9.2, 2f12.0)
   40 continue
      write (2,47)
   47 format (1x,//) 
c
c  Convert width and transports to km and Sv for screen output.
c
      do 723 k = 1,4
      wkm(k) = wms(k)/1000.
      transv(k) = tranms(k)/1.e6
  723 continue
c
c  Write out the inflow and marginal sea outflow data to the terminal.
c
      write (6,70)
   70 format (1x,/,'  Output from the water type factory is:',/,)
c
      write (6,935) 
  935 format (1x,10x,
     c  '      z     U       T       S    sigpot  sigins',
     c '     h     width   trans')
c
      write (6,21) dms(1), vms(1), tms(1), sms(1), sigpms(1),
     c   sigims(1), hms(1), wkm(1), transv(1)
   21 format (1x,/, '  inflow   ', f7.0,9f8.2)
c
      write (6,22) dms(2), vms(2), tms(2), sms(2), sigpms(2),
     c   sigims(2), hms(2), wkm(2), transv(2)
   22 format (1x,' source   ',f7.0,9f8.2)
c 
      wkm = 999.
      write (6,23) dms(3), vms(3), tms(3), sms(3), sigpms(3),
     c   sigims(3), hms(3), wkm(3), transv(3)
   23 format (1x,' entrained',f7.0,9f8.2)
c
      write (6,24) dms(4), vms(4), tms(4), sms(4), sigpms(4),
     c   sigims(4), hms(4), wkm(4), transv(4)
   24 format (1x,' product  ',f7.0,9f8.2)
c
      write (6,931) msflag
  931 format (1x,/,1x,' The msflag for this case was',i3,', where',/,
     c  '   = 0 means that all is normal,',/,
     c  '   = 1 there was no entrainment (which is OK)',/,
     c  '   = 2 the overflow went to the bottom (also OK)',/,
     c  '   < 0 there is big trouble and you may want to stop.') 

      iagain = 0
      write (6,90)
   90 format (1x,//,'  Do you want to rerun this case? (1 = yes)')
      read (5,*) iagain
      if(iagain.eq.1) go to 100
   99 continue
c         
      stop
      end
c 
c _____________________ end of the main program _________________________
c 
c
      subroutine mseabc
c
c
c  Written by Jim Price, WHOI, started in November 1994. 
c
c  This version 2.1 in September 1996. 
c
c  This subroutine computes inflow and outflow boundary conditions
c  for an OGCM as if a marginal sea were on that boundary.  Air-sea
c  interaction over the marginal sea and the exchange with the ocean
c  are computed by the Bryden and Kinder 1991 (Deep-Sea Research, 38,
c  S1, 445-463) model of exchange by a marginal sea.  Entrainment 
c  and descent of the resulting outflow are treated by the Price
c  and Baringer 1994 model (Prog. in Occeanography, 33, 161-200).
c
c
c  To run this little model you have to specify a lot of stuff about
c  the marginal sea and the continental shelf.   
c
c  Air-sea fluxes and surface area of the marginal sea:
c  1) ams is the surface area of the marginal sea, 
c  2) hflux is the heat flux from the sea to the atmosphere (a
c       negative value indicates heat loss from the sea),
c  3) emp is the evaporation minus precipitation rate  (a positive
c       value indicates that there is more evaporation),
c 
c  Geometry of the marginal sea strait and adjacent continental
c       slope:
c  4) width is the width of the connecting strait,
c  5) dsill is the depth of the sill (all depths are positive),
c  6) f is the Coriolis parameter of the strait,
c  7) dssb is the depth of the shelf-slope break,
c  8) slope is the maximum slope just over the shelf-slope break,
c  9) xssb is the distance downstream of the shelf-slope break.
c 
c  Temperature and salinity profile of the ocean just outside the 
c       marginal sea:
c  10 - 12) zocn, tocn, and socn are the depth, temperature and salinity of 
c       oceanic water just over the shelf-slope break (z is positive
c       downwards), and finally,
c  13) nocn is the number of elements in zocn, etc.
c
c
c  The output comes back in the following arrays:
c
c  tms(1) is the temperature of the oceanic water that inflows 
c    to the marginal sea. All temperatures are potential temps.
c  tms(2) is the temperature of the water that outflows from
c    the marginal sea into the ocean. This is the initial condition
c    for the outflow model. 
c  tms(3) is the temperature of the oceanic water that is entrained 
c    by the descending outflow.
c  tms(4) is the final temperature of the equilibrated outflow.
c
c  sms is the salinity
c  hms is the thickness
c  dms is the mid-depth of a layer (z is positive down as if pressure)
c  wms is the width
c  vms is the current, negative means flowing out of the ocean
c  tranms is the transport, signs as above. 
c
c  If some field does not make sense, then a default 999. 
c    is inserted into the array.
c
c  Other variables used internally in the subroutine are:
c
c  t1 and s1 are the inflow temperature and salinity (same as tms(1)). 
c  t2 and s2 are the outlfow temperature and salinity.  
c  tssb and sssb are the oceanic t,s values at the shelf slope break
c  tf and sf are the final temperature, salinity etc. of the equilibrated
c    outflow.
c  h1 and h2 are the thicknesses of the inflow and outflow layers (for 
c   now assumed equal).  Same for u1 and u2.
c
c  So far as the marginal sea goes, the only fluxes of importance are
c  the 'inflow' and the 'source' (which flows out of the marginal sea). So
c  far as the OGCM goes, the only fluxes of importance are the inflow
c  and the 'entrained' water (which are both flowing out of the ocean), and 
c  'product' water which is the equilibrated outflow and is a source to 
c  the ocean.  Note that the 'source' water is not relevent to the GCM. 
c
c  The status of the calculation is given by the flag, msflag, as follows:
c   msflag = 0 means that the subroutine seemed to function OK.
c   msflag = 1 indicates that there was no entrainment.
c   msflag = 2 means that the product water went to the bottom of the ocean.
c   msflag < 0 means something was badly wrong and you should stop the model. 
c
c  All of the data in this subroutine is in MKS units. 
c
c
c __________________ Data Storarge and Common Allocation _________________
c
c  This common is used to pass data back and forth to the marginal sea
c  subroutine, mseabc.  This should appear also in the main program. 
c
      dimension zocn(100), tocn(100), socn(100)
      dimension tms(4),sms(4),sigims(4),sigpms(4),dms(4),hms(4),
     c wms(4),vms(4),tranms(4)
c
      common/marsea/
     & ams,hflux,emp,dsill,width,f,dssb,slope,xssb,
     & zocn,tocn,socn,nocn,
     & tms,sms,sigims,sigpms,dms,hms,wms,vms,tranms,
     & msflag
c    
c
      msflag = 0
c
c _______________ Exchange Between the Marginal Sea and the Ocean ____________
c
c  The first task is to evaluate inflow and outflow from the marginal sea
c  using the Bryden and Kinder (1991, DSR) results for maximal two-layer 
c  exchange between a marginal sea and the open ocean.  This version
c  assumes that the exchange is driven entirely by the salinity difference,
c  and so is suitable for the Mediterranean only.
c 
c  g is the acceleration of gravity
c  rho is the nominal density of sea water
c  cp is the heat capacity of sea water
c
      g = 9.8
      rho = 1026.
c
c  Find the oceanic t,s at the surface and at the outflow interface. 
c
      h1 = dsill/2.
c
      call intx(zocn,tocn,nocn,0.,tsurf)
      call intx(zocn,tocn,nocn,h1,tmid)
      call intx(zocn,socn,nocn,0.,ssurf)
      call intx(zocn,socn,nocn,h1,smid)
c
c  Average these to estimate the inflow t,s.
c
      s1 = (ssurf + smid)/2.
      t1 = (tsurf + tmid)/2. 
c      
      sig01 = sigma(s1,t1,0.)
      cp = cpsw(s1,t1,0.)
c
c  Evaluate alpha and beta (the expansion coefficients) at the sea surface. 
c
      alpha = (1./rho)*(sigma(s1,t1+0.5,0.) - sigma(s1,t1-0.5,0.))
      beta = (1./rho)*(sigma(s1+0.5,t1,0.) - sigma(s1-0.5,t1,0.))
c 
c  Compute the net evaporation.
c
      enet = ams*emp
c
c  This version assumes maximal exchange, and thus that h1 = h2 = dsill/2.  
c
      h2 = h1
c
c  Find the salinity change of the outflowing water, deltas, and
c  from that compute the speed and transport of the outflow.  This 
c  requires finding the root of a polynomial.
c
      deltas = froot(beta,h2,width,s1,enet,irflag)
      if(irflag.eq.1.or.deltas.lt.0.) then
      msflag = -1
      return                		!  just quit if deltas < 0
      end if
c
      s2 = s1 + deltas
c
c  Compute the overflow transport.
c
      tran2 = enet*s2/deltas
      write (6,811) tran2
  811 format (1x, /////, 1x, '  tran2 is ', e12.3)
c
c  The width used in the exchange theory above is the width of the strait
c  at the surface.  Assuming a triangular cross-section, then the width of 
c  the lower layer is about 1/root2 as much, as is layer thickness.
c
      width2 = width/sqrt(2.)
      h2 = h2/sqrt(2.)
      u2 = tran2/(width2*h2)
c
c  Given the known exchange, compute the temperature change (here
c  temperature is assumed to be a passive tracer - OK for the Med only).
c
      deltat = ams*(hflux/(rho*cp))/(u2*width*h2)
      t2 = t1 + deltat 
      sig02 = sigma(s2,t2,0.)
c
c  Subject to the h1 = h2 assumption, the upper layer current and
c  transport are the same as the lower layer's  (the sign flip is to show
c  water leaving the ocean).  This ignores the small barotropic flow into
c  the marginal sea required to balance the volume loss due to evaporation.
c
      u1 = -u2
      tran1 = -tran2
c
c  Similarly, the layer mid-depths are:
c
      d1 = h1/2.
      d2 = h1 + h2/2.
c
c  Compute the insitu densities at these levels.
c
      ti1 = theta(s1, t1, 0., d1)
      sigi1 = sigma(s1, ti1, d1)
      ti2 = theta(s2, t2, 0., d2)
      sigi2 = sigma(s2, ti2, d2)
c
c  Load the data into the output arrays.
c
      tms(1) = t1
      tms(2) = t2
      sms(1) = s1
      sms(2) = s2
      sigpms(1) = sig01
      sigpms(2) = sig02
      hms(1) = h1
      hms(2) = h2
      dms(1) = d1
      dms(2) = d2
      sigims(1) = sigi1
      sigims(2) = sigi2
      vms(1) = u1
      vms(2) = u2
      wms(1) = width
      wms(2) = width2
      tranms(1) = tran1
      tranms(2) = tran2
c
c ________________________  Descent and Entrainment  _________________
c   
c
c  Given the outflow from the marginal sea, we now compute its
c  evolution as it descends the continental slope.  This is done with the
c  Price and Baringer (1994) overflow model. 
c
c  Fine the oceanic t,s and density at the depth of the shelf-slope break,
c  which is at the depth z = dssb.
c
      call intx(zocn,tocn,nocn,dssb,tssb)
      call intx(zocn,socn,nocn,dssb,sssb)
c
c  The notation sig0x indicates the potential density of layer x, while
c  sigx is the insitu density.  Insitu density should be used in density
c  comparisons that need to account for pressure dependent expansion
c  coefficients.  However, some ocean GCMs use potential density, and
c  to be consistent with a GCM it may be desirable to use potential
c  density throughout. This can be done by setting the flag ipot = 1  
c
      ipot = 0
c      if(ipot.eq.1) then
c      write (6,2245) 
c 2245 format (1x,//,1x,
c     c ' Watch out big guy, you are using potential density!',//)
c      end if
c
      sig0ssb = sigma(sssb,tssb,0.)
c
      zssb = dssb
      if(ipot.eq.1) zssb = 0.
c
      t2i = theta(s2,t2,0.,zssb)
      sig2 = sigma(s2,t2i,zssb)
c
      tssbi = theta(sssb,tssb,0.,zssb)
      sigssb = sigma(sssb,tssbi,zssb)
c
      deltad = sig2 - sigssb
c  
c  deltad is assumed positive for a marginal sea producing dense source water. 
c
c  Check to make sure that this marginal sea will indeed produce a
c  deep outflow. If not, jump down a bit.
c
      if(deltad.lt.0.) go to 8811
c 
c  Given the density difference, compute the gesotrophic speed.
c     
      ugeo = g*deltad*slope/(f*rho)     
      bflux = g*h2*(deltad/rho)*u2
c 
c  Apply some broadening to the outflow. First, compute the Ekman number, 
c  based upon the average speed up to the shelf-slope break.
c
      cd = 5.e-3    !  a biggish value, Baringer and Price, 1997 JPO
c
      uavg = (ugeo + u2)/2.
      h2geo = h2*u2/uavg
      E = cd*uavg/(h2geo*f)
      width0 = width
c
c  Compute how much the buoyancy flux per unit width would decrease due
c  to broadening of the outflow as it moves to the shelf-slope break. This
c  is the same as the width increase.
c
      width = width0 + 2.*E*xssb
      broad = width0/width
      bflux = bflux*broad
      wkm  = width/1000.
c
c  Compute phi and check to see if it is positive.
c
      phi = 1. - (bflux**0.333)/ugeo
c
c      write (6,3333) E, deltad, ugeo, bflux, broad, phi
c 3333 format (1x,/,1x,
c     1  'E, deltad, ugeo, bflux, broad, phi', 6e10.3,/)
c 
 8811 continue           
      if(phi.lt.0..or.deltad.lt.0.) then
c      write (6,110) deltad, bflux, ugeo, phi
c  110 format 
c     c (1x,/,1x,' deltad or phi less than 0.; bf, ug, phi',4e11.3,/)
      phi = 0.
      msflag = 1
      end if      
c
c  tf, sf, and sigf refer to the product water t,s etc. trane is the 
c  volume transport of the water that is entrained.  The assumption is
c  that this entrained water comes from the oceanic water column at
c  the depth of the shelf-slope break.
c
      tf = t2 - phi*(t2 - tssb)
      sf = s2 - phi*(s2 - sssb)
      sigf = sig2 - phi*deltad 
      sig0f = sigma(sf,tf,0.)
c
      tranf = tran2/(1. - phi)
      trane = -(tranf - tran2)
c 
c      write (6,2256) phi, tran2, trane, tranf
c 2256 format (1x, '  phi, 3tran', 3e12.3)
c
c  Now figure out where this final outflow water would become neutrally
c  buoyant in the oceanic water column.  Use insitu density for this unless
c  ipot = 1, in which case use potential density.
c
      zbot = zocn(nocn)
      do 20 j=1,400
      ztest = float(j)*25.
      if(ztest.gt.zbot) then
      if(msflag.ge.0) msflag = 2
      go to 21
      end if
c
      zdens = ztest
      if(ipot.eq.1) zdens = 0.
c
      tfi = theta(sf,tf,0.,zdens)
      sigfi = sigma(sf,tfi,zdens)
      call intx(zocn,tocn,nocn,zdens,toz)
      call intx(zocn,socn,nocn,zdens,soz)
      tocni = theta(soz, toz, 0., zdens)
      sigocni = sigma(soz,tocni,zdens)
      if(sigfi.lt.sigocni) go to 21
   20 continue
c
c  Come here to break out of the ztest.
c
   21 continue
      zneut = ztest - 12. 
c
      ti3 = theta(sssb, tssb, 0., dssb)
      sigissb = sigma(sssb, tissb, dssb)
      ti4 = theta(sf, tf, 0., zneut)
      sigif = sigma(sf, ti4, zneut) 
c
c  Write out data on the entrained water and final (or product) water.
c
      vms(3) = 999.
      vms(4) = 999.
      h3 = 999.
      h4 = 999.
c
c  Load these data into the output arrays.  Insert a default 
c  999. if an answer would be meaningless.  Specifically, the 
c  height and width of the entrained water and the final outflow 
c  depend entirely upon the bottom slope, and so can not and need not be
c  specified here.  The transports must all be specified, however.
c
      tms(3) = tssb
      tms(4) = tf
      sms(3) = sssb
      sms(4) = sf
      sigpms(3) = sig0ssb
      sigpms(4) = sig0f
      hms(3) = 999.
      hms(4) = 999.
      dms(3) = dssb
      dms(4) = zneut
      sigims(3) = sigissb 
      sigims(4) = sigif 
      vms(3) = 999.
      vms(4) = 999.
      wms(3) = width
      wms(4) = width
      tranms(3) = trane
      tranms(4) = tranf
c
      return
      end
c __________________ End of the Marginal Sea Model _______________
c
c
      subroutine intx(r,s,n,ri,si)
c
c  This subroutine interploates the array s to find the value 
c  at the point ri, where r is the independent variable.
c  The interpolation is linear, and assumes that r increases
c  montonically, and that ri is within the range of r. If
c  ri is not within the range of r, si is set equal to the
c  appropriate endpoint of s, and a warning is sent out.  

c
      dimension r(100), s(100)
c 
c  Check for out of range ri.
c
      if(ri.lt.r(1)) then
      si = s(1)
c      write (6,20) ri, r(1)
c  20  format (1x,/,1x,'warning, ri was below r(1) in call to intx',/,
c     c 'ri and r(1) were', 2e12.4)
      end if
c
c
      if(ri.gt.r(n)) then
      si = s(n)
c      write (6,21) ri, r(n)
c  21  format (1x,/,1x,'warning, ri was above r(n) in call to intx',/,
c     c 'ri and r(n) were', 2e12.4)
      end if    
c 
c  The data are within bounds for an interpolation.
c
      nm = n - 1
      do 2 k=1,nm
      if(ri.ge.r(k).and.ri.le.r(k+1)) go to 1
    2 continue
c
      go to 9
c
    1 continue
      si = (s(k)*(r(k+1) - ri) + s(k+1)*(ri - r(k)))/
     x (r(k+1) - r(k))
c
    9 continue
      return
      end
c
c_____________________________________________________________________
c
      function sigma(s,t,p)
c
c  Computes sigma from the svan function of seaprop.
c
      qq = svan(s,t,p,sig)
      sigma = sig
      return
      end
c
c  This subroutine finds a root (i.e., it is a root finder,
c  albeit a very dumb one that makes very special assumptions
c  about the polynomial). 
c 
      function froot(beta,h2,width,s1,enet,irflag)
c
      irflag = 0
      g = 9.8
c
      do 1 j=1,300
      ds = float(j)*0.02
      term1 = 0.283*sqrt(g*beta*ds*2.*h2)*width*h2
      term2 = enet*(2.*s1 + ds)/ds
      sumold = sum
      sum = term1 - term2
c      write(6,2) j, ds, term1, term2, sum
c    2 format (1x,'   in froot', i4, 4e12.4)
      if(sum.ge.0.) go to 3
    1 continue
      write (6,4) 
    4 format (1x,'   Root not found by froot',//)
      froot = 0.
      irflag = 1
      return
c
c  At the zero crossing, come here.
c
    3 continue
      deriv = (sum - sumold)/0.02
      dspart = -sumold/deriv
      ds = ds + dspart - 0.02 
      froot = ds
c
c  Make sure the ds (the root) is positive.
c
      if(ds.lt.0.) then
      ds = 0.
      irflag = 1
      end if
c
      return
      end
c 
c_______________________________ Seaprop ____________________________
c
c  FORTRAN below is from SEAPROP by Bob Millard.  It is included here 
c  to make this program entirely self-contained.
c
      REAL FUNCTION SVAN(S,T,P0,SIGMA)
C  MODIFIED RCM
C ******************************************************
C SPECIFIC VOLUME ANOMALY (STERIC ANOMALY) BASED ON 1980 EQUATION
C OF STATE FOR SEAWATER AND 1978 PRACTICAL SALINITY SCALE.
C REFERENCES
C MILLERO, ET AL (1980) DEEP-SEA RES.,27A,255-264
C MILLERO AND POISSON 1981,DEEP-SEA RES.,28A PP 625-629.
C BOTH ABOVE REFERENCES ARE ALSO FOUND IN UNESCO REPORT 38 (1981)
C UNITS:
C       PRESSURE        P0       DECIBARS
C       TEMPERATURE     T        DEG CELSIUS (IPTS-68)
C       SALINITY        S        (IPSS-78)
C       SPEC. VOL. ANA. SVAN     M**3/KG *1.0E-8
C       DENSITY ANA.    SIGMA    KG/M**3
C ******************************************************************
C CHECK VALUE: SVAN=981.3021 E-8 M**3/KG.  FOR S = 40 (IPSS-78) ,
C T = 40 DEG C, P0= 10000 DECIBARS.
C CHECK VALUE: SIGMA = 59.82037  KG/M**3 FOR S = 40 (IPSS-78) ,
C T = 40 DEG C, P0= 10000 DECIBARS.
C *******************************************************
      REAL P,T,S,SIG,SR,R1,R2,R3,R4
      REAL A,B,C,D,E,A1,B1,AW,BW,K,K0,KW,K35
C EQUIV
      EQUIVALENCE (E,D,B1),(BW,B,R3),(C,A1,R2)
      EQUIVALENCE (AW,A,R1),(KW,K0,K)
C ********************
C DATA
      DATA R3500,R4/1028.1063,4.8314E-4/
      DATA DR350/28.106331/
C   R4 IS REFERED TO AS  C  IN MILLERO AND POISSON 1981
C CONVERT PRESSURE TO BARS AND TAKE SQUARE ROOT SALINITY.
      P=P0/10.
      SR = SQRT(ABS(S))
C *********************************************************
C PURE WATER DENSITY AT ATMOSPHERIC PRESSURE
C   BIGG P.H.,(1967) BR. J. APPLIED PHYSICS 8 PP 521-537.
C
      R1 = ((((6.536332E-9*T-1.120083E-6)*T+1.001685E-4)*T
     X-9.095290E-3)*T+6.793952E-2)*T-28.263737
C SEAWATER DENSITY ATM PRESS.
C  COEFFICIENTS INVOLVING SALINITY
C  R2 = A   IN NOTATION OF MILLERO AND POISSON 1981
      R2 = (((5.3875E-9*T-8.2467E-7)*T+7.6438E-5)*T-4.0899E-3)*T
     X+8.24493E-1
C  R3 = B  IN NOTATION OF MILLERO AND POISSON 1981
      R3 = (-1.6546E-6*T+1.0227E-4)*T-5.72466E-3
C  INTERNATIONAL ONE-ATMOSPHERE EQUATION OF STATE OF SEAWATER
      SIG = (R4*S + R3*SR + R2)*S + R1
C SPECIFIC VOLUME AT ATMOSPHERIC PRESSURE
      V350P = 1.0/R3500
      SVA = -SIG*V350P/(R3500+SIG)
      SIGMA=SIG+DR350
C  SCALE SPECIFIC VOL. ANAMOLY TO NORMALLY REPORTED UNITS
      SVAN=SVA*1.0E+8
      IF(P.EQ.0.0) RETURN
C ******************************************************************
C ******  NEW HIGH PRESSURE EQUATION OF STATE FOR SEAWATER ********
C ******************************************************************
C        MILLERO, ET AL , 1980 DSR 27A, PP 255-264
C               CONSTANT NOTATION FOLLOWS ARTICLE
C********************************************************
C COMPUTE COMPRESSION TERMS
      E = (9.1697E-10*T+2.0816E-8)*T-9.9348E-7
      BW = (5.2787E-8*T-6.12293E-6)*T+3.47718E-5
      B = BW + E*S
C
      D = 1.91075E-4
      C = (-1.6078E-6*T-1.0981E-5)*T+2.2838E-3
      AW = ((-5.77905E-7*T+1.16092E-4)*T+1.43713E-3)*T
     X-0.1194975
      A = (D*SR + C)*S + AW
C
      B1 = (-5.3009E-4*T+1.6483E-2)*T+7.944E-2
      A1 = ((-6.1670E-5*T+1.09987E-2)*T-0.603459)*T+54.6746
      KW = (((-5.155288E-5*T+1.360477E-2)*T-2.327105)*T
     X+148.4206)*T-1930.06
      K0 = (B1*SR + A1)*S + KW
C EVALUATE PRESSURE POLYNOMIAL
C ***********************************************
C   K EQUALS THE SECANT BULK MODULUS OF SEAWATER
C   DK=K(S,T,P)-K(35,0,P)
C  K35=K(35,0,P)
C ***********************************************
      DK = (B*P + A)*P + K0
      K35  = (5.03217E-5*P+3.359406)*P+21582.27
      GAM=P/K35
      PK = 1.0 - GAM
      SVA = SVA*PK + (V350P+SVA)*P*DK/(K35*(K35+DK))
C  SCALE SPECIFIC VOL. ANAMOLY TO NORMALLY REPORTED UNITS
      SVAN=SVA*1.0E+8
      V350P = V350P*PK
C  ****************************************************
C COMPUTE DENSITY ANAMOLY WITH RESPECT TO 1000.0 KG/M**3
C  1) DR350: DENSITY ANAMOLY AT 35 (IPSS-78), 0 DEG. C AND 0 DECIBARS
C  2) DR35P: DENSITY ANAMOLY 35 (IPSS-78), 0 DEG. C ,  PRES. VARIATION
C  3) DVAN : DENSITY ANAMOLY VARIATIONS INVOLVING SPECFIC VOL. ANAMOLY
C ********************************************************************
C CHECK VALUE: SIGMA = 59.82037  KG/M**3 FOR S = 40 (IPSS-78),
C T = 40 DEG C, P0= 10000 DECIBARS.
C *******************************************************
      DR35P=GAM/V350P
      DVAN=SVA/(V350P*(V350P+SVA))
      SIGMA=DR350+DR35P-DVAN
      RETURN
      END

C
C      
      REAL FUNCTION TF(S,P)
C
C

C   FUNCTION TO COMPUTE THE FREEZING POINT OF SEAWATER
C
C   REFERENCE: UNESCO TECH. PAPERS IN THE MARINE SCIENCE NO. 28. 1978
C   EIGHTH REPORT JPOTS
C   ANNEX 6 FREEZING POINT OF SEAWATER F.J. MILLERO PP.29-35.
C
C  UNITS:
C         PRESSURE      P          DECIBARS
C         SALINITY      S          PSS-78
C         TEMPERATURE   TF         DEGREES CELSIUS
C         FREEZING PT.
C************************************************************
C  CHECKVALUE: TF= -2.588567 DEG. C FOR S=40.0, P=500. DECIBARS
      TF=(-.0575+1.710523E-3*SQRT(ABS(S))-2.154996E-4*S)*S-7.53E-4*P
      RETURN
      END
      REAL FUNCTION CPSW(S,T,P0)
C ****************************
C UNITS:
C       PRESSURE        P0       DECIBARS
C       TEMPERATURE     T        DEG CELSIUS (IPTS-68)
C       SALINITY        S        (IPSS-78)
C       SPECIFIC HEAT   CPSW     J/(KG DEG C)
C********************************************************
C REF: MILLERO ET AL,1973,JGR,78,4499-4507
C       MILLERO ET AL, UNESCO REPORT NO. 38 1981 PP. 99-188.
C PRESSURE VARIATION FROM LEAST SQUARES POLYNOMIAL
C DEVELOPED BY FOFONOFF 1980.
C CHECK VALUE: CPSW = 3849.500 J/(KG DEG. C) FOR S = 40 (IPSS-78),
C T = 40 DEG C, P0= 10000 DECIBARS
C   SCALE PRESSURE TO BARS
      P=P0/10.
C**************************
C SQRT SALINITY FOR FRACTIONAL TERMS
      SR = SQRT(ABS(S))
C SPECIFIC HEAT CP0 FOR P=0 (MILLERO ET AL ,UNESCO 1981)
      A = (-1.38385E-3*T+0.1072763)*T-7.643575
      B = (5.148E-5*T-4.07718E-3)*T+0.1770383
      C = (((2.093236E-5*T-2.654387E-3)*T+0.1412855)*T
     X    -3.720283)*T+4217.4
      CP0 = (B*SR + A)*S + C
C CP1 PRESSURE AND TEMPERATURE TERMS FOR S = 0
      A = (((1.7168E-8*T+2.0357E-6)*T-3.13885E-4)*T+1.45747E-2)*T
     X   -0.49592
      B = (((2.2956E-11*T-4.0027E-9)*T+2.87533E-7)*T-1.08645E-5)*T
     X   +2.4931E-4
      C = ((6.136E-13*T-6.5637E-11)*T+2.6380E-9)*T-5.422E-8
      CP1 = ((C*P+B)*P+A)*P
C CP2 PRESSURE AND TEMPERATURE TERMS FOR S > 0
      A = (((-2.9179E-10*T+2.5941E-8)*T+9.802E-7)*T-1.28315E-4)*T
     X   +4.9247E-3
      B = (3.122E-8*T-1.517E-6)*T-1.2331E-4
      A = (A+B*SR)*S
      B = ((1.8448E-11*T-2.3905E-9)*T+1.17054E-7)*T-2.9558E-6
      B = (B+9.971E-8*SR)*S
      C = (3.513E-13*T-1.7682E-11)*T+5.540E-10
      C = (C-1.4300E-12*T*SR)*S
      CP2 = ((C*P+B)*P+A)*P
C SPECIFIC HEAT RETURN
      CPSW = CP0 + CP1 + CP2
      RETURN
      END
C
C
C      
      REAL FUNCTION ATG(S,T,P)
C ****************************
C ADIABATIC TEMPERATURE GRADIENT DEG C PER DECIBAR
C REF: BRYDEN,H.,1973,DEEP-SEA RES.,20,401-408
C UNITS:
C       PRESSURE        P        DECIBARS
C       TEMPERATURE     T        DEG CELSIUS (IPTS-68)
C       SALINITY        S        (IPSS-78)
C       ADIABATIC       ATG      DEG. C/DECIBAR
C CHECKVALUE: ATG=3.255976E-4 C/DBAR FOR S=40 (IPSS-78),
C T=40 DEG C,P0=10000 DECIBARS
      DS = S - 35.0
      ATG = (((-2.1687E-16*T+1.8676E-14)*T-4.6206E-13)*P
     X+((2.7759E-12*T-1.1351E-10)*DS+((-5.4481E-14*T
     X+8.733E-12)*T-6.7795E-10)*T+1.8741E-8))*P
     X+(-4.2393E-8*T+1.8932E-6)*DS
     X+((6.6228E-10*T-6.836E-8)*T+8.5258E-6)*T+3.5803E-5
      RETURN
      END
C
C
C      
      REAL FUNCTION THETA(S,T0,P0,PR)
C ***********************************
C TO COMPUTE LOCAL POTENTIAL TEMPERATURE AT PR
C USING BRYDEN 1973 POLYNOMIAL FOR ADIABATIC LAPSE RATE
C AND RUNGE-KUTTA 4-TH ORDER INTEGRATION ALGORITHM.
C REF: BRYDEN,H.,1973,DEEP-SEA RES.,20,401-408
C FOFONOFF,N.,1977,DEEP-SEA RES.,24,489-491
C UNITS:
C       PRESSURE        P0       DECIBARS
C       TEMPERATURE     T0       DEG CELSIUS (IPTS-68)
C       SALINITY        S        (IPSS-78)
C       REFERENCE PRS   PR       DECIBARS
C       POTENTIAL TMP.  THETA    DEG CELSIUS
C CHECKVALUE: THETA= 36.89073 C,S=40 (IPSS-78),T0=40 DEG C,
C P0=10000 DECIBARS,PR=0 DECIBARS
C
C      SET-UP INTERMEDIATE TEMPERATURE AND PRESSURE VARIABLES
      P=P0
      T=T0
C**************
      H = PR - P
      XK = H*ATG(S,T,P)
      T = T + 0.5*XK
      Q = XK
      P = P + 0.5*H
      XK = H*ATG(S,T,P)
      T = T + 0.29289322*(XK-Q)
      Q = 0.58578644*XK + 0.121320344*Q
      XK = H*ATG(S,T,P)
      T = T + 1.707106781*(XK-Q)
      Q = 3.414213562*XK - 4.121320344*Q
      P = P + 0.5*H
      XK = H*ATG(S,T,P)
      THETA = T + (XK-2.0*Q)/6.0
      RETURN
      END
C
