subroutine model

	include 'common_blocks.h'	

        real*8 bbc, wadv(50),vda(50),ds,dconst
	real*8 ar(50), br(50), cr(50), rr(50)
c
cccccccccccccccccc    START MAIN LOOP OVER TIME    cccccccccccccccccccc
c
 	do jmain=2,isteps
	   istep = jmain 
	   do n=1,nz+1
	      wadv(n)=wvel(n,istep)
	   enddo
c
cccccccc BIOLOGICAL PROCESSES ccccccccccccccccccccccccccccccccccccccc
c
 	   call biosub(istep)
c
cccc PHYSICAL PROCESSES: VERTICAL ADVECTION, DIFFUSION, SINKING ccccc
c       
cccccccccc coefficients for tridiagonal matrix (diffusion) ccccccccccc
c
           dconst = delt/(2.0d0*(dz**2.d0))
           do n=2,nz-1
	      rkzm = (rkz(n-1,istep) + rkz(n,istep))/2.d0
	      rkzp = (rkz(n+1,istep) + rkz(n,istep))/2.d0
	      ar(n) = -dconst*rkzm
	      br(n) = 1.d0 + (dconst*(rkzm+rkzp))
	      cr(n) = -dconst*rkzp
	   enddo
	   rkzp = (rkz(1,istep) + rkz(2,istep))/2.d0
	   br(1) = 1.d0 + (dconst*rkzp)
	   cr(1) = -dconst*rkzp
	   rkzm = (rkz(nz,istep) + rkz(nz-1,istep))/2.d0
	   rkznzp1 = rkz(nz,istep)+(rkz(nz,istep)-rkz(nz-1,istep))
	   rkzp = (rkz(nz,istep) + rkznzp1)/2.d0
	   ar(nz)= -dconst*rkzm
	   br(nz)= 1.d0 + (dconst*(rkzm+rkzp))
c       
ccccccccccccccccccc STEP THROUGH ALL STATE VARIABLES cccccccccccccccc
c
	   do inb=1,nsv
c       
ccccccccccccccccccccccc   boundary conditions   ccccccccccccccccccccc
c
	      if(bbcnsv(inb).eq.-9.) then
		 bbc=bio(nz,istep,inb)-
     &                (bio(nz-1,istep,inb)-bio(nz,istep,inb))
	      else
		 bbc=bbcnsv(inb)
	      endif
c       
ccccccccccccccccccccc   optional surface aeolian flux   ccccccccccccc
c
              if(aeonsv(inb).ne.0.) 
     &           bio(1,istep,inb)=bio(1,istep,inb)+
     &            aeoflux(istep)*delt/(dz*86400.d0)           
c
ccccccccccccccccccccccc    vertical advection    cccccccccccccccccccc
c     
	      vda(1) = 0.5d0*(wadv(2)*
     &               (bio(1,istep,inb)-bio(2,istep,inb)))
     &               *delt/dz
	      bio(1,istep,inb) = bio(1,istep,inb) + vda(1)
c
	      do j=2,nz-1
		 vda(j) = 0.5d0*
     &           ( (wadv(j)*(bio(j-1,istep,inb)-bio(j,istep,inb)))     
     &           + (wadv(j+1)*(bio(j,istep,inb)-bio(j+1,istep,inb))) )
     &           * delt/dz
		 bio(j,istep,inb) = bio(j,istep,inb) + vda(j)
	      enddo
c
	      vda(nz) = 0.5d0*
     &           (wadv(nz)*(bio(nz-1,istep,inb)-bbc))
     &           * delt/dz
	      bio(nz,istep,inb) = bio(nz,istep,inb) + vda(nz)
