!  Unified model deck SOLANG, containing only routine SOLANG.         
!    This is part of logical component P233, performing the           
!  calculations of the earth's orbit described in the second page of  
!  the "Calculation of incoming insolation" section of UMDP 23, i.e.  
!  from the sin of the solar  declination, the position of each point 
!  and the time limits it calculates how much sunlight, if any, it    
!  receives.                                                          
!    Written in FORTRAN 77 with the addition of "!" comments and      
!  underscores in variable names.                                     
!    Written to comply with 12/9/89 version of UMDP 4 (meteorological 
!  standard).                                                         
!    Author:    William Ingram  21/3/89                               
!                      Reviewer: Clive Wilson Winter 1989/90          
!  First version.                                                     


!*********************************************************************  
! A01_2A selected for SCM use equivalent to frozen version              
! physics 1/3/93 J.Lean                                                 
!*********************************************************************  
  SUBROUTINE SOLANG(SINDEC,T,DT,SINLAT,LONGIT,K,LIT,COSZ)

    USE c_pi

    IMPLICIT NONE

    INTEGER, INTENT(IN) ::                                        &
      K     ! Number of points

    REAL, INTENT(IN) ::                                           &
      SINDEC,                                                     &
            ! Sin(solar declination)
      T,DT,                                                       &
            ! Start time (GMT) & timestep
      SINLAT(K),LONGIT(K)
            ! sin(latitude) & longitude of each point

    REAL, INTENT(OUT) ::                                          &
      LIT(K),                                                     &
            ! Sunlit fraction of the timestep
      COSZ(K)
            ! Mean cos(solar zenith angle) during the sunlit
            ! fraction

! This routine has no dynamically allocated work areas.  It calls the  
! intrinsic functions SQRT, ACOS & SIN, but no user functions or       
! subroutines.  The only structure is a loop over all the points to be 
! dealt with, with IF blocks nested inside to cover the various        
! possibilities.                                                       

    INTEGER :: J  ! Loop counter over points

    REAL ::                                                       &
      TWOPI,                                                      &
            ! 2*pi
      S2R,                                                        &
            ! Seconds-to-radians converter
      SINSIN,COSCOS,                                              &
            ! Products of the sines and of the cosines
            ! of solar declination and of latitude.
      HLD,COSHLD,                                                 &
            ! Half-length of the day in radians (equal to the
            ! hour-angle of sunset, and minus the hour-angle of
            ! sunrise) & its cosine.
      HAT,                                                        &
            ! Local hour angle at the start time.
      OMEGAB,OMEGAE,OMEGA1,OMEGA2,OMEGAS,                         &
            ! Beginning and end of the timestep and of the
            ! period over which cosz is integrated, and sunset
            ! - all measured in radians after local sunrise,
            ! not from local noon as the true hour angle is.
      DIFSIN,DIFTIM,                                              &
            ! A difference-of-sines intermediate value and the
            ! corresponding time period
      TRAD, DTRAD
            ! These are the start-time and length of the
            ! timestep (T & DT) converted to radians after
            ! midday GMT, or equivalently, hour angle of the
            ! sun on the Greenwich meridian.

!-----------------------------------------------------------------
! Set up parameter values
!-----------------------------------------------------------------
    PARAMETER( TWOPI = 2. * PI, S2R = PI / 43200.)



    TRAD = T * S2R - PI
    DTRAD = DT * S2R

!DIR$ IVDEP
    DO J = 1,K                          ! Loop over points

! Logically unnecessary statement without which the CRAY compiler will not vectorize this code
      HLD = 0.0

      SINSIN = SINDEC * SINLAT(J)
      COSCOS = SQRT( (1.-SINDEC**2) * (1.-SINLAT(J)**2) )
      COSHLD = SINSIN / COSCOS
      IF(COSHLD < -1.0) THEN             ! Perpetual night
        LIT(J)  = 0.0
        COSZ(J) = 0.0
      ELSE
        HAT = LONGIT(J) + TRAD           ! (3.2.2)                
        IF(COSHLD > 1.0) THEN            !   Perpetual day - hour 
          OMEGA1 = HAT                   ! angles for (3.2.3) are 
          OMEGA2 = HAT + DTRAD           ! start & end of timestep
        ELSE                                !   At this latitude some
