text:
!===============================================================================
! One of Andy Nowacki's Fortran utility modules for dealing with seismic
! anisotropy and other problems.
!
! Andy Nowacki <andy.nowacki@bristol.ac.uk>
!
! See the file LICENCE for licence details.
!===============================================================================
module statistical
!===============================================================================
!  Contains statistical functions
!  Andy Nowacki, University of Bristol
!
!  2011-06-13   + First iteration.  Includes F-distribution, Beta- (B), regularized
!                 Beta- (I) and incomplete Beta-functions
!  2011-06-17   + Added functions from Alan Miller's website to compute percentage
!                 values for the chi2 distribution.
!  2011-08-09   + Added circular statistical functions.

!  Declare public functions
!   public f_dist

!  Declare some internal constants
!  ** size constants
      integer, parameter, private :: i4 = selected_int_kind(9) ; ! long int
      integer, parameter, private :: r4 = selected_real_kind(6,37) ; ! SP
      integer, parameter, private :: r8 = selected_real_kind(15,307) ; ! DP
      integer, parameter, private :: c4 = r4
      integer, parameter, private :: c8 = r8

!  ** precision selector
      integer, parameter, private :: rs = r8
      integer, parameter, private :: cs = c8

!  ** maths constants and other useful things
      real(rs), parameter, private :: pi = 3.141592653589793238462643_rs

!  ** IO units
      integer, parameter, private :: lu_stdout = 5
      integer, parameter, private :: lu_stdin  = 6
      integer, parameter, private :: lu_stderr = 0

contains

!===============================================================================
function fact12(n)
!===============================================================================
!  Gives true integer value for n!, n <= 12
   implicit none
   integer,intent(in) :: n
   integer            :: fact12,i

   if (n < 0) then
      write(lu_stderr,'(a)') 'statistical: fact: error: n must be >= 0.'
      stop
   else if (n == 0) then
      fact12 = 1
      return
   else if (n > 12) then
      write(lu_stderr,'(a)') 'statistical: fact: error: n cannot be larger than 12 for 4-byte integer calculations.'
      stop
   endif
   fact12 = 1
   do i=1,n
      fact12 = i*fact12
   enddo
   return
end function fact12
!-------------------------------------------------------------------------------

!===============================================================================
function fact(n)
!===============================================================================
!  Gives approximate value of n! up to n<170
!  fact(n) ≈ exp{n.ln(n) - n + ln(n(1 + 4n(1 + 2n)))/6 + ln(pi)/2}
!  Approximation given by Ramanujan (1988)
!  Good to within ~1e-5 of integer value up to n=12.
   implicit none
   integer,intent(in) :: n
   real(rs)           :: fact,rn

   if (n == 0) then
      fact = 1._rs
      return
   else if (n < 0) then
      write(lu_stderr,'(a)') 'statistical: fact: error: n must be >= 0.'
      stop
   else if (n > 0 .and. n <= 12) then
      fact = real(fact12(n))
      return
   else if (n > 170) then
      write(lu_stderr,'(a)') 'statistical: fact: error: cannot represent n! in double precision for n > 170.'
      stop
   endif

   rn = real(n)
   fact = exp(n*log(rn) - rn + log(rn*(1._rs + 4._rs*rn*(1._rs + 2._rs*rn)))/6._rs &
                                                         + log(pi)/2._rs)
   return
end function fact
!-------------------------------------------------------------------------------

!===============================================================================
function stat_poisson_pmf(lambda, k)
!===============================================================================
! Return the probability mass function for the Poisson distribution with mean
! lambda at integer point k
   implicit none
   real(rs), intent(in) :: lambda
   integer, intent(in) :: k
   real(rs) :: stat_poisson_pmf
   stat_poisson_pmf = lambda**k*exp(-lambda)/fact(k)
end function stat_poisson_pmf
!-------------------------------------------------------------------------------

!===============================================================================
function beta_func(p,q)
!===============================================================================
!  Returns the Beta function B(p,q) = Gamma(p)*Gamma(q) / Gamma(p+q)
   implicit none
   real(rs),intent(in) :: p,q
   real(rs)            :: beta_func,rp,rq

   rp = real(p)  ;  rq = real(q)
   beta_func = gamma(rp) * gamma(rq) / gamma(rp + rq)
   return
end function beta_func
!-------------------------------------------------------------------------------

!===============================================================================
function incomp_beta_func(z,a,b)
!===============================================================================
!  Returns the incomplete Beta function B(z;a,b) = ∫_0^z(u^(a-1)*(1-u)^(b-1))du
   implicit none
   real(rs),intent(in) :: z,a,b
   real(rs)            :: incomp_beta_func,ra,rb,rn
   real(rs)            :: u,du   ! Dummy variable of integration
   integer             :: n

   ra = real(a)
   rb = real(b)
   if (z <= 0. .or. z > 1.) then
      write(lu_stderr,'(a)') 'statistical: incomp_beta_func: error: z must be in range [0,1].'
      stop
   endif
   du = 0.000001_rs
   u = 0._rs
   incomp_beta_func = 0._rs
   do while (u <= z)
      incomp_beta_func = incomp_beta_func + (u**(a-1._rs))*((1._rs-u)**(b-1._rs))*du
      u = u + du
   enddo
   incomp_beta_func = incomp_beta_func / beta_func(a,b)
   return
end function incomp_beta_func
!-------------------------------------------------------------------------------

!===============================================================================
function reg_beta_func(z,a,b)
!===============================================================================
!  Returns the regularised Beta function B(a,b) = B(z;a,b)/B(a,b)
   implicit none
   real(rs),intent(in) :: z,a,b
   real(rs)            :: reg_beta_func

   reg_beta_func = incomp_beta_func(z,a,b) / beta_func(a,b)
   return
end function reg_beta_func
!-------------------------------------------------------------------------------

!===============================================================================
function pochhammer(x,n)
!===============================================================================
!  Give the Pochhhammer symbol p(x)_n == gamma(x+n)/gamma(x)
   implicit none
   real(rs),intent(in) :: x
   integer,intent(in)  :: n
   real(rs)            :: pochhammer,rn

   rn = real(n)
   pochhammer = gamma(x+rn)/gamma(x)
   return
end function pochhammer
!-------------------------------------------------------------------------------

!===============================================================================
function f_dist(n,m,x)
!===============================================================================
!  Returns the value of the f distribution f(n,m,x)
   implicit none
   integer             :: m,n
   real(rs)            :: x
   real(rs)            :: f_dist
   real(rs)            :: rn,rm,half,one

   if (x <= 0.) then
      write(lu_stderr,'(a)') 'statistical: f_dist: error: x must be > 0.'
      stop
   endif

!  Convert input integers into real values and set value of half
   rn = real(n)
   rm = real(m)
   half = 0.5_rs
   one = 1.0_rs

   f_dist = gamma(half*(rn+rm)) * rn**(half*rn) * rm**(half*rm) * x**(half*rn - one) / &
            (gamma(half*rn) * gamma(half*rm) * (rm+rn*x)**(half*(rn + rm)))
   return
end function f_dist
!-------------------------------------------------------------------------------

!===============================================================================
function f_dist_cum(n,m,x)
!===============================================================================
!  Returns the value of the cumulative F distribution F(n,m,x)
   implicit none
   integer,intent(in)  :: m,n
   real(rs),intent(in) :: x
   real(rs)            :: f_dist_cum,u,du
   real(rs)            :: rm,rn
!  Convert input integers into real values
   rm = real(m)
   rn = real(n)

!   f_dist_cum = reg_beta_func(rn*x/(rm+rn*x), rn*0.5_rs, rm*0.5_rs)
   f_dist_cum = 0._rs
   du = 0.0001_rs
   u = 0._rs + du
   do while (u <= x)
      f_dist_cum = f_dist_cum + f_dist(n,m,u)*du
      u = u + du
   enddo
   return
end function f_dist_cum
!-------------------------------------------------------------------------------

!===============================================================================
!===============================================================================
! The following functions are taken from Alan Miller's page at:
! http://jblevins.org/mirror/amiller/
! All are transcriptions into Fortran of algorithms plublished in the Royal Statistical
! Society's Applied Statitstics journal.

!===============================================================================
function ppchi2(p, v, g) RESULT(fn_val)
!===============================================================================
! N.B. Argument IFAULT has been removed.
! This version by Alan Miller
! amiller @ bigpond.net.au
! Latest revision - 27 October 2000
!  Algorithm AS 91   Appl. Statist. (1975) Vol.24, P.35
!
!  To evaluate the percentage points of the chi-squared
!  probability distribution function.
!
!  p must lie in the range 0.000002 to 0.999998,
!  v must be positive,
!  g must be supplied and should be equal to ln(gamma(v/2.0))
!
!  Incorporates the suggested changes in AS R85 (vol.40(1), pp.233-5, 1991)
!  which should eliminate the need for the limited range for p above,
!  though these limits have not been removed from the routine.
!
!  If IFAULT = 4 is returned, the result is probably as accurate as
!  the machine will allow.
!
!  Auxiliary routines required: PPND = AS 111 (or AS 241) and GAMMAD = AS 239.

   IMPLICIT NONE
   INTEGER, PARAMETER    :: dp = rs

   REAL (dp), INTENT(IN)  :: p
   REAL (dp), INTENT(IN)  :: v
   REAL (dp), INTENT(IN)  :: g
   REAL (dp)              :: fn_val

! AJN: don't need interface as we're in a module.
!   INTERFACE
!     FUNCTION gammad(x, p) RESULT(fn_val)
!       IMPLICIT NONE
!       INTEGER, PARAMETER    :: dp = rs
!       REAL (dp), INTENT(IN) :: x, p
!       REAL (dp)             :: fn_val
!     END FUNCTION gammad
!
!     SUBROUTINE ppnd16 (p, normal_dev, ifault)
!       IMPLICIT NONE
!       INTEGER, PARAMETER      :: dp = rs
!       REAL (dp), INTENT(IN)   :: p
!       INTEGER, INTENT(OUT)    :: ifault
!       REAL (dp), INTENT(OUT)  :: normal_dev
!     END SUBROUTINE ppnd16
!   END INTERFACE

! Local variables

   REAL (dp)  :: a, b, c, p1, p2, q, s1, s2, s3, s4, s5, s6, t, x, xx
   INTEGER    :: i, if1

   INTEGER, PARAMETER    :: maxit = 20
   REAL (dp), PARAMETER  :: aa = 0.6931471806_dp, e = 0.5e-06_dp,         &
                            pmin = 0.000002_dp, pmax = 0.999998_dp,       &
                            zero = 0.0_dp, half = 0.5_dp, one = 1.0_dp,   &
                            two = 2.0_dp, three = 3.0_dp, six = 6.0_dp,   &
                            c1 = 0.01_dp, c2 = 0.222222_dp, c3 = 0.32_dp, &
                            c4 = 0.4_dp, c5 = 1.24_dp, c6 = 2.2_dp,       &
                            c7 = 4.67_dp, c8 = 6.66_dp, c9 = 6.73_dp,     &
                            c10 = 13.32_dp, c11 = 60.0_dp, c12 = 70.0_dp, &
                            c13 = 84.0_dp, c14 = 105.0_dp, c15 = 120.0_dp, &
                            c16 = 127.0_dp, c17 = 140.0_dp, c18 = 175.0_dp, &
                            c19 = 210.0_dp, c20 = 252.0_dp, c21 = 264.0_dp, &
                            c22 = 294.0_dp, c23 = 346.0_dp, c24 = 420.0_dp, &
                            c25 = 462.0_dp, c26 = 606.0_dp, c27 = 672.0_dp, &
                            c28 = 707.0_dp, c29 = 735.0_dp, c30 = 889.0_dp, &
                            c31 = 932.0_dp, c32 = 966.0_dp, c33 = 1141.0_dp, &
                            c34 = 1182.0_dp, c35 = 1278.0_dp, c36 = 1740.0_dp, &
                            c37 = 2520.0_dp, c38 = 5040.0_dp

!  Test arguments and initialise
   fn_val = -one
   IF (p < pmin .OR. p > pmax) THEN
      WRITE(lu_stderr,'(a)') 'statistical: PPCHI2: error: p must be between 0.000002 & 0.999998'
      stop !RETURN
   END IF
   IF (v <= zero) THEN
      WRITE(lu_stderr,'(a)') 'statistical: PPCHI2: error: Number of deg. of freedom <= 0'
      stop !RETURN
   END IF

   xx = half * v
   c = xx - one

!  Starting approximation for small chi-squared
   IF (v < -c5 * LOG(p)) THEN
      fn_val = (p * xx * EXP(g + xx * aa)) ** (one/xx)
      IF (fn_val < e) GO TO 6
      GO TO 4
   END IF

!  Starting approximation for v less than or equal to 0.32
   IF (v > c3) GO TO 3
   fn_val = c4
   a = LOG(one-p)

   2 q = fn_val
   p1 = one + fn_val * (c7+fn_val)
   p2 = fn_val * (c9 + fn_val * (c8 + fn_val))
   t = -half + (c7 + two * fn_val) / p1 - (c9 + fn_val * (c10 + three * fn_val)) / p2
   fn_val = fn_val - (one - EXP(a + g + half * fn_val + c * aa) * p2 / p1) / t
   IF (ABS(q / fn_val - one) > c1) GO TO 2
   GO TO 4

!  Call to algorithm AS 241 - note that p has been tested above.
   3 CALL ppnd16(p, x, if1)

!  Starting approximation using Wilson and Hilferty estimate
   p1 = c2 / v
   fn_val = v * (x * SQRT(p1) + one - p1) ** 3

!  Starting approximation for p tending to 1
   IF (fn_val > c6 * v + six) fn_val = -two * (LOG(one-p) - c * LOG(half * fn_val) + g)

