!----------------------------------------------------------------------
! This routine calculates sub-daily variability
! C. Huntingford (April 2001) - based on earlier version by P. Cox
!----------------------------------------------------------------------
  SUBROUTINE DAY_CALC(                                            &
    POINTSM,STEP_DAY,DAY_MON,SW_L,PRECIP_L,T_L,DTEMP_DAY_L,LW_L,  &
    PSTAR_L,WIND_L,RH15M_L,SW_SD,T_SD,LW_SD,CONV_RAIN_SD,         &
    LS_RAIN_SD,LS_SNOW_SD,PSTAR_SD,WIND_SD,QHUM_SD,MONTH,IDAY,    &
    LAT,LONG,SEC_DAY,NSDMAX,SEED_RAIN                             &
  )

    USE c_pi, ONLY : pi

    IMPLICIT NONE

    INTEGER ::                                                    &
      DAY_MON,                                                    &
              !IN "Number of days in a month" (day)
      POINTSM,                                                    &
              !IN Maximum number of points in grid.
      STEP_DAY,                                                   &
              !IN WORK Calculated number of timesteps per day
      SEC_DAY,                                                    &
              !IN Number of seconds in day (sec)
      NSDMAX,                                                     &
              !IN Maximum possible number of sub-daily timesteps.
      SEED_RAIN(4),                                               &
              !IN Seeding numbers required to disaggregate
              !   the rainfall
      N_TALLY !WORK Counting up number of precipitation periods. 

!-----------------------------------------------------------------------
! Single day arrays of incoming variables 
!-----------------------------------------------------------------------
    REAL ::                                                       &
      SW_L(POINTSM),                                              &
              ! IN Daily values for downward shortwave
              !    radiation (W/m2)
      PRECIP_L(POINTSM),                                          &
              ! In Daily values of rain+snow (mm/day)
      T_L(POINTSM),                                               &
              ! IN Daily values of temperature (K)
      DTEMP_DAY_L(POINTSM),                                       &
              ! IN Daily values of diurnal temperature range (K)
      LW_L(POINTSM),                                              &
              ! IN Daily values of surface longwave
              !    radiation (W/m**2). 
      PSTAR_L(POINTSM),                                           &
              ! IN Daily values of pressure at level 1 (Pa).
      WIND_L(POINTSM),                                            &
              ! IN Daily values of wind speed (m/s)
      QHUM_L(POINTSM),                                            &
              ! Humidity deficit calculated from RH15M (kg/kg) 
      RH15M_L(POINTSM)
              ! Humidity (%)

