!!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  module condition
   implicit none
   real*8 :: obsvconst(0:5)
   real*8 :: elements(0:27)
   real*8 :: c1(0:39),c2(0:39),c3(0:39),c4(0:39),mid(0:39)
   real*8 :: duration, V,P, azi, alt, t,cover
   integer :: year,month,day,typepe
  end module condition
!!***********************************************************
  subroutine eclip(UT,SLT,longitude,latitude,height,newmag, &
                  & yearin,monthin,dayin,t1_SLT,t4_SLT,tm_SLT,mag)
   USE condition,only: obsvconst,elements,duration,V,P,azi,alt,t, &
                & c1,c2,c3,c4,mid,cover,year,month,day

       use spmd_utils,          only: masterproc
       use cam_logfile,       only : iulog

   implicit none
   real*8, parameter  :: PI =3.14151926
   real*8 :: tmp,mag,ratio,vv(1:5),pp(1:5),azimuth(1:5),altitude(1:5),timee(1:5)
   integer :: i,j, eclitype  !! 0 - none, 1 - partial, 2 - annular, 3 - total
   integer :: yearin,monthin,dayin
   real*8 :: UT,sUT,longitude,latitude,height,t1,t4,tm
   real*8 :: newmag
   real*8 :: t1_SLT,t4_SLT,tm_SLT,SLT,x,SLTTemp
! Dang, 2017
   real :: R,M,D,S1 
    x=longitude

!    if (masterproc) write(iulog,*) 'eclip: Top SLT ', SLT

!! Observer constants -
!! (0) North Latitude (degree)
!! (1) East Longitude (degree)
!! (2) Altitude (metres)
!! (3) West time zone (hours)
!! (4) rho dsin O'
!! (5) rho dcos O'
!! (6) index into the elements array for the eclipse in question
!!
! Dang, 2017
!   year=2017
!   month=8
!   day=21
!   year=2019
!   month=12
!   day=26
!   year=2020
!   month=12
!   day=14
!   year=2021
!   month=6
!   day=10
!   year=2021
!   month=12
!   day=4

   year = yearin
   month = monthin
   day = dayin

!! Note that correcting for refraction will involve creating a "virtual" altitude
!! for each contact, and hence a different value of rho and O' for each contact!
   obsvconst(0)= latitude   !! degree
   obsvconst(1)= - longitude;   !! degree
   obsvconst(2)=height*1000.0; !! meter
   obsvconst(3)=0.0;  
!!! Eclipse Elements------------------------------
!
! First line -
!  (0) Julian date of maximum eclipse
!  (1) t0 - the TDT hour at which t=0
!  (2) tmin - the lowest allowed value of t
!  (3) tmax - the highest allowed value of t
!  (4) dUTC - the difference between the civilian "GMT" timescale and TDT
!  (5) dT - the difference between UT and TDT
! Second line -
!  (6) X0, X1, X2, X3 - X elements
! Third line -
! (10) Y0, Y1, Y2, Y3 - Y elements
! Fourth line -
! (14) D0, D1, D2 - D elements
! Fifth line -
! (17) M0, M1, M2 - mu elements
! Sixth line -
! (20) L10, L11, L12 - L1 elements
! Seventh line -
! (23) L20, L21, L22 - L2 elements
! Eigth line -
! (26) tan f1
! (27) tan f2------------------------------------------------------------
! Dang, 2017
         if (masterproc) write(iulog,*) 'Calling findEclipse for UT', UT

      call findEclipse(year,month,day)

!         if (masterproc) write(iulog,*) 'Done with findEclipse '

      tmp=atan(0.99664719*tan(obsvconst(0)*Pi/180.0));
      obsvconst(4)=0.99664719*dsin(tmp)+(obsvconst(2)/6378140.0)*dsin(obsvconst(0)*Pi/180.0);
      obsvconst(5)=dcos(tmp)+(obsvconst(2)/6378140.0*dcos(obsvconst(0)*Pi/180.0));
!! initial the values 
      timee=0; altitude=0; azimuth=0; pp=0; vv=0; eclitype=0; mag=0; ratio=0; duration=0; cover=0;
	  t1=0; t4=0; tm=0
      call getall()
!!==============================================================================================

!! Eclipse circumstances
!!  (0) Event type (C1=-2, C2=-1, Mid=0, C3=1, C4=2)
!!  (1) t
!! -- time-only dependent circumstances (and their per-hour derivatives) follow --
!!  (2) x  (3) y  (4)  d (5) dsin d  (6) dcos d  (7) mu  (8) l1  (9) l2  (10) dx
!! (11) dy (12) dd  (13) dmu  (14) dl1  (15) dl2
!! -- time and location dependent circumstances follow --
!! (16) h  (17) dsin h  (18) dcos h  (19) xi  (20) eta  (21) zeta  (22) dxi  (23) deta  (24) u
!! (25) v  (26) a  (27) b  (28) l1' (29) l2' (30) n^2
!! -- observational circumstances follow --
!! (31) p  (32) alt  (33) q  (34) v  (35) azi
!! (36) m (mid eclipse only) or limb correction applied (where available!)
!! (37) magnitude (mid eclipse only)
!! (38) moon/sun (mid eclipse only)
!! (39) calculated local event type for a transparent earth (0 = none, 1 = partial, 2 = annular, 3 = total)

         if (masterproc .and. mid(39) > 0) write(iulog,*) 'eclip: mid(39)', mid(39)

	  if (mid(39)==0) then
	     newmag=0.0  !! there is no eclipse
		 goto 30
	  endif

