suews_ctrl_input.f95

Go to the documentation of this file.

!Reading one line of meteorological forcing data in.
!Latest change:
!  Feb 2012, LJ:  Input fluxes qh and qe changed _obs as well as qn1_obs ending
!  Oct 2014, LJ:  Variables changed only be used in this part of code and these are passed to calling
!                 function in MetArray.
!  Jan 2015, HCW: Precip_hr, wuh and LAI_hr changed for generic timesteps
!  Jan 2016, LJ:  Removal of tabs
!  Feb 2017, HCW: Added file unit as argument so MetRead can be used for original met forcing file too
!  May 2018, TS: imporved the ability in downscaling datetime by introducing a datetime library
! To Do:
!       - Check observed SM calculation
!---------------------------------------------------------------------------------------------------
SUBROUTINE metread(lfn, MetArray, InputmetFormat, ldown_option, NetRadiationMethod, &
                   SnowUse, SMDMethod, SoilDepthMeas, SoilRocks, SoilDensity, SmCap)
 
   USE defaultnotused
 
   IMPLICIT NONE
 
   !INPUT
   REAL(KIND(1D0)), DIMENSION(24) :: MetArray !Array leaving the subroutine within
   !each INTERVAL (defined in RunControl.nml)
   ! - Met data now provided at a resolution of tstep, HCW Jan 2015
 
   REAL(KIND(1D0)) :: SmCap, & !Volumetric/gravimetric soil moisture capacity
                      SoilDepthMeas, & !Measured soil depth
                      SoilRocks, & !Rocks on ground
                      SoilDensity !Density of soil
 
   INTEGER :: InputmetFormat, & !Format of the meteorological forcing file
              ldown_option, & !Method of calculating Ldown
              NetRadiationMethod, & !Method of calculating Q*
              SMDMethod, & !Method of measured soil moisture
              SnowUse
 
   ! Variables read in
   REAL(KIND(1D0)) :: avkdn, & !Average downwelling shortwave radiation
                      avrh, & !Average relative humidity
                      avu1, & !Average wind speed
                      dectime, & !Decimal time
                      fcld_obs, & !Cloud fraction observed
                      iy, & !Year
                      id, & !Day
                      it, & !Hour
                      imin, & !Minute
                      kdiff, & !Diffuse shortwave radiation
                      kdir, & !Direct shortwave radiation
                      LAI_obs, & !Overall LAI of the study area
                      ldown_obs, & !Downwelling longwave radiation
                      Precip, & !Rainfall [mm]
                      Pres_hPa, & !Station air pressure in hPa
                      Pres_kPa, & !Station air pressure in kPa
                      snowFrac_obs, & !Observed surface fraction of snow (between 0 and 1)
                      qe_obs, & !Observed latent heat flux
                      qf_obs, & !Observed antrhropogeni heat flux
                      qh_obs, & !Observed sensible heat flux
                      qn1_obs, & !Observed net all-wave radiation
                      qs_obs, & !Observed storage heat flux
                      Temp_C, & !Air temperature
                      wdir, & !Wind direction
                      wu_m3, & !Water use provided in met forcing file [m3]
                      xsmd !Measured soil moisture deficit
 
   INTEGER :: iostat_var, lfn
 
   !-----------------------------------------------------------------------------------
   !-----------------------------------------------------------------------------------
 
   IF (inputmetformat == 0) THEN !Default format using LUMPS only
 
      READ (lfn, *, iostat=iostat_var) iy, id, it, imin, qn1_obs, avu1, avrh, &
         temp_c, wdir, pres_kpa, precip, avkdn, snowfrac_obs, ldown_obs, fcld_obs
 
      !Set other variables needed while running SUEWS to zero
      qf_obs = nan
      qs_obs = nan
      qh_obs = nan
      qe_obs = nan
      xsmd = -99999
      kdiff = nan
      kdir = nan
      wdir = nan
 
   ELSEIF (inputmetformat == 10) THEN !SUEWS reading
      READ (lfn, *, iostat=iostat_var) iy, id, it, imin, qn1_obs, qh_obs, qe_obs, qs_obs, qf_obs, avu1, avrh, &
         temp_c, pres_kpa, precip, avkdn, snowfrac_obs, ldown_obs, fcld_obs, &
         wu_m3, xsmd, lai_obs, kdiff, kdir, wdir
 
      !write(*,*) 'In LUMPS_MetRead (1)'
      !write(*,*) 'imin',imin
      !write(*,*) 'it',it
      !write(*,*) 'id',id
      !write(*,*) 'iy',iy
 
      !Calculate observed soil moisture deficits from either volumetric or gravimetric SoilStates
      IF (smdmethod == 1 .AND. xsmd /= -999) THEN !Soil moisture - volumetric
         xsmd = (smcap - xsmd)*soildepthmeas*soilrocks
      ELSEIF (smdmethod == 2 .AND. xsmd /= -999) THEN !Soil moisture -gravimetric
         xsmd = (smcap - xsmd)*soildensity*soildepthmeas*soilrocks
      ELSE
         xsmd = -999
      END IF
 
   ELSE
      CALL errorhint(55, 'RunControl.nml, InputMetFormat not usable.', notused, notused, inputmetformat)
   END IF
 
   !===============Meteorological variables reading done==========================
   pres_hpa = pres_kpa*10. ! convert to hPa
 
   ! get decimal time
   dectime = id + it + imin/60.0
 
   IF (iostat_var < 0) THEN
      iostat_var = 0
      CLOSE (lfn)
      RETURN
   END IF
 
   IF (avkdn < 0) THEN
      CALL errorhint(27, 'Met Data: avKdn - needed for StoreDrainPrm. resistance, If present, check file not tab delimited', &
                     avkdn, dectime, notusedi)
      !sg removed this is causing the problems with resistances
      !  AvKdn=0 !Solar radiation cannot be lower than 1
   END IF
 
   IF ((ldown_option == 1) .AND. (ldown_obs < 0)) THEN
      CALL errorhint(27, 'Met Data: LWdn (ldown_obs) - impact Q* calc', ldown_obs, dectime, notusedi)
 
   ELSEIF (ldown_option == 2) THEN
      IF (fcld_obs == -999.0 .OR. fcld_obs < 0 .OR. fcld_obs > 1) THEN
         CALL errorhint(27, 'Met Data: flcd_obs - impacts LW & Q* radiation', fcld_obs, dectime, notusedi)
      END IF
   END IF
 
   IF (qn1_obs == -999 .AND. netradiationmethod == 0) THEN !If measured Q* is used and it is -999
      CALL errorhint(27, 'Met Data: Q* - will impact everything', qn1_obs, dectime, notusedi)
   END IF
 
   IF (avu1 <= 0) THEN !If wind speed is negative
      CALL errorhint(27, 'Met Data: avU1 - impacts aeroydnamic resistances', avu1, dectime, notusedi)
   END IF
 
   IF (temp_c < -50 .OR. temp_c > 60) THEN !If temperature unrealistic
      CALL errorhint(27, 'Met Data: Temp_C - beyond what is expected', temp_c, dectime, notusedi)
   END IF
 
   IF (avrh > 100 .OR. avrh < 1) THEN !If relative humidity larger than 100%
      CALL errorhint(27, 'Met Data: avRH - beyond what is expected', avrh, dectime, notusedi)
   END IF
 
   IF (pres_kpa < 80) THEN !If pressure too low
      CALL errorhint(27, 'Met Data: Pres_kPa - too low - this could be fixed in model', pres_kpa, dectime, notusedi)
   END IF
 
   IF (precip < 0) THEN !If rain in negative, set it to zero
      CALL errorhint(27, 'Met Data: Precip - less than 0', precip, dectime, notusedi)
   END IF
 
   IF (snowfrac_obs == nan) snowfrac_obs = 0
 
   IF (snowuse == 0 .AND. (snowfrac_obs < 0 .OR. snowfrac_obs > 1)) THEN
      CALL errorhint(27, 'Met Data: snow not between [0  1]', snowfrac_obs, dectime, notusedi)
   END IF
 
   IF (xsmd < 0 .AND. smdmethod == 1) THEN !If soil moisture deficit is zero
      CALL errorhint(27, 'Met Data: xsmd - less than 0', xsmd, dectime, notusedi)
   END IF
 
   !Create an array to be printed out.
   metarray(1:24) = (/iy, id, it, imin, qn1_obs, qh_obs, qe_obs, qs_obs, qf_obs, avu1, &
                      avrh, temp_c, pres_hpa, precip, avkdn, snowfrac_obs, ldown_obs, &
                      fcld_obs, wu_m3, xsmd, lai_obs, kdiff, kdir, wdir/)
 
   !write(*,*) 'In LUMPS_MetRead (2)'
   !write(*,*) 'imin',imin
   !write(*,*) 'it',it
   !write(*,*) 'id',id
   !write(*,*) 'iy',iy
 
   RETURN
 
END SUBROUTINE metread
 
! sg Feb 2012 -- moving all related suborutine together
! subroutine run_control - checks that value is reasonable
! subrotuine skipHeader  - jumps over header rows in input files
! Last modified:
! LJ 27 Jan 2016 - Removal of tabs, cleaning of the code
!------------------------------------------------------------
 
!Information for the run
MODULE run_info
   IMPLICIT NONE
   CHARACTER(len=90), DIMENSION(14) :: text
   INTEGER :: lim0 = 0, lim1 = 1, lim2 = 2, lim4 = 4, lim3 = 3, lim6 = 6, lim8 = 8, lim12 = 12, lfn_us
   LOGICAL :: file_qs
END MODULE run_info
 
! run_control
! called from: LUMPS_initial
 
SUBROUTINE run_control(eval, LowerLimit, Upperlimit)
   ! ver - determines if value to be read is an integer or real and returns the value
   ! if ver=-9 - then use integer
   USE run_info
   IMPLICIT NONE
   INTEGER :: eval, i, lowerlimit, upperlimit
   CHARACTER(len=4) :: check
 
   IF (file_qs) THEN
101   READ (lfn_us, *) check
      WRITE (*, *) check
      DO i = 1, 3
         IF (check(i:i) == "#") THEN
            !   write(*,*)check(i:i),i,"check"
            GOTO 101
         ELSE
            backspace(lfn_us)
            READ (lfn_us, *) eval
            !write(*,*)eval, TEXT(1)
            EXIT
         END IF
      END DO
   END IF
 
   WRITE (12, 120) eval, text(1)
 
   IF (eval < lowerlimit .OR. eval > upperlimit) THEN
      WRITE (*, *) "Value out of range"
      WRITE (*, *) eval, text(1)
      stop
   END IF
 
   WRITE (*, 120) eval, text(1)
120 FORMAT(i4, 2x, a90)
 
   RETURN
END SUBROUTINE run_control
 
! Skip Headers------------------------------
SUBROUTINE skipheader(lfn, skip)
   USE defaultnotused
   IMPLICIT NONE
 
   INTEGER :: skip, lfn, i
   DO i = 1, skip
      READ (lfn, *, err=201, iostat=ios_out)
   END DO
 
   RETURN
 
201 reall = real(skip)
   CALL errorhint(20, 'In SkipHeader subroutine.', reall, notused, ios_out)
END SUBROUTINE skipheader
 
