!  larval_cod_SDP.f90 
!
!  FUNCTIONS:
!	larval_cod_SDP      - Entry point of console application.
!

!****************************************************************************
!
!  PROGRAM: Fortran F90 code. Written by Øyvind Fiksen (University of Bergen, Norway) 
!           and Christian Jørgensen, Uni Research, Norway. email:oyvind.fiksen@bio.uib.no)
!
!  PURPOSE:  Model of fish habitat selection, including stomach dynamics
!
!****************************************************************************

	program larval_cod_SDP
	!use AVVIEWER
	!use AVDEF
	implicit none

!Program execution
integer, parameter :: RunOptimization = 1   !1: include SDP, 0: skip SDP and read files

!Model resolution
integer, parameter :: Lmin = 5, Lmax = 20   !Length limits (mm) for SDP solution 
integer, parameter :: Smax = 10				!Number of stomach filling categories (+1)
integer, parameter :: zmax = 70, dz = 2		!MaxDepth, vertical resolution
integer, parameter :: MaxAct = 10, dAct = 2	!Activity
integer, parameter :: DayMax = 100			!Number of day iterations for fitness stabilization
integer, parameter :: Horizon = 24			!The number of time steps throughout a day
integer, parameter :: dt = 3600

!Iterations
integer, parameter :: MultSim = 0			!If 0, then multiple loops
integer, parameter :: MultNum = 30			!Number of iterations (30) in MultSim mode, =1 (or 0) in single mode
integer, parameter :: dSimNum = 1
!integer, parameter :: lesfil = 1			!Read (1) or make (0) binfiles of z and act
real*8,	 parameter :: MultRes = 1.			!Range in Multsim=MultNum*MultRes

