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

!****************************************************************************
!
!  PROGRAM: larval_cod_SDP
!
!  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 = 50, 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 (20) 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 = 6., 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 = 50.	!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 = 1,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
 end function m_OtherPred

 