!-----------------------------------------------------------------------
! Single day (global) arrays of variables above, but for up to hourly
! periods. NOTE: DTEMP_DAY_SD does not exist as T_SD combines T_L 
! and DTEMP_DAY_L 
!-----------------------------------------------------------------------

    REAL ::                                                       &
      SW_SD(POINTSM,NSDMAX),                                      &
              ! OUT Sub-daily values for downward
              !     shortwave radiation (W/m2)
      T_SD(POINTSM,NSDMAX),                                       &
              ! OUT Sub-daily values of temperature (K)
      LW_SD(POINTSM,NSDMAX),                                      &
              ! OUT Sub-daily values of surface 
              !     longwave radiation (W/m**2). 
      PSTAR_SD(POINTSM,NSDMAX),                                   &
              ! OUT Sub-daily values of pressure at level 1 (Pa).
      WIND_SD(POINTSM,NSDMAX),                                    &
              ! OUT Sub-daily values of wind speed (m/s)
      QHUM_SD(POINTSM,NSDMAX),                                    &
              ! OUT Sub_daily humidity deficit calculated
              !     from RH15M (kg/kg) 
      CONV_RAIN_SD(POINTSM,NSDMAX),                               &
              ! OUT Sub-daily convective rain (mm/day)
      LS_RAIN_SD(POINTSM,NSDMAX),                                 &
              ! OUT Sub-daily large scale rain (mm/day)
      LS_SNOW_SD(POINTSM,NSDMAX)
              ! OUT Sub-daily large scale snow (mm/day)

    INTEGER ::                                                    &
      N_EVENT(POINTSM,NSDMAX),                                    &
              ! 1: if rains/snows during timestep period
              ! 0: otherwise
      N_EVENT_LOCAL(NSDMAX)
              ! As N_EVENT, but for each gridpoint

    REAL ::                                                       &
      PREC_LOC(NSDMAX)
              ! Temporary value of rainfall for each gridbox.

    REAL ::                                                       &
      T_SD_LOCAL(POINTSM),                                        &
              ! WORK Intermediate calculation of temperature (K)
      PSTAR_SD_LOCAL(POINTSM),                                    &
              ! WORK Intermediate calculation of pressure (Pa)
      QS_SD_LOCAL(POINTSM)
              ! WORK Saturated humidity deficit associated with
              !      T_SD_LOCAL and PSTAR_SD_LOCAL
 
    INTEGER ::                                                    &
      ISTEP,                                                      &
              ! WORK Loop over sub-daily periods
      MONTH,                                                      &
              ! IN Month of interest
      IDAY,                                                       &
              ! IN Day number since beginning of month
      DAYNUMBER,                                                  &
              ! WORK Daynumber since beginning of year
      L       ! WORK Loop over land points

    REAL ::                                                       &
      SUN(POINTSM,NSDMAX),                                        &
              ! WORK Normalised solar radiation
      TIME_MAX(POINTSM),                                          &
              ! WORK GMT at which temperature is maximum (hours)
      LAT(POINTSM),                                               &
              ! IN Latitude (degrees)
      LONG(POINTSM),                                              &
              ! IN Longitude (degrees)
      TIME_DAY,                                                   &
              ! WORK Time of day (hours)
      TIMESTEP,                                                   &
              ! WORK Timestep (seconds)
      RANDOM_NUM_SD
              ! WORK Random number associated with rai

    REAL ::                                                       &
      TEMP_CONV,                                                  &
              ! WORK Temperature above which rainfall is convective (K)
      TEMP_SNOW
              ! WORK Temperature below which snow occurs (K)
    PARAMETER(TEMP_CONV = 293.15) 
    PARAMETER(TEMP_SNOW = 275.15) 

    REAL ::                                                       &
      INIT_HOUR_CONV_RAIN,                                        &
              !WORK Start of convective rain event (hour)
      INIT_HOUR_LS_RAIN,                                          &
              !WORK Start of large scale rain event (hour)
      INIT_HOUR_LS_SNOW,                                          &
              !WORK Start of large scale snow event (hour)
      DUR_CONV_RAIN,                                              &
              !WORK Start of convective rain event (hour)
      DUR_LS_RAIN,                                                &
              !WORK Duration of large scale rain event (hour)
      DUR_LS_SNOW,                                                &
              !WORK Start of large scale snow event (hour)
      END_HOUR_CONV_RAIN,                                         &
              !WORK End of convective rain event (hour)
      END_HOUR_LS_RAIN,                                           &
              !WORK End of large scale rain event (hour)
      END_HOUR_LS_SNOW,                                           &
              !WORK End of large scale snow event (hour)
      HOUREVENT,                                                  &
              !WORK Local variable giving hours during diurnal
              !     period for checking if precip. event occurs
      MAX_PRECIP_RATE,                                            &
              !WORK Maximum allowed precip. rate allowed within
              !     each sub-daily timestep (mm/day). (This only
              !     applies when STEP_DAY >= 2)
      PERIOD_LEN
              !WORK Length of period (hr)


!-----------------------------------------------------------------------
! Calculate the maximum precipitation rate. It is noted that 58 mm/day
! over 8 timesteps, and where all fell within a single 3 hour period
! caused numerical issues for MOSES. This corresponded to a rate of
! 464 mm/day during the 3-hour period. Hence, place a limit of 350 
! mm/day.
!-----------------------------------------------------------------------
    PARAMETER(MAX_PRECIP_RATE = 350.0)

!-----------------------------------------------------------------------
! First check whether sub-daily calculations are required. 
!-----------------------------------------------------------------------