!  Call to algorithm AS 239 and calculation of seven term Taylor series
   4 DO i = 1, maxit
      q = fn_val
      p1 = half * fn_val
      p2 = p - gammad(p1, xx)

      t = p2 * EXP(xx * aa + g + p1 - c * LOG(fn_val))
      b = t / fn_val
      a = half * t - b * c
      s1 = (c19 + a * (c17 + a * (c14 + a * (c13 + a * (c12 + c11 * a))))) / c24
      s2 = (c24 + a * (c29 + a * (c32 + a * (c33 + c35 * a)))) / c37
      s3 = (c19 + a * (c25 + a * (c28 + c31 * a))) / c37
      s4 = (c20 + a * (c27 + c34 * a) + c * (c22 + a * (c30 + c36 * a))) / c38
      s5 = (c13 + c21 * a + c * (c18 + c26 * a)) / c37
      s6 = (c15 + c * (c23 + c16 * c)) / c38
      fn_val = fn_val + t * (one + half * t * s1 - b * c * (s1 - b *   &
               (s2 - b * (s3 - b * (s4 - b * (s5 - b * s6))))))
      IF (ABS(q / fn_val - one) > e) RETURN
   END DO

   WRITE(lu_stderr,'(a)') 'statistical: PPCHI2: error: Max. number of iterations exceeded'
   stop

   6 RETURN
END FUNCTION ppchi2
!-------------------------------------------------------------------------------

!===============================================================================
function gammad(x, p) RESULT(fn_val)
!===============================================================================
!  ALGORITHM AS239  APPL. STATIST. (1988) VOL. 37, NO. 3

!  Computation of the Incomplete Gamma Integral

!  Auxiliary functions required: ALOGAM = logarithm of the gamma
!  function, and ALNORM = algorithm AS66

! ELF90-compatible version by Alan Miller
! Latest revision - 27 October 2000

! N.B. Argument IFAULT has been removed
! AJN: Remove need for alogam by using log_gamma, intrinsic in Fortran 2008 and
!      implemented in most recent compilers.

   IMPLICIT NONE
   INTEGER, PARAMETER    :: dp = rs
   REAL (dp), INTENT(IN) :: x, p
   REAL (dp)             :: fn_val

   ! Local variables
   REAL (dp)             :: pn1, pn2, pn3, pn4, pn5, pn6, arg, c, rn, a, b, an
   REAL (dp), PARAMETER  :: zero = 0.d0, one = 1.d0, two = 2.d0, &
                            oflo = 1.d+37, three = 3.d0, nine = 9.d0, &
                            tol = 1.d-14, xbig = 1.d+8, plimit = 1000.d0, &
                            elimit = -88.d0
! EXTERNAL alogam, alnorm

   fn_val = zero

!  Check that we have valid values for X and P
   IF (p <= zero .OR. x < zero) THEN
     WRITE(lu_stderr,'(a)') 'statistical: gammad(AS239): error: Either p <= 0 or x < 0'
     stop !RETURN
   END IF
   IF (x == zero) RETURN

!  Use a normal approximation if P > PLIMIT
   IF (p > plimit) THEN
     pn1 = three * SQRT(p) * ((x / p) ** (one / three) + one /(nine * p) - one)
     fn_val = alnorm(pn1, .false.)
     RETURN
   END IF

!  If X is extremely large compared to P then set fn_val = 1
   IF (x > xbig) THEN
     fn_val = one
     RETURN
   END IF

   IF (x <= one .OR. x < p) THEN

!  Use Pearson's series expansion.
!  (Note that P is not large enough to force overflow in ALOGAM).
!  No need to test IFAULT on exit since P > 0.
!     arg = p * LOG(x) - x - alogam(p + one, ifault)
     arg = p * LOG(x) - x - log_gamma(p + one)
     c = one
     fn_val = one
     a = p
     40   a = a + one
     c = c * x / a
     fn_val = fn_val + c
     IF (c > tol) GO TO 40
     arg = arg + LOG(fn_val)
     fn_val = zero
     IF (arg >= elimit) fn_val = EXP(arg)

   ELSE

!  Use a continued fraction expansion
!     arg = p * LOG(x) - x - alogam(p, ifault)
     arg = p * LOG(x) - x - log_gamma(p)
     a = one - p
     b = a + x + one
     c = zero
     pn1 = one
     pn2 = x
     pn3 = x + one
     pn4 = x * b
     fn_val = pn3 / pn4
     60   a = a + one
     b = b + two
     c = c + one
     an = a * c
     pn5 = b * pn3 - an * pn1
     pn6 = b * pn4 - an * pn2
     IF (ABS(pn6) > zero) THEN
       rn = pn5 / pn6
       IF (ABS(fn_val - rn) <= MIN(tol, tol * rn)) GO TO 80
       fn_val = rn
     END IF

     pn1 = pn3
     pn2 = pn4
     pn3 = pn5
     pn4 = pn6
     IF (ABS(pn5) >= oflo) THEN

!  Re-scale terms in continued fraction if terms are large
       pn1 = pn1 / oflo
       pn2 = pn2 / oflo
       pn3 = pn3 / oflo
       pn4 = pn4 / oflo
     END IF
     GO TO 60
     80   arg = arg + LOG(fn_val)
     fn_val = one
     IF (arg >= elimit) fn_val = one - EXP(arg)
   END IF

   RETURN
END FUNCTION gammad
!-------------------------------------------------------------------------------

!===============================================================================
subroutine ppnd16 (p, normal_dev, ifault)
!===============================================================================
! ALGORITHM AS241  APPL. STATIST. (1988) VOL. 37, NO. 3

! Produces the normal deviate Z corresponding to a given lower
! tail area of P; Z is accurate to about 1 part in 10**16.

! The hash sums below are the sums of the mantissas of the
! coefficients.   They are included for use in checking
! transcription.

! This ELF90-compatible version by Alan Miller - 20 August 1996
! N.B. The original algorithm is as a function; this is a subroutine

   IMPLICIT NONE

   INTEGER, PARAMETER      :: dp = rs
   REAL (dp), INTENT(IN)   :: p
   INTEGER, INTENT(OUT)    :: ifault
   REAL (dp), INTENT(OUT)  :: normal_dev

! Local variables
   REAL (dp) :: zero = 0.d0, one = 1.d0, half = 0.5d0,  &
                split1 = 0.425d0, split2 = 5.d0, const1 = 0.180625d0, &
                const2 = 1.6d0, q, r

! Coefficients for P close to 0.5
   REAL (dp) :: a0 = 3.3871328727963666080D0, &
                a1 = 1.3314166789178437745D+2, &
                a2 = 1.9715909503065514427D+3, &
                a3 = 1.3731693765509461125D+4, &
                a4 = 4.5921953931549871457D+4, &
                a5 = 6.7265770927008700853D+4, &
                a6 = 3.3430575583588128105D+4, &
                a7 = 2.5090809287301226727D+3, &
                b1 = 4.2313330701600911252D+1, &
                b2 = 6.8718700749205790830D+2, &
                b3 = 5.3941960214247511077D+3, &
                b4 = 2.1213794301586595867D+4, &
                b5 = 3.9307895800092710610D+4, &
                b6 = 2.8729085735721942674D+4, &
                b7 = 5.2264952788528545610D+3
! HASH SUM AB    55.8831928806149014439

! Coefficients for P not close to 0, 0.5 or 1.
   REAL (dp) :: c0 = 1.42343711074968357734D0, &
                c1 = 4.63033784615654529590D0, &
                c2 = 5.76949722146069140550D0, &
                c3 = 3.64784832476320460504D0, &
                c4 = 1.27045825245236838258D0, &
                c5 = 2.41780725177450611770D-1, &
                c6 = 2.27238449892691845833D-2, &
                c7 = 7.74545014278341407640D-4, &
                d1 = 2.05319162663775882187D0, &
                d2 = 1.67638483018380384940D0, &
                d3 = 6.89767334985100004550D-1, &
                d4 = 1.48103976427480074590D-1, &
                d5 = 1.51986665636164571966D-2, &
                d6 = 5.47593808499534494600D-4, &
                d7 = 1.05075007164441684324D-9
! HASH SUM CD    49.33206503301610289036

! Coefficients for P near 0 or 1.
   REAL (dp) :: e0 = 6.65790464350110377720D0, &
                e1 = 5.46378491116411436990D0, &
                e2 = 1.78482653991729133580D0, &
                e3 = 2.96560571828504891230D-1, &
                e4 = 2.65321895265761230930D-2, &
                e5 = 1.24266094738807843860D-3, &
                e6 = 2.71155556874348757815D-5, &
                e7 = 2.01033439929228813265D-7, &
                f1 = 5.99832206555887937690D-1, &
                f2 = 1.36929880922735805310D-1, &
                f3 = 1.48753612908506148525D-2, &
                f4 = 7.86869131145613259100D-4, &
                f5 = 1.84631831751005468180D-5, &
                f6 = 1.42151175831644588870D-7, &
                f7 = 2.04426310338993978564D-15
! HASH SUM EF    47.52583317549289671629

   ifault = 0
   q = p - half
   IF (ABS(q) <= split1) THEN
      r = const1 - q * q
      normal_dev = q * (((((((a7*r + a6)*r + a5)*r + a4)*r + a3)*r + a2)*r + a1)*r + a0) / &
               (((((((b7*r + b6)*r + b5)*r + b4)*r + b3)*r + b2)*r + b1)*r + one)
      RETURN
   ELSE
      IF (q < zero) THEN
         r = p
      ELSE
         r = one - p
      END IF
      IF (r <= zero) THEN
         ifault = 1
         normal_dev = zero
         RETURN
      END IF
      r = SQRT(-LOG(r))
      IF (r <= split2) THEN
         r = r - const2
         normal_dev = (((((((c7*r + c6)*r + c5)*r + c4)*r + c3)*r + c2)*r + c1)*r + c0) / &
                  (((((((d7*r + d6)*r + d5)*r + d4)*r + d3)*r + d2)*r + d1)*r + one)
      ELSE
         r = r - split2
         normal_dev = (((((((e7*r + e6)*r + e5)*r + e4)*r + e3)*r + e2)*r + e1)*r + e0) / &
                  (((((((f7*r + f6)*r + f5)*r + f4)*r + f3)*r + f2)*r + f1)*r + one)
      END IF
      IF (q < zero) normal_dev = - normal_dev
      RETURN
   END IF
END SUBROUTINE ppnd16
!-------------------------------------------------------------------------------

!===============================================================================
function alnorm( x, upper ) RESULT( fn_val )
!===============================================================================
!  Algorithm AS66 Applied Statistics (1973) vol.22, no.3

!  Evaluates the tail area of the standardised normal curve
!  from x to infinity if upper is .true. or
!  from minus infinity to x if upper is .false.

! ELF90-compatible version by Alan Miller
! Latest revision - 29 November 2001

   IMPLICIT NONE
   INTEGER, PARAMETER     ::  dp = rs
   REAL(DP), INTENT(IN)   ::  x
   LOGICAL,   INTENT(IN)  ::  upper
   REAL(DP)               ::  fn_val

!  Local variables
   REAL(DP), PARAMETER   ::  zero=0.0_DP, one=1.0_DP, half=0.5_DP, con=1.28_DP
   REAL(DP)              ::  z, y
   LOGICAL               ::  up

!  Machine dependent constants
   REAL(DP), PARAMETER  ::  ltone = 7.0_DP, utzero = 18.66_DP
   REAL(DP), PARAMETER  ::  p = 0.398942280444_DP, q = 0.39990348504_DP,   &
                            r = 0.398942280385_DP, a1 = 5.75885480458_DP,  &
                            a2 = 2.62433121679_DP, a3 = 5.92885724438_DP,  &
                            b1 = -29.8213557807_DP, b2 = 48.6959930692_DP, &
                            c1 = -3.8052E-8_DP, c2 = 3.98064794E-4_DP,     &
                            c3 = -0.151679116635_DP, c4 = 4.8385912808_DP, &
                            c5 = 0.742380924027_DP, c6 = 3.99019417011_DP, &
                            d1 = 1.00000615302_DP, d2 = 1.98615381364_DP,  &
                            d3 = 5.29330324926_DP, d4 = -15.1508972451_DP, &
                            d5 = 30.789933034_DP

   up = upper
   z = x
   IF( z < zero ) THEN
      up = .NOT. up
      z = -z
   END IF
   IF( z <= ltone  .OR.  (up  .AND.  z <= utzero) ) THEN
      y = half*z*z
      IF( z > con ) THEN
         fn_val = r*EXP( -y )/(z+c1+d1/(z+c2+d2/(z+c3+d3/(z+c4+d4/(z+c5+d5/(z+c6))))))
      ELSE
         fn_val = half - z*(p-q*y/(y+a1+b1/(y+a2+b2/(y+a3))))
      END IF
   ELSE
      fn_val = zero
   END IF

   IF( .NOT. up ) fn_val = one - fn_val
   RETURN
END FUNCTION alnorm
!-------------------------------------------------------------------------------

!===============================================================================
!-------------------------------------------------------------------------------
!  Circular statistical functions
!-------------------------------------------------------------------------------
!===============================================================================

!===============================================================================
function circ_mean(angle,degrees)
!===============================================================================
!  Returns the circular mean of a set of angles.  Input is a column vector
!  of arbitrary length.  Default is for input in degrees.
   implicit none
   real(rs),intent(in) :: angle(:)
   real(rs)            :: circ_mean
   logical,intent(in),optional :: degrees
   logical                     :: radians

   radians = .false.
   if (present(degrees)) radians = .not.degrees

   if (radians) then
      circ_mean = atan2(sum( sin( angle(1:size(angle)) ) ), &
                        sum( cos( angle(1:size(angle)) ) ) )
   else
      circ_mean = atan2(sum( sin( angle(1:size(angle))*pi/180._rs ) ),  &
                        sum( cos( angle(1:size(angle))*pi/180._rs ) ) )
      circ_mean = circ_mean*180._rs/pi
   endif

   return
end function circ_mean
!-------------------------------------------------------------------------------

!===============================================================================
subroutine circ_mean_bootstrap(angle,mu,mu1,mu2,B,degrees)
!===============================================================================
!  Calculate a sample mean using a bootstrap technique.  Assumes a symmetric
!  distribution.  See Fisher, Statistical analysis of circular data, §4.4
!  Optionally specify number of bootstrap samples to take.
   implicit none
   real(rs),intent(in)  :: angle(:)
   real(rs),intent(out) :: mu
   real(rs),intent(out),optional :: mu1,mu2
   integer,intent(in),optional   :: B
   logical,intent(in),optional   :: degrees
   logical :: degrees_in

   degrees_in = .true.
   if (present(degrees)) degrees_in = degrees

   if (size(angle) <= 9) then  ! For small N, calculate all n**n samples
      call circ_mean_bootstrap_smallN(angle,mu,mu1,mu2,degrees=degrees_in)
   else
      write(lu_stderr,'(a)') 'statistical: circ_mean_bootstrap: only implemented for n<=9 at the moment.'
      stop
   endif

   return