!-------------------------------------------------------------------------
SUBROUTINE inputheadercheck(FileName)
   ! Checks columns in input files match the columns expected by model code
   ! Model code columns are defined here
   ! Latest update:
   !   MH 21 Jun 2017 - Added parameters to SUEWS_AnthropogenicEmissions.txt
   !   MH 16 Jun 2017 - Added SUEWS_BiogenCO2.txt
   !   TS 02 Mar 2016  - AnOHM related variables added
   !   LJ 27 Jan 2016  - Removal of tabs
   !   LJ 07 July 2015 - snow albedo removed
   !   HCW 12 Nov 2014
   !-------------------------------------------------------------------------
 
   USE allocatearray
   USE colnamesinputfiles
   USE defaultnotused
 
   IMPLICIT NONE
 
   CHARACTER(len=50) :: FileName
 
   ! ========== Define expected column names here ==========
   ! =======================================================
 
   ! ========== SUEWS_NonVeg.txt =============
   headernonveg_reqd(ci_code) = "Code"
   headernonveg_reqd(ci_albmin) = "AlbedoMin"
   headernonveg_reqd(ci_albmax) = "AlbedoMax"
   headernonveg_reqd(ci_emis) = "Emissivity"
   headernonveg_reqd(ci_stormin) = "StorageMin"
   headernonveg_reqd(ci_stormax) = "StorageMax"
   headernonveg_reqd(ci_wetthresh) = "WetThreshold"
   headernonveg_reqd(ci_statelimit) = "StateLimit"
   headernonveg_reqd(ci_dreq) = "DrainageEq"
   headernonveg_reqd(ci_drcoef1) = "DrainageCoef1"
   headernonveg_reqd(ci_drcoef2) = "DrainageCoef2"
   headernonveg_reqd(ci_soiltcode) = "SoilTypeCode"
   headernonveg_reqd(ci_snowlimpat) = "SnowLimPatch"
   headernonveg_reqd(ci_snowlimrem) = "SnowLimRemove"
   headernonveg_reqd(ci_ohmcode_swet) = "OHMCode_SummerWet"
   headernonveg_reqd(ci_ohmcode_sdry) = "OHMCode_SummerDry"
   headernonveg_reqd(ci_ohmcode_wwet) = "OHMCode_WinterWet"
   headernonveg_reqd(ci_ohmcode_wdry) = "OHMCode_WinterDry"
   headernonveg_reqd(ci_ohmthresh_sw) = "OHMThresh_SW"
   headernonveg_reqd(ci_ohmthresh_wd) = "OHMThresh_WD"
   headernonveg_reqd(ci_estmcode) = "ESTMCode"
   headernonveg_reqd(ci_cpanohm) = "AnOHM_Cp" ! AnOHM TS
   headernonveg_reqd(ci_kkanohm) = "AnOHM_Kk" ! AnOHM TS
   headernonveg_reqd(ci_chanohm) = "AnOHM_Ch" ! AnOHM TS
 
   ! ========== SUEWS_Veg.txt ===============
   headerveg_reqd(cp_code) = "Code"
   headerveg_reqd(cp_albmin) = "AlbedoMin"
   headerveg_reqd(cp_albmax) = "AlbedoMax"
   headerveg_reqd(cp_emis) = "Emissivity"
   headerveg_reqd(cp_stormin) = "StorageMin"
   headerveg_reqd(cp_stormax) = "StorageMax"
   headerveg_reqd(cp_wetthresh) = "WetThreshold"
   headerveg_reqd(cp_statelimit) = "StateLimit"
   headerveg_reqd(cp_dreq) = "DrainageEq"
   headerveg_reqd(cp_drcoef1) = "DrainageCoef1"
   headerveg_reqd(cp_drcoef2) = "DrainageCoef2"
   headerveg_reqd(cp_soiltcode) = "SoilTypeCode"
   headerveg_reqd(cp_snowlimpat) = "SnowLimPatch"
   headerveg_reqd(cp_baset) = "BaseT"
   headerveg_reqd(cp_basete) = "BaseTe"
   headerveg_reqd(cp_gddfull) = "GDDFull"
   headerveg_reqd(cp_sddfull) = "SDDFull"
   headerveg_reqd(cp_laimin) = "LAIMin"
   headerveg_reqd(cp_laimax) = "LAIMax"
   headerveg_reqd(cp_porositymin) = "PorosityMin"
   headerveg_reqd(cp_porositymax) = "PorosityMax"
   headerveg_reqd(cp_gsmax) = "MaxConductance"
   headerveg_reqd(cp_laieq) = "LAIEq"
   headerveg_reqd(cp_leafgp1) = "LeafGrowthPower1"
   headerveg_reqd(cp_leafgp2) = "LeafGrowthPower2"
   headerveg_reqd(cp_leafop1) = "LeafOffPower1"
   headerveg_reqd(cp_leafop2) = "LeafOffPower2"
   headerveg_reqd(cp_ohmcode_swet) = "OHMCode_SummerWet"
   headerveg_reqd(cp_ohmcode_sdry) = "OHMCode_SummerDry"
   headerveg_reqd(cp_ohmcode_wwet) = "OHMCode_WinterWet"
   headerveg_reqd(cp_ohmcode_wdry) = "OHMCode_WinterDry"
   headerveg_reqd(cp_ohmthresh_sw) = "OHMThresh_SW"
   headerveg_reqd(cp_ohmthresh_wd) = "OHMThresh_WD"
   headerveg_reqd(cp_estmcode) = "ESTMCode"
   headerveg_reqd(cp_cpanohm) = "AnOHM_Cp" ! AnOHM TS
   headerveg_reqd(cp_kkanohm) = "AnOHM_Kk" ! AnOHM TS
   headerveg_reqd(cp_chanohm) = "AnOHM_Ch" ! AnOHM TS
   headerveg_reqd(cp_biogenco2code) = "BiogenCO2Code"
 
   ! ========== SUEWS_Water.txt ==================
   headerwater_reqd(cw_code) = "Code"
   headerwater_reqd(cw_albmin) = "AlbedoMin"
   headerwater_reqd(cw_albmax) = "AlbedoMax"
   headerwater_reqd(cw_emis) = "Emissivity"
   headerwater_reqd(cw_stormin) = "StorageMin"
   headerwater_reqd(cw_stormax) = "StorageMax"
   headerwater_reqd(cw_wetthresh) = "WetThreshold"
   headerwater_reqd(cw_statelimit) = "StateLimit"
   headerwater_reqd(cw_waterdepth) = "WaterDepth"
   headerwater_reqd(cw_dreq) = "DrainageEq"
   headerwater_reqd(cw_drcoef1) = "DrainageCoef1"
   headerwater_reqd(cw_drcoef2) = "DrainageCoef2"
   headerwater_reqd(cw_ohmcode_swet) = "OHMCode_SummerWet"
   headerwater_reqd(cw_ohmcode_sdry) = "OHMCode_SummerDry"
   headerwater_reqd(cw_ohmcode_wwet) = "OHMCode_WinterWet"
   headerwater_reqd(cw_ohmcode_wdry) = "OHMCode_WinterDry"
   headerwater_reqd(cw_ohmthresh_sw) = "OHMThresh_SW"
   headerwater_reqd(cw_ohmthresh_wd) = "OHMThresh_WD"
   headerwater_reqd(cw_estmcode) = "ESTMCode"
   headerwater_reqd(cw_cpanohm) = "AnOHM_Cp" ! AnOHM TS
   headerwater_reqd(cw_kkanohm) = "AnOHM_Kk" ! AnOHM TS
   headerwater_reqd(cw_chanohm) = "AnOHM_Ch" ! AnOHM TS
 
   ! ========== SUEWS_Snow.txt ===================
   headersnow_reqd(cs_code) = "Code"
   headersnow_reqd(cs_snowrmfactor) = "RadMeltFactor"
   headersnow_reqd(cs_snowtmfactor) = "TempMeltFactor"
   headersnow_reqd(cs_snowalbmin) = "AlbedoMin"
   headersnow_reqd(cs_snowalbmax) = "AlbedoMax"
   headersnow_reqd(cs_snowemis) = "Emissivity"
   headersnow_reqd(cs_snowtau_a) = "tau_a"
   headersnow_reqd(cs_snowtau_f) = "tau_f"
   headersnow_reqd(cs_snowplimalb) = "PrecipLimAlb"
   headersnow_reqd(cs_snowsdmin) = "SnowDensMin"
   headersnow_reqd(cs_snowsdmax) = "SnowDensMax"
   headersnow_reqd(cs_snowtau_r) = "tau_r"
   headersnow_reqd(cs_snowcrwmin) = "CRWMin"
   headersnow_reqd(cs_snowcrwmax) = "CRWMax"
   headersnow_reqd(cs_snowplimsnow) = "PrecipLimSnow"
   headersnow_reqd(cs_ohmcode_swet) = "OHMCode_SummerWet"
   headersnow_reqd(cs_ohmcode_sdry) = "OHMCode_SummerDry"
   headersnow_reqd(cs_ohmcode_wwet) = "OHMCode_WinterWet"
   headersnow_reqd(cs_ohmcode_wdry) = "OHMCode_WinterDry"
   headersnow_reqd(cs_ohmthresh_sw) = "OHMThresh_SW"
   headersnow_reqd(cs_ohmthresh_wd) = "OHMThresh_WD"
   headersnow_reqd(cs_estmcode) = "ESTMCode"
   headersnow_reqd(cs_cpanohm) = "AnOHM_Cp" ! AnOHM TS
   headersnow_reqd(cs_kkanohm) = "AnOHM_Kk" ! AnOHM TS
   headersnow_reqd(cs_chanohm) = "AnOHM_Ch" ! AnOHM TS
 
   ! ========== SUEWS_Soil.txt ===================
   headersoil_reqd(cso_code) = "Code"
   headersoil_reqd(cso_soildepth) = "SoilDepth"
   headersoil_reqd(cso_soilstcap) = "SoilStoreCap"
   headersoil_reqd(cso_ksat) = "SatHydraulicCond"
   headersoil_reqd(cso_soildens) = "SoilDensity"
   headersoil_reqd(cso_soilinfrate) = "InfiltrationRate"
   headersoil_reqd(cso_obssmdepth) = "OBS_SMDepth"
   headersoil_reqd(cso_obssmmax) = "OBS_SMCap"
   headersoil_reqd(cso_obssnrfrac) = "OBS_SoilNotRocks"
 
   ! ========== SUEWS_Conductance.txt ============
   headercond_reqd(cc_code) = "Code"
   headercond_reqd(cc_gsg1) = "G1"
   headercond_reqd(cc_gsg2) = "G2"
   headercond_reqd(cc_gsg3) = "G3"
   headercond_reqd(cc_gsg4) = "G4"
   headercond_reqd(cc_gsg5) = "G5"
   headercond_reqd(cc_gsg6) = "G6"
   headercond_reqd(cc_gsth) = "TH"
   headercond_reqd(cc_gstl) = "TL"
   headercond_reqd(cc_gss1) = "S1"
   headercond_reqd(cc_gss2) = "S2"
   headercond_reqd(cc_gskmax) = "Kmax"
   headercond_reqd(cc_gsmodel) = "gsModel"
 
   ! ========== SUEWS_OHMCoefficients.txt ========
   headerohmcoefficients_reqd(co_code) = "Code"
   headerohmcoefficients_reqd(co_a1) = "a1"
   headerohmcoefficients_reqd(co_a2) = "a2"
   headerohmcoefficients_reqd(co_a3) = "a3"
 
   ! ========== SUEWS_ESTMCoefficients.txt ========
   headerestmcoefficients_reqd(ce_code) = "Code"
   headerestmcoefficients_reqd(ce_surf_thick1) = "Surf_thick1"
   headerestmcoefficients_reqd(ce_surf_k1) = "Surf_k1"
   headerestmcoefficients_reqd(ce_surf_rhocp1) = "Surf_rhoCp1"
   headerestmcoefficients_reqd(ce_surf_thick2) = "Surf_thick2"
   headerestmcoefficients_reqd(ce_surf_k2) = "Surf_k2"
   headerestmcoefficients_reqd(ce_surf_rhocp2) = "Surf_rhoCp2"
   headerestmcoefficients_reqd(ce_surf_thick3) = "Surf_thick3"
   headerestmcoefficients_reqd(ce_surf_k3) = "Surf_k3"
   headerestmcoefficients_reqd(ce_surf_rhocp3) = "Surf_rhoCp3"
   headerestmcoefficients_reqd(ce_surf_thick4) = "Surf_thick4"
   headerestmcoefficients_reqd(ce_surf_k4) = "Surf_k4"
   headerestmcoefficients_reqd(ce_surf_rhocp4) = "Surf_rhoCp4"
   headerestmcoefficients_reqd(ce_surf_thick5) = "Surf_thick5"
   headerestmcoefficients_reqd(ce_surf_k5) = "Surf_k5"
   headerestmcoefficients_reqd(ce_surf_rhocp5) = "Surf_rhoCp5"
   headerestmcoefficients_reqd(ce_wall_thick1) = "Wall_thick1"
   headerestmcoefficients_reqd(ce_wall_k1) = "Wall_k1"
   headerestmcoefficients_reqd(ce_wall_rhocp1) = "Wall_rhoCp1"
   headerestmcoefficients_reqd(ce_wall_thick2) = "Wall_thick2"
   headerestmcoefficients_reqd(ce_wall_k2) = "Wall_k2"
   headerestmcoefficients_reqd(ce_wall_rhocp2) = "Wall_rhoCp2"
   headerestmcoefficients_reqd(ce_wall_thick3) = "Wall_thick3"
   headerestmcoefficients_reqd(ce_wall_k3) = "Wall_k3"
   headerestmcoefficients_reqd(ce_wall_rhocp3) = "Wall_rhoCp3"
   headerestmcoefficients_reqd(ce_wall_thick4) = "Wall_thick4"
   headerestmcoefficients_reqd(ce_wall_k4) = "Wall_k4"
   headerestmcoefficients_reqd(ce_wall_rhocp4) = "Wall_rhoCp4"
   headerestmcoefficients_reqd(ce_wall_thick5) = "Wall_thick5"
   headerestmcoefficients_reqd(ce_wall_k5) = "Wall_k5"
   headerestmcoefficients_reqd(ce_wall_rhocp5) = "Wall_rhoCp5"
   headerestmcoefficients_reqd(ce_internal_thick1) = "Internal_thick1"
   headerestmcoefficients_reqd(ce_internal_k1) = "Internal_k1"
   headerestmcoefficients_reqd(ce_internal_rhocp1) = "Internal_rhoCp1"
   headerestmcoefficients_reqd(ce_internal_thick2) = "Internal_thick2"
   headerestmcoefficients_reqd(ce_internal_k2) = "Internal_k2"
   headerestmcoefficients_reqd(ce_internal_rhocp2) = "Internal_rhoCp2"
   headerestmcoefficients_reqd(ce_internal_thick3) = "Internal_thick3"
   headerestmcoefficients_reqd(ce_internal_k3) = "Internal_k3"
   headerestmcoefficients_reqd(ce_internal_rhocp3) = "Internal_rhoCp3"
   headerestmcoefficients_reqd(ce_internal_thick4) = "Internal_thick4"
   headerestmcoefficients_reqd(ce_internal_k4) = "Internal_k4"
   headerestmcoefficients_reqd(ce_internal_rhocp4) = "Internal_rhoCp4"
   headerestmcoefficients_reqd(ce_internal_thick5) = "Internal_thick5"
   headerestmcoefficients_reqd(ce_internal_k5) = "Internal_k5"
   headerestmcoefficients_reqd(ce_internal_rhocp5) = "Internal_rhoCp5"
   headerestmcoefficients_reqd(ce_nroom) = "nroom"
   headerestmcoefficients_reqd(ce_alb_ibld) = "Internal_albedo"
   headerestmcoefficients_reqd(ce_em_ibld) = "Internal_emissivity"
   headerestmcoefficients_reqd(ce_ch_iwall) = "Internal_CHwall"
   headerestmcoefficients_reqd(ce_ch_iroof) = "Internal_CHroof"
   headerestmcoefficients_reqd(ce_ch_ibld) = "Internal_CHbld"
 
   ! ========== SUEWS_AnthropogenicEmission.txt ======
   headeranthropogenic_reqd(ca_code) = "Code"
   headeranthropogenic_reqd(ca_baset_hc) = "BaseT_HC"
   headeranthropogenic_reqd(ca_qf_a1) = "QF_A_WD"
   headeranthropogenic_reqd(ca_qf_b1) = "QF_B_WD"
   headeranthropogenic_reqd(ca_qf_c1) = "QF_C_WD"
   headeranthropogenic_reqd(ca_qf_a2) = "QF_A_WE"
   headeranthropogenic_reqd(ca_qf_b2) = "QF_B_WE"
   headeranthropogenic_reqd(ca_qf_c2) = "QF_C_WE"
   headeranthropogenic_reqd(ca_ahmin_wd) = "AHMin_WD"
   headeranthropogenic_reqd(ca_ahmin_we) = "AHMin_WE"
   headeranthropogenic_reqd(ca_ahslopeheating_wd) = "AHSlope_Heating_WD"
   headeranthropogenic_reqd(ca_ahslopeheating_we) = "AHSlope_Heating_WE"
   headeranthropogenic_reqd(ca_ahslopecooling_wd) = "AHSlope_Cooling_WD"
   headeranthropogenic_reqd(ca_ahslopecooling_we) = "AHSlope_Cooling_WE"
   headeranthropogenic_reqd(ca_tcriticheating_wd) = "TCritic_Heating_WD"
   headeranthropogenic_reqd(ca_tcriticheating_we) = "TCritic_Heating_WE"
   headeranthropogenic_reqd(ca_tcriticcooling_wd) = "TCritic_Cooling_WD"
   headeranthropogenic_reqd(ca_tcriticcooling_we) = "TCritic_Cooling_WE"
   headeranthropogenic_reqd(ca_enprofwd) = "EnergyUseProfWD"
   headeranthropogenic_reqd(ca_enprofwe) = "EnergyUseProfWE"
   headeranthropogenic_reqd(ca_co2mwd) = "ActivityProfWD"
   headeranthropogenic_reqd(ca_co2mwe) = "ActivityProfWE"
   headeranthropogenic_reqd(ca_traffprofwd) = "TraffProfWD"
   headeranthropogenic_reqd(ca_traffprofwe) = "TraffProfWE"
   headeranthropogenic_reqd(ca_popprofwd) = "PopProfWD"
   headeranthropogenic_reqd(ca_popprofwe) = "PopProfWE"
   headeranthropogenic_reqd(ca_minqfmetab) = "MinQFMetab"
   headeranthropogenic_reqd(ca_maxqfmetab) = "MaxQFMetab"
   headeranthropogenic_reqd(ca_minfcmetab) = "MinFCMetab"
   headeranthropogenic_reqd(ca_maxfcmetab) = "MaxFCMetab"
   headeranthropogenic_reqd(ca_frpddwe) = "FrPDDwe"
   headeranthropogenic_reqd(ca_frfossilfuel_heat) = "FrFossilFuel_Heat"
   headeranthropogenic_reqd(ca_frfossilfuel_nonheat) = "FrFossilFuel_NonHeat"
   headeranthropogenic_reqd(ca_ef_umolco2perj) = "EF_umolCO2perJ"
   headeranthropogenic_reqd(ca_enef_v_jkm) = "EnEF_v_Jkm"
   headeranthropogenic_reqd(ca_fcef_v_kgkmwd) = "FcEF_v_kgkmWD"
   headeranthropogenic_reqd(ca_fcef_v_kgkmwe) = "FcEF_v_kgkmWE"
   headeranthropogenic_reqd(ca_co2pointsource) = "CO2PointSource"
   headeranthropogenic_reqd(ca_trafficunits) = "TrafficUnits"
 
   ! ========== SUEWS_Irrigation.txt =============
   headerirrigation_reqd(cir_code) = "Code"
   headerirrigation_reqd(cir_iestart) = "Ie_start"
   headerirrigation_reqd(cir_ieend) = "Ie_end"
   headerirrigation_reqd(cir_intwu) = "InternalWaterUse"
   headerirrigation_reqd(cir_faut) = "Faut"
   headerirrigation_reqd(cir_h_maintain) = "H_maintain"
   headerirrigation_reqd(cir_ie_a1) = "Ie_a1"
   headerirrigation_reqd(cir_ie_a2) = "Ie_a2"
   headerirrigation_reqd(cir_ie_a3) = "Ie_a3"
   headerirrigation_reqd(cir_ie_m1) = "Ie_m1"
   headerirrigation_reqd(cir_ie_m2) = "Ie_m2"
   headerirrigation_reqd(cir_ie_m3) = "Ie_m3"
   headerirrigation_reqd(cir_daywat1) = "DayWat(1)"
   headerirrigation_reqd(cir_daywat2) = "DayWat(2)"
   headerirrigation_reqd(cir_daywat3) = "DayWat(3)"
   headerirrigation_reqd(cir_daywat4) = "DayWat(4)"
   headerirrigation_reqd(cir_daywat5) = "DayWat(5)"
   headerirrigation_reqd(cir_daywat6) = "DayWat(6)"
   headerirrigation_reqd(cir_daywat7) = "DayWat(7)"
   headerirrigation_reqd(cir_daywatper1) = "DayWatPer(1)"
   headerirrigation_reqd(cir_daywatper2) = "DayWatPer(2)"
   headerirrigation_reqd(cir_daywatper3) = "DayWatPer(3)"
   headerirrigation_reqd(cir_daywatper4) = "DayWatPer(4)"
   headerirrigation_reqd(cir_daywatper5) = "DayWatPer(5)"
   headerirrigation_reqd(cir_daywatper6) = "DayWatPer(6)"
   headerirrigation_reqd(cir_daywatper7) = "DayWatPer(7)"
 
   ! ========== SUEWS_Profiles.txt ===============
   headerprofiles_reqd(cpr_code) = "Code"
   headerprofiles_reqd(cpr_hours(1)) = "0"
   headerprofiles_reqd(cpr_hours(2)) = "1"
   headerprofiles_reqd(cpr_hours(3)) = "2"
   headerprofiles_reqd(cpr_hours(4)) = "3"
   headerprofiles_reqd(cpr_hours(5)) = "4"
   headerprofiles_reqd(cpr_hours(6)) = "5"
   headerprofiles_reqd(cpr_hours(7)) = "6"
   headerprofiles_reqd(cpr_hours(8)) = "7"
   headerprofiles_reqd(cpr_hours(9)) = "8"
   headerprofiles_reqd(cpr_hours(10)) = "9"
   headerprofiles_reqd(cpr_hours(11)) = "10"
   headerprofiles_reqd(cpr_hours(12)) = "11"
   headerprofiles_reqd(cpr_hours(13)) = "12"
   headerprofiles_reqd(cpr_hours(14)) = "13"
   headerprofiles_reqd(cpr_hours(15)) = "14"
   headerprofiles_reqd(cpr_hours(16)) = "15"
   headerprofiles_reqd(cpr_hours(17)) = "16"
   headerprofiles_reqd(cpr_hours(18)) = "17"
   headerprofiles_reqd(cpr_hours(19)) = "18"
   headerprofiles_reqd(cpr_hours(20)) = "19"
   headerprofiles_reqd(cpr_hours(21)) = "20"
   headerprofiles_reqd(cpr_hours(22)) = "21"
   headerprofiles_reqd(cpr_hours(23)) = "22"
   headerprofiles_reqd(cpr_hours(24)) = "23"
 
   ! ========== SUEWS_WithinGridWaterDist.txt ====
   headerwgwaterdist_reqd(cwg_code) = "Code"
   headerwgwaterdist_reqd(cwg_topaved) = "ToPaved"
   headerwgwaterdist_reqd(cwg_tobldgs) = "ToBldgs"
   headerwgwaterdist_reqd(cwg_toevetr) = "ToEveTr"
   headerwgwaterdist_reqd(cwg_todectr) = "ToDecTr"
   headerwgwaterdist_reqd(cwg_tograss) = "ToGrass"
   headerwgwaterdist_reqd(cwg_tobsoil) = "ToBSoil"
   headerwgwaterdist_reqd(cwg_towater) = "ToWater"
   headerwgwaterdist_reqd(cwg_torunoff) = "ToRunoff"
   headerwgwaterdist_reqd(cwg_tosoilstore) = "ToSoilStore"
 
   ! ========== SUEWS_BiogenCO2.txt ======
   headerbiogen_reqd(cb_code) = "Code"
   headerbiogen_reqd(cb_alpha) = "alpha"
   headerbiogen_reqd(cb_beta) = "beta"
   headerbiogen_reqd(cb_theta) = "theta"
   headerbiogen_reqd(cb_alpha_enh) = "alpha_enh"
   headerbiogen_reqd(cb_beta_enh) = "beta_enh"
   headerbiogen_reqd(cb_resp_a) = "resp_a"
   headerbiogen_reqd(cb_resp_b) = "resp_b"
   headerbiogen_reqd(cb_min_r) = "min_respi"
 
   ! =======================================================
 
   !write(*,*) 'Checking header for ', FileName
   ! Check columns in input files match model code
 
   IF (filename == 'SUEWS_NonVeg.txt') THEN
      IF (any(headernonveg_file /= headernonveg_reqd)) THEN
         WRITE (*, *) headernonveg_file == headernonveg_reqd
         WRITE (*, *) headernonveg_file
         WRITE (*, *) headernonveg_reqd
         CALL errorhint(56, 'Names or order of columns in SUEWS_NonVeg.txt does not match model code.', notused, notused, notusedi)
      END IF
 
   ELSEIF (filename == 'SUEWS_Veg.txt') THEN
      IF (any(headerveg_file /= headerveg_reqd)) THEN
         WRITE (*, *) headerveg_file == headerveg_reqd
         CALL errorhint(56, 'Names or order of columns in SUEWS_Veg.txt does not match model code.', notused, notused, notusedi)
      END IF
 
   ELSEIF (filename == 'SUEWS_Water.txt') THEN
      IF (any(headerwater_file /= headerwater_reqd)) THEN
         WRITE (*, *) headerwater_file == headerwater_reqd
         WRITE (*, *) headerwater_file
         WRITE (*, *) headerwater_reqd
         CALL errorhint(56, 'Names or order of columns in SUEWS_Water.txt does not match model code.', notused, notused, notusedi)
      END IF
 
   ELSEIF (filename == 'SUEWS_Snow.txt') THEN
      IF (any(headersnow_file /= headersnow_reqd)) THEN
         WRITE (*, *) headersnow_file == headersnow_reqd
         CALL errorhint(56, 'Names or order of columns in SUEWS_Snow.txt does not match model code.', notused, notused, notusedi)
      END IF
 
   ELSEIF (filename == 'SUEWS_Soil.txt') THEN
      IF (any(headersoil_file /= headersoil_reqd)) THEN
         WRITE (*, *) headersoil_file == headersoil_reqd
         CALL errorhint(56, 'Names or order of columns in SUEWS_Soil.txt does not match model code.', notused, notused, notusedi)
      END IF
 
   ELSEIF (filename == 'SUEWS_Conductance.txt') THEN
      IF (any(headercond_file /= headercond_reqd)) THEN
         WRITE (*, *) headercond_file == headercond_reqd
         CALL errorhint(56, 'Names or order of columns in SUEWS_Cond.txt does not match model code.', notused, notused, notusedi)
      END IF
 
   ELSEIF (filename == 'SUEWS_OHMCoefficients.txt') THEN
      IF (any(headerohmcoefficients_file /= headerohmcoefficients_reqd)) THEN
         WRITE (*, *) headerohmcoefficients_file == headerohmcoefficients_reqd
         CALL errorhint(56, 'Names or order of columns in SUEWS_OHMCoefficients.txt does not match model code.', &
                        notused, notused, notusedi)
      END IF
 
   ELSEIF (filename == 'SUEWS_ESTMCoefficients.txt') THEN
      IF (any(headerestmcoefficients_file /= headerestmcoefficients_reqd)) THEN
         WRITE (*, *) headerestmcoefficients_file == headerestmcoefficients_reqd
         CALL errorhint(56, 'Names or order of columns in SUEWS_ESTMCoefficients.txt does not match model code.', &
                        notused, notused, notusedi)
      END IF
 
   ELSEIF (filename == 'SUEWS_AnthropogenicEmission.txt') THEN
      IF (any(headeranthropogenic_file /= headeranthropogenic_reqd)) THEN
         WRITE (*, *) headeranthropogenic_file == headeranthropogenic_reqd
         CALL errorhint(56, 'Names or order of columns in SUEWS_AnthropogenicEmission.txt does not match model code.', &
                        notused, notused, notusedi)
      END IF
 
   ELSEIF (filename == 'SUEWS_Irrigation.txt') THEN
      IF (any(headerirrigation_file /= headerirrigation_reqd)) THEN
         WRITE (*, *) headerirrigation_file == headerirrigation_reqd
         WRITE (*, *) headerirrigation_file
         WRITE (*, *) headerirrigation_reqd
         CALL errorhint(56, 'Names or order of columns in SUEWS_Irrigation.txt does not match model code.', &
                        notused, notused, notusedi)
      END IF
 
   ELSEIF (filename == 'SUEWS_Profiles.txt') THEN
      IF (any(headerprofiles_file /= headerprofiles_reqd)) THEN
         WRITE (*, *) headerprofiles_file == headerprofiles_reqd
         CALL errorhint(56, 'Names or order of columns in SUEWS_Profiles.txt does not match model code.', &
                        notused, notused, notusedi)
      END IF
 
   ELSEIF (filename == 'SUEWS_WithinGridWaterDist.txt') THEN
      IF (any(headerwgwaterdist_file /= headerwgwaterdist_reqd)) THEN
         WRITE (*, *) headerwgwaterdist_file == headerwgwaterdist_reqd
         CALL errorhint(56, 'Names or order of columns in SUEWS_WithinGridWaterDist.txt does not match model code.', &
                        notused, notused, notusedi)
      END IF
 
   ELSEIF (filename == 'SUEWS_BiogenCO2.txt') THEN
      IF (any(headerbiogen_file /= headerbiogen_reqd)) THEN
         WRITE (*, *) headerbiogen_file == headerbiogen_reqd
         CALL errorhint(56, 'Names or order of columns in SUEWS_BiogenCO2.txt does not match model code.', &
                        notused, notused, notusedi)
      END IF
 
   ELSE
      WRITE (*, *) 'Problem in subroutine InputHeaderCheck. File header not specified in model code for ', filename
      CALL errorhint(58, filename, notused, notused, notusedi)
 
   END IF
 
