      Subroutine DINOGROW(GROADD,GROW,ZEROADD,CW,II)
C
      Include 'comdeck'
      Include 'dino.inc'
C-----------------------------------------------------------------------
C This subroutine calculates an exponential growth rate as
C a function of temperature, salinity and light
C
C Variables from dino.inc
C DIN_OBS: array of values of dissolved inorganic
C          nitrogen in the surface layer
C DASWRA: array of the day-averaged shortwave irradiance
C
C Passed Variables:
C GROADD: tracks the total number of cells added via growth
C         summed over the water column (2D)
C GROW: tracks the growth rate (day-1) for each grid cell
C ZEROADD: tracks the number of cells added when gaurding against
C          negative concentration.  Ideally, this is 0, and it should 
C          be if you use a positive definite advection scheme.  sporadic
C          cfl violations can also cause some problems here.
C CW: the cell concentration at each grid cell (cells/m3)
C II: index for multiple biological settings
C
C Local Variables
C CWOLD: Save previous value of concentrations (cells/m3)
C        (used for calculation of the number of cysts added via germination)
C C_DEPTH: Compensation depth (0 @ surface, negative down)
C
C variables recieving biological model parameters:
C G_EFF: Growth Efficiency (day-1 m2 w-1)
C MORT: Mortality (days-1)
C GMAX: maximum growth rate at optimal T,S (days-1)
C GR: The maintenance respiration rate (days-1)
C AKSN: Half saturation for nutrients (micro-molar)
C ATTENK: diffuse attenuation coefficient (m-1, (+))
C DASWR_AVG: mean day averaged short-wave irradiance
C            (used for approximation of the compensation depth)
C 
C local 1x1 variables included in loops/formulas
C RATE: final growth rate days-1 for a grid cell 
C GROWTHFAC: multiplier of old concentration to calculate new
C TTT: dummy variable for temperature(deg. C)
C SSS: dummy variable for salinity (ppt)
C DASWR: dummy variable for day-averaged shortwave radiations (watts/m2)
C T_FAC: Temperature scaling factor
C S_FAC: Salinity scaling factor
C G_FAC: Total growth adjustment due to temperature/salinity
C RAD: Stores irradiance (watts/m2) at mid-point of each grid cell
C DIN: Concentration of dissolved iorganic nitrogen (micro-molar)
C G_LIGHT: Light limited growth rate (days-1)
C G_DIN: DIN limited growth rate (days-1)
C TDAY: Floating point variable storing the day of the simulation 
C ITDAY: Integer variable based on TDAY which references the day-averaged
C        irradiance
C
C from hydrodynamic model
C DT: The total water column depth (+, meters) 
C ZZ: normalized decimal position of the mid-point of each sigma layer (-) 
C     e.g. -0.005 -0.025 -0.07 -0.1 -0.2 -0.3 ....
C DZ: with of each sigma layer as decimal fraction of water column
C ART: area of each grid cell (m2)
C FSM: free surface mask (1 if part of computational domain, 0 not) 
C THOUR: number of hours since the start of the run
C T: temperature (deg. C)
C S: salinity (ppt)
C
C-----------------------------------------------------------------------
      REAL CW(IM,JM,KB), CWOLD(IM,JM,KB)
      REAL C_DEPTH,G_EFF,MORT,GMAX,GR,AKSN,ATTENK
      REAL DASWR_AVG,RATE,GROWTHFAC,DASWR
      REAL TTT,SSS,T_FAC,S_FAC,G_FAC
      REAL RAD,DIN,G_LIGHT,G_DIN
      REAL GROADD(IM,JM),GROW(IM,JM,KB),ZEROADD(IM,JM)
      REAL TDAY
      INTEGER ITDAY

C reference day averaged irradiance to noon, local time
C THOUR = 0 at midnight GMT = 16 hours before noon EDT
C at 1600 GMT the first day = 1200 EDT, ITDAY = 1.
      TDAY = THOUR / 24.
      ITDAY = NINT(TDAY-0.67) + 1
      IF (ITDAY.EQ.0) THEN
        ITDAY = 1
      END IF
      DASWR = DASWRA(ITDAY)
    
C set parameter values 
      G_EFF = BIOMOD(II,5)
      MORT = BIOMOD(II,4) 
      GMAX = BIOMOD(II,1)
      GR = BIOMOD(II,6)
      AKSN = BIOMOD(II,3)
      ATTENK = BIOMOD(II,7)

