[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Biotic HOWTO



- ----------
X-Sun-Data-Type: text
X-Sun-Data-Description: text
X-Sun-Data-Name: text
X-Sun-Charset: us-ascii
X-Sun-Content-Lines: 23

OCMIP folks,

Attached to this message is code for constructing the nutrient restoring 
model.  Everything is there except for CO2 equilibria routines.  These
will become available with the C14 HOWTO.

Look at main.f first.

To access the phosphate maps:

ftp ftp.essc.psu.edu
login anonymous
cd /pub/meteo/ferial
mget po4map.dat.gz README readmap.fr nutrient-doc

Some details need to be specified as well, such as how long to run the model
and what to save.  These will be forthcoming shortly.

The HOWTO is preliminary; we'll finalize it in one week based on your
comments.  I don't expect it to change very much, though, since we've been
through much of the details already.

Ray
- ----------
X-Sun-Data-Type: h-file
X-Sun-Data-Description: h-file
X-Sun-Data-Name: bio.h
X-Sun-Charset: us-ascii
X-Sun-Content-Lines: 5

      parameter (km=??, nt=?, kc=2, kcp1=kc+1)
      common /bio/ jpo4(km), jdop(km), jo2(km), jdic(km), jalk(km),
     &  po4star(kc), zforg(km), zfca(km), kmt, tb(km,5), dz(km),
     &  z(0:km), tau, sigma, kappa, rnp, rcp, ro2p, o2min, r, zc, a, d
      real kappa, jpo4, jdop, jo2, jdic, jalk
- ----------
X-Sun-Data-Type: default
X-Sun-Data-Description: default
X-Sun-Data-Name: co2flux.f
X-Sun-Charset: us-ascii
X-Sun-Content-Lines: 67

      subroutine co2flux(t,s,kw660,ppo,dic,alk,po4,si,xco2,dz1,co2ex)
c
c**********************************************************************
c
c  Computes the time rate of change of DIC in the surface
c  layer due to air-sea gas exchange in mol/m^3/s.
c
c  Inputs:
c    t        model surface temperature (deg C)
c    s        model surface salinity (permil)
c    kw660    gas transfer velocity at a Schmidt number of 660, accounting
c               for sea ice fraction (cm/hr)
c    ppo      surface pressure divided by 1 atm.
c    dic      surface DIC concentration (mol/m^3)
c    alk      surface alkalinity (eq/m^3)
c    po4      surface phosphate concentration (mol/m^3)
c    si       surface silicate concentration (mol/m^3)
c    xco2     atmospheric CO2 mixing ratio (ppm)
c    dz1      surface grid box thickness (m)
c  Output:
c    co2ex    time rate of change of DIC in the surface layer due
c               to air-sea exchange (mol/m^3/s)
c
c  Two functions are called:
c    scco2    Schmidt number of CO2
c    co2sato  CO2 saturation concentration at 1 atm (mol/m^3)
c  One subroutine is called:
c    co2calc  Computes actual aqueous CO2 concentration in surface
c               layer (mol/m^3). 
c
c  Numbers in brackets refer to equation numbers in simulation design
c  document.
c
c  Ray Najjar 1/29/99
c
c**********************************************************************
c
      real kwco2
c
c  Compute the transfer velocity for CO2 in m/s [4]
c
      kwco2 = (scco2(t)/660)**-0.5*0.01/3600.0
c
c  Do the carbonate equilibrium calculations.  Pick values of
c  ph range.  You may want to change these per suggestions of Sabine
c  and Key.  Note also that co2calc subroutine expects concentrations in
c  mol/kg and eq/kg, so average surface density is used to convert.
c
      phlo = 6.0
      phhi = 9.0
      dicin = dic/rhoavg
      alkin = alk/rhoavg
      po4in = po4/rhoavg
      siin = si/rhoavg
c
      call co2calc(t,s,dicin,alkin,po4in,siin,phlo,phi,ph,co2star)
c
c  Compute the saturation concentration for CO2 [3]
c
      co2sat = co2sato(t,s,xco2)*ppo
c
c  Compute time rate of change of CO2 due to gas exchange [1]
c
      co2ex = kwco2*(co2sat-co2star)/dz1
c
      return
      end
- ----------
X-Sun-Data-Type: default
X-Sun-Data-Description: default
X-Sun-Data-Name: jbio.f
X-Sun-Charset: us-ascii
X-Sun-Content-Lines: 115

      subroutine jbio
c
c*********************************************************************
c
c  Subroutine for computing jpo4, jdop, jo2, jdic and jalk in units of
c  mol/m^3/s.  Numbers in brackets refer to equations in simulation
c  design document.
c
c*********************************************************************
c
c  Variables in common with main program
c
      include 'bio.h'
c
      real jprod(kc), jca(km), f(km), fca(km)
c
c  Compute phosphate uptake in surface layers as a result of
c  restoring towards observations.  Allow for likelihood
c  that grid box interface will not lie exactly at compensation
c  depth. [24]
c
      do k=1,kc
        if (tb(k,1) .gt. po4star(k)) then
          jprod(k) = (tb(k,1) - po4star(k))/tau
        else
          jprod(k) = 0.0
        endif
      enddo
      jprod(kc) = jprod(kc)*(zc - z(kc-1))/dz(kc)
c
c  Compute DOP source/sink term everywhere. [25]
c
      do k=1,kc
        jdop(k) = sigma*jprod(k) - kappa*tb(k,2)
      enddo
c
      do k=kcp1,kmt
        jdop(k) = -kappa*tb(k,2)
      enddo
c
c  Compute particulate P flux at compensation depth. [26]
c
      fc = 0.0
      do k=1,kc
        fc = fc + (1.0-sigma)*jprod(k)*dz(k)
      enddo
c
c  Do the same for the bottom of all grid boxes below compensation
c  depth. [27]
c 
      do k=kc,kmt
        f(k) = fc*zforg(k)
      enddo
c
c  Compute phosphate source/sink term.  Make sure material hitting
c  bottom is remineralized in bottom box. [28] (Note that the sigma
c  in equation 28a in the document should be removed.)
c
      do k=1,kc
        jpo4(k) = -jprod(k) + kappa*tb(k,2)
      enddo
c
      jpo4(kc) = jpo4(kc) + (fc - f(kc))/dz(kc)
c
      do k=kcp1,kmt
        jpo4(k) = (f(k-1) - f(k))/dz(k) + kappa*tb(k,2)
      enddo
      jpo4(kmt) = jpo4(kmt) + f(kmt)/dz(kmt)
c
c  Compute oxygen source/sink term. [30]
c
      do k=1,kmt
        if (tb(k,3) .gt. o2min) then
          jo2(k) = -ro2p*jpo4(k)
        else
          jo2(k) = 0.0
        endif
      enddo
c
c  Compute formation of calcium carbonate above compensation depth [31a]
c  Note that a minus sign is missing in front of the right hand side of
c  Equation 31a in the document.
c
      do k=1,kc
        jca(k) = -r*rcp*(1.0-sigma)*jprod(k)
      enddo
c
c  Compute flux of calcium carbonate at compensation depth and at 
c  bottom of all grid boxes below compesation depth [32]
c
      fcca = r*rcp*fc
c
      do k=kc,kmt
        fca(k) = fcca*zfca(k)
      enddo
c
c  Compute dissolution of calcium carbonate below compensation
c  depth. [31b]
c
      jca(kc) = jca(kc) + (fcca - fca(kc))/dz(kc)
c
      do k=kcp1,kmt
        jca(k) = (fca(k-1) - fca(k))/dz(k)
      enddo
      jca(kmt) = jca(kmt) + fca(kmt)/dz(kmt)
c
c  Compute source/sink terms for DIC and Alk.  [33] and [34].
c
      do k=1,kmt
        jdic(k) = rcp*jpo4(k) + jca(k)
        jalk(k) = -rnp*jpo4(k) + 2.0*jca(k)
      enddo
c
      return
      end
- ----------
X-Sun-Data-Type: default
X-Sun-Data-Description: default
X-Sun-Data-Name: main.f
X-Sun-Charset: us-ascii
X-Sun-Content-Lines: 164

      program main
c
c****************************************************************************
c
c  This code can serve as a template for the OCMIP nutrient-restoring model.
c  The code is based on Section 7 of the document "Design of OCMIP-2
c  simulations of chlorofluorocarbons, the solubility pump and common
c  biogeochemistry," which can be found on the OCMIP web site.  Variable
c  names here are very similar to those in the document, except where
c  noted.  You will probably want to change some variable names because
c  some of them (e.g., r and f) are probably already names of other
c  variables in your model.
c
c  I have chose clarity over efficiency in writing this code.  You should
c  modify the code as you see fit, or merely use it to guide your own
c  algorithms.
c
c  Units are all mks.  Concentrations are in mol/m^3.  rhoavg is the
c  average surface density and is used to convert from kg to m^3
c
c  The code here is only for a single water column.  You will have to
c  generalize the code to work at any latitude and longitude.
c  The maximum number of vertical grid points is equal to km.  kmt is the
c  actual number of ocean grid points in the vertical, and will vary with
c  latitude and longitude.
c
c  Tracer concentrations of the previous timestep are carried in the
c  array tb(km,nt), where nt = 5.  Indexing is as follows:
c    n = 1, PO4
c    n = 2, DOP
c    n = 3, O2
c    n = 4, DIC
c    n = 5, Alk
c
c  The code presented here will:
c    (1) Compute the time rate of change of phosphate, DOP, oxygen,
c        DIC and alkalinity due to biological processes:  jpo4, jdop,
c        jo2, jdic, jalk (mol/m^3/s)
c    (2) Compute the time rate of change of oxygen and DIC in the
c        surface box due to air-sea gas exchange:  o2ex and co2ex
c        (mol/m^3/s)
c
c  kc is the index of the grid box that includes the compensation depth
c  (75 m).  You will have to specify kc yourself.  If you have an
c  interface that lies exactly at 75 m, then choose kc to correspond to
c  the grid box just above the compensation depth.
c
c  Calling sequence is:
c
c  main ---> jbio
c    \
c     \----> o2flux ---> sco2 and o2sato
c      \
c       ---> co2flux ---> sco2, co2sato and co2calc
c
c  There is also an include file, bio.h, for main and jbio
c
c  Ray Najjar 1/29/99
c
c**************************************************************************
c
c  Variables in common with subroutine jbio
c
      include 'bio.h'
c
      real kw660
c
c  Set constants
c
      rhoavg = 1024.5
      sperd = 24.0*3600.0
      spery = 365.25*sperd
c
      tau = 30.0*sperd
      sigma = 0.67
      kappa = 1.0/0.5/spery
c
      rnp = 16.0
      rcp = 117.0
      ro2p = 170.0
c
      o2min = 4.0*rhoavg*1.0e-6
      d = 3500.0
      r = 0.07
      a = 0.9
      zc = 75.0
      xco2 = 278.0
      siavg = 7.5*rhoavg*1.0e-6
c
c  Convert initial concentrations from umol/kg to mol/m^3
c
      po4avg = 2.17*rhoavg*1.0e-6
      dopavg = 0.02*rhoavg*1.0e-6
      o2avg = 170.0*rhoavg*1.0e-6
      dicavg = 2230.0*rhoavg*1.0e-6
      alkavg = 2370.0*rhoavg*1.0e-6
c
c  dz(k) is the grid box thickness in your model.  k = 1 is the surface
c  box, k = km is the deepest box.  z(k) is the depth of the bottom of
c  grid box k.
c
      z(0) = 0.0
      do k=1,km
        z(k) = z(k-1) + dz(k)
      enddo
c
c  Compute arrays used for computing sinking fluxes.  zforg(k) and zfca(k)
c  are the sinking fluxes of organic matter and calcium carbonate, 
c  respectively, at the bottom of grid box k, divided by the flux at the
c  compensation depth.
c
      do k=1,km
        zforg(k) = (z(k)/zc)**-a
        zfca(k) = exp(-(z(k)-zc)/d)
      enddo
c
c  Initialization
c
      do k=1,km
        tb(k,1) = po4avg
        tb(k,2) = dopavg
        tb(k,3) = o2avg
        tb(k,4) = dicavg
        tb(k,5) = alkavg
      enddo
c
c  Need code here for reading in monthly maps of climatological fields
c  interpolated to your model grid:
c    ppo     surface pressure normalized to 1 atm
c    kw660   gas transfer velocity at a Schmidt number of 660, accounting
c              for sea ice fraction 
c    po4star  observed phosphate distribution
c
c  Start timestepping.  Need code here to linearly interpolate above maps
c  to model timestep
c
c  Compute time rates of change of PO4, DOP, O2, DIC and Alk due to
c  biological processes.  Set to zero if there are no ocean grid boxes
c  completely below the compensation depth.
c
      if (kmt .gt. kc) then
        call jbio
      else
        do k=1,km
          jpo4(k) = 0.0
          jdop(k) = 0.0
          jo2(k) = 0.0
          jdic(k) = 0.0
          jalk(k) = 0.0
        enddo
      endif
c
c  Compute time rates of change of oxygen and DIC in the surface layer
c  due to air-sea gas exchange (o2ex and co2ex).
c
      call o2flux(t,s,kw660,ppo,tb(1,3),dz(1),o2ex) 
c
      call co2flux(t,s,kw660,ppo,tb(1,4),tb(1,5),tb(1,1),siavg,
     &  xco2,dz(1),co2ex)
c
c  Update tracers (your code)
c
      stop
      end
- ----------
X-Sun-Data-Type: default
X-Sun-Data-Description: default
X-Sun-Data-Name: o2flux.f
X-Sun-Charset: us-ascii
X-Sun-Content-Lines: 46

      subroutine o2flux(t,s,kw660,ppo,o2,dz,o2ex)
c
c**********************************************************************
c
c  Computes the time rate of change of oxygen in the surface
c  layer due to air-sea gas exchange in mol/m^/s.
c
c  Inputs:
c    t       model surface temperature (deg C)
c    s       model surface salinity (permil)
c    kw660   gas transfer velocity at a Schmidt number of 660, accounting
c              for sea ice fraction (cm/hr)
c    ppo     surface pressure divided by 1 atm.
c    o2      surface ocean O2 concentration (mol/m^3)
c    dz      thickness of surface grid box (m)
c  Output:
c    o2ex    time rate of change of oxygen in the surface layer due
c              to air-sea exchange (mol/m^3/s)
c
c  Two functions are called:
c    sco2    Schmidt number of oxygen
c    o2sato  oxygen saturation concentration at 1 atm (mol/m^3)
c
c  Numbers in brackets refer to equation numbers in simulation design
c  document.
c
c  Ray Najjar 1/29/99
c
c**********************************************************************
c
      real kwo2
c
c  Compute the transfer velocity for O2 in m/s [4]
c
      kwo2 = (sco2(t)/660)**-0.5*0.01/3600.0
c
c  Compute the saturation concentrations for O2 [3]
c
      o2sat = o2sato(t,s)*ppo
c
c  Compute time rate of change of O2 due to gas exchange [1]
c
      o2ex = kwo2*(o2sat-o2)/dz(1)
c
      return
      end
- ----------
X-Sun-Data-Type: default
X-Sun-Data-Description: default
X-Sun-Data-Name: o2sato.f
X-Sun-Charset: us-ascii
X-Sun-Content-Lines: 46

      function o2sato(T,S)
c
C ********************************************************************
C                                     
C Computes the oxygen saturation concentration at 1 atm total pressure
c in mol/m^3 given the temperature (t, in deg C) and the salinity (s,
c in permil). 
C
C FROM GARCIA AND GORDON (1992), LIMNOLOGY and OCEANOGRAPHY.
C THE FORMULA USED IS FROM PAGE 1310, EQUATION (8).
c
C *** NOTE: THE "A3*TS^2" TERM (IN THE PAPER) IS INCORRECT. ***
C *** IT SHOULDN'T BE THERE.                                ***
C
C o2sato IS DEFINED BETWEEN T(freezing) <= T <= 40(deg C) AND
c 0 permil <= S <= 42 permil
C C
C CHECK VALUE:  T = 10.0 deg C, S = 35.0 permil, 
c o2sato = 0.282015 mol/m^3
C
C ********************************************************************
c
      DATA A0/ 2.00907   /,A1/ 3.22014   /, A2/ 4.05010 /,
     $     A3/ 4.94457   /,A4/-2.56847E-1/, A5/ 3.88767 /
      DATA B0/-6.24523E-3/,B1/-7.37614E-3/,
     $     B2/-1.03410E-2/,B3/-8.17083E-3/
      DATA C0/-4.88682E-7/
C
      TT  = 298.15-T
      TK  = 273.15+T
      TS  = LOG(TT/TK)
      TS2 = TS**2
      TS3 = TS**3
      TS4 = TS**4
      TS5 = TS**5
      CO  = A0 + A1*TS + A2*TS2 + A3*TS3 + A4*TS4 + A5*TS5
     $     + S*(B0 + B1*TS + B2*TS2 + B3*TS3)
     $     + C0*(S*S)
      o2sato = EXP(CO)
c
c  Convert from ml/l to mol/m^3
c
      o2sato = o2sato/22391.6*1000.0
C
      RETURN
      END
- ----------
X-Sun-Data-Type: default
X-Sun-Data-Description: default
X-Sun-Data-Name: scco2.f
X-Sun-Charset: us-ascii
X-Sun-Content-Lines: 14

      function scco2(t)
c
c*********************************************************************
c
c  Computes the Schmidt number of CO2 in seawater using the 
c  formulation presented by Wanninkhof (1992, J. Geophys. Res., 97,
c  7373-7382).  Input is temperature in deg C.
c
c*********************************************************************
c
      scco2 = 2073.1 - 125.62*t + 3.6276*t**2 - 0.043219*t**3
c
      return
      end
- ----------
X-Sun-Data-Type: default
X-Sun-Data-Description: default
X-Sun-Data-Name: sco2.f
X-Sun-Charset: us-ascii
X-Sun-Content-Lines: 14

      function sco2(t)
c
c*********************************************************************
c
c  Computes the Schmidt number of oxygen in seawater using the
c  formulation proposed by Keeling et al. (1998, Global Biogeochem.
c  Cycles, 12, 141-163).  Input is temperature in deg C.
c
c*********************************************************************
c
      sco2 = 1638.0 - 81.83*t + 1.483*t**2 - 0.008004*t**3
c
      return
      end