!! displaymid--------------------------------------------------  
      call gettime(mid);  timee(5)=t; tm=t;
      call getalt(mid);   altitude(5)=alt;
      call getazi(mid);   azimuth(5)=azi;
      call getp(mid);     pp(5)=p;
      call getv(mid);     vv(5)=v;
!! Display the information about 1st and 4th contact--------------
      call gettime(c1); timee(1)=t; t1=t;
      call getalt(c1);  altitude(1)=alt;
      call getazi(c1);  azimuth(1)=azi;
      call getp(c1);    pp(1)=p;
      call getv(c1);    vv(1)=v;
      call gettime(c4); timee(4)=t; t4=t;
      call getalt(c4);  altitude(4)=alt;
      call getazi(c4);  azimuth(4)=azi;
      call getp(c4);    pp(4)=p;
      call getv(c4);    vv(4)=v;
      mag=floor(100000.0*mid(37)+0.5)/100000.0;
      ratio= floor(100000.0*mid(38)+0.5)/100000.0;
!	  t1_SLT=mod(t1+longitude/15.0+24.0,24.0);
!	  t4_SLT=mod(t4+longitude/15.0+24.0,24.0);
!	  tm_SLT=mod(tm+longitude/15.0+24.0,24.0);
          t1_SLT=mod(t1+x/15.0+24.0,24.0);
         t4_SLT=mod(t4+x/15.0+24.0,24.0);
          tm_SLT=mod(tm+x/15.0+24.0,24.0);

!         if (masterproc) write(iulog,*) 'In eclip after gets, UT, SLT, t1_SLT, t4_SLT, timee(5) ', &
!	                                       UT, SLT, t1_SLT, t4_SLT, timee(5)

	  if (SLT<t1_SLT .OR. SLT>t4_SLT) then
!         if (masterproc) write(iulog,*) 'In eclip setting newmag to zero for night '
	      newmag=0;
		  goto 30
	  endif

!! linear interpolation     
!	  if (SLT<tm_SLT) then
!		  newmag=mag*(SLT-t1_SLT)/(tm_SLT-t1_SLT);
!	  else
!		  newmag=mag*(SLT-t4_SLT)/(tm_SLT-t4_SLT);
!	  endif

! Calculate the actual ratio of the obscured diameter of the solar photosphere   
! Dang, 2017

      R=0.5
      M=mag   
      D=sqrt(4*M*R-M*M)
      if (UT<=timee(5)) then 
        d=(timee(5)-UT)/(timee(5)-timee(1))*D
      else
        d=(timee(5)-UT)/(timee(5)-timee(4))*D
      endif
      S1=sqrt(d*d+(2*R-M)*(2*R-M))
      newmag=2*R-S1
!         if (masterproc) write(iulog,*) 'eclip: daytime - calculated newmag ', newmag
!         if (masterproc) write(iulog,*) 'eclip: daytime - R, S1 ', R, S1
30    return
  End subroutine eclip

!!-----------------------------Calculate a limb correction------------------------------------
   subroutine limbcorrectionc2(p, q)

      integer ::  n
	  real*8 :: i,C2limb2003May(0:317),p,q
	  real*8, parameter  :: PI =3.14151926

