subroutine calccost

	include 'common_blocks.h'

	real*8 totdat,totN,totZ,totpp
        real*8 Dflux(500),TavgD400(500)
	real*8 davg1,davg2,davg3,davg4,ST1,ST2,ST3,ST4
	tms=24.d0
	if(ntb.eq.1) then
	   open(unit=16,file=homedir//'cost_as.out')
	elseif(ntb.eq.2) then
	   open(unit=16,file=homedir//'cost_ep.out')
	endif
c       
ccccccc   COMPUTE model equivalents of SED FLUX cccccccccccccccc
c
	if(ntb.eq.1) then 
	   do nscount=1,nSTdat
	      Dflux(nscount)=0.d0
	   enddo
	   nscount=1
	   do jmain=1,isteps
	      if(dayST(nscount).ge.cdays(jmain).
     &          and.dayST(nscount).lt.cdays(jmain+1)) then
		 do inb=1,nsv
		    if(Cnsv(inb).eq.3.) then
		       tmpD400=0.d0
		       do nt=1,nint(timeST(nscount)*tms)
			  tmpD400=tmpD400+bio(nz,jmain+nt-1,inb)
     &                      *unitnsv(inb)
		       enddo
		       TavgD400(nscount)=14.d0*
     &                  tmpD400/(real(nint(timeST(nscount)*tms)))
		       Dflux(nscount)=Dflux(nscount)+
     &                 wnsv(inb)*TavgD400(nscount)*
     &                 dexp((-zST(nscount)+real(nz*dz))*
     &                 reminrate/wnsv(inb))
		    endif
		 enddo
		 nscount=nscount+1
	      endif
	   enddo
	elseif(ntb.eq.2) then
	   davg1=0.d0
	   davg2=0.d0
	   davg3=0.d0
	   davg4=0.d0
c Shallow EqPac trap data (units mmolC/m2/day)
	   ST1=5.3d0
	   ST2=6.7d0
	   ST3=12.1d0
	   ST4=12.8d0
	   do inb=1,nsv
	      if(Cnsv(inb).eq.3) then
		 do j=2784,2784+30*24-1
		    davg1=davg1+bio(12,j,inb)*wnsv(inb)
     &                    *(117./16.)*unitnsv(inb)
		 enddo
	      endif
	   enddo
	   davg1=davg1/(24.*30.)
	   do inb=1,nsv
	      if(Cnsv(inb).eq.3) then
		 do j=3648,3648+30*24-1
		    davg2=davg2+bio(12,j,inb)*wnsv(inb)
     &                    *(117./16.)*unitnsv(inb)
		 enddo
	      endif
	   enddo
	   davg2=davg2/(24.*30.)
	   do inb=1,nsv
	      if(Cnsv(inb).eq.3) then
		 do j=7248,7248+30*24-1
		    davg3=davg3+bio(12,j,inb)*wnsv(inb)
     &                *(117./16.)*unitnsv(inb)
		 enddo
	      endif
	   enddo
	   davg3=davg3/(24.*30.)
	   do inb=1,nsv
	      if(Cnsv(inb).eq.3) then
		 do j=8280,8280+30*24-1
		    davg4=davg4+bio(12,j,inb)*wnsv(inb)
     &                *(117./16.)*unitnsv(inb)
		 enddo
	      endif
	   enddo
	   davg4=davg4/(24.*30.)
	endif
c
ccccccccccccccccc  CALCULATE COST FUNCTION *** THREE TIMES *** ccccc
c
	do icalccost=1,3
	   if (icalccost.eq.1.and.ntb.eq.1) then
	      w(1)=1.d0/.06d0	         ! weight for chlorophyll
	      w(2)=1.d0/.6d0	         ! weight for nitrate+ammonium
	      w(3)=1.d0/(dlog10(1.25d0)) ! weight for sediment Trap Flux
	      w(4)=1.d0/(0.3d0)          ! weight for zooplankton
	      w(5)=1.d0/5.d0	         ! weight for primary production
	   elseif (icalccost.eq.1.and.ntb.eq.2) then
	      w(1)=1.d0/.06d0	         ! weight for chlorophyll
	      w(2)=1.d0/.6d0	         ! weight for nitrate+ammonium
	      w(3)=1.d0/1.25d0           ! weight for sediment Trap Flux
	      w(4)=1.d0/(0.3d0)          ! weight for zooplankton
	      w(5)=1.d0/5.d0	         ! weight for primary production
	   elseif(icalccost.eq.2.and.ntb.eq.1) then
	      w(1)=1.d0/0.2d0	         ! weight for chlorophyll
	      w(2)=1.d0/2.4d0	         ! weight for nitrate+ammonium
	      w(3)=1.d0/0.3d0            ! weight for sediment Trap Flux
	      w(4)=1.d0/0.4d0            ! weight for zooplankton
	      w(5)=1.d0/19.2d0	         ! weight for primary production
	   elseif(icalccost.eq.2.and.ntb.eq.2) then
	      w(1)=1.d0/0.09d0	         ! weight for chlorophyll
	      w(2)=1.d0/2.2d0	         ! weight for nitrate+ammonium
	      w(3)=1.d0/3.8d0            ! weight for sediment Trap Flux
	      w(4)=1.d0/0.2d0            ! weight for zooplankton
	      w(5)=1.d0/8.5d0	         ! weight for primary production
	   elseif(icalccost.eq.3.and.ntb.eq.1) then
	      w(1)=1.d0/0.2d0	         ! weight for chlorophyll
	      w(2)=1.d0/2.4d0	         ! weight for nitrate+ammonium
	      w(3)=1.d0/0.3d0            ! weight for sediment Trap Flux
	      w(4)=1.d0/3.0d0            ! weight for zooplankton
	      w(5)=1.d0/19.2d0	         ! weight for primary production
	   elseif(icalccost.eq.3.and.ntb.eq.2) then
	      w(1)=1.d0/0.09d0	         ! weight for chlorophyll
	      w(2)=1.d0/2.2d0	         ! weight for nitrate+ammonium
	      w(3)=1.d0/3.8d0            ! weight for sediment Trap Flux
	      w(4)=1.d0/3.0d0            ! weight for zooplankton
	      w(5)=1.d0/8.5d0	         ! weight for primary production
	   endif

	costP=0.d0
	costZ=0.d0
	costN=0.d0
	costST=0.d0
	costPP=0.d0
	if(ntb.eq.1)aNp=48.d0/36.d0
	if(ntb.eq.2)aNp=356.d0/321.d0
	if(ntb.eq.1)aNz=48.d0/48.d0
	if(ntb.eq.2)aNz=356.d0/326.d0
	if(ntb.eq.1)aNpr=48.d0/34.d0
	if(ntb.eq.2)aNpr=356.d0/321.d0
	if(ntb.eq.1)aNn=48.d0/37.d0
	if(ntb.eq.2)aNn=356.d0/356.d0
	if(ntb.eq.1)aNst=48.d0/30.d0
	if(ntb.eq.2)aNst=356.d0/4.d0
	nzcount=1
	nccount=1
        nncount=1
        nscount=1
        npcount=1
	do jmain=1,isteps
	   istep = jmain 
           do n=1,nz
              if(dayPrPr(npcount).ge.cdays(jmain).
     &              and.dayPrPr(npcount).lt.cdays(jmain+1)) then
                 if(zPrPr(npcount).ge.dz*(real(n)-1.d0).
     &              and.zPrPr(npcount).lt.(dz*real(n)))then
		    totpp=0.d0
		    do ii=1,24
		       totpp=totpp+bio(n,istep+ii-12,ipp)
		    enddo
		    costPP=costPP+aNpr*(w(5)*(aPrPr(npcount)
     &                   -totpp))**2.d0
                    npcount=npcount+1
                 endif
              endif
           enddo
	   do n=1,nz
	      if(dayZ(nzcount).ge.cdays(jmain).
     &              and.dayZ(nzcount).lt.cdays(jmain+1)) then
		 if(zZ(nzcount).ge.dz*(real(n)-1.d0).
     &                   and.zZ(nzcount).lt.(dz*real(n)))then
		    totZ=0.d0
		    do inb=1,nsv
		       if(Cnsv(inb).eq.2.) 
     &                   totZ=totZ+bio(n,istep,inb)*unitnsv(inb)
		    enddo
		    if(totz.ne.0.d0)
     &                costZ=costZ+aNz*(w(4)*(aZoo(nzcount)-totZ))**2.d0
		    nzcount=nzcount+1
		 endif
	      endif
	   enddo
	   do n=1,nz
	      if(dayChl(nccount).ge.cdays(jmain).
     &              and.dayChl(nccount).lt.cdays(jmain+1)) then
		 if(zChl(nccount).ge.dz*(real(n)-1.d0).
     &                   and.zChl(nccount).lt.(dz*real(n)))then
		    costP=costP+aNp*(w(1)*(Chl(nccount)
     &                          -bio(n,istep,ichl)))**2.d0
		    nccount=nccount+1
		 endif
	      endif
	   enddo
	   do n=1,nz
	      if(dayN(nncount).ge.cdays(jmain).
     &              and.dayN(nncount).lt.cdays(jmain+1)) then
		 if(zN(nncount).ge.dz*(real(n)-1.d0).
     &                   and.zN(nncount).lt.(dz*real(n))) then
		    totN=0.d0
		    do inb=1,nsv
		       if(Cnsv(inb).eq.1.) 
     &                      totN=totN+bio(n,istep,inb)*unitnsv(inb)
		    enddo
		    costN=costN+aNn*(w(2)*(aNut(nncount)-totN))**2.d0
		    nncount=nncount+1
		 endif
	      endif
	   enddo
	   if(ntb.eq.1) then 
	      if(dayST(nscount).ge.cdays(jmain).
     &             and.dayST(nscount).lt.cdays(jmain+1)) then
		 costST=costST+aNst*(w(3)*(dlog10(ST(nscount))-
     &             dlog10(Dflux(nscount))))**2.d0
		 nscount=nscount+1
	      endif          
	   endif
	enddo
c	
	if(ntb.eq.2) then
	   costST=aNst*( (w(3)*(ST1-davg1))**2.d0 
     &                    + (w(3)*(ST2-davg2))**2.d0 
     &                    + (w(3)*(ST3-davg3))**2.d0 
     &                    + (w(3)*(ST4-davg4))**2.d0 )
	endif

	if(ntb.eq.1)totdat=48*5
	if(ntb.eq.2)totdat=356*5
	write(16,*)'Total cost = ', 
     &         sngl((costP+costZ+costN+costST+costPP)/totdat)
	write(16,*)'P cost = ',sngl(costP/totdat)
	write(16,*)'Z cost = ',sngl(costZ/totdat)
	write(16,*)'N cost = ',sngl(costN/totdat)
	write(16,*)'ST cost =',sngl(costST/totdat)
	write(16,*)'PP cost =',sngl(costPP/totdat)
	write(*,*)'Total cost = ', 
     &         sngl((costP+costZ+costN+costST+costPP)/totdat)
	write(*,*)'P cost = ',sngl(costP/totdat)
	write(*,*)'Z cost = ',sngl(costZ/totdat)
	write(*,*)'N cost = ',sngl(costN/totdat)
	write(*,*)'ST cost =',sngl(costST/totdat)
	write(*,*)'PP cost =',sngl(costPP/totdat)

	enddo

	end