END SUBROUTINE inputheadercheck
 
!-------------------------------------------------------------------------
!
! TS 05 Jul 2018: No longer needed as interpolation is done through specific subroutines at each required instant
!Interpolates hourly profiles provided in SUEWS_Profiles.txt
! to resolution of the model timestep
! HCW 06 Feb 2015
!===================================================================================
! SUBROUTINE SUEWS_InterpHourlyProfiles(Gridiv,TstepP_ID,SurfChar_HrProf)
!
!   USE allocateArray
!   USE ColNamesInputFiles
!   USE sues_data
!
!   IMPLICIT NONE
!
!   INTEGER:: i,j, ii   !Used to count over hours and sub-hourly timesteps
!   INTEGER:: Gridiv, TstepP_ID
!   INTEGER,DIMENSION(24):: SurfChar_HrProf
!   REAL(KIND(1d0)):: deltaProf   !Change in hourly profiles per model timestep
!
!   ! Copy value for first hour
!   TstepProfiles(Gridiv,TstepP_ID,1) = SurfaceChar(Gridiv,SurfChar_HrProf(1))
!   DO i=1,24
!      j = (i+1)
!      IF(i == 24) j = 1   !If last hour of day, loop round to first hour of day for interpolation
!      deltaProf = ((SurfaceChar(Gridiv,SurfChar_HrProf(j)) - SurfaceChar(Gridiv,SurfChar_HrProf(i))))/nsh_real
!      DO ii=1,nsh
!         IF((nsh*(i-1)+ii+1) < (23*nsh+nsh+1))  THEN
!            TstepProfiles(Gridiv,TstepP_ID,(nsh*(i-1)+ii+1)) = SurfaceChar(Gridiv,SurfChar_HrProf(i)) + deltaProf*ii
!         ENDIF
!      ENDDO
!   ENDDO
!
! endsubroutine SUEWS_InterpHourlyProfiles
!===================================================================================
 
