c  Program layercode.f

c  This program is a 1-dimensional biology model taken from the 3-d physical
c  biological model of McCreary, et. al 

c  --------------------------------------------------------------------------
c  Input files

c  Diurnal forcing

c  Unit 10	layercode_data.asc	Contains the required physical fields
c  Unit 11	bioadv_layer.asc	Contains the horizontal advection 
c					terms for the biological fields 
c					(optional)

c  Output file
c  Unit 12	layercode_out.dat	Contains the calculated biological
c					fields at each timestep

c  No diurnal forcing

c  Unit 10	layercode_nod_data.asc	Contains the required physical fields
c  Unit 11	bioadv_nod_layer.asc	Contains the horizontal advection 
c					terms for the biological fields 
c					(optional)

c  Output file
c  Unit 12	layercode_nod_out.dat	Contains the calculated biological
c					fields at each timestep
c  --------------------------------------------------------------------------

c  Parameters for the biological model are the following:

c  az                Growth efficiency of zooplankton (dimensionless) 
c  an                Nitrogen generation coefficient due to grazing (dimensionless) 
c  deathp            Death rate of phytoplankton (s^-1)
c  deathz            Death rate of zooplankton (s^-1)
c  grazp             Zooplankton grazing rate of phytoplankton (s^-1 (umolesN Kg^-1)^-1)
c  grazz             Zooplankton grazing rate of zooplankton (s^-1 (umolesN Kg^-1)^-1)
c  eps               Nitrogen generation coefficient due to zoo death (dimensionless) 
c  pgrowth           Phytoplankton growth rate (s^-1 (umolesN Kg^-1)^-1)
c  e                 Remineralization rate of detritus (s^-1)
c  rkn               N-specific phytoplankton light attenuation coefficient 
c                       (m^-1 (umolesN Kg^-1)^-1)
c  rkw               Diffuse attenuation coefficient for clear seawater (m^-1)
c  q0                Normalization light intensity at the ocean surface (watts m^-2) 
c  wsink             Sinking rate of detritus (m/s)
c  gamma             Cross-thermocline mixing (s-1)
c  rntss             Nitrogen concentration below the mixed layer

      parameter (imax=162)
      parameter (jmax=110)

      logical rntadv,phyadv,zooadv,detadv,diurnal,bioadvect

      dimension alightn(4)
      dimension hkm(4)
      dimension rnt(4),phy(4),zoo(4),det(4)
      dimension rntk(4),phyk(4),zook(4),detk(4)
      dimension rntkm(4),phykm(4),zookm(4),detkm(4)
      dimension rntadva(4),phyadva(4),zooadva(4),detadva(4)

      dimension div(4)

      common/ecovar/az,an,deathp,deathz,phip,phiz,grazz,f0,eps,pgrowth,
     .              e,rkn,rkw,q0,wsink,gamma,gammap,gammad,gamma3,rnt0,
     .              rnt4,rntst,altnrm,parfac,convh,convv,convq,gamh,
     .              gam3h,gamph,gamdh,gamprh,gampprh,gamdprh,wsinkp

c  Specify biological advection logicals
      data rntadv,phyadv,zooadv,detadv/.true.,.true.,.true.,.true./

c  Specify diurnal or no diurnal forcing data
      data diurnal/.true./

c  Calculate bioadvect logical 
      if(rntadv.or.phyadv.or.zooadv.or.detadv) then
        bioadvect=.true.
      else
        bioadvect=.false. 
      endif

c  Open input files
      if(diurnal) then
        open(10,file='layercode_data.asc',status='old')
        if(bioadvect) then
          open(11,file='bioadv_layer.asc',status='old')
        endif
c       Open output file
        open(12,file='layercode_out.dat',status='new',
     .       form='unformatted')
      else
        open(10,file='layercode_nod_data.asc',status='old')
        if(bioadvect) then
          open(11,file='bioadv_nod_layer.asc',status='old')
        endif