!Parameters
real*8,  parameter :: PredationCostOfActivity = 1.	!Scaling predation effects of activity
real*8,  parameter :: AssimilationEfficiency=0.75
real*8,	 parameter :: PreyDensDryMass = 10., PDmean_z = 30., PDsd = 10.	!mg dw prey mass per m3, depth distribution
real*8,	 parameter :: FishDens = 1d-4				!number of fish predators per m3
real*8,  parameter :: GutSize = 0.06				!Gut capacity as a fraction of body mass (may be a function..)
real*8,	 parameter :: aPred = 2.5d-9, bPred = -0.85 !(2.78d-5 & -1.3 on Length! in Fiksen&al 2002. 
													!For McGurk it is a=2.2X10-4/86400 and b=-0.85 on g dw!
real*8,	 parameter :: Ambush = 0.5					!Fraction of invertebrate predation that is ambush (0.5?)
real*8,  parameter :: S_mortality = 0.d-6			!Starvation rate (s-1) at S=0
real*8,	 parameter :: SwimLim = 1.					!Migration range in fractions of body lengths per second
real*8,	 parameter :: SwimCost = 1.					!Scaling of migration costs in terms of basic metabolism when swimming SwimLim BL/dt
real*8,  parameter :: MAXLIG=500.,JulianDay=100.,B_latitude=70.		!Diel light cycle
real*8,  parameter :: MixedLayerChla = 2., MixedLayerDepth = 70.	!Chlorophyll a in mixed layer
real*8,  parameter :: MixedLayerTemp = 8.,MixedLayerSlope = 20.		!Mixed layer temp


!Variables
!Optimal depth decision and reproductive value matrixes:
real*8,  dimension(Lmin:Lmax,0:Smax,0:Zmax,0:Horizon) :: zOpt, V, ActOpt, zOpt_check,ActOpt_check

real*8,  dimension(Lmin:Lmax,0:Smax,0:Zmax) :: Vz
real*8,	 dimension(0:zmax,0:Horizon) :: Eb, PME, FE !Light, prey- and pred enc matrixes
real*8,  dimension(0:zmax) :: Chl_a, BeamAtt, PreyDensDM, Temp
integer  Day, L, S, t, MultVar					!Counters for different loops
integer  Depth, Dep, z, OptimalDepth, Activity
integer  z_up, z_down, dzswim, MeanDepth
integer  intL, intS								!integer part of length for interpolations
real*8   newL, newS								!length and S next time-step
real*8   dL, dS									!decimal part of length for interpolations
real*8   L_m									!length in meters
real*8	 act, OptimalAct, PauseDuration, swimDist
real*8	 Stomach
real*8   InitialS, delta
real*8   Fitness, F_Max
real*8   PreyMassEnc, PD
real*8   Weight, DeltaMass
real*8   Survival, Amb
real*8	 FishEnc, FishD
real*8	 zdiff,actdiff,zdiff1,actdiff1
real*8   MLDep,MLChla,MLTemp,MLTempSlope

character(12) filename

!Functions
real*8   SGR			!Specific growth rate from Folkvord model
real*8   R				!Respiration 30(
real*8   Length_Weight	!Length as a function of weight
real*8	 Weight_Length	!Weight as a function of length
real*8	 m_starv		!Mortality rate when S = 0
real*8	 m_OtherPred	!Predators other than fish

open(2, file = 'parameters.txt')
write(2,'(13I6,20E12.3)')Lmin,Lmax,Smax,zmax,dz,MaxAct,dAct,DayMax,Horizon,dt,MultSim,MultNum,dSimNum, &
	& AssimilationEfficiency,PreyDensDryMass,PDmean_z,PDsd,FishDens,GutSize,aPred,bPred,Ambush,S_mortality, &
	& SwimLim,SwimCost,MAXLIG,JulianDay,B_latitude,MixedLayerChla,MixedLayerDepth,MixedLayerTemp,MixedLayerSlope
close(2)

do Multvar = 26,MultNum,dSimNum
  PD = PreyDensDryMass	
  FishD = FishDens						!If only one sim
  MLChla = MixedLayerChla
  MLDep = MixedLayerDepth
  Amb = Ambush
  MLTemp = MixedLayerTemp
  MLTempSlope = MixedLayerSlope
	
 if(MultSim==0) then			 
  PD = dble(Multvar*MultRes*1.)				!If many sims - set the one relevant
  !FishD = MultVar*FishDens
  !MLChla = MultVar*MixedLayerChla
  !MLDep = dble(Multvar)
  !Amb =dble(Multvar)*0.1
  !MLTemp = MixedLayerTemp + dble(Multvar)
 endif

call SetEnvironment(MLDep,MLChla,MLTemp,MLTempSlope,Horizon,zmax,MAXLIG,JulianDay,&
& dz,B_latitude,PD,PDmean_z,PDsd,PreyDensDM,Chl_a,BeamAtt,Eb,Temp)

!stop

Print*
print('(A30,F10.4)'), 'prey abundance:', PreyDensDM(0)
print('(A30,E10.2)'), 'Fish abundance:', FishD
print('(A30,F10.4)'), 'Chla concentration:', MLChla
print('(A30,F10.4)'), 'Mixed layer depth:', MLDep
print('(A30,F10.4)'), 'Fraction ambush pred:', Amb

zDiff1=1000000.
actDiff1=100000.

if (RunOptimization .eq. 1) then  !Kjør SDP

!Solve the SDP equation:
print*
print*, 'Solving SDP..'

do Day = 0,DayMax						!Run enough days to avoid terminal effects
 if (Day .eq. 0) then
  zOpt(:,:,:,Horizon) = 0				!Reset
  V(:,:,:,Horizon) = 0.					!Reset
 else
  V(:,:,:,Horizon) = V(:,:,:,0)			!Fitness transferred between days - iterations will stabilize
  zOpt(:,:,:,Horizon) = zOpt(:,:,:,0)	!Optimal depth between days
 endif

 do t = Horizon-1, 0, -1   !Cycle over all timesteps in a day (hours or similar)
  V(Lmax,:,:,t+1) = 1.

  do L = Lmin, Lmax		!Cycle over length categories
   L_m = dble(L)*1.d-3		!Larval length in m
   !Weight = f(L)		!Expected body mass

   do S = 0, Smax		!Stomach fullness (number of fullness categories)
    InitialS = (dble(S)/dble(Smax))*GutSize*Weight_Length(L_m)  !Initial stomach content in mg dWeight

    do Dep = Zmax,0,-dz			!Depth loop
	   F_max = 0.					!Reset fitness maximum

	   !Behavioural options:
    do Activity = 0,MaxAct,dAct
     act = max(0.,min(dble(Activity)/dble(MaxAct),1.)) !Activity 0-1 scale

     dzswim = dz*idnint(act*SwimLim*L_m*dt/dble(dz))  !max vertical migration rounded off to nearest depth category
      z_up = max(dep-dzswim,0)		!Limit migration to nearby habitats (z(t+1)+-dz), always possible to migrate at least one step
	  z_down= min(dep+dzswim,Zmax)	!that can be reached by directional swimming
	    do z = z_up,z_down,dz		!Depth range (limited by previous)
!			swimDist = abs(dble(Dep-z))/dble(SwimLim*L_m*dt) !Swimming relative to max..

			MeanDepth = idnint(dble(Dep+z)/2.)
			!Mass of prey ingested over the timestep
			call ingestion(act,Eb(MeanDepth,t),BeamAtt(MeanDepth),PreyDensDM(MeanDepth),L_m,PauseDuration,PreyMassEnc)
			delta = dt*SGR(Weight_Length(L_m),Temp(MeanDepth))*Weight_Length(L_m)/AssimilationEfficiency       !Stomach content required...
			delta = delta + dt*R(Temp(MeanDepth),Weight_Length(L_m))/AssimilationEfficiency !...to sustain max growth
     
			if(delta .le. InitialS)then	!Growth if stomach can supply enough mass
				Weight = Weight_Length(L_m)*(1.+ (exp(dt*SGR(Weight_Length(L_m),Temp(MeanDepth)))-1.))  
				Stomach = InitialS-delta
			else							!New weight if stomach is limiting
				Weight = Weight_Length(L_m) + InitialS*AssimilationEfficiency - dt*R(Temp(MeanDepth),Weight_Length(L_m))  
				Stomach = 0.
			endif

			Stomach = min(Stomach + PreyMassEnc*dt,GutSize*Weight_Length(L_m))  !Stomach filling, constrained
			Stomach = max(0.,Stomach)

			Weight = Weight - act*dt*R(dble(7),Weight_Length(L_m)) !Energetic activity cost

			call FishPred(Eb(MeanDepth,t),BeamAtt(MeanDepth),FishD,L_m,FishEnc)	!Expected encounters with predators

			!Local predation rate from fish, other predators and starvation
			Survival = Amb*(PredationCostOfActivity*act)*m_OtherPred(Weight_Length(L_m),aPred,bPred) + (1.-Amb)*m_OtherPred(Weight_Length(L_m),aPred,bPred) !Inv pred during feeding
		    !Survival = Survival + Amb*swimDist*m_OtherPred(Weight_Length(L_m),aPred,bPred)	!Invertebrate predation during migration
			Survival = Survival + FishEnc + m_starv(dble(S),S_mortality)
			Survival = exp(-Survival*dt)  !Survival probability in dt

			!State interpolations
			newL = 1d3*Length_Weight(Weight)					!New length in state intervals (mm?)  
			intL = max(Lmin,min(idint(newL),Lmax-1))      		!Interpolation - stay within array
			dL   = max(0.,min(newL-dble(intL),1.))	      		!Interpolation - stay within array

			newS = dble(Smax)*Stomach/(GutSize*Weight_Length(L_m))  !New stomach state
			intS = max(0,min(idint(newS),Smax-1))
			dS   = max(0.,min(newS - dble(intS),1.))			!Interpolation

			!The SDP equation
			Fitness = (1.-dL)*(1.-dS)*V(intL  ,intS  ,z,t+1) + &   !Length down, stomach down
                (1.-dL)*    dS *V(intL  ,intS+1,z,t+1) + &   !Length down, stomach up
                    dL *(1.-dS)*V(intL+1,intS  ,z,t+1) + &   !Length up, stomach down
                    dL *    dS *V(intL+1,intS+1,z,t+1)       !Length up, stomach up 
      
			Fitness = Fitness * Survival		!Survival probability

			if(Fitness .gt. F_max) then		!Pick the optimal depth
			 F_max = Fitness
			 OptimalDepth = z
     		 OptimalAct = act
			endif

     enddo      !Loop over available depths
	   enddo	   !Activity loop

    V(L,S,Dep,t) = F_max
    zOpt(L,S,Dep,t) = OptimalDepth
	ActOpt(L,S,Dep,t) = OptimalAct

   enddo !All depths
  enddo !Stomach fullness
 enddo !Length classes
   
enddo !Day

 zDiff = sum(dabs(zOpt-zOpt_check))				!Difference between days in behavioural strategies 
 zOpt_check = zOpt
 actDiff = sum(dabs(ActOpt-ActOpt_check))		!Difference between days in behavioural strategies 
 ActOpt_check = ActOpt
 
 print('(A6,I4,5E12.4)'),'Day=', day, zDiff,actDiff,V(Lmin,Smax,10,Horizon),V(Lmax-5,Smax,10,Horizon)

 if((zDiff==0).and.(actDiff==0))exit				!Exit difference between days is zero
 zDiff1=zDiff
 actDiff1=actDiff

enddo !Day

!store optimal strategies
write(filename, '(a8,f4.1)') 'PreyDens',PD
open(1, file = filename//'zopt.bin', form='binary')
write(1)zopt
close(1)
open(1, file = filename//'actopt.bin', form='binary')
write(1)actopt
close(1)
open(1, file = filename//'V.bin', form='binary')
write(1) V
close(1)

!*******
!or read from  optimal strategies from .bin files:
elseif (RunOptimization .eq. 0) then
 write(filename, '(a8,f4.1)') 'PreyDens',PD
 open(1, file =   filename // 'zopt.bin', form='binary')
 read(1)zopt
 close(1)
 open(1, file = filename//'actopt.bin', form='binary')
 read(1)actopt
 close(1)
 open(1, file = filename // 'V.bin', form='binary')
 read(1)V
 close(1)
endif

call LarvalSDPForward(zOpt,ActOpt,V,Eb,PME,FE,Lmin,Lmax,Smax,zmax,DayMax,MaxAct,Horizon,dt,dz,Gutsize,AssimilationEfficiency, &
& BeamAtt,PreyDensDM,PD,Amb,aPred,bPred,S_mortality,FishD,Chl_a,Temp,MultSim,SwimCost,SwimLim,Day,zDiff,actDiff,PredationCostOfActivity)

enddo !Prey loop
!if((MultSim==0).and.(PreyDensDM.lt.5))goto 111

end program larval_cod_SDP

!*****************
function SGR(Weight,Temp) !Arild Folkvord's growth model, no food lim, valid for Weight [0.03-350] mgdw
 implicit none
 
 real*8 Temp,Weight,SGR
 if(Weight<0.03)print*,'Weight below limit of 0.03 mgdwgt!',Weight
 if(Weight>350)print*,'Weight above limit of 350 mgdwgt!',Weight

 SGR = 1.2 + 1.8*Temp - 0.078*Temp*dlog(Weight) - 0.0946*Temp*(dlog(Weight))**2. &
    + 0.0105*Temp*(dlog(Weight))**3.
 SGR = dlog(SGR/100. + 1)/86400. !Specific temperature-limited growth rate g/g/dt (in sec; /sec)

end function SGR

!**********
function R(Temp,Weight)  !Respiration in mass per second
 implicit none
 !Inherited values
  real*8 Weight, Temp
 !Parameters
  real*8,  parameter :: R0 = 2.38E-7, Finn_b=0.9, T_resp=0.088 !Parameters for respiration - Finn et al. 2002
 !Returned value
  real*8 R

 R  = R0*dexp(T_resp*Temp)   !Respiration from Finn & al. (2002), g dWeight/dt
 R  = (Weight**Finn_b)*R
end function R

function Length_Weight(Weight)  !length in m as from dry body mass in mg
 implicit none
 !Inherited values
 real*8 Weight
 !Returned value
 real*8 Length_Weight
 !Function from Otterlei & al 1999
 Length_Weight = dexp((dlog(Weight)+9.2)/3.9)
 Length_Weight = 1.d-3*Length_Weight  !From mm to m
end function Length_Weight    

function Weight_Length(L_m)  !Dry body mass (mg) from length in m
 implicit none
 !Inherited values
 real*8 L_m
 !Returned value
 real*8 Weight_Length
 !Function from Otterlei & al 1999
 Weight_Length = dexp(-9.2 + 3.9*dlog(1.d3*L_m))
 end function Weight_Length   
 
function m_starv(Stomach,S_mortality) !Death risk from starvation
 implicit none
 !Inherited values
 real*8 Stomach,S_mortality
 !Return value
 real*8 m_starv
 if(Stomach==0)then
  m_starv = S_mortality
 else
  m_starv = 0.
 endif
 end function m_starv


function m_OtherPred(Weight,aPred,bPred) !Death risk from other sources
 implicit none
 !Inherited values
 real*8 Weight,aPred,bPred
 !Return value
 real*8 m_OtherPred
 m_OtherPred = aPred*(Weight*1.d-3)**bPred
 !m_OtherPred=1.2E-6 !Fixed 0.01 death rate per hour (no size-dependence).
 end function m_OtherPred

 
 
subroutine SetEnvironment(MLDep,MLChla,MLTemp,MLTempSlope,Horizon,zmax,MAXLIG,JulianDay,&
& dz,B_latitude,PD,PDmean_z,PDsd,PreyDensDM,Chl_a,BeamAtt,Eb,Temp)
implicit none

integer Horizon,t,zmax,z,dz,L
real*8  MAXLIG,JulianDay,B_latitude,PD,BeamAtt(0:zmax),SLIG,Eb(0:zmax,0:Horizon)
real*8  Chl_a(0:zmax),DiffAtt(0:zmax),PreyDensDM(0:zmax),Temp(0:zmax)
!real*8, parameter :: k0 = 0.14, BeamDiffCoefficient = 3.
real*8  MLDep,MLChla,MLTemp,MLTempSlope,PreyPer_m3
real*8	PDsd,PDmean_z,normal

open(28, file = 'PreyDens_z.txt')
open(30, file = 'DepthEnv.txt')

!Set Chl_a profile
do z = 0,zmax

 if(dble(z)<=MLDep) then
  Chl_a(z) = MLChla
  PreyDensDM(z) = PD
  Temp(z) = MLTemp  !- 0.05*dble(z)
 else
  Chl_a(z) = 0.
  PreyDensDM(z) = PD
  Temp(z) = MLTemp  !- 0.05*dble(z) !-MLDep)
 endif
 !To generate a normal prey distribution:
PreyDensDM(z)= PD*normal(PDsd,PDmean_z,z)/normal(PDsd,PDmean_z,nint(PDmean_z))
Diffatt(z) = 0.0567 + Chl_a(z)*(0.0506*exp(-0.606*Chl_a(z))+0.0285)+(0.068*exp((-0.014)*(550-550))) !Coastal type 2, equations 13 and 15 in review
BeamAtt(z) = 0.0657 + (0.39*(Chl_a(z)**0.57))*(1.563 - 0.001149*550) ! equations 3.47 and 3.48 in Mobley (Voss 1992)

enddo

!Set light profile in depth and time
do t = 0,Horizon
 call SURLIG(dble(t),MAXLIG,JulianDay,B_latitude,SLIG)
 Eb(0,t) = SLIG !max(SLIG,dble(1E-4)) !Min surface light at night

do z = 1,zmax
 Eb(z,t) = Eb(z-1,t)*exp((-DiffAtt(z))) !
enddo
enddo

do z = 0,Zmax,dz
   Write(30,'(100E12.3)') dble(z),Temp(z),Chl_a(z),PreyDensDM(z),Eb(z,12),Eb(z,6),Eb(z,0)
enddo
close(30)

!Print the prey field in numbers per L over depth and for all size classes
do L=5,15
do z = 0,zmax
 call preydens(dble(0),Eb(0,0),BeamAtt(z),PreyDensDM(z),dble(0.001*L),dble(2),PreyPer_m3)
 write(28,'(100E12.3)')dble(L),dble(z),PreyPer_m3,PreyPer_m3*1E-3
enddo
enddo
close(28)

return
end subroutine


 REAL*8 FUNCTION normal(sd,mean_depth,depth)

 IMPLICIT NONE 	
 REAL*8, INTENT(in) :: sd               !standard deviation of normal distribution
 REAL*8, INTENT(in) :: mean_depth		!mean of normal distribution
 INTEGER, INTENT(in) :: depth			!depth 
 
 REAL*8 expon,pi
 pi = acos(-1.)

 !Normal distribution = 1/(sd*sqrt(2*pi)) * exp(-(x-mean)**2/2*sd**2)
 expon = -(((dble(depth) - mean_depth))**2)/(2.*(sd**2))
 normal = exp(expon)/(sqrt(2.*pi)*sd)
 
 END FUNCTION normal 

 
 
!Just to get prey density in numbers:
subroutine preydens(act,Eb,BeamAtt,Prey_mgdw_per_m3,Larval_m,PauseDuration,PreyPer_m3)
implicit none

!Parameters
real*8, parameter:: PreyRelLength=0.05
real*8, parameter:: PreyWidth = 0.5					!Prey diam:length ratio (for mass)
real*8, parameter:: PreyContrast=0.3
real*8, parameter:: VisRtoSDPreyinBL=0.8
real*8, parameter:: Ke = 1.
real*8, parameter:: omega = 0.
real*8, parameter:: PreySwimVel = 0.
real*8,  parameter :: Maxpause = 10.
real*8,  parameter :: MinPause = 2.
!real*8, parameter:: PauseFrequency = 0.43			!
!real*8, parameter:: PauseDuration = 2.				!Seconds per pause

real*8, parameter:: VisFieldShape = 0.5
real*8, parameter:: zp_mgdw_per_m3 = 1025*0.25*1.d6

!Variables
integer	IER
real*8 Eb,BeamAtt,Prey_mgdw_per_m3,Larval_m,PreyMassEnc
real*8 Prey_m,Prey_mgdw,PreyPer_m3,PreyImageArea,Swim_m_per_s
real*8 VisualRange,EyeSens,pi
real*8 act,PauseFrequency,PauseDuration

!Functions
real*8 Weight_Length
real*8 Length_Weight

pi=acos(-1.)

!Prey properties
Prey_m = Larval_m*PreyRelLength	!Prey length (prosome?)
!Prey_mgdw = (4/3.)*pi*((0.5*Prey_m)**3.)*zp_mgdw_per_m3 !Spherical prey
Prey_mgdw = pi*((0.5*PreyWidth*Prey_m)**2.)*Prey_m*zp_mgdw_per_m3 !Cylindric prey (diam=Pw*L)

PreyPer_m3 = Prey_mgdw_per_m3/prey_mgdw					!Prey density

return
end subroutine preydens

subroutine ingestion(act,Eb,BeamAtt,Prey_mgdw_per_m3,Larval_m,PauseDuration,PreyMassEnc)
implicit none

!Parameters
real*8, parameter:: PreyRelLength=0.05				!Prey length relative to fish
real*8, parameter:: PreyWidth = 0.5					!Prey diam:length ratio (for mass)
real*8, parameter:: PreyContrast=0.3
real*8, parameter:: VisRtoSDPreyinBL=1.			!=1 BL for a 0.015 mm2 prey (Fiksen & MacKenzie 2002)
real*8, parameter:: Ke = 1.
real*8, parameter:: omega = 0.
real*8, parameter:: PreySwimVel = 0.
real*8,  parameter :: Prey_handling_s = 1.			!Prey handling time (s)
real*8,  parameter :: rho = 0.3						!Time to move between search volumes
real*8,  parameter :: Pc = 0.2						!Prey capture probability (Fiksen & MacKenzie 2002: 0.2)

real*8, parameter:: VisFieldShape = 0.5
real*8, parameter:: zp_mgdw_per_m3 = 1025*0.25*1.d6 !Mauchline..

!Variables
integer	IER
real*8 Eb,BeamAtt,Prey_mgdw_per_m3,Larval_m,PreyMassEnc
real*8 Prey_m,Prey_mgdw,PreyPer_m3,PreyImageArea,Swim_m_per_s
real*8 VisualRange,EyeSens,pi
real*8 act,PauseFrequency,PauseDuration

!Functions
real*8 Weight_Length
real*8 Length_Weight

pi=acos(-1.)

!Prey properties
Prey_m = Larval_m*PreyRelLength	!Prey length (prosome?)
!Prey_mgdw = (4/3.)*pi*((0.5*Prey_m)**3.)*zp_mgdw_per_m3 !Spherical prey
Prey_mgdw = pi*((0.5*PreyWidth*Prey_m)**2.)*Prey_m*zp_mgdw_per_m3 !Cylindric prey (diam=Pw*L)

PreyPer_m3 = Prey_mgdw_per_m3/prey_mgdw					!Prey density
!PreyImageArea = pi*(Prey_m*0.5)**2						!Spherical prey
PreyImageArea = pi*(Prey_m*PreyWidth)*Prey_m*0.75		!Elongate prey 

!Assume larvae detect 0.1x0.2 mm prey (VisRtoSDPreyinBL x Larval_m) away in clear water 
!(from Fiksen&MacKenzie 2002):
EyeSens = ((Larval_m*VisRtoSDPreyinBL)**2.)/(PreyContrast*(0.1*0.2*0.75*1.d-6))

Swim_m_per_s = Larval_m

!Finding the visual range in m (from AU97)
VisualRange = dsqrt(EyeSens*PreyContrast*PreyImageArea*(Eb/(Ke+Eb))) !Approximation
if(VisualRange>0.05)then	!Exact visual range (above ca 5 cm)
 call getr(VisualRange,BeamAtt,PreyContrast,PreyImageArea,EyeSens,Ke,Eb,IER)
endif

!Cruising larvae:
PreyMassEnc = act*VisFieldShape*pi*(VisualRange**2)*Swim_m_per_s

!Holling disk foraging process
PreyMassEnc = (Pc*PreyMassEnc*PreyPer_m3*Prey_mgdw)/(1. + PreyMassEnc*PreyPer_m3*Prey_handling_s)

if(act.lt.1.d-10)PreyMassEnc = 0. !No feeding if not moving

return
end subroutine ingestion
!******

subroutine FishPred(Eb,BeamAtt,FishDens,Larval_m,FishEnc)
implicit none

!Parameters
real*8, parameter:: LarvalShape = 0.2		!Larval width:length ratio
real*8, parameter:: PreyContrast=0.3
real*8, parameter:: EyeSens = 5.d4
real*8, parameter:: Ke = 1.
real*8, parameter:: VisFieldShape = 0.5
real*8, parameter:: FishSwimVel = 0.10		!Fish cruising velocity (m/s)


!Variables
integer	IER
real*8 Eb,BeamAtt,LarvalWidth,Larval_m,FishDens
real*8 Prey_m,Prey_mgdw,PreyPer_m3,PreyImageArea
real*8 VisualRange,pi,FishEnc

pi=acos(-1.)
LarvalWidth = LarvalShape*Larval_m
PreyImageArea = LarvalWidth*Larval_m*0.75		!Larval image area 

!Finding the visual range in m (from AU97)
VisualRange = dsqrt(EyeSens*PreyContrast*PreyImageArea*(Eb/(Ke+Eb))) !Approximation
if(VisualRange>0.05)then	!Exact visual range (above ca 5 cm)
 call getr(VisualRange,BeamAtt,PreyContrast,PreyImageArea,EyeSens,Ke,Eb,IER)
endif

!Calculate lethal encounter rate with fish:
FishEnc = VisFieldShape*pi*(VisualRange**2)*FishSwimVel*FishDens

return
end subroutine FishPred

subroutine LarvalSDPForward(zOpt,actOpt,V,Eb,PME,FE,Lmin,Lmax,Smax,zmax,DayMax,MaxAct,Horizon,dt,dz,Gutsize,AssimilationEfficiency, &
& BeamAtt,PreyDensDM,PD,Ambush,aPred,bPred,S_mortality,FishDens,Chl_a,Temp,MultSim,SwimCost,SwimLim,ConvDay,zDiff,actDiff,PredationCostOfActivity)
implicit none
!Model resolution
integer  L,Lmin,Lmax,Smax,zmax,Horizon,MaxAct,Activity,S,z,newz,dt

!Variables
real*8,  dimension(Lmin:Lmax,0:Smax,0:zmax,0:Horizon) :: zOpt, ActOpt, V	!Optimal depth decision and reproductive value matrix
real*8,	 dimension(0:zmax,0:Horizon) :: Eb, PME, FE			!Light, prey- and pred enc matrixes
real*8,	 dimension(0:zmax) :: BeamAtt, Chl_a, PreyDensDM, Temp
integer  Day,t,i,j,ij,dz, MultSim, DayMax, Dep, MeanDepth					!Counters for different loops
integer  intS, intL, ConvDay, kk, jjj
real*8   L_m												!length in meters
real*8   delta, dS, dL, NewS, newL, intZ, Fitness, act, PauseDuration
real*8   PreyMassEnc,Lys, PreyPer_m3, PD, PredationCostOfActivity
real*8	 SwimCost,SwimLim,swimDist
real*8   Weight,Growth,age
real*8   Stomach,GutSize,AssimilationEfficiency
real*8   Psurv,aPred,bPred,S_mortality,DeathRate,Ambush,InvPred
real*8	 DayPred,DayInvPred,DayFishPred,DayInvpredSize,DayInvPredAct,DayInvPredMigr
real*8	 FishEnc, FishDens
real*8	 Time10mm,Time15mm,L30,zDiff,actDiff,Ps30,Ps10mm,Ps15mm
real*8	 TotPred,TotInvPred,TotInvpredSize,TotInvPredAct,TotInvPredMigr,TotFishPred

character(12) filename

!Functions
real*8   SGR					!Specific growth rate from Folkvord model
real*8   R						!Respiration
real*8   Length_Weight			!Length as a function of weight
real*8	 Weight_Length			!Weight as a function of length
real*8	 m_starv				!Mortality rate when S = 0
real*8	 m_OtherPred			!Predators other than fish

if(j==0)then
 open(21, file = 'MultSim.txt')
 open(20, file = 'Environment.txt')
  do z = 0,Zmax,dz
   call FishPred(Eb(z,12),BeamAtt(z),FishDens,dble(0.005),FishEnc)	!Expected encounters with predators
   Write(20,'(100E12.3)') dble(z),Temp(z),Chl_a(z),PreyDensDM(z),Eb(z,12),Eb(z,6),Eb(z,0),FishEnc*dt
  enddo
 j=1
endif

!Filename
write(filename, '(a8,f4.1)') 'PreyDens', PD
!Open output files
Open(15, file = filename//'Forward.txt')
Open(16, file = filename//'DayPred.txt')


L_m = dble(Lmin)*0.001					!Initial length
z = 10							!Habitat
Weight = Weight_Length(L_m)		!Initial body mass
Psurv = 1.
Stomach = 0.5*GutSize*Weight_Length(L_m)
kk=0
jjj=0
Day = -1
t = 0			!Starting point in day during forward simulation

do while (Day<99)
Day = Day + 1

do while (t<Horizon)
t = t + 1
if(L_m>=0.019) exit !Stop at large size
if(Weight<0.001)exit !Stop if no growth
age = dble(Day)+dble(t)/dble(Horizon)

!Interpolate states to find optimal depth:
NewS = (Stomach*dble(Smax))/(GutSize*Weight_Length(L_m))			 !Initial stomach content in mg dWeight
intS = max(0,min(idint(NewS),Smax-1))
dS = max(0.,min(NewS - dble(intS),1.))
newL = 1.d3*L_m										!New length in state intervals (mm?)
intL = max(Lmin,min(idint(newL),Lmax-1))
dL = max(0.,min(newL - dble(intL),1.))

!Optimal depth
intZ =  (1.-dL)*(1.-dS)*  zOpt(intL  ,intS  ,z,t) + &   !Length down, stomach down
	         (1.-dL)*    dS *  zOpt(intL  ,intS+1,z,t) + &   !Length down, stomach up
		            dL *(1.-dS)*  zOpt(intL+1,intS  ,z,t) + &   !Length up, stomach down
              dL *    dS *  zOpt(intL+1,intS+1,z,t)	      !Length up, stomach up

!Optimal activity
act =   (1.-dL)*(1.-dS)*ActOpt(intL  ,intS  ,z,t) + &   !Length down, stomach down
	         (1.-dL)*    dS *ActOpt(intL  ,intS+1,z,t) + &   !Length down, stomach up
		            dL *(1.-dS)*ActOpt(intL+1,intS  ,z,t) + &	  !Length up, stomach down
              dL *    dS *ActOpt(intL+1,intS+1,z,t)	     	!Length up, stomach up


Fitness = (1.-dL)*(1.-dS)*     V(intL  ,intS  ,z,t) + &   !Length down, stomach down
          (1.-dL)*    dS *     V(intL  ,intS+1,z,t) + &   !Length down, stomach up
              dL *(1.-dS)*     V(intL+1,intS  ,z,t) + &   !Length up, stomach down
              dL *    dS *     V(intL+1,intS+1,z,t)       !Length up, stomach up 

newz = dz*idnint(intZ/dble(dz))					
swimDist = z - newz
swimDist = abs(swimDist)/(SwimLim*L_m*dble(dt)) !Swimming relative to max..
act = min(1.,max(act,swimDist))					!Including swimming in activity

MeanDepth = idnint(dble(z+newz)/2.)  !Mean depth for environmental exposure

z = newz

!Mass of prey ingested over the timestep
call ingestion(min(1.,act),Eb(MeanDepth,t),BeamAtt(MeanDepth),PreyDensDM(MeanDepth),L_m,PauseDuration,PreyMassEnc)

delta = dt*SGR(weight,Temp(MeanDepth))*Weight/AssimilationEfficiency       !Stomach content required...
delta = delta + dt*R(Temp(MeanDepth),Weight)/AssimilationEfficiency		!...to sustain max growth

	if(delta .le. Stomach)then	!Growth if stomach can supply enough mass
		Growth = Weight*(exp(dt*SGR(Weight,Temp(MeanDepth)))-1.)
		Stomach = Stomach-delta
	else							!New weight if stomach is limiting
		Growth= Stomach*AssimilationEfficiency - dt*R(Temp(MeanDepth),Weight)  
		Stomach = 0.
	endif

	!Stomach = InitialS-delta     !Stomach evacuation
	Stomach = min(Stomach + PreyMassEnc*dt,GutSize*Weight)  !Stomach filling, constrained
	Stomach = max(0.,Stomach)

call FishPred(Eb(MeanDepth,t),BeamAtt(MeanDepth),FishDens,L_m,FishEnc)	!Expected encounters with predators

!Local predation rate from fish, other predators and starvation
InvPred = Ambush*(PredationCostOfActivity*act)*m_OtherPred(Weight,aPred,bPred) + (1.-Ambush)*m_OtherPred(Weight,aPred,bPred)
!InvPred = InvPred + Ambush*swimDist*m_OtherPred(Weight,aPred,bPred)	!Invertebrate predation during migration
DeathRate = InvPred + FishEnc + m_starv(Stomach,S_mortality)
Psurv = Psurv*exp(-DeathRate*dt)						!Risk of dying in dt

!DayInvPredMigr = DayInvPredMigr + Ambush*swimDist*m_OtherPred(Weight,aPred,bPred)*dt
DayInvPredAct = DayInvPredAct + Ambush*PredationCostOfActivity*act*m_OtherPred(Weight,aPred,bPred)*dt
DayInvpredSize = DayInvpredSize + (1.-Ambush)*m_OtherPred(Weight,aPred,bPred)*dt
DayFishPred = DayFishPred + FishEnc*dt

Write(15,'(100E12.4)')age,Psurv,FishEnc*dt,Weight,L_m,&
& Stomach/(GutSize*Weight_Length(L_m)),dble(MeanDepth), &
& Fitness,m_starv(Stomach,S_mortality),InvPred*dt,PreyMassEnc*dt/(GutSize*Weight_Length(L_m)), &
& Eb(MeanDepth,t),growth/Weight,act,Eb(0,t),(FishEnc+InvPred)*dt

!New states
Weight = Weight + Growth - act*dt*R(dble(7),Weight_Length(L_m))		 !Body mass may go up and down
!Weight = Weight - swimDist*SwimCost*dt*R(dble(7),Weight_Length(L_m)) !Energetic migration cost

L_m = max(dble(Lmin)*0.001,Length_Weight(Weight))			       		!Increase in length

if(L_m>0.007.and.kk==0)then
 Time10mm = age
 Ps10mm = Psurv
 kk=1
endif
if(age==30.)then
 L30 = L_m
 Ps30 = Psurv
endif

if(L_m>0.015.and.jjj==0)then
 time15mm = age
 Ps15mm = Psurv
 jjj=1
endif
 
enddo !Diel cycle
DayInvPred = DayInvpredSize + DayInvPredAct + DayInvPredMigr
DayPred = DayFishPred + DayInvpredSize + DayInvPredAct + DayInvPredMigr

Write(16,'(15E12.4)')dble(Day),L_m,DayPred,DayFishPred,DayInvPred,DayInvpredSize,DayInvPredAct, &
  & DayInvPredMigr,exp(-DayPred),DayInvpredSize+DayInvPredAct,DayInvpredSize+DayInvPredAct+DayInvPredMigr

TotPred = TotPred+DayPred								!Get total mortality from different sources
TotInvPred = TotInvPred + DayInvPred
TotFishPred = TotFishPred + DayFishPred
TotInvpredSize = TotInvpredSize + DayInvpredSize
TotInvPredAct = TotInvPredAct + DayInvPredAct
TotInvPredMigr = TotInvPredMigr + DayInvPredMigr
DayPred=0.
DayInvPred=0.
DayInvpredSize=0.
DayInvPredAct=0.
DayInvPredMigr=0.
DayFishPred = 0.
t = 0
enddo !Days
close(15)

!Print results:
write(21,'(100E18.8)')PD,FishDens,Time15mm,Time10mm,L30,Ps30,L_m,Weight,Psurv,V(Lmin,0,0,0),dble(ConvDay),zDiff,actDiff, &
 & TotPred,TotInvPred,TotInvpredSize,TotInvPredAct,TotInvPredMigr,TotFishPred,Ps10mm,Ps15mm

TotPred = 0.								!Reset
TotInvPred = 0.
TotFishPred = 0.
TotInvpredSize = 0.
TotInvPredAct = 0.
TotInvPredMigr = 0.

Time15mm = 0.
Time10mm = 0.
Ps15mm = 0.
Ps10mm = 0.

open(31, file = 'Behaviours.txt')
do t = 1,Horizon
 do Dep = 0,Zmax,dz
   write(31,'(F10.4,2I5,100E15.4)')PD,t,dep,(ActOpt(L,0,Dep,t), Dep - zOpt(L,0,Dep,t), L=5,15)
 enddo
enddo
close(31) 
return
end subroutine

      SUBROUTINE GETR(r,c,C0,Ap,Vc,Ke,Eb,IER)
!     -------------------------------------------------------------

!     Obtain visual range by solving the non-linear equation
!     by means of Newton-Raphson iteration and derivation in
!     subroutine DERIV. Initial value is calculated in EASYR.
!	  The calculation is based on the model described in Aksnes &
!	  Utne (1997) Sarsia 83:137-147.
!
!	  Programmed and tested 29.01.01 Dag L Aksnes
!			
	  implicit none
      REAL*8     r,c,C0,Vc,Ap,Ke,Eb
      REAL*8     EPS,RST,TOL,TOLF,F1,FDER,DX
      INTEGER    I,IEND,AS,IER

!     Input parameters
!		RST		: start value of r calculated by EASYR
!		c		: beam attenuation koefficient (m-1)
!		C0		: prey inherent contrast
!		Ap		: prey area (m^2)
!		Vc		: parameter characterising visual capacity (d.l.) 
!				  this parameter is denoted E' in Aksnes & Utne
!	   	Ke	    : saturation parameter (uE m-2 s-1) 
!	    Eb	    : background irradiance at depth DEPTH

!    Output parameters
!        F1     : function value of equation in DERIV
!        FDER   : the derivative of the function
!        r      : the predator's visual range (when F1=0)
!	   IER	  : = 1, No convergence after IEND steps. Error return.
!				= 2, Return in case of zero divisor.
!				= 3, r out of allowed range (negative)
!				= 0, valid r returned




!     Initial guess of visual range (RST)
      CALL EASYR(RST,C0,Ap,Vc,Ke,Eb)


!     Upper boundary of allowed error of visual range
      EPS = .0000001
!     Maximum number of iteration steps
      IEND = 100

!     Prepare iteration
      r = RST
      TOL = r
  	  CALL DERIV(r,F1,FDER,c,C0,Ap,Vc,Ke,Eb)
      TOLF = 100. * EPS

!     Start iteration loop
      DO 6 I = 1, IEND
         IF (F1) 1, 7, 1

!        Equation is not satisfied by r
1        IF (FDER) 2, 8, 2

!        Iteration is possible
2        DX = F1/FDER
         r = r-DX

!        Test on allowed range
         IF (r .LT. 0.) GOTO 9

         TOL = r
  	     CALL DERIV(r,F1,FDER,c,C0,Ap,Vc,Ke,Eb)

!        Test on satisfactory accuracy
         TOL = EPS
         AS = ABS(r)
         IF (AS-1.) 4, 4, 3
3        TOL = TOL*AS
4        IF (ABS(DX)-TOL) 5, 5, 6
5        IF (ABS(F1)-TOLF) 7, 7, 6
6     CONTINUE

!     No convergence after IEND steps. Error return.
      IER = 1
7     RETURN
!     Return in case of zero divisor
8     IER = 2
      RETURN
!     r out of allowed range (negative)
9     IER = 3
      RETURN

      END
 
!     -------------------------------------------------------------------
      SUBROUTINE EASYR(r,C0,Ap,Vc,Ke,Eb)
!     -------------------------------------------------------------------
!     Obtain a first estimate of visual range by using a simplified
!     expression of visual range
	  implicit none
      REAL*8       r,C0,Ap,Vc,Ke,Eb
      REAL*8       R2
!
!     Se calling routine for explanation of parameters
!
	  R2= ABS(C0)*Ap*Vc*(Eb/(Ke+Eb))
      r = DSQRT(R2)
      RETURN
      END


!     -------------------------------------------------------
      SUBROUTINE DERIV(r,F1,FDER,c,C0,Ap,Vc,Ke,Eb)
!     -------------------------------------------------------
!
!     Derivation of equation for visual range of a predator
	  implicit none
      REAL*8    r,c,C0,Ap,Vc,Ke,Eb
	  REAL*8	FR1,FR2,F1,FDER
     	
!     Input parameters
!        Se explanation in calling routine
!
!    Output parameters
!        F1     : function value of equation in DERIV
!        FDER   : the derivative of the function
!        r      : the predator's visual range (when F1=0)
!
!	   The function and the derivative is calculated on the basis of the 
!	   log-transformed expression


	   FR2=DLOG(ABS(C0)*Ap*Vc)
	   FR1=DLOG(((Ke+Eb)/Eb)*r*r*EXP(c*r))
	   F1 = FR1-FR2
       FDER = c + 2./r 

      RETURN
      END



!     -------------------------------------------------
      SUBROUTINE SURLIG(H,MAXLIG,D,B,SLIG)
!     -------------------------------------------------
!	  Surface irradiance after Skartveit & Olseth 1988
	  implicit none
      REAL*8     B,D,DELTA,HEIGHT,H,P,V,MAXLIG
      REAL*8     SLIG,H12,TWLIGHT

!     B:Degrees north
!     DELTA: sun declination
!     D: day of the year
!     H: hour of day
!     HEIGHT: sin(sunheight)
!     IRR. irradiance above sea surface uEm-2s-1
!     P: Pi
!     R: factor for distance variations between sun-earth 
!     SLIG:surface light
!     V: sunheight in degrees
!     TWLIGHT: light at 0-degree sun
!	  MAXLIG: level of irradiance at midday

      P = 3.1415927
      TWLIGHT = 5.76
		
        DELTA = .3979*DSIND(.9856*(D-80)+		&
            1.9171*(DSIND(.9856*D)-.98112))
		H12 = DELTA*DSIND(B*1.)-					&
              DSQRT(1.-DELTA**2)*DCOSD(B*1.)*COSD(15.*12)

        HEIGHT = DELTA*DSIND(B*1.)-					&
                 DSQRT(1.-DELTA**2)*DCOSD(B*1.)*COSD(15.*H)
		
        V = DASIND(HEIGHT)
		 
          IF (V .GE. 0.) THEN                  
            SLIG = MAXLIG*(HEIGHT/H12) + TWLIGHT
          ELSE IF (V .GE. -6.) THEN
            SLIG = ((TWLIGHT - .048)/6.)*(6.+V)+.048 
          ELSE IF (V .GE. -12.) THEN
            SLIG = ((.048 - 1.15E-4)/6.)*(12.+V)+1.15E-4
          ELSE IF (V .GE. -18) THEN
            SLIG = (((1.15E-4)-1.15E-5)/6.)*(18.+V)+1.15E-5
          ELSE 
            SLIG = 1.15E-5
          ENDIF

        !WRITE(*,800)D,H,HEIGHT/H12,HEIGHT

800   FORMAT(I4,3X,F5.2,3X,2(F15.10,2X))


      RETURN
      END

