!==============================================================================
! This source code is part of the
! Australian Community Atmosphere Biosphere Land Exchange (CABLE) model.
! This work is licensed under the CSIRO Open Source Software License
! Agreement (variation of the BSD / MIT License).
!
! You may not use this file except in compliance with this License.
! A copy of the License (CSIRO_BSD_MIT_License_v2.0_CABLE.txt) is located
! in each directory containing CABLE code.
!
! ==============================================================================
! Purpose: defines ranges to verify validity of inputs and outputs
!          checks mass balance and energy balance
!          switched on/off through namelist variables: check%*
!
! Contact: Bernard.Pak@csiro.au
!
! History: Small change to energy balance equation relative to 1.4b
!          Additional variables from 1.4b for range checking
!
! March 2014: Modifications to ebal
! and wbal (new wbal for sli) ! (Vanessa Haverd)
!==============================================================================

MODULE cable_checks_module
  ! Ranges_type in the module sets the acceptable ranges for all variables
  ! coming in or going out of the offline netcdf driver. The mass_balance
  ! and energy_balance subroutines calculate cumulative and per-timestep
  ! balances, as well as allow user to scrutinise balances in
  ! particular sections of the code - largely for diagnostics/fault finding.
  ! rh_sh - converts relative to sensible humidity if met file units require it
  !
  USE cable_IO_vars_module, ONLY: patch
  USE cable_abort_module, ONLY: range_abort
  USE cable_def_types_mod
  USE cable_common_module, ONLY: cable_user

  IMPLICIT NONE

  PRIVATE
  PUBLIC constant_check_range, check_range, ranges_type, ranges, mass_balance, energy_balance, rh_sh

  TYPE units_type
    CHARACTER(LEN=1) :: Rainf ! 's' (mm/s) or 'h' (mm/h)
    CHARACTER(LEN=1) :: PSurf  ! 'h'(hPa or mbar) or 'P'(Pa)
    CHARACTER(LEN=1) :: Tair  ! 'C' or 'K'
    CHARACTER(LEN=1) :: Qair  ! '%' or 'g' (spec hum)
    CHARACTER(LEN=1) :: CO2air ! 'p' (ppmv)
    CHARACTER(LEN=1) :: Wind ! 'm'(m/s)
  END TYPE units_type
  TYPE(units_type) :: units

  TYPE ranges_type
    REAL, DIMENSION(2) :: &
      nav_lon = [-360.0, 360.0], &
      nav_lat = [-90.0, 90.0], &
      time, &
      timestp, &
      ! possible forcing variables for CABLE
      SWdown = [0.0, 1360.0], & ! W/m^2
      LWdown = [0.0, 950.0], & ! W/m^2
      Rainf = [0.0, 0.1], & ! mm/s
      Snowf = [0.0, 0.1], & ! mm/s
      PSurf = [500.0, 1100.0], & ! mbar/hPa
      Tair = [200.0, 333.0], & ! K
      Qair = [0.0, 0.1], & ! g/g
      Tscrn = [-70.0, 70.0], & ! oC - INH
      Qscrn = [0.0, 0.1], & ! kg/kg
      CO2air = [160.0, 2000.0], & ! ppmv
      Wind = [0.0, 75.0], & ! m/s
      Wind_N = [-75.0, 75.0], & ! m/s
      Wind_E = [-75.0, 75.0], & ! m/s
      ! possible output variables
      Qmom = [-10.0, 8000.0], & ! kg/m/s2 - (INH generous range)
      Qh = [-2000.0, 2000.0], & ! W/m^2
      Qle = [-2500.0, 2500.0], & ! W/m^2
      Qg = [-4000.0, 4000.0], & ! W/m^2
      SWnet = [0.0, 1350.0], & ! W/m^2 (YP oct07)
      ! SWnet = [0.0, 1250.0],            & ! W/m^2
      LWnet = [-500.0, 510.0], & ! W/m^2
      Rnet = [-500.0, 1250.0], & ! W/m^2
      Evap = [-0.0045, 0.0045], &  ! note this is also used for snow melt !
      Ewater = [-0.0005, 0.0005], &
      ESoil = [-0.0015, 0.0015], &
      TVeg = [-0.0003, 0.0003], &
      ECanop = [-0.0003, 0.0003], &
      PotEvap = [-0.005, 0.005], &  !note should encompass Evap
      ACond = [0.0, 1.0], &
      SoilWet = [-0.4, 1.2], &
      Albedo = [0.0, 1.0], &
      visAlbedo = [0.0, 1.0], & ! vars intro for Ticket #27
      nirAlbedo = [0.0, 1.0], & ! vars intro for Ticket #27
      VegT = [213.0, 353.0], &
      SoilTemp = [213.0, 353.0], &
      SoilMoist = [0.0, 2000.0], &
      Qs = [0.0, 15.0], &
      Qsb = [0.0, 15.0], &
      DelSoilMoist = [-2000.0, 2000.0], &
      DelSWE = [-2000.0, 2000.0], &
      DelIntercept = [-100.0, 100.0], &
      SnowT = [213.0, 280.0], &
      BaresoilT = [213.0, 343.0], &
      AvgSurfT = [213.0, 333.0], &
      RadT = [200.0, 373.0], &
      ! vh_js !
      SWE = [0.0, 4000.0], &
      RootMoist = [0.0, 2000.0], &
      CanopInt = [0.0, 100.0], &
      NEE = [-70.0, 50.0], & ! umol/m2/s
      NPP = [-20.0, 75.0], & ! umol/m2/s
      GPP = [-20.0, 100.0], & ! umol/m2/s
      PAR = [-1000.0, 5000.0], & ! umol/m2/s
      AutoResp = [-50.0, 20.0], & ! umol/m2/s
      LeafResp = [-50.0, 20.0], & ! umol/m2/s
      HeteroResp = [-50.0, 20.0], & ! umol/m2/s
      HSoil = [-1000.0, 1000.0], &
      HVeg = [-1000.0, 1000.0], &
      SnowDepth = [0.0, 50.0], & ! EK nov07
      SnowMelt = [-99999.0, 999999.0], &
      Wbal = [-999999.0, 999999.0], &
      Ebal = [-999999.0, 999999.0], &
      ! vh_js !
      CanT = [213.0, 333.0], &
      Fwsoil = [0.0, 1.0], &
      ! parameters:
      albsoil = [0.0, 0.9], &
      isoil = [1.0, 30.0], &
      iveg = [1.0, 30.0], &
      bch = [2.0, 15.0], &
      latitude = [-90.0, 90.0], &
      c3 = [0.0, 1.0], & ! EK nov07
      clay = [0.0, 1.0], &
      css = [700.0, 2200.0], &
      rhosoil = [300.0, 3000.0], &
      hyds = [5.0E-7, 8.5E-3], & ! vh_js ! sep14
      rs20 = [0.0, 10.0], &
      sand = [0.0, 1.0], &
      sfc = [0.1, 0.5], &
      silt = [0.0, 1.0], &
      ssat = [0.35, 0.5], &
      sucs = [-0.8, -0.03], &
      swilt = [0.05, 0.4], &
      froot = [0.0, 1.0], &
      zse = [0.0, 5.0], &
      canst1 = [0.05, 0.15], &
      dleaf = [0.005, 0.4], &
      ejmax = [1.0E-5, 3.0E-4], &
      frac4 = [0.0, 1.0], &
      hc = [0.0, 100.0], &
      lai = [0.0, 8.0], &
      rp20 = [0.0, 10.0], &
      vbeta = [-999999.0, 999999.0], &
      g0 = [-0.5, 0.5], & ! Ticket #56 (must find better range)
      g1 = [0.0, 20.0], & ! Ticket #56 (must find better range)
      xalbnir = [0.0, 1.5], &
      meth = [0.0, 1.0], &
      za = [0.0, 150.0], &
      rpcoef = [0.05, 1.5], &
      shelrb = [1.0, 3.0], &
      vcmax = [5.0E-6, 1.5E-4], &
      xfang = [-1.0, 0.5], &
      ratecp = [0.01, 3.0], &
      ratecs = [0.01, 3.0], &
      refsbare = [0.0, 0.5], &
      taul = [0.0, 0.3], &
      refl = [0.0, 0.5], &
      tauw = [0.0, 0.1], &
      refw = [0.0, 0.5], &
      extkn = [0.0, 10.0], & ! YP oct07
      wai = [0.0, 5.0], & ! YP oct07
      vegcf = [0.0, 100.0], & ! YP oct07
      tminvj = [-20.0, 15.0], &
      tmaxvj = [-15.0, 30.0], &
      rootbeta = [0.7, 1.0], & ! YP oct07
      veg_class = [1.0, 20.0], &
      soil_class = [1.0, 20.0], &
      TotLivBiomass = [0.0, 1000.], &
      TotSoilCarb = [0.0, 1000.], &
      TotLittCarb = [0.0, 1000.], &
      Area = [0.0, 5000.], &
      !MD
      WatTable = [0.0, 1.0e10], &
      GWwb = [0.0, 1.0], &
      SatFrac = [0.0, 1.0], &
      Qrecharge = [-9999.0, 9999.0], &
      patchfrac = [0.0, 1.0], &
      ! soil
      slope = [0.0, 1.0], &
      slope_std = [0.0, 1.0], &
      GWdz = [0.0, 1.0], &
      ! ssnow
      ssdn = [0.0, 9999.0], &
      smass = [0.0, 9999.0], &
      sdepth = [0.0, 9999.0], &
      tggsn = [100.0, 300.0], &
      albsoiln = [0.0, 1.0], &
      S = [0.0, 1.5], &
      Tsoil = [-100.0, 100.0], &
      ! gw model
      gw_default = [0.0, 100000000.0], &
      default_l = [-99999.0, 9999999.0], & ! default large range
      default_s = [-99999.9, 99999.0] ! default smaller range
  END TYPE ranges_type

  TYPE(ranges_type), SAVE :: ranges

  INTERFACE check_range
    MODULE PROCEDURE :: check_range_d1
    MODULE PROCEDURE :: check_range_d2
    MODULE PROCEDURE :: check_range_d3
  END INTERFACE check_range