!      DUR_CONV_RAIN = 2.0
!      DUR_LS_RAIN = 5.0
!      DUR_LS_SNOW = 5.0

    DUR_CONV_RAIN = 6.0
    DUR_LS_RAIN = 1.0
    DUR_LS_SNOW = 1.0

    PERIOD_LEN = 24.0 / FLOAT(STEP_DAY)

!-----------------------------------------------------------------------
! First check whether sub-daily calculations are required. 
!-----------------------------------------------------------------------
    IF(STEP_DAY >= 2) THEN

!-----------------------------------------------------------------------
! Ensure that the durations are at least as long as a time period for
! the model to prevent solution "falling through gaps"
!-----------------------------------------------------------------------
      IF(DUR_CONV_RAIN <= PERIOD_LEN)                             &
        DUR_CONV_RAIN = PERIOD_LEN+1.0E-6
      IF(DUR_LS_RAIN <= PERIOD_LEN)                               &
        DUR_LS_RAIN = PERIOD_LEN+1.0E-6 
      IF(DUR_LS_SNOW <= PERIOD_LEN)                               &
        DUR_LS_SNOW = PERIOD_LEN+1.0E-6

      TIMESTEP = FLOAT(SEC_DAY) / FLOAT(STEP_DAY)

!-----------------------------------------------------------------------
! Calculate the diurnal cycle in the SW radiation
!-----------------------------------------------------------------------
      DAYNUMBER = INT((MONTH-1.0) * REAL(DAY_MON)) + IDAY
      CALL SUNNY(                                                 &
        DAYNUMBER,STEP_DAY,POINTSM,1990,LAT,LONG,SUN,TIME_MAX     &
      )

!-----------------------------------------------------------------------
! Loop over timesteps                                               
!-----------------------------------------------------------------------
      DO ISTEP=1,STEP_DAY                                             

!-----------------------------------------------------------------------
! Calculate timestep values of the driving data.                       
!-----------------------------------------------------------------------
        TIME_DAY = (REAL(ISTEP) - 0.5) * TIMESTEP

        DO L=1,POINTSM
          T_SD(L,ISTEP) = T_L(L) + 0.5 * DTEMP_DAY_L(L) *         &
                   COS(2*PI*(TIME_DAY-3600.0*TIME_MAX(L))/SEC_DAY)
          LW_SD(L,ISTEP) = LW_L(L)                                &
                         * (4.0 * T_SD(L,ISTEP) / T_L(L) - 3.0)
          SW_SD(L,ISTEP) = SW_L(L) * SUN(L,ISTEP)
        ENDDO

!-----------------------------------------------------------------------
! Calculate timestep values of the driving data that is not split up
! into diurnal behaviour.                       
!-----------------------------------------------------------------------
        DO L=1,POINTSM
          PSTAR_SD(L,ISTEP) = PSTAR_L(L)
          WIND_SD(L,ISTEP)  = WIND_L(L)
        ENDDO	  

!-----------------------------------------------------------------------
! Check that humidity value is not greater than QSAT (but otherwise, QHU
! is not split up into diurnal behaviour).                       
!-----------------------------------------------------------------------
        DO L=1,POINTSM
          PSTAR_SD_LOCAL(L) = PSTAR_SD(L,ISTEP)
          T_SD_LOCAL(L)     = T_SD(L,ISTEP)
        ENDDO

        CALL QSAT(QS_SD_LOCAL,T_SD_LOCAL,PSTAR_SD_LOCAL,POINTSM)
        
        DO L=1,POINTSM
          QHUM_SD(L,ISTEP) = 0.01 * RH15M_L(L) * QS_SD_LOCAL(L) 
        ENDDO

      ENDDO   ! End of timestep loop within the individual days.

!-----------------------------------------------------------------------
! Calculate daily rainfall disaggregation
!-----------------------------------------------------------------------
      DO L=1,POINTSM 
 
!-----------------------------------------------------------------------
! Precipitation is split into four components,
! these being large scale rain, convective rain, large scale snow,
! convective snow. Call random number generator for different
! durations.
!-----------------------------------------------------------------------
        CALL RNDM(RANDOM_NUM_SD,SEED_RAIN)