c       Open output file
        open(12,file='layercode_nod_out.dat',status='new',
     .       form='unformatted')
      endif

      pi=4.0*abs(atan(1.))
      rearth=6370.e5
      dypdeg=pi/180.*rearth

c  Set basin and resolution parameters
      dt=86400./48.
      tdt=2.*dt
      ymin=-29.*dypdeg
      ymax=+25.*dypdeg
      xmin=0.*dypdeg
      xmax=+80.*dypdeg
      xbas=xmax-xmin
      dx=(xmax-xmin)/(imax-2)
      dy=(ymax-ymin)/(jmax-2)
      recdx=1./dx
      recdx2=recdx*recdx
      recdy=1./dy
      recdy2=recdy*recdy

c  Set up biology parameters.
      call ecoparams

      q02=q0*q0

      indm=2
      ind2=2

      ntime=0
5     continue

c  Read in 1-d fields
      read(10,*,end=900)djd
      read(10,*)w12a,w23a,w34a,sol,qsol
      read(10,*)h1b,h2b,h3b,h4b
      read(10,*)rntk,rntkm
      read(10,*)phyk,phykm
      read(10,*)zook,zookm
      read(10,*)detk,detkm
      read(10,*)div

      if(bioadvect)then
        read(11,*)h1b,h2b,h3b,h4b
        read(11,*)rntadva,phyadva,zooadva,detadva
      endif

      ntime = ntime+1

C  Calculate biology fields
      w12b = w12a*convv
      w12e = .5*(w12b+abs(w12b))
      w12d = .5*(w12b-abs(w12b))
      w23b = w23a*convv
      w23e = .5*(w23b+abs(w23b))
      w23d = .5*(w23b-abs(w23b))
      w34b = w34a*convv
      w34e = .5*(w34b+abs(w34b))
      w34d = .5*(w34b-abs(w34b))

c  Limited light growth function
      qtop = sol*parfac*convq

      rk = rkw
      rlinm = qtop
      rz0 = rlinm+sqrt(rlinm*rlinm+q02)
      expzhm = exp(-rk*h1b)
      rzhm = rlinm*expzhm+sqrt(rlinm*rlinm*expzhm*expzhm+q02)
      alightn(1) = alog(rz0/rzhm)/(h1b*rk)

      rk = rkw
      rlinf = rlinm*expzhm
      rzhm = rlinf+sqrt(rlinf*rlinf+q02)
      expzhf = exp(-rk*h2b)
      rzhf = rlinf*expzhf+sqrt(rlinf*rlinf*expzhf*expzhf+q02)
      if(h2b.ne.0.) then
        alightn(2) = alog(rzhm/rzhf)/(h2b*rk)
      else
        alightn(2) = 0.
      endif

      rk = rkw
      rlin3 = rlinf*expzhf
      rzhf = rlin3+sqrt(rlin3*rlin3+q02)
      expzh3 = exp(-rk*h3b)
      rzh3 = rlin3*expzh3+sqrt(rlin3*rlin3*expzh3*expzh3+q02)
      alightn(3) = alog(rzhf/rzh3)/(h3b*rk)

      rk = rkw
      rlin4 = rlin3*expzh3
      rzhf = rlin4+sqrt(rlin4*rlin4+q02)
      expzh4 = exp(-rk*h4b)
      rzh4 = rlin4*expzh4+sqrt(rlin4*rlin4*expzh4*expzh4+q02)
      alightn(4) = alog(rzhf/rzh4)/(h4b*rk)
      
C  Layer 1 Nitrogen (implicit scheme)

      fi = phip*phykm(1)+phiz*zookm(1)
      fdenom = fi + f0
      sfi = fi/fdenom
      spi = phip*phykm(1)/fdenom
      szi = phiz*zookm(1)/fdenom

       a = rntkm(1) + tdt * (
     . + an * grazz * sfi * zookm(1)
     . + e * detkm(1)
     . - pgrowth * alightn(1) * phykm(1) * rntkm(1)/(rntkm(1)+rnt0) )

      if(rntadv) a = a + tdt * rntadva(1)

      b = tdt*(w12e+gamh)
      rnum = a*h1b+b*rntkm(indm)
      denom = (h1b+b)
      rnt(1) = rnum/denom