! C2 limb corrections for the 2003 May 31 annular eclipse in seconds--------------
!
! The first 3 elements of the array tell us that the remaining data starts at contact angle
! 197.37 degrees in 0.4 degree increments for 315 data points, which means that the last
! element is for angle 322.97 degrees
!
! These limb corrections were calculated by Fred Espenak, NASA/GSFC---------------
     C2limb2003May = [197.37, 0.4, 315.0, &
 &  9.85, 9.35, 8.85, 8.33, 7.80, 7.26, 6.71, 6.14, 5.56, 5.22, 5.06, &
 &  4.89, 4.70, 4.51, 4.31, 4.10, 3.88, 3.82, 3.91, 4.00, 4.07, 4.14, &
 &  4.20, 4.26, 4.31, 4.35, 4.38, 4.41, 4.43, 4.44, 4.44, 4.44, 4.43, &
 &  4.42, 4.39, 4.36, 4.33, 4.28, 4.23, 4.17, 4.11, 4.04, 3.96, 3.87, &
 &  3.84, 3.85, 3.86, 3.86, 3.86, 3.85, 3.83, 3.81, 3.78, 3.74, 3.70, &
 &  3.65, 3.59, 3.53, 3.47, 3.39, 3.31, 3.22, 3.13, 3.03, 2.93, 2.81, &
 &  2.70, 2.57, 2.44, 2.31, 2.16, 2.01, 1.86, 1.70, 1.59, 1.64, 1.69, &
 &  1.72, 1.76, 1.78, 1.81, 1.87, 1.92, 1.97, 2.01, 2.05, 2.08, 2.11, &
 &  2.13, 2.15, 2.16, 2.23, 2.37, 2.51, 2.64, 2.77, 2.89, 3.01, 3.12, &
 &  3.23, 3.33, 3.43, 3.52, 3.61, 3.69, 3.77, 3.84, 3.91, 3.97, 4.03, &
 &  4.08, 4.13, 4.17, 4.21, 4.24, 4.27, 4.29, 4.31, 4.32, 4.33, 4.33, &
 &  4.33, 4.32, 4.31, 4.29, 4.27, 4.24, 4.20, 4.17, 4.12, 4.08, 4.02, &
 &  3.96, 3.90, 3.83, 3.76, 3.68, 3.60, 3.51, 3.42, 3.32, 3.22, 3.11, &
 &  3.00, 2.88, 2.76, 2.63, 2.53, 2.50, 2.46, 2.41, 2.36, 2.31, 2.25, &
 &  2.19, 2.12, 2.04, 1.96, 1.88, 1.79, 1.70, 1.60, 1.49, 1.39, 1.27, &
 &  1.15, 1.03, 0.90, 0.77, 0.63, 0.49, 0.34, 0.19, 0.17, 0.21, 0.24, &
 &  0.26, 0.29, 0.30, 0.31, 0.32, 0.32, 0.31, 0.30, 0.29, 0.26, 0.24, &
 &  0.21, 0.17, 0.13, 0.09, 0.05, 0.18, 0.31, 0.43, 0.55, 0.66, 0.77, &
 &  0.87, 0.97, 1.05, 1.14, 1.22, 1.29, 1.36, 1.42, 1.48, 1.53, 1.57, &
 &  1.61, 1.65, 1.68, 1.70, 1.72, 1.73, 1.73, 1.73, 1.73, 1.72, 1.70, &
 &  1.69, 1.75, 1.81, 1.87, 1.92, 1.96, 2.00, 2.03, 2.06, 2.08, 2.10, &
 &  2.11, 2.12, 2.12, 2.12, 2.13, 2.16, 2.18, 2.19, 2.20, 2.21, 2.25, &
 &  2.35, 2.44, 2.52, 2.59, 2.66, 2.71, 2.77, 2.81, 2.85, 2.89, 2.91, &
 &  2.93, 2.94, 2.95, 2.94, 3.07, 3.22, 3.37, 3.51, 3.64, 3.77, 3.88, &
 &  3.99, 4.09, 4.18, 4.27, 4.34, 4.41, 4.47, 4.52, 4.57, 4.60, 4.63, &
 &  4.65, 4.66, 4.66, 4.66, 4.64, 4.62, 4.59, 4.55, 4.50, 4.45, 4.39, &
 &  4.31, 4.23, 4.15, 4.05, 3.94, 3.83, 3.71, 3.58, 3.44, 3.29, 3.14, &
 &  2.97, 2.80, 2.62, 2.88, 3.47, 4.05, 4.62, 5.16, 5.70, 6.21, 6.71, &
 &  7.19, 7.66, 8.10, 8.54, 8.95, 9.35, 9.73 ];

      i = (p * 180 / PI) - C2limb2003May(0)
      do while (i >= 360.0)
	    i=i- 360.0
	  end do
      do while (i < 0.0) 
	    i=i+360.0
	  end do
      i = i / C2limb2003May(1)
      if (i >= (C2limb2003May(2) - 1)) then 
         q = 999.0;
		 goto 10
      end if
      n = floor(i)
      q = (((C2limb2003May(n+4) - C2limb2003May(n+3)) * (i - n)) + C2limb2003May(n+3))  
10   end subroutine limbcorrectionc2
!!------------------------------------------------------------------
   subroutine limbcorrectionc3(p, q)

      integer ::  n
	  real*8 :: i,C3limb2003May(0:342),p,q
	  real*8, parameter  :: PI =3.14151926