CONTAINS

  SUBROUTINE check_range_d1(vname, parameter_r1, parameter_range, ktau, met)

    CHARACTER(LEN=*) :: vname
    INTEGER, INTENT(IN) :: ktau
    TYPE(met_type), INTENT(IN), OPTIONAL :: met

    REAL(4), INTENT(IN) :: parameter_r1(:)
    REAL, INTENT(IN) :: parameter_range(2)

    INTEGER :: index

    DO index = 1, SIZE(parameter_r1)
      IF (parameter_r1(index) < parameter_range(1) .OR. parameter_r1(index) > parameter_range(2)) THEN
        CALL range_abort(vname, ktau, parameter_r1(index),               &
        parameter_range, index, patch(index)%latitude, patch(index)%longitude, met)
      END IF
    END DO

  END SUBROUTINE check_range_d1

  SUBROUTINE check_range_d2(var_name, var, var_range, ktau, met)

    CHARACTER(LEN=*), INTENT(IN) :: var_name
    INTEGER, INTENT(IN) :: ktau
    TYPE(met_type), INTENT(IN), OPTIONAL :: met

    REAL(4), INTENT(IN) :: var(:, :)
    REAL, INTENT(IN) :: var_range(2)

    REAL :: max_val, min_val
    INTEGER :: index(2)

    index = 0
    max_val = MAXVAL(var)
    min_val = MINVAL(var)

    IF (min_val < var_range(1)) THEN
      index = FINDLOC(var, min_val)
    END IF

    IF (max_val > var_range(2)) THEN
      index = FINDLOC(var, max_val)
    END IF

    IF (index(1) > 0) THEN
      CALL range_abort(var_name, ktau, var(index(1), index(2)), &
                       var_range, index(1), &
                       patch(index(1))%latitude, patch(index(1))%longitude, met)
    END IF

  END SUBROUTINE check_range_d2

  SUBROUTINE check_range_d3(var_name, var, var_range, ktau, met)

    CHARACTER(LEN=*), INTENT(IN) :: var_name
    INTEGER, INTENT(IN) :: ktau
    TYPE(met_type), INTENT(IN), OPTIONAL :: met

    REAL(4), INTENT(IN) :: var(:, :, :)
    REAL, INTENT(IN) :: var_range(2)
    INTEGER :: index(3)
    REAL :: max_val, min_val

    index = 0
    max_val = MAXVAL(var)
    min_val = MINVAL(var)

    IF (min_val < var_range(1)) THEN
      index = FINDLOC(var, min_val)
    END IF

    IF (max_val > var_range(2)) THEN
      index = FINDLOC(var, max_val)
    END IF

    IF (index(1) > 0) THEN
      CALL range_abort(var_name, ktau, var(index(1), index(2), index(3)), var_range, &
                       index(1), patch(index(1))%latitude, patch(index(1))%longitude, met)
    END IF

  END SUBROUTINE check_range_d3

  SUBROUTINE constant_check_range(soil, veg, ktau, met)
    ! Notes
    ! Done for uncommented types in cable_def_types_mod
    ! Commented ones in cable_def_types need to be checked
    ! Commented ones here need to be range-specified
    ! In future, it's better to have a callback for a derived type having both ranges and values
    ! This also clashes with range checks in the output module

    TYPE(soil_parameter_type), INTENT(IN) :: soil
    TYPE(veg_parameter_type), INTENT(IN) :: veg

    INTEGER, INTENT(IN) :: ktau
    TYPE(met_type), INTENT(IN), OPTIONAL :: met

    ! Soil
    !~ 1D

    !~~ Integer
    ! CALL  check_range(soil%isoilm, ranges%isoilm) ! Review: isoilm vs isoil
    ! CALL  check_range(soil%nhorizons, ranges%nhorizons)

    !~~ r1
    CALL check_range("bch", soil%bch, ranges%bch, ktau, met)
    CALL check_range("c3", soil%c3, ranges%c3, ktau, met)
    CALL check_range("clay", soil%clay, ranges%clay, ktau, met)
    CALL check_range("css", soil%css, ranges%css, ktau, met)
    ! CALL  check_range(soil%hsbh, ranges%hsbh, ktau, met)
    CALL check_range("hyds", soil%hyds, ranges%hyds, ktau, met)
    ! CALL  check_range(soil%i2bp3, ranges%i2bp3, ktau, met)
    ! CALL  check_range(soil%ibp2, ranges%ibp2, ktau, met)
    CALL check_range("rhosoil", soil%rhosoil, ranges%rhosoil, ktau, met)
    CALL check_range("sand", soil%sand, ranges%sand, ktau, met)
    CALL check_range("sfc", soil%sfc, ranges%sfc, ktau, met)
    CALL check_range("silt", soil%silt, ranges%silt, ktau, met)
    CALL check_range("ssat", soil%ssat, ranges%ssat, ktau, met)
    CALL check_range("sucs", soil%sucs, ranges%sucs, ktau, met)
    CALL check_range("swilt", soil%swilt, ranges%swilt, ktau, met)
    CALL check_range("zse", soil%zse, ranges%zse, ktau, met)
    ! CALL  check_range(soil%zshh, ranges%zshh, ktau, met)
    ! CALL  check_range(soil%soilcol, ranges%soilcol, ktau, met)
    ! CALL  check_range(soil%albsoilf, ranges%albsoilf, ktau, met) ! Review: albsoil vs adding an f ??

    ! CALL  check_range(soil%cnsd, ranges%cnsd, ktau, met)
    ! CALL  check_range(soil%pwb_min, ranges%pwb_min, ktau, met)

    ! CALL  check_range(soil%hkrz, ranges%hkrz, ktau, met)
    ! CALL  check_range(soil%zdepth, ranges%zdepth, ktau, met)
    ! CALL  check_range(soil%srf_frac_ma, ranges%srf_frac_ma, ktau, met)
    ! CALL  check_range(soil%edepth_ma, ranges%edepth_ma, ktau, met)
    ! CALL  check_range(soil%qhz_max, ranges%qhz_max, ktau, met)
    ! CALL  check_range(soil%qhz_efold, ranges%qhz_efold, ktau, met)
    ! CALL  check_range(soil%drain_dens, ranges%drain_dens, ktau, met)
    ! CALL  check_range(soil%elev, ranges%elev, ktau, met)
    ! CALL  check_range(soil%elev_std, ranges%elev_std, ktau, met)
    ! CALL  check_range(soil%slope, ranges%slope, ktau, met)
    ! CALL  check_range(soil%slope_std, ranges%slope_std, ktau, met)
    ! CALL  check_range(soil%GWsucs_vec, ranges%GWsucs_vec, ktau, met)
    ! CALL  check_range(soil%GWhyds_vec, ranges%GWhyds_vec, ktau, met)
    ! CALL  check_range(soil%GWbch_vec, ranges%GWbch_vec, ktau, met)
    ! CALL  check_range(soil%GWssat_vec, ranges%GWssat_vec, ktau, met)
    ! CALL  check_range(soil%GWwatr, ranges%GWwatr, ktau, met)
    ! CALL  check_range(soil%GWz, ranges%GWz, ktau, met)
    ! CALL  check_range(soil%GWdz, ranges%GWdz, ktau, met)
    ! CALL  check_range(soil%GWrhosoil_vec, ranges%GWrhosoil_vec, ktau, met)

    ! CALL  check_range(soil%clitt, ranges%clitt, ktau, met)
    ! CALL  check_range(soil%zeta, ranges%zeta, ktau, met)
    ! CALL  check_range(soil%fsatma, ranges%fsatma, ktau, met)

    !~ 2D
    !~~ Integer
    ! CALL  check_range(soil%ishorizon, ranges%ishorizon, ktau, met)

    !~~ Real
    ! CALL  check_range(soil%heat_cap_lower_limit, ranges%heat_cap_lower_limit, ktau, met)
    CALL check_range("albsoil", soil%albsoil, ranges%albsoil, ktau, met)

    !~~ Real r2
    !~~~ REVIEW: soil%field_vec vs ranges%field
    !~~~ without vec failing on 1x1 - zse and css
    !~~~ CALL  check_range(soil%zse_vec, ranges%zse_vec, ktau, met)
    !~~~ CALL  check_range(soil%css_vec, ranges%css_vec, ktau, met)
    !~~~ CALL  check_range(soil%cnsd_vec, ranges%cnsd, ktau, met)

    !~~~ CALL  check_range(soil%sucs_vec, ranges%sucs, ktau, met)
    !~~~ CALL  check_range(soil%hyds_vec, ranges%hyds, ktau, met)
    !~~~ CALL  check_range(soil%bch_vec, ranges%bch, ktau, met)
    !~~~ CALL  check_range(soil%clay_vec, ranges%clay, ktau, met)
    !~~~ CALL  check_range(soil%sand_vec, ranges%sand, ktau, met)
    !~~~ CALL  check_range(soil%silt_vec, ranges%silt, ktau, met)
    ! CALL  check_range(soil%org_vec, ranges%org, ktau, met)
    !~~~ CALL  check_range(soil%rhosoil_vec, ranges%rhosoil, ktau, met)
    !~~~ CALL  check_range(soil%ssat_vec, ranges%ssat, ktau, met)
    ! CALL  check_range(soil%watr, ranges%watr, ktau, met)
    !~~~ CALL  check_range(soil%sfc_vec, ranges%sfc, ktau, met)
    !~~~ CALL  check_range(soil%swilt_vec, ranges%swilt, ktau, met)

    ! Vegetation
    !~ 1D
    !~~ Integer
    CALL check_range("iveg", REAL(veg%iveg, 4), ranges%iveg, ktau, met)
    ! CALL  check_range(veg%iLU, ranges%iLU, ktau, met)

    !~~ Real

    CALL check_range("canst1", veg%canst1, ranges%canst1, ktau, met)
    CALL check_range("dleaf", veg%dleaf, ranges%dleaf, ktau, met)
    CALL check_range("ejmax", veg%ejmax, ranges%ejmax, ktau, met)
    CALL check_range("meth", veg%meth, ranges%meth, ktau, met)
    CALL check_range("frac4", veg%frac4, ranges%frac4, ktau, met)
    CALL check_range("hc", veg%hc, ranges%hc, ktau, met)
    ! CALL  check_range("vlai", veg%vlai, ranges%vlai, ktau, met)
    CALL check_range("xalbnir", veg%xalbnir, ranges%xalbnir, ktau, met)
    CALL check_range("rp20", veg%rp20, ranges%rp20, ktau, met)
    CALL check_range("rpcoef", veg%rpcoef, ranges%rpcoef, ktau, met)
    CALL check_range("rs20", veg%rs20, ranges%rs20, ktau, met)
    CALL check_range("shelrb", veg%shelrb, ranges%shelrb, ktau, met)
    CALL check_range("vegcf", veg%vegcf, ranges%vegcf, ktau, met)
    CALL check_range("tminvj", veg%tminvj, ranges%tminvj, ktau, met)
    ! CALL  check_range("toptvj", veg%toptvj, ranges%toptvj, ktau, met)
    CALL check_range("tmaxvj", veg%tmaxvj, ranges%tmaxvj, ktau, met)
    CALL check_range("vbeta", veg%vbeta, ranges%vbeta, ktau, met)
    CALL check_range("vcmax", veg%vcmax, ranges%vcmax, ktau, met)
    CALL check_range("xfang", veg%xfang, ranges%xfang, ktau, met)
    CALL check_range("extkn", veg%extkn, ranges%extkn, ktau, met)
    ! CALL  check_range("vlaimax", veg%vlaimax, ranges%vlaimax, ktau, met)
    CALL check_range("wai", veg%wai, ranges%wai, ktau, met)
    ! CALL  check_range("a1gs", veg%a1gs, ranges%a1gs, ktau, met)
    ! CALL  check_range("d0gs", veg%d0gs, ranges%d0gs, ktau, met)
    ! CALL  check_range("alpha", veg%alpha, ranges%alpha, ktau, met)
    ! CALL  check_range("convex", veg%convex, ranges%convex, ktau, met)
    ! CALL  check_range("cfrd", veg%cfrd, ranges%cfrd, ktau, met)
    ! CALL  check_range("gswmin", veg%gswmin, ranges%gswmin, ktau, met)
    ! CALL  check_range("conkc0", veg%conkc0, ranges%conkc0, ktau, met)
    ! CALL  check_range("conko0", veg%conko0, ranges%conko0, ktau, met)
    ! CALL  check_range("ekc", veg%ekc, ranges%ekc, ktau, met)
    ! CALL  check_range("eko", veg%eko, ranges%eko, ktau, met)
    CALL check_range("g0", veg%g0, ranges%g0, ktau, met)
    CALL check_range("g1", veg%g1, ranges%g1, ktau, met)

    !~~ Real r_2

    CALL check_range("rootbeta", REAL(veg%rootbeta, 4), ranges%rootbeta, ktau, met)
    ! CALL  check_range("gamma", veg%gamma, ranges%gamma, ktau, met)
    ! CALL  check_range("ZR", veg%ZR, ranges%ZR, ktau, met)
    ! CALL  check_range("F10", veg%F10, ranges%F10, ktau, met)
    ! CALL  check_range("clit", veg%clit, ranges%clit, ktau, met)

    !~ 2D

    !~~ logical,        public,        DIMENSION(:, ktau, met), POINTER        ::        deciduous                ! Review: No need of this since T/F

    !~~ Integer
    ! CALL  check_range("disturbance_interval", veg%disturbance_interval, ranges%disturbance_interval, ktau, met)

    !~~ Real
    CALL check_range("refl", veg%refl, ranges%refl, ktau, met)
    ! CALL check_range("taul", veg%taul, ranges%taul, ktau, met)
    CALL check_range("froot", veg%froot, ranges%froot, ktau, met)

    !~~ Real r_2
    ! CALL  check_range("disturbance_intensity", veg%disturbance_intensity, ranges%disturbance_intensity, ktau, met)

  END SUBROUTINE constant_check_range

  !==============================================================================
  !
  ! Name: mass_balance
  !
  ! Purpose: Calculate cumulative and per-timestep balance, as well as allow user
  !          to scrutinise balance in particular sections of the code - largely
  !          for diagnostics/fault finding.
  !
  !==============================================================================

  SUBROUTINE mass_balance(dels,ktau, ssnow,soil,canopy,met,                            &
       air,bal)
    !* Calculates the energy and water balances at each time step as well as
    ! the cumulative balance over the simulation.

    ! Input arguments
    REAL,INTENT(IN)                           :: dels        
      !! Time step length \(s\)
    INTEGER, INTENT(IN)                       :: ktau        
      !! Number of time steps since the start of the simulation \(-\)
    TYPE (soil_snow_type),INTENT(IN)          :: ssnow       
      !! Soil and snow variables data
    TYPE (soil_parameter_type),INTENT(IN)     :: soil        
      !! Soil parameters
    TYPE (canopy_type),INTENT(IN)             :: canopy      
      !! Canopy variables data
    TYPE(met_type),INTENT(IN)                 :: met         
      !! Met forcing data
    TYPE (air_type),INTENT(IN)                :: air
      !! Air variables data
    TYPE (balances_type),INTENT(INOUT)        :: bal         
      !! Water balance variables data

    ! Local variables
    REAL(r_2), DIMENSION(:),POINTER, SAVE     :: owb         ! volumetric soil moisture 
                                                             ! at previous time step
    REAL(r_2), DIMENSION(mp)                  :: delwb       ! change in soilmoisture
    ! b/w tsteps
    REAL, DIMENSION(mp)                       :: canopy_wbal !canopy water balance
    INTEGER                                   :: j, k        ! do loop counter

    IF(ktau==1) THEN
       ALLOCATE( owb(mp) )
       ! initial vlaue of soil moisture
       owb=ssnow%wbtot
    END IF

    ! Calculate change in soil moisture b/w timesteps:
    delwb = ssnow%wbtot - owb

    ! Update old soil moisture to this timestep value.
    owb = ssnow%wbtot

    ! net water into soil (precip-(change in canopy water storage)
    !  - (change in snow depth) - (surface runoff) - (deep drainage)
    !  - (evaporated water from vegetation and soil(excluding fevw, since
    !      it's included in change in canopy storage calculation))
    ! rml 28/2/11 ! BP changed rnof1+rnof2 to ssnow%runoff which also included rnof5
    ! which is used when nglacier=2 in soilsnow routines (BP feb2011)
    bal%wbal = REAL(met%precip - canopy%delwc - ssnow%snowd+ssnow%osnowd        &
         - ssnow%runoff-(canopy%fevw+canopy%fevc                                &
         + canopy%fes/ssnow%cls)*dels/air%rlam - delwb - ssnow%qrecharge)

    ! Canopy water balance: precip-change.can.storage-throughfall-evap+dew
    canopy_wbal = REAL(met%precip-canopy%delwc-canopy%through                   &
         - (canopy%fevw+MIN(canopy%fevc,0.0_r_2))*dels/air%rlam)

    IF (cable_user%soil_struc=='sli') THEN  ! vh March 2014 !
       ! delwcol includes change in soil water, pond and snowpack
       bal%wbal = canopy_wbal + REAL(canopy%through - ssnow%delwcol-ssnow%runoff &
            - ssnow%evap - MAX(canopy%fevc,0.0)*dels/air%rlam, r_2)

    ENDIF


    IF(ktau==1) THEN
       bal%wbal_tot = 0.; bal%precip_tot = 0.
       bal%rnoff_tot = 0.; bal%evap_tot = 0.
    ENDIF

    IF(ktau>10) THEN ! Avoid wobbly balances for ktau<10 pending later fix
       ! Add to accumulation variables:
       bal%wbal_tot = bal%wbal_tot + bal%wbal
       bal%precip_tot = bal%precip_tot + met%precip
       bal%rnoff_tot = bal%rnoff_tot + ssnow%rnof1 + ssnow%rnof2
       bal%evap_tot = bal%evap_tot                                              &
            + (canopy%fev+canopy%fes/ssnow%cls) * dels/air%rlam
    END IF

  END SUBROUTINE mass_balance

  !==============================================================================
  !
  ! Name: energy_balance
  !
  ! Purpose: Calculate cumulative and per-timestep balance, as well as allow user
  !          to scrutinise balance in particular sections of the code - largely
  !          for diagnostics/fault finding.
  !
  ! MODULEs used: cable_data (inherited)
  !
  !==============================================================================

  SUBROUTINE energy_balance( dels,ktau,met,rad,canopy,bal,ssnow,                    &
       SBOLTZ,EMLEAF, EMSOIL )

    ! Input arguments
    REAL,INTENT(IN)              :: dels   ! time step size
    INTEGER, INTENT(IN)              :: ktau   ! time step size
    TYPE (canopy_type),INTENT(IN)     :: canopy ! canopy variable data
    TYPE(met_type),INTENT(IN)         :: met    ! met data
    TYPE(radiation_type),INTENT(IN)   :: rad    ! met data
    TYPE (balances_type),INTENT(INOUT):: bal
    TYPE (soil_snow_type),INTENT(IN)  :: ssnow  ! soil data
    REAL, INTENT(IN)::                                                         &
         SBOLTZ,  & !Stefan-Bolzman constant
         EMLEAF,  & !leaf emissivity
         EMSOIL     !leaf emissivity

    ! vh_js ! note changes to this subroutine. Need to use ssnow%otss (not ssnow%tss) in these calculations.

    ! vh ! March 2014

    bal%Radbal = met%fsd(:,1) + met%fsd(:,2) + met%fld  - rad%albedo(:,1)*met%fsd(:,1) - rad%albedo(:,2)*met%fsd(:,2)  &
         - (emsoil*sboltz*rad%transd*ssnow%otss**4) - &
         (emleaf*sboltz*(1-rad%transd)*canopy%tv**4) &
         - canopy%fnv - canopy%fns

    !  soil energy - INH Ticket #133 corrected for consistency with %Ebal
    !  this includes the correction terms
    bal%EbalSoil =canopy%fns - canopy%fes & ! *ssnow%cls &
         & -canopy%fhs -canopy%ga

    ! canopy energy balance
    bal%Ebalveg = canopy%fnv - canopy%fev -canopy%fhv

    ! soil + canopy energy balance
    ! SW absorbed + LW absorbed - (LH+SH+ghflux) should = 0
    bal%Ebal = SUM(rad%qcan(:,:,1),2)+SUM(rad%qcan(:,:,2),2)+rad%qssabs &  ! vh ! March 2014
         & +met%fld-sboltz*emleaf*canopy%tv**4*(1-rad%transd) &
         -rad%flws*rad%transd &
         & -canopy%fev-canopy%fes & ! *ssnow%cls &
         & -canopy%fh -canopy%ga

    !REV_CORR - likely testing and offline-as-online cases only
    !need to add on the correction to the soil net radiation
    !as %fes, %fhs and %ga include the correction terms but %Ebal doesn't
    ! %fns_cor=0 in standard offline runs
    IF (cable_user%L_REV_CORR) THEN
       bal%Ebal = bal%Ebal + canopy%fns_cor
    ENDIF

    ! Add to cumulative balance:
    bal%ebal_tot = bal%ebal_tot + bal%ebal
    bal%RadbalSum = bal%RadbalSum + bal%Radbal

  END SUBROUTINE energy_balance

  !==============================================================================
  !
  ! Name: rh_sh
  !
  ! Purpose: Converts relative humidity to specific humidity
  !
  ! CALLed from: units_in
  !              get_met_data
  !
  ! CALLs: svp
  !
  !==============================================================================

  SUBROUTINE rh_sh (relHum,tk,psurf,specHum)

    ! Input arguments
    REAL, INTENT (IN)  ::                                                  &
         psurf,  & ! surface pressure (hPa)
         relHum, & ! relative humidity (%)
         tk        ! air temp (K)
    REAL, INTENT (OUT) :: specHum ! specific humidity (kg/kg)

    ! Local variables
    REAL ::                                                                &
         es,     & ! saturation vapour pressure
         ws        ! specific humidity at saturation

    es = svp (tk) ! saturation vapour pressure
    ws = 0.622 * es / (psurf - es) ! specific humidity at saturation
    specHum = (relHum/100.0) * ws ! specific humidity

  END SUBROUTINE rh_sh

  !==============================================================================
  !
  ! Name: svp
  !
  ! Purpose: Calculates saturation vapour pressure
  !
  ! CALLed from: rh_sh
  !
  !==============================================================================

  FUNCTION svp(tk) RESULT (F_Result)

    ! Local variables
    REAL ::                  &
         eilog,                   &
         ewlog,                   &
         ewlog2,                  &
         ewlog3,                  &
         ewlog4,                  &
         F_Result,                &
         temp,                    &
         tk,                      &
         toot,                    &
         toto,                    &
         tsot

    temp = tk - 273.15
    IF (temp < -20.0) THEN
       ! ice saturation
       toot = 273.15 / tk
       toto = 1. / toot
       eilog = -9.09718 * (toot-1) - 3.56654 * (LOG (toot) / LOG (10.0))        &
            + 0.876793 * (1-toto) + (LOG (6.1071) / LOG (10.0))
       F_Result = 10.0**eilog
    ELSE
       tsot = 373.15 / tk
       ewlog = -7.90298 * (tsot-1) + 5.02808 * (LOG (tsot) / LOG (10.0))
       ewlog2 = ewlog - 1.3816e-07 * (10**(11.344 * (1 - (1/tsot))) - 1)
       ewlog3 = ewlog2 + 0.0081328 * (10**(-3.49149 * (tsot-1)) - 1)
       ewlog4 = ewlog3 + (LOG (1013.246) / LOG (10.0))
       F_Result = 10.0**ewlog4
    END IF

  END FUNCTION svp
  !==============================================================================
END MODULE cable_checks_module
