SUBROUTINE COS_SZA(KSTART,KEND,KCOL,PGEMU,PGELAM,LDRADIATIONTIMESTEP,PMU0)

!**** *COS_SZA*   
!
! (C) Copyright 2015- ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
!
! In applying this licence, ECMWF does not waive the privileges and immunities
! granted to it by virtue of its status as an intergovernmental organisation
! nor does it submit to any jurisdiction.
!
!     Purpose.
!     --------
!        Compute the cosine of the solar zenith angle.  Note that this
!        is needed for three different things: (1) as input to the
!        radiation scheme in which it is used to compute the path
!        length of the direct solar beam through the atmosphere, (2)
!        every timestep to scale the solar fluxes by the incoming
!        solar radiation at top-of-atmosphere, and (3) to compute the
!        albedo of the ocean.  For (1) we ideally want an average
!        value for the duration of a radiation timestep while for (2)
!        we want an average value for the duration of a model
!        timestep.

!**   Interface.
!     ----------
!        *CALL* *COS_SZA(...)

!        Explicit arguments : 
!        ------------------
!            PGEMU - Sine of latitude
!            PGELAM - Geographic longitude in radians
!            LDRadiationTimestep - Is this for a radiation timestep?
!            PMU0 - Output cosine of solar zenith angle

!        Implicit arguments :
!        --------------------
!            YRRIP%RWSOVR, RWSOVRM - Solar time for model/radiation timesteps
!            RCODECM, RSIDECM - Sine/cosine of solar declination
!            YRERAD%LAverageSZA - Average solar zenith angle in time interval?
!            YRRIP%TSTEP - Model timestep in seconds
!            YRERAD%NRADFR - Radiation frequency in timesteps

!     Method.
!     -------
!        Compute cosine of the solar zenith angle, mu0, from lat, lon
!        and solar time using standard formula.  If
!        YRERAD%LAverageSZA=FALSE then this is done at a single time,
!        which is assumed to be the mid-point of either the model or
!        the radiation timestep.  If YRERAD%LAverageSZA=TRUE then we
!        compute the average over the model timestep exactly by first
!        computing sunrise/sunset times. For radiation timesteps, mu0
!        is to be used to compute the path length of the direct solar
!        beam through the atmosphere, and the fluxes are subsequently
!        weighted by mu0.  Therefore night-time values are not used,
!        so we average mu0 only when the sun is above the horizon.

!     Externals.
!     ----------

!     Reference.
!     ----------
!        ECMWF Research Department documentation of the IFS
!
!        See also: Zhou, L., M. Zhang, Q. Bao, and Y. Liu (2015), On
!        the incident solar radiation in CMIP5
!        models. Geophys. Res. Lett., 42, 1930–1935. doi:
!        10.1002/2015GL063239.

!     Author.
!     -------
!      Robin Hogan, ECMWF, May 2015

!     Modifications:
!     --------------

USE PARKIND1 , ONLY : JPIM, JPRB
USE YOMHOOK  , ONLY : LHOOK, DR_HOOK
USE YOMCST   , ONLY : RPI, RDAY
USE YOMRIP   , ONLY : YRRIP
USE YOERIP   , ONLY : YRERIP
USE YOERAD   , ONLY : YRERAD
USE YOMLUN   , ONLY : NULOUT

!     ------------------------------------------------------------------

IMPLICIT NONE

INTEGER(KIND=JPIM),INTENT(IN) :: KSTART      ! Start column to process
INTEGER(KIND=JPIM),INTENT(IN) :: KEND        ! Last column to process
INTEGER(KIND=JPIM),INTENT(IN) :: KCOL        ! Number of columns in arrays
REAL(KIND=JPRB),   INTENT(IN) :: PGEMU(KCOL) ! Sine of latitude
REAL(KIND=JPRB),   INTENT(IN) :: PGELAM(KCOL)! Longitude in radians
LOGICAL, INTENT(IN) :: LDRADIATIONTIMESTEP   ! Is this for a radiation timestep?
REAL(KIND=JPRB),  INTENT(OUT) :: PMU0(KCOL)  ! Cosine of solar zenith angle

! Solar time at the start and end of the time interval
REAL(KIND=JPRB) :: ZSOLARTIMESTART, ZSOLARTIMEEND