!===================================================================================
! get interpolated profile values at specified time
! NO normalisation performed
FUNCTION get_prof_spectime_inst(Hour, Min, Sec, Prof_24h) RESULT(Prof_CurrTime)
 
   IMPLICIT NONE
 
   INTEGER :: i, j !Used to count over hours and sub-hourly timesteps
   INTEGER, INTENT(IN) :: hour, min, sec
   INTEGER :: total_sec, secperhour
   REAL(kind(1d0)), DIMENSION(0:23), INTENT(IN) :: prof_24h
   REAL(kind(1d0)) :: deltaprof !Change in hourly profiles per model timestep
   REAL(kind(1d0)) :: prof_currtime
 
   total_sec = min*60 + sec
   secperhour = 3600
 
   i = hour
   j = i + 1
   IF (j == 24) j = 0
 
   deltaprof = (prof_24h(j) - prof_24h(i))/secperhour
   prof_currtime = prof_24h(hour) + deltaprof*total_sec
 
END FUNCTION get_prof_spectime_inst
 
!===================================================================================
! get interpolated profile values at specified time
! normalise so the AVERAGE of the multipliers is equal to 1
FUNCTION get_prof_spectime_mean(Hour, Min, Sec, Prof_24h) RESULT(Prof_CurrTime)
 
   IMPLICIT NONE
 
   INTEGER :: i, j !Used to count over hours and sub-hourly timesteps
   INTEGER, INTENT(IN) :: hour, min, sec
   INTEGER :: total_sec, secperhour
   REAL(kind(1d0)), DIMENSION(0:23), INTENT(IN) :: prof_24h
   REAL(kind(1d0)), DIMENSION(0:23) :: prof_24h_mean
   REAL(kind(1d0)) :: deltaprof !Change in hourly profiles per model timestep
   REAL(kind(1d0)) :: prof_currtime
 
   total_sec = min*60 + sec
   secperhour = 3600
 
   prof_24h_mean = merge(prof_24h/(sum(prof_24h)/SIZE(prof_24h, dim=1)), 0.d0, sum(prof_24h) /= 0) ! prevent zero-division
   ! print*, Prof_24h_mean
 
   i = hour
   j = i + 1
   IF (j == 24) j = 0
 
   deltaprof = (prof_24h_mean(j) - prof_24h_mean(i))/secperhour
 
   ! print*, deltaProf,total_sec
   prof_currtime = prof_24h_mean(i) + deltaprof*total_sec
 
END FUNCTION get_prof_spectime_mean
 
!===================================================================================
! get interpolated profile values at specified time
! normalise so the SUM of the multipliers is equal to 1
FUNCTION get_prof_spectime_sum(Hour, Min, Sec, Prof_24h, dt) RESULT(Prof_CurrTime)
 
   IMPLICIT NONE
 
   INTEGER :: i, j !Used to count over hours and sub-hourly timesteps
   INTEGER, INTENT(IN) :: hour, min, sec, dt
   INTEGER :: total_sec, secperhour
   REAL(kind(1d0)), DIMENSION(0:23), INTENT(IN) :: prof_24h
   REAL(kind(1d0)), DIMENSION(0:23) :: prof_24h_sum
   REAL(kind(1d0)) :: deltaprof !Change in hourly profiles per model timestep
   REAL(kind(1d0)) :: prof_currtime
 
   total_sec = min*60 + sec
   secperhour = 3600
 
   prof_24h_sum = merge(prof_24h/(sum(prof_24h)), 0.d0, sum(prof_24h) /= 0) ! prevent zero-division
 
   i = hour
   j = i + 1
   IF (j == 24) j = 0
 
   deltaprof = (prof_24h_sum(j) - prof_24h_sum(i))/secperhour
   prof_currtime = prof_24h_sum(hour) + deltaprof*total_sec
   prof_currtime = prof_currtime*dt/secperhour
 
END FUNCTION get_prof_spectime_sum
!===================================================================================
! TODO: #79 all the parameter tables can be stored using a custom dictionary for easier and more consistent access:
! this package looks promising: https://github.com/zerothi/fdict
 
! Subroutines for matching codes in the input files
!  could re-write as a generic function later...
 
SUBROUTINE codematchohm(Gridiv, is, SWWD)
   ! Matches OHM coefficients a1, a2, a3 via OHM codes
   ! for summer/winter wet/dry conditions
   ! HCW 03 Nov 2014
   ! ---------------------------------------------------------
 
   USE allocatearray
   USE initial
   USE colnamesinputfiles
   USE defaultnotused
 
   IMPLICIT NONE
 
   INTEGER :: gridiv
   INTEGER :: is
   CHARACTER(len=4) :: SWWD
 
   iv5 = 0 ! Reset iv5 to zero
 
   IF (swwd == 'SWet') THEN
 
      DO iv5 = 1, nlinesohmcoefficients
         IF (ohmcoefficients_coeff(iv5, co_code) == surfacechar(gridiv, c_ohmcode_swet(is))) THEN
            EXIT
         ELSEIF (iv5 == nlinesohmcoefficients) THEN
            WRITE (*, *) 'Program stopped! OHM code (summer wet)', surfacechar(gridiv, c_ohmcode_swet(is)), &
               'not found in OHM_Coefficients.txt for surface', is, '.'
            CALL errorhint(57, 'Cannot find OHM code (summer wet)', surfacechar(gridiv, c_ohmcode_swet(is)), notused, notusedi)
         END IF
      END DO
 
   ELSEIF (swwd == 'SDry') THEN
 
      DO iv5 = 1, nlinesohmcoefficients
         IF (ohmcoefficients_coeff(iv5, co_code) == surfacechar(gridiv, c_ohmcode_sdry(is))) THEN
            EXIT
         ELSEIF (iv5 == nlinesohmcoefficients) THEN
            WRITE (*, *) 'Program stopped! OHM code (summer dry)', surfacechar(gridiv, c_ohmcode_sdry(is)), &
               'not found in OHM_Coefficients.txt for surface', is, '.'
            CALL errorhint(57, 'Cannot find OHM code (summer dry)', surfacechar(gridiv, c_ohmcode_sdry(is)), notused, notusedi)
         END IF
      END DO
 
   ELSEIF (swwd == 'WWet') THEN
 
      DO iv5 = 1, nlinesohmcoefficients
         IF (ohmcoefficients_coeff(iv5, co_code) == surfacechar(gridiv, c_ohmcode_wwet(is))) THEN
            EXIT
         ELSEIF (iv5 == nlinesohmcoefficients) THEN
            WRITE (*, *) 'Program stopped! OHM code (winter wet)', surfacechar(gridiv, c_ohmcode_wwet(is)), &
               'not found in OHM_Coefficients.txt for surface', is, '.'
            CALL errorhint(57, 'Cannot find OHM code (winter wet)', surfacechar(gridiv, c_ohmcode_wwet(is)), notused, notusedi)
         END IF
      END DO
 
   ELSEIF (swwd == 'WDry') THEN
 
      DO iv5 = 1, nlinesohmcoefficients
         IF (ohmcoefficients_coeff(iv5, co_code) == surfacechar(gridiv, c_ohmcode_wdry(is))) THEN
            EXIT
         ELSEIF (iv5 == nlinesohmcoefficients) THEN
            WRITE (*, *) 'Program stopped! OHM code (winter dry)', surfacechar(gridiv, c_ohmcode_wdry(is)), &
               'not found in OHM_Coefficients.txt for surface', is, '.'
            CALL errorhint(57, 'Cannot find OHM code (winter dry)', surfacechar(gridiv, c_ohmcode_wdry(is)), notused, notusedi)
         END IF
      END DO
 
   ELSE
      WRITE (*, *) 'Problem with CodeMatchOHM (in SUEWS_CodeMatch.f95). ', swwd, ' not recognised. Needs to be one of: ', &
         'SWet = Summer Wet, SDry = Summer Dry, WWet = WinterWet, WDry = Winter Dry. N.B. Case sensitive. '
      stop
   END IF
 
   RETURN
END SUBROUTINE codematchohm
! ---------------------------------------------------------
 
SUBROUTINE codematchestm(Gridiv, is)
   ! Matches ESTM coefficients via ESTM code
   ! Modified HCW 16 Jun 2016 - for SUEWS surface types
   !                          - removed summer/winter wet/dry option
   ! S.O. 04 Feb 2016
   ! ---------------------------------------------------------
 
   USE allocatearray
   USE initial
   USE colnamesinputfiles
   USE defaultnotused
 
   IMPLICIT NONE
 
   INTEGER :: gridiv
   INTEGER :: is
 
   iv5 = 0 ! Reset iv5 to zero
 
   DO iv5 = 1, nlinesestmcoefficients
      IF (estmcoefficients_coeff(iv5, ce_code) == surfacechar(gridiv, c_estmcode(is))) THEN
         EXIT
      ELSEIF (iv5 == nlinesestmcoefficients) THEN
         WRITE (*, *) 'Program stopped! ESTM code', surfacechar(gridiv, c_estmcode(is)), &
            'not found in ESTM_Coefficients.txt for surface', is, '.'
         CALL errorhint(57, 'Cannot find ESTM code', surfacechar(gridiv, c_estmcode(is)), notused, notusedi)
      END IF
   END DO
 
   RETURN
END SUBROUTINE codematchestm
! ---------------------------------------------------------
 
SUBROUTINE codematchestm_class(Gridiv, is, ii)
   ! Matches ESTM coefficients via ESTM codes in SiteSelect for Paved and Bldgs ESTM classes
   ! HCW 16 Jun 2016
   ! ---------------------------------------------------------
 
   USE allocatearray
   USE initial
   USE colnamesinputfiles
   USE defaultnotused
 
   IMPLICIT NONE
 
   INTEGER :: gridiv
   INTEGER :: is, ii
 
   iv5 = 0 ! Reset iv5 to zero
 
   IF (is == bldgsurf) THEN
      DO iv5 = 1, nlinesestmcoefficients
         IF (estmcoefficients_coeff(iv5, ce_code) == surfacechar(gridiv, c_code_estmclass_bldgs(ii))) THEN
            EXIT
         ELSEIF (iv5 == nlinesestmcoefficients) THEN
            WRITE (*, *) 'Program stopped! ESTM code', surfacechar(gridiv, c_code_estmclass_bldgs(ii)), &
               'not found in ESTM_Coefficients.txt for surface', is, '.'
            CALL errorhint(57, 'Cannot find ESTM code', surfacechar(gridiv, c_code_estmclass_bldgs(ii)), notused, notusedi)
         END IF
      END DO
   ELSEIF (is == pavsurf) THEN
      DO iv5 = 1, nlinesestmcoefficients
         IF (estmcoefficients_coeff(iv5, ce_code) == surfacechar(gridiv, c_code_estmclass_paved(ii))) THEN
            EXIT
         ELSEIF (iv5 == nlinesestmcoefficients) THEN
            WRITE (*, *) 'Program stopped! ESTM code', surfacechar(gridiv, c_code_estmclass_paved(ii)), &
               'not found in ESTM_Coefficients.txt for surface', is, '.'
            CALL errorhint(57, 'Cannot find ESTM code', surfacechar(gridiv, c_code_estmclass_paved(ii)), notused, notusedi)
         END IF
      END DO
   ELSE
      WRITE (*, *) 'Problem with CodeMatchESTM_Class (in SUEWS_ctrl_input.f95). ', is, ' not correct. Needs to be either ', &
         '1 = Paved surfaced, 2 = Bldgs surfaces.'
      stop
   END IF
   RETURN
END SUBROUTINE codematchestm_class
! ---------------------------------------------------------
 
SUBROUTINE codematchprof(Gridiv, SurfaceCharCodeCol)
   ! Matches Soil characteristics via codes in *SurfaceChar*
   ! for energy use/water use/snow clearing
   ! HCW 20 Nov 2014
   ! ---------------------------------------------------------
 
   USE allocatearray
   USE initial
   USE colnamesinputfiles
   USE defaultnotused
 
   IMPLICIT NONE
 
   INTEGER :: Gridiv
   INTEGER :: SurfaceCharCodeCol
 
   iv5 = 0 ! Reset iv5 to zero
 
   DO iv5 = 1, nlinesprofiles
      IF (profiles_coeff(iv5, cpr_code) == surfacechar(gridiv, surfacecharcodecol)) THEN
         EXIT
      ELSEIF (iv5 == nlinesprofiles) THEN
         WRITE (*, *) 'Program stopped! Profile code ', surfacechar(gridiv, surfacecharcodecol), 'not found in SUEWS_Profiles.txt.'
         CALL errorhint(57, 'Cannot find code in SUEWS_Profiles.txt', surfacechar(gridiv, surfacecharcodecol), notused, notusedi)
      END IF
   END DO
 
   RETURN
END SUBROUTINE codematchprof
! ---------------------------------------------------------
 
SUBROUTINE codematchdist(rr, CodeCol, codeColSameSurf)
   ! Matches within-grid water distribution via codes
   ! Checks water cannot flow from one surface to the same surface
   ! Checks water distribution fractions sum to 1.
   ! HCW 10 Nov 2014
   ! ---------------------------------------------------------
 
   USE allocatearray
   USE initial
   USE colnamesinputfiles
   USE defaultnotused
 
   IMPLICIT NONE
 
   INTEGER :: rr
   INTEGER :: codeCol, codeColSameSurf
 
   iv5 = 0 ! Reset iv5 to zero
 
   DO iv5 = 1, nlineswgwaterdist
      IF (wgwaterdist_coeff(iv5, cwg_code) == siteselect(rr, codecol)) THEN
         EXIT
      ELSEIF (iv5 == nlineswgwaterdist) THEN
         WRITE (*, *) 'Program stopped! Within-grid water distribution code ', siteselect(rr, codecol), &
            'not found in SUEWS_WaterDistWithinGrid.txt.'
         CALL errorhint(57, 'Cannot find code in SUEWS_WaterDistWithinGrid.txt', siteselect(rr, codecol), notused, notusedi)
      END IF
   END DO
 
   ! Check water flow to same surface is zero (previously in RunControlByGridByYear in SUEWS_Initial.f95)
   IF (wgwaterdist_coeff(iv5, codecolsamesurf) /= 0) THEN
      CALL errorhint(8, 'Diagonal elements should be zero as water cannot move from one surface to the same surface.', &
                     wgwaterdist_coeff(iv5, codecolsamesurf), notused, notusedi)
   END IF
 
   !! QUESTION: MODIFY THIS?
   ! Check water either moves to runoff or soilstore, but not to both
   ! Model returns an error if both ToRunoff and ToSoilStore are non-zero.
   !! - Probably should remove this...
   !! - Also look at SUEWS_translate, as the non-zero value goes into WaterDist
   IF (wgwaterdist_coeff(iv5, cwg_torunoff) /= 0 .AND. wgwaterdist_coeff(iv5, cwg_tosoilstore) /= 0) THEN
      CALL errorhint(9, 'One of these (ToRunoff,ToSoilStore) should be zero.', &
                     wgwaterdist_coeff(iv5, cwg_torunoff), wgwaterdist_coeff(iv5, cwg_tosoilstore), notusedi)
   END IF
 
   !! Also do for water surface once implemented
   IF (codecol /= c_wgwatercode) THEN ! Except for Water surface
      ! Check total water distribution from each surface adds up to 1
      IF (sum(wgwaterdist_coeff(iv5, cwg_topaved:cwg_tosoilstore)) > 1.0000001 &
          .OR. sum(wgwaterdist_coeff(iv5, cwg_topaved:cwg_tosoilstore)) < 0.9999999) THEN
         CALL errorhint(8, 'Total water distribution from each surface should add up to 1.', &
                        sum(wgwaterdist_coeff(iv5, cwg_topaved:cwg_tosoilstore)), notused, notusedi)
      END IF
   END IF
 
   RETURN