!-----------------------------------------------------------------------
! Calculate type of precipitation. The decision is based purely up
! mean daily temperature, T_L. The cutoffs are:
!
! Convective scale rain (duration CONV_RAIN_DUR): T_L > 20.0oC
! Large scale rain (duration LS_RAIN_DUR) : 20.0oC > T_L > 2oC 
! Large scale snow (duration LS_SNOW_DUR) : T_L < 2oC
! Convective snow - IGNORED
!
!-----------------------------------------------------------------------
! Initialise arrays
        DO ISTEP=1,STEP_DAY       
          CONV_RAIN_SD(L,ISTEP) = 0.0 
          LS_RAIN_SD(L,ISTEP)   = 0.0 
          LS_SNOW_SD(L,ISTEP)   = 0.0 
          N_EVENT(L,ISTEP)      = 0
          N_EVENT_LOCAL(ISTEP)  = 0
          PREC_LOC(ISTEP)       = 0.0
        ENDDO

! Calculate rainfall disaggregation. 

! Start with convective rain 
! (temperatures based upon mean daily temperature)

! First check if warm enough for convective rain
        IF(T_L(L) >= TEMP_CONV) THEN

          INIT_HOUR_CONV_RAIN = RANDOM_NUM_SD*(24.0-DUR_CONV_RAIN) 
          END_HOUR_CONV_RAIN = INIT_HOUR_CONV_RAIN+DUR_CONV_RAIN

          N_TALLY = 0

          DO ISTEP=1,STEP_DAY
            HOUREVENT = (REAL(ISTEP)-0.5)*PERIOD_LEN
            IF(HOUREVENT >= INIT_HOUR_CONV_RAIN .AND.             &
               HOUREVENT < END_HOUR_CONV_RAIN) THEN
              N_EVENT(L,ISTEP) = 1
              N_TALLY = N_TALLY + 1
	    ENDIF
          ENDDO

          DO ISTEP=1,STEP_DAY
            IF(N_EVENT(L,ISTEP).EQ.1) THEN      !Rains on this day
	      CONV_RAIN_SD(L,ISTEP) =                             &
                        (REAL(STEP_DAY)/REAL(N_TALLY))*PRECIP_L(L)
              PREC_LOC(ISTEP) = CONV_RAIN_SD(L,ISTEP)
              N_EVENT_LOCAL(ISTEP) = N_EVENT(L,ISTEP) 
            ENDIF
          ENDDO

! Check that no convective rain periods 
! exceed MAX_PRECIP_RATE, or if so, 
! then redistribute. The variable that is redistributed is local
! variable PREC_LOC - CONV_RAIN_SD is then set to this after the
! call to REDIS.

          CALL REDIS(                                             &
            NSDMAX,STEP_DAY,MAX_PRECIP_RATE,PREC_LOC,             &
            N_EVENT_LOCAL,N_TALLY                                 &
          )

          DO ISTEP=1,STEP_DAY
            CONV_RAIN_SD(L,ISTEP) = PREC_LOC(ISTEP)
          ENDDO

! Now look at large scale rainfall components
        ELSE IF(T_L(L) < TEMP_CONV .AND. T_L(L) >= TEMP_SNOW) THEN

          INIT_HOUR_LS_RAIN = RANDOM_NUM_SD*(24.0-DUR_LS_RAIN) 
          END_HOUR_LS_RAIN = INIT_HOUR_LS_RAIN+DUR_LS_RAIN 

          N_TALLY = 0

          DO ISTEP=1,STEP_DAY
            HOUREVENT = (REAL(ISTEP)-0.5)*PERIOD_LEN
            IF(HOUREVENT >= INIT_HOUR_LS_RAIN .AND.               &
               HOUREVENT < END_HOUR_LS_RAIN) THEN
              N_EVENT(L,ISTEP) = 1
              N_TALLY = N_TALLY + 1
            ENDIF
          ENDDO

          DO ISTEP=1,STEP_DAY
            IF(N_EVENT(L,ISTEP) == 1) THEN      !Rains on this day
              LS_RAIN_SD(L,ISTEP) =                               &
                        (REAL(STEP_DAY)/REAL(N_TALLY))*PRECIP_L(L) 
              PREC_LOC(ISTEP) = LS_RAIN_SD(L,ISTEP)
              N_EVENT_LOCAL(ISTEP) = N_EVENT(L,ISTEP) 
            ENDIF
          ENDDO