! C3 limb corrections for the 2003 May 31 annular eclipse in seconds
!
! The first 3 elements of the array tell us that the remaining data starts at contact angle
! 353.37 degrees in 0.4 degree increments for 340 data points, which means that the last
! element is for angle 128.97 degrees
!
! These limb corrections were calculated by Fred Espenak, NASA/GSFC

 C3limb2003May = [353.37, 0.4, 340.0, &  
 &  -8.67, -6.49, -6.14, -6.36, -6.62, -6.87, -7.09, -7.29, -7.47, -7.64, -7.78, &  
 &  -7.90, -8.01, -8.09, -8.15, -8.20, -8.22, -8.23, -8.22, -8.18, -8.13, -8.06, &  
 &  -7.97, -7.86, -7.74, -7.59, -7.43, -7.25, -7.05, -6.83, -6.60, -6.34, -6.07, &  
 &  -5.78, -5.48, -5.15, -4.81, -4.79, -4.82, -4.84, -4.85, -4.84, -4.83, -4.80, &  
 &  -4.76, -4.70, -4.64, -4.56, -4.48, -4.38, -4.26, -4.14, -4.01, -3.86, -3.81, &  
 &  -3.79, -3.78, -3.92, -4.05, -4.17, -4.29, -4.40, -4.50, -4.59, -4.67, -4.74, &  
 &  -4.81, -4.87, -4.92, -4.96, -5.00, -5.02, -5.04, -5.05, -5.05, -5.04, -5.03, &  
 &  -5.00, -4.97, -4.93, -4.89, -4.83, -4.77, -4.70, -4.62, -4.53, -4.44, -4.33, &  
 &  -4.22, -4.20, -4.26, -4.32, -4.37, -4.42, -4.45, -4.49, -4.51, -4.53, -4.54, &  
 &  -4.55, -4.54, -4.54, -4.52, -4.50, -4.47, -4.44, -4.39, -4.35, -4.29, -4.23, &  
 &  -4.16, -4.09, -4.01, -3.92, -3.82, -3.72, -3.62, -3.50, -3.38, -3.37, -3.40, &  
 &  -3.42, -3.44, -3.46, -3.47, -3.47, -3.46, -3.45, -3.44, -3.42, -3.39, -3.36, &  
 &  -3.32, -3.27, -3.22, -3.18, -3.12, -3.06, -3.00, -2.93, -2.91, -2.90, -2.87, &  
 &  -2.85, -2.81, -2.77, -2.73, -2.68, -2.62, -2.57, -2.52, -2.63, -2.75, -2.87, &  
 &  -2.99, -3.10, -3.21, -3.31, -3.40, -3.49, -3.58, -3.66, -3.73, -3.81, -3.87, &  
 &  -3.93, -3.99, -4.04, -4.08, -4.12, -4.16, -4.19, -4.21, -4.23, -4.25, -4.26, &  
 &  -4.26, -4.26, -4.25, -4.24, -4.23, -4.23, -4.23, -4.22, -4.21, -4.20, -4.18, &  
 &  -4.15, -4.12, -4.09, -4.05, -4.04, -4.09, -4.14, -4.19, -4.23, -4.26, -4.31, &  
 &  -4.37, -4.44, -4.49, -4.54, -4.59, -4.63, -4.67, -4.76, -4.88, -4.99, -5.11, &  
 &  -5.21, -5.31, -5.41, -5.50, -5.58, -5.66, -5.74, -5.81, -5.88, -5.94, -5.99, &  
 &  -6.04, -6.09, -6.13, -6.16, -6.20, -6.22, -6.24, -6.26, -6.27, -6.27, -6.28, &  
 &  -6.28, -6.27, -6.26, -6.24, -6.22, -6.19, -6.16, -6.12, -6.08, -6.03, -5.98, &  
 &  -5.92, -5.86, -5.79, -5.72, -5.64, -5.56, -5.47, -5.47, -5.45, -5.44, -5.41, &  
 &  -5.39, -5.35, -5.31, -5.27, -5.22, -5.17, -5.11, -5.05, -4.98, -4.90, -4.82, &  
 &  -4.74, -4.67, -4.60, -4.53, -4.46, -4.43, -4.61, -4.78, -4.95, -5.11, -5.27, &  
 &  -5.42, -5.57, -5.70, -5.84, -5.96, -6.08, -6.20, -6.30, -6.40, -6.50, -6.59, &  
 &  -6.67, -6.75, -6.82, -6.88, -6.96, -7.13, -7.29, -7.44, -7.58, -7.72, -7.86, &  
 &  -7.98, -8.10, -8.22, -8.32, -8.42, -8.52, -8.60, -8.68, -8.75, -8.82, -8.88, &  
 &  -8.93, -8.98, -9.02, -9.06, -9.08, -9.11, -9.12, -9.13, -9.13, -9.12, -9.11, &  
 &  -9.09, -9.06, -9.03, -8.98, -8.94, -8.90, -8.87, -8.83, -8.79, -8.74, -8.68, &  
 &  -8.62, -8.54, -8.47, -8.38, -8.29, -8.18, -8.09, -8.59, -9.08, -9.56 ];
      i = (p * 180 / PI) - C3limb2003May(0)
      do while (i >= 360.0)
	    i=i- 360.0
	  end do
      do while (i < 0.0) 
	    i=i+360.0
	  end do
      i = i / C3limb2003May(1)
      if (i >= (C3limb2003May(2) - 1)) then 
         q = 999.0; 
		 goto 10
      end if
      n = floor(i)
      q = (((C3limb2003May(n+4) - C3limb2003May(n+3)) * (i - n)) + C3limb2003May(n+3))  
 10  end subroutine limbcorrectionc3
