README:  Ecosystem Modeling Testbed 
      
To use this one-dimensional ecosystem modeling testbed, you will need to 

(1) provide a subroutine (derivs.f) that computes the biological sources and
    sinks of each state variable, and 
(2) edit a subroutine (readinit.f) that assigns initial conditions to each
    state variable. I do not anticipate that you will need to edit any of
    the other subroutines; however if you do, please let me know. 

Main Program:

driver.f

This is the main program that drives all required subroutines.  You will
need to edit this to replace my home and data directory paths with your
own. The directory "homedir" should contain all model code (testbed.tar),
and the directory "datdir" should contain all the forcing and data files
(forcing.tar). The lengths of the corresponding character statements in
common_blocks.h will also need to be changed.

The vertical domain is 150 m, with a resolution of 10 m resulting in a total
of 15 depth intervals. The time step is 3600 seconds (one hour.) The Arabian
Sea simulation begins on Year Day 274 of 1994 and extends into December 1995
(YD 718, 10660 time steps) and the Equatorial Pacific simulation begins on
Year Day 274 of 1991 and extends into December 1992 (YD718, 10660 time
steps.) 

When running the program you will be required to answer four questions: 

(1) Are you applying the Arabian Sea testbed (ntb=1) or the Equatorial
    Pacific testbed (ntb=2)?

(2) How many state variables will independently undergo the physical
    processes of diffusion, advection and mixing (nsv)?  

(3) What is the rate of remineralization of detrital nitrogen (reminrate,
    [d-1]), between 150m (the bottom boundary of the model) and 800m (the
    depth of the sediment traps)?

(4) What is the average sinking velocity of detrital nitrogen (wtrap, [m
    d-1]), between 150m (the bottom boundary of the model) and 800m (the
    depth of the sediment traps)?

You will also be required to enter six vectors, each of length nsv. These
vectors will describe the sinking rate, bottom boundary condition, unit
conversion, etc for each state variable.  The example given below is for
a simple P,Z,D,N model, where the first element in each vector applies to
the phytoplankton component, the second element applies to the zooplankton
component, the third element applies to the detrital component, and the
fourth element applies to the nitrate component. 

Vertical sinking rates of each state variable:

   wnsv=[0.,0.,32.,0.]    

Note: Units must be m d-1.

Bottom boundary conditions of each state variable:

   bbcnsv=[0.,0.,-9.,16.]   

Note: Units must be consistent with derivs.f subroutine.  A value of -9
in this vector indicates that the following statement will define the bottom
boundary condition for state variable inb:   

   bbc = bio(nz,t,inb)-(bio(nz-1,t,inb)-bio(nz,istep,inb))

Horizontal advection of each state variable:

   HAnsv=[0.,0.,0.,1.] for EqPac and 
   HAnsv=[0.,0.,0.,0.] for Arabian Sea

Note: When running the testbed code in the Arabian Sea, no horizontal
advection should be applied, and each element of this vector should be
zero. When running the code in the Equatorial Pacific, assume that
horizontal advection of plankton and detritus is negligible, and set these
elements to zero. Horizontal advection of nutrients is generally not
insignificant in the Equatorial Pacific. Enter a value of "1", for nitrate,
"2" for ammonium, and "3" for iron.  Other nutrients (PO4, SiO4) can be
assumed to have a characteristic time scale for horizontal advection that is
similar to that of nitrate, and therefore a value of "1" can be used for
these nutrients as well.

Aeolian flux source of each model component:

   aeonsv=[0.,0.,0.,0.]

If an aeolian flux of iron is desired, set element for iron equal to "1".
The subroutine "readinit.f" will need to be modified to read in a time
series of aeolian flux, or the flux can be set to a constant value.

Define which state variables will be used in the cost function:

   Cnsv=[0.,2.,3.,1.]

