SUBROUTINE GETR(r,c,C0,Ap,Vc,Ke,Eb,IER)
 !     -------------------------------------------------------------

 !     Obtain visual range by solving the non-linear equation
 !     by means of Newton-Raphson iteration and derivation in
 !     subroutine DERIV. Initial value is calculated in EASYR.
 !	  The calculation is based on the model described in Aksnes &
 !	  Utne (1997) Sarsia 83:137-147.
 !
 !	  Programmed and tested 29.01.01 Dag L Aksnes
 !			
 implicit none
 REAL*8     r,c,C0,Vc,Ap,Ke,Eb
 REAL*8     EPS,RST,TOL,TOLF,F1,FDER,DX
 INTEGER    I,IEND,AS,IER

 !     Input parameters
 !		RST		: start value of r calculated by EASYR
 !		c		: beam attenuation koefficient (m-1)
 !		C0		: prey inherent contrast
 !		Ap		: prey area (m^2)
 !		Vc		: parameter characterising visual capacity (d.l.) 
 !				  this parameter is denoted E' in Aksnes & Utne
 !	   	Ke	    : saturation parameter (uE m-2 s-1) 
 !	    Eb	    : background irradiance at depth DEPTH

 !    Output parameters
 !        F1     : function value of equation in DERIV
 !        FDER   : the derivative of the function
 !        r      : the predator's visual range (when F1=0)
 !	   IER	  : = 1, No convergence after IEND steps. Error return.
 !				= 2, Return in case of zero divisor.
 !				= 3, r out of allowed range (negative)
 !				= 0, valid r returned




 !     Initial guess of visual range (RST)
 CALL EASYR(RST,C0,Ap,Vc,Ke,Eb)


 !     Upper boundary of allowed error of visual range
 EPS = .0000001
 !     Maximum number of iteration steps
 IEND = 100

 !     Prepare iteration
 r = RST
 TOL = r
 CALL DERIV(r,F1,FDER,c,C0,Ap,Vc,Ke,Eb)
 TOLF = 100. * EPS

 !     Start iteration loop
 DO 6 I = 1, IEND
  IF (F1) 1, 7, 1

  !        Equation is not satisfied by r
1 IF (FDER) 2, 8, 2

  !        Iteration is possible
2 DX = F1/FDER
  r = r-DX

  !        Test on allowed range
  IF (r .LT. 0.) GOTO 9

  TOL = r
  CALL DERIV(r,F1,FDER,c,C0,Ap,Vc,Ke,Eb)

  !        Test on satisfactory accuracy
  TOL = EPS
  AS = ABS(r)
  IF (AS-1.) 4, 4, 3
3 TOL = TOL*AS
4 IF (ABS(DX)-TOL) 5, 5, 6
5 IF (ABS(F1)-TOLF) 7, 7, 6
6 CONTINUE

 !     No convergence after IEND steps. Error return.
 IER = 1
7 RETURN
 !     Return in case of zero divisor
8 IER = 2
 RETURN
 !     r out of allowed range (negative)
9 IER = 3
 RETURN

 END

 !     -------------------------------------------------------------------
 SUBROUTINE EASYR(r,C0,Ap,Vc,Ke,Eb)
 !     -------------------------------------------------------------------
 !     Obtain a first estimate of visual range by using a simplified
 !     expression of visual range
 implicit none
 REAL*8       r,C0,Ap,Vc,Ke,Eb
 REAL*8       R2
 !
 !     Se calling routine for explanation of parameters
 !
 R2= ABS(C0)*Ap*Vc*(Eb/(Ke+Eb))
 r = DSQRT(R2)
 RETURN
 END


 !     -------------------------------------------------------
 SUBROUTINE DERIV(r,F1,FDER,c,C0,Ap,Vc,Ke,Eb)
 !     -------------------------------------------------------
 !
 !     Derivation of equation for visual range of a predator
 implicit none
 REAL*8    r,c,C0,Ap,Vc,Ke,Eb
 REAL*8	FR1,FR2,F1,FDER

 !     Input parameters
 !        Se explanation in calling routine
 !
 !    Output parameters
 !        F1     : function value of equation in DERIV
 !        FDER   : the derivative of the function
 !        r      : the predator's visual range (when F1=0)
 !
 !	   The function and the derivative is calculated on the basis of the 
 !	   log-transformed expression


 FR2=DLOG(ABS(C0)*Ap*Vc)
 FR1=DLOG(((Ke+Eb)/Eb)*r*r*EXP(c*r))
 F1 = FR1-FR2
 FDER = c + 2./r 

 RETURN
 END
 
 !     -------------------------------------------------
      SUBROUTINE SURLIG(H,MAXLIG,D,B,SLIG)
