      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