end subroutine circ_mean_bootstrap
!-------------------------------------------------------------------------------

!===============================================================================
subroutine circ_mean_bootstrap_smallN(angle,mu,mu1,mu2,degrees,force)
!===============================================================================
!  For a small sample, calculate the bootstrap mean and confidence interval
!  by taking all possible subsamples of the data (= n**n).
!  If n==9, then n**n real*4s takes up ~1.5 GB.  Use this routine with care on
!  lesser or shared machines!  Hence a warning is displayed, unless the force=.true.
!  option is employed on the command line.
   implicit none
   real(rs),intent(in)  :: angle(:)
   real(rs),intent(out) :: mu,mu1,mu2
   logical,intent(in),optional :: degrees,force
   logical :: force_in
   real(rs) :: conversion
   real(r4) :: Csum,Ssum
   real(r4),allocatable :: mean(:)  ! Use single preision for storage of bootstrap samples
   real(rs),allocatable :: rangle(:)
   integer :: n,i,i1,i2,i3,i4,i5,i6,i7,i8,i9

!  We're potentially allocating <=1.5 GB of memory, so make sure we know what
!  we're doing.
   force_in = .false.
   if (present(force)) force_in = force
   if (.not.force_in) then
      if (size(angle) >= 8) then
         write(*,'(a,f0.0,a)') 'circ_mean_bootstrap_smallN is about to allocate ',&
                  8.*size(angle)**size(angle)/(1024.**2),' MB of memory.'
         write(*,'(a)') 'Hit return to continue, or ^C to abort.'
         read(*,*)
      endif
   endif

!  Check we're only testing with small sample size
   n = size(angle)
   if (n > 9 .or. n == 1) then
      write(lu_stderr,'(a)') &
      'statistical: circ_mean_bootstrap_smallN: number of samples too large for small N subroutine, or n is 1.'
      stop
   endif

!  Allocate memory for all n**n sample means: here we go!
   allocate(mean(n**n))

!  Convert angles into radians if necesary
   allocate(rangle(n))
   conversion = pi/180._rs
   if (present(degrees)) then
      if (.not.degrees) conversion = 1._rs
   endif
   rangle = conversion*angle

!  Calculate sample mean for all possible bootstrap samples
   i = 1
   if (n == 2) then
      do i1=1,n
       do i2=1,n
               Csum = cos(rangle(i1)) + cos(rangle(i2))
               Ssum = sin(rangle(i1)) + sin(rangle(i2))
               mean(i) = atan2(Ssum,Csum)
               i = i + 1
       enddo
      enddo

   else if (n == 3) then
      do i1=1,n
       do i2=1,n
        do i3=1,n
               Csum = cos(rangle(i1)) + cos(rangle(i2)) + cos(rangle(i3))
               Ssum = sin(rangle(i1)) + sin(rangle(i2)) + sin(rangle(i3))
               mean(i) = atan2(Ssum,Csum)
               i = i + 1
        enddo
       enddo
      enddo

   else if (n == 4) then
      do i1=1,n
       write(*,'(a,i0,a)') 'Executing ',i1,' of 8 passes.'
       do i2=1,n
        do i3=1,n
         do i4=1,n
                 Csum = cos(rangle(i1)) + cos(rangle(i2)) + cos(rangle(i3)) + &
                        cos(rangle(i4))
                 Ssum = sin(rangle(i1)) + sin(rangle(i2)) + sin(rangle(i3)) + &
                        sin(rangle(i4))
                 mean(i) = atan2(Ssum,Csum)
                 i = i + 1
         enddo
        enddo
       enddo
      enddo

   else if (n == 5) then
      do i1=1,n
       write(*,'(a,i0,a)') 'Executing ',i1,' of 8 passes.'
       do i2=1,n
        do i3=1,n
         do i4=1,n
          do i5=1,n
                 Csum = cos(rangle(i1)) + cos(rangle(i2)) + cos(rangle(i3)) + &
                        cos(rangle(i4)) + cos(rangle(i5))
                 Ssum = sin(rangle(i1)) + sin(rangle(i2)) + sin(rangle(i3)) + &
                        sin(rangle(i4)) + sin(rangle(i5))
                 mean(i) = atan2(Ssum,Csum)
                 i = i + 1
          enddo
         enddo
        enddo
       enddo
      enddo

   else if (n == 6) then
      do i1=1,n
       write(*,'(a,i0,a)') 'Executing ',i1,' of 8 passes.'
       do i2=1,n
        do i3=1,n
         do i4=1,n
          do i5=1,n
           do i6=1,n
                 Csum = cos(rangle(i1)) + cos(rangle(i2)) + cos(rangle(i3)) + &
                        cos(rangle(i4)) + cos(rangle(i5)) + cos(rangle(i6))
                 Ssum = sin(rangle(i1)) + sin(rangle(i2)) + sin(rangle(i3)) + &
                        sin(rangle(i4)) + sin(rangle(i5)) + sin(rangle(i6))
                 mean(i) = atan2(Ssum,Csum)
                 i = i + 1
           enddo
          enddo
         enddo
        enddo
       enddo
      enddo

   else if (n == 7) then
      do i1=1,n
       write(*,'(a,i0,a)') 'Executing ',i1,' of 8 passes.'
       do i2=1,n
        do i3=1,n
         do i4=1,n
          do i5=1,n
           do i6=1,n
            do i7=1,n
                 Csum = cos(rangle(i1)) + cos(rangle(i2)) + cos(rangle(i3)) + &
                        cos(rangle(i4)) + cos(rangle(i5)) + cos(rangle(i6)) + &
                        cos(rangle(i7))
                 Ssum = sin(rangle(i1)) + sin(rangle(i2)) + sin(rangle(i3)) + &
                        sin(rangle(i4)) + sin(rangle(i5)) + sin(rangle(i6)) + &
                        sin(rangle(i7))
                 mean(i) = atan2(Ssum,Csum)
                 i = i + 1
            enddo
           enddo
          enddo
         enddo
        enddo
       enddo
      enddo

   else if (n == 8) then
      do i1=1,n
       write(*,'(a,i0,a)') 'Executing ',i1,' of 8 passes.'
       do i2=1,n
        do i3=1,n
         do i4=1,n
          do i5=1,n
           do i6=1,n
            do i7=1,n
             do i8=1,n
                 Csum = cos(rangle(i1)) + cos(rangle(i2)) + cos(rangle(i3)) + &
                        cos(rangle(i4)) + cos(rangle(i5)) + cos(rangle(i6)) + &
                        cos(rangle(i7)) + cos(rangle(i8))
                 Ssum = sin(rangle(i1)) + sin(rangle(i2)) + sin(rangle(i3)) + &
                        sin(rangle(i4)) + sin(rangle(i5)) + sin(rangle(i6)) + &
                        sin(rangle(i7)) + sin(rangle(i8))
                 mean(i) = atan2(Ssum,Csum)
                 i = i + 1
             enddo
            enddo
           enddo
          enddo
         enddo
        enddo
       enddo
      enddo

   else if (n == 9) then
      do i1=1,n
       write(*,'(a,i0,a)') 'Executing ',i1,' of 9 passes.'
       do i2=1,n
        do i3=1,n
         do i4=1,n
          do i5=1,n
           do i6=1,n
            do i7=1,n
             do i8=1,n
              do i9=1,n
                 Csum = cos(rangle(i1)) + cos(rangle(i2)) + cos(rangle(i3)) + &
                        cos(rangle(i4)) + cos(rangle(i5)) + cos(rangle(i6)) + &
                        cos(rangle(i7)) + cos(rangle(i8)) + cos(rangle(i9))
                 Ssum = sin(rangle(i1)) + sin(rangle(i2)) + sin(rangle(i3)) + &
                        sin(rangle(i4)) + sin(rangle(i5)) + sin(rangle(i6)) + &
                        sin(rangle(i7)) + sin(rangle(i8)) + sin(rangle(i9))
                 mean(i) = atan2(Ssum,Csum)
                 i = i + 1
              enddo
             enddo
            enddo
           enddo
          enddo
         enddo
        enddo
       enddo
      enddo
   endif

   mu = sum(mean)/(conversion*n**n)

   return
end subroutine circ_mean_bootstrap_smallN
!-------------------------------------------------------------------------------

!===============================================================================
function circ_res_length(angle,degrees)
!===============================================================================
!  Returns the resultant length of the angles.  Input is a column vector of
!  arbitrary length.  Default is for input in degrees.
   implicit none
   real(rs),intent(in) :: angle(:)
   real(rs)            :: circ_res_length
   logical,intent(in),optional :: degrees
   logical                     :: radians

   radians = .false.
   if (present(degrees)) radians = .not.degrees

   if (radians) then
      circ_res_length = sqrt((sum(sin(angle(1:size(angle)))))**2 + &
                             (sum(cos(angle(1:size(angle)))))**2) / real(size(angle))
   else
      circ_res_length = sqrt((sum( sin(angle(1:size(angle))*pi/180._rs)) )**2 + &
                (sum( cos(angle(1:size(angle))*pi/180._rs)) )**2) / real(size(angle))
   endif

   return
end function circ_res_length
!-------------------------------------------------------------------------------

!===============================================================================
function circ_sd(angle,degrees)
!===============================================================================
!  Returns the circular standard deviation.  Input is a column vector of
!  arbitrary length.  Default is for input in degrees.
   implicit none
   real(rs),intent(in) :: angle(:)
   real(rs)            :: circ_sd
   logical,intent(in),optional :: degrees
   logical                     :: degrees_in

   degrees_in = .true.
   if (present(degrees)) degrees_in = degrees

   circ_sd = sqrt(-2._rs*log(circ_res_length(angle,degrees=degrees_in)))

   return
end function circ_sd
!===============================================================================

!===============================================================================
function circ_variance(angle,degrees)
!===============================================================================
!  Returns the circular variance of a set of angles.  Input is a column vector
!  of arbitrary length.  Default is for input in degrees.
   implicit none
   real(rs),intent(in) :: angle(:)
   real(rs)            :: circ_variance,mean
   logical,intent(in),optional :: degrees
   logical                     :: radians

   radians = .false.
   if (present(degrees)) radians = .not.degrees

   if (radians) then
      mean = circ_mean(angle,degrees=.false.)
      circ_variance = 1._rs - sum(cos(angle - mean))/real(size(angle))
   else
      mean = circ_mean(angle,degrees=.true.) * pi/180._rs
      circ_variance = 1._rs - sum(cos(angle*pi/180._rs - mean))/real(size(angle))
   endif

   return
end function circ_variance
!-------------------------------------------------------------------------------

!===============================================================================
function circ_uncent_trig_moment(angle,p,degrees)
!===============================================================================
!  Returns the uncentered pth trigonometric moment.  Input is a column vector of
!  arbitrary length.  Default is for input in degrees.
   implicit none
   real(rs),intent(in) :: angle(:)
   integer,intent(in)  :: p
   complex(cs)         :: circ_uncent_trig_moment
   real(rs) :: Cp,Sp
   integer  :: n
   logical,intent(in),optional :: degrees
   real(rs) :: conversion

!  Check for valid p
   if (p < 1) then
      write(lu_stderr,'(a)') 'statistical: circ_uncent_trig_moment: p must be > 0.'
      stop
   endif

!  Convert from degrees to radians unless otherwise
   conversion = pi/180._rs
   if (present(degrees)) then
      if (.not.degrees) conversion = 1._rs
   endif

   n = size(angle)
   Cp = sum(cos(real(p)*angle(1:n)*conversion))/real(n)
   Sp = sum(sin(real(p)*angle(1:n)*conversion))/real(n)

   circ_uncent_trig_moment = cmplx(Cp, Sp)

   return
end function circ_uncent_trig_moment
!-------------------------------------------------------------------------------

!===============================================================================
function circ_cent_trig_moment(angle,p,degrees)
!===============================================================================
!  Returns the centered pth trigonometric moment.  Input is a column vector of
!  arbitrary length.  Default is for input in degrees.
   implicit none
   real(rs),intent(in) :: angle(:)
   integer, intent(in) :: p
   complex(cs)         :: circ_cent_trig_moment
   real(rs) :: Cp,Sp,mean,conversion
   integer  :: n
   logical,intent(in),optional :: degrees

!  Check for valid p
   if (p < 1) then
      write(lu_stderr,'(a)') 'statistical: circ_cent_trig_moment: P must be > 0.'
      stop
   endif

!  Convert from degrees to radians unless otherwise desired
   conversion = pi/180._rs
   if (present(degrees)) then;  if (.not.degrees) conversion = 1._rs; endif

   n = size(angle)
   mean = circ_mean(angle,degrees=degrees)
   Cp = sum(cos(real(p)*conversion*(angle(1:n) - mean)))/real(n)
   Sp = sum(sin(real(p)*conversion*(angle(1:n) - mean)))/real(n)

   circ_cent_trig_moment = cmplx(Cp, Sp)

   return
end function circ_cent_trig_moment
!-------------------------------------------------------------------------------

!===============================================================================
function circ_dispers(angle,degrees)
!===============================================================================
!  Returns the circular dispersion of a set of angles.  Input is a column vector
!  of arbitrary length.  Default is for input in degrees.
   implicit none
   real(rs),intent(in) :: angle(:)
   real(rs)            :: circ_dispers
   real(rs) :: mean,rho2,R,conversion
   integer :: n
   logical,intent(in),optional :: degrees

!  By default, convert from degrees to radians for calculations
   conversion = pi/180._rs
   if (present(degrees)) then
      if (.not.degrees) conversion = 1._rs
   endif

   mean = circ_mean(angle,degrees=degrees)
   n = size(angle)
   R = sum(cos(conversion*(angle(1:n) - mean)))/real(n)
   rho2 = sum(cos(2._rs*conversion*(angle(1:n) - mean)))/real(n)
   circ_dispers = (1._rs - rho2)/(2._rs*R**2)

   return