C zero out some some arrays
      DO 31 I = 2,IM
        DO 21 J = 2,JM
            GROADD(I,J) = 0.
            ZEROADD(I,J) = 0.
            DO 11 K = 1,KB
              GROW(I,J,K) = 0. 
11	    CONTINUE
21      CONTINUE
31    CONTINUE

C assume that deep newly germinated cells can get up to the
C compensation depth without losses. tanh(x) --> x as x gets
C small - use mean day-averaged  irradiance (watts) to get
C the compensation depth (C_DEPTH)

      DASWR_AVG = 345.5
      C_DEPTH = 1.0/ATTENK*LOG(GR/(G_EFF*DASWR_AVG))

      Do 30 I = 2,IMM1 
        Do 20 J = 2,JMM1 
          Do 10 K = 1,KBM1
            If (DT(I,J).GT.0.0) then
             TTT=t(i,j,k)
             SSS=s(i,j,k)

C Temperature factor, polynomial for < 5, 1 term Taylor Expansion < 5
C 0.343 = first derivative at 5
             IF (TTT.GE.5.0) THEN 
               T_FAC=-0.000513*TTT**3+0.0160*TTT**2-0.0867*TTT+0.382
             ELSE
               T_FAC=-0.000513*5.0**3+0.0160*5.0**2-0.0867*5.0+0.382
               T_FAC=T_FAC - 0.0343*(5.0-TTT)
             END IF
C salinity adjustment
             S_FAC=0.0000882*SSS**3-0.00808*SSS**2+0.220*SSS-0.872
             G_FAC = S_FAC*T_FAC

             RAD = DASWR*exp(ATTENK*ZZ(K)*DT(I,J))
             IF (ZZ(K)*DT(I,J).GT.C_DEPTH) THEN
	       G_LIGHT = (GMAX*G_FAC+GR)*
     &                 TANH(G_EFF*RAD/(GMAX*G_FAC+GR)) - GR
             ELSE
               G_LIGHT = 0
             END IF

C setting DIN from observations

	     IF ((ZZ(K)*DT(I,J)).GT.C_DEPTH) THEN
	       DIN = DIN_OBS(I,J)
	     ELSE
	       DIN = 15.0
	     END IF

             IF (AKSN.GT.0.0) THEN
	       G_DIN = GMAX*G_FAC*DIN/(AKSN+DIN)  
             ELSE
               G_DIN = GMAX*G_FAC
             ENDIF

C growth is either DIN or light limited

	     RATE = MIN(G_LIGHT,G_DIN)

C apply loss terms if above mean critical depth
            IF (ZZ(K)*DT(I,J).GT.C_DEPTH) THEN
              RATE = RATE - MORT 
            END IF

C note that the two is included here because the version of ECOM used
C has a leap-frog scheme.  Adjust accordingly if this is not the case
C with your code.
            GROWTHFAC = 1 + 2.0*dti * (RATE/(3600.*24.))
            GROW(I,J,K) = RATE

C track concentration before growth for calculation of cells added during
C growth.
            CWOLD(I,J,K) = max(CW(I,J,K),0.0) 

C guard against negative values - the array zeroadd tracks any additions
C due to this fix.  Should be 0 unless 1) you are not using a positive
C definite advection scheme, or 2) you have many CFL violations
C track the number of cells added with zeroadd - note that 0.5 is
C included to track over one time step for output

            IF ((CW(I,J,K)*GROWTHFAC).LT.0) THEN
              ZEROADD(I,J) = ZEROADD(I,J)+ART(I,J)*DZ(K)*DT(I,J)*
     &                       FSM(I,J)*0.5*ABS(CW(I,J,K)*GROWTHFAC)
            END IF
            
            CW(I,J,K) = max(CW(I,J,K) * GROWTHFAC,0.0)

C track the number of cells added with zeroadd - note that 0.5 is
C included to track over one time step for output

            GROADD(I,J)=GROADD(I,J)+(CW(I,J,K)-CWOLD(I,J,K))*
     &        ART(I,J)*DZ(K)*DT(I,J)*FSM(I,J)*0.5
            
           end if
  
   10     Continue
   20   Continue
   30 Continue
      RETURN
      End