c  Layer 1 Phytoplankton (implicit scheme)
      a = phykm(1) + tdt * (
     .   pgrowth * alightn(1) * phykm(1) * rntkm(1)/(rntkm(1)+rnt0)
     . - grazz * zookm(1) * spi
     . - deathp * phykm(1) )

      if(phyadv) a = a + tdt * phyadva(1)

      b = tdt*(w12e+gamph)
      rnum = a*h1b+b*phykm(indm)
      denom = (h1b+b)
      phy(1) = rnum/denom

c  Layer 1 Zooplankton (implicit scheme)
      a = zookm(1) + tdt * (
     . + az * grazz * sfi * zookm(1)
     . - grazz * szi * zookm(1)
     . - deathz * zookm(1) )

      if(zooadv) a = a + tdt * zooadva(1)

      b = tdt*(w12e+gamh)
      rnum = a*h1b+b*zookm(indm)
      denom = (h1b+b)
      zoo(1) = rnum/denom

c  Layer 1 Detritus (implicit scheme)
      a = detkm(1) + tdt * (
     . + (1 - az - an) * grazz * sfi * zookm(1)
     . + deathp * phykm(1)
     . + deathz * zookm(1)
     . - e * detkm(1) )

      if(detadv) a = a + tdt * detadva(1)

      b = tdt*(w12e+gamdh)
      rnum = a*h1b+b*detkm(indm)
      denom = (h1b+b+wsink*tdt)
      det(1) = rnum/denom

c  Layer 2 Nitrogen (implicit scheme)
      fi = phip*phykm(2)+phiz*zookm(2)
      fdenom = fi+f0
      sfi = fi/fdenom
      spi = phip*phykm(2)/fdenom
      szi = phiz*zookm(2)/fdenom

       a = rntkm(2) + tdt*(
     . + an * grazz * sfi * zookm(2)
     . +e * detkm(2)
     . -pgrowth * alightn(2) * phykm(2) * rntkm(2)/(rntkm(2)+rnt0) )

      if(rntadv) a = a + tdt * rntadva(2)

      b = tdt*(w12d-gamh)
      c = tdt*(w23e+gamprh)
      rnum = a * h2b - b * rntkm(1) + c * rntkm(3)
      denom = (h2b-b+c)
      rnt(2) = rnum/denom

c  Layer 2 Phytoplankton (implicit scheme)
      a = phykm(2) + tdt*(
     .   pgrowth * alightn(2) * phykm(2) * rntkm(2)/(rntkm(2)+rnt0)
     . - grazz * zookm(2) * spi
     . - deathp * phykm(2) )

      if(phyadv) a = a + tdt * phyadva(2)

      b = tdt * (w12d-gamph)
      c = tdt * (w23e+gampprh)
      rnum = a * h2b - b * phykm(1) + c * phykm(3)
      denom = (h2b - b + c)
      phy(2) = rnum/denom

c  Layer 2 Zooplankton (implicit scheme)
      a = zookm(2) + tdt * (
     . + az * grazz * sfi * zookm(2)
     . - grazz * szi * zookm(2)
     . - deathz * zookm(2) )

      if(zooadv) a = a + tdt * zooadva(2)

      b = tdt * (w12d - gamh)
      c = tdt * (w23e + gamprh)
      rnum = a * h2b - b * zookm(1) + c * zookm(3)
      denom = (h2b-b+c)
      zoo(2) = rnum/denom

c  Layer 2 Detritus (implicit scheme)
      a = detkm(2)+tdt*(
     . + (1 - az - an) * grazz * sfi * zookm(2)
     . + deathp * phykm(2)
     . + deathz * zookm(2)
     . - e * detkm(2) )

      if(detadv) a = a + tdt * detadva(2)

      b = tdt * (w12d - gamdh)
      c = tdt * (w23e + gamdprh)
      rnum = a * h2b - b * detkm(1) + c * detkm(3) 
     .       + wsink * tdt * detkm(1)
      denom = (h2b-b+c+wsinkp*tdt)
      det(2) = rnum/denom