end function circ_dispers
!-------------------------------------------------------------------------------

!===============================================================================
function circ_correl(a,b,degrees)
!===============================================================================
!  Returns the circular correlation between two sets of angles.  Input is two
!  column vectors of arbitrary (but the same) length.  Default is for input
!  in dgerees.
   implicit none
   real(rs),intent(in) :: a(:), b(:)
   real(rs)            :: circ_correl, amean, bmean
   logical,intent(in),optional :: degrees
   logical                     :: radians

!  Check arrays are same length
   if (size(a) /= size(b)) then
      write(0,'(a)') 'Error: statistical: circ_correl: Input arrays are not same length.'
      stop
   endif

   radians = .false.
   if (present(degrees)) radians = .not.degrees

   if (radians) then
      amean = circ_mean(a,degrees=.false.)
      bmean = circ_mean(a,degrees=.false.)
      circ_correl = sum( sin(a - amean) * sin(b - bmean) )/ &
            sqrt( sum( sin(a-amean)*sin(a-amean)*sin(b-bmean)*sin(b-bmean) ) )
   else
      amean = circ_mean(a,degrees=.true.) * pi/180._rs
      bmean = circ_mean(b,degrees=.true.) * pi/180._rs
      circ_correl = sum( sin(a*pi/180._rs - amean) * sin(b*pi/180._rs - bmean) )/ &
            sqrt( sum( sin(a*pi/180._rs-amean)*sin(a*pi/180._rs-amean) * &
                  sin(b*pi/180._rs-bmean)*sin(b*pi/180._rs-bmean) ) )
   endif

   return
end function circ_correl
!-------------------------------------------------------------------------------

!===============================================================================
function circ_von_mieses(theta,mu,kappa,degrees)
!===============================================================================
!  Returns the value of the von Mieses distribution for given mean and 'pointiness'
   implicit none

   real(rs),intent(in) :: theta,mu,kappa
   real(rs)            :: circ_von_mieses,conversion
   logical,intent(in),optional :: degrees

   conversion = pi/180._rs
   if (present(degrees)) then; if (.not.degrees) conversion = 1._rs; endif

   circ_von_mieses = exp(kappa*cos(conversion*(theta - mu)))/(2._rs*pi*BessI0(kappa))

   return

   contains
      function BessI0(x)
      !  Evaluates the modified Bessel function of the first kind at x.
      !  Polynomial approximation taken from Abramowitz & Stegun, 1964, Handbook
      !   of mathematical functions.  See also Numerical Recipes §6.6
         implicit none
         real(rs),intent(in) :: x
         real(rs)            :: BessI0
         real(r8)            :: t

         if (abs(x) < 3.75_r8) then
            t = abs(x/3.75_r8)
            BessI0 = 1._r8 + &
                     3.5156229_r8*t**2 + 3.0899424_r8*t**4 +  1.2067492_r8*t**6 + &
                     0.2659732_r8*t**8 + 0.0360768_r8*t**10 + 0.0045813_r8*t**12
         else
            t = abs(3.75_r8/x)
            BessI0 = (exp(abs(x))/sqrt(abs(x)))*(0.39894228_rs + &
                     0.01328592_r8*t    + 0.00225319_r8*t**2 - 0.00157565_r8*t**3 + &
                     0.00916281_r8*t**4 - 0.02057706_r8*t**5 + 0.02635537_r8*t**6 - &
                     0.01647633_r8*t**7 + 0.00392377_r8*t**8)
         endif
         return
      end function BessI0

end function circ_von_mieses
!-------------------------------------------------------------------------------

!===============================================================================
function circ_test_random_orient(theta,theta0,alpha,degrees) result(pass)
!===============================================================================
! Test whether a set of orientations (i.e., 180-degrees ambiguous), theta_i,
! point in a certain direction, theta0, statistically differently from random,
! using the V-test as modified by Rayleigh
! (test 95 in: 100 Statistical Tests. G.K. Kanji, SAGE Publications.).
   implicit none
   real(rs), intent(in) :: theta(:), theta0
   real(rs), intent(in), optional :: alpha
   logical, optional, intent(in) :: degrees
   real(rs) :: conversion, alpha_in
   real(rs) :: r, mean, V
   logical :: radians, pass
   integer :: n

   n = size(theta)
   if (n < 5) then
      write(0,'(a)') 'statistical: circ_test_random_orient: Error: Cannot ' // &
         'calculate significance for samples < 5'
      stop
   endif

   radians = .false.
   alpha_in = 0.05_rs
   if (present(alpha)) alpha_in = alpha

   if (present(degrees)) radians = .not.degrees
   conversion = pi/180._rs
   if (radians) conversion = 1._rs

   mean = circ_mean(2._rs*theta, degrees=.not.radians)
   r = circ_res_length(2._rs*theta, degrees=.not.radians)
   V = 2._rs*sqrt(2._rs*size(theta))*r*cos(conversion*mean - 2._rs*conversion*theta0)

   ! Lookup values in table and decide on significance
   pass = .false.
   if (V > v_test_table()) pass = .true.

   contains
      function v_test_table() result(value)
         integer :: i
         real(rs) :: value
         integer, parameter, dimension(33) :: table_n = &
            (/ 5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26, &
            27,28,29,30,40,50,60,70,100,500,1000 /)
         real(rs), parameter, dimension(6) :: table_alpha = &
            (/ 0.1, 0.05, 0.01, 0.005, 0.001, 0.0001 /)
         real(rs), parameter :: table(33,6) = transpose(reshape( (/ &
            1.3051_rs, 1.6524_rs, 2.2505_rs, 2.4459_rs, 2.7938_rs, 3.0825_rs, &
            1.3009_rs, 1.6509_rs, 2.2640_rs, 2.4695_rs, 2.8502_rs, 3.2114_rs, &
            1.2980_rs, 1.6499_rs, 2.2734_rs, 2.4858_rs, 2.8886_rs, 3.2970_rs, &
            1.2958_rs, 1.6492_rs, 2.2803_rs, 2.4978_rs, 2.9164_rs, 3.3578_rs, &
            1.2942_rs, 1.6484_rs, 2.2856_rs, 2.5070_rs, 2.9375_rs, 3.4034_rs, &
            1.2929_rs, 1.6482_rs, 2.2899_rs, 2.5143_rs, 2.9540_rs, 3.4387_rs, &
            1.2918_rs, 1.6479_rs, 2.2933_rs, 2.5201_rs, 2.9672_rs, 3.4669_rs, &
            1.2909_rs, 1.6476_rs, 2.2961_rs, 2.5250_rs, 2.9782_rs, 3.4899_rs, &
            1.2902_rs, 1.6474_rs, 2.2985_rs, 2.5290_rs, 2.9873_rs, 3.5091_rs, &
            1.2895_rs, 1.6472_rs, 2.3006_rs, 2.5325_rs, 2.9950_rs, 3.5253_rs, &
            1.2890_rs, 1.6470_rs, 2.3023_rs, 2.5355_rs, 3.0017_rs, 3.5392_rs, &
            1.2885_rs, 1.6469_rs, 2.3039_rs, 2.5381_rs, 3.0075_rs, 3.5512_rs, &
            1.2881_rs, 1.6467_rs, 2.3052_rs, 2.5404_rs, 3.0126_rs, 3.5617_rs, &
            1.2877_rs, 1.6466_rs, 2.3064_rs, 2.5424_rs, 3.0171_rs, 3.5710_rs, &
            1.2874_rs, 1.6465_rs, 2.3075_rs, 2.5442_rs, 3.0211_rs, 3.5792_rs, &
            1.2871_rs, 1.6464_rs, 2.3085_rs, 2.5458_rs, 3.0247_rs, 3.5866_rs, &
            1.2868_rs, 1.6464_rs, 2.3093_rs, 2.5473_rs, 3.0279_rs, 3.5932_rs, &
            1.2866_rs, 1.6463_rs, 2.3101_rs, 2.5486_rs, 3.0308_rs, 3.5992_rs, &
            1.2864_rs, 1.6462_rs, 2.3108_rs, 2.5498_rs, 3.0335_rs, 3.6047_rs, &
            1.2862_rs, 1.6462_rs, 2.3115_rs, 2.5509_rs, 3.0359_rs, 3.6096_rs, &
            1.2860_rs, 1.6461_rs, 2.3121_rs, 2.5519_rs, 3.0382_rs, 3.6142_rs, &
            1.2858_rs, 1.6461_rs, 2.3127_rs, 2.5529_rs, 3.0402_rs, 3.6184_rs, &
            1.2856_rs, 1.6460_rs, 2.3132_rs, 2.5538_rs, 3.0421_rs, 3.6223_rs, &
            1.2855_rs, 1.6460_rs, 2.3136_rs, 2.5546_rs, 3.0439_rs, 3.6258_rs, &
            1.2853_rs, 1.6459_rs, 2.3141_rs, 2.5553_rs, 3.0455_rs, 3.6292_rs, &
            1.2852_rs, 1.6459_rs, 2.3145_rs, 2.5560_rs, 3.0471_rs, 3.6323_rs, &
            1.2843_rs, 1.6456_rs, 2.3175_rs, 2.5610_rs, 3.0580_rs, 3.6545_rs, &
            1.2837_rs, 1.6455_rs, 2.3193_rs, 2.5640_rs, 3.0646_rs, 3.6677_rs, &
            1.2834_rs, 1.6454_rs, 2.3205_rs, 2.5660_rs, 3.0689_rs, 3.6764_rs, &
            1.2831_rs, 1.6453_rs, 2.3213_rs, 2.5674_rs, 3.0720_rs, 3.6826_rs, &
            1.2826_rs, 1.6452_rs, 2.3228_rs, 2.5699_rs, 3.0775_rs, 3.6936_rs, &
            1.2818_rs, 1.6449_rs, 2.3256_rs, 2.5747_rs, 3.0877_rs, 3.7140_rs, &
            1.2817_rs, 1.6449_rs, 2.3260_rs, 2.5752_rs, 3.0890_rs, 3.7165_rs /), (/6,33/)))
         ! The function just checks that the value you asked for alpha
         ! is near to one of those we offer--if not, you get an error
         do i = 1, size(table_alpha)
            if (abs(table_alpha(i) - alpha_in)/table_alpha(i) < 0.1_rs) then
               value = table(minloc(abs(n-table_n),1), i)
               return
            endif
         enddo
         write(0,'(a,f0.6,a)') 'statistical: circ_test_random_orient: v_test_table: ' &
            //'Error: Cannot find a value for alpha near that requested (',alpha_in,')'
         stop
      end function v_test_table
end function circ_test_random_orient
!-------------------------------------------------------------------------------

!===============================================================================
!-------------------------------------------------------------------------------
!  Spherical statistical functions
!-------------------------------------------------------------------------------
!===============================================================================

!===============================================================================
function sph_mean(x, y, z) result(mean)
!===============================================================================
! Compute the mean direction for a set of vectors representing the x, y and z
! directions of points on the unit sphere.  Each input value is normalised to
! the unit sphere even if it isn't on input.
! Output is a three-vector giving the direction as (x,y,z)
   real(rs), intent(in), dimension(:) :: x, y, z
   real(rs) :: mean(3)
   if (size(x) /= size(y) .or. size(y) /= size(z) .or. size(z) /= size(x)) then
      write(0,'(a)') 'statistical: sph_mean: Input arrays must be same length'
      stop
   endif
   mean = (/sum(x), sum(y), sum(z)/)/sph_resultant(x, y, z)
end function sph_mean
!-------------------------------------------------------------------------------

!===============================================================================
function sph_resultant(x, y, z) result(resultant)
!===============================================================================
! Compute the spherical resultant
   real(rs), intent(in), dimension(:) :: x, y, z
   real(rs) :: resultant(3)
   if (size(x) /= size(y) .or. size(y) /= size(z) .or. size(z) /= size(x)) then
      write(0,'(a)') 'statistical: sph_resultant: Input arrays must be same length'
      stop
   endif
   resultant = sqrt(sum(x)**2 + sum(y)**2 + sum(z)**2)
end function sph_resultant
!-------------------------------------------------------------------------------

!===============================================================================
function sph_mean_resultant(x, y, z) result(mean_resultant)
!===============================================================================
   real(rs), intent(in), dimension(:) :: x, y, z
   real(rs) :: mean_resultant(3)
   mean_resultant = sph_resultant(x, y, z)/size(x)
end function sph_mean_resultant
!-------------------------------------------------------------------------------

!===============================================================================
!-------------------------------------------------------------------------------
!  Linear statistical functions
!-------------------------------------------------------------------------------
!===============================================================================

!===============================================================================
function stat_mean(x)
!===============================================================================
! Compute the mean
   real(rs), intent(in), dimension(:) :: x
   real(rs) :: stat_mean
   stat_mean = sum(x)/size(x)
end function stat_mean
!-------------------------------------------------------------------------------

!===============================================================================
function stat_variance(x)
!===============================================================================
! Compute the variance
   real(rs), intent(in), dimension(:) :: x
   real(rs) :: stat_variance
   stat_variance = sum((x - stat_mean(x))**2)/(size(x) - 1)
end function stat_variance
!-------------------------------------------------------------------------------

!===============================================================================
function stat_sd(x)
!===============================================================================
! Compute the standard deviation
   real(rs), intent(in), dimension(:) :: x
   real(rs) :: stat_sd
   stat_sd = sqrt(stat_variance(x))
end function stat_sd
!-------------------------------------------------------------------------------

!===============================================================================
function stat_covariance(x, y)
!===============================================================================
! Compute the covariance between two arrays of the same length
   real(rs), intent(in), dimension(:) :: x, y
   real(rs) :: stat_covariance
   if (size(x) /= size(y)) then
      write(0,'(a)') 'statistical: stat_covariance: Input arrays must be same length'
      stop
   endif
   stat_covariance = sum((x - stat_mean(x))*(y - stat_mean(y)))/size(x)
end function stat_covariance
!-------------------------------------------------------------------------------
!_______________________________________________________________________________
end module statistical

max_stars_repo_name: JackWalpole/seismo-fortran-fork