Note: Enter a value of "1" for state variables that will be summed to
compare to the N.dat data file: nitrate/nitrite/ammonium. Enter a value of
"2" for state variables that will be summed to compare to the Z.dat data
file: all zooplankton/heterotrophs. Enter a "3" for state variables that
will be summed to compare to the ST.dat data file: all detrital nitrogen
components.  All remaining state variables should be given values of
"0". (Phytoplankton should have a value of "0" since chlorophyll, not
phytoplankton, is included in the cost function.)

Define how your specific model units should be converted to standard units: 

   unitnsv=[0.,6.625,1.,1.]

Note: Any state variable that has a "Cnsv" element that equals zero,
i.e. that is not involved in the cost function computation, should have a
corresponding "unitnsv" value equal to zero. Non-zero elements of "Cnsv"
must correspond to non-zero elements of "unitnsv".  The standard unit for N
is mmolN m-3. If this is not the unit used in derivs.f, the element of
"unitnsv" corresponding to N should be equal to this conversion factor. The
standard unit for Z is mmolC m-3. If this is not the currency used for Z in
your derivs.f subroutine, the element of "unitnsv" corresponding to Z should
be equal to this conversion factor. (In this example the model currency for
Z is mmolN m-3, and therefore unitnsv(2)= 6.625.) The standard unit for
detrital nitrogen is mmolN m-3. If this is not the currency used for D in
your derivs.f subroutine, the element of "unitnsv" corresponding to D should
be equal to this conversion factor.

Header Files:

common_blocks.h

This header file contains common statements, variable definitions, array
definitions and parameter assignments that are universally applied.

Input Subroutines:

readdata.f

This subroutine reads in the data files that are required for the cost
function computation. Data were obtained from the JGOFS website, and include
particulate organic nitrogen flux (sediment traps at 800 m),
nitrate+nitrite+ammonium (bottle data), zooplankton carbon (MOCNESS net
tows), primary production (carbon assimilation from dawn to dawn) and
chlorophyll. The sediment trap data are lagged as a function of the value
supplied for "wtrap", accounting for the time difference due to sinking from
the bottom of the model domain (150 m), to the depth of the sediment traps
(roughly 800 m).

readforcing.f

This subroutine reads in data from a sequence of forcing files, which are
required to run an ecosystem model in the testbed. These data include: PAR
[W m-2], mixed-layer depth [m], vertical velocity [m d-1], vertical
diffusivity [m2 s-1], and temperature [oC]. Temperature and PAR are
obtained from data (when available).  Vertical diffusivity is computed
following Pacanowki and Philander (1981).  Mixed-layer depth and vertical
velocity are computed from a 3-D coupled biological-physical model
(Murtugudde, pers. comm.).  PAR is converted to total solar radiation (QSW)
by dividing through by 0.4.  To conform with the testbed model, the units
for all the applied data are transformed to MKS. 


For the equatorial Pacific testbed, horizontal advection information is also
read within this subroutine.  For example, the file EP_HA_Ni.dat contains a
time series of the inverse of the characteristic time scale for horizontal
advection of ammonium (HN [s-1]), again computed from a 3-D coupled
biological-physical model (Murtugudde, pers. comm.) over a one-degree length
scale:

To complete this example, the effect of horizontal advection on nitrate
concentration is computed as (in horizadv.f):

Analogous equations hold for ammonium and iron concentration.  For other
major nutrients such as PO4 and SiO4 for which EP_HA_xx.dat data files do
not exist, we assume that horizontal advection can be approximated as:

readinit.f

This subroutine reads in initial conditions for seven state variables (eight
for the Equatorial Pacific). Initial nutrient conditions are obtained from
JGOFS Arabian Sea and EqPac Process Study Cruises: nitrate [mmolN m-3],
phosphate (reactive phosphorus) [mmolP m-3], silicate (silicic acid/reactive
silica) [mmolS m-3], and iron [umolFe m-3].  Initial chlorophyll [mg chl
m-3], zooplankton [mmolN m-3] and detrital nitrogen [mmolN m-3]
concentrations are obtained from output from a 3-D biological-physical model
(Murtugudde, pers. comm.)  