END SUBROUTINE codematchdist
! ---------------------------------------------------------
 
SUBROUTINE codematchnonveg(rr, CodeCol)
   ! Matches Impervious characteristics via codes in SiteSelect
   ! HCW 20 Nov 2014
   ! ---------------------------------------------------------
 
   USE allocatearray
   USE initial
   USE colnamesinputfiles
   USE defaultnotused
 
   IMPLICIT NONE
 
   INTEGER :: rr
   INTEGER :: codeCol
 
   iv5 = 0 ! Reset iv5 to zero
 
   DO iv5 = 1, nlinesnonveg
      IF (nonveg_coeff(iv5, ci_code) == siteselect(rr, codecol)) THEN
         EXIT
      ELSEIF (iv5 == nlinesnonveg) THEN
         WRITE (*, *) 'Program stopped! NonVeg code ', siteselect(rr, codecol), 'not found in SUEWS_NonVeg.txt.'
         CALL errorhint(57, 'Cannot find code in SUEWS_NonVeg.txt', siteselect(rr, codecol), notused, notusedi)
      END IF
   END DO
 
   RETURN
END SUBROUTINE codematchnonveg
! ---------------------------------------------------------
 
SUBROUTINE codematchveg(rr, CodeCol)
   ! Matches Pervious characteristics via codes in SiteSelect
   ! HCW 20 Nov 2014
   ! ---------------------------------------------------------
 
   USE allocatearray
   USE initial
   USE colnamesinputfiles
   USE defaultnotused
 
   IMPLICIT NONE
 
   INTEGER :: rr
   INTEGER :: codeCol
 
   iv5 = 0 ! Reset iv5 to zero
 
   DO iv5 = 1, nlinesveg
      IF (veg_coeff(iv5, cp_code) == siteselect(rr, codecol)) THEN
         EXIT
      ELSEIF (iv5 == nlinesveg) THEN
         WRITE (*, *) 'Program stopped! Veg code ', siteselect(rr, codecol), 'not found in SUEWS_Vegs.txt.'
         CALL errorhint(57, 'Cannot find code in SUEWS_Veg.txt', siteselect(rr, codecol), notused, notusedi)
      END IF
   END DO
 
   RETURN
END SUBROUTINE codematchveg
! ---------------------------------------------------------
 
SUBROUTINE codematchwater(rr, CodeCol)
   ! Matches Water characteristics via codes in SiteSelect
   ! HCW 20 Nov 2014
   ! ---------------------------------------------------------
 
   USE allocatearray
   USE initial
   USE colnamesinputfiles
   USE defaultnotused
 
   IMPLICIT NONE
 
   INTEGER :: rr
   INTEGER :: codeCol
 
   iv5 = 0 ! Reset iv5 to zero
 
   DO iv5 = 1, nlineswater
      IF (water_coeff(iv5, cw_code) == siteselect(rr, codecol)) THEN
         EXIT
      ELSEIF (iv5 == nlineswater) THEN
         WRITE (*, *) 'Program stopped! Water code ', siteselect(rr, codecol), 'not found in SUEWS_Water.txt.'
         CALL errorhint(57, 'Cannot find code in SUEWS_Water.txt', siteselect(rr, codecol), notused, notusedi)
      END IF
   END DO
 
   RETURN
END SUBROUTINE codematchwater
! ---------------------------------------------------------
 
SUBROUTINE codematchsnow(rr, CodeCol)
   ! Matches Snow characteristics via codes in SiteSelect
   ! HCW 20 Nov 2014
   ! ---------------------------------------------------------
 
   USE allocatearray
   USE initial
   USE colnamesinputfiles
   USE defaultnotused
 
   IMPLICIT NONE
 
   INTEGER :: rr
   INTEGER :: codeCol
 
   iv5 = 0 ! Reset iv5 to zero
 
   DO iv5 = 1, nlinessnow
      IF (snow_coeff(iv5, cs_code) == siteselect(rr, codecol)) THEN
         EXIT
      ELSEIF (iv5 == nlinessnow) THEN
         WRITE (*, *) 'Program stopped! Snow code ', siteselect(rr, codecol), 'not found in SUEWS_Snow.txt.'
         CALL errorhint(57, 'Cannot find code in SUEWS_Snow.txt', siteselect(rr, codecol), notused, notusedi)
      END IF
   END DO
 
   RETURN
END SUBROUTINE codematchsnow
! ---------------------------------------------------------
 
SUBROUTINE codematchconductance(rr, CodeCol)
   ! Matches Conductance characteristics via codes in SiteSelect
   ! HCW 20 Nov 2014
   ! ---------------------------------------------------------
 
   USE allocatearray
   USE initial
   USE colnamesinputfiles
   USE defaultnotused
 
   IMPLICIT NONE
 
   INTEGER :: rr
   INTEGER :: codeCol
 
   iv5 = 0 ! Reset iv5 to zero
 
   DO iv5 = 1, nlinesconductance
      IF (conductance_coeff(iv5, cc_code) == siteselect(rr, codecol)) THEN
         EXIT
      ELSEIF (iv5 == nlinesconductance) THEN
         WRITE (*, *) 'Program stopped! Conductance code ', siteselect(rr, codecol), 'not found in SUEWS_Conductance.txt.'
         CALL errorhint(57, 'Cannot find code in SUEWS_Conductance.txt', siteselect(rr, codecol), notused, notusedi)
      END IF
   END DO
 
   RETURN
END SUBROUTINE codematchconductance
! ---------------------------------------------------------
 
SUBROUTINE codematchanthropogenic(rr, CodeCol)
   ! Matches AnthropogenicHeat characteristics via codes in SiteSelect
   ! HCW 20 Nov 2014
   ! MH 21 Jun 2017
   ! ---------------------------------------------------------
 
   USE allocatearray
   USE initial
   USE colnamesinputfiles
   USE defaultnotused
 
   IMPLICIT NONE
 
   INTEGER :: rr
   INTEGER :: codeCol
 
   iv5 = 0 ! Reset iv5 to zero
 
   DO iv5 = 1, nlinesanthropogenic
      IF (anthropogenic_coeff(iv5, ca_code) == siteselect(rr, codecol)) THEN
         EXIT
      ELSEIF (iv5 == nlinesanthropogenic) THEN
         WRITE (*, *) 'Program stopped! Anthropogenic code ', siteselect(rr, codecol), &
            'not found in SUEWS_AnthropogenicEmission.txt.'
         CALL errorhint(57, 'Cannot find code in SUEWS_AnthropogenicEmission.txt', &
                        siteselect(rr, codecol), notused, notusedi)
      END IF
   END DO
 
   RETURN
END SUBROUTINE codematchanthropogenic
! ---------------------------------------------------------
 
SUBROUTINE codematchirrigation(rr, CodeCol)
   ! Matches Irrigation characteristics via codes in SiteSelect
   ! HCW 20 Nov 2014
   ! ---------------------------------------------------------
 
   USE allocatearray
   USE initial
   USE colnamesinputfiles
   USE defaultnotused
 
   IMPLICIT NONE
 
   INTEGER :: rr
   INTEGER :: codeCol
 
   iv5 = 0 ! Reset iv5 to zero
 
   DO iv5 = 1, nlinesirrigation
      IF (irrigation_coeff(iv5, cir_code) == siteselect(rr, codecol)) THEN
         EXIT
      ELSEIF (iv5 == nlinesirrigation) THEN
         WRITE (*, *) 'Program stopped! Irrigation code ', siteselect(rr, codecol), 'not found in SUEWS_Irrigation.txt.'
         CALL errorhint(57, 'Cannot find code in SUEWS_Irrigation.txt', siteselect(rr, codecol), notused, notusedi)
      END IF
   END DO
 
   RETURN
END SUBROUTINE codematchirrigation
! ---------------------------------------------------------
 
SUBROUTINE codematchsoil(Gridiv, SurfaceCharCodeCol)
   ! Matches Soil characteristics via codes in *SurfaceChar*
   ! HCW 20 Nov 2014
   ! ---------------------------------------------------------
 
   USE allocatearray
   USE initial
   USE colnamesinputfiles
   USE defaultnotused
 
   IMPLICIT NONE
 
   INTEGER :: Gridiv
   INTEGER :: SurfaceCharCodeCol
 
   iv5 = 0 ! Reset iv5 to zero
 
   DO iv5 = 1, nlinessoil
      IF (soil_coeff(iv5, cso_code) == surfacechar(gridiv, surfacecharcodecol)) THEN
         EXIT
      ELSEIF (iv5 == nlinessoil) THEN
         WRITE (*, *) 'Program stopped! Soil code ', surfacechar(gridiv, surfacecharcodecol), 'not found in SUEWS_Soil.txt.'
         CALL errorhint(57, 'Cannot find code in SUEWS_Soil.txt', surfacechar(gridiv, surfacecharcodecol), notused, notusedi)
      END IF
   END DO
 
   RETURN
END SUBROUTINE codematchsoil
! ---------------------------------------------------------
 
SUBROUTINE codematchbiogen(Gridiv, SurfaceCharCodeCol)
   ! Matches Biogen characteristics via codes in *SuraceChar*
   ! MH 16 Jun 2017
   ! ---------------------------------------------------------
 
   USE allocatearray
   USE initial
   USE colnamesinputfiles
   USE defaultnotused
 
   IMPLICIT NONE
 
   INTEGER :: Gridiv
   INTEGER :: SurfaceCharCodeCol
 
   iv5 = 0 ! Reset iv5 to zero
 
   DO iv5 = 1, nlinesbiogen
      IF (biogen_coeff(iv5, cb_code) == surfacechar(gridiv, surfacecharcodecol)) THEN
         EXIT
      ELSEIF (iv5 == nlinesbiogen) THEN
         WRITE (*, *) 'Program stopped! Biogen code ', surfacechar(gridiv, surfacecharcodecol), 'not found in SUEWS_BiogenCO2.txt.'
         CALL errorhint(57, 'Cannot find code in SUEWS_BiogenCO2.txt', surfacechar(gridiv, surfacecharcodecol), notused, notusedi)
      END IF
   END DO
 
   RETURN
END SUBROUTINE codematchbiogen
! ---------------------------------------------------------
 
MODULE metdisagg
   !========================================================================================
   ! Disaggregation of meteorological forcing data
   !  Code to disaggregate met forcing data from resolution provided to the model time-step
   ! Created: HCW 10 Feb 2017
   !
   ! ---- Key for MetDisaggMethod settings ----
   !10  ->  linear disaggregation for timestamps representing the end of the averaging period
   !20  ->  linear disaggregation for instantaneous variables (e.g. air temperature, humidity, pressure, wind speed in WFDEI dataset)
   !100 -> evenly distribute rainfall among all subintervals in a rainy interval
   !101 -> evenly distribute rainfall among RainAmongN subintervals in a rainy interval
   !        - requires RainAmongN to be set in RunControl.nml
   ! If KdownZen = 1 -> include additional zenith check in kdown disaggregation
   !
   ! N.B. wdir downscaling is currently not implemented
   !
   !========================================================================================
 
   USE allocatearray
   USE colnamesinputfiles
   USE data_in
   USE initial
   USE narp_module, ONLY: narp_cal_sunposition
 
   IMPLICIT NONE
 