c  Layer 3 Nitrogen
      fi = phip * phykm(3) + phiz * zookm(3)
      fdenom = fi+f0
      sfi = fi/fdenom
      spi = phip * phykm(3)/fdenom
      szi = phiz * zookm(3)/fdenom

      rnt(3) = rntkm(3)+ tdt * (
     . + an * grazz * sfi * zookm(3)
     . + e * detkm(3)
     . - pgrowth * alightn(3) * phykm(3) * rntkm(3)/(rntkm(3)+rnt0)
     . + (w23d - gamh) * (rntkm(3) - rntkm(ind2))/h3b 
     . - (w34e + gamprh) * (rntkm(3) - rntkm(4))/h3b)

      if(rntadv) rnt(3) = rnt(3) + tdt * rntadva(2)

c  Layer 3 Phytoplankton
      phy(3) = phykm(3) + tdt * ( 
     .   pgrowth * alightn(3) * phykm(3) * rntkm(3)/(rntkm(3)+rnt0)
     . - grazz * zookm(3) * spi
     . - deathp * phykm(3)
     . + (w23d-gamph)*(phykm(3)-phykm(ind2))/h3b 
     . - (w34e+gampprh)*(phykm(3)-phykm(4))/h3b)

      if(phyadv) phy(3) = phy(3) + tdt * phyadva(3)

c  Layer 3 Zooplankton
      zoo(3) = zookm(3)+tdt * (
     . + az * grazz * sfi * zookm(3)
     . - grazz * szi * zookm(3)
     . - deathz * zookm(3)
     . + (w23d - gamh) * (zookm(3) - zookm(ind2))/h3b
     . - (w34e + gamprh) * (zookm(3) - zookm(4))/h3b) 

      if(zooadv) zoo(3) = zoo(3) + tdt * zooadva(3)

c  Layer 3 Detritus
      det(3) = detkm(3)+tdt*(
     . + (1 - az - an) * grazz * sfi * zookm(3)
     . + deathz * zookm(3)
     . + deathp * phykm(3)
     . - e * detkm(3)
     . + wsinkp/h3b * (detkm(ind2) - detkm(3))
     . + (w23d - gamdh) * (detkm(3) - detkm(ind2))/h3b  
     . - (w34e + gamdprh) * (detkm(3) - detkm(4))/h3b)

      if(detadv) det(3) = det(3) + tdt * detadva(3)

c  Layer 4 Nitrogen
      fi = phip * phykm(4) + phiz * zookm(4)
      fdenom = fi + f0
      sfi = fi/fdenom
      spi = phip * phykm(4)/fdenom
      szi = phiz * zookm(4)/fdenom

      rnt(4) = rntkm(4)+ tdt*(
     . + an * grazz * sfi * zookm(4)
     . + e * detkm(4)
     . - pgrowth * alightn(4) * phykm(4) * rntkm(4)/(rntkm(4)+rnt0)
     . + (w34d - gamh) * (rntkm(4) - rntkm(3))/h4b 
     . - gam3h * (rntkm(4) - rnt4)/h4b)

      if(rntadv) rnt(4) = rnt(4) +tdt*rntadva(4)

c  Layer 4 Phytoplankton
      phy(4) = phykm(4) + tdt * (
     .   pgrowth * alightn(4) * phykm(4) * rntkm(4)/(rntkm(4)+rnt0)
     . - grazz * zookm(4) * spi
     . - deathp * phykm(4)
     . + (w34d-gamph)*(phykm(4)-phykm(3))/h4b)

      if(phyadv) phy(4) = phy(4) + tdt * phyadva(4)

c  Layer 4 Zooplankton
      zoo(4) = zookm(4)+tdt * (
     . + az * grazz * sfi * zookm(4)
     . - grazz * szi * zookm(4)
     . - deathz * zookm(4)
     . + (w34d-gamh)*(zookm(4)-zookm(3))/h4b)

      if(zooadv) zoo(4) = zoo(4) + tdt * zooadva(4)