The subroutine "readinitExample.f" will have to be adapted for use
with your specific ecosystem model. Two examples are provided, both of which
use mmolN m-3 as model currency.  If alternative model currencies are used,
appropriate unit conversions will have to be made.

In the four-component example, a constant Chl:N ratio = 1.59 is used to
convert initial chlorophyll concentrations [mg m-3] to phytoplankton
concentrations  [mmolN m-3]. In the eight-component example, the initial
plankton concentrations are split 70/30 between the small and large size
classes, respectively.  Initial ammonium concentrations are set to zero in
the Arabian Sea, but are non-zero in the Equatorial Pacific.  Aeolian iron
flux, if used, is also defined within this subroutine.

Forward Model Routines:

model.f

This subroutine contains the main time loop in which the physical processes
(advection, diffusion, sinking and mixing) and the biological processes are
carried out. All (nsv) biological state variables are contained in the array
bio(z,t,inb), where z is depth, t is time, and inb is the specific number
assigned to each state variable.  In addition, two diagnostics are tracked
within this array: bio(z,t,nsv+1) is chlorophyll concentration, and
bio(z,t,nsv+2) is the rate of primary production.

For each time step, the biological sources and sinks are computed in
biosub.f.  

The subroutine then loops through all nsv state variables, applying vertical
advection, vertical diffusion, horizontal advection, sinking and mixing. 

Bottom boundary conditions are obtained from values in "bbcnsv". Aeolian
flux is added, if applicable. Vertical advection is computed using a simple
centered, finite difference, gradient-form advection scheme. Vertical
diffusion is applied using a Crank-Nicholson vertically variable diffusion
operation, with a closed upper boundary and an open bottom boundary. The
coefficients for the tridiagonal solver are set outside the state variable
loop.  At each time step, all state variables are mixed down to the bottom
of the mixed layer.
 
biosub.f

Within the main depth loop of this subroutine, there are two primary
computations. (1) The subsurface light field is calculated at each depth
using a simple two-wavelength optical model, and (2) bio(z,t,inb) is
computed from bio(z,t-1,inb) using a second-order Runge Kutta scheme where
the first time step is Eulerian. This involves two calls to the ecosystem
model subroutine, derivs.f, which computes the time change in each
biological compartment.

derivs.f

This subroutine is provided by the user of the testbed. In the example
provided (derivs_4.f), the array y(1) = bio(z,t-1,1)= phytoplankton
concentration at the previous time step and current depth and is passed from
the biosub.f subroutine. The array dydtt(1) contains the time rate of change
of phytoplankton (units mmolN m-3 s-1) and is passed from derivs.f to
biosub.f. Additionally, chlorophyll concentration and the rate of primary
production must be computed within this subroutine. The array dydtt(ichl) is
chlorophyll concentration (units must be mg C m-3) and dydtt(ipp) is the
rate of primary production occurring over the last time step (units must be
mg C m-3 s-1).

tridag.f

This subroutine solves the Crank-Nicholson scheme for vertical diffusion.

horizadv.f

Horizontal advection is applied as described above (see readforcing.f).  

sink.f

This subroutine applies a modified second order flux-corrected transport
scheme (Smolarkiewics, 1983) to calculate the effect of sinking on each
model component with a nonzero element of the vector "wnsv". 

Output Routines:

calccost.f

This subroutine computes the cost function, and writes the result to the
screen. 

output.f

This subroutine computes model equivalents of the observed data (chlorophyll
[mg Chl m-3], zooplankton [mmolC m-3], nitrate + ammonium [mmolN m-3],
primary production [mgC m-3 d-1] and export flux [mg N m-2 d-1]).  These are
written to the files PrPr.out, Z.out, etc… Time series of each state
variable that is included in the cost function are output into two separate
output files, e.g. vertical profiles are included in
Vprof_ArabSea.out/Vprof_EqPac.out and vertical averages are given in
Vavg_ArabSea.out/Vavg_EqPac.out.