CONTAINS
 
   !======================================================================================
   SUBROUTINE disaggregatemet(iBlock, igrid)
      ! Subroutine to disaggregate met forcing data to model time-step
      ! HCW 10 Feb 2017
      !======================================================================================
 
      USE sues_data
      USE defaultnotused
 
      IMPLICIT NONE
 
      INTEGER :: lunit = 100
      INTEGER :: tdiff !Time difference (in minutes) between first and second rows of original met forcing file
      INTEGER :: i, ii !counter
      INTEGER :: iBlock, igrid
      INTEGER, DIMENSION(NperTstepIn) :: seq1Nper
      INTEGER, DIMENSION(nsd) :: seq1nsd
      INTEGER, DIMENSION(nColumnsMetForcingData) :: MetDisaggMethod ! Stores method to use for disaggregating met data
      REAL(KIND(1D0)), DIMENSION(nColumnsMetForcingData) :: MetArrayOrig
      REAL(KIND(1D0)), DIMENSION(ReadLinesOrigMetData*NperTstepIn, ncolumnsMetForcingData) :: Met_tt
      REAL(KIND(1D0)), DIMENSION(ReadLinesOrigMetData*NperTstepIn) :: Met_tt_kdownAdj
      CHARACTER(LEN=9), DIMENSION(ncolumnsMetForcingData) :: HeaderMet
      CHARACTER(LEN=10*ncolumnsMetForcingData) :: HeaderMetOut
      ! REAL(KIND(1d0)),DIMENSION(ReadLinesOrigMetData):: dectimeOrig
      ! REAL(KIND(1d0)),DIMENSION(ReadLinesOrigMetData*Nper):: dectimeDscd, dectimeFast
      REAL(KIND(1D0)), DIMENSION(ReadLinesOrigMetData*NperTstepIn) :: dectimeFast
      REAL(KIND(1D0)), DIMENSION(ReadLinesOrigMetData*NperTstepIn) :: idectime ! sun position at middle of time-step before
 
      ! INTEGER, DIMENSION(Nper):: temp_iy, temp_id, temp_ih, temp_im, temp_ihm ! temorary varaibles for disaggragation
      ! REAL(KIND(1d0)), DIMENSION(Nper):: temp_dectime ! temorary varaibles for disaggragation
      ! INTEGER :: Days_of_Year
      ! INTEGER, DIMENSION(Nper)::ndays_iy ! number of days in iy
 
      ! Allocate and initialise arrays to receive original forcing data --------------------
      ALLOCATE (metfordisagg(readlinesorigmetdata, ncolumnsmetforcingdata))
      ALLOCATE (metfordisaggprev(ncolumnsmetforcingdata))
      ALLOCATE (metfordisaggnext(ncolumnsmetforcingdata))
      metfordisagg(:, :) = -999
      metfordisaggprev(:) = -999
      metfordisaggnext(:) = -999
      ! Initialise array to receive disaggregated data
      met_tt = -999
      ! Intialise disaggregation method
      metdisaggmethod = -999
 
      ! Generate useful sequences
      seq1nper = (/(i, i=1, npertstepin, 1)/)
      seq1nsd = (/(i, i=1, nsd, 1)/)
 
      ! Get methods to use for disaggregation from RunControl
      IF (diagnosedisagg == 1) WRITE (*, *) 'DisaggMethod: ', disaggmethod, 'RainDisaggMethod:', raindisaggmethod
      IF (disaggmethod == 1) THEN
         metdisaggmethod(:) = 10 !linear disaggregation of averages
      ELSEIF (disaggmethod == 2) THEN
         metdisaggmethod(:) = 20 !linear disaggregation of instantaneous values
      ELSEIF (disaggmethod == 3) THEN !WFDEI set up, where T, Q, pres, U are instantaneous
         metdisaggmethod(:) = 10 !linear disaggregation of averages
         metdisaggmethod(10:13) = 20 !linear disagg instantaneous values for U, RH, Tair, pres
      ELSE
         CALL errorhint(2, 'Problem in SUEWS_MetDisagg: DisaggMethod value should be 1, 2, or 3', &
                        notused, notused, disaggmethod)
      END IF
      ! Set rainfall
      metdisaggmethod(14) = raindisaggmethod
 
      ! Read data ---------------------------------------------------------------------
      IF (diagnosedisagg == 1) WRITE (*, *) 'Reading file: ', trim(fileorigmet)
      OPEN (lunit, file=trim(fileorigmet), status='old')
      ! CALL skipHeader(lunit,SkipHeaderMet)  !Skip header -> read header instead
      READ (lunit, *) headermet
      !write(*,*) HeaderMet
      ! Skip over lines that have already been read and downscaled
      IF (skippedlinesorig > 0) THEN
         DO i = 1, skippedlinesorig - 1 ! minus 1 here because last line of last block needs to be read again
            READ (lunit, *)
         END DO
         ! Read in last line of previous block
         CALL metread(lunit, metarrayorig, inputmetformat, ldown_option, netradiationmethod, &
                      snowuse, smdmethod, soildepthmeas, soilrocks, soildensity, smcap)
         metfordisaggprev(1:ncolumnsmetforcingdata) = metarrayorig
      END IF
      ! print*, 'MetForDisagg',MetForDisagg(1:3,1:4)
      ! print*, 'ReadLinesOrigMetDataMax',ReadLinesOrigMetDataMax
      ! Read in current block
      DO i = 1, readlinesorigmetdatamax
         CALL metread(lunit, metarrayorig, inputmetformat, ldown_option, netradiationmethod, &
                      snowuse, smdmethod, soildepthmeas, soilrocks, soildensity, smcap)
         metfordisagg(i, 1:ncolumnsmetforcingdata) = metarrayorig
      END DO
      ! print*, 'MetForDisagg',MetForDisagg(1:3,1:4)
      ! Read in first line of next block (except for last block)
      IF (iblock /= readblocksorigmetdata) THEN
         CALL metread(lunit, metarrayorig, inputmetformat, ldown_option, netradiationmethod, &
                      snowuse, smdmethod, soildepthmeas, soilrocks, soildensity, smcap)
         metfordisaggnext(1:ncolumnsmetforcingdata) = metarrayorig
      END IF
      CLOSE (lunit)
 
      ! Check resolution of original met forcing data -------------------------------------
      ! Find time difference (in minutes) between first and second row
      tdiff = int(metfordisagg(2, 4) - metfordisagg(1, 4)) !Try using minutes
      IF (tdiff == 0) tdiff = int(metfordisagg(2, 3) - metfordisagg(1, 3))*60 !If no difference in minutes, try using hours
      IF (tdiff < 0) THEN !If time difference is negative (e.g. change of day), instead use second and third row
         tdiff = int(metfordisagg(3, 4) - metfordisagg(2, 4))
         IF (tdiff == 0) tdiff = int(metfordisagg(3, 3) - metfordisagg(2, 3))*60 !If no difference in minutes, try using hours
      END IF
      ! Check actual resolution matches specified input resolution
      IF (tdiff /= resolutionfilesin/60) THEN
         CALL errorhint(2, 'Problem in SUEWS_MetDisagg: timestamps in met forcing file inconsistent with ResolutionFilesIn', &
                        REAL(ResolutionFilesIn, KIND(1D0)), NotUsed, tdiff*60)
      END IF
 
      ! Check file only contains a single year --------------------------------------------
      ! Very last data point is allowed to be (should be) timestamped with following year
      IF (any(metfordisagg(1:(readlinesorigmetdatamax - 1), 1) /= metfordisagg(1, 1))) THEN
         CALL errorhint(3, 'Problem in SUEWS_MetDisagg: multiple years found in original met forcing file.', &
                        metfordisagg(1, 1), notused, notusedi)
      END IF
 
      ! Disaggregate time columns ---------------------------------------------------------
      IF (diagnose == 1) WRITE (*, *) 'Disaggregating met forcing data (', trim(fileorigmet), ') to model time-step...'
 
      CALL disaggregatedatetime(metfordisagg(:, 1:4), tstep, npertstepin, readlinesorigmetdatamax, met_tt(:, 1:4))
 
      ! Disaggregate other columns --------------------------------------------------------
      DO ii = 5, ncolumnsmetforcingdata
         IF (ii == 14) THEN !Do something different for rainfall and snowfall (if present)
            IF (metdisaggmethod(14) == 100) THEN
               met_tt(:, 14) = disaggp_amongn(metfordisagg(:, 14), &
                                              npertstepin, npertstepin, readlinesorigmetdata, readlinesorigmetdatamax)
               IF (all(metfordisagg(:, 16) == -999)) THEN
                  met_tt(:, 16) = -999
               ELSE
                  met_tt(:, 16) = disaggp_amongn(metfordisagg(:, 16), &
                                                 npertstepin, npertstepin, readlinesorigmetdata, readlinesorigmetdatamax)
               END IF
            ELSEIF (metdisaggmethod(14) == 101) THEN
               IF (rainamongn == -999) THEN
                  CALL errorhint(2, 'Problem in SUEWS_MetDisagg: RainDisaggMethod requires RainAmongN', &
                                 REAL(RainAmongN, KIND(1D0)), NotUsed, RainDisaggMethod)
               ELSEIF (rainamongn > npertstepin) THEN
                  CALL errorhint(2, 'Problem in SUEWS_MetDisagg: RainAmongN > Nper', &
                                 REAL(NperTstepIn, KIND(1D0)), NotUsed, RainAmongN)
               ELSE
                  met_tt(:, 14) = disaggp_amongn(metfordisagg(:, 14), &
                                                 rainamongn, npertstepin, readlinesorigmetdata, readlinesorigmetdatamax)
                  IF (all(metfordisagg(:, 16) == -999)) THEN
                     met_tt(:, 16) = -999
                  ELSE
                     met_tt(:, 16) = disaggp_amongn(metfordisagg(:, 16), &
                                                    rainamongn, npertstepin, readlinesorigmetdata, readlinesorigmetdatamax)
                  END IF
               END IF
            ELSEIF (metdisaggmethod(14) == 102) THEN
               IF (all(multrainamongn == -999)) THEN
                  CALL errorhint(2, 'Problem in SUEWS_MetDisagg: RainDisaggMethod requires MultRainAmongN', &
                                 REAL(MultRainAmongN(1), KIND(1D0)), NotUsed, RainDisaggMethod)
               ELSEIF (all(multrainamongnupperi == -999)) THEN
                  CALL errorhint(2, 'Problem in SUEWS_MetDisagg: RainDisaggMethod requires MultRainAmongNUpperI', &
                                 multrainamongnupperi(1), notused, raindisaggmethod)
               ELSEIF (any(multrainamongn > npertstepin)) THEN
                  CALL errorhint(2, 'Problem in SUEWS_MetDisagg: MultRainAmongN > Nper', real(npertstepin, kind(1d0)), notused, &
                                 maxval(multrainamongn))
               ELSE
                  met_tt(:, 14) = disaggp_amongnmult(metfordisagg(:, 14), multrainamongnupperi, multrainamongn, npertstepin, &
                                                     readlinesorigmetdata, readlinesorigmetdatamax)
                  IF (all(metfordisagg(:, 16) == -999)) THEN
                     met_tt(:, 16) = -999
                  ELSE
                     met_tt(:, 16) = disaggp_amongnmult(metfordisagg(:, 16), multrainamongnupperi, multrainamongn, npertstepin, &
                                                        readlinesorigmetdata, readlinesorigmetdatamax)
                  END IF
               END IF
            ELSE
               WRITE (*, *) 'Disaggregation code for rain not recognised'
            END IF
         ELSEIF (ii == 24) THEN !wind direction disaggregation not coded yet...
            IF (any(metfordisagg(:, ii) /= -999)) THEN
               WRITE (*, *) 'Disaggregation of wind direction not currently implemented!'
            END IF
         ELSE
            IF (all(metfordisagg(:, ii) == -999)) THEN
               !IF(DiagnoseDisagg==1) write(*,*) 'No data for col.', ii
               met_tt(:, ii) = -999
            ELSE
               met_tt(:, ii) = disagg_lin(metfordisagg(:, ii), metfordisaggprev(ii), metfordisaggnext(ii), metdisaggmethod(ii), &
                                          npertstepin, readlinesorigmetdata, readlinesorigmetdatamax, iblock)
            END IF
         END IF
      END DO
 
      ! Adjust kdown disaggregation using zenith angle
      IF (kdownzen == 1) THEN
         IF (diagnosedisagg == 1) WRITE (*, *) 'Adjusting disaggregated kdown using zenith angle'
         met_tt_kdownadj(:) = met_tt(:, 15)
         ! Translate location data from SurfaceChar to find solar angles
         lat = surfacechar(igrid, c_lat)
         lng = surfacechar(igrid, c_lng)
         timezone = surfacechar(igrid, c_tz)
         alt = surfacechar(igrid, c_alt)
         ! Calculate dectime at downscaled time-step
         dectimefast(:) = met_tt(:, 2) + met_tt(:, 3)/24.0 + met_tt(:, 4)/(60.0*24.0)
         idectime = dectimefast - halftimestep ! sun position at middle of timestep before
         DO i = 1, (readlinesorigmetdatamax*npertstepin)
            CALL narp_cal_sunposition(met_tt(i, 2), idectime(i), timezone, lat, lng, alt, azimuth, zenith_deg)
            ! If sun below horizon, set disaggregated kdown to zero
            IF (zenith_deg > 90) THEN
               !write(*,*) Met_tt(i,1:4)
               met_tt_kdownadj(i) = 0.0
            END IF
         END DO
         ! Redistribute kdown over each day
         DO i = 1, (readlinesorigmetdatamax*npertstepin/nsd) ! Loop over each day
            met_tt_kdownadj((i - 1)*nsd + seq1nsd) = &
               met_tt_kdownadj((i - 1)*nsd + seq1nsd) &
               *sum(met_tt((i - 1)*nsd + seq1nsd, 15)) &
               /sum(met_tt_kdownadj((i - 1)*nsd + seq1nsd))
         END DO
         ! Copy adjusted kdown back to Met_tt array
         met_tt(:, 15) = met_tt_kdownadj(:)
      END IF
 
      ! Copy disaggregated data to MetForcingDataArray
      metforcingdata(:, 1:24, gridcounter) = met_tt(:, 1:24)
 
      ! If snow is -999, set to zero (also in LUMPS_metRead.f95)
      IF (all(metforcingdata(:, 16, gridcounter) == -999)) metforcingdata(:, 16, gridcounter) = 0
 
      ! Undo pressure conversion again for writing out
      met_tt(:, 13) = met_tt(:, 13)/10.0
 
      ! Write out disaggregated file ------------------------------------------------------
      IF (keeptstepfilesin == 1) THEN
         IF (iblock == 1) THEN
            ! Prepare header
            DO i = 1, ncolumnsmetforcingdata
               IF (i == 1) THEN
                  headermetout = adjustl(headermet(i))
               ELSE
                  headermetout = trim(headermetout)//' '//adjustl(headermet(i))
               END IF
            END DO
            ! Write out header
            OPEN (78, file=trim(filedscdmet), err=112)
            WRITE (78, '(a)') headermetout
         ELSE
            OPEN (78, file=trim(filedscdmet), position='append') !,err=112)
         END IF
         ! Write out data
         DO i = 1, (readlinesorigmetdatamax*npertstepin)
            WRITE (78, 303) (int(met_tt(i, ii)), ii=1, 4), met_tt(i, 5:ncolumnsmetforcingdata)
         END DO
         !IF(iBlock == ReadBlocksOrigMetData) THEN
         !   WRITE(78,'(i2)') -9
         !   WRITE(78,'(i2)') -9
         !ENDIF
         CLOSE (78) !Close output file
      END IF
 
303   FORMAT((i4, 1x), 3(i3, 1x), 9(f12.6, 1x), (f9.4, 1x), 10(f9.4, 1x)) !Allows 4 dp for rainfall
 
      ! Deallocate arrays -----------------------------------------------------------------
      DEALLOCATE (metfordisagg)
      DEALLOCATE (metfordisaggprev)
      DEALLOCATE (metfordisaggnext)
 
      RETURN
 