!     -------------------------------------------------
!	  Surface irradiance after Skartveit & Olseth 1988
	  implicit none
      REAL*8     B,D,DELTA,HEIGHT,H,P,V,MAXLIG
      REAL*8     SLIG,H12,TWLIGHT

!     B:Degrees north
!     DELTA: sun declination
!     D: day of the year
!     H: hour of day
!     HEIGHT: sin(sunheight)
!     IRR. irradiance above sea surface uEm-2s-1
!     P: Pi
!     R: factor for distance variations between sun-earth 
!     SLIG:surface light
!     V: sunheight in degrees
!     TWLIGHT: light at 0-degree sun
!	  MAXLIG: level of irradiance at midday

      P = 3.1415927
      TWLIGHT = 5.76
		
        DELTA = .3979*DSIND(.9856*(D-80)+		&
            1.9171*(DSIND(.9856*D)-.98112))
		H12 = DELTA*DSIND(B*1.)-					&
              DSQRT(1.-DELTA**2)*DCOSD(B*1.)*COSD(15.*12)

        HEIGHT = DELTA*DSIND(B*1.)-					&
                 DSQRT(1.-DELTA**2)*DCOSD(B*1.)*COSD(15.*H)
		
        V = DASIND(HEIGHT)
		 
          IF (V .GE. 0.) THEN                  
            SLIG = MAXLIG*(HEIGHT/H12) + TWLIGHT
          ELSE IF (V .GE. -6.) THEN
            SLIG = ((TWLIGHT - .048)/6.)*(6.+V)+.048 
          ELSE IF (V .GE. -12.) THEN
            SLIG = ((.048 - 1.15E-4)/6.)*(12.+V)+1.15E-4
          ELSE IF (V .GE. -18) THEN
            SLIG = (((1.15E-4)-1.15E-5)/6.)*(18.+V)+1.15E-5
          ELSE 
            SLIG = 1.15E-5
          ENDIF

        !WRITE(*,800)D,H,HEIGHT/H12,HEIGHT

800   FORMAT(I4,3X,F5.2,3X,2(F15.10,2X))


      RETURN
      END
      !end subroutine SURLIG




subroutine qsw0(radfl0,radmax,ssurf,clouds,rlat,time,npart,daysinyear) 
!---------------------------------------------------------------------------------------------

!Anders.
!Modify qsw0 to accept time of day (npart) as an input variable and respond by returning
!an hour value of radiafion(ssurf) used in source file "Environment.f90".
!---------------------------------------------------------------------------------------------

	  implicit none
! radfl0 (out)		- 24 hrs mean solar irrradiance at the marine surface layer --- (unit: w/m^2)
! radmax (out)		- daily maximum solar irrradiance at the marine surface layer --- (unit: w/m^2)
! ssurf (out)		- Anders: adding ssurf as an output variable
! clouds (in)		- cloud input (fraction of cloudy sky)
! rlat (in)			- latitude (radians)
! time (in)			- julian day
! daysinyear (in)	- number of days in the year
! npart (in)		- Anders: adding npart as an input variable 

      real, intent(out) :: radfl0,radmax,ssurf	!Anders: adding ssurf
      real, intent(in)  ::  clouds,rlat
      real, intent(in) :: time
      integer, intent(in)  :: daysinyear, npart	!Anders: adding npart 

      real pi2,deg,rad,eepsil
      integer ifrac !npart									!Anders: removing npart
      real fraci,absh2o,s0,day,dangle,decli,sundv
      real cc,sin2,cos2,cosz,scosz,stot,bioday,biohr,hangle
      real srad,sdir,sdif,altdeg,cfac !ssurf				!Anders: removing ssurf




! --- -------------------------------------------------------------------
! --- compute 24 hrs mean solar irrradiance at the marine surface layer
! --- (unit: w/m^2)
! --- -------------------------------------------------------------------