c
ccccccccccccccccccccccc    vertical diffusion    ccccccccccccccccccc
c       
	      do n=2,nz-1
		 rkzm = (rkz(n-1,istep)+rkz(n,istep))/2.d0
		 rkzp = (rkz(n+1,istep)+rkz(n,istep))/2.d0
		 rr(n) = dconst*rkzm*bio(n-1,istep,inb)+
     &                (1.d0-(dconst*(rkzp+rkzm)))*bio(n,istep,inb)
     &                +dconst*rkzp*bio(n+1,istep,inb)
	      enddo
	      rkzp=(rkz(1,istep)+rkz(2,istep))/2.d0
	      rr(1)=(1.d0-(dconst*rkzp))*bio(1,istep,inb)
     &                +dconst*rkzp*bio(2,istep,inb)
	      rkzm=(rkz(nz,istep)+rkz(nz-1,istep))/2.d0
	      rkznzp1=rkz(nz,istep)+(rkz(nz,istep)-rkz(nz-1,istep))
	      rkzp=(rkz(nz,istep)+rkznzp1)/2.d0 

	      rr(nz)=dconst*rkzm*bio(nz-1,istep,inb)
     &             +(1.d0-(dconst*(rkzp+rkzm)))*bio(nz,istep,inb)
     &             + 2.d0*dconst*rkzp*bbc
	      call tridag(ar,br,cr,rr,bio(1,istep,inb),nz)
c
ccccccccccccccccccccc    horizontal advection    cccccccccccccccccccc
c
	      if(HAnsv(inb).ne.0) call horizadv(istep,inb)
c
cccccccccccccccccccccccc    sinking    cccccccccccccccccccccccccccccc
c
	      if(wnsv(inb).ne.0.) call sink(istep,inb)
c       
cccccccccccccc  step forward to next state variable  ccccccccccccccc
c
	   enddo
c       
cccccccccccccccccccccccc MIXING DOWN TO MLD cccccccccccccccccccccccc
c
	   ndml=int(dml(istep)/dz) !calculate no. of cells in mixed layer
	   do inb = 1,nsv
	      ds = 0.d0
	      do j=1,ndml
		 ds = ds + bio(j,istep,inb)/real(ndml)
	      enddo
	      do j=1,ndml
		 bio(j,istep,inb) = ds
	      enddo
	   enddo
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
	enddo
	end

C---------------------------------------------------------------------+
C       SUBROUTINE BIOSUB
C---------------------------------------------------------------------+
        subroutine biosub(istep)

        include 'common_blocks.h'  

        real*8 kw,kchl,hh,kb,riob(50)
        real*8 yt(nbio),dyt(nbio)
        real*8 tempbio(nbio),dydtemp(nbio)
        integer istep

	kw=0.05d0
	kchl=0.1d0
	hh=delt*0.5d0
	riob(1) = qi(istep)*0.4d0 !convert back from total radiation to PAR

c  *************** MAIN LOOP OVER DEPTH *********************

        do k = 1,nz

	   ndep=dint(dz*k-dz/2.)
	   if(ndep.le.10) then
	      Temp(k)=Tdat(istep,1)
	   elseif(ndep.gt.10.and.ndep.le.30) then
	      Temp(k)=Tdat(istep,2)
	   elseif(ndep.gt.30.and.ndep.le.50) then
	      Temp(k)=Tdat(istep,3)
	   elseif(ndep.gt.50.and.ndep.le.70) then
	      Temp(k)=Tdat(istep,4)
	   elseif(ndep.gt.70.and.ndep.le.90) then
	      Temp(k)=Tdat(istep,5)
	   elseif(ndep.gt.90.and.ndep.le.110) then
	      Temp(k)=Tdat(istep,6)
	   elseif(ndep.gt.110.and.ndep.le.130) then
	      Temp(k)=Tdat(istep,7)
	   elseif(ndep.gt.130.) then
	      Temp(k)=Tdat(istep,8)
	   endif

           kb = kw+kchl*bio(k,istep-1,ichl)
	   rlt = riob(k)*(1.d0 - exp(-kb*dz))/(kb*dz)
           riob(k+1) = riob(k)*dexp(-kb*dz)
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccx
c       
	   do inb = 1,nbio
	      tempbio(inb)=bio(k,istep-1,inb)
	      dydtemp(inb)=0.d0
	   enddo