c  Layer 4 Detritus
      det(4) = detkm(4)+tdt*(
     . + (1 - az - an) * grazz * sfi * zookm(4)
     . + deathz * zookm(4)
     . + deathp * phykm(4)
     . - e * detkm(4)
     . + wsinkp/h4b * (detkm(3)-detkm(4))
     . + (w34d-gamdh)*(detkm(4)-detkm(3))/h4b )

      if(detadv) det(4) = det(4) + tdt * detadva(4)

c  set all biological variables nonnegative
      do kk=1,4
       rnt(kk) = amax1(rnt(kk),0.)
       phy(kk) = amax1(phy(kk),0.)
       zoo(kk) = amax1(zoo(kk),0.)
       det(kk) = amax1(det(kk),0.)
      enddo 

c  Write out current time step data
      write(12)rnt,phy,zoo,det

c  Increment time level values
      do n=1,4
       rntkm(n) = rntk(n)
       phykm(n) = phyk(n)
       zookm(n) = zook(n)
       detkm(n) = detk(n)

       rntk(n) = rnt(n)
       phyk(n) = phy(n)
       zook(n) = zoo(n)
       detk(n) = det(n)
      enddo

      goto 5
900   close(10)

      end

      subroutine ecoparams
c  This subroutine sets the parameters used in the biological model

      common/ecovar/az,an,deathp,deathz,phip,phiz,grazz,f0,eps,pgrowth,
     .              e,rkn,rkw,q0,wsink,gamma,gammap,gammad,gamma3,rnt0,
     .              rnt4,rntst,altnrm,parfac,convh,convv,convq,gamh,
     .              gam3h,gamph,gamdh,gamprh,gampprh,gamdprh,wsinkp

c  Set parameters for the biological model

      data  az,  an,   deathp,    deathz, phip,   phiz,   grazz
     .     /.1,  .4,  1.1574E-6,   0.,   .8333333, .1666666, 2.3148E-5/

      data f0/1./ 
      data parfac/.4/   

      data  eps, pgrowth,   e,          rkn,          rkw
     .     /0.0, 2.8935E-5, 5.7870E-7,  0.,          3.0E-2/ 

      data  q0,   wsink,       gamma,      gammap,    gammad
     .     /40.  , 4.6296E-5 ,  1.1574E-8,  1.1574E-8, 1.1574E-8/  

      data gamma3/5.7870e-9/ 

      data rntst,rnt0,rnt4/20.0,  1.0,  20.0/

      data hfac/100./   !"Typical" depth for gamma mixing (m)

      gamh=gamma*hfac
      gam3h=gamma3*hfac
      gamph=gammap*hfac
      gamdh=gammad*hfac

      gamprh=gamh
      gampprh=gamph
      gamdprh=gamdh

      wsinkp=wsink

c  Define conversion factors for use with biological model
c  the biological model expects q, h, and v in units of
c  watts/(meters)**2, meters, and meters/sec, respectively.
      convh=1./100.
      convv=1./100.
      convq=1./2.39e-5

c  Print out values
      print *,'az= ',az
      print *,'an= ',an
      print *,'deathp= ',deathp
      print *,'deathz= ',deathz
      print *,'phip= ',phip
      print *,'phiz= ',phiz
      print *,'grazz= ',grazz
      print *,'f0= ',f0
      print *,'eps= ',eps
      print *,'pgrowth= ',pgrowth
      print *,'e= ',e
      print *,'rkn= ',rkn
      print *,'rkw= ',rkw
      print *,'q0= ',q0
      print *,'wsink= ',wsink
      print *,'gamma= ',gamma
      print *,'gammap= ',gammap
      print *,'gammad= ',gammad
      print *,'gamma3= ',gamma3
      print *,'rnt0= ',rnt0
      print *,'rnt4= ',rnt4
      print *,'rntst= ',rntst
      print *,'altnrm= ',altnrm
      print *,'parfac= ',parfac

      return        
      end