! $Header$
!
MODULE simplehydrol_mod

!*******************************************************************************************
! This module contains a simple hydrology model to compute the soil water content,
! the melting and accumulation of snow as well as ice sheet "calving" (rough assumptions)
! It is especially used over land and landice surfaces when the coupling with ORCHIDEE
! is not active, and over sea ice (especially for snow) when the coupling with NEMO
! is not active.
! contact: F. Cheruy, frederique.cheruy@lmd.ipsl.fr ; E. Vignon, etienne.vignon@lmd.ipsl.fr
!*******************************************************************************************
   USE dimphy, ONLY: klon
   USE indice_sol_mod

   IMPLICIT NONE
   SAVE

! run_off_ter and run_off_lic are the runoff at the compressed grid knon for
! land and land-ice respectively
! Note: run_off_lic is used in mod_landice and therfore not private
   REAL, ALLOCATABLE, DIMENSION(:), PRIVATE    :: run_off_ter
   !$OMP THREADPRIVATE(run_off_ter)
   REAL, ALLOCATABLE, DIMENSION(:)             :: run_off_lic
   !$OMP THREADPRIVATE(run_off_lic)

! run_off_lic_0 is the runoff at land-ice a time-step earlier, on the global 1D array grid
   REAL, ALLOCATABLE, DIMENSION(:), PRIVATE    :: run_off_lic_0
   !$OMP THREADPRIVATE(run_off_lic_0)

   REAL, PRIVATE                               :: tau_calv
   !$OMP THREADPRIVATE(tau_calv)
   REAL, ALLOCATABLE, DIMENSION(:, :)           :: ffonte_global
   !$OMP THREADPRIVATE(ffonte_global)
   REAL, ALLOCATABLE, DIMENSION(:, :)           :: fqfonte_global
   !$OMP THREADPRIVATE(fqfonte_global)
   REAL, ALLOCATABLE, DIMENSION(:, :)           :: fqcalving_global
   !$OMP THREADPRIVATE(fqcalving_global)
   REAL, ALLOCATABLE, DIMENSION(:)             :: runofflic_global
   !$OMP THREADPRIVATE(runofflic_global)
#ifdef ISO
   REAL, ALLOCATABLE, DIMENSION(:, :), PRIVATE  :: xtrun_off_ter
   !$OMP THREADPRIVATE(xtrun_off_ter)
   REAL, ALLOCATABLE, DIMENSION(:, :)           :: xtrun_off_lic
   !$OMP THREADPRIVATE(xtrun_off_lic)
   REAL, ALLOCATABLE, DIMENSION(:, :), PRIVATE  :: xtrun_off_lic_0
   !$OMP THREADPRIVATE(xtrun_off_lic_0)
   REAL, ALLOCATABLE, DIMENSION(:, :, :), PRIVATE:: fxtfonte_global
   !$OMP THREADPRIVATE(fxtfonte_global)
   REAL, ALLOCATABLE, DIMENSION(:, :, :), PRIVATE:: fxtcalving_global
   !$OMP THREADPRIVATE(fxtcalving_global)
   REAL, ALLOCATABLE, DIMENSION(:, :), PRIVATE  :: xtrunofflic_global
   !$OMP THREADPRIVATE(xtrunofflic_global)
#endif

CONTAINS
!
!****************************************************************************************
   SUBROUTINE simplehydrol_init(restart_runoff)

! This subroutine allocates and initialize variables in the module.
! The variable run_off_lic_0 is initialized to the field read from
! restart file. The other variables are initialized to zero.
!
!****************************************************************************************
! Input argument
      REAL, DIMENSION(klon), INTENT(IN) :: restart_runoff

! Local variables
      INTEGER                           :: error
      CHARACTER(len=80)              :: abort_message
      CHARACTER(len=20)              :: modname = 'fonte_neige_init'

! Allocate run-off at landice and initilize with field read from restart
!****************************************************************************************

      ALLOCATE (run_off_lic_0(klon), stat=error)
      IF (error /= 0) THEN
         abort_message = 'Pb allocation run_off_lic'
         CALL abort_physic(modname, abort_message, 1)
      END IF
      run_off_lic_0(:) = restart_runoff(:)