! The time of half a model/radiation timestep, in radians
REAL(KIND=JPRB) :: ZHALFTIMESTEP

! For efficiency we precompute sin(solar declination)*sin(latitude)
REAL(KIND=JPRB) :: ZSINDECSINLAT(KSTART:KEND)
!...and cos(solar declination)*cos(latitude)
REAL(KIND=JPRB) :: ZCOSDECCOSLAT(KSTART:KEND)
! ...and cosine of latitude
REAL(KIND=JPRB) :: ZCOSLAT(KSTART:KEND)

! Tangent of solar declination
REAL(KIND=JPRB) :: ZTANDEC

! Hour angles (=local solar time in radians plus pi)
REAL(KIND=JPRB) :: ZHOURANGLESTART, ZHOURANGLEEND
REAL(KIND=JPRB) :: ZHOURANGLESUNSET, ZCOSHOURANGLESUNSET

INTEGER(KIND=JPIM) :: JCOL        ! Column index

REAL(KIND=JPRB) :: ZHOOK_HANDLE

!     ------------------------------------------------------------------
IF (LHOOK) CALL DR_HOOK('COS_SZA',0,ZHOOK_HANDLE)

! An average solar zenith angle can only be computed if the solar time
! is centred on the time interval
IF (YRERAD%LAVERAGESZA .AND. .NOT. YRERAD%LCENTREDTIMESZA) THEN
  WRITE(NULOUT,*) 'ERROR IN COS_SZA: LAverageSZA=TRUE but LCentredTimeSZA=FALSE'
  CALL ABOR1('COS_SZA: ABOR1 CALLED')
ENDIF

DO JCOL = KSTART,KEND
  ZCOSLAT(JCOL) = SQRT(1.0_JPRB - PGEMU(JCOL)**2)
ENDDO