112   CALL errorhint(52, trim(filedscdmet), notused, notused, notusedi)
 
   END SUBROUTINE disaggregatemet
   !======================================================================================
 
   !======================================================================================
   SUBROUTINE disaggregateestm(iBlock)
      ! Subroutine to disaggregate met forcing data to model time-step
      ! HCW 10 Feb 2017
      !======================================================================================
 
      USE sues_data
      USE defaultnotused
 
      IMPLICIT NONE
 
      INTEGER :: lunit = 101
      INTEGER :: tdiff !Time difference (in minutes) between first and second rows of original met forcing file
      INTEGER :: i, ii !counter
      INTEGER :: iBlock
      INTEGER, DIMENSION(NperESTM) :: seq1NperESTM
      INTEGER, DIMENSION(nsd) :: seq1nsd
      INTEGER, DIMENSION(ncolsESTMdata) :: ESTMDisaggMethod ! Stores method to use for disaggregating met data
      REAL(KIND(1D0)), DIMENSION(ncolsESTMdata) :: ESTMArrayOrig
      REAL(KIND(1D0)), DIMENSION(ReadLinesOrigESTMData*NperESTM, ncolsESTMdata) :: ESTM_tt
      CHARACTER(LEN=9), DIMENSION(ncolsESTMdata) :: HeaderESTM
      CHARACTER(LEN=10*ncolsESTMdata) :: HeaderESTMOut
      ! REAL(KIND(1d0)),DIMENSION(ReadLinesOrigESTMData):: dectimeOrig
      ! REAL(KIND(1d0)),DIMENSION(ReadLinesOrigESTMData*NperESTM):: dectimeDscd
      INTEGER :: iostat_var
 
      ! INTEGER, DIMENSION(NperESTM):: temp_iy, temp_id, temp_ih, temp_im, temp_ihm
 
      ! Allocate and initialise arrays to receive original forcing data --------------------
      ALLOCATE (estmfordisagg(readlinesorigestmdata, ncolsestmdata))
      ALLOCATE (estmfordisaggprev(ncolsestmdata))
      ALLOCATE (estmfordisaggnext(ncolsestmdata))
      estmfordisagg(:, :) = -999
      estmfordisaggprev(:) = -999
      estmfordisaggnext(:) = -999
      ! Initialise array to receive disaggregated data
      estm_tt = -999
      ! Intialise disaggregation method
      estmdisaggmethod = -999
 
      ! Generate useful sequences
      seq1nperestm = (/(i, i=1, nperestm, 1)/)
      seq1nsd = (/(i, i=1, nsd, 1)/)
 
      ! Get methods to use for disaggregation from RunControl
      ! (N.B.DisaggMethodESTM is set as 1 or 2 in RunControl; ESTMDisaggMethod is array of ncolsESTMdata used here)
      IF (disaggmethodestm == 1) THEN
         estmdisaggmethod(:) = 10 !linear disaggregation of averages
      ELSEIF (disaggmethodestm == 2) THEN
         estmdisaggmethod(:) = 20 !linear disaggregation of instantaneous values
      ELSE
         CALL errorhint(2, 'Problem in SUEWS_ESTMDisagg: DisaggMethodESTM value should be 1 or 2', &
                        notused, notused, disaggmethodestm)
      END IF
 
      ! Read data ---------------------------------------------------------------------
      IF (diagnosedisaggestm == 1) WRITE (*, *) 'Reading file: ', trim(fileorigestm)
      OPEN (lunit, file=trim(fileorigestm), status='old')
      ! CALL skipHeader(lunit,SkipHeaderMet)  !Skip header -> read header instead
      READ (lunit, *) headerestm
      !write(*,*) HeaderMet
      ! Skip over lines that have already been read and downscaled
      IF (skippedlinesorigestm > 0) THEN
         DO i = 1, skippedlinesorigestm - 1 ! minus 1 here because last line of last block needs to be read again
            READ (lunit, *)
         END DO
         ! Read in last line of previous block
         READ (lunit, *, iostat=iostat_var) estmarrayorig
         estmfordisaggprev(1:ncolsestmdata) = estmarrayorig
      END IF
      ! Read in current block
      DO i = 1, readlinesorigestmdatamax
         READ (lunit, *, iostat=iostat_var) estmarrayorig
         estmfordisagg(i, 1:ncolsestmdata) = estmarrayorig
      END DO
      ! Read in first line of next block (except for last block)
      IF (iblock /= readblocksorigmetdata) THEN
         READ (lunit, *, iostat=iostat_var) estmarrayorig
         estmfordisaggnext(1:ncolsestmdata) = estmarrayorig
      END IF
      CLOSE (lunit)
 
      ! Check resolution of original met forcing data -------------------------------------
      ! Find time difference (in minutes) between first and second row
      tdiff = int(estmfordisagg(2, 4) - estmfordisagg(1, 4)) !Try using minutes
      IF (tdiff == 0) tdiff = int(estmfordisagg(2, 3) - estmfordisagg(1, 3))*60 !If no difference in minutes, try using hours
      IF (tdiff < 0) THEN !If time difference is negative (e.g. change of day), instead use second and third row
         tdiff = int(estmfordisagg(3, 4) - estmfordisagg(2, 4))
         IF (tdiff == 0) tdiff = int(estmfordisagg(3, 3) - estmfordisagg(2, 3))*60 !If no difference in minutes, try using hours
      END IF
      ! Check actual resolution matches specified input resolution
      IF (tdiff /= resolutionfilesinestm/60) THEN
         CALL errorhint(2, 'Problem in SUEWS_ESTMDisagg: timestamps in ESTM forcing file inconsistent with ResolutionFilesInESTM', &
                        REAL(ResolutionFilesInESTM, KIND(1D0)), NotUsed, tdiff*60)
      END IF
 
      ! Disaggregate time columns ---------------------------------------------------------
      IF (diagnose == 1) THEN
         WRITE (*, *) 'Disaggregating ESTM forcing data (', trim(fileorigestm), ') to model time-step...'
      END IF
      CALL disaggregatedatetime(estmfordisagg(:, 1:4), tstep, nperestm, readlinesorigmetdatamax, estm_tt(:, 1:4))
      ! ! Convert to dectime
      ! dectimeOrig = ESTMForDisagg(:,2) + ESTMForDisagg(:,3)/24.0 + ESTMForDisagg(:,4)/(60.0*24.0)
      !
      ! DO i=1,ReadLinesOrigESTMDataMax
      !    ! Downscale dectime using dectimeOrig(i) [becomes timestamp of last subinterval]
      !    dectimeDscd(NperESTM*(i-1)+Seq1NperESTM) = dectimeOrig(i) - (tstep/60.0)/(60.0*24.0)*(/(ii, ii=(NperESTM-1),0, -1)/)
      !    ! Convert to required formats
      !    temp_iy   = INT(ESTMForDisagg(i,1))   !Copy year
      !    temp_id   = FLOOR(dectimeDscd(NperESTM*(i-1)+Seq1NperESTM))   !DOY
      !    ! To avoid precision errors, round here
      !    !  - this should not be a problem as a difference of 1 = 1 min, so a difference of 0.001 << 1 min
      !    temp_ihm  = NINT(((dectimeDscd(NperESTM*(i-1)+Seq1NperESTM) - temp_id/1.0)*60.0*24.0)*1000.0)/1000   !Minutes of the day (1440 max)
      !    temp_ih = (temp_ihm-MOD(temp_ihm,60))/60   !Hours
      !    temp_im = MOD(temp_ihm,60)   !Minutes
      !
      !    IF(dectimeOrig(i) == 1.0000 .AND. i > 1) THEN   !If year changes and it is not the beginning of the dataset
      !       WRITE(*,*) 'Year change encountered: ', dectimeOrig(i), dectimeOrig(i-1)
      !       ! Re-downscale dectime using dectimeOrig(i-1)
      !       dectimeDscd(NperESTM*(i-1)+Seq1NperESTM) = dectimeOrig(i-1) + (tstep/60.0)/(60.0*24.0)*Seq1NperESTM
      !       ! Convert to required formats
      !       temp_iy   = INT(ESTMForDisagg(i,1))   !Copy year
      !       temp_id   = FLOOR(dectimeDscd(NperESTM*(i-1)+Seq1NperESTM))   !DOY
      !       temp_ihm  = NINT(((dectimeDscd(NperESTM*(i-1)+Seq1NperESTM) - temp_id/1.0)*60.0*24.0)*1000.0)/1000   !Mins of the day (1440 max)
      !       temp_ih = (temp_ihm-MOD(temp_ihm,60))/60   !Hours
      !       temp_im = MOD(temp_ihm,60)   !Minutes
      !       ! Adjust year and DOY to account for year change
      !       temp_iy(1:(NperESTM-1)) = temp_iy(1:(NperESTM-1)) - 1  !Subtract 1 from year for all except final timestamp
      !       temp_id(NperESTM) = 1  !Set day for final timestamp to 1
      !    ENDIF
      !
      !    !IF(i==1 .or. i== ReadlinesOrigESTMDataMax) THEN
      !    !   write(*,*) temp_iy
      !    !   write(*,*) temp_id
      !    !   !write(*,*) temp_ihm
      !    !   write(*,*) temp_ih
      !    !   write(*,*) temp_im
      !    !ENDIF
      !
      !    ! Copy to ESTM_tt array
      !    ESTM_tt(NperESTM*(i-1)+Seq1NperESTM,1) = temp_iy
      !    ESTM_tt(NperESTM*(i-1)+Seq1NperESTM,2) = temp_id
      !    ESTM_tt(NperESTM*(i-1)+Seq1NperESTM,3) = temp_ih
      !    ESTM_tt(NperESTM*(i-1)+Seq1NperESTM,4) = temp_im
      !
      ! ENDDO
 
      ! Disaggregate other columns --------------------------------------------------------
      ! All other columns are temperatures
      DO ii = 5, ncolsestmdata
         IF (all(estmfordisagg(:, ii) == -999)) THEN
            !IF(DiagnoseDisaggESTM==1) write(*,*) 'No data for col.', ii
            estm_tt(:, ii) = -999
         ELSE
            estm_tt(:, ii) = disagg_lin(estmfordisagg(:, ii), estmfordisaggprev(ii), estmfordisaggnext(ii), estmdisaggmethod(ii), &
                                        nperestm, readlinesorigestmdata, readlinesorigestmdatamax, iblock)
         END IF
      END DO
 
      ! Copy disaggregated data to MetForcingDataArray
      estmforcingdata(:, 1:ncolsestmdata, gridcounter) = estm_tt(:, 1:ncolsestmdata)
 
      ! Write out disaggregated file ------------------------------------------------------
      IF (keeptstepfilesin == 1) THEN
         IF (iblock == 1) THEN
            ! Prepare header
            DO i = 1, ncolsestmdata
               IF (i == 1) THEN
                  headerestmout = adjustl(headerestm(i))
               ELSE
                  headerestmout = trim(headerestmout)//' '//adjustl(headerestm(i))
               END IF
            END DO
            ! Write out header
            OPEN (78, file=trim(filedscdestm), err=113)
            WRITE (78, '(a)') headerestmout
         ELSE
            OPEN (78, file=trim(filedscdestm), position='append') !,err=113)
         END IF
         ! Write out data
         DO i = 1, (readlinesorigestmdatamax*nperestm)
            WRITE (78, 304) (int(estm_tt(i, ii)), ii=1, 4), estm_tt(i, 5:ncolsestmdata)
         END DO
         !IF(iBlock == ReadBlocksOrigMetData) THEN
         !   WRITE(78,'(i2)') -9
         !   WRITE(78,'(i2)') -9
         !ENDIF
         CLOSE (78) !Close output file
      END IF
 
304   FORMAT((i4, 1x), 3(i3, 1x), 9(f9.4, 1x))
 
      ! Deallocate arrays -----------------------------------------------------------------
      DEALLOCATE (estmfordisagg)
      DEALLOCATE (estmfordisaggprev)
      DEALLOCATE (estmfordisaggnext)
 
      RETURN
 