! Allocate other variables and initilize to zero
!****************************************************************************************
      ALLOCATE (run_off_ter(klon), stat=error)
      IF (error /= 0) THEN
         abort_message = 'Pb allocation run_off_ter'
         CALL abort_physic(modname, abort_message, 1)
      END IF
      run_off_ter(:) = 0.

      ALLOCATE (run_off_lic(klon), stat=error)
      IF (error /= 0) THEN
         abort_message = 'Pb allocation run_off_lic'
         CALL abort_physic(modname, abort_message, 1)
      END IF
      run_off_lic(:) = 0.

      ALLOCATE (ffonte_global(klon, nbsrf))
      IF (error /= 0) THEN
         abort_message = 'Pb allocation ffonte_global'
         CALL abort_physic(modname, abort_message, 1)
      END IF
      ffonte_global(:, :) = 0.0

      ALLOCATE (fqfonte_global(klon, nbsrf))
      IF (error /= 0) THEN
         abort_message = 'Pb allocation fqfonte_global'
         CALL abort_physic(modname, abort_message, 1)
      END IF
      fqfonte_global(:, :) = 0.0

      ALLOCATE (fqcalving_global(klon, nbsrf))
      IF (error /= 0) THEN
         abort_message = 'Pb allocation fqcalving_global'
         CALL abort_physic(modname, abort_message, 1)
      END IF
      fqcalving_global(:, :) = 0.0

      ALLOCATE (runofflic_global(klon))
      IF (error /= 0) THEN
         abort_message = 'Pb allocation runofflic_global'
         CALL abort_physic(modname, abort_message, 1)
      END IF
      runofflic_global(:) = 0.0

! Read tau_calv
!***************
      CALL conf_interface(tau_calv)

   END SUBROUTINE simplehydrol_init
!************************************************************************************

#ifdef ISO
!************************************************************************************
   SUBROUTINE simplehydrol_init_iso(xtrestart_runoff)

! This subroutine allocates and initialize variables in the module for water isotopes.
! The variable run_off_lic_0 is initialized to the field read from
! restart file. The other variables are initialized to zero.
!************************************************************************************

      USE infotrac_phy, ONLY: niso
#ifdef ISOVERIF
      USE isotopes_mod, ONLY: iso_eau, iso_HDO
      USE isotopes_verif_mod
#endif

! Declarations
!****************************************************************************************

! Input argument
      REAL, DIMENSION(niso, klon), INTENT(IN) :: xtrestart_runoff

! Local variables
      INTEGER                           :: error
      CHARACTER(len=80)              :: abort_message
      CHARACTER(len=20)              :: modname = 'simplehydrol_init'
      INTEGER                           :: i

! Allocate run-off at landice and initilize with field read from restart
!****************************************************************************************

      ALLOCATE (xtrun_off_lic_0(niso, klon), stat=error)
      IF (error /= 0) THEN
         abort_message = 'Pb allocation run_off_lic'
         CALL abort_physic(modname, abort_message, 1)
      END IF

      xtrun_off_lic_0(:, :) = xtrestart_runoff(:, :)

#ifdef ISOVERIF
      IF (iso_eau > 0) THEN
         CALL iso_verif_egalite_vect1D( &
      &           xtrun_off_lic_0, run_off_lic_0, 'simplehydrol 100', &
      &           niso, klon)
      END IF !IF (iso_eau > 0) THEN
#endif

