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