!!==============================================================================================

!! Eclipse circumstances
!!  (0) Event type (C1=-2, C2=-1, Mid=0, C3=1, C4=2)
!!  (1) t
!! -- time-only dependent circumstances (and their per-hour derivatives) follow --
!!  (2) x  (3) y  (4)  d (5) dsin d  (6) dcos d  (7) mu  (8) l1  (9) l2  (10) dx
!! (11) dy (12) dd  (13) dmu  (14) dl1  (15) dl2
!! -- time and location dependent circumstances follow --
!! (16) h  (17) dsin h  (18) dcos h  (19) xi  (20) eta  (21) zeta  (22) dxi  (23) deta  (24) u
!! (25) v  (26) a  (27) b  (28) l1' (29) l2' (30) n^2
!! -- observational circumstances follow --
!! (31) p  (32) alt  (33) q  (34) v  (35) azi
!! (36) m (mid eclipse only) or limb correction applied (where available!)
!! (37) magnitude (mid eclipse only)
!! (38) moon/sun (mid eclipse only)
!! (39) calculated local event type for a transparent earth (0 = none, 1 = partial, 2 = annular, 3 = total)

!!------------------------------------------------------------------------
!! Populate the circumstances array with the time-only dependent circumstances (x, y, d, m, ...)
   subroutine timedependent(circumstances) 
     use condition, only: elements,typepe
     implicit none
	 real*8 circumstances(0:39)
	 real*8, parameter  :: PI =3.14151926
     real*8 :: t, ans

  t = circumstances(1)
  !! Calculate x
  ans = elements(9) * t + elements(8)
  ans = ans * t + elements(7)
  ans = ans * t + elements(6)
  circumstances(2) = ans
  !! Calculate dx
  ans = 3.0 * elements(9) * t + 2.0 * elements(8)
  ans = ans * t + elements(7)
  circumstances(10) = ans
  !! Calculate y
  ans = elements(13) * t + elements(12)
  ans = ans * t + elements(11)
  ans = ans * t + elements(10)
  circumstances(3) = ans
  !! Calculate dy
  ans = 3.0 * elements(13) * t + 2.0 * elements(12)
  ans = ans * t + elements(11)
  circumstances(11) = ans
  !! Calculate d
  ans = elements(16) * t + elements(15)
  ans = ans * t + elements(14)
  ans = ans * PI / 180.0
  circumstances(4) = ans
  !! dsin d and dcos d
  circumstances(5) = dsin(ans)
  circumstances(6) = dcos(ans)
  !! Calculate dd
  ans = 2.0 * elements(16) * t + elements(15)
  ans = ans * PI / 180.0
  circumstances(12) = ans
  !! Calculate m
  ans = elements(19) * t + elements(18)
  ans = ans * t + elements(17)
  if (ans >= 360.0) then 
    ans = ans - 360.0
  end if
  ans = ans * PI / 180.0
  circumstances(7) = ans
  !! Calculate dm
  ans = 2.0 * elements(19) * t + elements(18)
  ans = ans * PI / 180.0
  circumstances(13) = ans
  !! Calculate l1 and dl1
  typepe = circumstances(0)
  if ((typepe == -2) .OR. (typepe == 0) .OR. (typepe == 2)) then
    ans = elements(22) * t + elements(21)
    ans = ans * t + elements(20)
    circumstances(8) = ans
    circumstances(14) = 2.0 * elements(22) * t + elements(21)
  end if
  !! Calculate l2 and dl2
  if ((typepe == -1) .OR. (typepe == 0) .OR. (typepe == 1)) then
    ans = elements(25) * t + elements(24)
    ans = ans * t + elements(23)
    circumstances(9) = ans
    circumstances(15) = 2.0 * elements(25) * t + elements(24)
  end if
   end subroutine timedependent
!!====================================================================
!!------------------------------------------------------------------------
!! Populate the circumstances array with the time and location dependent circumstances
  subroutine timelocdependent(circumstances)
     use condition
	 implicit none
	 real*8, parameter  :: PI =3.14151926
	 real*8 :: circumstances(0:39)

  call timedependent(circumstances)
  !! Calculate h, dsin h, dcos h
  circumstances(16) = circumstances(7) - obsvconst(1)*PI/180.0 - (elements(5) / 13713.44)
  circumstances(17) = dsin(circumstances(16))
  circumstances(18) = dcos(circumstances(16))
  !! Calculate xi
  circumstances(19) = obsvconst(5) * circumstances(17)
  !! Calculate eta
  circumstances(20) = obsvconst(4) * circumstances(6) - obsvconst(5) * circumstances(18) * circumstances(5)
  !! Calculate zeta
  circumstances(21) = obsvconst(4) * circumstances(5) + obsvconst(5) * circumstances(18) * circumstances(6)
  !! Calculate dxi
  circumstances(22) = circumstances(13) * obsvconst(5) * circumstances(18)
  !! Calculate deta
  circumstances(23) = circumstances(13) * circumstances(19) * circumstances(5) - circumstances(21) * circumstances(12)
  !! Calculate u
  circumstances(24) = circumstances(2) - circumstances(19)
  !! Calculate v
  circumstances(25) = circumstances(3) - circumstances(20)
  !! Calculate a
  circumstances(26) = circumstances(10) - circumstances(22)
  !! Calculate b
  circumstances(27) = circumstances(11) - circumstances(23)
  !! Calculate l1'
  typepe = circumstances(0)
  if ((typepe == -2) .OR. (typepe == 0) .OR. (typepe == 2)) then
    circumstances(28) = circumstances(8) - circumstances(21) * elements(26)
  endif
  !! Calculate l2'
  if ((typepe == -1) .OR. (typepe == 0) .OR. (typepe == 1)) then
    circumstances(29) = circumstances(9) - circumstances(21) * elements(27)
  endif
  !! Calculate n^2
  circumstances(30) = circumstances(26) * circumstances(26) + circumstances(27) * circumstances(27)

 end subroutine timelocdependent
 !!==========================================================================================