113   CALL errorhint(52, trim(filedscdestm), notused, notused, notusedi)
 
   END SUBROUTINE disaggregateestm
   !======================================================================================
 
   !======================================================================================
   SUBROUTINE disaggregatedatetime(DateTimeForDisagg, tstep, Nper, ReadLinesOrigMetDataMax, DateTimeDscd)
      USE datetime_module, ONLY: daysinyear
      IMPLICIT NONE
      INTEGER, INTENT(in) :: tstep, Nper, ReadLinesOrigMetDataMax
      REAL(KIND(1D0)), DIMENSION(ReadLinesOrigMetData, 4), INTENT(in) :: DateTimeForDisagg
      REAL(KIND(1D0)), DIMENSION(ReadLinesOrigMetData*Nper, 4), INTENT(out) :: DateTimeDscd
 
      REAL(KIND(1D0)), DIMENSION(ReadLinesOrigMetData) :: dectimeOrig
      REAL(KIND(1D0)), DIMENSION(Nper) :: temp_dectime !temorary varaibles for disaggragation
      INTEGER, DIMENSION(Nper) :: temp_iy, temp_id, temp_ih, temp_im, temp_ihm ! temorary varaibles for disaggragation
      INTEGER, DIMENSION(Nper) :: ndays_iy ! number of days in iy
 
      INTEGER :: i, ii
      INTEGER, DIMENSION(Nper) :: seq1Nper
 
      ! Generate useful sequences
      seq1nper = (/(i, i=1, nper, 1)/)
      ! Convert to dectime
      ! dectimeOrig = MetForDisagg(:,2) + MetForDisagg(:,3)/24.0 + MetForDisagg(:,4)/(60.0*24.0)
      ! correct to dectime(year_start)=0 and dectime(year_end)=days of year (i.e., 365 or 366 if leap year) ! TS 09 May 2018
      dectimeorig = (datetimefordisagg(:, 2) - 1) + datetimefordisagg(:, 3)/24.0 + datetimefordisagg(:, 4)/(60.0*24.0)
 
      DO i = 1, readlinesorigmetdatamax
 
         ! Downscale dectime using dectimeOrig(i) [becomes timestamp of last subinterval]
         temp_dectime = dectimeorig(i) - (tstep/60.0)/(60.0*24.0)*(/(ii, ii=(nper - 1), 0, -1)/)
         temp_dectime = nint(temp_dectime*60*60*24)/(60*60*24*1.)
 
         ! Convert to required formats
         ! get year: year-1 if dectime <0, copy `year` from metforcing otherwise
         temp_iy = merge(int(datetimefordisagg(i, 1)) - 1, int(datetimefordisagg(i, 1)), temp_dectime < 0)
 
         ! get day of year:
         ndays_iy = daysinyear(temp_iy) ! get days of year for DOY correction when temp_dectime<0 during year-crossing
         temp_dectime = merge(temp_dectime + ndays_iy, temp_dectime, temp_dectime < 0) ! correct minus temp_dectime to positive values
         temp_id = floor(temp_dectime) + 1 !DOY
 
         temp_ihm = nint((temp_dectime + 1 - temp_id/1.0)*60.0*24.0) !Minutes of the day (1440 max)
         temp_ih = (temp_ihm - mod(temp_ihm, 60))/60 !Hours
         temp_ih = merge(temp_ih, 0, mask=(temp_ih < 24))
         temp_im = mod(temp_ihm, 60) !Minutes
 
         ! Copy to Met_tt array
         datetimedscd(nper*(i - 1) + seq1nper, 1) = temp_iy
         datetimedscd(nper*(i - 1) + seq1nper, 2) = temp_id
         datetimedscd(nper*(i - 1) + seq1nper, 3) = temp_ih
         datetimedscd(nper*(i - 1) + seq1nper, 4) = temp_im
 
      END DO
 
   END SUBROUTINE disaggregatedatetime
   !======================================================================================
 
   ! Define functions here:
   FUNCTION disagg_lin(Slow, SlowPrev, SlowNext, DisaggType, Nper_loc, ReadLinesOrig_loc, ReadLinesOrigMax_loc, iBlock) RESULT(Fast)
 
      USE defaultnotused
      USE sues_data
 
      IMPLICIT NONE
 
      INTEGER :: disaggtype !Type of disaggregation: 10 for averaged variables; 20 for instantaneous variables
      INTEGER :: nper_loc !Number of subintervals per interval (local Nper)
      INTEGER :: readlinesorig_loc, readlinesorigmax_loc !Number of lines to read in original file (local)
      INTEGER :: iblock
      REAL(kind(1d0)), DIMENSION(ReadLinesOrig_loc*Nper_loc) :: fast !Array to receive disaggregated data
      REAL(kind(1d0)), DIMENSION(ReadLinesOrig_loc) :: slow !Array to disaggregate
      REAL(kind(1d0)) :: slowprev, slownext
      INTEGER, DIMENSION(Nper_loc) :: fastrows !Group of rows that are filled with each iteration
      INTEGER, DIMENSION(FLOOR(Nper_loc/2.0)) :: firstrows10 !Rows at the beginning that are not filled during iteration (for averages)
      INTEGER, DIMENSION(Nper_loc - FLOOR(Nper_loc/2.0)) :: lastrows10 !Rows at the end that are not filled during iteration
      INTEGER, DIMENSION(Nper_loc) :: firstrows20 !Rows at the beginning that are not filled during iteration (for instantaneous)
      INTEGER, DIMENSION(Nper_loc) :: seq1nper_loc !1 to Nper_loc
      INTEGER :: xnper_loc !XNper_loc = 2 for even Nper_loc; XNper_loc=1 for odd Nper_loc
      INTEGER :: i, ii !counters
 
      ! Calculate XNper_loc (differentiates between disaggregations with odd and even Nper_loc)
      IF (mod(nper_loc, 2) == 0) xnper_loc = 2
      IF (mod(nper_loc, 2) == 1) xnper_loc = 1
 
      seq1nper_loc = (/(i, i=1, nper_loc, 1)/)
 
      ! Setup counters for iteration
      IF (disaggtype == 10) THEN
         fastrows = floor(nper_loc/2.0) + seq1nper_loc ! Rows to create at model time-step
         firstrows10 = (/(i, i=1, (fastrows(1) - 1), 1)/) !For start of dataset
         lastrows10 = &
            (/(i, i=nper_loc*(readlinesorigmax_loc - 1 - 1) + fastrows(nper_loc) + 1, &
               (readlinesorigmax_loc*nper_loc), 1)/) ! For end of dataset
      ELSEIF (disaggtype == 20) THEN
         fastrows = nper_loc + seq1nper_loc !Rows to create at model time-step
         firstrows20 = (/(i, i=1, (fastrows(1) - 1), 1)/) !For start of dataset
      END IF
 
      ! Initialise fast array to -999
      fast = -999
      ! Linearly disaggregate
      IF (disaggtype == 10) THEN !Averaged variables
         IF (diagnosedisagg == 1) WRITE (*, *) 'Linearly disaggregating averaged variable'
         DO i = 1, (readlinesorigmax_loc - 1)
            fast(nper_loc*(i - 1) + fastrows) = slow(i) - &
                                                (slow(i + 1) - slow(i))/(xnper_loc*nper_loc) + &
                                                (slow(i + 1) - slow(i))/nper_loc*(/(ii, ii=1, nper_loc, 1)/)
         END DO
 
         ! For first few rows, use previous met block
         IF (iblock == 1) THEN
            fast(firstrows10) = fast(fastrows(1)) !Use repeat values at the start of the year
         ELSE
            fast(firstrows10) = slowprev - &
                                (slow(1) - slowprev)/(xnper_loc*nper_loc) + &
                                (slow(1) - slowprev)/nper_loc* &
                                (/(ii, ii=(nper_loc - SIZE(firstrows10) + 1), nper_loc, 1)/)
         END IF
         ! For last few rows, use next met block
         IF (iblock == readblocksorigmetdata) THEN
            fast(lastrows10) = fast(nper_loc*(readlinesorigmax_loc - 1 - 1) + fastrows(nper_loc)) !Use repeat values at the end of the year
         ELSE
            fast(lastrows10) = slow(readlinesorigmax_loc) - &
                               (slownext - slow(readlinesorigmax_loc))/(xnper_loc*nper_loc) + &
                               (slownext - slow(readlinesorigmax_loc))/nper_loc* &
                               (/(ii, ii=1, SIZE(lastrows10), 1)/)
         END IF
      ELSEIF (disaggtype == 20) THEN !Instantaneous variables
         IF (diagnosedisagg == 1) WRITE (*, *) 'Linearly disaggregating instantaneous variable'
         DO i = 1, (readlinesorigmax_loc - 1)
            fast(nper_loc*(i - 1) + fastrows) = (slow(i) + &
                                                 (slow(i + 1) - slow(i))/nper_loc*2*(seq1nper_loc - 1) + &
                                                 slow(i))/2
         END DO
         ! For first few rows, use previous met block
         IF (iblock == 1) THEN
            fast(firstrows20) = fast(fastrows(1)) !Use repeat values at the start of the year
         ELSE
            fast(firstrows20) = (slowprev + &
                                 (slow(1) - slowprev)/nper_loc*2* &
                                 ((/(ii, ii=(nper_loc - SIZE(firstrows20) + 1), nper_loc, 1)/) - 1) + &
                                 slowprev)/2
         END IF
         !! Last few rows are already filled for the instantaneous value disaggregation
         !IF(iBlock==ReadBlocksOrigMetData) THEN
         !   Fast(LastRows20) = Fast(Nper_loc*(ReadLinesOrigMax_loc-1-1)+FastRows(Nper_loc))   !Use repeat values at the end of the year
         !ELSE
         !   Fast(LastRows20) = (Slow(ReadLinesOrigMax_loc) + &
         !                     (SlowNext-Slow(ReadLinesOrigMax_loc))/Nper_loc*2 * &
         !                        ((/(ii, ii=1,SIZE(LastRows20), 1)/)-1) + &
         !                       Slow(ReadLinesOrigMax_loc))/2
         !ENDIF
      END IF
 
      IF (any(fast(1:readlinesorigmax_loc*nper_loc) == -999)) THEN
         WRITE (*, *) 'Problem: -999s (', count(fast(1:readlinesorigmax_loc*nper_loc) == -999), ') in disaggregated data.'
         CALL errorhint(13, 'Problem in SUEWS_MetDisagg: -999 values in disaggregated data.', notused, notused, notusedi)
      END IF
 
   END FUNCTION disagg_lin
   !======================================================================================
 
   !======================================================================================
   FUNCTION disaggp_amongn(Slow, amongN, Nper_loc, ReadLinesOrig_loc, ReadLinesOrigMax_loc) RESULT(Fast)
      ! Subroutine to disaggregate precipitation by evenly distributing among N subintervals
      !  (i.e. equal intensity in N subintervals)
      ! See Ward et al. (in review), meanN, 0.5N or 0.25N approach
      ! HCW 10 Feb 2017
      !======================================================================================
 
      USE defaultnotused
      USE sues_data
 
      IMPLICIT NONE
 
      INTEGER :: amongn !Number of subintervals over which rain will be distributed
      INTEGER :: nper_loc !Number of subintervals per interval (local Nper)
      INTEGER :: readlinesorig_loc, readlinesorigmax_loc !Number of lines to read in original file (local)
      REAL(kind(1d0)), DIMENSION(ReadLinesOrig_loc*Nper_loc) :: fast !Array to receive disaggregated data
      REAL(kind(1d0)), DIMENSION(ReadLinesOrig_loc) :: slow !Array to disaggregate
      INTEGER, DIMENSION(:), ALLOCATABLE :: subintervals !Array of subintervals that contain rain
      INTEGER, DIMENSION(Nper_loc) :: seq1nper_loc !1 to Nper_loc
      INTEGER :: i
 
      ! For each averaging period, get subintervals which will receive rain
      ALLOCATE (subintervals(amongn))
      subintervals(:) = -999
 
      seq1nper_loc = (/(i, i=1, nper_loc, 1)/)
 
      IF (diagnosedisagg == 1) WRITE (*, *) 'Distributing over ', amongn, ' subintervals for variable'
      ! If all subintervals are to contain rain, don't need to generate random numbers
      IF (amongn == nper_loc) THEN
         subintervals(:) = seq1nper_loc
      END IF
      IF (amongn > nper_loc) &
         CALL errorhint(2, 'Problem in SUEWS_MetDisagg: no. of rainy periods cannot exceed number of subintervals', &
                        REAL(nper_loc, kind(1d0)), notused, amongn)
 
      ! Initialise fast array to -999
      fast = -999
      DO i = 1, readlinesorigmax_loc
         fast(nper_loc*(i - 1) + seq1nper_loc) = 0 !Fill all subintervals with zeros initially
         IF (slow(i) > 0) THEN !If there is some rainfall during this interval...
            IF (amongn < nper_loc) THEN
               subintervals(:) = -999
               subintervals = randomsamples(amongn, nper_loc)
            END IF
            fast(nper_loc*(i - 1) + subintervals) = slow(i)/amongn
         END IF
      END DO
 
      IF (any(fast(1:readlinesorigmax_loc*nper_loc) == -999)) THEN
         WRITE (*, *) 'Problem: -999s (', count(fast(1:readlinesorigmax_loc*nper_loc) == -999), ') in disaggregated data'
         CALL errorhint(13, 'Problem in SUEWS_MetDisagg: -999 values in disaggregated data.', notused, notused, notusedi)
      END IF
 
   END FUNCTION disaggp_amongn
   !======================================================================================
 
   !======================================================================================
   FUNCTION disaggp_amongnmult(Slow, multupperI, multamongN, Nper_loc, ReadLinesOrig_loc, ReadLinesOrigMax_loc) RESULT(Fast)
      ! Subroutine to disaggregate precipitation by evenly distributing among N subintervals
      !  (i.e. equal intensity in N subintervals) for different intensity bins
      ! Based on analsysis by Wen Gu
      ! HCW 21 Apr 2017
      !======================================================================================
 
      USE defaultnotused
      USE sues_data
 
      IMPLICIT NONE
 
      REAL(kind(1d0)), DIMENSION(5) :: multupperi !Upper bound of intensity bin
      INTEGER, DIMENSION(5) :: multamongn !Number of subintervals over which rain will be distributed (array)
      INTEGER :: thisamongn !Number of subintervals over which rain will be distributed
      INTEGER :: nper_loc !Number of subintervals per interval (local Nper)
      INTEGER :: readlinesorig_loc, readlinesorigmax_loc !Number of lines to read in original file (local)
      REAL(kind(1d0)), DIMENSION(ReadLinesOrig_loc*Nper_loc) :: fast !Array to receive disaggregated data
      REAL(kind(1d0)), DIMENSION(ReadLinesOrig_loc) :: slow !Array to disaggregate
      INTEGER, DIMENSION(:), ALLOCATABLE :: subintervals !Array of subintervals that contain rain
      INTEGER, DIMENSION(Nper_loc) :: seq1nper_loc !1 to Nper_loc
      INTEGER :: i
 
      seq1nper_loc = (/(i, i=1, nper_loc, 1)/)
 
      IF (diagnosedisagg == 1) WRITE (*, *) 'Distributing over variable subintervals depending on intensity for variable'
 
      ! Initialise fast array to -999
      fast = -999
      DO i = 1, readlinesorigmax_loc
         fast(nper_loc*(i - 1) + seq1nper_loc) = 0 !Fill all subintervals with zeros initially
         IF (slow(i) > 0) THEN !If there is some rainfall during this interval...
            !Use intensity in this interval to decide number of subintervals to fill with rain
            IF (slow(i) <= multupperi(1)) THEN
               thisamongn = multamongn(1)
            ELSEIF (slow(i) > multupperi(1) .AND. slow(i) <= multupperi(2)) THEN
               thisamongn = multamongn(2)
            ELSEIF (slow(i) > multupperi(2) .AND. slow(i) <= multupperi(3)) THEN
               thisamongn = multamongn(3)
            ELSEIF (slow(i) > multupperi(3) .AND. slow(i) <= multupperi(4)) THEN
               thisamongn = multamongn(4)
            ELSEIF (slow(i) > multupperi(4) .AND. slow(i) <= multupperi(5)) THEN
               thisamongn = multamongn(5)
            ELSEIF (slow(i) > multupperi(5)) THEN
               thisamongn = multamongn(5)
               CALL errorhint(4, 'Precip in met forcing file exceeds maxiumum MultRainAmongNUpperI', &
                              slow(i), multrainamongnupperi(5), notusedi)
            END IF
 
            ! For each averaging period, get subintervals which will receive rain
            ALLOCATE (subintervals(thisamongn))
            subintervals(:) = -999
 
            IF (thisamongn > nper_loc) CALL errorhint(2, 'Problem in SUEWS_MetDisagg: no. of rainy periods cannot exceed '// &
                                                      'number of subintervals', real(nper_loc, kind(1d0)), notused, thisamongn)
 
            IF (thisamongn == nper_loc) THEN ! If all subintervals are to contain rain, don't need to generate random numbers
               subintervals(:) = seq1nper_loc
            ELSEIF (thisamongn < nper_loc) THEN
               subintervals = randomsamples(thisamongn, nper_loc)
            END IF
            fast(nper_loc*(i - 1) + subintervals) = slow(i)/thisamongn
            !write(*,*) Slow(i), thisamongN
            DEALLOCATE (subintervals)
         END IF
      END DO
 
      IF (any(fast(1:readlinesorigmax_loc*nper_loc) == -999)) THEN
         WRITE (*, *) 'Problem: -999s (', count(fast(1:readlinesorigmax_loc*nper_loc) == -999), ') in disaggregated data'
         CALL errorhint(13, 'Problem in SUEWS_MetDisagg: -999 values in disaggregated data.', notused, notused, notusedi)
      END IF
 
   END FUNCTION disaggp_amongnmult
   !======================================================================================
 
   !======================================================================================
   FUNCTION randomsamples(N, OutOf) RESULT(Samples)
      ! Generates N/OutOf random samples without repeats
      !   e.g. for N = 3 and OutOf = 12, a possibility for Samples = 7,3,11
      ! HCW 10 Feb 2017
      !======================================================================================
 
      IMPLICIT NONE
 
      INTEGER :: i !counter
      INTEGER :: n !number of samples to return
      INTEGER :: outof !number to sample from
      INTEGER :: x !next sample to be added
      REAL(kind(1d0)) :: r !random number
      INTEGER, DIMENSION(:), ALLOCATABLE :: samples !Array to receive random samples
 
      ! Allocate and initialise Samples
      ALLOCATE (samples(n))
      samples(:) = -999
 
      ! Generate random sample (no repeats)
      i = 0 !Set counter to zero initially
      DO WHILE (any(samples == -999))
         CALL random_number(r)
         x = int(r*outof) + 1
         !write(*,*) X
         !write(*,*) COUNT(Samples == X)
         IF (count(samples == x) == 0) THEN
            ! Only keep if this subinterval has not already been selected
            i = i + 1
            samples(i) = x
         END IF
         !write(*,*) Samples
      END DO
 
   END FUNCTION randomsamples
   !======================================================================================
 
END MODULE metdisagg
!========================================================================================