max_issues_repo_name: JackWalpole/seismo-fortran-fork

max_forks_repo_name: JackWalpole/seismo-fortran-fork


program main

  use heat
  implicit none

  type (heat_model) :: model

  write(*,"(a)") "Start"
  call initialize_from_defaults(model)

  write(*,"(a)") "Model initial values..."
  call print_values(model)

  do while (model%t < model%t_end)
     call advance_in_time(model)
     write(*,"(a, f6.1)") "Model values at time = ", model%t
     call print_values(model)
  end do

  call cleanup(model)
  write(*,"(a)") "Finish"

end program main


text:
      SUBROUTINE REGRID
     &           (NSD, NNW_SD, NGP_SD, NW_IPAR, NW_RPAR, CRD)
C ---------------------------------------------------------------------------
C
C ---------------------------------------------------------------------------
      IMPLICIT NONE
      INCLUDE 'knd_params.inc'
      INCLUDE 'net_params.inc'
C
      INTEGER(KIND=IK) NSD
      INTEGER(KIND=IK) NNW_SD (*)
      INTEGER(KIND=IK) NGP_SD (*)
      INTEGER(KIND=IK) NW_IPAR (N_IP, *)
      REAL   (KIND=RK) NW_RPAR (N_RP, *)
      REAL   (KIND=RK) CRD (4, *)
C
      INTEGER(KIND=IK) ISD, INW, IGP, NNW, NGP
C
      INW = 1
      IGP = 1
C
      DO ISD = 1, NSD
         NNW = NNW_SD (ISD)
         NGP = NGP_SD (ISD)
         CALL REGRID_SD
     &        (NNW, NW_IPAR (1, INW), NW_RPAR (1, INW), CRD (1, IGP))
         INW = INW + NNW
         IGP = IGP + NGP
      END DO
C
      RETURN
      END
      SUBROUTINE REGRID_SD
     &           (NNW, NW_IPAR, NW_RPAR, CRD)
C ---------------------------------------------------------------------------
C
C ---------------------------------------------------------------------------
      IMPLICIT NONE
      INCLUDE 'knd_params.inc'
      INCLUDE 'net_params.inc'
      INCLUDE 'net_funcs.inc'
C
      INTEGER(KIND=IK) NNW
      INTEGER(KIND=IK) NW_IPAR (N_IP, *)
      REAL   (KIND=RK) NW_RPAR (N_RP, *)
      REAL   (KIND=RK) CRD (4, *)
C
      INTEGER(KIND=IK) INW, IGP, NGP
C
      IGP = 1
C
      DO INW = 1, NNW
         NGP = GET_NGP (NW_IPAR (1, INW))
         CALL REGRID_NW
     &        (NGP, NW_IPAR (1, INW), NW_RPAR (1, INW), CRD (1, IGP))
         IGP = IGP + NGP
      END DO
C
      RETURN
      END
      SUBROUTINE REGRID_NW
     &           (NGP, NW_IPAR, NW_RPAR, CRD)
C ---------------------------------------------------------------------------
C
C ---------------------------------------------------------------------------
      IMPLICIT NONE
      INCLUDE 'knd_params.inc'
      INCLUDE 'net_params.inc'
      INCLUDE 'anl_params.inc'
      INCLUDE 'tme_params.inc'
      INCLUDE 'net_funcs.inc'
C
      INTEGER(KIND=IK) NGP
      INTEGER(KIND=IK) NW_IPAR (N_IP)
      REAL   (KIND=RK) NW_RPAR (N_RP)
      REAL   (KIND=RK) CRD (4, *)
C
      INTEGER(KIND=IK) BCT
C
      IF (FIRST_TIME .AND. ANL_TYPE .EQ. ANL_RFWAVE) THEN
         BCT = GET_BCT (NW_IPAR)
         CALL REGRID_RF (BCT, NGP, ANL_PAR, ANL_PAR (67), CRD)
      END IF
C
      RETURN
      END
      SUBROUTINE REGRID_RF (BCT, N, A, B, CRD)
C ---------------------------------------------------------------------------
C
C ---------------------------------------------------------------------------
      IMPLICIT NONE
      INCLUDE 'knd_params.inc'
      INCLUDE 'bct_params.inc'
      INCLUDE 'tme_params.inc'
      INCLUDE 'tme_funcs.inc'
C
      INTEGER(KIND=IK) BCT
      INTEGER(KIND=IK) N
      REAL   (KIND=RK) A (*)
      REAL   (KIND=RK) B (*)
      REAL   (KIND=RK) CRD (4, N)
C
      INTEGER(KIND=IK) I, J, NRF
      REAL   (KIND=RK) L, H, X0, W, K
      REAL   (KIND=RK) X, Z, ETA, ETA_X, T, C
      REAL   (KIND=RK) Z_XI (N)
      INTEGER(KIND=IK) BC (2)
      DATA             BC /0, 0/
C
      T   = GET_TIME ()
      L   = A (2)
      H   = A (3)
      X0  = A (4)
      W   = A (5)
      K   = A (6)
      C   = W / K
      NRF = INT (A (8))
C
      DO I = 1, N
         X     = CRD (1, I) - X0 - C * T
         Z     = CRD (2, I)
         ETA   = 0.5D0 * B (1)
         ETA_X = 0.0D0
         DO J = 1, NRF
           ETA   = ETA    +         B (J + 1) * COS (J * K * X)
           ETA_X = ETA_X  - J * K * B (J + 1) * SIN (J * K * X)
         END DO
         CRD (2, I) = -H + (H + ETA)  * (Z + H) / H
         Z_XI   (I) = CRD (4, I)
      END DO
C
      IF (BCT .EQ. BCT_BERNOULLI) THEN
         DO I = 1, N
            CRD (4, I) = Z_XI (I)
         END DO
      ELSE
         CALL SPLINE (BC, N, 1, 4, CRD (2, 1), CRD (4, 1))
C        WRITE (*, '(1X,4E14.6)') CRD
      END IF
C
      RETURN
      END

ext: f

max_stars_repo_name: meggermo/waterwaves

max_issues_repo_name: meggermo/waterwaves

max_forks_repo_name: meggermo/waterwaves


text:
module interface_f_mod

  use, intrinsic :: iso_c_binding
  use, intrinsic :: iso_fortran_env
  use :: flcl_mod

  implicit none

  public

  interface
    subroutine f_kokkos_initialize() &
      & bind(c, name='c_kokkos_initialize')
      use, intrinsic :: iso_c_binding
      implicit none
    end subroutine f_kokkos_initialize
  end interface

  interface
    subroutine f_kokkos_finalize() &
      & bind(c, name='c_kokkos_finalize')
      use, intrinsic :: iso_c_binding
      implicit none
    end subroutine f_kokkos_finalize
  end interface

  interface
    subroutine f_refine_mesh_status( nd_new_mesh_status, nd_mesh_status, dimension, refine_count ) &
      & bind(c, name="c_refine_mesh_status")
      use, intrinsic :: iso_c_binding
      use :: flcl_mod
      implicit none
      type(nd_array_t), intent(inout) :: nd_mesh_status
      type(nd_array_t), intent(inout) :: nd_new_mesh_status
      integer(c_int32_t), intent(inout) :: dimension
      integer(c_int32_t), intent(inout) :: refine_count
    end subroutine f_refine_mesh_status
  end interface

  interface
    subroutine f_axpy_full( v_x, v_y, s_val, c_levs ) &
      & bind(c, name="c_axpy_full")
      use, intrinsic :: iso_c_binding
      use :: flcl_mod
      use :: levels_f_mod
      implicit none
      type(c_ptr), intent(inout) :: v_x
      type(c_ptr), intent(inout) :: v_y
      real(c_double), intent(inout) :: s_val
      type(c_levels_t), intent(inout) :: c_levs
    end subroutine f_axpy_full
  end interface

  interface
    subroutine f_axpy_ltop( v_x, v_y, s_val, c_levs ) &
      & bind(c, name="c_axpy_ltop")
      use, intrinsic :: iso_c_binding
      use :: flcl_mod
      use :: levels_f_mod
      implicit none
      type(c_ptr), intent(inout) :: v_x
      type(c_ptr), intent(inout) :: v_y
      real(c_double), intent(inout) :: s_val
      type(c_levels_t), intent(inout) :: c_levs
    end subroutine f_axpy_ltop
  end interface

  interface
    subroutine f_axpy_ltop_chunks( v_x, v_y, s_val, c_levs ) &
      & bind(c, name="c_axpy_ltop_chunks")
      use, intrinsic :: iso_c_binding
      use :: flcl_mod
      use :: levels_f_mod
      implicit none
      type(c_ptr), intent(inout) :: v_x
      type(c_ptr), intent(inout) :: v_y
      real(c_double), intent(inout) :: s_val
      type(c_levels_t), intent(inout) :: c_levs
    end subroutine f_axpy_ltop_chunks
  end interface

  contains

    subroutine kokkos_initialize()
      use, intrinsic :: iso_c_binding
      implicit none
      call f_kokkos_initialize()
    end subroutine kokkos_initialize

    subroutine kokkos_finalize()
      use, intrinsic :: iso_c_binding
      implicit none
      call f_kokkos_finalize()
    end subroutine kokkos_finalize

    subroutine refine_mesh_status( p_mesh_status, p_new_mesh_status, dimension, refine_count )

      use, intrinsic :: iso_c_binding
      use :: flcl_mod
      implicit none
      integer(c_int32_t), dimension(:), intent(inout), pointer :: p_mesh_status
      integer(c_int32_t), dimension(:), intent(inout), pointer :: p_new_mesh_status
      integer(c_int32_t), intent(inout) :: dimension
      integer(c_int32_t), intent(inout) :: refine_count

      ! integer(c_int32_t), dimension(:), pointer :: p_temp_mesh_status
      type(nd_array_t) :: nd_mesh_status
      type(nd_array_t) :: nd_new_mesh_status

      nd_mesh_status = to_nd_array( p_mesh_status )
      nd_new_mesh_status = to_nd_array( p_new_mesh_status )
      ! write(*,*)'refine_mesh_status: mesh_status_dims(1) ', mesh_status_dims(1)
      call f_refine_mesh_status( nd_mesh_status, nd_new_mesh_status, dimension, refine_count )

      ! p_temp_mesh_status => p_mesh_status
      ! p_mesh_status => p_new_mesh_status
      ! nullify(p_temp_mesh_status)

    end subroutine refine_mesh_status

    subroutine axpy_full( v_x, v_y, s_val, c_levs )
      use, intrinsic :: iso_c_binding
      use :: flcl_mod
      use :: levels_f_mod
      implicit none
      real(c_double), intent(inout) :: s_val
      type(c_levels_t), intent(inout) :: c_levs
      type(c_ptr), intent(inout) :: v_x
      type(c_ptr), intent(inout) :: v_y

      call f_axpy_full( v_x, v_y, s_val, c_levs )
    end subroutine axpy_full

    subroutine axpy_ltop( v_x, v_y, s_val, c_levs )
      use, intrinsic :: iso_c_binding
      use :: flcl_mod
      use :: levels_f_mod
      implicit none
      real(c_double), intent(inout) :: s_val
      type(c_levels_t), intent(inout) :: c_levs

      type(c_ptr), intent(inout) :: v_x
      type(c_ptr), intent(inout) :: v_y

      call f_axpy_ltop( v_x, v_y, s_val, c_levs )
    end subroutine axpy_ltop

    subroutine axpy_ltop_chunks( v_x, v_y, s_val, c_levs )
      use, intrinsic :: iso_c_binding
      use :: flcl_mod
      use :: levels_f_mod
      implicit none
      real(c_double), intent(inout) :: s_val
      type(c_levels_t), intent(inout) :: c_levs

      type(c_ptr), intent(inout) :: v_x
      type(c_ptr), intent(inout) :: v_y

      call f_axpy_ltop_chunks( v_x, v_y, s_val, c_levs )
    end subroutine axpy_ltop_chunks

end module interface_f_mod

max_stars_repo_name: lanl/xkt

max_issues_repo_name: lanl/xkt

max_forks_repo_name: lanl/xkt


text:
Module mem_grid

real   , allocatable, dimension(:,:) :: lat, lon, topo, ldf
integer, allocatable, dimension(:,:) :: IstartIend
integer, allocatable, dimension(:)   :: recvcounts,displs
integer                              :: nPEs,rank

end Module

max_stars_repo_name: maxinye/laps-mirror

max_stars_repo_licenses: Intel Unlicense OLDAP-2.2.1 NetCDF

max_issues_repo_name: longwosion/laps-mirror

max_issues_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense

max_forks_repo_name: longwosion/laps-mirror

max_forks_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense


text:
      SUBROUTINE TRITRA(IND,CONV,NPH,NPV,LINE,SAMP,MODE)
C	THIS IS THE UCL VERSION OF TRITRA
C  TRANSFORM COORDINATES USING GEOMA PARAMETERS
C  CONV=GEOMA PARAMETERS STARTING AFTER 'TIEP' KEYWORD
C  NPH,NPV NUMBER POINTS HORIZONTAL,VERTICAL
C  LINE,SAMP  R*4 REPLACED BY OUTPUT VALUES
C  MODE=0 NL,NS TO OL,OS  MODE=1 OL,OS TO NL,NS
      REAL LINE,CONV(4,NPH,NPV),CORNER(2,4),POINT(2),A(16),B(4)
      CHARACTER*80 MSG1
      LOGICAL   INSIDE
C==================================================================
      II=2
      LL=1
      IITO=4
      LLTO=3
      IF(MODE.EQ.0) GO TO 5
      II=4
      LL=3
      IITO=2
      LLTO=1
5     CONTINUE
      POINT(1)=SAMP
      POINT(2)=LINE
      IND=0
      NAH=NPH-1
      NAV=NPV-1
C
C  FIND CLOSEST TIEPOINT TO LINE,SAMP
C  LOOP EACH ROW POINTS
      RMAX=1.0E+20
      SAMOLD=RMAX
      DO 90 L=1,NPV