! Allocate other variables and initialize to zero
!****************************************************************************************

      ALLOCATE (xtrun_off_ter(niso, klon), stat=error)
      IF (error /= 0) THEN
         abort_message = 'Pb allocation xtrun_off_ter'
         CALL abort_physic(modname, abort_message, 1)
      END IF
      xtrun_off_ter(:, :) = 0.

      ALLOCATE (xtrun_off_lic(niso, klon), stat=error)
      IF (error /= 0) THEN
         abort_message = 'Pb allocation xtrun_off_lic'
         CALL abort_physic(modname, abort_message, 1)
      END IF
      xtrun_off_lic(:, :) = 0.

      ALLOCATE (fxtfonte_global(niso, klon, nbsrf))
      IF (error /= 0) THEN
         abort_message = 'Pb allocation fxtfonte_global'
         CALL abort_physic(modname, abort_message, 1)
      END IF
      fxtfonte_global(:, :, :) = 0.0

      ALLOCATE (fxtcalving_global(niso, klon, nbsrf))
      IF (error /= 0) THEN
         abort_message = 'Pb allocation fxtcalving_global'
         CALL abort_physic(modname, abort_message, 1)
      END IF
      fxtcalving_global(:, :, :) = 0.0

      ALLOCATE (xtrunofflic_global(niso, klon))
      IF (error /= 0) THEN
         abort_message = 'Pb allocation xtrunofflic_global'
         CALL abort_physic(modname, abort_message, 1)
      END IF
      xtrunofflic_global(:, :) = 0.0

   END SUBROUTINE simplehydrol_init_iso
#endif
!****************************************************************************************

!****************************************************************************************
   SUBROUTINE simplehydrol(knon, nisurf, knindex, dtime, &
                           tsurf, precip_rain, precip_snow, &
                           snow, qsol, tsurf_new, evap, ice_sub &
#ifdef ISO
                           , fq_fonte_diag, fqfonte_diag, snow_sub_diag, fqcalving_diag &
                           , max_eau_sol_diag, runoff_diag, run_off_lic_diag, coeff_rel_diag &
#endif
                           )
!$gpum horizontal knon klon
      USE indice_sol_mod
#ifdef ISO
      USE infotrac_phy, ONLY: niso
      !use isotopes_mod, ONLY: ridicule_snow,iso_eau,iso_HDO
#ifdef ISOVERIF
      USE isotopes_verif_mod
#endif
#endif
      USE yoethf_mod_h
      USE clesphys_mod_h
      USE yomcst_mod_h

!**********************************************************************************************
! This routines is a simple hydrology model to compute the soil water content,
! the melting and accumulation of snow as well as ice sheet "calving" terms (rough assumptions)
! It is especially used over land and landice surfaces when the coupling with ORCHIDEE
! is not active, and over sea ice (especially for snow above it) when the coupling with NEMO
! is not active.
! contact: F. Cheruy, frederique.cheruy@lmd.ipsl.fr ; E. Vignon, etienne.vignon@lmd.ipsl.fr
!**********************************************************************************************

      INCLUDE "FCTTRE.h"

! Declaration
!****************************************************************************************

! Input variables
!----------------
      INTEGER, INTENT(IN)                  :: knon  ! number of horizontal grid points
      INTEGER, INTENT(IN)                  :: nisurf ! index for surface type that is considered
      INTEGER, DIMENSION(knon), INTENT(IN) :: knindex ! list of horizontal indices on the native
      ! horizontal grid for the considered surface type

      REAL, INTENT(IN)                  :: dtime ! time step [s]
      REAL, DIMENSION(knon), INTENT(IN)    :: tsurf ! surface temperature [K]
      REAL, DIMENSION(knon), INTENT(IN)    :: precip_rain ! rainfall flux [kg/m2/s]
      REAL, DIMENSION(knon), INTENT(IN)    :: precip_snow ! snowfall flux [kg/m2/s]

! Input/Output variables
!-----------------------

      REAL, DIMENSION(knon), INTENT(INOUT) :: snow ! snow amount on ground [kg/m2]
      REAL, DIMENSION(knon), INTENT(INOUT) :: qsol ! amount of water in the soil [kg/m2]
      REAL, DIMENSION(knon), INTENT(INOUT) :: tsurf_new ! updated surface temperature [K]
      REAL, DIMENSION(knon), INTENT(INOUT) :: evap ! evaporation flux [kg/m2]

! Output variables
!-----------------

      REAL, DIMENSION(knon), INTENT(OUT)   :: ice_sub ! sublimation flux from ice over landice surfaces [kg/m2/s]
