c program norsea c c This main program is used to exercise the marginal sea model, c mseabc, which computes the inflows and outflows to a marginal c sea. No links required to compile; needs a data file irmzts.dat c to define hydrographic conditions just south of Denmark Strait. c c Written by Jim Price, WHOI, jprice@whoi.edu, 508-289-2526. c c Version 1.0 was complete in November 1994. c Version 2.1 in September, 1996. c c This version 2.3 of March 1998 is for Nordic Seas only. 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 Whitehead et al (1974, GAFD) 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, irmzts.dat, can be used to specify c conditions found in the Irminger Sea (Denmark St. 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,iwhich 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 [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.3 (Nordic Seas 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='irmzts.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 write (6,111) 111 format (1x,/,1x, c ' Set default to DSt or Combined Nordic Sea (1 or 2)') read (5,*) mordst c c 100 continue c 120 continue c c_________________ a Norwegian-Greenland Sea-like outflow __________ c c The values set below look like the Norwegian-Greenland Sea and c give inflow/outflow transports and T/S that are very roughly the c Denmark St. and Faroe Bank Channel outflows added together. c c The area of the marginal sea should be the ice-free area, presumably. c ams = 1000.e3*2700.e3/1.e12 c hflux = -40. c emp = -0.1 hflux = -30. ! this is about half, since that goes thru DSt. emp = -0.08 rlat = 65. width = 40.e3/1.e3 c dssb = 1000. slope = 2.8e-2 xssb = 100.e3 c dsill = 600. hupper = 200. c if(mordst.eq.2) then c c to make a combined Nordic Sea (DSt plus FBC) case you can use c the following (but are you still using the irmzts.dat ocean profile?) c (this will not be the same as the combined Nordic Sea case in c Price and Yang (1998), which was coupled to MOM). c dsill = 700. ! sill depth more like Faroe Bank Channel width = 50. hflux = -60. emp = -0.15 c end if c 129 continue 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 [km] c vms is the current, (negative means flowing out of the ocean) [m/sec], c tranms is the transport, sign as above [Sv] 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),wms(j),vms(j),tranms(j) 41 format (1x, 5f9.2,2f9.0,2f9.2) 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 overflow went to a mid-depth, (which is OK)',/, c ' = 1 there was no entrainment (which is OK)',/, c ' = 2 the overflow went to the bottom (which is OK)',/, c ' < 0 there is big trouble and you may want to stop.') c 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, in November 1994. c c This version 2.3 in April 1998. 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 Whitehead et al 1974 (Geophys Fluid Dynamics, c (101-125) 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,iwhich 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. As of V 2.0, this model 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 = (sigma(s1,3.0,0.) - sigma(s1,2.0,0.)) beta = (sigma(s1+0.5,t1,0.) - sigma(s1-0.5,t1,0.)) c c set alpha to be a cold water value c c c write (6,700) alpha, beta c 700 format (1x,////,1x,'alpha, beta', 2e12.3) c c Use the Whitehead et al 1974 version of geostrophic exchange c hupper = 200. ! this is the thickness of the surface layer S0 = 36. c c compute the buoyancy flux, B c bt = alpha*hflux/(rho*cp) bs = beta*emp*S0 B = alpha*hflux/(rho*cp) + beta*emp*S0 c c compte the deltarho of the outflow for this B c deltar = sqrt(rho*2*f*B*ams/(g*(dsill - hupper)**2)) c c compte the volume transport c q = (g*deltar/rho)*((dsill-hupper)**2)/(2.*f) c deltat = (1/q)*ams*hflux/(rho*cp) deltas = ams*emp*S0/q c c write (6,201) deltar, q, deltat, deltas c 201 format (1x,///,1x,' deltar, q, delt, dels are ...', 7e12.2) c c This version assumes maximal exchange. c h2 = (dsill - hupper)/2. h1 = dsill - h2 c t2 = t1 + deltat s2 = s1 + deltas c u2 = q/(width*h2) tran2 = q c sig02 = sigma(s2,t2,0.) c c compute speed and thickness of the source (depends upon slope1) c u1 = -u2*(h2/h1) tran1 = -q tran2 = q c c Similarly, the layer mid-depths are: c d1 = h1/2. d2 = dsill - 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) = width 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 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 geostrophic 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 ! big, but based upon 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 write (6,3333) E, deltad, ugeo, bflux, broad, phi 3333 format (1x,/,1x, 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 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 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