! SUMMA - Structure for Unifying Multiple Modeling Alternatives
! Copyright (C) 2014-2020 NCAR/RAL; University of Saskatchewan; University of Washington
!
! This file is part of SUMMA
!
! For more information see: http://www.ral.ucar.edu/projects/summa
!
! This program is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! along with this program.  If not, see <http://www.gnu.org/licenses/>.

module check_icond4chm_module
USE nrtype

! access missing values
USE globalData,only:integerMissing  ! missing integer
USE globalData,only:realMissing     ! missing double precision number

! define modeling decisions
USE mDecisions_module,only:  &
 moisture,                   & ! moisture-based form of Richards' equation
 mixdform                      ! mixed form of Richards' equation

implicit none
private
public::check_icond4chm
contains

 ! ************************************************************************************************
 ! public subroutine check_icond4chm: read model initial conditions
 ! ************************************************************************************************
 subroutine check_icond4chm(&
                        indxGRU,                           & ! index of GRU in gru_struc
                        indxHRU,                           & ! index of HRU in gru_struc
                        progData,                          & ! model prognostic (state) variables
                        mparData,                          & ! model parameters
                        indxData,                          & ! layer index data
                        err,message)                         ! error control
 ! --------------------------------------------------------------------------------------------------------
 ! modules
 USE nrtype
 USE var_lookup,only:iLookParam                          ! variable lookup structure
 USE var_lookup,only:iLookProg                           ! variable lookup structure
 USE var_lookup,only:iLookIndex                          ! variable lookup structure
 USE globalData,only:gru_struc                           ! gru-hru mapping structures
 USE data_types,only:var_dlength                         ! actual data
 USE data_types,only:var_ilength                         ! actual data
 USE globaldata,only:iname_soil,iname_snow               ! named variables to describe the type of layer
 USE multiconst,only:&
                       LH_fus,    &                      ! latent heat of fusion                (J kg-1)
                       iden_ice,  &                      ! intrinsic density of ice             (kg m-3)
                       iden_water,&                      ! intrinsic density of liquid water    (kg m-3)
                       gravity,   &                      ! gravitational acceleration           (m s-2)
                       Tfreeze                           ! freezing point of pure water         (K)
 USE snow_utils_module,only:fracliquid                   ! compute volumetric fraction of liquid water in snow based on temperature
 USE updatState_module,only:updateSnow                   ! update snow states
 USE updatState_module,only:updateSoil                   ! update soil states
 implicit none

 ! --------------------------------------------------------------------------------------------------------
 ! variable declarations
 ! dummies
 integer(i4b),intent(in)            :: indxGRU        ! index of GRU in gru_struc
 integer(i4b),intent(in)            :: indxHRU        ! index of HRU in gru_struc
 type(var_dlength),intent(inout)    :: progData       ! prognostic vars
 type(var_dlength),intent(in)       :: mparData       ! parameters
 type(var_ilength),intent(in)       :: indxData       ! layer indexes
 integer(i4b),intent(out)           :: err            ! error code
 character(*),intent(out)           :: message        ! returned error message

 ! locals
 character(len=256)                :: cmessage       ! downstream error message

 ! temporary variables for realism checks
 integer(i4b)                      :: iLayer         ! index of model layer
 integer(i4b)                      :: iSoil          ! index of soil layer
 real(dp)                          :: fLiq           ! fraction of liquid water on the vegetation canopy (-)
 real(dp)                          :: vGn_m          ! van Genutchen "m" parameter (-)
 real(dp)                          :: tWat           ! total water on the vegetation canopy (kg m-2)
 real(dp)                          :: scalarTheta    ! liquid water equivalent of total water [liquid water + ice] (-)
 real(dp)                          :: h1,h2          ! used to check depth and height are consistent
 integer(i4b)                      :: nLayers        ! total number of layers
 real(dp)                          :: kappa          ! constant in the freezing curve function (m K-1)
 integer(i4b)                      :: nSnow          ! number of snow layers
 real(dp),parameter                :: xTol=1.e-10_dp ! small tolerance to address precision issues
 ! --------------------------------------------------------------------------------------------------------

 ! Start procedure here
 err=0; message="check_icond4chm/"

 ! --------------------------------------------------------------------------------------------------------
 ! Check that the initial conditions do not conflict with parameters, structure, etc.
 ! --------------------------------------------------------------------------------------------------------
 ! ensure the spectral average albedo is realistic
 if(progData%var(iLookPROG%scalarSnowAlbedo)%dat(1) > mparData%var(iLookPARAM%albedoMax)%dat(1)) &
    progData%var(iLookPROG%scalarSnowAlbedo)%dat(1) = mparData%var(iLookPARAM%albedoMax)%dat(1)
 if(progData%var(iLookPROG%scalarSnowAlbedo)%dat(1) < mparData%var(iLookPARAM%albedoMinWinter)%dat(1)) &
    progData%var(iLookPROG%scalarSnowAlbedo)%dat(1) = mparData%var(iLookPARAM%albedoMinWinter)%dat(1)
 ! ensure the visible albedo is realistic
 if(progData%var(iLookPROG%spectralSnowAlbedoDiffuse)%dat(1) > mparData%var(iLookPARAM%albedoMaxVisible)%dat(1)) &
    progData%var(iLookPROG%spectralSnowAlbedoDiffuse)%dat(1) = mparData%var(iLookPARAM%albedoMaxVisible)%dat(1)
 if(progData%var(iLookPROG%spectralSnowAlbedoDiffuse)%dat(1) < mparData%var(iLookPARAM%albedoMinVisible)%dat(1)) &
    progData%var(iLookPROG%spectralSnowAlbedoDiffuse)%dat(1) = mparData%var(iLookPARAM%albedoMinVisible)%dat(1)
 ! ensure the nearIR albedo is realistic
 if(progData%var(iLookPROG%spectralSnowAlbedoDiffuse)%dat(2) > mparData%var(iLookPARAM%albedoMaxNearIR)%dat(1)) &
    progData%var(iLookPROG%spectralSnowAlbedoDiffuse)%dat(2) = mparData%var(iLookPARAM%albedoMaxNearIR)%dat(1)
 if(progData%var(iLookPROG%spectralSnowAlbedoDiffuse)%dat(2) < mparData%var(iLookPARAM%albedoMinNearIR)%dat(1)) &
    progData%var(iLookPROG%spectralSnowAlbedoDiffuse)%dat(2) = mparData%var(iLookPARAM%albedoMinNearIR)%dat(1)

 ! ensure the initial conditions are consistent with the constitutive functions
 ! associate local variables with variables in the data structures
 associate(&
 ! state variables in the vegetation canopy
 scalarCanopyTemp  => progData%var(iLookPROG%scalarCanopyTemp)%dat(1)   , & ! canopy temperature
 scalarCanopyIce   => progData%var(iLookPROG%scalarCanopyIce)%dat(1)    , & ! mass of ice on the vegetation canopy (kg m-2)
 scalarCanopyLiq   => progData%var(iLookPROG%scalarCanopyLiq)%dat(1)    , & ! mass of liquid water on the vegetation canopy (kg m-2)
 ! state variables in the snow+soil domain
 mLayerTemp        => progData%var(iLookPROG%mLayerTemp)%dat            , & ! temperature (K)
 mLayerVolFracLiq  => progData%var(iLookPROG%mLayerVolFracLiq)%dat      , & ! volumetric fraction of liquid water in each snow layer (-)
 mLayerVolFracIce  => progData%var(iLookPROG%mLayerVolFracIce)%dat      , & ! volumetric fraction of ice in each snow layer (-)
 mLayerMatricHead  => progData%var(iLookPROG%mLayerMatricHead)%dat      , & ! matric head (m)
 mLayerLayerType   => indxData%var(iLookINDEX%layerType)%dat            , & ! type of layer (ix_soil or ix_snow)
 ! depth varying soil properties
 vGn_alpha         => mparData%var(iLookPARAM%vGn_alpha)%dat            , & ! van Genutchen "alpha" parameter (m-1)
 vGn_n             => mparData%var(iLookPARAM%vGn_n)%dat                , & ! van Genutchen "n" parameter (-)
 theta_sat         => mparData%var(iLookPARAM%theta_sat)%dat            , & ! soil porosity (-)
 theta_res         => mparData%var(iLookPARAM%theta_res)%dat            , & ! soil residual volumetric water content (-)
 ! snow parameters
 snowfrz_scale     => mparData%var(iLookPARAM%snowfrz_scale)%dat(1)     , & ! scaling parameter for the snow freezing curve (K-1)
 FCapil            => mparData%var(iLookPARAM%FCapil)%dat(1)              & ! fraction of pore space in tension storage (-)
 )  ! (associate local variables with model parameters)

 ! compute the constant in the freezing curve function (m K-1)
 kappa  = (iden_ice/iden_water)*(LH_fus/(gravity*Tfreeze))  ! NOTE: J = kg m2 s-2

 ! modify the liquid water and ice in the canopy
 if(scalarCanopyIce > 0._dp .and. scalarCanopyTemp > Tfreeze)then
  message=trim(message)//'canopy ice > 0 when canopy temperature > Tfreeze'
  err=20; return
 end if
 fLiq = fracliquid(scalarCanopyTemp,snowfrz_scale)  ! fraction of liquid water (-)
 tWat = scalarCanopyLiq + scalarCanopyIce           ! total water (kg m-2)
 scalarCanopyLiq = fLiq*tWat                        ! mass of liquid water on the canopy (kg m-2)
 scalarCanopyIce = (1._dp - fLiq)*tWat              ! mass of ice on the canopy (kg m-2)

 ! number of layers
 nLayers = gru_struc(indxGRU)%hruInfo(indxHRU)%nSnow + gru_struc(indxGRU)%hruInfo(indxHRU)%nSoil
 nSnow   = gru_struc(indxGRU)%hruInfo(indxHRU)%nSnow

 ! loop through all layers
 do iLayer=1,nLayers

  ! compute liquid water equivalent of total water (liquid plus ice)
  if (iLayer>nSnow) then ! soil layer = no volume expansion
   iSoil       = iLayer - nSnow
   vGn_m       = 1._dp - 1._dp/vGn_n(iSoil)
   scalarTheta = mLayerVolFracIce(iLayer) + mLayerVolFracLiq(iLayer)
  else ! snow layer = volume expansion allowed
   iSoil       = integerMissing
   vGn_m       = realMissing
   scalarTheta = mLayerVolFracIce(iLayer)*(iden_ice/iden_water) + mLayerVolFracLiq(iLayer)
  end if

  ! *****
  ! * check that the initial volumetric fraction of liquid water and ice is reasonable...
  ! *************************************************************************************
  
  select case(mlayerLayerType(iLayer))

   ! ***** snow
   case(iname_snow)
    ! (check liquid water)
    if(mLayerVolFracLiq(iLayer) < 0._dp)then; write(message,'(a,1x,i0)') trim(message)//'cannot initialize the model with volumetric fraction of liquid water < 0: layer = ',iLayer; err=20; return; end if
    if(mLayerVolFracLiq(iLayer) > 1._dp)then; write(message,'(a,1x,i0)') trim(message)//'cannot initialize the model with volumetric fraction of liquid water > 1: layer = ',iLayer; err=20; return; end if
    ! (check ice)
    if(mLayerVolFracIce(iLayer) > 0.80_dp)then; write(message,'(a,1x,i0)') trim(message)//'cannot initialize the model with volumetric fraction of ice > 0.80: layer = ',iLayer; err=20; return; end if
    if(mLayerVolFracIce(iLayer) < 0.05_dp)then; write(message,'(a,1x,i0)') trim(message)//'cannot initialize the model with volumetric fraction of ice < 0.05: layer = ',iLayer; err=20; return; end if
    ! check total water
    if(scalarTheta > 0.80_dp)then; write(message,'(a,1x,i0)') trim(message)//'cannot initialize the model with total water fraction [liquid + ice] > 0.80: layer = ',iLayer; err=20; return; end if
    if(scalarTheta < 0.05_dp)then; write(message,'(a,1x,i0)') trim(message)//'cannot initialize the model with total water fraction [liquid + ice] < 0.05: layer = ',iLayer; err=20; return; end if

   ! ***** soil
   case(iname_soil)

    ! (check liquid water)
    if(mLayerVolFracLiq(iLayer) < theta_res(iSoil)-xTol)then; write(message,'(a,1x,i0)') trim(message)//'cannot initialize the model with volumetric fraction of liquid water < theta_res: layer = ',iLayer; err=20; return; end if
    if(mLayerVolFracLiq(iLayer) > theta_sat(iSoil)+xTol)then; write(message,'(a,1x,i0)') trim(message)//'cannot initialize the model with volumetric fraction of liquid water > theta_sat: layer = ',iLayer; err=20; return; end if
    ! (check ice)
    if(mLayerVolFracIce(iLayer) < 0._dp                )then; write(message,'(a,1x,i0)') trim(message)//'cannot initialize the model with volumetric fraction of ice < 0: layer = '        ,iLayer; err=20; return; end if
    if(mLayerVolFracIce(iLayer) > theta_sat(iSoil)+xTol)then; write(message,'(a,1x,i0)') trim(message)//'cannot initialize the model with volumetric fraction of ice > theta_sat: layer = ',iLayer; err=20; return; end if
    ! check total water
    if(scalarTheta < theta_res(iSoil)-xTol)then; write(message,'(a,1x,i0)') trim(message)//'cannot initialize the model with total water fraction [liquid + ice] < theta_res: layer = ',iLayer; err=20; return; end if
    if(scalarTheta > theta_sat(iSoil)+xTol)then; write(message,'(a,1x,i0)') trim(message)//'cannot initialize the model with total water fraction [liquid + ice] > theta_sat: layer = ',iLayer; err=20; return; end if

   case default
    write(*,*) 'Cannot recognize case in initial vol water/ice check: type=', mlayerLayerType(iLayer)
    err=20; message=trim(message)//'cannot identify layer type'; return
  end select

  ! *****
  ! * check that the initial conditions are consistent with the constitutive functions...
  ! *************************************************************************************
  select case(mLayerLayerType(iLayer))

   ! ** snow
   case(iname_snow)

    ! check that snow temperature is less than freezing
    if(mLayerTemp(iLayer) > Tfreeze)then
     message=trim(message)//'initial snow temperature is greater than freezing'
     err=20; return
    end if

    ! ensure consistency among state variables
    call updateSnow(&
                    ! input
                    mLayerTemp(iLayer),                        & ! intent(in): temperature (K)
                    scalarTheta,                               & ! intent(in): mass fraction of total water (-)
                    snowfrz_scale,                             & ! intent(in): scaling parameter for the snow freezing curve (K-1)
                    ! output
                    mLayerVolFracLiq(iLayer),                  & ! intent(out): volumetric fraction of liquid water (-)
                    mLayerVolFracIce(iLayer),                  & ! intent(out): volumetric fraction of ice (-)
                    fLiq,                                      & ! intent(out): fraction of liquid water (-)
                    err,cmessage)                                ! intent(out): error control
    if(err/=0)then; message=trim(message)//trim(cmessage); return; end if  ! (check for errors)

   ! ** soil
   case(iname_soil)

    ! ensure consistency among state variables
    call updateSoil(&
                    ! input
                    mLayerTemp(iLayer),                                                    & ! intent(in): layer temperature (K)
                    mLayerMatricHead(iLayer-nSnow),                                        & ! intent(in): matric head (m)
                    vGn_alpha(iSoil),vGn_n(iSoil),theta_sat(iSoil),theta_res(iSoil),vGn_m, & ! intent(in): van Genutchen soil parameters
                    ! output
                    scalarTheta,                                                           & ! intent(out): volumetric fraction of total water (-)
                    mLayerVolFracLiq(iLayer),                                              & ! intent(out): volumetric fraction of liquid water (-)
                    mLayerVolFracIce(iLayer),                                              & ! intent(out): volumetric fraction of ice (-)
                    err,cmessage)                                                            ! intent(out): error control
    if(err/=0)then; message=trim(message)//trim(cmessage); return; end if  ! (check for errors)

   case default; err=10; message=trim(message)//'unknown case for model layer'; return
  end select

 end do  ! (looping through layers)

   ! end association to variables in the data structures
 end associate

 ! if snow layers exist, compute snow depth and SWE
 if(nSnow > 0)then
  progData%var(iLookPROG%scalarSWE)%dat(1)       = sum( (progData%var(iLookPROG%mLayerVolFracLiq)%dat(1:nSnow)*iden_water + &
                                                                             progData%var(iLookPROG%mLayerVolFracIce)%dat(1:nSnow)*iden_ice)        * &
                                                                             progData%var(iLookPROG%mLayerDepth)%dat(1:nSnow) )
 end if  ! if snow layers exist

 ! check that the layering is consistent
 do iLayer=1,nLayers
  h1 = sum(progData%var(iLookPROG%mLayerDepth)%dat(1:iLayer)) ! sum of the depths up to the current layer
  h2 = progData%var(iLookPROG%iLayerHeight)%dat(iLayer) - progData%var(iLookPROG%iLayerHeight)%dat(0)  ! difference between snow-atm interface and bottom of layer
  if(abs(h1 - h2) > 1.e-6_dp)then
   write(message,'(a,1x,i0)') trim(message)//'mis-match between layer depth and layer height; layer = ', iLayer, '; sum depths = ',h1,'; height = ',h2
   err=20; return
  end if
 end do
 end subroutine check_icond4chm

end module check_icond4chm_module