#ifdef ISO
      ! diagnostics for isotopes
      REAL, DIMENSION(knon), INTENT(OUT) :: fq_fonte_diag
      REAL, DIMENSION(knon), INTENT(OUT) :: fqfonte_diag
      REAL, DIMENSION(knon), INTENT(OUT) ::  snow_sub_diag
      REAL, DIMENSION(knon), INTENT(OUT) ::  fqcalving_diag
      REAL, INTENT(OUT) :: max_eau_sol_diag
      REAL, DIMENSION(knon), INTENT(OUT) ::  runoff_diag
      REAL, DIMENSION(knon), INTENT(OUT) :: run_off_lic_diag
      REAL, INTENT(OUT) :: coeff_rel_diag
#endif

! Local variables
!----------------

      INTEGER               :: i, j
      REAL                  :: fq_fonte ! quantify of snow that is melted [kg/m2]
      REAL                  :: coeff_rel
      REAL, PARAMETER       :: snow_max = 3000. ! maximum snow amount over ice sheets [kg/m2]
      REAL, PARAMETER       :: max_eau_sol = 150.0 ! maximum water amount in the soil [kg/m2]
      REAL, PARAMETER       :: chasno = 3.334E+05/(2.3867E+06*0.15) ! Latent heat of ice melting / (cp water) / tuning param=0.15
      REAL, DIMENSION(knon) :: ffonte    ! flux of energy associated with snow melting [W/m2]
      REAL, DIMENSION(knon) :: fqcalving ! flux of water associated with calving [kg/m2]
      REAL, DIMENSION(knon) :: fqfonte   ! flux of water associated with snow melting [kg/s/m2]
      REAL, DIMENSION(knon) :: d_ts      ! increment surface temperature [K]
      REAL, DIMENSION(knon) :: bil_eau_s ! water budget in soil [kg/m2/s]
      REAL, DIMENSION(knon) :: snow_sub ! snow sublimation flux [kg/m2/s]

      LOGICAL               :: is_snow_melting ! Is snow melting?

#ifdef ISO
      max_eau_sol_diag = max_eau_sol
#endif

! initial calculations
!****************************************************************************************
      coeff_rel = dtime/(tau_calv*rday)
      bil_eau_s(:) = 0.

! Snow increment snow due to precipitation and sublimation
!****************************************************************************************
      WHERE (precip_snow > 0.)
         snow = snow + (precip_snow*dtime)
      END WHERE

      snow_sub(:) = 0.
      ice_sub(:) = 0.

      IF (.NOT. ok_lic_cond) THEN
!---only positive sublimation has an impact on snow
!---note that this could create a bit of water
!---this was the default until CMIP6
!---Note that evap includes BOTH liquid water evaporation AND snow+ice sublimation
         WHERE (evap(:) > 0.)
            snow_sub(:) = MIN(snow(:)/dtime, evap(:))    !---one cannot sublimate more than the amount of snow
            snow(:) = snow(:) - snow_sub(:)*dtime         !---snow that remains on the ground
            snow(:) = MAX(0.0, snow(:))                     !---just in case
         END WHERE
      ELSE
!---now considers both positive and negative sublimation (so surface condensation) in the budget of snow
         snow_sub(:) = MIN(snow(:)/dtime, evap(:))    !---one cannot evaporate more than the amount of snow
         snow(:) = snow(:) - snow_sub(:)*dtime         !---snow that remains or deposits on the ground
         snow(:) = MAX(0.0, snow(:))                     !---just in case
      END IF

!---diagnostics of sublimation/condensation of ice over landice surfaces (when all the snow above has been sublimated)
!---in principle it should be 0 when ok_lic_cond that is when surface water condensation over landice was not allowed
      IF (nisurf == is_lic) THEN
         DO i = 1, knon
            ice_sub(i) = evap(i) - snow_sub(i)
         END DO
      END IF

!---diagnostics for isotopes
#ifdef ISO
      snow_sub_diag(:) = snow_sub(:)
      coeff_rel_diag = coeff_rel
#endif

! total water flux that goes into the soil (liquid precipitation - "liquid" evaporation)
!****************************************************************************************
      bil_eau_s(:) = (precip_rain(:)*dtime) - (evap(:) - snow_sub(:))*dtime