! points are sunlit, some not.  Different ones need different treatment.
          HLD = ACOS(-COSHLD)               ! (3.2.4)
! The logic seems simplest if one takes all "times" - actually hour     
! angles - relative to sunrise (or sunset), but they must be kept in the
! range 0 to 2pi for the tests on their orders to work.                 
          OMEGAB = HAT + HLD                                         
          IF(OMEGAB < 0.0)    OMEGAB = OMEGAB + TWOPI
          IF(OMEGAB >= TWOPI) OMEGAB = OMEGAB - TWOPI
          IF(OMEGAB >= TWOPI) OMEGAB = OMEGAB - TWOPI
!            !  Line repeated - otherwise could have failure if         
!            !  longitudes W are > pi rather than < 0.                  
          OMEGAE = OMEGAB + DTRAD                                    
          IF(OMEGAE > TWOPI)  OMEGAE = OMEGAE - TWOPI
          OMEGAS = 2. * HLD

! Now that the start-time, end-time and sunset are set in terms of hour 
! angle, can set the two hour-angles for (3.2.3).  The simple cases are 
! start-to-end-of-timestep, start-to-sunset, sunrise-to-end and sunrise-
! -to-sunset, but two other cases exist and need special treatment.     
          IF(OMEGAB <= OMEGAS .OR. OMEGAB < OMEGAE) THEN
            OMEGA1 = OMEGAB - HLD
          ELSE
            OMEGA1 = - HLD
          ENDIF

          IF(OMEGAE <= OMEGAS) THEN
            OMEGA2 = OMEGAE - HLD
          ELSE
            OMEGA2 = OMEGAS - HLD
          ENDIF

!  Put in an arbitrary marker for the case when the sun does not rise   
!  during the timestep (though it is up elsewhere at this latitude).    
!  (Cannot set COSZ & LIT within the ELSE ( COSHLD < 1 ) block          
!  because 3.2.3 is done outside this block.)                           
          IF(OMEGAE > OMEGAB .AND. OMEGAB > OMEGAS)               &
            OMEGA2=OMEGA1

        ENDIF           ! This finishes the ELSE (perpetual day) block

        DIFSIN = SIN(OMEGA2) - SIN(OMEGA1)             ! Begin (3.2.3)
        DIFTIM = OMEGA2 - OMEGA1

! Next, deal with the case where the sun sets and then rises again      
! within the timestep.  There the integration has actually been done    
! backwards over the night, and the resulting negative DIFSIN and DIFTIM
! must be combined with positive values representing the whole of the   
! timestep to get the right answer, which combines contributions from   
! the two separate daylit periods.  A simple analytic expression for the
! total sun throughout the day is used.  (This could of course be used  
! alone at points where the sun rises and then sets within the timestep)
        IF(DIFTIM < 0.0) THEN
          DIFSIN = DIFSIN + 2. * SQRT(1.-COSHLD**2)                   
          DIFTIM = DIFTIM + 2. * HLD                                  
        ENDIF                                                         

        IF (DIFTIM == 0.0) THEN
! Pick up the arbitrary marker for night points at a partly-lit latitude
          COSZ(J) = 0.0
          LIT(J)  = 0.0
        ELSE
          COSZ(J) = DIFSIN*COSCOS/DIFTIM + SINSIN     ! (3.2.3)      
          LIT(J) = DIFTIM / DTRAD
        ENDIF
      ENDIF            ! This finishes the ELSE (perpetual night) block
    ENDDO

    RETURN                                                            

 END SUBROUTINE SOLANG
