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