c
	   if(istep.eq.2) then
	      call derivs(tempbio,dydtemp,istep,k)
	      do inb=1,nsv
		 bio(k,istep,inb)=bio(k,istep-1,inb)
     &                      +delt*dydtemp(inb) !Euler step 
	      enddo
	      bio(k,istep,ipp)=delt*dydtemp(ipp)
              bio(k,istep,ichl)=dydtemp(ichl)         
	   else
              call derivs(tempbio,dydtemp,istep,k)
	      do inb=1,nsv
		 yt(inb)=bio(k,istep-1,inb)+delt
     &             *dydtemp(inb) 
	      enddo
	      call derivs(yt,dyt,istep,k) 

	      do inb=1,nsv
		 bio(k,istep,inb)=bio(k,istep-1,inb)+0.5d0*delt*
     &              (dydtemp(inb)+dyt(inb))
	      enddo
	      bio(k,istep,ipp)=0.5d0*delt*(dydtemp(ipp)+dyt(ipp))
              bio(k,istep,ichl)=0.5d0*(dydtemp(ichl)+dyt(ichl))
	   endif
	enddo

        return
        end 

C---------------------------------------------------------------------+
C       SUBROUTINE TRIDAG: Solves for a vector U of length N the      |
C       tridiagonal linear set given by the Crank-Nicholson scheme.    |
C       A, B, C and R are input vectors and are not modified.  This   |
C       subroutine is from Press et al., Numerical Recipes, p. 40.
C---------------------------------------------------------------------+
        subroutine tridag(a,b,c,r,u,n)
        implicit double precision (a-h,o-z)
        parameter (nmax=1000)
        real*8 gam(nmax), a(n), b(n), c(n), r(n), u(n), 
     &                bet,xu(n),bL(n)
c
        if(b(1).eq.0.d0)pause
         bet=b(1)
         xu(1)=r(1)/bet
         
         do j=2,n
            gam(j) = c(j-1)/bet
            bet=b(j)-a(j)*gam(j)
            if(bet.eq.0.d0)pause
            xu(j)=(r(j)-a(j)*xu(j-1))/bet
	 enddo
c
        bL(n)=xu(n)
        do j=n-1,1,-1
           bL(j)=xu(j)-gam(j+1)*xu(j+1)
	   xu(j)=bL(j)
	enddo

	do j=1,n
	   u(j)=bL(j)
	enddo

        return
        end
C---------------------------------------------------------------------+
C       SUBROUTINE SINK: 
C---------------------------------------------------------------------+
        subroutine sink(istep,inb)
        include 'common_blocks.h' 
c
        real*8 astar(nz),fdiff,fbelow,fabove,wnew
c
c       Calculate ASTAR using a simple upstream scheme at the top boundary: 
c
	wnew=wnsv(inb)/86400.d0
        fabove = 0.0d0         !No flux from above the top grid point
        fbelow = wnew*bio(1,istep,inb)*delt/dz  
        astar(1) = bio(1,istep,inb) - (fbelow - fabove) 
c
c       Calculate ASTAR values using a standard upstream scheme 
c
        do j=2,nz
           fabove = wnew*bio(j-1,istep,inb)*delt/dz
           fbelow = wnew*bio(j,istep,inb)*delt/dz
	   fdiff=fbelow-fabove
	   astar(j) = bio(j,istep,inb) - fdiff
	enddo
c
        do j=1,nz
	   bio(j,istep,inb)=astar(j)              
	enddo
c
        return
        end
C---------------------------------------------------------------------+
C       SUBROUTINE HORIZADV: 
C---------------------------------------------------------------------+
        subroutine horizadv(istep,inb)
c
        include 'common_blocks.h' 

	real*8 hadv

	do j=1,nz
	   if(HAnsv(inb).eq.2.) then 
	      hadv=HAAm(j,istep)
	   elseif(HAnsv(inb).eq.3.) then
	      hadv=HAFe(j,istep)
	   elseif(HAnsv(inb).eq.1.) then
	      hadv=HANi(j,istep)
	   endif
	   bio(j,istep,inb) = 
     &        bio(j,istep,inb)-hadv*bio(j,istep,inb)*delt
	enddo

	return
	end