!! Iterate on C1 or C4-----------------------------------------------------------------------
  subroutine c1c4iterate(circumstances) 
  use condition
  implicit none
  real*8 :: circumstances(0:39) 
  real*8 :: sign, tmp, n,tt
  integer :: iter

  call timelocdependent(circumstances)
  if (circumstances(0) < 0) then
    sign=-1.0
  else 
    sign=1.0
  end if
  tmp=1.0
  iter=0;
  do while (((tmp > 0.000001) .OR. (tmp < -0.000001)) .AND. (iter < 50)) 
    n = sqrt(circumstances(30))
    tmp = circumstances(26) * circumstances(25) - circumstances(24) * circumstances(27)
    tmp = tmp / n / circumstances(28)
	tt=1.0 - tmp * tmp;
	if(tt<0) tt=0;
    tmp = sign * sqrt(tt) * circumstances(28) / n
    tmp = (circumstances(24) * circumstances(26) + circumstances(25) * circumstances(27)) / circumstances(30) - tmp
    circumstances(1) = circumstances(1) - tmp
    call timelocdependent(circumstances)
    iter = iter +1
  end do
  end subroutine c1c4iterate
!!======================================================================================
!! Get C1 and C4 data--------------------------------------------------------------------
!!   Entry conditions -
!!   1. The mid array must be populated
!!   2. The magnitude at mid eclipse must be > 0.0
  subroutine getc1c4() 
  USE condition, only : c1,c4,mid
  real*8 :: tmp, n,tt
  n = sqrt(mid(30))
  tmp = mid(26) * mid(25) - mid(24) * mid(27)
  tmp = tmp / n / mid(28)
  tt=1.0 - tmp * tmp;
  if(tt<0) tt=0;
  tmp = sqrt(tt) * mid(28) / n
  c1(0) = -2;
  c4(0) = 2;
  c1(1) = mid(1) - tmp;
  c4(1) = mid(1) + tmp;
  call c1c4iterate(c1)
  call c1c4iterate(c4)
  end subroutine getc1c4

!!===============================================================

!! Iterate on C2 or C3--------------------------------------------------------
  subroutine c2c3iterate(circumstances) 
    Use condition
    real*8 :: sign,  tmp, n,circumstances(0:39),tt
	integer :: iter

  call timelocdependent(circumstances)
  if (circumstances(0) < 0) then
    sign=-1.0
   else 
    sign=1.0
  endif
  if (mid(29) < 0.0) then
    sign = -sign
  endif
  tmp=1.0
  iter=0
  do while (((tmp > 0.000001) .OR. (tmp < -0.000001)) .AND. (iter < 50)) 
    n = sqrt(circumstances(30))
    tmp = circumstances(26) * circumstances(25) - circumstances(24) * circumstances(27)
    tmp = tmp / n / circumstances(29)
    tt=1.0 - tmp * tmp;
	if(tt<0) tt=0;
    tmp = sign * sqrt(tt) * circumstances(29) / n
    tmp = (circumstances(24) * circumstances(26) + circumstances(25) * circumstances(27)) / circumstances(30) - tmp
    circumstances(1) = circumstances(1) - tmp
    call timelocdependent(circumstances)
    iter= iter+1
  end do
end subroutine c2c3iterate
!!===============================================================================
!!--------------------------------------------------------------------------------
!! Get C2 and C3 data
!!   Entry conditions -
!!   1. The mid array must be populated
!!   2. There mut be either a total or annular eclipse at the location!
  subroutine getc2c3() 
  use condition, only : c2,c3,mid
  real*8 :: tmp, n

  n = sqrt(mid(30))
  tmp = mid(26) * mid(25) - mid(24) * mid(27)
  tmp = tmp / n / mid(29)
  tt=1.0 - tmp * tmp;
  if(tt<0) tt=0;
  tmp = sqrt(tt) * mid(29) / n
  c2(0) = -1
  c3(0) = 1
  if (mid(29) < 0.0) then
    c2(1) = mid(1) + tmp
    c3(1) = mid(1) - tmp
   else 
    c2(1) = mid(1) - tmp
    c3(1) = mid(1) + tmp
  endif
  call c2c3iterate(c2)
  call c2c3iterate(c3)
  end subroutine getc2c3