C  LOOP EACH POINT IN ROW L
      DO 91 I=1,NPH
      SAMNEW=CONV(II,I,L)
      IF(SAMNEW.EQ.SAMOLD) GO TO 91
      R=(LINE-CONV(LL,I,L))**2+(SAMP-SAMNEW)**2
      IF(R.GT.RMAX) GO TO 92
      RMAX=R
      ISAV=I
      LSAV=L
92    SAMOLD=SAMNEW
91    CONTINUE
90    CONTINUE
C
C  SEARCH TWENTY QUADRILATERALS AROUND CLOSEST POINT
      ILEFT=ISAV-2
      IF(ILEFT.LT.1) ILEFT=1
      IRIGHT=ISAV+2
      IF(IRIGHT.GT.NAH) IRIGHT=NAH
      LTOP=LSAV-2
      IF(LTOP.LT.1) LTOP=1
      LBOT=LSAV+1
      IF(LBOT.GT.NAV) LBOT=NAV
C  AREA ROW LOOP
      DO 10 L=LTOP,LBOT
C  AREA COLUMN LOOP
      DO 20 I=ILEFT,IRIGHT
      N=1
C  LOAD UPPER LEFT CORNER
      CORNER(1,N)=CONV(II,I,L)
      CORNER(2,N)=CONV(LL,I,L)
C  CHECK IF POINTS ARE SAME
      IF(CORNER(1,N).EQ.CONV(II,I+1,L)) GO TO 40
      N=N+1
      CORNER(1,N)=CONV(II,I+1,L)
      CORNER(2,N)=CONV(LL,I+1,L)
40    CONTINUE
C  LOAD LOWER RIGHT POINT
      N=N+1
      CORNER(1,N)=CONV(II,I+1,L+1)
      CORNER(2,N)=CONV(LL,I+1,L+1)
C  CHECK IF POINTS ARE SAME
      IF(CONV(II,I,L+1).EQ.CORNER(1,N)) GO TO 41
      N=N+1
      CORNER(1,N)=CONV(II,I,L+1)
      CORNER(2,N)=CONV(LL,I,L+1)
41    CONTINUE
C  CHECK FOR DEGENERACY
      IF(N.GT.2) GO TO 43
      IND=1
      write (MSG1, 9000) I, L
9000  FORMAT('DEGENERATE QUADRILATERAL AT AH=',I3,'    AV=',I3)
      CALL XVMESSAGE(MSG1,' ')
      RETURN
43    CONTINUE
C  CHECK IF INSIDE CURRENT QUADRILATERAL
      IF(INSIDE(POINT,CORNER,N)) GO TO 50
20    CONTINUE
10    CONTINUE
C
C  POINT OUTSIDE QUADRILATERAL ARRAY - SET TO ONE OF 4 QUADS.
      I=ISAV
      L=LSAV
      IF(LINE.GT.CONV(LL,I,L)) GO TO 85
      IF(SAMP.GT.CONV(II,I,L)) GO TO 83
      I=I-1
      L=L-1
      IF(L.LT.1) L=1
      IF(I.LT.1) I=1
      GO TO 89
C  SAMP IS TO RIGHT OF ISAV
83    L=L-1
      IF(L.LT.1) L=1
      IF(I.EQ.NPH) GO TO 89
      IF(CONV(II,I,L).EQ.CONV(II,I+1,L)) I=I+1
      GO TO 89
C  LINE IS BELOW LSAV
85    IF(SAMP.GT.CONV(II,I,L)) GO TO 87
      I=I-1
      IF(I.LT.1) I=1
      GO TO 89
87    IF(I.EQ.NPH) GO TO 89
      IF(CONV(II,I,L).EQ.CONV(II,I+1,L)) I=I+1
89    IF(I.GT.NAH) I=NAH
      IF(L.GT.NAV) L=NAV
C  CHECK IF IS QUADRILATERAL OR REALLY A TRIANGLE (WHAT A PAIN)
      N=3
      IF(CONV(II,I,L).EQ.CONV(II,I+1,L)) GO TO 82
      IF(CONV(II,I,L+1).EQ.CONV(II,I+1,L+1)) GO TO 82
      N=4
82    CONTINUE
C
C  SOLVE SIMULTANEOUS EQUATIONS X,Y=A+BX+CY+DXY
50    IF(N.LT.4) GO TO 60
C  QUADRILATERALS ONLY
      DO 62 M=LLTO,IITO
      DO 61 J=1,4
61    A(J)=1.0
      A(5)=CONV(II,I,L)
      A(6)=CONV(II,I+1,L)
      A(7)=CONV(II,I,L+1)
      A(8)=CONV(II,I+1,L+1)
      A(9)=CONV(LL,I,L)
      A(10)=CONV(LL,I+1,L)
      A(11)=CONV(LL,I,L+1)
      A(12)=CONV(LL,I+1,L+1)
      A(13)=A(5)*A(9)
      A(14)=A(6)*A(10)
      A(15)=A(7)*A(11)
      A(16)=A(8)*A(12)
      B(1)=CONV(M,I,L)
      B(2)=CONV(M,I+1,L)
      B(3)=CONV(M,I,L+1)
      B(4)=CONV(M,I+1,L+1)
      CALL SIMQ(A,B,N,IND)
      IF(IND.NE.0) CALL XVMESSAGE('SIMQ MATRIX IS SINGULAR',' ')
      IF(IND.NE.0) RETURN
      IF(M.EQ.IITO) GO TO 63
      XLINE=B(1)+B(2)*SAMP+B(3)*LINE+B(4)*LINE*SAMP
      GO TO 62
63    SAMP=B(1)+B(2)*SAMP+B(3)*LINE+B(4)*LINE*SAMP
      LINE=XLINE
62    CONTINUE
      RETURN
C
C  TRIANGLES ONLY
60    DO 65 M=LLTO,IITO
      DO 66 J=1,3
66    A(J)=1.0
      A(4)=CONV(II,I,L)
      A(7)=CONV(LL,I,L)
      B(1)=CONV(M,I,L)
      IF(CONV(II,I+1,L).EQ.A(4)) GO TO 67
      A(5)=CONV(II,I+1,L)
      A(8)=CONV(LL,I+1,L)
      B(2)=CONV(M,I+1,L)
67    A(6)=CONV(II,I,L+1)
      A(9)=CONV(LL,I,L+1)
      B(3)=CONV(M,I,L+1)
      IF(CONV(II,I+1,L+1).EQ.A(6)) GO TO 68
      A(5)=CONV(II,I+1,L+1)
      A(8)=CONV(LL,I+1,L+1)
      B(2)=CONV(M,I+1,L+1)
68    CALL SIMQ(A,B,N,IND)
      IF(IND.NE.0) CALL XVMESSAGE('SIMQ MATRIX IS SINGULAR',' ')
      IF(IND.NE.0) RETURN
      IF(M.EQ.IITO) GO TO 69
      XLINE=B(1)+B(2)*SAMP+B(3)*LINE
      GO TO 65
69    SAMP=B(1)+B(2)*SAMP+B(3)*LINE
      LINE=XLINE
65    CONTINUE
      RETURN
      END

max_stars_repo_name: NASA-AMMOS/VICAR

max_issues_repo_name: NASA-AMMOS/VICAR

max_forks_repo_name: NASA-AMMOS/VICAR


text:
C$Procedure      SETERR ( Set Error Status )

      LOGICAL FUNCTION SETERR ( STATUS )

C$ Abstract
C
C     Set the SPICELIB error status.  DO NOT CALL THIS ROUTINE.
C
C$ Disclaimer
C
C     THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE
C     CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S.
C     GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE
C     ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE
C     PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS"
C     TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY
C     WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A
C     PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC
C     SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE
C     SOFTWARE AND RELATED MATERIALS, HOWEVER USED.
C
C     IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA
C     BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT
C     LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND,
C     INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS,
C     REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE
C     REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY.
C
C     RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF
C     THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY
C     CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE
C     ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE.
C
C$ Required_Reading
C
C     ERROR
C
C$ Keywords
C
C     ERROR
C
C$ Declarations

       LOGICAL               STATUS

C$ Brief_I/O
C
C      VARIABLE  I/O  DESCRIPTION
C      --------  ---  --------------------------------------------------
C      STATUS     I   Status indicator.
C
C
C      The function takes an UNSPECIFIED (and meaningless) value
C      on exit.
C
C$ Detailed_Input
C
C      STATUS  Indicates the new status.  When .TRUE., it
C              means that an error condition exists.
C
C$ Detailed_Output
C
C      None.
C
C      This purpose of this routine is to set status; the
C      function takes an UNSPECIFIED value on exit.  The
C      assigned value does not have any meaning.
C
C$ Parameters
C
C      None.
C
C$ Exceptions
C
C      None.
C
C$ Files
C
C      None.
C
C$ Particulars
C
C      DO NOT CALL THIS ROUTINE.
C
C      This is a data structure access routine for the
C      SPICELIB status.  This routine should be used for no
C      other purpose; in particular, it should not be used
C      to signal errors.  Use SIGERR or FAILED for that.
C
C      This routine assigns a value to SETERR on exit.
C      However, the value is not meaningful.
C
C$ Examples
C
C      None.  DON'T CALL THIS ROUTINE.
C
C      No examples.  If you don't know EXACTLY what a
C      ``data structure access routine'' is, don't call
C      this routine.  If you do know, you don't need an
C      example.
C
C$ Restrictions
C
C      DON'T CALL THIS ROUTINE.
C
C$ Literature_References
C
C      None.
C
C$ Author_and_Institution
C
C      N.J. Bachman    (JPL)
C
C$ Version
C
C-     SPICELIB Version 1.0.1, 10-MAR-1992 (WLT)
C
C         Comment section for permuted index source lines was added
C         following the header.
C
C-     SPICELIB Version 1.0.0, 31-JAN-1990 (NJB)
C
C-&

C$ Index_Entries
C
C     None.
C
C-&

C$ Revisions
C
C-     Beta Version 1.0.1, 08-FEB-1989 (NJB)
C
C         Warnings added to discourage use of this routine in
C         non-error-handling code.
C
C-&

C
C     Local Variables:
C

C
C     The SPICELIB status:
C

      LOGICAL               SVSTAT

      SAVE                  SVSTAT

C
C     Declaration of the entry point, FAILED:
C

      LOGICAL FAILED

C
C     Initial values:
C

      DATA    SVSTAT      / .FALSE. /

C
C     Executable Code:
C

      SVSTAT = STATUS

C
C     Give SETERR a value; the value does not have any
C     meaning, but it appears standard FORTRAN requires this.
C

      SETERR = .TRUE.

      RETURN

C$Procedure      FAILED ( Error Status Indicator )

       ENTRY FAILED ()

C$ Abstract
C
C      True if an error condition has been signalled via SIGERR.
C      FAILED is the SPICELIB status indicator.
C
C$ Disclaimer
C
C     THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE
C     CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S.
C     GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE
C     ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE
C     PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS"
C     TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY
C     WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A
C     PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC
C     SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE
C     SOFTWARE AND RELATED MATERIALS, HOWEVER USED.
C
C     IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA
C     BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT
C     LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND,
C     INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS,
C     REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE
C     REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY.
C
C     RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF
C     THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY
C     CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE
C     ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE.
C
C$ Required_Reading
C
C      ERROR
C
C$ Keywords
C
C      ERROR
C
C$ Declarations
C
C      None.
C
C$ Brief_I/O
C
C      VARIABLE  I/O  DESCRIPTION
C      --------  ---  --------------------------------------------------
C
C      The function takes the value .TRUE. if an error condition
C      was detected; it is .FALSE. otherwise.
C
C
C$ Detailed_Input
C
C      None.
C
C$ Detailed_Output
C
C      Please read the required reading file before reading this!
C
C      The value taken by FAILED indicates status.
C
C      The status value applies to the SPICELIB routines,
C      and to any other routines which call the status-setting
C      routine, SIGERR.
C
C      When FAILED has the value, .TRUE., an error condition
C      exists.   .FALSE. means "no error."
C
C      More specifically, when FAILED has the value .TRUE.,
C      some routine has indicated an error by calling the
C      SPICELIB routine, SIGERR.  All SPICELIB routines
C      which can detect errors do this.  Non-SPICELIB
C      routines may also reference SIGERR if desired.
C
C      When FAILED has the value .FALSE., either no routine
C      has yet signalled an error via SIGERR, or the status
C      has been reset using, what else, RESET.
C
C      FAILED is initialized to have the value, .FALSE.
C      This indicates a  "no error" status.
C
C      See "particulars" below for (slightly) more information.
C
C
C
C$ Parameters
C
C      None.
C
C$ Exceptions
C
C      None.
C
C      However, this routine is part of the SPICELIB error
C      handling mechanism.
C
C$ Files
C
C      None.
C
C$ Particulars
C
C      See the required reading file for details of error
C      processing.  However, here are some notes:
C
C      When any SPICELIB routine detects an error, the
C      status is set to indicate an error condition via
C      a call to SIGERR.   After SIGERR
C      returns, further calls to FAILED will return the
C      value, .TRUE., indicating an error condition.
C
C      Non-SPICELIB routines may also call SIGERR to indicate
C      an error condition; FAILED will reflect such calls
C      as well.
C
C      It is possible to re-set the error status to indicate
C      "no error" using the SPICELIB routine, RESET (see).
C
C      The effect on FAILED of resetting the status is
C      that FAILED will again return the value .FALSE.,
C      indicating "no error."
C
C      One of the main virtues of the SPICELIB error
C      handling mechanism is that you don't HAVE to test the
C      error status after every call to a SPICELIB routine.
C      If you set the error handling mode to 'RETURN', using
C      the routine, ERRACT, SPICELIB routines won't crash
C      when an error occurs; following the detection of the
C      error, each routine will return immediately upon entry.
C      Therefore, you call several SPICELIB routines in a
C      row, and just test status at the end of the sequence
C      of calls, if you wish.  See "examples" below.
C
C
C$ Examples
C
C      1.  Here's an example of a simple call to RDTEXT, followed
C          by a test of the status.
C
C
C      C
C      C     We read a line of text from file SPUD.DAT:
C      C
C
C            CALL RDTEXT ( 'SPUD.DAT', LINE, EOF )
C
C            IF ( FAILED() ) THEN
C
C      C        An error occurred during the read.
C
C               [respond to error here]
C
C            END IF
C
C
C      2.    Here's an example in which we don't want to
C            put the error test inside our loop.  We just
C            test the error status after the loop terminates.
C            We can do this because we (that is, you, the user)
C            have made the call,
C
C                   CALL ERRACT ( 'RETURN' )
C
C            prior to execution of the following code.  If an
C            error does occur, the remaining calls to RDTEXT
C            will have no effect.  Here's the example:
C
C      C
C      C     We read the first 5000 lines of a file, or until
C      C     EOF is reached, whichever comes first:
C      C
C      C     Note:  the "DO WHILE" construct is available in
C      C     VAX FORTRAN.
C      C
C
C            LCOUNT = 0
C            DO WHILE (  ( .NOT. EOF ) .AND. ( LCOUNT .LE. 5000 )  )
C
C               CALL RDTEXT ( 'SPUD.DAT', LINE(LCOUNT), EOF )
C
C               LCOUNT = LCOUNT + 1
C
C            END DO
C
C            IF ( FAILED() ) THEN
C      C
C      C        An error occurred during the read
C      C
C               [respond to error here]
C
C            END IF
C
C
C
C$ Restrictions
C
C      This routine automatically detects errors occurring in
C      the SPICELIB code.  To make this routine work
C      for your own routines, your routines must call SIGERR
C      to report errors.
C
C$ Literature_References
C
C      None.
C
C$ Author_and_Institution
C
C      N.J. Bachman    (JPL)
C      H.A. Neilan     (JPL)
C
C$ Version
C
C-     SPICELIB Version 1.0.1, 10-MAR-1992 (WLT)
C
C         Comment section for permuted index source lines was added
C         following the header.
C
C-     SPICELIB Version 1.0.0, 31-JAN-1990 (NJB)
C
C-&