! Check that no large scale rain periods 
! exceed MAX_PRECIP_RATE, or if so, 
! then redistribute.

          CALL REDIS(                                             &
            NSDMAX,STEP_DAY,MAX_PRECIP_RATE,PREC_LOC,             &
            N_EVENT_LOCAL,N_TALLY                                 &
          )

          DO ISTEP=1,STEP_DAY
            LS_RAIN_SD(L,ISTEP) = PREC_LOC(ISTEP)
          ENDDO

! Now look at large scale snow components
        ELSE          

          INIT_HOUR_LS_SNOW = RANDOM_NUM_SD*(24.0-DUR_LS_SNOW) 
          END_HOUR_LS_SNOW = INIT_HOUR_LS_SNOW+DUR_LS_SNOW 

          N_TALLY = 0

          DO ISTEP=1,STEP_DAY
            HOUREVENT = (REAL(ISTEP)-0.5)*PERIOD_LEN
            IF(HOUREVENT >= INIT_HOUR_LS_SNOW .AND.               &
               HOUREVENT < END_HOUR_LS_SNOW) THEN
              N_EVENT(L,ISTEP) = 1
              N_TALLY = N_TALLY + 1
            ENDIF
          ENDDO

          DO ISTEP=1,STEP_DAY
            IF(N_EVENT(L,ISTEP) == 1) THEN       !Rains on this da
              LS_SNOW_SD(L,ISTEP) =                               &
                        (REAL(STEP_DAY)/REAL(N_TALLY))*PRECIP_L(L) 
              PREC_LOC(ISTEP) = LS_SNOW_SD(L,ISTEP)
              N_EVENT_LOCAL(ISTEP) = N_EVENT(L,ISTEP) 
            ENDIF
          ENDDO

! Check that no large scale snow periods exceed MAX_PRECIP_RATE, 
! or if so, then redistribute.
          CALL REDIS(                                             &
            NSDMAX,STEP_DAY,MAX_PRECIP_RATE,PREC_LOC,             &
            N_EVENT_LOCAL,N_TALLY                                 &
          )

          DO ISTEP=1,STEP_DAY
            LS_SNOW_SD(L,ISTEP) = PREC_LOC(ISTEP)
          ENDDO

        ENDIF

      ENDDO        ! End of large loop over different land points.
!                    in calculation of different rainfall behaviours

    ELSE           ! Now case where no subdaily variation (TSTEP=1)

      DO L=1,POINTSM
        PSTAR_SD_LOCAL(L) = PSTAR_SD(L,ISTEP)
        T_SD_LOCAL(L) = T_SD(L,ISTEP)
      ENDDO

      CALL QSAT(QS_SD_LOCAL,T_SD_LOCAL,PSTAR_SD_LOCAL,POINTSM)
 
      DO L = 1,POINTSM
        SW_SD(L,1) = SW_L(L)
        T_SD(L,1) = T_L(L)
        LW_SD(L,1) = LW_L(L)
        PSTAR_SD(L,1) = PSTAR_L(L)
        WIND_SD(L,1) = WIND_L(L)
        QHUM_SD(L,1) = 0.01*RH15M_L(L)*QS_SD_LOCAL(L) 

        IF(T_L(L) >= TEMP_CONV) THEN
          CONV_RAIN_SD(L,1) = PRECIP_L(L)
        ELSE IF(T_L(L) < TEMP_CONV .AND. T_L(L) >= TEMP_SNOW) THEN  
          LS_RAIN_SD(L,1) = PRECIP_L(L)
        ELSE
          LS_SNOW_SD(L,1) = PRECIP_L(L)
        ENDIF

      ENDDO

    ENDIF !End of loop to chose whether sub-daily is required

    RETURN

  END SUBROUTINE DAY_CALC
