
 subroutine SetEnvironment(MLDep,MLChla,MLTemp,DLTemp,DeltaTemp,Horizon,zmax,JulianDay,&
 & dz,PD,PDmean_z,PDsd,PreyDensDM,Chl_a,BeamAtt,Eb,Temp,ssurf,CloudCover,LatR,LatD,DaysInYear, FixedLight, MAXLIG)
 implicit none

 integer Horizon,t,zmax,z,dz,L,JulianDay,DaysInYear, FixedLight
 real*8  PD,BeamAtt(0:zmax),MAXLIG, 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,DLTemp,DeltaTemp,PreyPer_m3
 real*8	PDsd,PDmean_z,normal, LatD
 real    radmax, radfl0, LatR, CloudCover, ssurf			!variables needed for qsw0
 
 open(32,file= 'light.txt')
open(28, file = 'PreyDens_z.txt')
 open(30, file = 'DepthEnv.txt')


 !Set Chl_a , Temp profile
 do z = 0,zmax
  ! make Chl_a a function of depth
  if(dble(z)<=MLDep*0.5) then !Half way to MLDep
   Chl_a(z) = MLChla
   PreyDensDM(z) = PD
   !Temp(z) = MLTemp  !- 0.05*dble(z)
  else
   Chl_a(z)=MLChla*MLDep*normal(MLDep*0.4,MLDep/2,z) ! Sloping function to approach zero
   !Chl_a(z) = 0.5
   PreyDensDM(z) = 0 !PD
   !Temp(z) = MLTemp  !- 0.05*dble(z) !-MLDep)
  endif

  ! make temperature a function of depth:
 if(dble(z)<=MLDep*0.5) then ! Half way to MLDep
   Temp(z)=DLTemp+DeltaTemp
 else
   Temp(z)=DLTemp+(DeltaTemp*MLDep*normal(MLDep*0.4,MLDep/2,z)) ! Sloping function (StdDev,mean,z) to approach zero
  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
  if (FixedLight .eq. 0) then ! use subroutine qsw10 from hycom
   call qsw0(radfl0,radmax,ssurf,CloudCover,LatR,float(JulianDay),t,DaysInYear)
   Eb(0,t)=ssurf*4.6		! Converting from w/m2 to uE/m2/s
   !Eb(0,t)=max(ssurf,0.0001)*4.6	
  else ! use fixed function
   call SURLIG(dble(t),MAXLIG,dble(JulianDay),LatD,SLIG)                      
   Eb(0,t) = SLIG !max(SLIG,dble(1E-4)) !Min surface light at night
   
  endif
  do z = 1,zmax
   Eb(z,t) = Eb(z-1,t)*exp((-DiffAtt(z))) !
   
  enddo !depth (z)
 enddo ! hour (t)
 do t = 0,Horizon
  do  z = 0,Zmax,dz
  Write(32,'(100E12.3)') dble(z),dble(t), log10(Eb(z,t))
  enddo
  enddo
 close(32)
 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 

 REAL*8 FUNCTION normal_dble(sd,mean,val)

 IMPLICIT NONE 	
 REAL*8, INTENT(in) :: sd               !standard deviation of normal distribution
 REAL*8, INTENT(in) :: mean		!mean of normal distribution
 real*8, INTENT(in) :: val			!depth 

 REAL*8 expon,pi
 pi = acos(-1.)

 !Normal distribution = 1/(sd*sqrt(2*pi)) * exp(-(x-mean)**2/2*sd**2)
 expon = -(((dble(val) - mean))**2)/(2.*(sd**2))
 normal_dble = exp(expon)/(sqrt(2.*pi)*sd)

 END FUNCTION normal_dble 

 

 !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