subroutine readinit
	include 'common_blocks.h'

	real*8 depth,binitC,binitZ,binitN,binitD,binitF,binitS
     &          ,binitP,binitA

	if(ntb.eq.1) then
	   open (unit=30,status='old',file=datdir//'AS_Chl.init')
	   open (unit=31,status='old',file=datdir//'AS_Z.init')
	   open (unit=32,status='old',file=datdir//'AS_N.init')
	   open (unit=33,status='old',file=datdir//'AS_D.init')
	   open (unit=34,status='old',file=datdir//'AS_Fe.init')
	   open (unit=35,status='old',file=datdir//'AS_SiO4.init')
	   open (unit=36,status='old',file=datdir//'AS_PO4.init')
	elseif(ntb.eq.2) then
	   open (unit=30,status='old',file=datdir//'EqPac_Chl.init')
	   open (unit=31,status='old',file=datdir//'EqPac_Z.init')
	   open (unit=32,status='old',file=datdir//'EqPac_N.init')
	   open (unit=33,status='old',file=datdir//'EqPac_D.init')
	   open (unit=34,status='old',file=datdir//'EqPac_Fe.init')
	   open (unit=35,status='old',file=datdir//'EqPac_SiO4.init')
	   open (unit=36,status='old',file=datdir//'EqPac_PO4.init')
	   open (unit=37,status='old',file=datdir//'EqPac_A.init')
	endif
	do n=1,nz+1
	   read(30,*)depth,binitC
	   read(31,*)depth,binitZ
	   read(32,*)depth,binitN
	   read(33,*)depth,binitD
	   read(34,*)depth,binitF
	   read(35,*)depth,binitS
	   read(36,*)depth,binitP
	   if(ntb.eq.2)read(37,*)depth,binitA
c
cccc  USER ADD INITIAL CONDITIONS HERE  ccccccccccccccccccccccccc
c
cccc Mandatory initial conditions for workshop comparison: ccccc
c
c Single phytoplankton: bio(n,1,?)=binitC/n2chl
c Small phytoplankton: bio(n,1,?)=binitC*0.7d0/n2chl
c Large phytoplankton: bio(n,1,?)=binitC*0.3d0/n2chl
c Single zooplankton: bio(n,1,?)=binitZ
c Small zooplankton: bio(n,1,?)=binitZ*0.7d0
c Large zooplankton: bio(n,1,?)=binitZ*0.3d0
c Single detritus: bio(n,1,?)=binitD
c Small detritus: bio(n,1,?)=binitD
c Large detritus: bio(n,1,?)=0.d0
c Nitrate: bio(n,1,?)=binitN
c ntb=1 (Arabian Sea) 
c     Ammonium: bio(n,1,?)=0.d0
c ntb=2 (EqPac) 
c     Ammonium: bio(n,1,?)=binitA
c
cccc Optional (suggested) initial conditions: ccccccccccccccccc
c
c Tricho:  bio(n,1,?)=0.d0
c DON:     bio(n,1,?)=.15d0
c Fe:      bio(n,1,?)=binitF
c SiO4:    bio(n,1,?)=binitS
c PO4:     bio(n,1,?)=binitP
c
cccccccccccccccc Example: four-component model  ccccccccccccccc
c
c	   bio(n,1,1)=binitC/n2chl
c	   bio(n,1,2)=binitZ
c	   bio(n,1,3)=binitD
c	   bio(n,1,4)=binitN
c
cccccccccccccccc Example: eight-component model  cccccccccccccc
c
c	   bio(n,1,1)=binitC*0.7d0/n2chl
c	   bio(n,1,2)=binitC*0.3d0/n2chl
c	   bio(n,1,3)=binitZ*0.7d0
c	   bio(n,1,4)=binitZ*0.3d0
c	   bio(n,1,5)=binitD
c	   bio(n,1,6)=0.d0
c	   bio(n,1,7)=binitN
c	   if(ntb.eq.1) bio(n,1,8)=0.d0
c	   if(ntb.eq.2) bio(n,1,8)=binitA
c
cccc Optional aeolian flux in units of nmol Fe m^-2 d^-1: cccccccc
c
ccccccccc  not provided and probably not necessary ccccccccccccccc
c
	   do inb=1,nsv
	      if(aeonsv(inb).eq.1) then
		 write(*,*)'Read Aeolian Flux Data in readinit.f

c       open(unit=38,status='old',file=datdir//'AeolianFeFlux')
c       do jmain=1,isteps
c          read(38,*)aeoflux(jmain)
c       enddo

	      endif
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
	enddo
	
	close(30)
	close(31)
	close(32)
	close(33)
	close(34)
	close(35)
	close(36)
	if(ntb.eq.2)close(37)
c
        return
        end