text: c=======================================================================
  c Soil nutrient uptake model with switches file processing. 
  c Modified from OurPrism2.f by Jonathan Rees.
  
  c User documentation is in a separate file.  Don't get caught without it.
  
  c Switch options: (default if not specified is "calculated")
  c   Fixed <value>
  c   Input
  c   Calculated
  c   Calculated <parameter>
  c   Linear <initial>
  c   Exponential <initial>
  c   Output
  
  c State variables and parameters:
  c  Ruth     Carol
  c   name     name    Units
  c  ----     ------   -----
  c  alpha    RtAbso   cm/s       Root absorbing power
  c  b        SolBuf   1          Soil buffer capacity (cm3/cm3)
  c  C0       CSurf    mol/cm3    Concentration at root surface
  c  Cav      SChem    mol/cm3    Average nutrient concentration
  c  De       DiffCo   cm2/s      Effective diffusion coefficient
  c  dt       Interval s          Duration of a timestep
  c  f        --       1          Fugacity (for computing v0) (NYI)
  c  f_slope  --       1/s        Change in f (NYI)
  c  growth   DlyGrwth 1/s        Growth rate (exponential; cm/cm/s)
  c  growth   --       cm/s       Growth rate (linear; cm/s)
  c  Km       Km       mol/cm3    Solution concentration at which net influx
  c  L        cm       cm         Root length
  c  r0       RtRad    cm         Root radius
  c  rx       RadAve   cm         Average radial distance to the next root's 
  c                                zone of influence
  c  source   --       mol/cm3/s  Constant nutrient supply/depletion (Cav pump)
  c  t        --       s          Time
  c  theta    --       1          Water content (cm3/cm3)   NYI
  c  U        UPTK     mol        Uptake (this is the main output of the model)
  c  Unew     Unew     mol        Uptake of new roots
  c  steps    --       1          Number of timesteps to run
  c  v0       VelIn    cm/s       Inward radial velocity of water at root surf
  c  Vmax     Imax     mol/cm2/s  Maximum influx at high concentration is 1/2 Vmax
  c  w0       RtWUp    cm3/s      Rate of water uptake
  c
  c Obsolescent (left over from OurPrism2):
  c  --       Cmin     mol/cm3    Concentration at which alpha goes to 0
  c  iter     Iteration --        Number of steps so far
  c  oalpha   RtAbso   cm/s       Root absorbing power (calculated iteratively)
  
  
  c To promote portability, I try to code to the FORTRAN 90 standard.  I
  c have checked conformance using MS Fortran's "Adhere to Fortran 90
  c Language Standard" option with array bounds and other checks enabled.
  c 
  c Language features used:
  c -- IF (...) THEN ... ELSEIF ... ENDIF
  c -- STATUS= in OPEN statements
  c -- Blank lines
  c -- "!" comments
  c -- Underscore as an alphabetic character.
  c 
  c Language features not used:
  c -- REAL*8 isn't official (this is an IBMish extension).  Replaced with
  c    DOUBLE PRECISION.
  c -- Statements not to go past column 72.
  c -- I assume that only the first 6 characters of an identifier are
  c    significant.
  c -- Fortran 90 probably has some structured looping constructs, but I
  c    have not bothered to learn them.
  c 
  c Style notes:
  c 
  c -- File extension ".fi" means "Fortran include" and is the extension
  c    conventionally used by Microsoft Fortran Powerstation.
  c 
  c -- I start statements in column 7, not in column 8 (as in many of
  c    the TREGRO files).  Less whitespace (up to a point) means less
  c    wasted space and therefore makes programs more readable, because
  c    you can take in more information at a glance.
  c 
  c -- Indentation: 2 characters for progressive indentation of IF's and
  c    loops.
  c 
  c -- Case:  Upper case is traditional for FORTRAN, but unnecessary
  c    these days.  I am undecided which is more readable.
  c    Use upper or lower case 'c' for comments?
  c 
  c -- At present all the source code is in one file, Nut.f, with one
  c    include file, Nut.fi.  This is better for development since it
  c    makes global searches a little easier.  Also it's silly to repeat
  c    the same long comments at the top of each file.  There is a
  c    glossary of variables at the top, but it is not as detailed as
  c    Carol's.
  
  
  c *------------------------------------------------------------
  c * Static information: Names of state variables.
  c *------------------------------------------------------------
  
        BLOCK DATA MUMBLE
        INCLUDE 'Nut.fi'
  
  c This list of variable name strings must be kept in sync with the
  c corresponding PARAMETER definitions in Nut.fi.
  
        DATA VARNAMES /
       >  'alpha',
       >  'b',
       >  'C0',
       >  'Cav',
       >  'Cmin',
       >  'De',
       >  'Dl',
       >  'dt',
       >  'f',
       >  'f_slope',
       >  'growth',
       >  'iter',
       >  'Kd',
       >  'Km',
       >  'L',
       >  'oalpha',
       >  'r0',
       >  'rho',
       >  'rx',
       >  'source',
       >  'steps',
       >  't',
       >  'theta',
       >  'U',
       >  'Unew',
       >  'v0',
       >  'Vmax',
       >  'w0', 
       >  'Lv' /
        DATA UNITNAMES /
       >  'cm/s',         ! alpha
       >  '1',            ! b
       >  'mol/cm3',      ! C0
       >  'mol/cm3',      ! Cav
       >  'mol/cm3',      ! Cmin
       >  'cm2/s',        ! De
       >  'cm3/g',        ! Dl
       >  's',            ! dt
       >  '1',            ! f
       >  '1/s',          ! f_slope
       >  '1/s',          ! growth
       >  '1',            ! iter
       >  'cm3/g',        ! Kd
       >  'mol/cm3',      ! Km
       >  'cm',           ! L
       >  'cm/s',         ! oalpha
       >  'cm',           ! r0
       >  'g/cm3',        ! rho
       >  'cm',           ! rx
       >  'mol/cm3/s',    ! source
       >  '1',            ! steps
       >  's',            ! t
       >  '1',            ! theta
       >  'mol',          ! U
       >  'mol',          ! Unew
       >  'cm/s',         ! v0
       >  'mol/cm2/s',    ! Vmax
       >  'cm3/s',        ! w0
       >  '1/cm2'/        ! Lv  (cm/cm3)
        END
  
  
  c *------------------------------------------------------------
  c * Nut: The main program.
  c *------------------------------------------------------------
  
        PROGRAM NUT
        IMPLICIT NONE
        INCLUDE 'Nut.fi'
  
  c Functions used:
        LOGICAL READIN
        INTEGER READFIELDS, RTRIM, VARNUM, IMEMQ
  
  
  c Internal variables
        LOGICAL EXISTP
        CHARACTER*60 CFNAME, RFNAME
        INTEGER I, ITER, STEPS, J, NI
        CHARACTER*80 COMMENT
        DOUBLE PRECISION VAL
        DOUBLE PRECISION f_init, theta_init
  
  c For input switches file processing: NF = number of fields
        INTEGER NF, MAXF
        PARAMETER (MAXF=NVARS+1)
        CHARACTER*15 FIELDS(MAXF)
  
        TAB = CHAR(9)
  
        WRITE (STDOUT, 901)
   901    FORMAT (' Soil Nutrient Uptake Model'/
       >          ' Original program by Carol Zollweg'/
       >          ' Modifications by Jonathan Rees, August 1997'/)
  
  c     *------------------------------------------------
  c     * Read name of the control file (which holds both switches and inputs).
  c     *------------------------------------------------
   10   WRITE (STDOUT, *) 'Please enter the name of the control file: '
        READ (STDIN, *) CFNAME
        J = RTRIM(CFNAME)
        CFNAME = CFNAME(1:j)
        OPEN (1, FILE=CFNAME, STATUS='OLD', ERR=15)
        GOTO 18
   15   WRITE (STDOUT, *) '!! Couldn''t find control file ', CFNAME
        GOTO 10
  
   18   WRITE (STDOUT, *) '...Reading control information and data from ',
       >                  CFNAME
  
  c     *------------------------------------
  c     * Open the results file 
  c     *------------------------------------
   20   WRITE (STDOUT, *) 'Please enter the name of the results file: '
        READ (STDIN, *) RFNAME
        J = RTRIM(RFNAME)
        RFNAME = RFNAME(1:J)
  
        INQUIRE (FILE=RFNAME, EXIST=EXISTP)
        IF (EXISTP) THEN
          WRITE (STDOUT, *) 'A file by that name already exists.  ',
       >                    'Replace it? (y/n)'
          READ (STDIN, *) COMMENT
          IF (COMMENT .EQ. 'y' .OR. COMMENT .EQ. 'yes' .OR.
       >      COMMENT .EQ. 'Y' .OR. COMMENT .EQ. 'YES') GOTO 23
          GOTO 20
        ENDIF
  
   23   OPEN (3, FILE=RFNAME, ERR=25)
        REWIND (3)      !relictual. should be unnecessary
        GOTO 28
   25   WRITE (STDOUT, *) '!! Couldn''t open results file ', RFNAME,
       >  ' for writing'
        GOTO 20
  
   28   WRITE (STDOUT, *) '...Writing results to ', RFNAME
  
  c     *---------------------------------------------
  c     * Process control information
  c     *---------------------------------------------
  
  c The first few lines are an identifying comment.  They are simply copied
  c to the result file.
  c     LOOP
   30     READ (1, '(A80)', END=777) COMMENT
          J = RTRIM(COMMENT)
          WRITE (3, 33) COMMENT(1:MAX(J,1))
   33       FORMAT (A)
          IF (J .GT. 0) GOTO 30
  c     END LOOP
  
  c Initialize the variable status table, etc.
        DO 45 I = 1, NVARS
          STATUS(I) = CALCULATED
          VALS(I) = 0.0
          SCALE(I) = 1.0
   45   CONTINUE
  c Oalpha is a special case.
        STATUS(Voalpha) = FIXED
        NUMINVARS = 0
  
  c Default values
        VALS(Vdt) = 24 * 60 * 60    !seconds in a day
  
  c Subsequent lines are parameter status specifiers:
  c   name units status value...
  c     LOOP
   50     NF = READFIELDS(1, FIELDS, MAXF)
          IF (NF .LT. 0) GOTO 777
          IF (NF .EQ. 0) GOTO 59
          IF (NF .LT. 3) 
       >    WRITE(STDOUT, *) '!! Invalid status line syntax -- ',
       >                     FIELDS(1)
  
  c Fields are:   <name> <units> <status> [<parameters>]
          I = VARNUM(FIELDS(1))
          IF (I .EQ. 0) THEN
            WRITE (STDOUT, *) '!! Unrecognized variable name ', 
       >                      FIELDS(1), ' in switches table'
            GOTO 50
          ENDIF
  
  c Perform unit normalization
          IF (FIELDS(2) .EQ. 'day') THEN
            FIELDS(2) = 's'
            SCALE(I) = 1.0/(24.0 * 60.0 * 60.0)
          ELSE IF (FIELDS(2) .EQ. '1/day') THEN
            FIELDS(2) = '1/s'
            SCALE(I) = (24.0 * 60.0 * 60.0)
          ELSE IF (FIELDS(2) .EQ. 'm') THEN
            FIELDS(2) = 'cm'
            SCALE(I) = 0.01
          ENDIF
  
  c Check unit against permissible units
          IF (FIELDS(2) .NE. UNITNAMES(I)) THEN
            WRITE (STDOUT, *) '!! Unrecognized unit ', FIELDS(2),
       >                      ' for variable ', FIELDS(1)
            WRITE (STDOUT, *) '!! (Expected e.g. ', UNITNAMES(I), ')'
          ENDIF
  
  c Read the initial value and scale it accordingly
          IF (NF .LT. 4 .OR. RTRIM(FIELDS(4)) .EQ. 0) THEN
            VALS(I) = 0.0
          ELSE
            READ(FIELDS(4), *) VAL
            VALS(I) = VAL / SCALE(I)
          ENDIF
  
  c The status field is "input", "fixed", "calculated", etc.
          IF (FIELDS(3) .EQ. 'input') THEN
            STATUS(I) = INPUT
            NUMINVARS = NUMINVARS + 1
            IF (NF .GT. 3)
       >      WRITE (STDOUT, *) '!! Extraneous parameter for ',
       >        FIELDS(1)
          ELSE IF (FIELDS(3) .EQ. 'calculated') THEN
            STATUS(I) = CALCULATED
          ELSE 
            IF (NF .LT. 4) WRITE (STDOUT, *) '!! Missing parameter for ',
       >      FIELDS(1)
            IF (FIELDS(3) .EQ. 'fixed') THEN
              STATUS(I) = FIXED
            ELSE IF (FIELDS(3) .EQ. 'linear') THEN
              STATUS(I) = LINEAR
            ELSE IF (FIELDS(3) .EQ. 'exponential') THEN
              STATUS(I) = EXPONENTIAL
            ELSE
              WRITE (STDOUT, *) '!! Invalid status: ', FIELDS(3)
              WRITE (STDOUT, *) '!! (Wanted input, fixed, calculated, ',
       >                        'linear, or exponential)'
            ENDIF
          ENDIF
          GOTO 50
  c     END LOOP
  
  c Blank line encountered.
   59   CONTINUE
  
  c     *--------------------------------
  c     * Process the column names line.  This determines what will be
  c     * written to the results file, and also specifies the order of
  c     * the fields in the input records.
  c     *--------------------------------
  
        NCOLS = READFIELDS(1, FIELDS, MAXF)
        NI = 0
        DO 60 J = 1, NCOLS
          I = VARNUM(FIELDS(J))
          IF (I .EQ. 0) THEN
            WRITE (STDOUT, *) '!! Unrecognized variable name in inputs ',
       >                      'file: [', FIELDS(J), ']'
          ELSE IF (STATUS(I) .EQ. INPUT) THEN
            NI = NI + 1
            INVARS(NI) = I
          ENDIF
          COLVARS(J) = I
   60   CONTINUE
  
        IF (NI .NE. NUMINVARS) THEN
          WRITE (STDOUT, *) '!! Wrong number of inputs: switches list',
       >        ' has', NUMINVARS, ', but column headings list has ', NI
        ENDIF
  
  c The results file should have at least one column, yes?
        IF (NCOLS .EQ. 0) THEN
          WRITE (STDOUT, *) '!! Warning: no output columns specified!'
          GOTO 98
        ENDIF
  
  c The results file should have an uptake column, yes?
        IF (IMEMQ(VU, COLVARS, NCOLS) .EQ. 0) THEN
          WRITE (STDOUT, *) '!! Warning: no column specified for ',
       >    'uptake (U)'
        ENDIF
  
        if(debugp) then
          print *, '| columns: ', 
       >    (varnames(colvars(j)), j = 1, ncols)
          if (numinvars .gt. 0) then
            print *, '| inputs:  ', 
       >      (varnames(invars(j)), j = 1, numinvars)
          else
            print *, '| no inputs'
          endif
        endif
  
  c     *---------------------------------
  c     * Initialize, etc.
  c     *---------------------------------
  
  c Begin writing the results file (comment has already been written, above)
  
        DO 70 I = 1, NCOLS
          FIELDS(I) = VARNAMES(COLVARS(I))
   70   CONTINUE
        CALL WRITEFIELDS(3, FIELDS, NCOLS)
  
  c Compute format to be used in output routine
        J = 1
        FORM(1:J) = '('
        DO 71 I = J, (NCOLS-1)*10 + J, 10
          FORM(I+1:I+10) = 'g12.6,1h-,'
          FORM(I+9:I+9) = TAB
     71 CONTINUE
        FORMLEN = NCOLS*10 + J + 1
        FORM(FORMLEN:FORMLEN) = ')'
  
        if(debugp)
       >  print *, '| Format = [', FORM(1:FORMLEN), ']'
  
  c  Sanity check
        IF (NUMINVARS .LE. 0 .AND. STATUS(Vsteps) .NE. FIXED)
       >  WRITE (STDOUT, *) '!! No termination condition ',
       >                    '(no inputs, and steps isn''t fixed)'
  
  c     *-----------------------------------------------------------------
  c     * Begin iterating
  c     *-----------------------------------------------------------------
  
  c  The following are conversions to integer
        ITER = 1
        STEPS = VALS(Vsteps)
        if(debugp) print *, '| init done'
  
  c     BEGIN LOOP
   90   CONTINUE
  
  c  Read in the ith parameter set (if there is one)
          IF (NUMINVARS .GT. 0) THEN
            IF (.NOT. READIN(1)) GOTO 98
          ENDIF
  
          VALS(Viter) = ITER
          if(debugp) print *, 'Iteration #', ITER
  
  c  Perform one time step (time t to time t+dt) of uptake model
          CALL RUNMODEL
  
  c  Get initial values of f and theta to establish linear relationship
          IF (ITER .EQ. 1) THEN
            f_init = VALS(Vf)
            theta_init = VALS(Vtheta)
          ENDIF
  
  c  Write results
          WRITE (3, FORM) (VALS(COLVARS(I)) * SCALE(COLVARS(I)),
       >                   I = 1, NCOLS)
  
  c  Stop if it's time to stop
          IF (STATUS(Vsteps) .EQ. FIXED .AND. ITER .GE. STEPS)
       >    GOTO 98
  
  c  Prepare for next iteration
          CALL UPDATE(f_init, theta_init)
          VALS(Vt) = VALS(Vt) + VALS(Vdt)
          ITER = ITER + 1
  
          GOTO 90
  c     END LOOP
  
  c     *-----------------------------------------------------------------
  c     * End iterating
  c     *-----------------------------------------------------------------
  
  c  Close control and results files
   98   CLOSE(1)
        CLOSE(3)
        WRITE (STDOUT, *) 'Simulation completed.'
        WRITE (STDOUT, *) 'Ran', ITER, ' iterations.'
  
        GOTO 778
  
   777  WRITE (STDOUT,*) '!! Not enough stuff in control file'
  
   778  CONTINUE
        END
  
  
  c *---------------------------------------------------------------------
  c * RUNMODEL oversees execution of one iteration of the uptake model.
  c *---------------------------------------------------------------------
  
        SUBROUTINE RUNMODEL
        IMPLICIT NONE
        INCLUDE 'Nut.fi'
  
  c Set up inputs to the uptake model.
        alpha     = VALS(Valpha)
        b         = VALS(Vb)
        Cav       = VALS(VCav)
        De        = VALS(VDe) 
        Dl        = VALS(VDl)
        dt        = VALS(Vdt)
        f         = VALS(Vf) 
        growth    = VALS(Vgrowth)
        Km        = VALS(VKm)
        L         = VALS(VL)
        r0        = VALS(Vr0)
        rho       = VALS(Vrho)
        rx        = VALS(Vrx)
        source    = VALS(Vsource)
        theta     = VALS(Vtheta)
        Unew      = VALS(VUnew)
        v0        = VALS(Vv0)
        Vmax      = VALS(VVmax)
        w0        = VALS(Vw0)
  
  c Obsolescent?
        Cmin      = VALS(VCmin)    ! Cmin (mol/cm3)
  
        if(debugp) print *, '| L:', L
        if(debugp) print *, '| r0:', r0
        if(debugp) print *, '| c0:', C0
  
  c MODEL does the real work of computing values of dependent variables.
        CALL MODEL
  
  c Copy calculated values back into the VALS vector.
        VALS(Valpha)  = alpha
        VALS(Vb)      = b
        VALS(VCav)    = Cav
        VALS(VC0)     = C0
        VALS(VDe)     = De
        VALS(VL)      = L
        VALS(Vrx)     = rx
        VALS(Vv0)     = v0
        VALS(VU)      = U
        VALS(VUnew)   = Unew
  
        RETURN
        END
  
  
  c=======================================================================
  c Derived from OurPrism2 file UpTk.f
  
  c *------------------------------------------------------------
  c * MODEL
  c * Do all computations for one step of the uptake model.
  c *
  c * Calculate the root length for both new and established roots, the
  c * root density, and the average radial distance to the next root's
  c * zone of influence.  The type of root used depends on the soil
  c * horizon.
  c * The root length (cm) = m/gC * gC for the root * 100cm/m     
  c * The root density (m/m**3) = m/gC * gC for the root / soil volume
  c * The average radial distance to the next root's zone of      
  c * influence, rx, (cm) =                                       
  c *      1/sqrt(Pi * RootDens (m/m**3)/10000cm**2/m**2)
  c *---------------------------------------------------------------------
  
        SUBROUTINE MODEL
        IMPLICIT NONE
  
  c Common variables
        INCLUDE 'Nut.fi'
  c Function
        DOUBLE PRECISION CALCPc
  c Internal variables
        DOUBLE PRECISION Deb, A, rav
        DOUBLE PRECISION GAMMA
        DOUBLE PRECISION DeltaL, Pc, RTRATIO
        DOUBLE PRECISION UEST
  c      DOUBLE PRECISION Lv
  c Temporaries
        DOUBLE PRECISION AA, AB, AC, AD, AE, AF, FOO
  
  c EVENSMLR is important.  v0 goes down as low as 10^-11 (for Robinia
  c  pseudoacacia), so the threshold should probably be smaller than that.
  c  If v0 is smaller than that but nonzero, we should probably
  c  generate a warning, because user should be informed that the
  c  the diffusion-only formula will be invoked.
        DOUBLE PRECISION CM3_L, EVENSMLR
          PARAMETER (CM3_L = 1000.0, EVENSMLR = 1E-12)
  
  c     *-----------------------------------------------------------*
  c     * If there are no roots, there is no uptake!  Exit routine. *
  c     *-----------------------------------------------------------*
        IF (L .LE. 0.0d0) THEN
          WRITE (STDOUT,*) '!! Root length <= 0'
          RETURN
        ENDIF
  
  c     *-----------------------------------------------------------------
  c     * Compute b if requested.
  c     *-----------------------------------------------------------------
        IF (STATUS(Vb) .EQ. CALCULATED) THEN
          b = theta + rho*Kd
        ENDIF
  
  c     *-----------------------------------------------------------------
  c     * Compute De if requested.
  c     *-----------------------------------------------------------------
        IF (STATUS(VDe) .EQ. CALCULATED) THEN
          Deb = Dl * theta * f
          De = Deb / b
        ELSE
          Deb = De * b
        ENDIF
  
  c     *-----------------------------------------------------------------
  c     * Compute rx (zone of influence), if required.
  c     * Lv is root length per unit volume (cm of root per cm^3 of soil).
  c     *-----------------------------------------------------------------
        IF (STATUS(Vrx) .EQ. CALCULATED) THEN
          rx = 1.0d0 / SQRT(PI * Lv)
        ENDIF
  
  c     *---------------------------------
  c     * Calculate the root ratio rx/r0. 
  c     *---------------------------------
        RTRATIO = rx / r0
        if(debugp) print *, '| RTRATIO:', RTRATIO
  
  c     *------------------------------------------------------------
  c     * Calculate v0 and gamma.  v0 equals
  c     * the water uptake w0 (cm3/s) divided by the root surface      
  c     * area (cm2).  Don't allow the inward velocity to equal zero.
  c     * Root surface area calculated with established roots       
  c     * only to make it consistent with root water uptake w0.
  c     * Root water uptake is calculated with established roots.    
  c     *------------------------------------------------------------
        IF (STATUS(Vv0) .EQ. CALCULATED) THEN
          v0 = MAX(EVENSMLR, w0 / (2 * PI * r0 * L))
          if(debugp) print *, '| v0 (calculated):', v0
        ENDIF
        gamma = (r0 * v0)/Deb
        if(debugp) print *, '| gamma:', gamma
  
  
        IF (STATUS(Valpha) .EQ. CALCULATED) THEN
  c     *-----------------------------------------------------------------
  c     * Ruth's method for calculating root absorbing power (alpha).
  c     *-----------------------------------------------------------------
  c       -- CALCALPHA computes alpha and C0, but not Pc.
          CALL CALCALPHA(gamma, RTRATIO)
          if(debugp) print *, '| C0 (direct):', C0
          if(debugp) print *, '| alpha (direct):', alpha
          Pc = C0 / Cav
          if(debugp) print *, '| Pc (direct):', Pc
        ENDIF
  
        IF (STATUS(Voalpha) .EQ. CALCULATED) THEN
  c     *-----------------------------------------------------------------
  c     * Calculate the root absorbing power (alpha) as a function of
  c     * Vmax and Km by iterating because alpha is a function of the
  c     * concentration at the root surface and that concentration is a
  c     * function of alpha.
  c     * 
  c     * N.b. UPALPHA calculates C0 and Pc in addition to oalpha.
  c     *-----------------------------------------------------------------
          CALL UPALPHA(EVENSMLR,      
       >               0.1d0,      ! OLDC0 (OLDCSURF)
       >               Vmax, Km, Cmin, gamma,             
       >               v0, RTRATIO, r0, Deb, b, Cav,
       >               .TRUE.,
       >               alpha,      ! This parameter appears to be ignored
       >               oalpha, C0, Pc, debugp)
          if(debugp) print *, '| C0 (iterative):', C0
          if(debugp) print *, '| alpha (iterative):', oalpha
          if(debugp) print *, '| Pc (iterative):', Pc
          alpha = oalpha
        ENDIF
  
  c     *-----------------------------------------------------------------
  c     * If alpha is calculated, then Pc and C0 are calculated as a side effect.
  c     * If alpha is fixed, then we have to compute them.
  c     *-----------------------------------------------------------------
        IF (STATUS(Valpha) .NE. CALCULATED .AND.
       >    STATUS(Voalpha) .NE. CALCULATED) THEN
  c       -- If alpha is fixed, use it to compute Pc and C0.
          Pc = CALCPc(alpha, gamma, RTRATIO, Deb, EVENSMLR)
          C0 = Pc * Cav
        ENDIF
  
  c     *-----------------------------------------------------------------
  c     * Calculate rav = the radius where the concentration equals
  c     * the average concentration. If v0 is close to zero, use a
  c     * different equation.
  c     *-----------------------------------------------------------------
        IF (v0 .LE. EVENSMLR) THEN
          AA = Deb/(alpha * r0)
          AB = (1 - Pc) / Pc
          FOO = EXP(AA * AB)
        ELSE
          AA = 2 / (2 - GAMMA)
          AC = -1 / GAMMA
          AD = RTRATIO**(2 - GAMMA) - 1
          AE = (RTRATIO * RTRATIO) - 1
          FOO = ((AA * (AD/AE))**AC)
        ENDIF
        rav = r0 * FOO
        if (debugp) print *, '| rav:', rav
  
  c     *------------------------------------------------------------
  c     * Calculate the uptake due to the established roots (U_est). 
  c     *------------------------------------------------------------
        Uest = L * 2 * PI * r0 * alpha * C0 * dt
        if (debugp) print *, '| Uest:', Uest
  
  c     *-----------------------------------------------------------
  c     * Calculate the uptake by the new roots (if any).
  c     * To do this, first calculate A = the amount of solute removed
  c     * from the depletion zone per unit length of new root.  If v0
  c     * is very close to zero, use a different equation.
  c     *-----------------------------------------------------
        IF (STATUS(VUnew) .EQ. CALCULATED) THEN
          IF (v0 .LE. EVENSMLR) THEN
            AB = r0 * r0
            AC = rav * rav
            AD = AC - AB
            AE = Pc * ALPHA * r0 / Deb
            AF = LOG(rav / r0)
            A = PI * b * Cav *
       >         (((1 - Pc) * AD) - (AE * ((AC * (AF - 0.5)) + AB*0.5)))
          ELSE
            AA = 1 - Pc*alpha/v0
            AB = rav*rav - r0*r0
            AC = 2*Pc/v0
            AD = (v0 - ALPHA)/(2 - GAMMA)
            AE = ((rav/r0) ** (2 - GAMMA)) - 1
            A = PI * b * Cav *
       >         (AA*AB - (AC * AD * r0*r0 * AE))
          ENDIF
          if (debugp) print *, '| A:', A
  
  c     *-----------------------------------------------------------
  c     * Length of roots that are brand new.
  c     *-----------------------------------------------------------
          DeltaL = L * (growth * dt)
  
  c     *-----------------------------------------------------------
  c     * Don't allow uptake by new roots to be negative.
  c     **** Ruth questions this.  Maybe it _should_ be allowed to go negative.
  c     *-----------------------------------------------------------
          Unew = MAX(A * DeltaL, 0.0d0)
          if (debugp) print *, '| Unew:', Unew
        ENDIF
  
  c     *---------------------------------------------------------
  c     * The total amount of nutrient taken up by the root is
  c     * the amount up by the established roots plus the amount
  c     * taken up by the new roots.
  c     *---------------------------------------------------------
        U = Uest + Unew
        if (debugp) print *, '| U:', U
  
        RETURN
        END
  
  c=======================================================================
  c Alpha computation copied mindlessly from Ruth's handwritten notes.
  c
  c Inputs:  gamma, RTRATIO
  c          v0, Cav, Km, Vmax
  c Outputs: alpha, C0
  
        SUBROUTINE CALCALPHA(gamma, RTRATIO)
        IMPLICIT NONE
        DOUBLE PRECISION gamma, RTRATIO
  
  c Common variables.
        INCLUDE 'Nut.fi'
  
  c Internal variables.
        DOUBLE PRECISION TEMP1, TEMP2
        DOUBLE PRECISION AA, BB
        DOUBLE PRECISION ROOT1, ROOT2
  
  c Compute intermediate quantities.
          TEMP1 = RTRATIO**(2 - GAMMA) - 1
          TEMP2 = (RTRATIO * RTRATIO) - 1
        AA = (2 / (2 - GAMMA)) * (TEMP1 / TEMP2)
        if(debugp) print *, '| AA:', AA
  
        BB = v0 * Cav
        if(debugp) print *, '| BB:', BB
  
        CALL QUADRATIC(AA * v0,
       >               AA*Km*v0 + (1-AA)*Vmax - BB,
       >               BB * Km,
       >               root1, root2)
        if(debugp) print *, '| roots are', root1, root2
  
        C0 = MAX(root1, root2)
        IF (C0 .LT. 0) THEN
          WRITE (STDOUT, *) '!! negative roots in C0 computation'
        ENDIF
  
        alpha = Vmax / (Km + C0)
  
        RETURN
        END
  
  c=======================================================================
  c Update state variables for next iteration.
  
        SUBROUTINE UPDATE(f_init, theta_init)
        IMPLICIT NONE
        DOUBLE PRECISION f_init, theta_init
  
  c Common and internal variables.
        INCLUDE 'Nut.fi'
  
  c Formula for Cav copied from UpComp.f
        IF (STATUS(VCav) .EQ. CALCULATED) THEN
          Cav = Cav - U / (L * PI * (rx * rx) * b)
        ENDIF
  
        IF (STATUS(VCav) .EQ. LINEAR)
       >  Cav = Cav + (source * dt)
  
  c Root growth can be linear or exponential
        IF (STATUS(VL) .EQ. LINEAR) THEN
          L = L + (growth * dt)
        ELSE IF (STATUS(VL) .EQ. EXPONENTIAL) THEN
          L = L * (1.0d0 + (growth * dt))
        ENDIF
  
  c Update fugacity, which is linear in theta
        IF (STATUS(Vf) .EQ. CALCULATED) THEN
          f = f_init + f_slope*(theta - theta_init)
        ENDIF
  
        VALS(VCav) = Cav
        VALS(VL) = L
        VALS(Vf) = f
  
        END
  
  
  c=======================================================================
  c     *----------------------------------------------------------*
  c     * Calculate the ratio Pc = C0/Cav.
  c     * If v0 <= EvenSmlr a different calculation for Pc is
  c     * used based on the diffusion-only case of Baldwin et al.
  c     * 1973 (eqn III).  This equation is independent of  v0
  c     * and gamma.
  c     *----------------------------------------------------------*
        DOUBLE PRECISION FUNCTION 
       >  CALCPc(myalpha, gamma, RTRATIO, Deb, EVENSMLR)
        IMPLICIT NONE
        DOUBLE PRECISION myalpha, gamma, RTRATIO, Deb, EVENSMLR
        INCLUDE 'Nut.fi'
  
        DOUBLE PRECISION Pc, AA, AB, AC, AD, AE
  
        IF (v0 .LE. EVENSMLR) THEN
  c       rx = RTRATIO * r0      -- this has already been computed
          AA = myalpha * r0 / (2 * Deb)
          AB = rx * rx
          AC = AB * myalpha * r0 / Deb
          AD = AB - r0 * r0 
          AE = LOG(rx/r0)
          Pc = 1/(1 - AA + AC/AD * AE)
          if(debugp) print *, '| Pc (small v0): ', Pc
        ELSE
          AA = v0 - myalpha
          AB = 2 / (2 - gamma)
          AC = RTRATIO ** (2 - gamma)
          AD = RTRATIO * RTRATIO
          AE = (AC - 1.0) / (AD - 1.0)
          Pc = v0 / (myalpha + AA * AB * AE)
          if(debugp) print *, '| Pc (normal): ', Pc
        ENDIF
  
        CALCPc = Pc
        RETURN
        END
  
  c=======================================================================
  c *-----------------------------------
  c * Quadratic equation solver.
  c * A numerical analyst would probably tell you that this is absolutely
  c * the worst way to solve a quadratic equation, due to its
  c * vulnerability to overflow, underflow, and roundoff errors.  A better
  c * way might be found in a reference such as "Numerical Recipes," which I
  c * haven't yet consulted.
  c *-----------------------------------
  
        SUBROUTINE QUADRATIC(A, B, C, R1, R2)
        IMPLICIT NONE
        DOUBLE PRECISION A, B, C, R1, R2
        DOUBLE PRECISION DISCRIM, S
  
        DISCRIM = ((B * B) - (4 * A * C))
        IF (DISCRIM .LT. 0) THEN
          PRINT *, '!! Quadratic equation has complex roots'
          R1 = 0
          R2 = 0
          RETURN
        ENDIF
        S = SQRT(DISCRIM)
        R1 = (-B + S)/(2*A)
        R2 = (-B - S)/(2*A)
        RETURN
        END
  
  
  c=======================================================================
  c I/O utilities.
  
  c *-----------------------------------
  c * Read a line of input parameters.
  c * Returns .true. normally, .false. at end of file.
  c *-----------------------------------
  
        LOGICAL FUNCTION READIN(UNIT)
        IMPLICIT NONE
        INTEGER UNIT
  
        INCLUDE 'Nut.fi'
        INTEGER I
        DOUBLE PRECISION INPUTS(NVARS)
  
        READ (UNIT, *, END=777) (INPUTS(I), I=1,NUMINVARS)
        DO 20 I = 1, NUMINVARS
          VALS(INVARS(I)) = INPUTS(I)
    20  CONTINUE
        READIN = .TRUE.
        RETURN
  c 
   777  WRITE (STDOUT, *) 'End of input records'
        READIN = .FALSE.
        END
  
  
  c *-------------------------------------------
  c * Separate a line into tab-separated fields.
  c * MAXF = maximum number of fields (i.e. size of output array).
  c * Returns the number of fields read, or -1 on EOF.
  c *-------------------------------------------
        INTEGER FUNCTION READFIELDS(UNIT, FIELDS, MAXF)
        IMPLICIT NONE
        INTEGER UNIT, MAXF
        CHARACTER*(*) FIELDS(MAXF)
  
        CHARACTER*320 BUFFER
        CHARACTER*1 TAB
        INTEGER POS, J, NF, FIN, RTRIM
  c         
        TAB = CHAR(9)
  c         
        POS = 1
        NF = 0
        READ (UNIT, '(A150)', END=9) BUFFER 
        FIN = RTRIM(BUFFER)
   1    IF (POS .LE. FIN .AND. NF .LT. MAXF) THEN
          J = INDEX(BUFFER(POS:FIN), TAB)
          IF (J .GT. 0) THEN
             NF = NF + 1
             FIELDS(NF) = BUFFER(POS:POS+J-2)
             POS = POS + J
             GOTO 1
          ELSE
             NF = NF + 1
             FIELDS(NF) = BUFFER(POS:FIN)
          ENDIF
        ENDIF
  
        READFIELDS = NF
        RETURN
  
  c end of file
   9    READFIELDS = -1
        RETURN
        END
  
  
  c *---------------------------------------------------------------------
  c * Write NF fields, separated by tabs. 
  c * [Trailing whitespace isn't written.  -- yes it is]
  c *---------------------------------------------------------------------
  
        SUBROUTINE WRITEFIELDS(UNIT, FIELDS, NF)
        IMPLICIT NONE
        INTEGER UNIT, NF
        CHARACTER*15 FIELDS(NF)
  
        CHARACTER*320 BUFFER
        INTEGER POS, F, I
  c      INTEGER RTRIM
        CHARACTER*1 TAB
  c   
        TAB = CHAR(9)
        POS = 1
        DO 10 F = 1, NF
  c        I = RTRIM(FIELDS(F))
          I = 15     !TEMPORARY CHEAT, TRIAL PERIOD ONLY
          BUFFER(POS:POS+I-1) = FIELDS(F)(1:I)
             POS = POS + I
             IF (F < NF) THEN
                BUFFER(POS:POS) = TAB
                POS = POS + 1
             ENDIF
   10     CONTINUE
        WRITE (UNIT, 11) BUFFER(1:POS-1)
   11     FORMAT (A)
        RETURN
        END
  
  
  c *---------------------------------------------------------------------
  c * Find position of rightmost non-whitespace character in string
  c *---------------------------------------------------------------------
        FUNCTION RTRIM(S)
        IMPLICIT NONE
        CHARACTER*(*) S
        INTEGER POS, RTRIM
        POS = LEN(S)
  c     LOOP
   10     IF (S(POS:POS) .NE. ' ') THEN
            RTRIM = POS
            RETURN
          ENDIF
          POS = POS - 1
          IF (POS .LE. 0) THEN
            RTRIM = 0
            RETURN
          ENDIF
          GOTO 10
  c     END LOOP
  
        END
  
  
  c *---------------------------------------------------------------------
  c * Determine whether a given integer occurs in a vector.
  c * (MEMQ is a similar function in the Lisp programming language.)
  c *---------------------------------------------------------------------
        INTEGER FUNCTION IMEMQ(N, V, LEN)
        INTEGER N, LEN, V(LEN)
        INTEGER I
        DO 1 I = 1, LEN
          IF (V(I) .EQ. N) THEN
            IMEMQ = I
            RETURN
          ENDIF
      1 CONTINUE
        IMEMQ = 0
        RETURN
        END
  
  
  c *------------------------------------------------------------
  c * Find the index of a variable given its name.
  c *------------------------------------------------------------
        INTEGER FUNCTION VARNUM(NAME)
        IMPLICIT NONE
        INCLUDE 'Nut.fi'
        CHARACTER*(*) NAME
        INTEGER I
  
        DO 1 I = 1, NVARS
          IF (NAME .EQ. VARNAMES(I)) THEN
            VARNUM = I
            RETURN
          ENDIF
      1 CONTINUE
        VARNUM = 0
        RETURN
        END
  
  
  c=======================================================================
  c From Carol's file UpAlpha.f
  
  c Old iterative method for computing alpha.
  
  c INPUTS:  EVENSMLR, INC0, Vmax, KM, CMIN, GAMMA, v0, 
  c          RTRATIO, r0, Deb, B, Cav,
  c          CALC_alpha,                  -- flag
  c          INALPHA,      -- observed only if CALC_alpha is .FALSE.
  c
  c OUTPUTS:
  c          ALPHA, C0, Pc
  
  
  c Pc = The ratio between the concentration of the 
  c substance at the root surface and the      
  c average solution concentration.
  
        SUBROUTINE UPALPHA(EVENSMLR,   INC0,    Vmax, 
       >                   Km,         Cmin,     GAMMA,
       >                   v0,         RTRATIO,  r0,    
       >                   Deb,        b,        Cav,   
       >                   CALC_ALPHA, INALPHA,  ALPHA,        
       >                   C0,         Pc,       DEBUGP)
        IMPLICIT NONE
  c Parameters
        DOUBLE PRECISION EVENSMLR, INC0, Vmax
        DOUBLE PRECISION Km, Cmin, GAMMA
        DOUBLE PRECISION RTRATIO, r0, Deb, b, Cav, v0
        LOGICAL   CALC_ALPHA
        DOUBLE PRECISION INALPHA
        DOUBLE PRECISION ALPHA, C0, Pc
        LOGICAL DEBUGP
  
  c Function
        DOUBLE PRECISION CALCPc
  
  c Internal variables
        DOUBLE PRECISION OLDC0
        INTEGER ITER
        LOGICAL DONE
  
  c Constants
        DOUBLE PRECISION CRITERION
          PARAMETER (CRITERION = 0.001d0)
        INTEGER MAXITER
          PARAMETER (MAXITER = 100)
  
  c Initialize
        ITER = 0 
        OLDC0 = INC0
  
  c     *--------------------------------------------------------------* 
  c     * If alpha is to be calculated, then check if concentration is * 
  c     * below minimum indicating that the root absorbing power is    *
  c     * zero.  If it is, indicate that the routine converged in zero *
  c     * iterations.  Otherwise, indicate that the routine has not    *
  c     * converged.                                                   *
  c     *--------------------------------------------------------------* 
        IF (CALC_alpha) THEN 
          IF (Cav .LE. Cmin) THEN
            ALPHA = 0.0
            DONE = .TRUE.
          ELSE
            DONE = .FALSE.
          ENDIF
  
  c     *---------------------------------------------------------*
  c     * If alpha is not to be calculated make sure that routine *
  c     * does not iterate.                                       *
  c     *---------------------------------------------------------*
        ELSE
          DONE = .TRUE.
        ENDIF
  
  c     LOOP
   10     CONTINUE
  c     *-----------------------------------------------*
  c     * Increment the number of iterations completed. *
  c     *-----------------------------------------------*
          ITER = ITER + 1
  c     *-------------------------------------------------------*
  c     * If alpha is to be calculated then calculate the root  *
  c     * absorbing power based on the old value of             *
  c     * concentration at the root surface. Otherwise, set it  *
  c     * to the input value.   Vmax / (Km + c0)                *
  c     *-------------------------------------------------------*
          IF (CALC_ALPHA) THEN
             ALPHA = Vmax / (Km + OLDC0)
          ELSE
             ALPHA = INALPHA
          ENDIF
          if(debugp) print *, '| alpha in loop: ', alpha
  
          Pc = CALCPc(alpha, gamma, RTRATIO, Deb, EVENSMLR)
  
  c     *----------------------------------------------------------*
  c     * Calculate the new concentration of the substance at the  *
  c     * root surface: C_0 = P_c * C_av
  c     *----------------------------------------------------------*
          C0 = Pc * Cav
  c     *----------------------------------------------------*
  c     * Is this new concentration close enough to the old  *
  c     * concentration to consider the solution converged?  *
  c     * (Only check if alpha is being calculated)          *
  c     *----------------------------------------------------*
          IF (CALC_ALPHA) THEN
            if(debugp) print *, '| criterion: ', criterion
            IF (ABS((C0 - OLDC0) / (OLDC0 + Km)) .LE.
       >          CRITERION) THEN
              DONE = .TRUE.
            ELSE
              OLDC0 = C0
            ENDIF
          ENDIF
  
  c     *---------------------------------------------------------*
  c     * If the maximum number of iterations has been run, allow *
  c     * routine to end.                                         *
  c     *---------------------------------------------------------*
          IF (ITER .GE. MAXITER) DONE = .TRUE.
  
          IF (.NOT.DONE) GOTO 10
  c     END LOOP
  
        RETURN
        END
  