C$ Index_Entries
C
C     error status indicator
C
C-&

C$ Revisions
C
C-     Beta Version 1.1.0, 18-DEC-1989 (HAN)
C
C         Empty parentheses added to the ENTRY statement in order to
C         comply with the ANSI Fortran 77 Standard.
C
C-&

C
C     Executable Code:
C

C
C     Grab saved status value:
C

      FAILED = SVSTAT

      END

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


text:
!==============================================================================!
  module Div_Mod
!------------------------------------------------------------------------------!
!   Variable definitions to be used in Divisor.  There are so few in fact      !
!   that there is very little reason for this module to exist at all.          !
!------------------------------------------------------------------------------!
  implicit none
!==============================================================================!

  ! Buffer send index and buffer receive index.
  ! Used for plotting dcomposed grids with links.
  integer, allocatable :: buf_send_ind(:), buf_recv_ind(:), buf_pos(:)

  end module

max_stars_repo_name: MassimoBorrelliPhysicist/T-Flows

max_issues_repo_name: MassimoBorrelliPhysicist/T-Flows

max_forks_repo_name: MassimoBorrelliPhysicist/T-Flows


module grid_mod
  use optical_properties_mod, only:nwl,nlayer
  implicit none
  save
  !contains
  integer, parameter :: nxg=101,nyg=101,nzg=101
  !integer, parameter :: nxg=200,nyg=200,nzg=200
  real*8 :: xface(-3:nxg+4),yface(-3:nyg+4),zface(-3:nzg+4)
  real*8 :: MASK(nxg,nyg,nzg,nlayer)
  real*8 :: rhokap(nxg,nyg,nzg)

  REAL*8 :: PL_SUM(nwl,NXG,NYG,NZG)

  real*8,parameter :: xmax=0.05,ymax=0.05,zmax=0.05
  !real*8,parameter :: xmax=1.d0,ymax=1.d0,zmax=1.d0
  real*8 :: grid_max(3)
  real*8 :: delta

end module grid_mod


text:
	SUBROUTINE RU_TTAA  ( report, lenr, itopwn, wnknot, iwndht,
     +			      ipoint, istnm, data, nlev, tdata, nlevt,
     +			      wxdata, nlevw, drpwnd, iret )
C************************************************************************
C* RU_TTAA								*
C*									*
C* This subroutine decodes TTAA reports.  These reports contain		*
C* surface, mandatory, tropopause and max wind data below 100 mb.       *
C* Surface/mandatory data are ordered PRES TEMP DWPT DRCT SPED HGHT.    *
C* Tropopause data are ordered PRES TEMP DWPT DRCT SPED.  Max wind data *
C* are ordered PRES DRCT SPED.                                          *
C*									*
C* RU_TTAA  ( REPORT, LENR, ITOPWN, WNKNOT, IWNDHT, IPOINT, ISTNM, 	*
C*            DATA, NLEV, TDATA, NLEVT, WXDATA, NLEVW, DRPWND, IRET )	*
C*									*
C* Input parameters:							*
C*	REPORT		CHAR*		TTAA report			*
C*	LENR		INTEGER		Length of station report	*
C*	ITOPWN		INTEGER		Pressure for last wind report	*
C*	WNKNOT		LOGICAL		Flag for wind in knots		*
C*	IWNDHT		INTEGER		Dropsonde wind cutoff value	*
C*									*
C* Input and output parameters:						*
C*	IPOINT		INTEGER		Pointer within report		*
C*	ISTNM		INTEGER		Station number			*
C*									*
C* Output parameters:							*
C*	DATA (6,NLEV)	REAL		Mandatory level data		*
C*	NLEV		INTEGER		Number of mandatory levels	*
C*	TDATA (5,NLEVT) REAL		Tropopause level data           *
C*	NLEVT		INTEGER		Number of tropopauses           *
C*	WXDATA(3,NLEVW)	REAL		Maximum wind level data         *
C*	NLEVW		INTEGER		Number of maximum wind levels   *
C*	DRPWND		LOGICAL		Flag to replace dropsonde wind	*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C*					 -2 = no data found		*
C*					 -4 = bad dropsonde id          *
C**									*
C* Log:									*
C* M. desJardins/GSFC	 6/86						*
C* D. Kidwell/NCEP	 2/01	Added tropopause and max wind data      *
C* D. Kidwell/NCEP	 3/01	Added RU_DRP1 call and argument istnm   *
C* D. Kidwell/NCEP	 2/05	Added argument drop to RU_MLVL call     *
C* m.gamazaychikov/SAIC	07/05	Added parms to CS and code for drpwnd	*
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
C*
	CHARACTER*(*)	report
	REAL		data (*), tdata (*), wxdata (*)
	LOGICAL		wnknot
C*
	LOGICAL		above, drop, drpwnd
        INCLUDE         'ERMISS.FNC'
C------------------------------------------------------------------------
	iret   = 0
	above  = .false.
	drpwnd = .false.
C
C*	Get the surface data.
C
	CALL RU_MSFC  ( report, lenr, wnknot, ipoint, data, ier1 )
C
C*	Check for a dropsonde.
C
	drop = .false.
        IF ( ( istnm .ge. 0 ) .and. ( istnm .le. 99 ) ) drop = .true.
C
C*	Read the mandatory level data.  Start the data after the
C*	surface data.
C
	CALL RU_MLVL  ( report, lenr, wnknot, itopwn, above, drop,
     + 			ipoint, data ( 7 ), nlev, ier2 )
C
C*	Check for neither mandatory nor surface data.
C
	IF  ( ( ier1 .ne. 0 ) .and. ( nlev .eq. 0 ) )  THEN
	    iret = -2
	  ELSE
	    nlev = nlev + 1
	END IF
C
C*	Read the tropopause data.
C
 	CALL RU_MTRP ( report, lenr, wnknot, above, ipoint, tdata,
     +		       nlevt, ier )
C
C*	Read the maximum wind data.
C
 	CALL RU_MMXW ( report, lenr, wnknot, above, ipoint, wxdata,
     +		       nlevw, ier )
C
C*      Check the station number.  If it is an hour value, this is a
C*      dropsonde report and number must be gotten from the 61616
C*      section.
C
        IF ( drop ) THEN
            CALL RU_DRP1 ( report, lenr, ipoint, istnm, lstwnd, ier )
	    IF ( iret .eq. 0 )  iret = ier
            IF ( iret .eq. 0 ) THEN
C
C*             Check if the last reported wind is less than or
C*             equal to the cutoff value
C
	       IF ( lstwnd .ne. IMISSD ) THEN
	          IF ( lstwnd .le. iwndht) THEN
	             IF ( ERMISS ( data ( 4 ) ) .and. .not.
     +                    ERMISS ( data ( 1 ) ) ) THEN
		        drpwnd = .true.
	             END IF
	          END IF
	       END IF
	    END IF
        END IF
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
!  SFM  04/06/2013  Code module added for WOlson
!
      SUBROUTINE RADTRAN(UMU, NLYR, TB, BTEMP, LYRTEMP, LYRHGT, KEXT,
     $                  SALB, ASYM, FISOT, EMIS, EBAR, LAMBERT, PRNT)
C
C     CHRIS KUMMEROW
C     INCLUDES ASYMPTOTIC EXPRESSIONS FOR TERM3, TERM4, AND TERM5 IN
C     THE LIMIT OF SMALL EFFECTIVE OPTICAL DEPTHS; BILL OLSON FEB, 1995.
C
!     WSO 04/07/2013  large layer optical depth numerics (double precision) on outgoing radiance
!     WSO 04/07/2013  delta-eddington option
!     WSO 04/07/2013  there are many code changes relative to radtran.f, here

      implicit none

      INTEGER *4 MAXLYR
      PARAMETER    ( MAXLYR = 80 )

      CHARACTER*1  POLN
      LOGICAL LAMBERT, PRNT(3)
      integer ier, imode, icode
      INTEGER NLYR, I, J, NANG, NN
      REAL  UMU, BTEMP, LYRTEMP(0:MAXLYR), LYRHGT(0:MAXLYR),
     $      KEXT(MAXLYR), SALB(MAXLYR), ASYM(MAXLYR)
      REAL  FISOT, EMIS, EBAR

      logical delta
      real *8 fdelta
      real *8 kext_delta(maxlyr)
      real *8 salb_delta(maxlyr)
      real *8 asym_delta(maxlyr)

      real *8 l(maxlyr)
      real *8 h(maxlyr)
      real *8 b0(maxlyr)
      real *8 b1(maxlyr)
      real *8 w(2*maxlyr, 2*maxlyr)
      real *8 btempd
      real *8 bb(2*maxlyr)
      real *8 dp(maxlyr)
      real *8 dm(maxlyr)
      real *8 z(0:maxlyr)
      real *8 mu
      real *8 nu
      real *8 x(2*maxlyr)
      real *8 fisotd
      real *8 emisd
      real *8 ebard

      real *8 xa
      real *8 xb
      real *8 xc
      real *8 xd
      real *8 ya
      real *8 yb
      real *8 yc
      real *8 dz
      real *8 xnu

      real *4 XIUP
      real *4 I_IN(maxlyr+1, 100)
      real *4 IOUT(0:maxlyr)
      real *4 TERM1
      real *4 TERM2
      real *4 TERM3
      real *4 TERM4
      real *4 TERM5
      real *4 TERM4prime
      real *4 TERM5prime
      real *4 tb

!delta-Eddington option
      delta = .false.
      if(.not. delta) then
        do i = 1, nlyr
          kext_delta(i) = dble(kext(i))
          salb_delta(i) = dble(salb(i))
          asym_delta(i) = dble(asym(i))
        end do
      else
        do i = 1, nlyr
          fdelta = dble(asym(i) * asym(i))
          kext_delta(i) = dble((1. - salb(i) * fdelta) * kext(i))
          salb_delta(i) = dble((1. - fdelta) * salb(i) /
     $     (1. - salb(i) * fdelta))
          asym_delta(i) = dble(asym(i) / (1. + asym(i)))
        end do
      end if

      fisotd = dble(fisot)
      emisd = dble(emis)
      ebard = dble(ebar)
      btempd = dble(btemp)

C
c      write(6,'(1x,"umu:  ",f7.2,"  nlyr: ",i5,"  tb: ",f7.2,
c     + "  btemp: ",f7.2,"  fisot: ",f7.2,
c     + /1x,"emis: ",f10.3,"  ebar: ",f10.3,"  lamb: ",l5,
c     + /1x,"  lyrtemp(0): ",f7.2,"  lyrhgt(0): ",f7.2)')
c     + umu,nlyr,tb,btemp,fisot,emis,ebar,lambert,lyrtemp(0),
c     + lyrhgt(0)
c      do j=nlyr,1,-1
c      write(6,'(1x,"lyrtemp: ",f7.2,"  lyrhgt: ",f7.2,
c     + "  kext: ",f10.4,"  salb: ",f10.4,"  asym: ",f10.4)')
c     + lyrtemp(j),lyrhgt(j),kext(j),salb(j),asym(j)
c      end do
C                                CALCULATE SOME CONSTANTS
      Z(0) = dble(LYRHGT(0))
      DO 20  J = 1,NLYR
        Z(J)  = dble(LYRHGT(J))
        B0(J) = dble(LYRTEMP(J-1))
        B1(J) = dble((LYRTEMP(J) - LYRTEMP(J-1))/
     $   (LYRHGT(J) - LYRHGT(J-1)))
        L(J) = dsqrt(3.*KEXT_delta(J)*KEXT_delta(J)*
     $         (1. - SALB_delta(J))*
     $         (1. - SALB_delta(J)*ASYM_delta(J)))
        H(J) = dble(1.5*KEXT_delta(J)*
     $   (1. - SALB_delta(J)*ASYM_delta(J)))
  20  CONTINUE

C     FILL IN THE MATRIX ELEMENTS WHICH FORM THE BOUNDARY CONDITIONS
C     AT THE TOP, BOTTOM AND LAYER INTERFACES OF THE CLOUD.  THERE ARE
C     TWO QUANTITIES, D+ "DP", AND D- "DM" AT EACH BOUNDARY, SO THE
C     MATRIX HAS DIMENSIONS  2*NLYR BY 2*NLYR.
C     ORDER OF D'S:  D+(1),D-(1),D+(2),D-(2), .... , D+(NLYR),D-(NLYR)

C     SET ALL MATRIX ELEMENTS TO ZERO
      DO 45  I = 1,2*NLYR
       DO 45  J = 1,2*NLYR
        W(I,J) = 0.d0
  45  CONTINUE

