c
program Tdmain ! version of 15 March, 2013
c
c Exercise the Tdsub subroutine by defining and supplying
c the required hurricane and ocean parameters: the ocean temperature
c and salinity profiles, bottom depth, and hurricane parameters, the
c translation speed, maximum stress and radius to maximum stress.
c
c Written by Jim Price, August 2009; this Fortran version 11 March 2012
c
c Revised 15 March 2013 to make the finite differencing internally
c consistent. Also added a lower limit on the allowed storm translation
c speed.
c
c
c Please direct comments or questions to J Price, WHOI, Woods Hole, MA
c 02543, USA, tel. 508-289-2526. best is email = jprice@whoi.edu
c
c Assumes the following logical units:
c 5 is the terminal for input,
c 6 is the terminal for output,
c 7 may be opened for data storage (optional),
c 8 may be opened to read hurricane data (optional),
c 9 may be opened to read the initial t, s profile (optional).
c
character*16 ifilein, ifileout
c
c The arrays below store the values of:
c temperature and salinity (t, s),
c density, d, and depth of the grid points, z.
c
dimension zo(500), to(500), so(500)
dimension t(500), s(500), z(500), d(500)
c
c nzmax should match the size of the arrays dimensioned above
c
nzmax = 500
c
write (6,800)
800 format (
x 1x, 'This program provides a means to run the Td,',/,
x 1x, '1-d mixing model. To run this model you must',/,
x 1x, 'prescribe the temperature and salinity profiles',/,
x 1x, 'and four parameters defining the site and the hurricane.',/,
x 1x, 'Default values and units are given below in parentheses.',/,
x 1x, 'To choose the default, enter a comma.',/)
c
write (6,457)
457 format (1x,//,
x 1x,'Define the ocean site and the hurricane by:',/,
x 1x,'zbot, the bottom depth, [500] m, ',/,
x 1x,'Uh, the translation speed, [5.5] m/sec, ',/,
x 1x,'Tau, the maximum stress, [5.5] Pa,' /,
x 1x,'R, the radius to maximum stress, [35.] km' /)
c
zbot = 500. ! this large value amounts to very deep water
Uh = 5.5
Tau = 5.5
Rmax = 35.
c
write (6,777)
777 format (1x,'Enter zbot, Uh, Tau, Rmax ')
read (5,*) zbot, Uh, Tau, Rmax
c
c Convert Rmax to meters
c
Rmax = Rmax*1000.
c
c Prescribe the z, t, s profiles.
c
c Read in the data used to define the initial t,s profiles
c from logical unit nunit (read in below). The data
c is presumed to be free format, with each record being a
c triple (z, t, s). The first record should be the surface
c value, i.e., 0., 28., 36. is a reasonable surface value,
c and the last true value should be deep enough to avoid
c entering into the caluclation (unless intentional).
c The end of file is defined either by the end
c of the data (if read from a disk file), or by a z value less
c than -10. These data are then linearly interpolated to give
c profiles at dz resolution needed by the Td model.
c
write (6,458)
458 format (1x,//,
x 1x,'The initial temperature and salinity profile can',/,
x 1x,'be defined in one of three ways:',/,
x 1x,'Use data hardwired into the program (enter 1, default), or',/,
x 1x,'enter z, t, s data from the terminal (enter 2), or',/,
x 1x,'read the data from a disk file (enter 3).')
c
mzts = 1
read (5,*) mzts
c
c The parameters/profiles below are those in Price (2009).
c
h0 = 30.
z3 = 300.
dtdz = -0.05
c
if(mzts.eq.1) then
no = 3
zo(1) = 0.
zo(2) = h0
zo(3) = z3
to(1) = 29.
to(2) = 29.
to(3) = to(2) + (z3 - ho)*dtdz
so(1) = 36.
so(2) = 36.
so(3) = 36.
end if
if(mzts.eq.1) go to 53
c
if(mzts.eq.3) then
write (6,3290)
3290 format (1x,/,1x,'Enter the name of the z,t,s data file.')
read (5,116) ifilein
116 format(a15)
open (unit=9,file=ifilein)
end if
c
c
do 50 j=1,nzmax
if(mzts.eq.2) then
write (6,876)
876 format (1x,'Enter a z,t,s (end = z < -10).')
read (5,*) zo(j),to(j),so(j)
no = j-1
end if
c
if(mzts.eq.3) then
read(9,*,end=53) zo(j),to(j),so(j)
end if
c
if(zo(j).lt.-9.) go to 53
no = no + 1
c if(mzts.eq.3) write (6,503) zo(j),to(j),so(j)
c 503 format (1x,'Initial z, t, s are:', 3f10.2)
50 continue
53 continue
c
c Interpolate the initial data to get regular arrays at dz resolution.
c
dz = 1.0 ! 1 - 10 should be OK
nz = 300 ! 300 m profiles = plenty
c
do 52 j=1,nz
z(j) = float(j-1)*dz + dz/2
call intx(zo,to,no,z(j),t(j))
call intx(zo,so,no,z(j),s(j))
d(j) = sgt(t(j),s(j),gg)
c write (6,54) j,z(j),t(j),s(j),d(j)
c 54 format (1x,' z, t, s, d:',i4,2x,6f9.2)
52 continue
c
c Definition of z, t, s is finished.
c
c
c At this point, we are ready to call Tdsub.
c
call Tdsub(t, s, z, dz, nz, zbot, Uh, Tau, Rmax, Td, dmix)
c
SST0 = to(1)
write (6, 743) zbot, Uh, Tau, Rmax, dmix, SST0, Td
743 format(1x,/,1x,
c ' zbot Uh tau Rmax dmix SST0 Td ' /,
c f7.0, 2f8.2, 2x, 2f8.0, 2f8.2, //,)
c
c Store data on disk, logical unit 7, if you told it to.
c
istor = 0 ! set this to 1 or 0 if you do/do not want to save data
if(istor.eq.1) then
write(6, 223)
223 format(1x, ' Enter the file name for data storage')
read (5,116) ifileout
open(unit=7, file=ifileout)
write (7,3188) SST0, h0, dtdz, zbot, Uh, Tau, Rmax, dmix, Td
3188 format (1x,20f15.3)
end if
c
stop
end
c
subroutine Tdsub(t, s, z, dz, nz, zbot, Uh, Tau, Rmax, Td, dmix)
c
c version of 15 March 2013
c
c Evaluate the variable-depth mixing model of Price (2009),
c Ocean Sci., 5, 351–368, www.ocean-sci.net/5/351/2009/
c Please see that paper and including the Supplementary Material
c for a discussion of this model, including its interpretation and
c its inherent limitations.
c
c Written by Jim Price in 2009, this Fortran starting on 11 March, 2013.
c Revised 15 March 2013 to make the finite differencing internally
c consistent. Also added an enforced lower limit on storm translation
c speed and a nominal air/sea heat exchange.
c
c Please direct questions and comments to Jim Price, WHOI,
c Woods Hole, MA, 02543, MA, tel. 508-548-1400, x2526.
c jprice@whoi.edu
c
c
c Input:
c t, s, z are regular arrays of temperature [C], salinity and
c depth [m], depth is positive down and zero at the sea surface
c dz is the grid interval [m]
c nz is the number of grid points in the vertical
c zbot is the bottom depth [m]
c Uh is the hurricane translation speed [m/sec]
c Tau is the maximum wind stress [Pa}
c Rmax is the radius of Tau [m]
c
c Output:
c Td, the post-hurricane SST [C]
c dmix, the depth of mixing [m]
C
c
dimension t(500), s(500), d(500), z(500)
dimension istab(500)
c
c Rb is the critical bulk richardson number (= 0.65,
c from the PWP model).
c
c Set a few constants/flags that will not change during this run.
c
g = 9.8
ro = 1.024e3
Rb = 0.65
cpw = 4.1e3
c
c Qnet is a nominal, net air/sea heat exchange.
c The factor Qon is multiplied onto Qnet, so you can set
c Qon = 0 to eliminate air/sea heat exchange, or = 1 for exchange on.
c
Qnet = -5.e7 ! [J/m^2] after Price (2009)
Qon = 1.0
c
c Check Uh for very small values, i.e., a nearly stationary storm.
c Enforce a limit if Uh is less than a small value, Uhlimit, because, .
c as discussed in Price (2009) Supplement, this 1-d model does
c not account for non-local process, e.g., vertical advection, and, it
c computes storm residence time as ~ 1/Uh.
c
Uhlimit = 1.5
if (abs(Uh).le.Uhlimit) then
write (6, 21) Uhlimit
21 format (1x,//,' Note ..... Uh was reset to Uhlimit = ', f6.1)
Uh = Uhlimit
endif
c
c Evaluate alpha and beta, the thermal and haline expansion
c coefficients used in a linear equation of state.
c dr is the reference density.
c
tr = t(1)
sr = s(1)
dr = 1000. + sgt(tr,sr,gg)
alpha = sgt(tr+.5,sr,gg) - sgt(tr-.5,sr,gg)
beta = sgt(tr,sr+.5,gg) - sgt(tr,sr-.5,gg)
c
c Compute the (linearized) density profile
c
do 71 j=1,nz
d(j) = dr + (t(j) - tr)*alpha + (s(j) - sr)*beta
71 continue
c
c Compute the momentum impulse caused by the hurricane
c
Simcon = 1.2
timeh = 4.0*Rmax/Uh
U = Simcon*Tau*timeh/ro
c
c Now do the bulk Richardson number form of mixing.
c
dsum = 0.
c
c Determine the maximum depth that mixing will be allowed to go, e.g.,
c not deeper than the bottom
c
nb = abs(zbot)/dz
nm = min(nz-1, nb) ! this is as far as the mixing will go
c
c write(6,23) U, zbot, nb, nm, dr
c 23 format(1x, 'U, zbot, nb, nm, dr', 2f14.2, 2i7, f10.4)
c
do 700 j=1,nz
z(j) = abs(z(j)) ! just to be sure
700 continue
c
c The code below implements Eqn. (5) of Price (2009)
c
do 800 j = 1,nm
dsum = dsum + dz*d(j)
davg = dsum/(z(j) + dz/2)
stab = g*(d(j+1) - davg)*(z(j) + dz/2.)
Rv = (stab*((z(j) + dz/2.)**2))/(ro*(U**2))
c
c write (6,7) davg, stab, z(j), Rv
c 7 format(1x, ' davg, stab, z(j), Rv are', 5f10.3)
c
istab(j) = 0
if(Rv.le.Rb) then
istab(j) = 1
end if
c
800 continue
c
c The thickness of the upper ocean layer that is mixed by this process
c
dmix = 0.
nmix = 0
do 801 j=1,nm
dmix = dmix + float(istab(j))*abs(dz)
nmix = nmix + istab(j)
801 continue
c
call mix1(t, nmix)
c
c Subtract a nominal heat loss to the hurricane (if Qon = 1.0)
c
delTQ = Qnet/(dmix*ro*cpw)
c
Td = t(1) + Qon*delTQ ! SST after the mixing and cooling process
c
return
end
c
subroutine mix1(b,k)
c
c This subroutine homogenizes the array b down to level k.
c
dimension b(500)
bs = 0.
do 1 j=1,k
bs = bs + b(j)
1 continue
bs = bs/float(k)
do 2 j=1,k
b(j) = bs
2 continue
return
end
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 printed.
c
dimension r(500), s(500)
c
c Check for out of range ri.
c
if(ri.lt.r(1)) then
si = s(1)
write (6,20) ri, r(1)
20 format (1x,/,1x,'warning, ri was below r(1) in call to intx',/,
c 'ri and r(1) were', 2e12.4)
end if
c
c
if(ri.gt.r(n)) then
si = s(n)
write (6,21) ri, r(n)
21 format (1x,/,1x,'warning, ri was above r(n) in call to intx',/,
c 'ri and r(n) were', 2e12.4)
end if
c
c OK, 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
function sg0(s)
c
c A sigma-0 subroutine neede by the sigma-t subroutine;
c taken from seaprop.
c
c sigma-0 knudsen
c
sg0 = ((6.76786136e-6*s-4.8249614e-4)*s+0.814876577)*s
x -0.0934458632
return
end
c
c
function sgt(t,s,sg)
c
c A sigma-t subroutine taken from seaprop;
c sigma-t knudsen
c
sg = sg0(s)
20 sgt = ((((-1.43803061e-7*t-1.98248399e-3)*t-0.545939111)*t
x +4.53168426)*t)/(t+67.26)+((((1.667e-8*t-8.164e-7)*t
x +1.803e-5)*t)*sg+((-1.0843e-6*t+9.8185e-5)*t-4.7867e-3)*t
x +1.0)*sg
return
end