! Snow melting and calving (we remove the excess of snow wrt snowmax over ice sheets)
! + update surface temperature
!****************************************************************************************

      ffonte(:) = 0.0
      fqcalving(:) = 0.0
      fqfonte(:) = 0.0

      ! snow melting
      DO i = 1, knon
         ! Is snow melting?
         is_snow_melting = (snow(i) > epsfra .OR. nisurf == is_sic .OR. nisurf == is_lic) .AND. tsurf_new(i) >= RTT

         IF (is_snow_melting) THEN
            ! quantity of snow that is melted
            ! it is based on the energy conservation equation
            ! Lm*Dq = cp*DT*tuning_param (tuning_param=0.15)
            fq_fonte = MIN(MAX((tsurf_new(i) - RTT)/chasno, 0.0), snow(i))
            ! flux of energy corresponding to snow melting
            ffonte(i) = fq_fonte*RLMLT/dtime
            ! flux of water corresponding to snow melting
            fqfonte(i) = fq_fonte/dtime
            ! update of snow amount on ground
            snow(i) = MAX(0., snow(i) - fq_fonte)
            ! flux of melted water goes into the soil
            bil_eau_s(i) = bil_eau_s(i) + fq_fonte
            ! surface temperature update
            tsurf_new(i) = tsurf_new(i) - fq_fonte*chasno
            ! diag for isotopes
#ifdef ISO
            fq_fonte_diag(i) = fq_fonte
#endif

            ! snow/ice melting over ice surfaces
            IF (nisurf == is_sic .OR. nisurf == is_lic) THEN
               ! pay attention, melting over sea ice and landice
               ! is not bounded by the amount of available snow (no MIN)
               ! so when snow has been completely melted, the ice below melts
               ! which is an infinite source of water for the model
               ! BUT:
               ! when snow has been fully melted, the flux due to ice melting should be explicitly computed
               ! why are we adding the flux to that previously computed (double counting).
               ! why ffonte and tsurf_new updates are not in ok_lic_melt?
               ! why over lic and sic we impose tsurf=RTT and not over lands when snow remains?
               ! now by default, ok_lic_melt = false which means ffonte and fqfonte are not consistent
               ! moreover, imposing tsurf_new=RTT means that the update in tsurf is not consistent
               ! with the quantity of melting snow.
               !
               ! Suggestion:
               ! -  compute separately fq_fonte over lic/sic only if ok_lic_melt (lower-bound by 0 and not snow)
               ! - add an output variable ice_melt = max(0,fq_fonte - snow)/dtime to quantify the melt of ice (net water source)
               !   and update snow with the melt of snow only i.e. fq_fonte - ice_melt
               ! - remove the tsurf_new = RTT over lic and sic but implies a loss of convergence 

               fq_fonte = MAX((tsurf_new(i) - RTT)/chasno, 0.0)
               ffonte(i) = ffonte(i) + fq_fonte*RLMLT/dtime

               IF (ok_lic_melt) THEN
                  fqfonte(i) = fqfonte(i) + fq_fonte/dtime
                  bil_eau_s(i) = bil_eau_s(i) + fq_fonte
               END IF
               tsurf_new(i) = RTT
            END IF
            d_ts(i) = tsurf_new(i) - tsurf(i)
         END IF

         ! so called 'calving', if there is an excess of snow wrt snowmax
         ! it is instantaneously removed
         fqcalving(i) = MAX(0., snow(i) - snow_max)/dtime
         snow(i) = MIN(snow(i), snow_max)
      END DO
#ifdef ISO
      DO i = 1, knon
         fqcalving_diag(i) = fqcalving(i)
         fqfonte_diag(i) = fqfonte(i)
      END DO !DO i = 1, knon
#endif

! Soil water content and runoff
!****************************************************************************************
      ! over land surfaces
      IF (nisurf == is_ter) THEN
         DO i = 1, knon
            j = knindex(i)
            ! qsol update with bil_eau_s
            qsol(i) = qsol(i) + bil_eau_s(i)
            ! water that exceeds max_eau_sol feeds the runoff
            run_off_ter(j) = run_off_ter(j) + MAX(qsol(i) - max_eau_sol, 0.0)
