C	Come here to calculate the derivatives:
C---------------------------------------------------------------------+
C	SUBROUTINE DERIVS_4:   This subroutine calculates time change |
C	in each biological compartment. Note that arrays y(nbio) and  |
C	dydtt(nbio) are dummy arrays which are passed different       |
C	values of bio and dydt depending on the step size in question.|
C	This version is for the 4-component (McCreary) ecosystem.     |
C---------------------------------------------------------------------+
	subroutine derivs(y,dydtt,istep,k)

	include 'common_blocks.h'

	real*8 y(nbio), dydtt(nbio)

C-----------------------------------------------------------+
C	Matrix configuration for the 4-component ecosystem: |
C	[1]: P                                              |
C	[2]: Z                                              |
C	[3]: D                                              |
C	[4]: N                                              |
C-----------------------------------------------------------+
C-----------------------------+
C	Local declarations    |
C-----------------------------+

	real*8 nlim, phar, fi, fdenom, sfi, spi, szi, q02
c	
C	Documentation for this set of ecosystem equations is 
C	available on the web at the following address:
C	http://www.nova.edu/ocean/arabsea/testbed/biolayermodel.pdf

        f0=1.0d0*1.0
        deathz=3.d-7
        az=0.1d0*.5
        rphip=0.83333333d0
        phiz=0.1666666d0
        e=1.3*reminrate/86400.d0
        pgrowth=2.8935d-5*1.5
        an=0.4d0
        deathp=1.1574d-6*0.
        grazz=4.6296d-05*3.0
        q0=40.d0*1.5
        rnt0=1.d0*.2

C	Calculate terms used in JPM ecosystem.

	nlim=y(4)/(y(4)+rnt0)		! Nutrient limitation

	fi = rphip*y(1)+phiz*y(2)	! Total food supply for Z
	fdenom = fi + f0
	sfi = fi/fdenom			! Grazing limitation

	spi = rphip*y(1)/fdenom		! Z herbivory
	szi = phiz*y(2)/fdenom		! Z carnivory

C	Light at depth. Note, rlt (in W/m**2) is passed from the 
C	subroutine CALCLIGHT.

	q02=q0*q0
	phar = rlt/((rlt*rlt+q02)**0.5d0)	! light limitation

C--------------------------------------+
C	Primary Production equation    |
C--------------------------------------+

	dydtt(ipp)=pgrowth*phar*nlim*y(1)*6.625d0*12.d0

C---------------------------------+
C	Chlorophyll equation      |
C---------------------------------+

	dydtt(ichl)=y(1)*n2chl

C---------------------------------+
C	Phytoplankton equation    |
C---------------------------------+

	dydtt(ipp)=pgrowth*phar*nlim*y(1)*6.625d0*12.d0
	dydtt(1) = pgrowth * phar * nlim * y(1)
     .	- grazz * spi * y(2)
     .  - deathp * y(1)

C---------------------------------+
C	Zooplankton equation      |
C---------------------------------+
	dydtt(2) = az * grazz * sfi * y(2)
     . - grazz * szi * y(2)
     . - deathz * y(2)

C---------------------------------+
C	Detritus equation         |
C---------------------------------+
	dydtt(3) = (1.d0 - az - an) * grazz * sfi * y(2)
     . + deathp * y(1)
     . + deathz * y(2)
     . - e * y(3)

C---------------------------------+
C	Nitrogen equation         |
C---------------------------------+
	dydtt(4) = an * grazz * sfi * y(2)
     . + e * y(3)
     . - pgrowth * phar * nlim * y(1)

	return
	end