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,DLTemp,DeltaTemp,MultSim,SwimCost,SwimLim,ConvDay,zDiff, &
 actDiff,PredationCostOfActivity,Tempname, PDname,x_sd,xnum,PreyVarMagnitude,IndNum,MaxForwardDays,LengthWriteInterval,MaxForwardLength,WriteForwardInds,MaxBin, &
 TimeVarMinLimit,TimeVarMaxLimit,LengthVarMinLimit,LengthVarMaxLimit,MaxTimeTraitNum,MaxLengthTraitNum)
 implicit none
 
 !Model resolution
 integer  L,Lmin,Lmax,Smax,zmax,Horizon,MaxAct,Activity,S,z,newz,dt,xnum,IndNum,TraitNo

 !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,Tempname,PDname, Ind_i,ForwardDays !Counters for different loops
 integer  intS, intL, ConvDay, kk, jjj, intx
 real*8   L_m, DLTemp, DeltaTemp,Eb12,Eb24,EbAv, TempCum !Temp12, Temp24, TempAv,Depth12,Depth24,Eb12,Eb24,EbAv												!length in meters
 real*8   delta, dS, dL, NewS, newL, intZ, Fitness, act, PauseDuration,ActivityAv, ActiveTime
 real*8   PreyMassEnc,PreyMassEnc_allx(0:xnum),Lys, PreyPer_m3, PD, PredationCostOfActivity
 real*8	  SwimCost,SwimLim,swimDist
 real*8   Weight,Growth,age
 real*8   Stomach,GutSize,AssimilationEfficiency,StomachAv,Stomach12,Stomach24,VolumeSearched
 real*8   Psurv,aPred,bPred,S_mortality,DeathRate,Ambush,InvPred, Survival_t,Pweight,TimeInLenIncr,InvPred1
 real*8	  DayPred,DayInvPred,DayFishPred,DayInvpredSize,DayInvPredAct,DayInvPredMigr
 real*8	  FishEnc, FishDens
 real*8	  Time10mm(1:IndNum),Time15mm,L30,zDiff,actDiff,Ps30,Ps10mm,Ps15mm,sumP
 real*8	  TotPred,TotInvPred,TotInvpredSize,TotInvPredAct,TotInvPredMigr,TotFishPred, SizeDep, Fish
 real*8   x_dble, varx, meanx, pi, x_sd, rand1, rand2, dx,yyy
 integer  RightBin,LengthBin
 real*8,  intent(in) :: PreyVarMagnitude


 !For writing
 integer MaxForwardDays,WriteForwardInds,MaxBin,Bin,MaxTimeTraitNum,MaxLengthTraitNum, Timestep
 real*8  LengthWriteInterval,MaxForwardLength, StomachAtBeginningOfTimestep
 real*8  DepthBinned(-2:zmax,0:(Horizon*(MaxForwardDays+1))-1)
 real*8  TimeBinned(1:MaxTimeTraitNum,-2:MaxBin,0:(Horizon*(MaxForwardDays+1))-1)
 real*8  LengthBinned(1:MaxLengthTraitNum,-2:MaxBin,0:idnint(((MaxForwardLength*1000.)-dble(Lmin))/LengthWriteInterval))
 real*8  TimeVarMinLimit(MaxTimeTraitNum),TimeVarMaxLimit(MaxTimeTraitNum),Pweigth
 real*8  LengthVarMinLimit(MaxLengthTraitNum),LengthVarMaxLimit(MaxLengthTraitNum)
 real*8  TraitValTime (1:MaxTimeTraitNum)
 real*8  TraitValLength(1:MaxLengthTraitNum)
 integer MaxLengthBin,OldLengthBin 
 
 character(16) filename2
 character(12) filename
 character(2) number(-2:99)

 !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

  
 do TraitNo=1, 99
  write(number(TraitNo),'(i2.2)') TraitNo
 enddo
 number(-2)='Vr'
 number(-1)='Av'
 number( 0)='Sm'
 MaxLengthBin = idnint((MaxForwardLength*1000.-dble(Lmin))/LengthWriteInterval)
 if(j==0)then
  
  open(21, file = 'MultSim.txt')
  open(20, file = 'Environment.txt')
  open(23,file='FishEnc.txt')
  open(24,file='FishPredsize.txt')
  open(25,file='InvPredsize.txt')
  
  do L=5,15
      call FishPred(Eb(50,19),BeamAtt(50),FishDens,dble(0.001*L),FishEnc)	!Expected encounters with predators
      Write(24,'(100E12.3)') dble(L),FishEnc*dt
  enddo !length
  close(24)

  do t = 0,Horizon
     do z = 0,Zmax,dz
      call FishPred(Eb(z,t),BeamAtt(z),FishDens,dble(0.005),FishEnc)	!Expected encounters with predators
      Write(23,'(100E12.3)') dble(z),dble(t),FishEnc*dt
     enddo !depth
  enddo !hour
  
  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
 write(filename2, '(a8,i2,a4,i2)') 'PreyDens',int(PD+10),'Temp',int(DLTemp+10) !for two loops
 !Open output files
 Open(15, file = 'Forward'//filename2//'.txt')
 write(15,'(100A12)')'ID','Age','Psurv','FishEnc','Weight','Length','Stomach','MeanDepth','Fitness','m_starv','InvPred','PreyMassEnc','EbMeanDep','Growth','act','Eb','Pred'
 Open(16, file = 'DayPred'//filename2//'.txt')

 
 Survival_t = 0.
 TimeBinned = 0.
 LengthBinned = 0.
 DepthBinned = 0.
 
 do Ind_i=1,IndNum,1 !loop over individuals

 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
 Timestep = -1
 TempCum=0
 !LengthBin=1
 ActiveTime = 0

 meanx = dble(xnum)/2.
 varx = x_sd**2
 pi = acos(-1.)
 
 do while (Day<MaxForwardDays)
  Day = Day + 1

  do while (t<Horizon)
   t = t + 1
   Timestep = Timestep + 1
   554 call RANDOM_NUMBER(rand1)
   call RANDOM_NUMBER(rand2)
   x_dble= meanx+dsqrt(-2.*dlog(max(0.00001,rand1)))*dcos(2.*pi*rand2)*dsqrt(varx) !mean +/- sd
   if (((x_dble > dble(xnum)) .or. (x_dble) < 0.)) goto 554  !If this brings newx outside [0,xnum] then try again
   intx = max(0,min(idint(x_dble),xnum-1))
   dx = max(0.,min(x_dble - dble(intx),1.))
   
   if(L_m>=MaxForwardLength) 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
   StomachAtBeginningOfTimestep = Stomach
   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))	
   !newz = 30
   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
   !act = 0.5
   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,xnum,PreyVarMagnitude,PreyMassEnc_allx)

   PreyMassEnc = dx * PreyMassEnc_allx(intx+1) + (1.-dx)*PreyMassEnc_allx(intx)  !interpolated Prey Mass Encountered to continuous environmental value
   
   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
   
   if (Ind_i<WriteForwardInds)then                
   Write(15,'(100E12.4)')dble(Ind_i),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
   endif
   
   !Writing every timestep
   TraitValTime(1) = log10(PSurv) 
   TraitValTime(2) = newL
   TraitValTime(3) = Weight
   TraitValTime(4) = Stomach/(GutSize*Weight_Length(L_m))
   TraitValTime(5) = act
   TraitValTime(6) = Growth
   TraitValTime(7) = dble(MeanDepth)
   TraitValTime(8) = log10(Eb(MeanDepth,t))
   TraitValTime(9) = FishEnc*dt
   TraitValTime(10)= (InvPred)*dt 
   TraitValTime(11)= (FishEnc+InvPred)*dt
   do TraitNo = 1,MaxTimeTraitNum
    RightBin = idint(((TraitValTime(TraitNo)-TimeVarMinLimit(TraitNo))/(TimeVarMaxLimit(TraitNo)-TimeVarMinLimit(TraitNo)))*dble(MaxBin)) + 1
    if ((RightBin .ge. 0) .and. (RightBin .le. MaxBin)) then
     if (TraitNo .eq. 1) then 
      Pweight=1.      !For survival
     else     
      Pweight=PSurv !For all other traits
     endif
     TimeBinned(TraitNo,RightBin,Timestep) = TimeBinned(TraitNo,RightBin,Timestep) + Pweight !Length at beginning of t
    endif
    TimeBinned(TraitNo,0,Timestep) = TimeBinned(TraitNo,0,Timestep) + Pweight  !Total larvae alive for normalization
    TimeBinned(TraitNo,-1,Timestep) = TimeBinned(TraitNo,-1,Timestep) + Pweight*TraitValTime(TraitNo)  !For weighted average
    TimeBinned(TraitNo,-2,Timestep) = TimeBinned(TraitNo,-2,Timestep) + Pweight*(TraitValTime(TraitNo)**2) !For single-pass variance
   enddo
   !Write depth distribution
   DepthBinned( z,Timestep) = DepthBinned(z,Timestep) + PSurv
   DepthBinned(-1,Timestep) = DepthBinned(-1,Timestep) + PSurv*dble(z)  !Weighted mean depth
   DepthBinned(-2,Timestep) = DepthBinned(-2,Timestep) + PSurv*dble(z)**2  !Weight sum of squares for variance calculation
   
   TempCum=Temp(MeanDepth)+Temp(MeanDepth) !cummulative temp exp

   !Writing for every Lengthstep
   LengthBin = idnint((newL-dble(Lmin))/LengthWriteInterval)
   
   if (act .gt. 0) then 
    ActiveTime = ActiveTime + 1
   endif
   TraitValLength( 1) = log10(PSurv)  
   TraitValLength( 2) = Age
   TraitValLength( 3) = Stomach/(GutSize*Weight_Length(L_m))  !StomachAtBeginningOfTimestep   !Av
   TraitValLength( 4) = Stomach/(GutSize*Weight_Length(L_m))  !AtBeginningOfTimestep   !12
   TraitValLength( 5) = Stomach/(GutSize*Weight_Length(L_m))  !AtBeginningOfTimestep   !24
   TraitValLength( 6) = act !ActivityAv
   TraitValLength( 7) = ActiveTime  !Summed over 24h
   TraitValLength( 8) = Growth
   TraitValLength( 9) = dble(MeanDepth)  !12
   TraitValLength(10) = dble(MeanDepth)  !24
   TraitValLength(11) = Temp(MeanDepth)  ! 12
   TraitValLength(12) = Temp(MeanDepth)  ! 24
   TraitValLength(13) = Temp(MeanDepth)  ! Average
   TraitValLength(14) = TempCum !cum
   TraitValLength(15) = log10(Eb(MeanDepth,12))
   TraitValLength(16) = log10(Eb(MeanDepth,24))
   TraitValLength(17) = log10(Eb(MeanDepth,t))
   TraitValLength(18) = Ambush*(PredationCostOfActivity*act)*m_OtherPred(Weight,aPred,bPred)*dt! Ambush
   TraitValLength(19) = (1.-Ambush)*m_OtherPred(Weight,aPred,bPred)*dt !SizeDep
   TraitValLength(20)=  FishEnc*dt !Fish
   TraitValLength(21) = dble(MeanDepth)
   
   
   do TraitNo = 1,MaxLengthTraitNum
    if ((t .ne. 12) .and. ((TraitNo .eq. 4) .or. (TraitNo .eq. 9) .or. (TraitNo .eq. 11) .or. (TraitNo .eq. 15))) cycle
    if ((t .ne. 24) .and. ((TraitNo .eq. 5) .or. (TraitNo .eq. 7) .or. (TraitNo .eq. 10) .or. (TraitNo .eq. 12) .or. (TraitNo .eq. 16))) cycle
    RightBin = idint(((TraitValLength(TraitNo)-LengthVarMinLimit(TraitNo))/(LengthVarMaxLimit(TraitNo)-LengthVarMinLimit(TraitNo)))*dble(MaxBin)) + 1
    if ((RightBin .ge. 0) .and. (RightBin .le. MaxBin)) then
     if (TraitNo .eq. 1) then 
      Pweight=1.      !For survival
     else     
      Pweight=PSurv !For all other traits
     endif
     LengthBinned(TraitNo,RightBin,LengthBin) = LengthBinned(TraitNo,RightBin,LengthBin) + Pweight !Length at beginning of t
     LengthBinned(TraitNo,0,LengthBin) = LengthBinned(TraitNo,0,LengthBin) + Pweight  !Total larvae alive for normalization
     LengthBinned(TraitNo,-1,LengthBin) = LengthBinned(TraitNo,-1,LengthBin) + Pweight*TraitValLength(TraitNo)  !For weighted average
     LengthBinned(TraitNo,-2,LengthBin) = LengthBinned(TraitNo,-2,LengthBin) + Pweight*(TraitValLength(TraitNo)**2) !For single-pass variance
    endif
   enddo
  if (t.eq. 24) ActiveTime = 0
 
 !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(Ind_i) = 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 !Hour
  DayInvPred = DayInvpredSize + DayInvPredAct + DayInvPredMigr
  DayPred = DayFishPred + 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
 
 enddo ! individuals
 
 !Normalize writing data
 do TraitNo = 1, MaxTimeTraitNum
  do Timestep = 1, ((MaxForwardDays+1)*Horizon)-1
   TimeBinned(TraitNo,1:MaxBin,Timestep) = TimeBinned(TraitNo,1:MaxBin,Timestep)/TimeBinned(TraitNo,0,Timestep)
  enddo
  TimeBinned(TraitNo,-1,:)    = TimeBinned(TraitNo,-1,:)/TimeBinned(TraitNo,0,:)  !Mean normalized according to weighting...
  TimeBinned(TraitNo,-2,:)    = TimeBinned(TraitNo,-2,:)/TimeBinned(TraitNo,0,:)  !...and variance.. then 
  TimeBinned(TraitNo,-2,:)    = TimeBinned(TraitNo,-2,:) - (TimeBinned(TraitNo,-1,:)**2)  !...sum of sqaures minus mean squared is variance
 enddo
 do TraitNo = 1, MaxLengthTraitNum
  
  do Bin = 1, MaxBin

  LengthBinned(TraitNo,Bin,:) = LengthBinned(TraitNo,Bin,:)/LengthBinned(TraitNo,0,:)
  enddo
  LengthBinned(TraitNo,-1,:)    = LengthBinned(TraitNo,-1,:)/LengthBinned(TraitNo,0,:)  !Mean normalized according to weighting...
  LengthBinned(TraitNo,-2,:)    = LengthBinned(TraitNo,-2,:)/LengthBinned(TraitNo,0,:)  !... and variance.. then
  LengthBinned(TraitNo,-2,:)    = LengthBinned(TraitNo,-2,:) - (LengthBinned(TraitNo,-1,:)**2)  !...sum of sqaures minus mean squared is variance
 enddo
 do Timestep = 1, ((MaxForwardDays+1)*Horizon)-1
  sumP = sum(DepthBinned(0:zmax,Timestep))
  DepthBinned(0:zmax,Timestep) = DepthBinned(0:zmax,Timestep)/sumP
  DepthBinned(-1,Timestep)     = DepthBinned(-1,Timestep)/sumP
  DepthBinned(-2,Timestep)     = DepthBinned(-2,Timestep)/sumP
 enddo
 DepthBinned(-2,:)     = DepthBinned(-2,:) - (DepthBinned(-1,:)**2)
 
 !writing
 
 open (17, file = 'Timestep'//filename2//'.txt')
 write (17,'(500a22)') 'Age',('LogSurv'//number(Bin),Bin=-2,MaxBin),('Length'//number(Bin),Bin=-2,MaxBin),('Weight'//number(Bin),Bin=-2,MaxBin),('StomachFullness'//number(Bin),Bin=-2,MaxBin),('Activity'//number(Bin),Bin=-2,MaxBin),('Growth'//number(Bin),Bin=-2,MaxBin),('MeanDepth'//number(Bin),Bin=-2,MaxBin),('logEb(ambient)'//number(Bin),Bin=-2,MaxBin),('FishEnc'//number(Bin),Bin=-2,MaxBin),('InvPred'//number(Bin),Bin=-2,MaxBin),('TotPred'//number(Bin),Bin=-2,MaxBin)
 do Timestep = 1, ((MaxForwardDays+1)*Horizon)-1
 write (17,'(500E13.5)') real(Timestep)/real(Horizon), (TimeBinned(TraitNo,:,Timestep),TraitNo = 1,MaxTimeTraitNum)
 enddo
 close (17)
 
 Open(18, file = 'Lengthstep'//filename2//'.txt')
 write (18,'(1000a22)') 'Length',('LogSurv'//number(Bin),Bin=-2,MaxBin),('Age'//number(Bin),Bin=-2,MaxBin),('StomachAv'//number(Bin),Bin=-2,MaxBin),('Stomach12'//number(Bin),Bin=-2,MaxBin),('Stomach24'//number(Bin),Bin=-2,MaxBin),('ActivityAv'//number(Bin),Bin=-2,MaxBin),('ActiveTime'//number(Bin),Bin=-2,MaxBin),('Growth'//number(Bin),Bin=-2,MaxBin),('Depth12'//number(Bin),Bin=-2,MaxBin),('Depth24'//number(Bin),Bin=-2,MaxBin),('Temp12'//number(Bin),Bin=-2,MaxBin),('Temp24'//number(Bin),Bin=-2,MaxBin),('TempAv'//number(Bin),Bin=-2,MaxBin),('TempCum'//number(Bin),Bin=-2,MaxBin),('LogLight12'//number(Bin),Bin=-2,MaxBin),('LogLight24'//number(Bin),Bin=-2,MaxBin),('LogLightAv'//number(Bin),Bin=-2,MaxBin),('Ambush'//number(Bin),Bin=-2,MaxBin),('Sizedep'//number(Bin),Bin=-2,MaxBin),('Fish'//number(Bin),Bin=-2,MaxBin),('MeanDepth'//number(Bin),Bin=-2,MaxBin)
 do intL = 0, MaxLengthBin
  write (18,'(1000E13.5)') real(Lmin)+ real(intL)*LengthWriteInterval, (LengthBinned(TraitNo,:,intL),TraitNo = 1,MaxLengthTraitNum)
 enddo
 close (18)
 
 open (19, file = 'Depthdist'//filename2//'.txt')
 write (19,'(105a22)') 'Age',('Depth'//number(Bin),Bin=-2,-1),'Depth00','Depth01',('Depth'//number(Bin),Bin=2,zmax)
  do Timestep = 1, ((MaxForwardDays+1)*Horizon)-1
  write (19,'(105es13.5)') real(Timestep)/real(Horizon), DepthBinned(:,Timestep)
 enddo
close(19) 

 !Print results:
 write(21,'(100E18.8)')PD,DLTemp,DeltaTemp,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