#ifdef ISO
            runoff_diag(i) = MAX(qsol(i) - max_eau_sol, 0.0)
#endif
            qsol(i) = MIN(qsol(i), max_eau_sol)
         END DO
         ! over landice surfaces
      ELSE IF (nisurf == is_lic) THEN
         DO i = 1, knon
            j = knindex(i)
            !--temporal filtering
            run_off_lic(j) = coeff_rel*fqcalving(i) + (1.-coeff_rel)*run_off_lic_0(j)
            run_off_lic_0(j) = run_off_lic(j)
            !--add melting snow and liquid precip to runoff over ice cap
            run_off_lic(j) = run_off_lic(j) + fqfonte(i) + precip_rain(i)
         END DO
      END IF

#ifdef ISO
      DO i = 1, knon
         run_off_lic_diag(i) = run_off_lic(knindex(i))
      END DO
#endif

! Save ffonte, fqfonte and fqcalving in global arrays for each
! sub-surface separately
!****************************************************************************************
      DO i = 1, knon
         j = knindex(i)
         ffonte_global(j, nisurf) = ffonte(i)
         fqfonte_global(j, nisurf) = fqfonte(i)
         fqcalving_global(j, nisurf) = fqcalving(i)
      END DO

      IF (nisurf == is_lic) THEN
         DO i = 1, knon
            runofflic_global(knindex(i)) = run_off_lic(knindex(i))
         END DO
      END IF

   END SUBROUTINE simplehydrol
!****************************************************************************************

!****************************************************************************************
   SUBROUTINE simplehydrol_final(restart_runoff &
#ifdef ISO
                                 , xtrestart_runoff &
#endif
                                 )
!
! This subroutine returns run_off_lic_0 for later writing to restart file.
!****************************************************************************************

#ifdef ISO
      USE infotrac_phy, ONLY: niso
#ifdef ISOVERIF
      USE isotopes_mod, ONLY: iso_eau
      USE isotopes_verif_mod
#endif
#endif

      REAL, DIMENSION(klon), INTENT(OUT) :: restart_runoff
#ifdef ISO
      REAL, DIMENSION(niso, klon), INTENT(OUT) :: xtrestart_runoff
#ifdef ISOVERIF
      INTEGER :: i
#endif
#endif

! Set the output variables
      restart_runoff(:) = run_off_lic_0(:)
#ifdef ISO
      xtrestart_runoff(:, :) = xtrun_off_lic_0(:, :)
#ifdef ISOVERIF
      IF (iso_eau > 0) THEN
         DO i = 1, klon
            IF (iso_verif_egalite_nostop(run_off_lic_0(i) &
         &                              , xtrun_off_lic_0(iso_eau, i) &
         &                              , 'simplehydrol 413') &
         &      == 1) then
               WRITE (*, *) 'i=', i
               STOP
            END IF
         END DO !DO i=1,klon
      END IF !IF (iso_eau > 0) then
#endif
#endif

! Deallocation of all varaibles in the module

      IF (ALLOCATED(run_off_lic_0)) DEALLOCATE (run_off_lic_0)
      IF (ALLOCATED(run_off_ter)) DEALLOCATE (run_off_ter)
      IF (ALLOCATED(run_off_lic)) DEALLOCATE (run_off_lic)
      IF (ALLOCATED(ffonte_global)) DEALLOCATE (ffonte_global)
      IF (ALLOCATED(fqfonte_global)) DEALLOCATE (fqfonte_global)
      IF (ALLOCATED(fqcalving_global)) DEALLOCATE (fqcalving_global)
      IF (ALLOCATED(runofflic_global)) DEALLOCATE (runofflic_global)
#ifdef ISO
      IF (ALLOCATED(xtrun_off_lic_0)) DEALLOCATE (xtrun_off_lic_0)
      IF (ALLOCATED(xtrun_off_ter)) DEALLOCATE (xtrun_off_ter)
      IF (ALLOCATED(xtrun_off_lic)) DEALLOCATE (xtrun_off_lic)
      IF (ALLOCATED(fxtfonte_global)) DEALLOCATE (fxtfonte_global)
      IF (ALLOCATED(fxtcalving_global)) DEALLOCATE (fxtcalving_global)
      IF (ALLOCATED(xtrunofflic_global)) DEALLOCATE (xtrunofflic_global)