!==============================================================================================
!----------------------------------------------------------------------------------------------
!! Get the observational circumstances
  subroutine observational(circumstances) 
  use condition, only : obsvconst,mid
  implicit none
  real*8, parameter  :: PI =3.14151926
  real*8 :: contacttype, dcoslat, dsinlat
  real*8 :: circumstances(0:39)
  !! We are looking at an "external" contact UNLESS this is a total eclipse AND we are looking at
  !! c2 or c3, in which case it is an INTERNAL contact! Note that if we are looking at mid eclipse,
  !! then we may not have determined the type of eclipse (mid[39]) just yet!
  if (circumstances(0) == 0) then
    contacttype = 1.0
   else if ((mid(39) == 3) .AND. ((circumstances(0) == -1) .OR. (circumstances(0) == 1))) then
      contacttype = -1.0
       else 
      contacttype = 1.0
   endif
  !! Calculate p
  circumstances(31) = atan2(contacttype*circumstances(24), contacttype*circumstances(25))
  !! Calculate alt
  dsinlat = dsin(obsvconst(0)*Pi/180.0)
  dcoslat = dcos(obsvconst(0)*Pi/180.0)
  circumstances(32) = dasin(circumstances(5) * dsinlat + circumstances(6) * dcoslat * circumstances(18))
  !! Calculate q
  circumstances(33) = dasin(dcoslat * circumstances(17) / dcos(circumstances(32)))
  if (circumstances(20) < 0.0) then
    circumstances(33) = PI - circumstances(33)
  endif
  !! Calculate v
  circumstances(34) = circumstances(31) - circumstances(33)
  !! Calculate azi
  circumstances(35) = atan2(-1.0*circumstances(17)*circumstances(6), circumstances(5)*dcoslat &
    &- circumstances(18)*dsinlat*circumstances(6))
  end subroutine observational
!!=======================================================================================
!!----------------------------------------------------------------------------------------
!! Calculate mid eclipse
   subroutine getmid() 
    use condition, only : mid
    real*8 :: tmp
    integer :: iter
    mid(0) = 0
    mid(1) = 0.0
    iter = 0
    tmp = 1.0
    call timelocdependent(mid)
    do while (((tmp > 0.000001) .OR. (tmp < -0.000001)) .AND. (iter < 50)) 
       tmp = (mid(24) * mid(26) + mid(25) * mid(27)) / mid(30)
       mid(1) = mid(1) - tmp
       iter = iter + 1
       call timelocdependent(mid)
    end do
  end subroutine getmid
!!====================================================================================
!! Populate the c1, c2, mid, c3 and c4 arrays
subroutine getall() 
  use condition, only: c1,c2,c3,c4,mid,year,month,day
  implicit none
  real*8 :: C2limb2003May,C3limb2003May
  call getmid()
  call observational(mid)
  !! Calculate m, magnitude and moon/sun
  mid(36) = sqrt(mid(24)*mid(24) + mid(25)*mid(25))
  mid(37) = (mid(28) - mid(36)) / (mid(28) + mid(29))
  mid(38) = (mid(28) - mid(29)) / (mid(28) + mid(29))
  if (mid(37) > 0.0) then
    call getc1c4()
    if ((mid(36) < mid(29)) .OR. (mid(36) < -mid(29))) then
      call getc2c3()
      if (mid(29) < 0.0) then
        mid(39) = 3 !! Total eclipse
       else 
        mid(39) = 2 !! Annular eclipse
      end if
      call observational(c2)
      call observational(c3)
      !! 2003 May 31 eclipse limb corrections -
      if (year==2003 .and. month==5 .and. day==31) then
        call limbcorrectionc2(c2(31), c2(36))
        call limbcorrectionc3(c3(31), c3(36))
        if (c2(36) < 990.0) then
          c2(1) = c2(1) + c2(36) / 3600.0
        endif
        if (c3(36) < 990.0) then
          c3(1) = c3(1) + c3(36) / 3600.0
        endif
      else 
        c2(36) = 999.9
        c3(36) = 999.9
      endif
    else 
      mid(39) = 1 !! Partial eclipse
    end if
    call observational(c1)
    call observational(c4)
 else 
    mid(39) = 0 !! No eclipse
 end if
end subroutine getall
!!===========================================================================================
!!----------------------Get the local time of an event---------------------------------------
  subroutine gettime(circumstances)
	 use condition, only: obsvconst,elements,t
     real*8 :: circumstances(0:39)
!! Calculate the local time. Add 0.05 seconds, as we will be rounding up to the nearest 0.1 sec
      t = circumstances(1) + elements(1) - obsvconst(3) - (elements(4) - 0.05) / 3600.0
      if (t < 0.0) then
       t = t + 24.0
      endif
      if (t >= 24.0) then
       t = t - 24.0
      endif
 end subroutine gettime