C     FILL IN THE NON-ZERO MATRIX ELEMENTS
      W(1,1)   = ((ebard - 2.d0)*L(1)/H(1)) + ebard
      W(1,2)   = ((2.d0 - ebard)*L(1)/H(1)) + ebard
      DO 50  I = 2,2*(NLYR-1),2
       W(I,I-1)   =
     $  (1.d0 - L(I/2)/H(I/2))*dexp(+L(I/2)*(Z(I/2)-Z(I/2-1)))
       W(I,I  )   =  (1.d0 + L(I/2)/H(I/2))*
     $  dexp(-L(I/2)*(Z(I/2)-Z(I/2-1)))
       W(I,I+1)   = -(1.d0 - L(I/2+1)/H(I/2+1))
       W(I,I+2)   = -(1.d0 + L(I/2+1)/H(I/2+1))

       W(I+1,I-1) =  (1.d0 + L(I/2)/H(I/2))*
     $  dexp(+L(I/2)*(Z(I/2)-Z(I/2-1)))
       W(I+1,I)   =  (1.d0 - L(I/2)/H(I/2))*
     $  dexp(-L(I/2)*(Z(I/2)-Z(I/2-1)))
       W(I+1,I+1) = -(1.d0 + L(I/2+1)/H(I/2+1))
       W(I+1,I+2) = -(1.d0 - L(I/2+1)/H(I/2+1))
c        write(*,*) (w(i,i+j),j=-2,2)
  50  CONTINUE
      W(2*NLYR,2*NLYR-1) =  (1.d0 + L(NLYR)/H(NLYR))*dexp(+L(NLYR)*
     $                       (Z(NLYR)-Z(NLYR-1)))
      W(2*NLYR,2*NLYR)   =  (1.d0 - L(NLYR)/H(NLYR))*dexp(-L(NLYR)*
     $                       (Z(NLYR)-Z(NLYR-1)))

C     FILL IN THE ROW OF CONSTANTS IN THE LINEAR EQUATIONS
      BB(1)    = ebard*btempd - ebard*B0(1) - (ebard - 2.d0)*B1(1)/H(1)
      DO 55  I = 2,2*(NLYR-1),2
       BB(I)   =  + B1(I/2)/H(I/2) - B1(I/2+1)/H(I/2+1)
       BB(I+1) =  - B1(I/2)/H(I/2) + B1(I/2+1)/H(I/2+1)
  55  CONTINUE
      BB(2*NLYR)  =  fisotd - B0(NLYR) - B1(NLYR)*(Z(NLYR)
     $               - Z(NLYR-1) + 1.d0/H(NLYR))
C
C
C     MATRIX INVERSION IN DONE IN SUBROUTINE LINPAK
c      write(*,*) bb
c      write(*,*) b1
c      stop
      imode=0
      call band(imode,2*nlyr,w,bb,icode)
!      call linsysa(w,bb,x,2*nlyr,2*maxlyr,ier)

!     CALL LINPAK(NLYR, W, BB, RCOND)
!      do i=1,2*nlyr
!         write(*,*) x(i), bb(i)
!      enddo
!      stop
!

      do i=1,2*nlyr
!         bb(i)=x(i)
      enddo
!      do i = 1, nlyr
!        write(*, '("i: ", i5, "  bbim1: ", e15.6,
!     $   "  bbi: ", e15.6)') i, bb(2*i-1), bb(2*i)
!      enddo
c     write(*,*) bb
c     stop
C
c     WRITE(*,*) ' '
c     IF ( PRNT(2) ) WRITE(4,*) ' '
      DO 60  I = 1,NLYR
       DP(I) = BB(2*I-1)
       DM(I) = BB(2*I)
c       WRITE(6,653) I, DP(I), I, DM(I)
c       IF ( PRNT(2) )  WRITE(4,653) I, DP(I), I, DM(I)
  60  CONTINUE
 653  FORMAT(10X,'D+ (',I2,') = ',E11.4,6X,'D-(',I2,') = ',E11.4)

C     AFTER D'S ARE KNOWN, CALCULATE SURFACE RADIANCE

      MU = dble(UMU)
      NU = dble(-UMU)
C     FOR THE FOLLOWING CALCULATIONS, REFER TO APPENDIX B OF THESIS
C     *********************************************************************

      IF ( LAMBERT ) THEN
C      CALCULATE THE DOWNWELLING FLUX THROUGH THE ATMOSPHERE AT 81 ANGLES
       NANG = 81
       DO 997  NN = 1,NANG
        XNU = -(2.d0*NN - 1.d0)/(NANG*2.d0)
        I_IN(NLYR+1,NN) = fisotd
C       LOOP THROUGH THE REMAINING LAYERS
        DO 100  J = NLYR,1,-1
C        CALCULATE RADIANCE FROM TOP OF LAYER "J"
         XA = B0(J) - 1.5d0*SALB_delta(J)*ASYM_delta(J)*XNU*B1(J)/H(J)
         XB = B1(J)
         XC = SALB_delta(J)*DP(J)*
     $    (1.d0 - 1.5d0*ASYM_delta(J)*XNU*L(J)/H(J))
         XD = SALB_delta(J)*DM(J)*
     $    (1.d0 + 1.5d0*ASYM_delta(J)*XNU*L(J)/H(J))
         YA = KEXT_delta(J)/XNU
         YB = YA + L(J)
         YC = YA - L(J)
         DZ = Z(J) - Z(J-1)

         TERM1 = I_IN(J+1,NN)*dexp(YA*DZ)
         TERM2 = XA*(1. - dexp(YA*DZ))
c     new tests
         if(dabs(ya*dz) .lt. 1.d-5) then
           term3=-xb*ya*dz*dz
         else
           TERM3 = XB/YA*(dexp(YA*DZ)*(1.d0 - YA*DZ) - 1.d0)
         end if
         if(dabs(yb*dz) .lt. 1.d-5) then
           term4=-xc*ya*dz
         else
           TERM4 = XC*YA/YB*(1.d0 - dexp(YB*DZ))
         end if
         if(dabs(yc*dz) .lt. 1.d-5) then
           term5=-xd*ya*dz
         else
           TERM5 = XD*YA/YC*(1.d0 - dexp(YC*DZ))
         end if
         I_IN(J,NN) = TERM1 + TERM2 + TERM3 + TERM4 + TERM5
!        write(6,'(1x,"lev: ",i5,"  ang: ",i5,"  term1: ",f10.4,
!     +            "  term2: ",f10.4,"  term3: ",f10.4,"  term4: ",f10.4,
!     +            "  term5: ",f10.4,"  tbdn: ",f7.2)')
!     +   j,nn,term1,term2,term3,term4,term5,i_in(j,nn)
  100   CONTINUE
  997  CONTINUE
C
C      CALCULATE THE TOTAL DOWNWELLING FLUX REACHING THE SURFACE
       XIUP = 0.
       DO 47  NN = 1,NANG
        XIUP = XIUP + I_IN(1,NN)*(1./NANG)*(2.*NN-1.)/(2.*NANG)
  47   CONTINUE
       XIUP = 2.*XIUP
c       write(6,'(1x,"xiup: ",f7.2)') xiup
C
      ELSE

C      CALCULATE THE DOWNWELLING FLUX THROUGH THE ATMOSPHERE AT ANGLE MU
       NN = 22                      ! THIS IS A DUMMY INDEX
       I_IN(NLYR+1,NN) = fisotd
C      LOOP THROUGH THE REMAINING LAYERS
       DO 110  J = NLYR,1,-1
C       CALCULATE RADIANCE FROM TOP OF LAYER "J"
        XA = B0(J) - 1.5d0*SALB_delta(J)*ASYM_delta(J)*NU*B1(J)/H(J)
        XB = B1(J)
        XC = SALB_delta(J)*DP(J)*
     $   (1.d0 - 1.5d0*ASYM_delta(J)*NU*L(J)/H(J))
        XD = SALB_delta(J)*DM(J)*
     $   (1.d0 + 1.5d0*ASYM_delta(J)*NU*L(J)/H(J))
        YA = KEXT_delta(J)/NU
        YB = YA + L(J)
        YC = YA - L(J)
        DZ = Z(J) - Z(J-1)

        TERM1 = I_IN(J+1,NN)*dexp(YA*DZ)
        TERM2 = XA*(1. - dexp(YA*DZ))
c     new tests
        if(dabs(ya*dz) .lt. 1.d-5) then
          term3=-xb*ya*dz*dz
        else
          TERM3 = XB/YA*(dexp(YA*DZ)*(1.d0 - YA*DZ) - 1.d0)
        end if
        if(dabs(yb*dz) .lt. 1.d-5) then
          term4=-xc*ya*dz
        else
          TERM4 = XC*YA/YB*(1. - dexp(YB*DZ))
        end if
        if(abs(yc*dz) .lt. 1.d-5) then
          term5=-xd*ya*dz
        else
          TERM5 = XD*YA/YC*(1.d0 - dexp(YC*DZ))
        end if
        I_IN(J,NN) = TERM1 + TERM2 + TERM3 + TERM4 + TERM5
!        write(*, '("down layer: ", i5, "  i_in: ", e15.6,
!     $   "  ya: ", e15.6, "  dz: ", e15.6, "  exp(ya*dz): ", e15.6,
!     $   "  yc: ", e15.6, "  exp(yc*dz): ", e15.6, "  term1: ", e15.6,
!     $   "  term2: ", e15.6, "  term3: ", e15.6,
!     $   "  term4: ", e15.6, "  term5: ", e15.6)')
!     $   J, I_IN(J+1,NN), YA, DZ, dexp(YA*DZ), YC, dexp(YC*DZ),
!     $   TERM1, TERM2, TERM3, TERM4, TERM5
  110  CONTINUE
       XIUP = I_IN(1,22)
C
      ENDIF

C
C
      IOUT(0) = emisd*btempd + (1.d0 - emisd)*XIUP
      DO 101  J = 1,NLYR
C      CALCULATE THE UPWELLING RADIANCES AT THE TOP OF EACH LAYER J
       XA = B0(J) - 1.5d0*SALB_delta(J)*ASYM_delta(J)*MU*B1(J)/H(J)
       XB = B1(J)
       XC = SALB_delta(J)*DP(J)*
     $  (1.d0 - 1.5d0*ASYM_delta(J)*MU*L(J)/H(J))
       XD = SALB_delta(J)*DM(J)*
     $  (1.d0 + 1.5d0*ASYM_delta(J)*MU*L(J)/H(J))
       YA = KEXT_delta(J)/MU
       YB = YA + L(J)
       YC = YA - L(J)
       DZ = Z(J) - Z(J-1)

       TERM1 = IOUT(J-1)*dexp(-YA*DZ)
       TERM2 = XA*(1.d0 - dexp(-YA*DZ))
c     new tests
       if(dabs(ya*dz) .lt. 1.d-5) then
         term3=0.
       else
         TERM3 = XB/YA*(YA*DZ - 1.d0 + dexp(-YA*DZ))
       end if
       if(dabs(yb*dz) .lt. 1.d-5) then
         TERM4=xc*ya*dz*dexp(-ya*dz)
       else
         TERM4 = XC*YA/YB*(dexp((YB-YA)*DZ) - dexp(-YA*DZ))
       end if
       if(dabs(yc*dz) .lt. 1.d-5) then
         TERM5=xd*ya*dz*dexp(-ya*dz)
         TERM5prime=xd*ya*dz*dexp(-ya*dz)
       else
!modification for large optical depth in layer
         TERM5 = XD*YA/YC*dexp(-YA*DZ)*(dexp(YC*DZ) - 1.d0)
         TERM5prime = XD*YA/YC*(dexp(-L(J)*DZ) - dexp(-YA*DZ))
       end if
!       IOUT(J) = TERM1 + TERM2 + TERM3 + TERM4 + TERM5
!        write(*, '("up   layer: ", i5, "  i_ou: ", e13.6,
!     $   "  ya: ", e13.6, "  exp: ", e13.6, "  L: ", e13.6,
!     $   "  yb: ", e13.6, "  exp: ", e13.6, "  term1: ", e13.6,
!     $   "  term2: ", e13.6, "  term3: ", e13.6,
!     $   "  term4: ", e13.6, "  term5: ", e13.6, "  term5prime: ", e13.6)')
!     $   J, IOUT(J-1), YA, dexp(-YA*DZ), L(J), YB, dexp(-L(J)*DZ),
!     $   TERM1, TERM2, TERM3, TERM4, TERM5, TERM5prime
       IOUT(J) = TERM1 + TERM2 + TERM3 + TERM4 + TERM5prime
  101 CONTINUE
C
C
c      WRITE(*,*) ' '
c      IF ( PRNT(2) )  WRITE(4,*) ' '
      DO 44  J = 0,NLYR
c       WRITE(*,77) NLYR+1-J, I_IN(NLYR+1-J,22), NLYR-J, IOUT(NLYR-J)
c       IF ( PRNT(2) ) WRITE(4,77) NLYR+1-J, I_IN(NLYR+1-J,22),
c     $                NLYR-J, IOUT(NLYR-J)
  44  CONTINUE
  77  FORMAT(10X,'Iin(',I2,') = ',F7.1,10X,'Iout(',I2,') = ',F7.1)

C
      TB = IOUT(NLYR)
!      write(*,*) tb
      RETURN
      END


text:
! Test code from Liao
!     if (N > 5 .or. N < 7) READ 1, N, N, N

!     READ 1, N
      READ 1, N
!     READ *, N
      READ *, N
!     READ (*,1) PI
      READ (*,1) PI

    1 FORMAT (I3)
      END

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
  30 Day minimum rental

  7x8 (smallest): $55/mo + $25 application fee + $20 key deposit

  They want you to pay the first and last months rents up front.

Bonus  This self storage storage facilty is freakishly close to campus (over by the Domes), a plus for any outgoing freshmen who need to leave some stuff in Davis for the summer.

   Stuff!  If you got too much and you want to keep it, you need one of the storage units in Davis!

Im pretty sure this place is (or at least used to be) called Russell Park Storage. I rented from them last summer and thats the name that they went by. Somebody want to double check? Users/AbbyLawson

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search