#endif

   END SUBROUTINE simplehydrol_final
!****************************************************************************************
   SUBROUTINE simplehydrol_get_vars(pctsrf, fqcalving_out, &
                                    fqfonte_out, ffonte_out, run_off_lic_out &
#ifdef ISO
                                    , fxtcalving_out, fxtfonte_out, xtrun_off_lic_out &
#endif
                                    )

! This routine cumulates ffonte, fqfonte and fqcalving respectively for
! all type of surfaces according to their fraction.
!
! This routine is called from physiq_mod before outputs' writting (histwrite)
!****************************************************************************************

      USE indice_sol_mod
#ifdef ISO
      USE infotrac_phy, ONLY: niso
#endif

! Input variables
!----------------
      REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: pctsrf ! fraction of subsurfaces [0-1]

! Output variables
!-----------------
      REAL, DIMENSION(klon), INTENT(OUT)      :: fqcalving_out ! flux of water associated with calving [kg/m2/s]
      REAL, DIMENSION(klon), INTENT(OUT)      :: fqfonte_out  ! flux of water associated with snow melting [kg/m2/s]
      REAL, DIMENSION(klon), INTENT(OUT)      :: ffonte_out   ! flux of energy associated with snow melting [W/m2]
      REAL, DIMENSION(klon), INTENT(OUT)      :: run_off_lic_out ! runoff flux [kg/m2/s]

#ifdef ISO
      REAL, DIMENSION(niso, klon), INTENT(OUT) :: fxtcalving_out
      REAL, DIMENSION(niso, klon), INTENT(OUT) :: fxtfonte_out
      REAL, DIMENSION(niso, klon), INTENT(OUT) :: xtrun_off_lic_out
      INTEGER   :: i, ixt
#endif

! Local variables
!----------------
      INTEGER   :: nisurf
!****************************************************************************************

      ffonte_out(:) = 0.0
      fqfonte_out(:) = 0.0
      fqcalving_out(:) = 0.0
#ifdef ISO
      fxtfonte_out(:, :) = 0.0
      fxtcalving_out(:, :) = 0.0
#endif

      DO nisurf = 1, nbsrf
         ffonte_out(:) = ffonte_out(:) + ffonte_global(:, nisurf)*pctsrf(:, nisurf)
         fqfonte_out(:) = fqfonte_out(:) + fqfonte_global(:, nisurf)*pctsrf(:, nisurf)
         fqcalving_out(:) = fqcalving_out(:) + fqcalving_global(:, nisurf)*pctsrf(:, nisurf)
      END DO

      run_off_lic_out(:) = runofflic_global(:)

#ifdef ISO
      DO nisurf = 1, nbsrf
         DO i = 1, klon
            DO ixt = 1, niso
               fxtfonte_out(ixt, i) = fxtfonte_out(ixt, i) + fxtfonte_global(ixt, i, nisurf)*pctsrf(i, nisurf)
               fxtcalving_out(ixt, i) = fxtcalving_out(ixt, i) + fxtcalving_global(ixt, i, nisurf)*pctsrf(i, nisurf)
            END DO
         END DO
      END DO
      xtrun_off_lic_out(:, :) = xtrunofflic_global(:, :)
#endif

   END SUBROUTINE simplehydrol_get_vars
!****************************************************************************************
!
!#ifdef ISO
!  subroutine simplehydrol_export_xtrun_off_lic_0(knon,xtrun_off_lic_0_diag)
!    use infotrac_phy, ONLY: niso
!
!    ! inputs
!    INTEGER, INTENT(IN)                      :: knon
!    real, INTENT(IN), DIMENSION(niso,klon)   :: xtrun_off_lic_0_diag
!
!    xtrun_off_lic_0(:,:)=xtrun_off_lic_0_diag(:,:)
!
!  end subroutine simplehydrol_export_xtrun_off_lic_0
!#endif