IF (LDRADIATIONTIMESTEP) THEN
  ! Compute the effective cosine of solar zenith angle for a radiation
  ! timestep

  ! Precompute quantities that may be used more than once
  DO JCOL = KSTART,KEND
    ZSINDECSINLAT(JCOL) = YRERIP%RSIDECM * PGEMU(JCOL)
    ZCOSDECCOSLAT(JCOL) = YRERIP%RCODECM * ZCOSLAT(JCOL)
  ENDDO

  IF (.NOT. YRERAD%LAVERAGESZA) THEN
    ! Original method: compute the value at the centre of the
    ! radiation timestep (assuming that LCentredTimeSZA=TRUE - see
    ! updtim.F90)
    DO JCOL = KSTART,KEND
      ! It would be more efficient to do it like this...
      ! PMU0(JCOL)=MAX(0.0_JPRB, ZSinDecSinLat(JCOL) &
      !      & - ZCosDecCosLat(JCOL) * COS(YRERIP%RWSOVRM + PGELAM(JCOL)))
      ! ...but for bit reproducibility with previous cycle we do it
      ! like this:
      PMU0(JCOL) = MAX(0.0_JPRB, ZSINDECSINLAT(JCOL) &
           & - YRERIP%RCODECM*COS(YRERIP%RWSOVRM)*ZCOSLAT(JCOL)*COS(PGELAM(JCOL)) &
           & + YRERIP%RCODECM*SIN(YRERIP%RWSOVRM)*ZCOSLAT(JCOL)*SIN(PGELAM(JCOL)))
    ENDDO

  ELSE
    ! Compute the average MU0 for the period of the radiation
    ! timestep, excluding times when the sun is below the horizon

    ! First compute the sine and cosine of the times of the start and
    ! end of the radiation timestep
    ZHALFTIMESTEP = YRRIP%TSTEP * REAL(YRERAD%NRADFR) * RPI / RDAY
    ZSOLARTIMESTART = YRERIP%RWSOVRM - ZHALFTIMESTEP
    ZSOLARTIMEEND   = YRERIP%RWSOVRM + ZHALFTIMESTEP

    ! Compute tangent of solar declination, with check in case someone
    ! simulates a planet completely tipped over
    ZTANDEC = YRERIP%RSIDECM / MAX(YRERIP%RCODECM, 1.0E-12)

    DO JCOL = KSTART,KEND
      ! Sunrise equation: cos(hour angle at sunset) =
      ! -tan(declination)*tan(latitude)
      ZCOSHOURANGLESUNSET = -ZTANDEC * PGEMU(JCOL) &
           &              / MAX(ZCOSLAT(JCOL), 1.0E-12)
      IF (ZCOSHOURANGLESUNSET > 1.0) THEN
        ! Perpetual darkness
        PMU0(JCOL) = 0.0_JPRB
      ELSE
        ! Compute hour angle at start and end of time interval,
        ! ensuring that the hour angle of the centre of the time
        ! window is in the range -PI to +PI (equivalent to ensuring
        ! that local solar time = solar time + longitude is in the
        ! range 0 to 2PI)
        IF (YRERIP%RWSOVRM + PGELAM(JCOL) < 2.0_JPRB*RPI) THEN
          ZHOURANGLESTART = ZSOLARTIMESTART + PGELAM(JCOL) - RPI
          ZHOURANGLEEND   = ZSOLARTIMEEND   + PGELAM(JCOL) - RPI 
        ELSE
          ZHOURANGLESTART = ZSOLARTIMESTART + PGELAM(JCOL) - 3.0_JPRB*RPI
          ZHOURANGLEEND   = ZSOLARTIMEEND   + PGELAM(JCOL) - 3.0_JPRB*RPI
        ENDIF

        IF (ZCOSHOURANGLESUNSET >= -1.0) THEN
          ! Not perpetual daylight or perpetual darkness, so we need
          ! to check for sunrise or sunset lying within the time
          ! interval
          ZHOURANGLESUNSET = ACOS(ZCOSHOURANGLESUNSET)
          IF (ZHOURANGLEEND <= -ZHOURANGLESUNSET &
               & .OR. ZHOURANGLESTART >= ZHOURANGLESUNSET) THEN
            ! The time interval is either completely before sunrise or
            ! completely after sunset
            PMU0(JCOL) = 0.0_JPRB
            CYCLE
          ENDIF

          ! Bound the start and end hour angles by sunrise and sunset
          ZHOURANGLESTART = MAX(-ZHOURANGLESUNSET, &
               &                MIN(ZHOURANGLESTART, ZHOURANGLESUNSET))
          ZHOURANGLEEND   = MAX(-ZHOURANGLESUNSET, &
               &                MIN(ZHOURANGLEEND,   ZHOURANGLESUNSET))
        ENDIF

        IF (ZHOURANGLEEND - ZHOURANGLESTART > 1.0E-8) THEN
          ! Compute average MU0 in the interval ZHourAngleStart to
          ! ZHourAngleEnd
          PMU0(JCOL) = ZSINDECSINLAT(JCOL) &
               & + (ZCOSDECCOSLAT(JCOL) &
               &    * (SIN(ZHOURANGLEEND) - SIN(ZHOURANGLESTART))) &
               & / (ZHOURANGLEEND - ZHOURANGLESTART)

          ! Just in case...
          IF (PMU0(JCOL) < 0.0_JPRB) THEN
            PMU0(JCOL) = 0.0_JPRB
          ENDIF
        ELSE
          ! Too close to sunrise/sunset for a reliable calculation
          PMU0(JCOL) = 0.0_JPRB
        ENDIF

      ENDIF
    ENDDO
  ENDIF

