
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
!PreyDensDM(z)= PD*normal(PDsd,PDmean_z,z)/normal(PDsd,PDmean_z,nint(PDmean_z))
!Temp(z) = MLTemp/(1. + exp(-(MLDep - dble(z))/MLTempSlope))

 !DiffAtt(z) = k0 + 0.054*Chl_a(z)**0.667 + 0.0088*Chl_a(z) !Riley 1956
 !BeamAtt(z) = DiffAtt(z)*BeamDiffCoefficient
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
!Diffatt(z) = (0.064 + 0.06*0.357*(Chl_a(z)**0.65))*(1.+0.2*exp(-0.014*(550-440))) !Oceanic, type 1, Morel 1991
BeamAtt(z) = 0.0657 + (0.39*(Chl_a(z)**0.57))*(1.563 - 0.001149*550) ! equations 3.47 and 3.48

enddo
!print*,sum(PreyDensDM(:))

!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 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
!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))

!or alternatively, the eye-sensitivity from Fiksen & al. 1998:
!Eyesens = 5.04/(1. + exp(-((Larval_m*1d3 - 5.33)/0.64)))
!Eyesens = 10**Eyesens

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 = VisFieldShape*pi*(VisualRange**2)*Swim_m_per_s*PreyPer_m3*Prey_mgdw

!Behavioural option: Jump frequency or pause duration
PauseDuration = MaxPause - act*(MaxPause-MinPause)
PauseFrequency = 1./(PauseDuration + 0.0)

!Saltatory search:
PreyMassEnc = 0.667*pi*(VisualRange**3.)*PauseFrequency + & 
&  pi*(VisualRange**2)*sqrt(PreySwimVel**2 + 2*omega**2)*PauseFrequency*PauseDuration
PreyMassEnc = PreyMassEnc*PreyPer_m3*Prey_mgdw

!print('(10E9.1)'),Larval_m,PreyPer_m3,Prey_m,Prey_mgdw,Eyesens,VisualRange/Larval_m,&
!& PreyMassEnc*3600./Weight_Length(Larval_m),Weight_Length(Larval_m)

return
end subroutine preydens