! --- set various quantities
 
      !print *,'calling qsw0'
      pi2=8.*atan(1.)          !        2 times pi
      deg=360./pi2             !        convert from radians to degrees
      rad=pi2/360.             !        convert from degrees to radians
      eepsil=1.e-9             !        small number

      
	  !ifrac=24                !        split each 12 hrs day into ifrac parts
      ifrac =12				   !		Anders: changing ifrac from 24 to 12
	  
	  fraci=1./real(ifrac)     !        1 over ifrac

      absh2o=0.09              ! ---    absorption of water and ozone
      s0=1365.                 ! w/m^2  solar constant

! --- -------------------------------------------------------------------
! --- compute 24 hrs mean solar radiation at the marine surface layer 
! --- -------------------------------------------------------------------

!KAL  ttime=time+dt/86400.
!KAL? day=aint(time*365./360.)          !accumulated day number (jan1=0,364,..)
      !day=aint(time/float(daysinyear))  !accumulated day number (jan1=0,364,..)
      day=amod(time,float(daysinyear))    !0 < day < 364
      day=floor(day)
      dangle=pi2*day/float(daysinyear)   !day-number-angle, in radians 
      if (day<0. .or. day>daysinyear+1) then
         print *,'qsw0: Error in day for day angle'
         print *,'Day angle is ',day,daysinyear
         stop
      end if

      


! --- compute astronomic quantities -- 
      decli=.006918+.070257*sin(dangle)   -.399912*cos(dangle) &
                  +.000907*sin(2.*dangle)-.006758*cos(2.*dangle) &
                  +.001480*sin(3.*dangle)-.002697*cos(3.*dangle)

      sundv=1.00011+.001280*sin(dangle)   +.034221*cos(dangle) &
                  +.000077*sin(2.*dangle)+.000719*cos(2.*dangle)

! --- compute cloudiness fraction

! KAL cc=clouds(i,j,l0)*w0+clouds(i,j,l1)*w1
! KAL.  +clouds(i,j,l2)*w2+clouds(i,j,l3)*w3
      cc = clouds

! --- compute astronomic quantities

      sin2=sin(rlat)*sin(decli)
      cos2=cos(rlat)*cos(decli)

! --- split each day into ifrac parts, and compute the solar radiance for 
! --- each part. by assuming symmetry of the irradiance about noon, it
! --- is sufficient to compute the irradiance for the first 12 hrs of
! --- the (24 hrs) day (mean for the first 12 hrs equals then the mean
! --- for the last 12 hrs)

      scosz=0.
      stot=0.
      radmax=0.0

!      do npart=1,ifrac									!Anders: removing do loop as npart is now set as an input variable
         bioday=day+(npart-.5)*fraci*.5
         biohr=bioday*86400.							!hour of day in seconds
         biohr=amod(biohr+43200.,86400.)				!hour of day;  biohr=0  at noon
         hangle=pi2*biohr/86400.						!hour angle, in radians
         cosz=amax1(0.,sin2+cos2*cos(hangle))			!cosine of the zenith angle
         scosz=scosz+cosz								 !  ..accumulated..

         srad =s0*sundv*cosz							 !extraterrestrial radiation

!         print *,i,j,npart,srad,cosz,eepsil
         ! obs: .7^100 = 3x10^-16 , an already ridicolously low number ...
         sdir=srad*0.7**(min(100.,1./(cosz+eepsil)))    !direct radiation component
!         sdir=srad*0.7**(1./(cosz+eepsil))				!direct radiation component
!         sdir=srad * exp(-0.356674943938732447/(cosz+eepsil))         

         sdif=((1.-absh2o)*srad-sdir)*.5				!diffusive radiation component

         altdeg=amax1(0.,asin(min(1.0,sin2+cos2)))*deg	!solar noon altitude in degrees


         cfac=(1.-0.62*cc+0.0019*altdeg)				!cloudiness correction 

         ssurf=(sdir+sdif)*cfac
         radmax=max(radmax,ssurf)
         stot=stot+ssurf
	

!      enddo											!Anders: removing do loop

      scosz=scosz*fraci									!24-hrs mean of  cosz
      radfl0=stot*fraci									!24-hrs mean shortw rad in w/m^2

!      cawdir=1.-amax1(0.15,0.05/(scosz+0.15)) 
!      print *,time,day,dangle,decli,sundv,radfl0(63,95)
!      co-albedo over water for dir light
	  print*,npart,ssurf
      
 end subroutine qsw0

 
 