MODULE CABLE_CRU

  USE netcdf                         ! Access to netcdf routines
       
  USE casa_ncdf_module, ONLY: HANDLE_ERR, GET_UNIT                       ! Finds an unused unit number for file opensYMDHMS2DOYSOD, DOYSOD2YMDHMS

  USE cable_IO_vars_module, ONLY: &  ! Selected cable_iovars.F90 variables:
       logn,            &             ! Log file unit number
       land_x, land_y,  &             ! Col (x) & row (y) indices of each land point in land mask (dimension mland)
       exists                         ! Only for exists%Snowf, which we will set to .FALSE. because there is no snow
  ! in CRU-NCEP. Setting this ensures snow will be determined in CABLE from temperature.

  IMPLICIT NONE

  ! Define a type for CRU-NCEP information, and the subtype METVALS

  TYPE CRU_MET_TYPE
     REAL, DIMENSION(:), ALLOCATABLE :: METVALS  ! Define a spatial vector of meteorology for one timestep
  END TYPE CRU_MET_TYPE

  TYPE CRU_TYPE
     INTEGER  :: mland                    ! Number of land cells
     INTEGER  :: NMET                     ! Number of met variable types (rain, lwdn etc) NOT INCLUDING prevTmax and nextTmin
     INTEGER  :: xdimsize, ydimsize       ! Landmask grid size dimensions (x=cols, y=rows)
     INTEGER  :: tdimsize                 ! Time dimension of metfiles (met data timesteps per annual file)
     INTEGER  :: CYEAR                    ! Current run year, same as CurYear, not necessarily the same as MetYear
     INTEGER  :: MetStart                 ! First year of met
     INTEGER  :: MetEnd                   ! Last year of met
     INTEGER  :: CTSTEP                   ! Current met data timestep (1 to tdimsize, i.e. 365 for CRU-NCEP annual daily files)
     INTEGER  :: DTsecs                   ! Model timestep in seconds, converted from namelist value in hours
     INTEGER  :: ktau                     ! Current model timestep, reset at the start of a new year of met
     INTEGER, DIMENSION(9) :: F_ID, V_ID  ! NetCDF object id's for files and variables (NetCDF bookkeeping stuff)
     REAL, DIMENSION(:) ,ALLOCATABLE :: AVG_LWDN ! Avg of one day's diurnal cycle of lwdn calculated by Swinbank. AVG_LWDN
     ! is used to rescale the diurnal cycle to match the day's CRUNCEP lwdn. (dim=mland)
     REAL, DIMENSION(:) ,ALLOCATABLE :: CO2VALS  ! Global annual CO2 values (dim is the number of years of data, or 1 if time-invariant)
     LOGICAL  :: DirectRead     ! Flag to do with reading small numbers of points efficiently. Set true for small numbers of points
     LOGICAL  :: LeapYears      ! Flag for whether leaps years occur, required by CABLE. Always false for CRUNCEP (no Feb 29th)
     LOGICAL,DIMENSION(:,:),ALLOCATABLE :: LandMask ! Logical landmask, true for land, false for non-land
     !
     CHARACTER(len=30)  :: Run            ! Where run type is      : "S0_TRENDY", "S1_TRENDY", "S2_TRENDY"
     CHARACTER(len=15)  :: CO2            ! CO2 takes value        : "static1860", "1860_1900", "1901_2015"
     CHARACTER(len=15)  :: Ndep            ! Ndep takes value        : "static1860", "1860_1900", "1901_2015"
     CHARACTER(len=15)  :: Forcing        ! Met Forcing takes value: "spinup",        "spinup", "1901_2015"
     !
     CHARACTER(len=200) :: BasePath       ! Full path for the location of data used for CRU runs "/x/y"
     CHARACTER(len=200) :: MetPath        ! Full path for the location of the met files "/x/y"
     CHARACTER(len=200) :: LandMaskFile   ! Land mask filename, without path
     CHARACTER(len=30) ,DIMENSION(9) :: VAR_NAME  ! Netcdf variable 'Name' for each type of met (dim=# of met vars). Note: Name, not 'Title'
     CHARACTER(len=200),DIMENSION(9) :: MetFile   ! Met file names incl metpath, constructed in CRU_GET_FILENAME (dim=# of met vars)
     TYPE(CRU_MET_TYPE), DIMENSION(11) :: MET  ! Met data vectors (METVALS) for one timestep, dim=# of met vars + 2 for prev Tmax and next Tmin
     REAL,  DIMENSION(:), ALLOCATABLE :: NdepVALS
     INTEGER :: NdepF_ID, NdepV_ID
     INTEGER  :: Ndep_CTSTEP   ! counter for Ndep in input file
  END TYPE CRU_TYPE

  TYPE (CRU_TYPE):: CRU  ! Define the variable CRU, of type CRU_TYPE

  ! Define local parameter names representing the position of each met var within variable MET.
  ! prevTmax and nextTmin are special cases of Tmax and Tmin that do not count as extra met variables per se.
  INTEGER, PRIVATE, PARAMETER :: &
       rain     =  1, &
       lwdn     =  2, &
       swdn     =  3, &
       pres     =  4, &
       qair     =  5, &
       tmax     =  6, &
       tmin     =  7, &
       uwind    =  8, &
       vwind    =  9, &
                                !
       prevTmax = 10, &
       nextTmin = 11

  ! Error status of various operations (mostly netcdf-related). Typically 0 means ok, > 0 means unexpected condition.
  INTEGER, PRIVATE :: ErrStatus

  REAL, PRIVATE, PARAMETER :: SecDay = 86400. ! Number of seconds in a day

  ! Filename prefix expected in the names of met files. Used by CRU_GET_FILENAME to construct met file names.
  CHARACTER(len=6), DIMENSION(9), PARAMETER, PRIVATE :: &
       PREF = (/ "rain  ", "lwdown", "swdown", "press ", "qair  ", "tmax  ", "tmin  ", "uwind ", "vwind " /)

CONTAINS

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

  SUBROUTINE CRU_INIT( CRU )

    ! Initialise the contents of the CRU defined type collection, from the CRU namelist file
    ! and by obtaining dimensions from the landmask

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

    USE cable_IO_vars_module, ONLY: &
         latitude, longitude, & ! (R) Lat and long of landcells only (?)
         nmetpatches,         & ! (I) Size of patch dimension in met file, if it exists
         mask,                & ! (I) Land/sea mask (1,0)
         metGrid,             & ! (C4) Either 'land' or 'mask' for whether the data are packed or not (?)
         sdoy, smoy, syear,   & ! (I) Start time day of year, month, year
         shod,                & ! (R) Start time hour of day
         xdimsize, ydimsize,  & ! (I) Size of grid dimensions
         lat_all, lon_all       ! (R) Grids with the lat or lon of each cell (i.e. repetition along rows/cols), for CABLE.

    USE cable_def_types_mod,  ONLY: mland  ! (I) Number of land cells

    IMPLICIT NONE

    TYPE (CRU_TYPE):: CRU

    INTEGER              :: ErrStatus  ! Error status returned by nc routines (zero=ok, non-zero=error)
    INTEGER              :: nmlunit    ! Unit number for reading namelist file
    INTEGER              :: FID        ! NetCDF id for the landmask file
    INTEGER              :: latID, lonID, timID  ! NetCDF ids for dimensions in the landmask file
    INTEGER              :: landID     ! NetCDF id for the landmask variable in the landmask file
    INTEGER              :: tdimsize   ! Time dimension in the met file, not used here apparently (??)
    INTEGER              :: landcnt    ! Manually incremented counter for the number of land cells
    INTEGER              :: xcol, yrow ! Column and row position in the data file grids
    INTEGER              :: imetvar    ! loop counter through met variables

    ! Temporary local names for CRU% variables as they are read from the namelist file.
    ! Note that CRU%CO2 and CRU%Forcing are assigned based on the value of Run, not read as options from the namelist file.
    LOGICAL              :: DirectRead = .FALSE.
    CHARACTER(len=30)    :: Run
    CHARACTER(len=200)   :: BasePath
    CHARACTER(len=200)   :: MetPath
    CHARACTER(len=200)   :: LandMaskFile
    REAL                 :: DThrs   ! CABLE timestep (hrs), converted immediately to integer seconds for CRU%DTsecs
    REAL,DIMENSION(:)     ,ALLOCATABLE :: CRU_lats, CRU_lons  ! Lat/long values for each grid rows/cols from landmask.
    INTEGER,DIMENSION(:,:),ALLOCATABLE :: landmask

    ! Flag for errors
    LOGICAL              :: ERR = .FALSE.

    NAMELIST /CRUNML/ BasePath, MetPath, LandMaskFile, Run, DThrs, DirectRead

    ! Read CRU namelist settings
    CALL GET_UNIT(nmlunit)  ! CABLE routine finds spare unit number
    OPEN (nmlunit,FILE="cru.nml",STATUS='OLD',ACTION='READ')
    READ (nmlunit,NML=CRUNML)
    CLOSE(nmlunit)

    ! Assign namelist settings to corresponding CRU defined-type elements
    CRU%BasePath     = BasePath
    CRU%MetPath      = MetPath
    CRU%LandMaskFile = LandMaskFile
    CRU%Run          = Run
    CRU%DTsecs       = INT(DThrs * 3600.)  ! in seconds
    CRU%DirectRead   = DirectRead

    ! Assign Forcing and CO2 labels based only on the value of CRU%Run
    SELECT CASE (TRIM(CRU%Run))
    CASE( "S0_TRENDY" )
       CRU%Forcing = "spinup"
       CRU%CO2     = "static1860"
       CRU%Ndep     = "static1860"
       WRITE(*   ,*)"Run = 'spinup': Therefore Forcing = 'spinup', CO2 = 'static1860'"
       WRITE(logn,*)"Run = 'spinup': Therefore Forcing = 'spinup', CO2 = 'static1860'"
    CASE( "S0_TRENDY_CO2" )
       CRU%Forcing = "spinup"
       CRU%CO2     = "1901_2015"
       CRU%Ndep     = "static1860"
       WRITE(*   ,*)"Run = 'S0_CO2': Therefore Forcing = 'S0_CO2', CO2 = '1860_2015'"
       WRITE(logn,*)"Run = 'S0_CO2': Therefore Forcing = 'S0_CO2', CO2 = '1860_2015'"
    CASE( "S0_TRENDY_Ndep" )
       CRU%Forcing = "spinup"
       CRU%CO2     = "static1860"
       CRU%Ndep     = "1901_2015"
       WRITE(*   ,*)"Run = 'S0_Ndep': Therefore Forcing = 'S0_Ndep', Ndep = '1860_2015'"
       WRITE(logn,*)"Run = 'S0_Ndep': Therefore Forcing = 'S0_Ndep', Ndep = '1860_2015'"
    CASE( "S0_TRENDY_Precip" )
       CRU%Forcing = "spinup"
       CRU%CO2     = "static1860"
       CRU%Ndep     = "static1860"
       WRITE(*   ,*)"Run = 'spinup': Therefore Forcing = 'spinup', CO2 = 'static1860'"
       WRITE(logn,*)"Run = 'spinup': Therefore Forcing = 'spinup', CO2 = 'static1860'"
    CASE( "S0_TRENDY_Temp" )
       CRU%Forcing = "spinup"
       CRU%CO2     = "static1860"
       CRU%Ndep     = "static1860"
       WRITE(*   ,*)"Run = 'spinup': Therefore Forcing = 'spinup', CO2 = 'static1860'"
       WRITE(logn,*)"Run = 'spinup': Therefore Forcing = 'spinup', CO2 = 'static1860'"
    CASE( "S0_TRENDY_Temp_Precip" )
       CRU%Forcing = "spinup"
       CRU%CO2     = "static1860"
       CRU%Ndep     = "static1860"
       WRITE(*   ,*)"Run = 'spinup': Therefore Forcing = 'spinup', CO2 = 'static1860'"
       WRITE(logn,*)"Run = 'spinup': Therefore Forcing = 'spinup', CO2 = 'static1860'"
    CASE( "S0_TRENDY_CO2_Temp" )
       CRU%Forcing = "spinup"
       CRU%CO2     = "1901_2015"
       CRU%Ndep     = "static1860"
       WRITE(*   ,*)"Run = 'spinup': Therefore Forcing = 'spinup', CO2 = '1860_2015'"
       WRITE(logn,*)"Run = 'spinup': Therefore Forcing = 'spinup', CO2 = '1860_2015'"
    CASE( "S0_TRENDY_CO2_Precip" )
       CRU%Forcing = "spinup"
       CRU%CO2     = "1901_2015"
       CRU%Ndep     = "static1860"
       WRITE(*   ,*)"Run = 'spinup': Therefore Forcing = 'spinup', CO2 = '1860_2015'"
       WRITE(logn,*)"Run = 'spinup': Therefore Forcing = 'spinup', CO2 = '1860_2015'"
    CASE( "S0_TRENDY_CO2_Temp_Precip" )
       CRU%Forcing = "spinup"
       CRU%CO2     = "1901_2015"
       CRU%Ndep     = "static1860"
       WRITE(*   ,*)"Run = 'spinup': Therefore Forcing = 'spinup', CO2 = '1860_2015'"
       WRITE(logn,*)"Run = 'spinup': Therefore Forcing = 'spinup', CO2 = '1860_2015'"
    CASE( "S1_TRENDY" )
       CRU%Forcing = "spinup"
       CRU%CO2     = "1860_1900"
       CRU%Ndep     = "1860_1900"
       WRITE(*   ,*)"Run = 'S1_TRENDY': Therefore Forcing = 'spinup', CO2 = '1860_1900'"
       WRITE(logn,*)"Run = 'S1_TRENDY': Therefore Forcing = 'spinup', CO2 = '1860_1900'"
    CASE( "S2_TRENDY" )
       CRU%Forcing = "1901_2015"
       CRU%CO2     = "1901_2015"
       CRU%Ndep     = "1901_2015"
       WRITE(*   ,*)"Run = 'S2_TRENDY': Therefore Forcing = 'spinup', CO2 = '1901_2015'"
       WRITE(logn,*)"Run = 'S2_TRENDY': Therefore Forcing = 'spinup', CO2 = '1901_2015'"
    CASE( "S2_TRENDY_precip" )
       CRU%Forcing = "1901_2015"
       CRU%CO2     = "1901_2015"
       CRU%Ndep     = "1901_2015"
       WRITE(*   ,*)"Run = 'S2_TRENDY': Therefore Forcing = 'spinup', CO2 = '1901_2015'"
       WRITE(logn,*)"Run = 'S2_TRENDY': Therefore Forcing = 'spinup', CO2 = '1901_2015'"
    CASE( "S2_TRENDY_precip0" )
       CRU%Forcing = "1901_2015"
       CRU%CO2     = "1901_2015"
       CRU%Ndep     = "1901_2015"
       WRITE(*   ,*)"Run = 'S2_TRENDY': Therefore Forcing = 'spinup', CO2 = '1901_2015'"
       WRITE(logn,*)"Run = 'S2_TRENDY': Therefore Forcing = 'spinup', CO2 = '1901_2015'"
    CASE  default
       WRITE(*   ,*)"Wrong CRU%Run: ",CRU%Run
       WRITE(*   ,*)"Use: S0_TRENDY, S1_TRENDY, or S2_TRENDY!"
       WRITE(logn,*)"Wrong CRU%Run: ",CRU%Run
       WRITE(logn,*)"Use: S0_TRENDY, S1_TRENDY, or S2_TRENDY!"
       ERR = .TRUE.
    END SELECT

    ! Print settings
    WRITE(*   ,*)"========================================= CRU ============"
    WRITE(*   ,*)"CRU settings chosen:"
    WRITE(*   ,*)" BasePath: ",TRIM(CRU%BasePath)
    WRITE(*   ,*)" LandMask: ",TRIM(CRU%LandMaskFile)
    WRITE(*   ,*)" Run                : ",TRIM(CRU%Run)
    WRITE(*   ,*)" Forcing (assigned) : ",TRIM(CRU%Forcing)
    WRITE(*   ,*)" CO2     (assigned) : ",TRIM(CRU%CO2)
    WRITE(*   ,*)" Ndep     (assigned) : ",TRIM(CRU%Ndep)
    WRITE(*   ,*)" DT(secs): ",CRU%DTsecs
    WRITE(logn,*)"========================================= CRU ============"
    WRITE(logn,*)"CRU settings chosen:"
    WRITE(logn,*)" BasePath: ",TRIM(CRU%BasePath)
    WRITE(logn,*)" LandMask: ",TRIM(CRU%LandMaskFile)
    WRITE(logn,*)" Run                : ",TRIM(CRU%Run)
    WRITE(logn,*)" Forcing (assigned) : ",TRIM(CRU%Forcing)
    WRITE(logn,*)" CO2     (assigned) : ",TRIM(CRU%CO2)
    WRITE(logn,*)" Ndep     (assigned) : ",TRIM(CRU%Ndep)
    WRITE(logn,*)" DT(secs): ",CRU%DTsecs

    ! Error trap for bad namelist.
    IF ( ERR ) THEN
       WRITE(logn,*)"Invalid settings in CRU_INIT"
       STOP "Invalid settings in CRU_INIT"
    ENDIF

    ! If this is a S0_TRENDY run look for met data in the spinup directory instead.
    !IF (TRIM(CRU%Run) .EQ. "S0_TRENDY") THEN
    !   CRU%MetPath = TRIM(CRU%MetPath)//"/spinup_data"
    !ENDIF

    ! Set variable names to their NetCDF 'Names' (i.e. not their 'Titles')
    CRU%NMET = 9
    CRU%VAR_NAME(rain)  = "Total_Precipitation"
    CRU%VAR_NAME(lwdn)  = "Incoming_Long_Wave_Radiation"
    CRU%VAR_NAME(swdn)  = "Incoming_Short_Wave_Radiation"
    CRU%VAR_NAME(pres)  = "Pression"
    CRU%VAR_NAME(qair)  = "Air_Specific_Humidity"
    ! CRU%VAR_NAME(tmax)  = "maximum_6h_air_temperature"
    ! CRU%VAR_NAME(tmin)  = "minimum_6h_air_temperature"
    CRU%VAR_NAME(tmax)  = "maximum_air_temperature"
    CRU%VAR_NAME(tmin)  = "minimum_air_temperature"
    CRU%VAR_NAME(uwind) = "U_wind_component"
    CRU%VAR_NAME(vwind) = "V_wind_component"

    WRITE(*   ,*)"========================================= CRU ============"
    WRITE(logn,*)"========================================= CRU ============"

    ! Now read landmask file
    ! Landmask file into init! Get LAt, LON etc. from there
    ! LMFILE = TRIM(CRU%LandMaskFile)
    WRITE(*   ,*) 'Opening CRU landmask file: ',TRIM(LandMaskFile)
    WRITE(logn,*) 'Opening CRU landmask file: ',TRIM(LandMaskFile)

    ! Open the land mask file
    ErrStatus = NF90_OPEN(TRIM(LandMaskFile), NF90_NOWRITE, FID)
    CALL HANDLE_ERR(ErrStatus, "Opening CRU Land-mask file"//TRIM(LandMaskFile))

    ! Latitude: Get the dimension ID, find the size of the dimension, assign it to CRU.
    ErrStatus = NF90_INQ_DIMID(FID,'latitude',latID)
    ErrStatus = NF90_INQUIRE_DIMENSION(FID,latID,len=ydimsize)
    CALL HANDLE_ERR(ErrStatus, "Inquiring 'lat'"//TRIM(LandMaskFile))
    CRU%ydimsize = ydimsize

    ! Collect the latitudes into CRU_lats
    ALLOCATE( CRU_lats ( ydimsize ) )
    ErrStatus = NF90_INQ_VARID(FID,'latitude',latID)
    CALL HANDLE_ERR(ErrStatus, "Inquiring 'latitudes'"//TRIM(LandMaskFile))
    ErrStatus = NF90_GET_VAR(FID,latID,CRU_lats)
    CALL HANDLE_ERR(ErrStatus, "Reading 'latitudes'"//TRIM(LandMaskFile))

    ! Longitude: Get the dimension ID, find the size of the dimension, assign it to CRU.
    ErrStatus = NF90_INQ_DIMID(FID,'longitude',lonID)
    ErrStatus = NF90_INQUIRE_DIMENSION(FID,lonID,len=xdimsize)
    CALL HANDLE_ERR(ErrStatus, "Inquiring 'lon'"//TRIM(LandMaskFile))
    CRU%xdimsize = xdimsize

    ! Collect the longitudes into CRU_lons
    ALLOCATE( CRU_lons ( xdimsize ) )
    ErrStatus = NF90_INQ_VARID(FID,'longitude',lonID)
    CALL HANDLE_ERR(ErrStatus, "Inquiring 'longitudes'"//TRIM(LandMaskFile))
    ErrStatus = NF90_GET_VAR(FID,lonID,CRU_lons)
    CALL HANDLE_ERR(ErrStatus, "Reading 'longitudes'"//TRIM(LandMaskFile))

    ! Allocate the landmask arrays for...
    ALLOCATE( CRU%landmask ( xdimsize, ydimsize) )  ! Passing out to other CRU routines (logical)
    ALLOCATE( landmask ( xdimsize, ydimsize) )      ! Local use in this routine (integer)
    ALLOCATE ( mask( xdimsize, ydimsize) )          ! Use by CABLE

    ! Check that the land mask variable is called "land" in the land mask file,
    ! and read it into local variable landmask
    ErrStatus = NF90_INQ_VARID(FID,'land',landID)
    CALL HANDLE_ERR(ErrStatus, "Inquiring 'land' "//TRIM(LandMaskFile))
    ErrStatus = NF90_GET_VAR(FID,landID,landmask)
    CALL HANDLE_ERR(ErrStatus, "Reading 'land' "//TRIM(LandMaskFile))

    ! Convert the integer landmask into the logical CRU%landmask
    WHERE ( landmask .GT. 0 )
       CRU%landmask = .TRUE.
       mask           = 1
    ELSEWHERE
       CRU%landmask = .FALSE.
       mask           = 0
    END WHERE

    ! Count the number of land cells -> mland
    CRU%mland = COUNT(CRU%landmask)

    ! Allocate CABLE land-only vectors for lat/long and row/col values/indices.
    ALLOCATE( latitude(CRU%mland), longitude(CRU%mland) )
    ALLOCATE( land_y  (CRU%mland), land_x   (CRU%mland) )

    ! Allocate vectors for each of the different met quantities, including extra
    ! prev/next temperatures for the Cesarracio temperature calculations in the
    ! weather generator.
    DO imetvar = 1, CRU%NMET
       ALLOCATE( CRU%MET(imetvar)%METVALS(CRU%mland) )
    END DO
    ALLOCATE( CRU%MET(prevTmax)%METVALS(CRU%mland) )
    ALLOCATE( CRU%MET(nextTmin)%METVALS(CRU%mland) )
    ! allocate array for Nitrogen deposition input data
    ALLOCATE( CRU%NdepVALS(CRU%mland) )

    ! Copy the col/row and lat/long positions of each land cell into the corresponding
    ! land only CABLE vectors. Q: We know mland at this point. Why not use landcnt to confirm
    ! the correct value of mland?
    landcnt = 1
    DO yrow = 1, ydimsize
       DO xcol = 1, xdimsize
          IF ( .NOT. CRU%landmask(xcol,yrow) ) CYCLE   ! Go to next iteration if not a land cell

          !          WRITE(6,FMT='(A15,I5,2(1X,F8.2),2(1x,I3))')"i, lo,la, xcol,yrow",landcnt,CRU_lons(xcol),CRU_lats(yrow),xcol, yrow

          ! C
          land_x(landcnt)    = xcol
          land_y(landcnt)    = yrow
          longitude(landcnt) = CRU_lons(xcol)
          latitude(landcnt)  = CRU_lats(yrow)
          landcnt = landcnt + 1
       END DO
    END DO

    ! Set global CABLE variables
    metGrid     = "mask"
    ALLOCATE( mask(xdimsize, ydimsize) )
    mask        = landmask
    mland       = CRU%mland
    nmetpatches = 1
    ALLOCATE( lat_all(xdimsize, ydimsize), lon_all(xdimsize, ydimsize) )
    DO xcol = 1, xdimsize
       lat_all(xcol,:) = CRU_lats
    END DO
    DO yrow = 1, ydimsize
       lon_all(:,yrow) = CRU_lons
    END DO

    ! CABLE TIME-UNITS needed by load-parameters (only on CABLE_init)
    shod        = 0.
    sdoy        = 1
    smoy        = 1
    syear       = CRU%CYEAR

    ! Used to rescale the diurnal cycle from Swinbank calculation to match CRU-NCEP provided value.
    ALLOCATE( CRU%AVG_LWDN(mland) )

    DEALLOCATE ( landmask, CRU_lats, CRU_lons )

    ErrStatus = NF90_CLOSE(FID)
    CALL HANDLE_ERR(ErrStatus, "Closing mask-file"//TRIM(LandMaskFile))

  END SUBROUTINE CRU_INIT

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

  SUBROUTINE CRU_GET_FILENAME ( CRU, cyear, par, FN )

    ! Build the filename FN: One annual file of daily met for one met quantity.

    IMPLICIT NONE

    TYPE(CRU_TYPE), INTENT(INOUT) :: CRU  ! Information about CRU
    INTEGER,              INTENT(IN)    :: cyear  ! Current year as an integer
    INTEGER,              INTENT(IN)    :: par    ! Index (1-9) of which met quantity will be sought
    CHARACTER(LEN=200),   INTENT(OUT)   :: FN     ! Met filename (outgoing)
    INTEGER   :: i, idx
    CHARACTER(4)   :: cy  ! Character representation of cyear
    CHARACTER(200) :: mp  ! Local repr of met path

    ! Create a character version of the year for building that part of the filename.
    WRITE(cy,FMT='(I4)')cyear

    ! Initialise the filename with the met path
    mp = CRU%MetPath
    FN = TRIM(mp)

    ! Build the rest of the filename according to the value of par, which references 11 possible
    ! types of met through the parameter names rain, lwdn, etc.
!    SELECT CASE ( par )
!    CASE(rain) ; FN = TRIM(FN)//"/rain/cruncep2015_1_rain_"//cy//".daymean.nc"
!    CASE(lwdn) ; FN = TRIM(FN)//"/lwdown/cruncep2015_1_lwdown_"//cy//".daymean.nc"
!    CASE(swdn) ; FN = TRIM(FN)//"/swdown/cruncep2015_1_swdown_"//cy//".daymean.nc"
!    CASE(pres) ; FN = TRIM(FN)//"/press/cruncep2015_1_press_"//cy//".daymean.nc"
!    CASE(qair) ; FN = TRIM(FN)//"/qair/cruncep2015_1_qair_"//cy//".daymean.nc"
!    CASE(tmax,PrevTmax) ; FN = TRIM(FN)//"/tmax/cruncep2015_1_tair_"//cy//".daymax.nc"
!    CASE(tmin,NextTmin) ; FN = TRIM(FN)//"/tmin/cruncep2015_1_tair_"//cy//".daymin.nc"
!    CASE(uwind) ; FN = TRIM(FN)//"/uwind/cruncep2015_1_uwind_"//cy//".daymean.nc"
!    CASE(vwind) ; FN = TRIM(FN)//"/vwind/cruncep2015_1_vwind_"//cy//".daymean.nc"
!    END SELECT


    SELECT CASE ( par )
    CASE(rain) ; FN = TRIM(FN)//"/rain/cruncepV8_rain_"//cy//".daytot.nc"
    CASE(lwdn) ; FN = TRIM(FN)//"/lwdown/cruncepV8_lwdown_"//cy//".daymean.nc"
    CASE(swdn) ; FN = TRIM(FN)//"/swdown/cruncepV8_swdown_"//cy//".daymean.nc"
    CASE(pres) ; FN = TRIM(FN)//"/press/cruncepV8_press_"//cy//".daymean.nc"
    CASE(qair) ; FN = TRIM(FN)//"/qair/cruncepV8_qair_"//cy//".daymean.nc"
    CASE(tmax,PrevTmax) ; FN = TRIM(FN)//"/tmax/cruncepV8_tmax_"//cy//".daymax.nc"
    CASE(tmin,NextTmin) ; FN = TRIM(FN)//"/tmin/cruncepV8_tmin_"//cy//".daymin.nc"
    CASE(uwind) ; FN = TRIM(FN)//"/uwind/cruncepV8_uwind_"//cy//".daymean.nc"
    CASE(vwind) ; FN = TRIM(FN)//"/vwind/cruncepV8_vwind_"//cy//".daymean.nc"
    END SELECT

  END SUBROUTINE CRU_GET_FILENAME

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

  SUBROUTINE GET_CRU_CO2( CRU, CO2air )

    ! Get CO2 values for use with a CRU-NCEP run. Assign a static 1860 value if specified otherwise
    ! on the first call read all the annual values from a file into the CRU%CO2VALS array. On the first
    ! and subsequent

    IMPLICIT NONE

    TYPE(CRU_TYPE) :: CRU           ! All the info needed for CRU met runs
    REAL, INTENT(OUT)    :: CO2air  ! A single annual value of CO2air in ppm for the current year.

    INTEGER              :: i, iunit, iyear, IOS = 0
    CHARACTER            :: CO2FILE*200
    LOGICAL,        SAVE :: CALL1 = .TRUE.  ! A *local* variable recording the first call of this routine

    ! For S0_TRENDY, use only static 1860 CO2 value and return immediately
    IF ( TRIM(CRU%CO2) .EQ. "static1860") THEN
       CO2air = 286.42   ! CO2 in ppm for 1860
       RETURN

       ! If not S0_TRENDY, varying CO2 values will be used...
    ELSE

       ! On the first call, allocate the CRU%CO2VALS array to store the entire history of annual CO2
       ! values, open the (ascii) CO2 file and read the values into the array.
       IF (CALL1) THEN
          ALLOCATE( CRU%CO2VALS( 1750:2016 ) )
          CO2FILE = TRIM(CRU%BasePath)//"/co2/1750_2015_globalCO2_time_series.csv"
          CALL GET_UNIT(iunit)
          OPEN (iunit, FILE=TRIM(CO2FILE), STATUS="OLD", ACTION="READ")
          DO WHILE( IOS .EQ. 0 )
             READ(iunit, FMT=*, IOSTAT=IOS) iyear, CRU%CO2VALS(iyear)
          END DO
          CLOSE(iunit)

          CALL1 = .FALSE.

       END IF

       ! In all varying CO2 cases, return the element of the array for the current year
       ! as a single CO2 value.
       !
       CO2air = CRU%CO2VALS( CRU%CYEAR )

    END IF

  END SUBROUTINE GET_CRU_CO2

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

  SUBROUTINE GET_CRU_Ndep( CRU )

    ! Get Ndep values for use with a CRU-NCEP run. Assign a static 1860 value if specified otherwise
    ! on the first call read all the annual values from a file into the CRU%CO2VALS array. On the first
    ! and subsequent

    IMPLICIT NONE

    TYPE(CRU_TYPE), INTENT(INOUT) :: CRU           ! All the info needed for CRU met runs
    REAL    :: tmparr(720,360)        ! Temporary array for reading one day of met before
    ! packing into CRU%NdepVALS(k)

    INTEGER              :: i, iunit, iyear, IOS = 0, k, t
    INTEGER :: xds, yds        ! Ndep file dimensions of long (x), lat (y)

    LOGICAL,        SAVE :: CALL1 = .TRUE.  ! A *local* variable recording the first call of this routine
    CHARACTER(200) :: NdepFILE

    ! Abbreviate dimensions for readability.
    xds = CRU%xdimsize
    yds = CRU%ydimsize

    ! For S0_TRENDY, use only static 1860 CO2 value and return immediately



    ! On the first call, allocate the CRU%CO2VALS array to store the entire history of annual CO2
    ! values, open the (ascii) CO2 file and read the values into the array.
    IF (CALL1) THEN

       NdepFILE = TRIM(CRU%BasePath)// &
            "/ndep/NOy_plus_NHx_dry_plus_wet_deposition_hist_1850_2015_annual.nc"

       ! Open the NDep and access the variables by their name and variable id.
       WRITE(*   ,*) 'Opening ndep data file: ', NdepFILE
       WRITE(logn,*) 'Opening ndep data file: ', NdepFILE


       ErrStatus = NF90_OPEN(TRIM(NdepFILE), NF90_NOWRITE, CRU%NdepF_ID)
       CALL HANDLE_ERR(ErrStatus, "Opening CRU file "//NdepFILE )
       ErrStatus = NF90_INQ_VARID(CRU%NdepF_ID,'N_deposition', CRU%NdepV_ID)
       CALL HANDLE_ERR(ErrStatus, "Inquiring CRU var "//"N_deposition"// &
            " in "//NdepFILE )

       ! Set internal counter
       CRU%Ndep_CTSTEP = 1

       IF ( TRIM(CRU%Ndep) .EQ. "static1860") THEN
          ! read Ndep at year 1860 (noting that file starts at 1850)
          CRU%Ndep_CTSTEP = 11
          t =  CRU%Ndep_CTSTEP
          ErrStatus = NF90_GET_VAR(CRU%NdepF_ID, CRU%NdepV_ID, tmparr, &
               start=(/1,1,t/),count=(/xds,yds,1/) )
          CALL HANDLE_ERR(ErrStatus, "Reading from "//NdepFILE )
          DO k = 1, CRU%mland
             CRU%NdepVALS(k) = tmparr( land_x(k), land_y(k) )
          END DO


       END IF
       CALL1 = .FALSE.
    END IF

    IF ( TRIM(CRU%Ndep) .NE. "static1860") THEN

       ! read Ndep at current year (noting that file starts at 1850 and ends in 2015)
       CRU%Ndep_CTSTEP = MIN(CRU%CYEAR, 2015) - 1850 + 1
       t =  CRU%Ndep_CTSTEP
       ErrStatus = NF90_GET_VAR(CRU%NdepF_ID, CRU%NdepV_ID, tmparr, &
            start=(/1,1,t/),count=(/xds,yds,1/) )
       CALL HANDLE_ERR(ErrStatus, "Reading from "//NdepFILE )
       DO k = 1, CRU%mland
          CRU%NdepVALS(k) = tmparr( land_x(k), land_y(k) )
       END DO

    END IF





  END SUBROUTINE GET_CRU_Ndep

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

  SUBROUTINE OPEN_CRU_MET( CRU )

    ! Opens each of the met files required for one year. This is where the distinction is made between
    ! the nominal run year (CYEAR) and the year of met required (MetYear), which is different for
    ! S0_TRENDY and S1_TRENDY than for a standard run (S2_TRENDY).

    USE cable_IO_vars_module, ONLY: timeunits ! (Char33) Name of time units read from nc file

    IMPLICIT NONE

    TYPE( CRU_TYPE ), INTENT(INOUT) :: CRU ! All CRU-NCEP related quantities and flags

    INTEGER             :: iVar            ! Loop counter through met variables
    INTEGER             :: tID             ! Numerical variable identifier returned by NetCDF routines,
    ! in this case for time. Needed to retrieve the time units.
    INTEGER             :: MetYear         ! Year of met to access. Equals CYEAR for normal runs, but
    ! must be calculated for S0_TRENDY and initialisation runs.
    INTEGER, SAVE       :: RunStartYear    ! The value of CRU%CYEAR on the first call, also equals syear.
    ! Allows the calculation of MetYear during S0_TRENDY and init runs.
    LOGICAL, SAVE       :: CALL1 = .TRUE.  ! A *local* variable recording the first call of this routine

    ! Keep the initial value of CYEAR for calculation of different MetYear if required.
    !IF (CALL1) RunStartYear = 1710 ! edit vh !
    IF (CALL1) RunStartYear = 1691 ! edit vh !
    DO iVar = 1, CRU%NMET  ! For each met variable

       ! For S0_TRENDY and initialisation, calculate the required met year for repeatedly cycling through the
       ! 30 years of 1901-1930 spinup meteorology. For normal runs 1901-2015, MetYear = CYEAR.
!    IF ( TRIM(CRU%Run) .EQ. 'S0_TRENDY' .OR.  ( TRIM(CRU%Run) .EQ. 'S1_TRENDY' )) THEN
!      MetYear = 1901 + MOD(CRU%CYEAR-RunStartYear,30)
!    ELSE IF ( TRIM(CRU%Run) .EQ. 'S2_TRENDY' ) THEN
!      MetYear = CRU%CYEAR
!    ENDIF
       IF ( TRIM(CRU%Run) .EQ. 'S0_TRENDY' .OR.  ( TRIM(CRU%Run) .EQ. 'S1_TRENDY' ) &
            .OR.  ( TRIM(CRU%Run) .EQ. 'S0_TRENDY_CO2') &
            .OR.  ( TRIM(CRU%Run) .EQ. 'S0_TRENDY_Ndep' )) THEN
          MetYear = 1901 + MOD(CRU%CYEAR-RunStartYear,30)
       ELSEIF  ( TRIM(CRU%Run) .EQ. 'S0_TRENDY_Precip' .OR. &
            TRIM(CRU%Run) .EQ. 'S0_TRENDY_CO2_Precip'.OR. &
            TRIM(CRU%Run) .EQ. 'S0_TRENDY_CO2_Temp_Precip'.OR. &
            TRIM(CRU%Run) .EQ. 'S0_TRENDY_Temp_Precip' ) THEN
          IF (iVar.EQ.1) THEN
             MetYear = CRU%CYEAR
          ELSE
             MetYear = 1901 + MOD(CRU%CYEAR-RunStartYear,30)
          ENDIF

          IF  ( TRIM(CRU%Run) .EQ. 'S0_TRENDY_CO2_Temp_Precip'.OR. &
               TRIM(CRU%Run) .EQ. 'S0_TRENDY_Temp_Precip' ) THEN
             IF (iVar.EQ.6 .OR. iVar.EQ.7) THEN
                MetYear = CRU%CYEAR
             ELSE
                MetYear = 1901 + MOD(CRU%CYEAR-RunStartYear,30)
             ENDIF
          ENDIF

       ELSEIF  ( TRIM(CRU%Run) .EQ. 'S0_TRENDY_Temp' .OR. &
            TRIM(CRU%Run) .EQ. 'S0_TRENDY_CO2_Temp' ) THEN

          IF (iVar.EQ.6 .OR. iVar.EQ.7) THEN
             MetYear = CRU%CYEAR
          ELSE
             MetYear = 1901 + MOD(CRU%CYEAR-RunStartYear,30)
          ENDIF
       ELSE IF ( TRIM(CRU%Run) .EQ. 'S2_TRENDY' ) THEN
          MetYear = CRU%CYEAR
       ELSE IF ( TRIM(CRU%Run) .EQ. 'S2_TRENDY_precip0' ) THEN
          IF (iVar.EQ.1) THEN
             ! special for baseline precip
             MetYear = 1901 + MOD(CRU%CYEAR-RunStartYear,30)
          ELSE
             MetYear = CRU%CYEAR
          ENDIF
       ELSE IF ( TRIM(CRU%Run) .EQ. 'S2_TRENDY_precip' ) THEN
          IF (iVar.NE.1) THEN
             ! special for baseline non-precip
             MetYear = 1901 + MOD(CRU%CYEAR-RunStartYear,30)
          ELSE
             MetYear = CRU%CYEAR
          ENDIF
       ENDIF

       CALL CRU_GET_FILENAME( CRU, MetYear, iVar, CRU%MetFile(iVar) ) ! Call routine to build the filenames.

       ! Open the new met files and access the variables by their name and variable id.
       WRITE(*   ,*) 'Opening met data file: ', CRU%MetFile(iVar)
       WRITE(logn,*) 'Opening met data file: ', CRU%MetFile(iVar)

       ErrStatus = NF90_OPEN(TRIM(CRU%MetFile(iVar)), NF90_NOWRITE, CRU%F_ID(iVar))
       CALL HANDLE_ERR(ErrStatus, "Opening CRU file "//CRU%MetFile(iVar) )
       ErrStatus = NF90_INQ_VARID(CRU%F_ID(iVar),TRIM(CRU%VAR_NAME(iVar)), CRU%V_ID(iVar))
       CALL HANDLE_ERR(ErrStatus, "Inquiring CRU var "//TRIM(CRU%VAR_NAME(iVar))// &
            " in "//CRU%MetFile(iVar) )
    END DO

    ! Set internal counter
    CRU%CTSTEP = 1

    CALL1 = .FALSE. ! No longer the first call (saved).

  END SUBROUTINE OPEN_CRU_MET

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

  SUBROUTINE CRU_GET_DAILY_MET( CRU, LastDayOfYear, LastYearOfMet )

    IMPLICIT NONE

    TYPE(CRU_TYPE) :: CRU
    LOGICAL, INTENT(IN)  :: LastDayOfYear, LastYearOfMet
    REAL    :: tmparr(720,360)        ! Temporary array for reading one day of met before
    ! packing into CRU%MET(iVar)%METVALS(k)
    REAL    :: tmp, stmp(365)
    INTEGER :: iVar, ii, k, x, y, realk
    INTEGER :: t, tplus1              ! The current and next timestep
    INTEGER :: fid, vid, tid          ! Netcdf id's for file, variable, and time
    INTEGER :: xds, yds, tds          ! Metfile dimensions of long (x), lat (y), and time (t)
    INTEGER :: MetYear                ! Year of meteorology currently in use
    INTEGER :: NextMetYear            ! Next met year: Where to look for the nextTmin on Dec 31st
    CHARACTER(LEN=200) :: filename

    INTEGER, SAVE       :: RunStartYear    ! The value of CRU%CYEAR on the first call, also equals syear.
    ! Allows the calculation of MetYear during S0_TRENDY and init runs.
    LOGICAL, SAVE :: CALL1 = .TRUE.   ! A *local* variable recording the first call of this routine


    ! If first call...
    ! Keep the initial value of CYEAR for calculation of different MetYear if required.
    IF (CALL1) THEN
       !RunStartYear = CRU%CYEAR
       RunStartYear = 1691
       ! If this is not the first call, capture the existing Tmax from the previous day as the
       ! 'previous Tmax' before reading another one. Move the existing next day's Tmin into the current
       ! Tmin before reading a new 'next Tmin'. These previous and next values are required for the
       ! Cesaraccio et al. algorithm used by the weather generator to interpolate daily into subdiurnal
       ! temperatures.
    ELSE
       CRU%MET(prevTmax)%METVALS(:) = CRU%MET(  Tmax  )%METVALS(:)
       CRU%MET(  Tmin  )%METVALS(:) = CRU%MET(NextTmin)%METVALS(:)
    ENDIF

    ! For S0_TRENDY and initialisation, calculate the year of meteorology as mod 50, so we repeatedly cycle
    ! through the 50 years of 1951-2000 spinup meteorology. For normal runs 1901-2015, MetYear = CYEAR.
    ! Stop with error for anything else.


!  IF ( TRIM(CRU%Run) .EQ. 'S0_TRENDY' .OR.  ( TRIM(CRU%Run) .EQ. 'S1_TRENDY' )) THEN
!    MetYear = 1901 + MOD(CRU%CYEAR-RunStartYear,30)
!  ELSE IF ( TRIM(CRU%Run) .EQ. 'S2_TRENDY' ) THEN
!    MetYear = CRU%CYEAR
!  ELSE
!    STOP 'Error in cable_cru.F90: CRU%Run not S0_TRENDY, S1_TRENDY, or 1901-2015'
!  ENDIF

    !print *, "runstartyear, metyear", runstartyear, metyear

    ! Abbreviate dimensions for readability.
    xds = CRU%xdimsize
    yds = CRU%ydimsize

    ! Loop through all 9 met variables (not including prevTmax and nextTmin, which are addressed
    ! separately as special cases of Tmax and Tmin)

    !print *,  "CRU%CTSTEP, LastDayOfYear, LastYearOfMet", CRU%CTSTEP, LastDayOfYear, LastYearOfMet
    DO iVar = 1, CRU%NMET

       IF ( TRIM(CRU%Run) .EQ. 'S0_TRENDY' .OR.  ( TRIM(CRU%Run) .EQ. 'S1_TRENDY' ) &
            .OR.  ( TRIM(CRU%Run) .EQ. 'S0_TRENDY_CO2') &
            .OR.  ( TRIM(CRU%Run) .EQ. 'S0_TRENDY_Ndep' )) THEN
          MetYear = 1901 + MOD(CRU%CYEAR-RunStartYear,30)
       ELSEIF  ( TRIM(CRU%Run) .EQ. 'S0_TRENDY_Precip' &
            .OR.  TRIM(CRU%Run) .EQ. 'S0_TRENDY_CO2_Precip'.OR. &
            TRIM(CRU%Run) .EQ. 'S0_TRENDY_CO2_Temp_Precip'.OR. &
            TRIM(CRU%Run) .EQ. 'S0_TRENDY_Temp_Precip' ) THEN
          IF (iVar.EQ.1) THEN
             MetYear = CRU%CYEAR
          ELSE
             MetYear = 1901 + MOD(CRU%CYEAR-RunStartYear,30)
          ENDIF
          IF  ( TRIM(CRU%Run) .EQ. 'S0_TRENDY_CO2_Temp_Precip'.OR. &
               TRIM(CRU%Run) .EQ. 'S0_TRENDY_Temp_Precip' ) THEN
             IF (iVar.EQ.6 .OR. iVar.EQ.7) THEN
                MetYear = CRU%CYEAR
             ELSE
                MetYear = 1901 + MOD(CRU%CYEAR-RunStartYear,30)
             ENDIF
          ENDIF

       ELSEIF  ( TRIM(CRU%Run) .EQ. 'S0_TRENDY_Temp' &
            .OR.  TRIM(CRU%Run) .EQ. 'S0_TRENDY_CO2_Temp') THEN
          IF (iVar.EQ.6 .OR. iVar.EQ.7) THEN
             MetYear = CRU%CYEAR
          ELSE
             MetYear = 1901 + MOD(CRU%CYEAR-RunStartYear,30)
          ENDIF
       ELSE IF ( TRIM(CRU%Run) .EQ. 'S2_TRENDY' ) THEN
          MetYear = CRU%CYEAR
       ELSE IF ( TRIM(CRU%Run) .EQ. 'S2_TRENDY_precip0' ) THEN
          IF (iVar.EQ.1) THEN
             ! special for baseline precip
             MetYear = 1901 + MOD(CRU%CYEAR-RunStartYear,30)
          ELSE
             MetYear = CRU%CYEAR
          ENDIF
       ELSE IF ( TRIM(CRU%Run) .EQ. 'S2_TRENDY_precip' ) THEN
          IF (iVar.NE.1) THEN
             ! special for baseline non-precip
             MetYear = 1901 + MOD(CRU%CYEAR-RunStartYear,30)
          ELSE
             MetYear = CRU%CYEAR
          ENDIF
       ENDIF


       SELECT CASE (iVar)

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

       CASE (Tmin) ! When Tmin comes up we will jump ahead to deal with nextTmin, since we
          ! are assigning Tmin from the previous value of nextTmin

          ! Bump the variable index ii to be the value for nextTmin and the time index ahead by 1 day.
          ii = nextTmin
          t  = CRU%CTSTEP
          tplus1 = CRU%CTSTEP + 1

          ! The first call exception for Tmin is to read both the first (t) and second (t+1) timesteps.
          ! The first goes into Tmin, the second into nextTmin.
          IF ( CALL1 ) THEN

             ! In the case of a small number of points, it is more efficient to pull the data for each land point
             ! directly from the file (DirectRead) rather than reading the whole array and extracting the land
             ! points from that.
             IF ( CRU%DirectRead ) THEN

                ! For each land cell, read the first day value of Tmin into tmp, and copy it to the METVALS vector.

                DO k = 1, CRU%mland

                   ! t (first value) into iVar (tmin)
                   ErrStatus = NF90_GET_VAR( CRU%F_ID(iVar), CRU%V_ID(iVar), tmp, &
                        start=(/land_x(k),land_y(k),t/) )
                   CALL HANDLE_ERR(ErrStatus, "Reading direct from "//CRU%MetFile(iVar) )
                   CRU%MET(iVar)%METVALS(k) = tmp

                   ! tplus1 (second value) into ii (nextTmin)
                   ErrStatus = NF90_GET_VAR( CRU%F_ID(iVar), CRU%V_ID(iVar), tmp, &
                        start=(/land_x(k),land_y(k),tplus1/) )
                   CALL HANDLE_ERR(ErrStatus, "Reading direct from "//CRU%MetFile(iVar) )
                   CRU%MET(ii)%METVALS(k) = tmp

                END DO

                ! Not DirectRead: Read the whole spatial grid of first day Tmin into tmparr, and copy it to the
                ! METVALS vector.
             ELSE

                ! t (first value) into iVar (tmin)
                ErrStatus = NF90_GET_VAR(CRU%F_ID(iVar), CRU%V_ID(iVar), tmparr, &
                     start=(/1,1,t/),count=(/xds,yds,1/) )
                CALL HANDLE_ERR(ErrStatus, "Reading from "//CRU%MetFile(iVar) )

                DO k = 1, CRU%mland
                   CRU%MET(iVar)%METVALS(k) = tmparr( land_x(k), land_y(k) )
                END DO

                ! tplus1 (second value) into ii (nextTmin)
                ErrStatus = NF90_GET_VAR(CRU%F_ID(iVar), CRU%V_ID(iVar), tmparr, &
                     start=(/1,1,tplus1/),count=(/xds,yds,1/) )
                CALL HANDLE_ERR(ErrStatus, "Reading from "//CRU%MetFile(iVar) )

                DO k = 1, CRU%mland
                   CRU%MET(ii)%METVALS(k) = tmparr( land_x(k), land_y(k) )
                END DO

             ENDIF  ! End of if DirectRead

          END IF   ! End of if CALL1

          ! If this is the last day of the year we will need to get the nextTmin
          ! value from the next year's met file.
          IF ( LastDayOfYear ) THEN

             IF ( LastYearOfMet ) THEN  ! If there is no more met to open...

                CONTINUE ! Do nothing. The current value of nextTmin will just be reused.

             ELSE  ! There is another year of Tmin data available

                t = 1   ! Time index is set to the first day of the next year

                ! Add one to the calculation of MetYear
!          IF ( TRIM(CRU%Run) .EQ. 'S0_TRENDY' .OR.  ( TRIM(CRU%Run) .EQ. 'S1_TRENDY' )) THEN
!            NextMetYear = 1901 + MOD(CRU%CYEAR + 1 - RunStartYear,30)
!          ELSE IF ( TRIM(CRU%Run) .EQ. 'S2_TRENDY' ) THEN
!            NextMetYear = CRU%CYEAR + 1
!          ENDIF
                IF ( TRIM(CRU%Run) .EQ. 'S0_TRENDY' .OR.  ( TRIM(CRU%Run) .EQ. 'S1_TRENDY' ) &
                     .OR.  ( TRIM(CRU%Run) .EQ. 'S0_TRENDY_CO2') &
                     .OR.  ( TRIM(CRU%Run) .EQ. 'S0_TRENDY_Ndep' )) THEN
                   NextMetYear = 1901 + MOD(CRU%CYEAR-RunStartYear,30)
                ELSEIF  ( TRIM(CRU%Run) .EQ. 'S0_TRENDY_Precip'&
                     .OR.  TRIM(CRU%Run) .EQ. 'S0_TRENDY_CO2_Precip' .OR. &
                     TRIM(CRU%Run) .EQ. 'S0_TRENDY_CO2_Temp_Precip'.OR. &
                     TRIM(CRU%Run) .EQ. 'S0_TRENDY_Temp_Precip' ) THEN
                   IF (iVar.EQ.1) THEN
                      NextMetYear = CRU%CYEAR
                   ELSE
                      NextMetYear = 1901 + MOD(CRU%CYEAR-RunStartYear,30)
                   ENDIF
                   IF  ( TRIM(CRU%Run) .EQ. 'S0_TRENDY_CO2_Temp_Precip'.OR. &
                        TRIM(CRU%Run) .EQ. 'S0_TRENDY_Temp_Precip' ) THEN
                      IF (iVar.EQ.6 .OR. iVar.EQ.7) THEN
                         NextMetYear = CRU%CYEAR
                      ELSE
                         NextMetYear = 1901 + MOD(CRU%CYEAR-RunStartYear,30)
                      ENDIF
                   ENDIF


                ELSEIF  ( TRIM(CRU%Run) .EQ. 'S0_TRENDY_Temp' &
                     .OR.  TRIM(CRU%Run) .EQ. 'S0_TRENDY_CO2_Temp' ) THEN
                   IF (iVar.EQ.6 .OR. iVar.EQ.7) THEN
                      NextMetYear = CRU%CYEAR
                   ELSE
                      NextMetYear = 1901 + MOD(CRU%CYEAR-RunStartYear,30)
                   ENDIF
                ELSE IF ( TRIM(CRU%Run) .EQ. 'S2_TRENDY' ) THEN
                   NextMetYear = CRU%CYEAR
                ELSE IF ( TRIM(CRU%Run) .EQ. 'S2_TRENDY_precip0' ) THEN
                   IF (iVar.EQ.1) THEN
                      ! special for baseline precip
                      NextMetYear = 1901 + MOD(CRU%CYEAR-RunStartYear,30)
                   ELSE
                      NextMetYear = CRU%CYEAR
                   ENDIF
                ELSE IF ( TRIM(CRU%Run) .EQ. 'S2_TRENDY_precip' ) THEN
                   IF (iVar.NE.1) THEN
                      ! special for baseline non-precip
                      NextMetYear = 1901 + MOD(CRU%CYEAR-RunStartYear,30)
                   ELSE
                      NextMetYear = CRU%CYEAR
                   ENDIF
                ENDIF

                ! Open the Tmin file for next year in preparation for reading Jan 1st...
                CALL CRU_GET_FILENAME( CRU, NextMetYear, Tmin, filename )
                ErrStatus = NF90_OPEN(TRIM(filename), NF90_NOWRITE, fid)
                CALL HANDLE_ERR(ErrStatus, "Opening CRU file "//filename )
                ErrStatus = NF90_INQ_VARID(fid,TRIM(CRU%VAR_NAME(iVar)), vid)
                CALL HANDLE_ERR(ErrStatus, "Inquiring CRU var "//filename )


                ! If DirectRead is specified (small domain) pull the points directly from the file.
                IF ( CRU%DirectRead ) THEN

                   DO k = 1, CRU%mland
                      ErrStatus = NF90_GET_VAR(fid, vid, CRU%MET(ii)%METVALS(k), &  ! Values to ii (nextTmin)
                           start=(/land_x(k), land_y(k),t/) )
                      CALL HANDLE_ERR(ErrStatus, "Reading directly from "//filename )
                   END DO

                   ! Else not DirectRead: Read full grid into temp array and extract domain from that.
                ELSE

                   ErrStatus = NF90_GET_VAR(fid, vid, tmparr, &
                        start=(/1,1,t/),count=(/xds,yds,1/) )
                   CALL HANDLE_ERR(ErrStatus, "Reading from "//filename )
                   DO k = 1, CRU%mland
                      CRU%MET(ii)%METVALS(k) = tmparr( land_x(k), land_y(k) )   ! Values to ii (nextTmin)
                   END DO

                ENDIF  ! End of If DirectRead

                ! Close the next year's Tmin met file
                ErrStatus = NF90_CLOSE(fid)
                CALL HANDLE_ERR(ErrStatus, "Closing CRU file "//filename)

             ENDIF ! End of If LastYearOfMet

          END IF  ! End of If LastDayOfYear

          ! If this is not the first call or last day of the year (i.e. somewhere in the middle)
          ! Just read the next timestep (tplus1) of Tmin data (iVar) into the nextTmin (ii) variable.

          IF ((.NOT. CALL1) .AND. (.NOT. LastDayOfYear)) THEN
             IF ( CRU%DirectRead ) THEN

                DO k = 1, CRU%mland
                   ErrStatus = NF90_GET_VAR(CRU%F_ID(iVar), CRU%V_ID(iVar), CRU%MET(ii)%METVALS(k), &
                        start=(/land_x(k),land_y(k),tplus1/) )
                   CALL HANDLE_ERR(ErrStatus, "Reading directly from "//CRU%MetFile(iVar) )
                END DO

             ELSE

                ErrStatus = NF90_GET_VAR(CRU%F_ID(iVar), CRU%V_ID(iVar), tmparr, &
                     start=(/1,1,tplus1/),count=(/xds,yds,1/) )
                CALL HANDLE_ERR(ErrStatus, "Reading from "//CRU%MetFile(iVar) )
                DO k = 1, CRU%mland
                   CRU%MET(ii)%METVALS(k) = tmparr( land_x(k), land_y(k) )
                END DO
             ENDIF
          END IF

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

       CASE (TMax)

          ! For Tmax, we need the previous day's Tmax for the Cesaraccio et al subdiurnal
          ! temperature interpolation algorithm.

          IF ( CALL1 ) THEN   ! On the first day of the run...

             ! Assign the met variable index to be prevTmax. iVar will still refer to Tmax.
             ii = prevTmax

             ! If the very first MetYear is more than 1901, there is a previous year of Tmax data available.
             ! such as during spinups, when the starting year is 1951.
             IF ( MetYear .GT. 1901 ) THEN

                ! Open the previous year's Tmax file (MetYear-1)
                CALL CRU_GET_FILENAME( CRU, MetYear-1, iVar, filename )
                ErrStatus = NF90_OPEN(TRIM(filename), NF90_NOWRITE, fid)
                CALL HANDLE_ERR(ErrStatus, "Opening CRU file "//filename )

                ! Obtain the size of the time dimension (tds), which locates the data for Dec 31st in the file.
                ErrStatus = NF90_INQ_DIMID(fid,'time',tid)
                ErrStatus = NF90_INQUIRE_DIMENSION(fid,tid,len=tds)
                CALL HANDLE_ERR(ErrStatus, "Inquiring 'time'"//TRIM(filename))

                ! Obtain the variable ID of the Tmax variable.
                ErrStatus = NF90_INQ_VARID(fid,TRIM(CRU%VAR_NAME(iVar)), vid)
                CALL HANDLE_ERR(ErrStatus, "Inquiring CRU var "//filename )

                ! If DirectRead is specified (small domain) pull the data for Dec 31st (tds) directly from the file.
                ! Otherwise, read the whole Dec 31st grid into tmparr and pack the land points from there.
                IF ( CRU%DirectRead ) THEN
                   !print *,"about to read last value of MetYear-1 ", MetYear - 1, trim(filename)
                   DO k = 1, CRU%mland
                      ErrStatus = NF90_GET_VAR(fid, vid, CRU%MET(ii)%METVALS(k), &
                           start=(/land_x(k), land_y(k),tds/) )
                      CALL HANDLE_ERR(ErrStatus, "Reading from "//filename )
                   END DO

                ELSE

                   ErrStatus = NF90_GET_VAR(fid, vid, tmparr, &
                        start=(/1,1,tds/),count=(/xds,yds,1/) )
                   CALL HANDLE_ERR(ErrStatus, "Reading from "//filename )

                   DO k = 1, CRU%mland
                      CRU%MET(ii)%METVALS(k) = tmparr( land_x(k), land_y(k) )
                   END DO

                ENDIF  ! End of If DirectRead

                ErrStatus = NF90_CLOSE(fid)
                CALL HANDLE_ERR(ErrStatus, "Closing CRU file "//filename)

                ! Having read the last day of the previous year's Tmax into prevTmax, and closed the file,
                ! now read the first day of Tmax from the current year's file (which is already open).

                t  = CRU%CTSTEP  ! CRU%CTSTEP here should be 1 (?)

                ! Read the current points directly into the vector for small domains, or read the whole grid into tmparr
                ! and extract them from there. Read the current timestep value into the iVar variable (Tmax).
                IF ( CRU%DirectRead ) THEN

                   DO k = 1, CRU%mland
                      ErrStatus = NF90_GET_VAR(CRU%F_ID(iVar), CRU%V_ID(iVar), CRU%MET(iVar)%METVALS(k), &
                           start=(/land_x(k),land_y(k),t/) )
                      CALL HANDLE_ERR(ErrStatus, "Reading directly from "//CRU%MetFile(iVar) )
                   END DO

                ELSE

                   ErrStatus = NF90_GET_VAR(CRU%F_ID(iVar), CRU%V_ID(iVar), tmparr, &
                        start=(/1,1,t/),count=(/xds,yds,1/) )
                   CALL HANDLE_ERR(ErrStatus, "Reading from "//CRU%MetFile(iVar) )
                   DO k = 1, CRU%mland
                      CRU%MET(iVar)%METVALS(k) = tmparr( land_x(k), land_y(k) )
                   END DO
                ENDIF

             ELSE ! MetYear = 1901

                ! Where MetYear is 1901, there is no previous Tmax value available, so we read the
                ! Jan 1 value into Tmax, then assign it to prevTmax as well

                ! Open the 1901 Tmax file (MetYear)
                CALL CRU_GET_FILENAME( CRU, 1901, iVar, filename )
                ErrStatus = NF90_OPEN(TRIM(filename), NF90_NOWRITE, fid)
                CALL HANDLE_ERR(ErrStatus, "Opening CRU file "//filename )

                t  = CRU%CTSTEP  ! CRU%CTSTEP here should be 1 (?)

                ! Read the current points directly into the vector for small domains, or read the whole grid into tmparr
                ! and extract them from there. Read the current timestep value into the iVar variable (Tmax).
                IF ( CRU%DirectRead ) THEN

                   DO k = 1, CRU%mland
                      ErrStatus = NF90_GET_VAR(CRU%F_ID(iVar), CRU%V_ID(iVar), CRU%MET(iVar)%METVALS(k), &
                           start=(/land_x(k),land_y(k),t/) )
                      CALL HANDLE_ERR(ErrStatus, "Reading directly from "//CRU%MetFile(iVar) )
                   END DO

                ELSE

                   ErrStatus = NF90_GET_VAR(CRU%F_ID(iVar), CRU%V_ID(iVar), tmparr, &
                        start=(/1,1,t/),count=(/xds,yds,1/) )
                   CALL HANDLE_ERR(ErrStatus, "Reading from "//CRU%MetFile(iVar) )
                   DO k = 1, CRU%mland
                      CRU%MET(iVar)%METVALS(k) = tmparr( land_x(k), land_y(k) )
                   END DO
                ENDIF

                ! Now set the value for prevTmax (ii) equal to the current value for Tmax.
                CRU%MET(ii)%METVALS(:) = CRU%MET(iVar)%METVALS(:)

             ENDIF  ! End of If MetYear > 1901

          ELSE  ! Not CALL1

             ! This is not the first call, so the only thing remaining to do for Tmax is read the new value.
             ! (Current Tmax was assigned to prevTmax at the top of the routine).

             ! Set the variable index to iVar (Tmax). The timestep is the current timestep.
             t  = CRU%CTSTEP

             ! Read the current points directly into the vector for small domains, or read the whole grid into tmparr
             ! and extract them from there.
             IF ( CRU%DirectRead ) THEN

                DO k = 1, CRU%mland
                   ErrStatus = NF90_GET_VAR(CRU%F_ID(iVar), CRU%V_ID(iVar), CRU%MET(iVar)%METVALS(k), &
                        start=(/land_x(k),land_y(k),t/) )
                   CALL HANDLE_ERR(ErrStatus, "Reading directly from "//CRU%MetFile(iVar) )
                END DO

             ELSE

                ErrStatus = NF90_GET_VAR(CRU%F_ID(iVar), CRU%V_ID(iVar), tmparr, &
                     start=(/1,1,t/),count=(/xds,yds,1/) )
                CALL HANDLE_ERR(ErrStatus, "Reading from "//CRU%MetFile(iVar) )
                DO k = 1, CRU%mland
                   CRU%MET(iVar)%METVALS(k) = tmparr( land_x(k), land_y(k) )
                END DO
             ENDIF

          END IF   ! End of If CALL1

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

       CASE DEFAULT   ! All variables other than Tmin or Tmax...

          ! iVar is not Tmin or Tmax so the variable index and timestep index are unchanged.
          ii = iVar
          t  = CRU%CTSTEP

          ! Standard read of the current variable, for the current timestep:
          ! Directly read the current points into the met vector (more efficient for small domains),
          ! or read the whole grid into tmparr and extract them from there.
          IF ( CRU%DirectRead ) THEN

             DO k = 1, CRU%mland
                ErrStatus = NF90_GET_VAR(CRU%F_ID(iVar), CRU%V_ID(iVar), CRU%MET(ii)%METVALS(k), &
                     start=(/land_x(k),land_y(k),t/) )
                CALL HANDLE_ERR(ErrStatus, "Reading directly from "//CRU%MetFile(iVar) )
             END DO

          ELSE

             ErrStatus = NF90_GET_VAR(CRU%F_ID(iVar), CRU%V_ID(iVar), tmparr, &
                  start=(/1,1,t/),count=(/xds,yds,1/) )
             CALL HANDLE_ERR(ErrStatus, "Reading from "//CRU%MetFile(iVar) )
             DO k = 1, CRU%mland
                CRU%MET(ii)%METVALS(k) = tmparr( land_x(k), land_y(k) )
             END DO
          ENDIF

       END SELECT

    END DO  ! End of loop through all met variables

    ! Convert pressure Pa -> hPa

    CRU%MET(pres)%METVALS(:) = CRU%MET(pres)%METVALS(:) / 100.

    !print *, 'CRU%CTSTEP ', CRU%CTSTEP

    ! Increment the internal timestep counter
    CRU%CTSTEP = CRU%CTSTEP + 1

    ! CALL1 can only happen once!
    IF (CALL1) CALL1 = .FALSE.

    !  print *, 'CRU%MET(tmin)%METVALS(:)'
    !  print *, CRU%MET(tmin)%METVALS(:)
    !  print *, 'CRU%MET(nexttmin)%METVALS(:)'
    !  print *, CRU%MET(nexttmin)%METVALS(:)
    !  print *, 'CRU%MET(prevtmax)%METVALS(:)'
    !  print *, CRU%MET(prevtmax)%METVALS(:)
    !  print *, 'CRU%MET(tmax)%METVALS(:)'
    !  print *, CRU%MET(tmax)%METVALS(:)
    !  print *, 'CRU%MET(pres)%METVALS(:)'
    !  print *, CRU%MET(pres)%METVALS(:)
    !  print *, 'CRU%MET(rain)%METVALS(:)'
    !  print *, CRU%MET(rain)%METVALS(:)

  END SUBROUTINE CRU_GET_DAILY_MET

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

  SUBROUTINE CRU_GET_SUBDIURNAL_MET(CRU, MET, CurYear, ktau, kend, LastYearOfMet )

    ! Obtain one day of CRU-NCEP meteorology, subdiurnalise it using a weather
    ! and return the result to the CABLE driver.

    USE cable_def_types_mod,   ONLY: MET_TYPE
    USE cable_IO_vars_module,  ONLY: LANDPT, latitude
  USE casa_ncdf_module, ONLY: DOYSOD2YMDHMS
    USE cable_weathergenerator,ONLY: WEATHER_GENERATOR_TYPE, WGEN_INIT, &
         WGEN_DAILY_CONSTANTS, WGEN_SUBDIURNAL_MET
    USE cable_checks_module,   ONLY: rh_sh

    IMPLICIT NONE

    INTEGER, INTENT(IN)  :: CurYear, ktau, kend
    LOGICAL, INTENT(IN)  :: LastYearOfMet

    TYPE(CRU_TYPE) :: CRU

    ! Define MET the CABLE version, different from the MET defined and used within the CRU variable.
    ! The structure of MET_TYPE is defined in cable_define_types.F90
    TYPE(MET_TYPE) :: MET

    ! Local variables
    LOGICAL   :: newday, LastDayOfYear  ! Flags for occurence of a new day (0 hrs) and the last day of the year.
    INTEGER   :: iland                  ! Loop counter through 'land' cells (cells in the spatial domain)
    INTEGER   :: itimestep              ! Loop counter through subdiurnal timesteps in a day
    INTEGER   :: imetvar                ! Loop counter through met variables
    INTEGER   :: dM, dD                 ! Met date as year, month, and day returned from DOYSOD2YMDHMS
    INTEGER   :: is, ie                 ! Starting and ending vegetation type per land cell
    REAL      :: dt                     ! Timestep in seconds
    REAL      :: CO2air                 ! CO2 concentration in ppm
    REAL      :: etime
    CHARACTER :: LandMaskFile*200       ! Name of the land mask file

    TYPE(WEATHER_GENERATOR_TYPE), SAVE :: WG
    LOGICAL,                      SAVE :: CALL1 = .TRUE.  ! A *local* variable recording the first call of this routine

    ! Purely for readability...
    dt = CRU%DTsecs

    ! On first step read and check CRU settings and read land-mask
    IF ( CALL1 ) CALL WGEN_INIT( WG, CRU%mland, latitude, dt )

    ! Pass time-step information to CRU
    CRU%CYEAR = CurYear


    CRU%ktau  = ktau     ! ktau is the current timestep in the year.

!  this only works with CANBERRA cable_serial, as ktau    !
!  restarts on Jan 1                                      !
    ! Based on the ktau timestep, calculate date and time information (the same for the whole spatial dimension.)
    met%hod (:) = REAL(MOD( (ktau-1) * NINT(dt), INT(SecDay)) ) / 3600.  ! Hour of the day
    met%doy (:) = INT(REAL(ktau-1) * dt / SecDay ) + 1                   ! Day of Year = days since 0 hr 1st Jan
    met%year(:) = CurYear                                                ! Current year

    ! Using the day-of-year and seconds-of-day calculate the month and day-of-month, using the time information
    ! for the first land cell only (because they will be the same across the domain).
    !                            In      In         In        Out       Out       Optional Out
    ! SUBROUTINE DOYSOD2YMDHMS( Year, Yearday, SecondsOfDay, Month, DayOfMonth, [Hour, Min, Sec])

    CALL DOYSOD2YMDHMS(CurYear, INT(met%doy(1)), INT(met%hod(1)) * 3600, dM, dD)

    met%moy (:) = dM     ! Record the month

    ! It's a new day if the hour of the day is zero.
    newday = ( met%hod(landpt(1)%cstart).EQ. 0 )

    ! Beginning-of-year accounting
    IF (ktau .EQ. 1) THEN  ! ktau is always reset to 1 at the start of the year.

       ! Read a new annual CO2 value and convert it from ppm to mol/mol
       CALL GET_CRU_CO2( CRU, CO2air )
       met%ca(:) = CO2air / 1.e+6  !

       CALL GET_CRU_Ndep( CRU )
       DO iland = 1, CRU%mland
          met%Ndep(landpt(iland)%cstart:landpt(iland)%cend) = &
               CRU%NdepVALS(iland)*86400000.  ! kg/m2/s > g/m2/d (1000.*3600.*24.)
       END DO

       ! Open a new annual CRU-NCEP met file.
       CALL OPEN_CRU_MET( CRU )

    ENDIF

    ! %%%%%% PRB to add his own comments from here down for this routine.
    ! Now get the Met data for this day

    IF ( newday ) THEN

       !print *, CRU%CTSTEP, ktau, kend
       !   CALL CPU_TIME(etime)
       !   PRINT *, 'b4 daily ', etime, ' seconds needed '

       LastDayOfYear = (ktau .EQ. kend-((SecDay/dt)-1))
       CALL CRU_GET_DAILY_MET( CRU, LastDayOfYear, LastYearOfMet )
       !   CALL CPU_TIME(etime)
       !   PRINT *, 'after daily ', etime, ' seconds needed '
       !   STOP

       ! Air pressure assumed to be constant over day
       DO iland = 1, CRU%mland
          met%pmb(landpt(iland)%cstart:landpt(iland)%cend) = CRU%MET(pres)%METVALS(iland) !CLN interpolation??
       END DO

       ! Convert wind from u and v components to wind speed by Pythagorean Theorem
       WG%WindDay        = SQRT( (CRU%MET(uwind)%METVALS * CRU%MET(uwind)%METVALS) +  &
            (CRU%MET(vwind)%METVALS * CRU%MET(vwind)%METVALS) )
       ! Convert all temperatures from K to C
       WG%TempMinDay     = CRU%MET(  Tmin  )%METVALS - 273.15
       WG%TempMaxDay     = CRU%MET(  Tmax  )%METVALS - 273.15
       WG%TempMinDayNext = CRU%MET(NextTmin)%METVALS - 273.15
       WG%TempMaxDayPrev = CRU%MET(PrevTmax)%METVALS - 273.15
       ! Convert solar radiation from J /m2/s to MJ/m2/d
       WG%SolarMJDay     = CRU%MET(  swdn  )%METVALS * 1.e-6 * SecDay ! ->[MJ/m2/d]
       ! Convert precip from mm to m/day
       WG%PrecipDay      = CRU%MET(  rain  )%METVALS  / 1000. ! ->[m/d]
       WG%SnowDay        = 0.0

       CALL WGEN_DAILY_CONSTANTS( WG, CRU%mland, INT(met%doy(1))+1 )

       ! To get the diurnal cycle for lwdn get whole day and scale with
       ! LWDN from file later

       CRU%AVG_LWDN(:) = 0.
       DO itimestep = 1, NINT(SecDay/dt)
          CALL WGEN_SUBDIURNAL_MET( WG, CRU%mland, itimestep-1 )
          CRU%AVG_LWDN = CRU%AVG_LWDN + WG%PhiLD
       END DO
       CRU%AVG_LWDN = CRU%AVG_LWDN / (SecDay/dt)
    END IF ! End of If newday

    ! Decision has been made, that first tstep of the day is at 0:01 am

    CALL WGEN_SUBDIURNAL_MET( WG, CRU%mland, NINT(met%hod(1)*3600./dt) )

    ! assign to cable variables

    ! strangely, met% is not save over Wait all in MPI...!
    !  met%pmb = CRU%MET(pres)%METVALS !CLN interpolation??

    ! Assign weather-generated data, or daily values, as required to CABLE variables.
    !
    DO iland = 1, CRU%mland
       is = landpt(iland)%cstart
       ie = landpt(iland)%cend

       met%precip    (is:ie)   = WG%Precip(iland)


       ! Cable's swdown is split into two components, visible and nir, which
       ! get half of the CRU-NCEP swdown each.
       met%fsd       (is:ie,1) = WG%PhiSD(iland) * 0.5    ! Visible
       met%fsd       (is:ie,2) = WG%PhiSD(iland) * 0.5  ! NIR

       ! Convert C to K for cable's tk
       met%tk        (is:ie)   = WG%Temp(iland) + 273.15

       met%ua        (is:ie)   = WG%Wind(iland)
       met%coszen    (is:ie)   = WG%coszen(iland)

       ! For longwave down, scale the diurnal series returned by the weather generator (WG%PhiLD(iland))
       ! to the daily value from CRU-NCEP.
       met%fld       (is:ie)   = CRU%MET(lwdn)%METVALS(iland) * WG%PhiLD(iland) / CRU%AVG_LWDN(iland)

       ! Specific humidity (qair g/g) was not sent to the weather generator. Here we assign the
       ! daily value to the whole diurnal cycle
       met%qv        (is:ie)   = CRU%MET(qair)%METVALS(iland)

       ! calculate snowfall based on total precip and air T
       !(ref Jin et al. Table II, Hyd Proc, 1999)

!    if (WG%Temp(iland) > 2.5) then
!       met%precip_sn(is:ie) = 0.0
!    elseif ((WG%Temp(iland) <= 2.5) .and. (WG%Temp(iland) > 2.0)) then
!       met%precip_sn(is:ie) = 0.6* met%precip(is:ie)
!    elseif ((WG%Temp(iland) <= 2.0) .and. (WG%Temp(iland) > 0.0)) then
!       met%precip_sn(is:ie) = (1.0 - (54.62 - 0.2 *(WG%Temp(iland) + 273.15)))* met%precip(is:ie) ! this facr can be > 1 !
!    elseif (WG%Temp(iland) <= 0.0) then
!       met%precip_sn(is:ie) = met%precip(is:ie)
!    endif

       IF (WG%Temp(iland) <= 0.0) THEN
          met%precip_sn(is:ie) = met%precip(is:ie)
       ELSE
          met%precip_sn(is:ie) = 0.0
       ENDIF
       !met%precip(is:ie) = met%precip(is:ie) -  met%precip_sn(is:ie)



    END DO

    ! initialise within canopy air temp
    met%tvair     = met%tk
    met%tvrad     = met%tk

    ! If this is the end of the year or the end of the met, close the current met files.
    !print *, "ktau, kend, LastYearOfMet as close test:", ktau, kend
    IF (ktau .EQ. kend) THEN
       DO imetvar=1, CRU%NMET
          !print *, 'Close CRU%MetFile(imetvar)', CRU%MetFile(imetvar)
          ErrStatus = NF90_CLOSE(CRU%F_ID(imetvar))
          CALL HANDLE_ERR(ErrStatus, "Closing CRU file"//CRU%MetFile(imetvar))
       END DO
    END IF

    ! CALL1 is over...
    CALL1 = .FALSE.

  END SUBROUTINE CRU_GET_SUBDIURNAL_MET

END MODULE CABLE_CRU