!****************************************************************************************
#ifdef ISO
   SUBROUTINE gestion_neige_besoin_varglob_simplehydrol(klon, knon, &
                                                        xtprecip_snow, xtprecip_rain, &
                                                        fxtfonte_neige, fxtcalving, &
                                                        knindex, nisurf, run_off_lic_diag, coeff_rel_diag)

      ! In this routine, we need global variables from simplehydrol_mod
      ! It must be included in simplehydrol_mod
      ! The other part of 'gestion_neige' is in insotopes_routines_mod because of circular
      ! dependencies

      USE infotrac_phy, ONLY: ntiso, niso
      USE isotopes_mod, ONLY: iso_eau
      USE indice_sol_mod
#ifdef ISOVERIF
      USE isotopes_verif_mod
#endif
      IMPLICIT NONE

      ! inputs
      INTEGER, INTENT(IN)                     :: klon, knon
      REAL, DIMENSION(ntiso, knon), INTENT(IN) :: xtprecip_snow, xtprecip_rain
      REAL, DIMENSION(niso, knon), INTENT(IN)  :: fxtfonte_neige, fxtcalving
      INTEGER, INTENT(IN)                     :: nisurf
      INTEGER, DIMENSION(knon), INTENT(IN)    :: knindex
      REAL, DIMENSION(klon), INTENT(IN)       :: run_off_lic_diag
      REAL, INTENT(IN)                        :: coeff_rel_diag

      ! locals
      INTEGER :: i, ixt, j

#ifdef ISOVERIF
      IF (nisurf == is_lic) THEN
         IF (iso_eau > 0) THEN
            DO i = 1, knon
               j = knindex(i)
               CALL iso_verif_egalite(xtrun_off_lic_0(iso_eau, j), &
         &             run_off_lic_0(j), 'gestion_neige_besoin_varglob_simplehydrol 625')
            END DO
         END IF
      END IF
#endif

! run_off_lic calculation
      IF (nisurf == is_lic) THEN

         DO i = 1, knon
            j = knindex(i)
            DO ixt = 1, niso
               xtrun_off_lic(ixt, i) = (coeff_rel_diag*fxtcalving(ixt, i)) &
          &                            + (1.-coeff_rel_diag)*xtrun_off_lic_0(ixt, j)
               xtrun_off_lic_0(ixt, j) = xtrun_off_lic(ixt, i)
               xtrun_off_lic(ixt, i) = xtrun_off_lic(ixt, i) + fxtfonte_neige(ixt, i) + xtprecip_rain(ixt, i)
            END DO !DO ixt=1,niso
#ifdef ISOVERIF
            IF (iso_eau > 0) THEN
               IF (iso_verif_egalite_choix_nostop(xtrun_off_lic(iso_eau, i), &
        &                  run_off_lic_diag(i), 'gestion_neige_besoin_varglob_simplehydrol 1201a', &
        &                  errmax, errmaxrel) == 1) THEN
                  WRITE (*, *) 'i,j=', i, j
                  WRITE (*, *) 'coeff_rel_diag=', coeff_rel_diag
                  STOP
               END IF
            END IF
#endif
         END DO
      END IF !IF (nisurf == is_lic) THEN

! Save ffonte, fqfonte and fqcalving in global arrays for each
! sub-surface separately
      DO i = 1, knon
         DO ixt = 1, niso
            fxtfonte_global(ixt, knindex(i), nisurf) = fxtfonte_neige(ixt, i)
            fxtcalving_global(ixt, knindex(i), nisurf) = fxtcalving(ixt, i)
         END DO !do ixt=1,niso
      END DO

      IF (nisurf == is_lic) THEN
         DO i = 1, knon
            DO ixt = 1, niso
               xtrunofflic_global(ixt, knindex(i)) = xtrun_off_lic(ixt, i)
            END DO ! DO ixt=1,niso
         END DO
      END IF

   END SUBROUTINE gestion_neige_besoin_varglob_simplehydrol
#endif

END MODULE simplehydrol_mod