c
program pwpmain
c
c
c This main program is used to drive the pwp subroutine with
c a daily cycle of heating, or, an arbitrary (observed) time
c series of surface flux data.
c
c The program 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 air/sea flux data (optional),
c 9 may be opened to read the initial t, s profile (optional).
c
c
c Written by Jim Price, april 27, 1989. Please direct any
c comments or questions to Jim Price, WHOI, Woods Hole, MA
c 02543, USA, tel. 508-289-2526.
c
c Version 2: changed the gradient Ri relaxation method.
c Version 3: clarified the bulk Ri relaxation.
c Version 4: fixed an error that disabled background diffusion.
c
c Last editted on 20 Jan, 1999 by JFP.
c
c
character*16 ifile
c
c The arrays below store the values of:
c east and north current (u, v),
c temperature and salinity (t, s),
c density (d, sigma units), and
c depth of the grid points, z.
c
c All units are mks, except density, d, which is in sigma
c units.
c
common t(500), s(500), u(500), v(500), d(500), absrb(500),
1 dz, dt, nz, g, ro, cpw, rg, rb, em1, em2, em3
dimension z(500)
c
c The arrays below are used to store the initial z,t,s.
c
dimension zo(500), to(500), so(500), do(500)
c
c The arrays below store times series of surface flux data.
c Depending upon the time step and length of the record you
c could easily run out of these.
c
dimension daya(2500), qia(2500), qla(2500), txa(2500),
x tya(2500), empa(2500)
c
c nzmax should match the size of the initial profile arrays
c dimensioned above, and nmet should match the flux arrays.
c
nzmax = 500
nmet = 2500
c
write (6,800)
800 format (
x 1x, 'This program provides a means to run the pwp',/,
x 1x, 'upper ocean model. To run this model you must',/,
x 1x, 'prescribe the air/sea fluxes, the initial',/,
x 1x, 'temperature and salinity profile, and several',/,
x 1x, 'parameters defining the site and the model.',/,
x 1x, 'The site and model parameters are listed below,',/,
x 1x, 'where default values and units given in parentheses.',/,
x 1x, 'To choose the default enter a comma.',/)
write (6,802)
802 format (
x 1x, 'beta1, longwave extinction coefficient (0.6 m)',/,
x 1x, 'beta2, shortwave extinction coefficient (20.0 m)',/,
x 1x, 'rlat, latitude (31.0)',/,
x 1x, 'udrag, decay time scale of current (9999. days)',/,
x 1x,/,
x 1x, 'dt, time step (1000. sec)',/,
x 1x, 'days, the number of days to run (1.0 day)',/,
x 1x, 'dz, grid interval (1.0 m)',/,
x 1x,/,
x 1x, 'rb, critical bulk richardson number (0.65)',/,
x 1x, 'rg, critical gradient richardson number (0.25)',/,
x 1x, 'em1, coefficient of buoyancy production (0.)',/,
x 1x, 'em2, coefficient of shear production (0.)',/,
x 1x, 'em3, coefficient of ustar**3 (0.)',/)
c
c
c Initialize rlat, the latitude.
c
rlat = 31.
c
c Initialize the solar radiation absorption;
c beta1 is the non-penetrating part of the insolation,
c beta2 is the penetrating (shortwave) part.
c values below are appropriate for fairly clear waters.
c
beta1 = 0.60
beta2 = 20.
c
c Initialize udrag, the e-folding time (days) of current
c due to unresolved processes.
c
udrag = 9999.
c
c
write (6,888)
888 format (1x,'Enter beta1, beta2, rlat, udrag.')
read (5,*) beta1, beta2, rlat, udrag
c
c
c
c dt and dz are the time step (sec) and the vertical grid
c interval (m), and days is the number of days to run
c the integration.
c
dt = 1.0e3
dz = 1.0
days = 1.
c
write (6,706)
706 format (1x,/,1x,'Enter dt, days, dz.')
read (5,*) dt, days, dz
c
dtd = dt/8.64e4
f = 2.*7.292e-5*sin(rlat*3.141/180.)
c
c
c Model constants are rb, the critical bulk richardson number,
c rg, the critical gradient richardson number, and emn, the
c coefficients that multiply the energy source terms in
c the turbulent energy budget. Values set below are the PWP
c model values.
c
rb = 0.65
rg = 0.25
em1 = 0.
em2 = 0.
em3 = 0.
c
write (6,772)
772 format (1x,/,1x,'Enter rb, rg, em1, em2, em3.')
read (5,*) rb,rg,em1,em2,em3
c
c
c
c Set up to input air/sea flux data.
c
write (6,457)
457 format (1x,//,
x 1x,'You can prescribe the air/sea flux data by:',/,
x 1x,'specifying the parameters of an idealized diurnal',/,
x 1x,'cycle from terminal input (enter 1 (default)), or, ',/,
x 1x,'read the air/sea flux data from a disk file (enter 2).')
c
metu = 1
read (5,*) metu
c
if(metu.eq.1) then
c
write (6,2790)
2790 format (1x,/,1x,'In order to define the air-sea fluxes ',/,
x 1x, 'for a diurnal cycle, the program will now ',/,
x 1x, 'query for the values of the following variables, all ',/,
x 1x, 'of which may be defaulted to the values in parens.',/,
x 1x, '(which are from day 131 of the PWP paper):',/)
write (6,801)
801 format (
x 1x, 'qimax, noon amplitude of insolation (978. w/m**2)',/,
x 1x, 'ql, steady heat loss (-126. w/m**2)',/,
x 1x, 'tx, steady east wind stress (0.07 pa)',/,
x 1x, 'emp, evaporation minus precipitation (0.0 m/sec)',/,
x 1x, 'pqfac, duration of insolation (10.0 hours)',/)
c
c Initialize the amplitude of solar insolation (qi), the
c heat loss (ql), and the wind stress (tx, uses east only).
c
qimax = 978.
ql = -126.
tx = 0.07
emp = 0.
c
c Initialize pqfac, the duration of daylight in hours.
c
pqfac = 10.
c
c
write (6,777)
777 format (1x,'Enter qimax, ql, tx, emp, pqfac.')
read (5,*) qimax,ql,tx,emp,pqfac
c
c Convert pqfac to a non-d form.
c
pqfac = pqfac/12.
if(pqfac.lt.0.00001) then
pqfac = 0.00001
end if
end if
c
c Come here if you intend to read a time series of flux data from
c a disk file.
c
if(metu.eq.2) then
c
write (6,117)
117 format (1x,/,1x,'Enter the name of the the air/sea flux file.')
read (5,116) ifile
116 format (a16)
open (unit=8,file=ifile)
c
c
do 233 j=1,nmet
read (8,*,end=459) daya(j), qia(j), qla(j), txa(j), tya(j),
x empa(j)
c
c
if(daya(j).lt.-9.) go to 459
c
nmeto = nmeto + 1
c
if(nmeto.lt.25) then
write (6,8838) nmeto, daya(j), qia(j), qla(j), txa(j), tya(j),
x empa(j)
8838 format (1x,'First 25 flux data:', i4, f7.2,2f7.0,2f7.2,e10.3)
end if
c
c
233 continue
459 continue
c
c Set nmet equal to actual number of observations
c
nmet = nmeto
c
write (6,8833)
8833 format (1x,/,1x,'Enter the starting time (day).')
read (5,*) day0
c
end if
c
c Finished with specifying air/sea flux data.
c
c
c
c Initialize the z, t, s, d 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., 20., 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 model.
c
write (6,458)
458 format (1x,//,
x 1x,'The initial temperature and salinity profile can',/,
x 1x,'be defined by:',/,
x 1x,'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
if(mzts.eq.1) then
nts = 4
zo(1) = 0.
zo(2) = 10.
zo(3) = 11.
zo(4) = 250.
to(1) = 20.
to(2) = 20.
to(3) = 10.
to(4) = 10.
so(1) = 36.
so(2) = 36.
so(3) = 36.
so(4) = 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) ifile
open (unit=9,file=ifile)
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)
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
nts = nts + 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 Compute the size of the arrays actually used.
c
nzo = int(zo(nts)/dz)
nz = nzo
if(nz.gt.nzmax) nz = nzmax
c
c Interpolate the initial data to get a profile at dz resolution.
c
do 52 j=1,nz
z(j) = float(j-1)*dz
call intx(zo,to,nzs,nts,z(j),t(j))
call intx(zo,so,nzs,nts,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)
54 format (1x,'Density prof.:',i4,2x,6f10.3)
52 continue
c
c Initialization of z, t, s, d is finished.
c
c
c
c Section below allows for application of a simple "background"
c diffusion to be applied to the profiles. this is not part of
c the pwp model per se, but is included here for the hell of it.
c
rkz = 0.
c
rkzmax = (dz**2)/(2.*dt)
write (6,3355) rkzmax
3355 format (1x,//,1x
x 'If you want to have a background diffusion, then',
x /,1x,'enter Kz (m**2/sec) (default = 0; max =',e10.3,')')
read (5,*) rkz
dstab = dt*rkz/dz**2
if(dstab.gt.0.5) then
write (6,8433) dstab
8433 format (1x,/,1x,'Warning, this value of kz will be unstable:',
x /,1x,'dstab = ', e10.3)
end if
c
c
write (6,8899)
8899 format (1x,/,1x,'Do you want a screen listing of data ?',/,
x 1x,'(0 = no, n = every nth step (default = 2)')
iwrt = 2
read (5,*) iwrt
c
c Section below stores some data onto a disk file for later
c analysis or plotting.
c
istor = 0
write (6,3186)
3186 format (1x,/,1x,'Do you want to store the data on disk ?',/,1x,
x '(0 = no (default), 1 = yes)')
read (5,*) istor
c
if(istor.eq.1) then
write (6,1188)
1188 format (1x,/,1x,'Enter the name of the file to store data.')
read (5,116) ifile
open (unit=7,file=ifile)
c
itfreq = 2
izfreq = 5
c
write (6,3187)
3187 format (1x,/,1x,'Enter the decimation rate for time and depth.',
x /,1x,'(default is 2 and 5)')
read (5,*) itfreq,izfreq
end if
c
c At this point, we are ready to call the initialization
c entry point in the pwp subroutine.
c
call pwpint(beta1,beta2,rlat,udrag)
c
c
c Time step the model for up to nnn steps.
c
nnn = 5000
c
c nnn is set to 5000 to avoid running forever
c in case the limit timed > days is missed.
c
do 40 m=1,nnn
c
c Compute the solar insolation as a function of the time,
c where timed is time since start in days.
c
c
timed = float(m-1)*dtd
if(timed.gt.days) go to 499
c
if(metu.eq.1) then
p2 = 3.141*2.
tims = amod(timed,1.)
tims = tims - 0.50
qi = 0.
if(cos(p2*tims).gt.0.) then
qi = qimax*cos(tims*p2/pqfac)
if(qi.lt.0.) qi = 0.
end if
emp = 0.
end if
c
if(metu.eq.2) then
c
day = timed + day0
call intx (daya,qia,nzs,nmet,day,qi)
call intx (daya,qla,nzs,nmet,day,ql)
call intx (daya,txa,nzs,nmet,day,tx)
call intx (daya,tya,nzs,nmet,day,ty)
call intx (daya,empa,nzs,nmet,day,emp)
c
end if
c
c nn is the number of time steps actually taken.
c
nn = nn + 1
c
c Time step the model by an increment of time, dt.
c
call pwpgo(qi,ql,emp,tx,ty,dml,dtl)
c
c
c Apply a "background" diffusion if rkz is non-zero.
c
if(rkz.gt.0) then
c
call diffus (rkz,nz,dz,dt,d)
call diffus (rkz,nz,dz,dt,s)
call diffus (rkz,nz,dz,dt,t)
call diffus (rkz,nz,dz,dt,u)
call diffus (rkz,nz,dz,dt,v)
c
end if
c
c Store the density profile if this is the first step.
c
if(m.eq.1) then
do 2278 j=1,nz
do(j) = d(j)
2278 continue
end if
c
c
c This is the place to compute diagnostics.
c
c Find max/min values on the first day.
c
if(m.eq.1) dmlmin = 100.
if(m.eq.1) tmin = t(1)
if(timed.lt.1.) then
if(t(1).lt.tmin) tmin = t(1)
if(t(1).gt.tmax) tmax = t(1)
spd = sqrt(u(1)**2 + v(1)**2)
if(spd.gt.spdmax) spdmax = spd
if(dml.lt.dmlmin) dmlmin = dml
end if
c
c
c Write out a few things on every iwrtth time step.
c
if(iwrt.lt.1) go to 1339
c
if(m.eq.1) write (6,444)
444 format (1x,/,1x,'model solution variables are:',/,1x,
x ' m timed qi ql tx t(1) s(1) u(1) ',
x ' v(1) dml dtl',/,1x)
c
if(mod(m,iwrt).ne.0) go to 1339
write (6,133) m,timed,qi,ql,tx,t(1),s(1),
c u(1),v(1),dml,dtl
133 format (1x,i4,f6.2,2f7.0,f7.2,4f7.2,2f6.0)
c
1339 continue
c
c
c Compute the net heat flux as a check on heat conservation.
c
qnet = qnet + dt*(qi + ql)
c
c
c Store data on disk, logical unit 7, if you told it to.
c
if(istor.eq.1) then
if(mod(m,itfreq).eq.0) then
write (7,3188) timed,qi,ql,tx,t(1),s(1),u(1),v(1),
x dml
3188 format (1x,20f10.3)
end if
end if
c
40 continue
499 continue
c
c
write (6,870) tmin,tmax,dmlmin,spdmax
870 format (1x,/,1x,'tmin, tmax, dmlmin, spdmax are',4f8.2)
c
c
c Compute heat and potential energy diagnostics.
c dpe will be the change in potential energy over the
c integration period, and heat is the change in heat
c content. heat should nearly equal the time-integrated
c surface heat flux, qnet, if heat has been conserved.
c
dpe = 0.
heat = 0.
c
do 8376 j=1,nz
call intx(zo,to,nzs,nts,z(j),ti)
tanom = t(j) - ti
heat = heat + tanom*dz
if(j.eq.1) tanoms = tanom
c
danom = d(j) - do(j)
dpe = dpe - g*danom*z(j)*dz
c write (6,4465) j,do(j),d(j),danom,dpe
4465 format (1x,'do,d,danom,dpe',i4,3f10.3,e12.2)
8376 continue
c
qnet = qnet/(ro*cpw)
write (6,3877) heat, qnet, dpe
3877 format (1x,/,1x,'heat, qnet, dpe are', 3f9.2)
c
c End of the diagnostics.
c
end
c
c
subroutine pwpint (beta1,beta2,rlat,udrag)
c
c
c This subroutine is an implementation of the Price, Weller,
c Pinkel upper ocean model described in JGR 91, C7 8411-8427
c July 15, 1986 (PWP). This version was coded and documented by
c Jim Price in April 1989.
c
c Edited on 20 September 1993 by JFP to allow for a critical
c gradient Richardson number other than 1/4, and to implement a
c different and a priori better means of achieving convergence
c of gradient ri mixing (see subroutine stir for details).
c The major difference is that the revised scheme gives a more
c smoothly varying mixed layer depth over a diurnal cycle.
c Edited on 14 December 1998 by JFP to clarify the bulk Ri
c relaxation.
c
c This model also implements an energy budget form of an
c entrainment parameterization that is very similar to that
c described in Price, Mooers, Van Leer, JPO 8, 4, 582-599 (and
c references therein). This part of the model should be
c treated as developmental only, as there are several features
c that are arbitrary to this model (i.e., the depth of the ml
c during times of heating can be the grid interval, dz). To
c use this parameterization set the bulk Richardson number to
c zero, and set em1, or em2, or em3 to non-zero. The gradient
c Richardson number can be zero or not.
c
c
c Please direct questions and comments to Jim Price, WHOI,
c Woods Hole, MA, 02543, MA, tel. 508-548-1400, x2526.
c
c All units within the model are mks, except density, d,
c which is in sigma units.
c
c
common t(500), s(500), u(500), v(500), d(500), absrb(500),
1 dz, dt, nz, g, ro, cpw, rg, rb, em1, em2, em3
c
c
c nz is the number of grid points in the vertical.
c dz is the grid interval (meters).
c beta1 and beta2 are the longwave and shortwave extinction
c coefficients for solar insolation (meters).
c rlat is the latitude (deg).
c tr, sr are the reference values of temperature and
c salinity used in the density equation (C, and ppt).
c udrag is the e-folding time scale of current due to
c unresolved processes (days) (set to infinite if > 100).
c rg is the critical gradient richardson number (= 1/4,
c in the pwp model).
c rb is the critical bulk richardson number (= 0.65,
c in the pwp model).
c em1, em2, em3 are coefficients that multiply energy
c source terms in the turbulent kinetic energy budget:
c buoyancy production, shear production, and ustar**3; all
c are zero in the pwp model, but can be used here by setting
c rb = 0, and then at least one of the emns to non-zero.
c There are some peculiar features of any energy budget model
c that arise when the ocean is heated, so be careful !
c
save d2r, hcon, ucon, f, tr, sr, dr, alpha, beta, energy
c
c Set a few constants that will not change during the run.
c
d2r = 2.*3.14159/360.
g = 9.8
ro = 1.024e3
cpw = 4183.3
hcon = dz*ro*cpw
ucon = 0.
if(udrag.lt.100.) ucon = 1./(8.64e4*udrag)
f = 2.*7.29e-5*sin(rlat*d2r)
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 = 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 absorption profile; absrb(j) is the fraction
c of the solar insolation that will be absorbed within
c the grid cell j.
c
call absorb(beta1,beta2)
c
c Finished with initialization.
c
c
return
c
c
c Come here to step ahead by an amount dt.
c
entry pwpgo (qi,ql,emp,tx,ty,dml,dtl)
c
c Input to the subroutine are the following:
c
c qi and ql are the insolation and heat loss (qi > 0,
c ql < 0 almost always) (W/m**2).
c emp is evaporation minus precipitation (m/sec)
c (emp > if evaporation exceeds precipitation).
c tx and ty are the east and north stress components (Pa).
c
c Apply heat and fresh water fluxes to the top most grid cell.
c
q1 = qi*absrb(1) + ql
t(1) = t(1) + dt*q1/(dz*ro*cpw)
c
c Absorb solar radiation at depth.
c
do 70 j=2,nz
t(j) = t(j) + qi*absrb(j)*dt/(dz*ro*cpw)
70 continue
s(1) = s(1) + s(1)*emp*dt/dz
c
c Compute the density (sigma units), and relieve static
c instability, if it occurs.
c
do 71 j=1,nz
d(j) = dr + (t(j) - tr)*alpha + (s(j) - sr)*beta
71 continue
c
call mldep(ds1)
call si(ml,sipe)
if(ml.gt.1) call mix5(ml)
call mldep(ds2)
c
c At this point the density proifile should be statically stable.
c
c Time step the momentum equation.
c
c Rotate the current throughout the water column through an
c angle equal to inertial rotation for half a time step.
c
ang = f*dt/2.
sa = sin(ang)
ca = cos(ang)
c
do 32 j=1,nz
call rot(u(j),v(j),sa,ca)
32 continue
c
c Apply the wind stress to the mixed layer as it now exists.
c
c Find the ml depth.
c
call mldep(dml)
ml = int(dml/dz)
c
du = (tx/(dml*ro))*dt
dv = (ty/(dml*ro))*dt
do 33 j=1,ml
u(j) = u(j) + du
v(j) = v(j) + dv
33 continue
c
c Apply drag to the current (this is a horrible parameterization of
c inertial-internal wave dispersion).
c
if(ucon.gt.1.e-10) then
do 37 j = 1,nz
u(j) = u(j)*(1. - dt*ucon)
v(j) = v(j)*(1. - dt*ucon)
37 continue
end if
c
c Rotate another half time step.
c
do 34 j = 1,nz
call rot(u(j),v(j),sa,ca)
34 continue
c
c Finished with the momentum equation for this time step.
c
c Cause the mixed-layer to deepen.
c
c
if(rb.le.0.00001) go to 61
c
c Come here to do the bulk Richardson number instability
c form of mixing (as in PWP).
c
nzm = nz - 1
c
rvc = rb
do 80 j = ml,nzm
delr = (d(j+1) - d(1))/ro
ds2 = (u(j+1) - u(1))**2 + (v(j+1) - v(1))**2
ds2 = ds2 + 1.e-8
h1 = float(j)*dz
rv = g*delr*h1/ds2
if(rv.gt.rvc) go to 81
c
jp = j + 1
call mix5(jp)
c
80 continue
81 continue
go to 65
c
c
61 continue
c
c
if(em1.le.0..and.em2.le.0..and.em3.le.0.) go to 65
c
c Come here to mix according to an energy balance closure;
c this implementation follows Price, Mooers, Van Leer, JPO,
c 8, 582-599, July 1978. This is not part of the PWP model
c per se, and has some rather unusual properties, so watch out !
c
c
call mldep(dml)
ml = dml/dz
taumag = sqrt(tx**2 + ty**2)
ustar = sqrt(taumag/ro)
ustarc = taumag*ustar*dt
c
c Compute the energy available to do entrainment.
c
energy = energy + em3*ustarc + em1*sipe
c
delkes = 0.
penes = 0.
c
nzm = nz - 1
do 60 l=ml,nzm
c
c Determine how much energy is required to entrain the next
c level below the base of the mixed-layer.
c
call pe(l+1,pene)
if(energy.lt.pene) go to 65
if(energy.gt.pene) then
c
c Come here if sufficient energy is avaliable to entrain at
c least one level into the mixed layer.
c
c Compute the change in kinetic energy due to entrainment.
c
eke1 = eke(l+1)
call mix5(l+1)
eke2 = eke(l+1)
delke = eke1 - eke2
c
penes = penes + pene
delkes = delkes + delke
c
c Add and subtract energy losses and gains due to entrainment.
c
energy = energy + em2*delke - pene
c
c
end if
60 continue
c
c End of the energy budget method for estimation entrainment.
c
65 continue
c
c Relieve gradient Richardson number instability if it occurs.
c
if(rg.gt.0.) call rimix(nmix,kcd)
c
call mldep(dml)
dtl = dz*float(kcd)
c
return
end
c
subroutine pe(m,pene)
c
c Compute the energy required to entrain density level
c d(m) into a mixed-layer above.
c
common t(500), s(500), u(500), v(500), d(500), absrb(500),
1 dz, dt, nz, g, ro, rg, rb, em1, em2, em3
c
dr = d(m) - d(m-1)
ho2 = dz*float(m)/2.
pene = g*dz*ho2*dr
return
end
c
subroutine rimix(nmix,kcd)
c
c This subroutine performs the gradient Richardson number
c relaxation by mixing adjacent cells just enough to bring
c them to a new richardson number = rnew (defined in
c subroutine stir). rg is the critical gradient Richardson number
c at which mixing is presumed to begin. rg = 0.25 is a nominal
c value.
c
c All of u, v, d, t, s, dz, nz are exactly as in the pwpri subroutine.
c nmix is the number of iterations required to achieve the relaxation,
c kcd is the deepest grid level reached (the base of the transition
c layer).
c
common t(500), s(500), u(500), v(500), d(500), absrb(500),
1 dz, dt, nz, g, ro, cwp, rg, rb, em1, em2, em3
c
dimension r(500)
c
rc = rg
nzm = nz - 1
kcd = 0
nmix = 0
c
j1 = 1
j2 = nzm
10 continue
c
c Compute the gradient Richardson number, taking care
c to avoid dividing by zero in the mixed-layer. The numerical
c values of the minimum allowable density and velocity differences
c are entirely arbitary, and should not effect the calculations
c (except that on some occasions thay evidently have!)
c
do 1 j = j1,j2
dd = d(j+1) - d(j)
if(dd.lt.1.e-3) dd = 1.e-3
dv = (u(j+1) - u(j))**2 + (v(j+1) - v(j))**2
if(dv.lt.1.e-6) dv = 1.e-6
r(j) = g*dz*dd/(dv*ro)
1 continue
c
c Find the smallest value of r in profile.
c
rs = r(1)
js = 1
do 2 j=2,nzm
if(r(j).lt.rs) go to 3
go to 2
3 continue
rs = r(j)
js = j
2 continue
c
c Check to see whether the smallest r is critical or not.
c
if(rs.gt.rc) go to 99
c
c Mix the cells js and js+1 that had the smallest Richardson number.
c
if(js.ge.kcd) kcd = js + 1
call stir(rc,rs,t,js)
call stir(rc,rs,s,js)
call stir(rc,rs,d,js)
call stir(rc,rs,u,js)
call stir(rc,rs,v,js)
nmix = nmix + 1
c
c Recompute the Richardson number over the part of the profile
c that has changed.
c
j1 = js - 2
if(j1.lt.1) j1 = 1
j2 = js + 2
if(j2.gt.nzm) j2 = nzm
go to 10
c
99 continue
return
end
c
c
subroutine stir (rc,r,a,j)
c
c
c This subroutine mixes cells j and j+1 just enough so that
c the Richardson number after the mixing is brought up to
c the value rnew. In order to have this mixing process
c converge, rnew must exceed the critical value of the
c richardson number where mixing is presumed to start. if
c r critical = rc = 0.25 (the nominal value), and r = 0.20, then
c rnew = 0.3 would be reasonable. If r were smaller, then a
c larger value of rnew - rc is used to hasten convergence.
c
c This subroutine was modified by JFP in Sep 93 to allow for an
c aribtrary rc and to achieve faster convergence.
c
c rc is the value of the critical gradient Richardson number,
c r is the Richardson number before this mixing,
c a is the array to be mixed,
c j and j+1 define the cell (grid level) where r was found to
c be critical.
c
dimension a(500)
c
c Set the convergence parameter.
c
rcon = 0.02 + (rc - r)/2.
rnew = rc + rcon/5. ! an arbitrary change
c
f = 1. - r/rnew
da = (a(j+1) - a(j))*f/2.
a(j+1) = a(j+1) - da
a(j) = a(j) + da
return
end
c
c
c
subroutine mldep(dml)
c
c This subroutine scans through the density array d to find
c the depth of the surface mixed-layer. The degree of density
c homogeneity, deps, is arbitary and should not effect the
c results.
c
c
common t(500), s(500), u(500), v(500), d(500), absrb(500),
1 dz, dt, nz, g, ro, cpw, rg, rb, em1, em2, em3
c
deps = 1.e-4
nzm = nz - 1
do 1 j=1,nzm
dd = abs(d(j+1) - d(j))
if(dd.gt.deps) go to 2
1 continue
2 continue
dml = float(j)*dz
return
end
c
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
c
subroutine mix5(k)
c
c This subroutine mixes arrays t, s, u, v, d down to level k.
c
common t(500), s(500), u(500), v(500), d(500), absrb(500),
1 dz, dt, nz, g, ro, cpw, rg, rb, em1, em2, em3
c
ts = 0.
ss = 0.
us = 0.
vs = 0.
ds = 0.
do 1 j=1,k
ts = ts + t(j)
ss = ss + s(j)
us = us + u(j)
vs = vs + v(j)
ds = ds + d(j)
1 continue
x = float(k)
do 2 j=1,k
u(j) = us/x
v(j) = vs/x
t(j) = ts/x
s(j) = ss/x
d(j) = ds/x
2 continue
return
end
c
c
subroutine rot(u,v,sa,ca)
c
c This subroutine rotates the vector (u,v) through an
c angle whose sine and cosine are sa and ca.
c
u0 = u
v0 = v
u = u0*ca + v0*sa
v = v0*ca - u0*sa
return
end
c
c
subroutine si(ml,sipe)
c
common t(500), s(500), u(500), v(500), d(500), absrb(500),
1 dz, dt, nz, g, ro, cpw, rg, rb, em1, em2, em3
c
dimension di(500)
c
c Find and relieve static instability that may occur in the
c density array d. This simulates free convection.
c ml is the depth of the surface mixed layer after adjustment,
c sipe is the change in potential energy due to free convection.
c
do 10 j=1,nz
di(j) = d(j)
10 continue
c
mml = 1
do 1 j=2,nz
if(d(j).gt.d(j-1)) go to 2
if(d(j).lt.d(j-1)) then
mml = j
call mix1(d,mml)
end if
1 continue
2 continue
ml = mml
c
nj = mml + 3
do 30 j=1,nj
30 continue
c
c
pesum = 0.
if(ml.ge.2) then
c
c Compute the change in potential energy.
c
do 11 j= ml,1,-1
z = float(j)*dz - dz/2.
pesum = pesum + (d(j) - di(j))*g*z*dz
11 continue
end if
sipe = pesum
c
return
end
c
c
function eke(k)
c
common t(500), s(500), u(500), v(500), d(500), absrb(500),
1 dz, dt, nz, g, ro, cpw, rg, rb, em1, em2, em3
c
c Compute the kinetic energy in the u,v profiles from the surface
c down to level k.
c
ekes = 0.
do 1 j=1,k
ekes = ekes + u(j)**2 + v(j)**2
1 continue
eke = ekes*ro*dz/2.
return
end
c
c
subroutine absorb(beta1,beta2)
c
common t(500), s(500), u(500), v(500), d(500), absrb(500),
1 dz, dt, nz, g, ro, cpw, rg, rb, em1, em2, em3
c
c Compute solar radiation absorption profile. This
c subroutine assumes two wavelengths, and a double
c exponential depth dependence for absorption.
c
c Subscript 1 is for red, non-penetrating light, and
c 2 is for blue, penetrating light. rs1 is the fraction
c assumed to be red.
c
abss = 0.
rs1 = 0.6
rs2 = 1.0 - rs1
c
do 3 j=1,nz
absrb(j) = 0.
3 continue
c
do 1 j=1,nz
z1 = float(j-1)*dz
z2 = z1 + dz
c
z1b1 = z1/beta1
z2b1 = z2/beta1
z1b2 = z1/beta2
z2b2 = z2/beta2
c
c The checks below are to avoid math overflows.
c
if(z2b1.lt.70.) absrb(j) = absrb(j) +
x rs1*(exp(-z1b1) - exp(-z2b1))
c
if(z2b2.lt.70.) absrb(j) = absrb(j) +
x rs2*(exp(-z1b2) - exp(-z2b2))
c
c Make sure that ab(z) integrates to 1 at large depth.
c
abss = abss + absrb(j)
c if(j.lt.20) write (6,2) j,z2,absrb(j),abss
c 2 format (1x,'absorption',i5,3f8.3)
1 continue
return
end
c
c
subroutine intx(r,s,nzs,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(1000), s(1000)
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
c
9 continue
return
end
c
c
c
subroutine diffus (rkz,nz,dz,dt,a)
dimension a(1000), work(1000)
c
c
c This subroutine applies a simple diffusion
c operation to the array a. It leaves the endpoints
c unchanged (assumes nothing about the
c boundary conditions).
c
c
dconst = dt*rkz/dz**2
nzm = nz - 1
c
do 1 j=2,nzm
work(j) = dconst*(a(j-1) + a(j+1) - 2.*a(j))
1 continue
c
do 2 j=2,nzm
2 a(j) = a(j) + work(j)
c
return
end
c
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