!!===========================================================================================
!!----------------------------- Get the altitude----------------------------------------------
subroutine getalt(circumstances) 
  use condition, only: alt
  implicit none
  real*8 :: circumstances(0:39)
  real*8, parameter  :: PI =3.14151926
  alt = circumstances(32) * 180.0 / PI
end subroutine getalt
!!==========================================================================================
!!------------------------------Get the azimuth--------------------------------------------
subroutine getazi(circumstances) 
  use condition, only: azi
  implicit none
  real*8 :: circumstances(0:39)
  real*8, parameter  :: PI =3.14151926

  azi = circumstances(35) * 180.0 / PI
  if (azi < 0.0) then
    azi = azi + 360.0
  endif
  if (azi >= 360.0) then
    azi = azi - 360.0
  endif
end subroutine getazi
!!============================================================================================
!!------------------------------Get the P-----------------------------------------------------
subroutine getp(circumstances) 
  use condition, only: P
  implicit none
  real*8 :: circumstances(0:39)
  real*8, parameter  :: PI =3.14151926
  p = circumstances(31) * 180.0 / PI
  if (p < 0.0) then
    p = p + 360.0
  endif
  if (p >= 360.0) then
    p = p - 360.0
  endif
end subroutine getp
!!===========================================================================================
!!---------------------------------------Get V-----------------------------------------------
subroutine getv(circumstances) 
  use condition, only: V
  implicit none
  real*8 :: circumstances(0:39)
  real*8, parameter  :: PI =3.14151926
  V = floor(120.5 - circumstances(34) * 60.0 / PI) / 10.0
  if (V > 13.0) then
    V = V - 12.0
  endif
  if (V > 13.0) then
    V = V - 12.0
  endif
  if (V < 1.0) then
    V = V + 12.0
  endif
end subroutine getV
!!============================================================================================
!!-----------------------Get duration------------------------------------------------------
subroutine getduration() 
  use condition, only : c2,c3,duration
  implicit none
  duration=c3(1)-c2(1);
  if (duration<0.0) then
    duration=duration+24.0
  else if (duration >= 24.0) then
    duration=duration-24.0
  endif
end subroutine getduration
!!========================================================================================

!!----------------------------- Get the coverage-------------------------------------------
subroutine getcoverage() 
  use condition, only: mid,cover
  implicit none
  real*8, parameter  :: PI =3.14151926
  real*8 :: a, b, c

  if (mid(37) <= 0.0) then
    cover=0
  end if
  if (mid(37) >= 1.0) then
    cover=1.0
  end if
  if (mid(39) == 2) then
    c = mid(38) * mid(38)
   else 
    c = dacos((mid(28)*mid(28) + mid(29)*mid(29) - 2.0*mid(36)*mid(36)) / (mid(28)*mid(28) - mid(29)*mid(29)))
    b = dacos((mid(28)*mid(29) + mid(36)*mid(36))/mid(36)/(mid(28)+mid(29)))
    a = PI - b - c
    c = ((mid(38)*mid(38)*a + b) - mid(38)*dsin(c))/PI
  end if
  cover = c+0.5/1000.0
end subroutine getcoverage
!--------------------------------------------------------------------------
  subroutine findEclipse(syear,smonth,sday)
     Use condition, only: elements

       use spmd_utils,          only: masterproc
       use cam_logfile,       only : iulog

     implicit none
	 integer :: iyear,imonth,iday,syear,smonth,sday,I
    iyear=1900;
!         if (masterproc) write(iulog,*) 'Opening Elements file '

     open(9,file='/glade/work/dangt/highres_tiegcm/tiegcm_res1.25_data/eclipse/Elements.dat',status='Old')
!         if (masterproc) write(iulog,*) 'Reading Elements file '

     do while(iyear<=syear)
!         do i=1,9
!            read(9,*)
!         end do
        read(9,*,end=10)iyear,imonth,iday,elements(0:27)
	    if (iyear==syear .and. imonth==smonth .and. iday==sday) then

!         if (masterproc) write(iulog,*) 'Year,month,day match eclipse '
!         if (masterproc) write(iulog,*) 'Year,month,day ',iyear,imonth,iday
!         if (masterproc) write(iulog,*) 'elements ', elements

!	       read(9,*) elements(0:5)
!               print *,'elements0-5',elements(0:5)
!               read(9,*) elements(6:9)
!	       read(9,*) elements(10:13)
!		   read(9,*) elements(14:16)
!		   read(9,*) elements(17:19)
!		   read(9,*) elements(20:22)
!		   read(9,*) elements(23:25)
!		   read(9,*) elements(26:27)
	       goto 20
        end if
      end do
10    close(9)
      write(*,*)'there is no eclipse at this time'
      return
!! Get the observer's geocentric position
20    close(9)
  end subroutine findEclipse
!--------------------------------------------------------------------------