ELSE
  ! Compute the cosine of solar zenith angle for a model timestep

  ! Precompute quantities that may be used more than once
  DO JCOL = KSTART,KEND
    ZSINDECSINLAT(JCOL) = YRRIP%RSIDEC * PGEMU(JCOL)
    ZCOSDECCOSLAT(JCOL) = YRRIP%RCODEC * ZCOSLAT(JCOL)
  ENDDO

  IF (.NOT. YRERAD%LAVERAGESZA) THEN
    ! Original method: compute the value at the centre of the
    ! model timestep
    DO JCOL = KSTART,KEND
      ! It would be more efficient to do it like this...
      ! PMU0(JCOL) = MAX(0.0_JPRB, ZSinDecSinLat(JCOL)        &
      !      & - ZCosDecCosLat(JCOL)*COS(YRRIP%RWSOVR + PGELAM(JCOL)))
      ! ...but for bit reproducibility with previous cycle we do it
      ! like this:
      PMU0(JCOL) = MAX(0.0_JPRB, ZSINDECSINLAT(JCOL) &
           & - YRRIP%RCODEC*COS(YRRIP%RWSOVR)*ZCOSLAT(JCOL)*COS(PGELAM(JCOL)) &
           & + YRRIP%RCODEC*SIN(YRRIP%RWSOVR)*ZCOSLAT(JCOL)*SIN(PGELAM(JCOL)))
    ENDDO

  ELSE
    ! Compute the average MU0 for the period of the model timestep

    ! First compute the sine and cosine of the times of the start and
    ! end of the model timestep
    ZHALFTIMESTEP   = YRRIP%TSTEP * RPI / RDAY
    ZSOLARTIMESTART = YRRIP%RWSOVR - ZHALFTIMESTEP
    ZSOLARTIMEEND   = YRRIP%RWSOVR + ZHALFTIMESTEP

    ! Compute tangent of solar declination, with check in case someone
    ! simulates a planet completely tipped over
    ZTANDEC = YRRIP%RSIDEC / MAX(YRRIP%RCODEC, 1.0E-12)

    DO JCOL = KSTART,KEND
      ! Sunrise equation: cos(hour angle at sunset) =
      ! -tan(declination)*tan(latitude)
      ZCOSHOURANGLESUNSET = -ZTANDEC * PGEMU(JCOL) &
           &              / MAX(ZCOSLAT(JCOL), 1.0E-12)
      IF (ZCOSHOURANGLESUNSET > 1.0) THEN
        ! Perpetual darkness
        PMU0(JCOL) = 0.0_JPRB
      ELSE
        ! Compute hour angle at start and end of time interval,
        ! ensuring that the hour angle of the centre of the time
        ! window is in the range -PI to +PI (equivalent to ensuring
        ! that local solar time = solar time + longitude is in the
        ! range 0 to 2PI)
        IF (YRRIP%RWSOVR + PGELAM(JCOL) < 2.0_JPRB*RPI) THEN
          ZHOURANGLESTART = ZSOLARTIMESTART + PGELAM(JCOL) - RPI
          ZHOURANGLEEND   = ZSOLARTIMEEND   + PGELAM(JCOL) - RPI 
        ELSE
          ZHOURANGLESTART = ZSOLARTIMESTART + PGELAM(JCOL) - 3.0_JPRB*RPI
          ZHOURANGLEEND   = ZSOLARTIMEEND   + PGELAM(JCOL) - 3.0_JPRB*RPI
        ENDIF

        IF (ZCOSHOURANGLESUNSET >= -1.0) THEN
          ! Not perpetual daylight or perpetual darkness, so we need
          ! to check for sunrise or sunset lying within the time
          ! interval
          ZHOURANGLESUNSET = ACOS(ZCOSHOURANGLESUNSET)
          IF (ZHOURANGLEEND <= -ZHOURANGLESUNSET &
               & .OR. ZHOURANGLESTART >= ZHOURANGLESUNSET) THEN
            ! The time interval is either completely before sunrise or
            ! completely after sunset
            PMU0(JCOL) = 0.0_JPRB
            CYCLE
          ENDIF

          ! Bound the start and end hour angles by sunrise and sunset
          ZHOURANGLESTART = MAX(-ZHOURANGLESUNSET, &
               &                MIN(ZHOURANGLESTART, ZHOURANGLESUNSET))
          ZHOURANGLEEND   = MAX(-ZHOURANGLESUNSET, &
               &                MIN(ZHOURANGLEEND,   ZHOURANGLESUNSET))
        ENDIF

        ! Compute average MU0 in the model timestep, although the
        ! numerator considers only the time from ZHourAngleStart to
        ! ZHourAngleEnd that the sun is above the horizon
        PMU0(JCOL) = (ZSINDECSINLAT(JCOL) * (ZHOURANGLEEND-ZHOURANGLESTART)   &
           & + ZCOSDECCOSLAT(JCOL)*(SIN(ZHOURANGLEEND)-SIN(ZHOURANGLESTART))) &
           & / (2.0_JPRB * ZHALFTIMESTEP)

        ! This shouldn't ever result in negative values, but just in
        ! case
        IF (PMU0(JCOL) < 0.0_JPRB) THEN
          PMU0(JCOL) = 0.0_JPRB
        ENDIF

      ENDIF
    ENDDO
  ENDIF

ENDIF


!     ------------------------------------------------------------------
IF (LHOOK) CALL DR_HOOK('COS_SZA',1,ZHOOK_HANDLE)
END SUBROUTINE COS_SZA
