
! $Header$

MODULE fonte_neige_mod

! This module will treat the process of snow, melting, accumulating, calving, in 
! case of simplified soil model.

!****************************************************************************************
  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 fonte_neige_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)
    ENDIF
    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)
    ENDIF
    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)
    ENDIF
    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)
    ENDIF
    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)
    ENDIF
    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)
    ENDIF
    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)
    ENDIF
    runofflic_global(:) = 0.0

!****************************************************************************************
! Read tau_calv

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


  END SUBROUTINE fonte_neige_init

#ifdef ISO
  SUBROUTINE fonte_neige_init_iso(xtrestart_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.

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

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

! Local variables
    INTEGER                           :: error
    CHARACTER (len = 80)              :: abort_message 
    CHARACTER (len = 20)              :: modname = 'fonte_neige_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_gcm(modname,abort_message,1)
    ENDIF    
    
    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,'fonte_neige 100', &
             niso,klon)
      ENDIF !IF (iso_eau > 0) THEN
#endif        

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

!****************************************************************************************

    ALLOCATE(xtrun_off_ter(niso,klon), stat = error)
    IF (error /= 0) THEN
       abort_message='Pb allocation xtrun_off_ter'
       CALL abort_gcm(modname,abort_message,1)
    ENDIF
    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_gcm(modname,abort_message,1)
    ENDIF
    xtrun_off_lic(:,:) = 0.

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

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

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

  END SUBROUTINE fonte_neige_init_iso
#endif

!****************************************************************************************

  SUBROUTINE fonte_neige( knon, nisurf, knindex, dtime, &
       tsurf, precip_rain, precip_snow, &
       snow, qsol, tsurf_new, evap &
#ifdef ISO    
   ,fq_fonte_diag,fqfonte_diag,snow_evap_diag,fqcalving_diag   &
   ,max_eau_sol_diag,runoff_diag,run_off_lic_diag,coeff_rel_diag   &
#endif
     )

    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
        
! Routine de traitement de la fonte de la neige dans le cas du traitement
! de sol simplifie!
! LF 03/2001
! input:
!   knon         nombre de points a traiter
!   nisurf       surface a traiter
!   knindex      index des mailles valables pour surface a traiter
!   dtime        
!   tsurf        temperature de surface
!   precip_rain  precipitations liquides
!   precip_snow  precipitations solides

! input/output:
!   snow         champs hauteur de neige
!   qsol         hauteur d'eau contenu dans le sol
!   tsurf_new    temperature au sol
!   evap

  INCLUDE "YOETHF.h"
  INCLUDE "YOMCST.h"
  INCLUDE "FCTTRE.h"
  INCLUDE "clesphys.h"

! Input variables
!****************************************************************************************
    INTEGER, INTENT(IN)                  :: knon
    INTEGER, INTENT(IN)                  :: nisurf
    INTEGER, DIMENSION(klon), INTENT(IN) :: knindex
    REAL   , INTENT(IN)                  :: dtime
    REAL, DIMENSION(klon), INTENT(IN)    :: tsurf
    REAL, DIMENSION(klon), INTENT(IN)    :: precip_rain
    REAL, DIMENSION(klon), INTENT(IN)    :: precip_snow

    ! Input/Output variables
!****************************************************************************************

    REAL, DIMENSION(klon), INTENT(INOUT) :: snow
    REAL, DIMENSION(klon), INTENT(INOUT) :: qsol
    REAL, DIMENSION(klon), INTENT(INOUT) :: tsurf_new
    REAL, DIMENSION(klon), INTENT(INOUT) :: evap

#ifdef ISO    
        ! sortie de quelques diagnostiques
    REAL, DIMENSION(klon), INTENT(OUT) :: fq_fonte_diag
    REAL, DIMENSION(klon), INTENT(OUT) :: fqfonte_diag
    REAL, DIMENSION(klon), INTENT(OUT) ::  snow_evap_diag 
    REAL, DIMENSION(klon), INTENT(OUT) ::  fqcalving_diag  
    REAL,                  INTENT(OUT) :: max_eau_sol_diag  
    REAL, DIMENSION(klon), INTENT(OUT) ::  runoff_diag   
    REAL, DIMENSION(klon), INTENT(OUT) :: run_off_lic_diag  
    REAL,                  INTENT(OUT) :: coeff_rel_diag
#endif

! Local variables
!****************************************************************************************

    INTEGER               :: i, j
    REAL                  :: fq_fonte
    REAL                  :: coeff_rel
    REAL, PARAMETER       :: snow_max=3000.
    REAL, PARAMETER       :: max_eau_sol = 150.0
!! PB temporaire en attendant mieux pour le modele de neige
! REAL, parameter :: chasno = RLMLT/(2.3867E+06*0.15)
    REAL, PARAMETER       :: chasno = 3.334E+05/(2.3867E+06*0.15)
!IM cf JLD/ GKtest
    REAL, PARAMETER       :: chaice = 3.334E+05/(2.3867E+06*0.15)
! fin GKtest
    REAL, DIMENSION(klon) :: ffonte
    REAL, DIMENSION(klon) :: fqcalving, fqfonte
    REAL, DIMENSION(klon) :: d_ts
    REAL, DIMENSION(klon) :: bil_eau_s, snow_evap

    LOGICAL               :: neige_fond

#ifdef ISO
        max_eau_sol_diag=max_eau_sol
#endif


!****************************************************************************************
! Start calculation
! - Initialization

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

!****************************************************************************************
! - Increment snow due to precipitation and evaporation
! - Calculate the water balance due to precipitation and evaporation (bil_eau_s)

!****************************************************************************************
    WHERE (precip_snow > 0.) 
       snow = snow + (precip_snow * dtime)
    END WHERE

    snow_evap = 0.
  
    IF (.NOT. ok_lic_cond) THEN
!---only positive evaporation has an impact on snow 
!---note that this could create a bit of water
!---this was the default until CMIP6 
      WHERE (evap > 0. )
         snow_evap = MIN (snow / dtime, evap)    !---one cannot evaporate more than the amount of snow
         snow = snow - snow_evap * dtime         !---snow that remains on the ground
         snow = MAX(0.0, snow)                   !---just in case
      END WHERE
    ELSE
!--now considers both positive and negative evaporation in the budget of snow 
      snow_evap = MIN (snow / dtime, evap)    !---one cannot evaporate more than the amount of snow
      snow = snow - snow_evap * dtime         !---snow that remains or deposits on the ground
      snow = MAX(0.0, snow)                   !---just in case
   ENDIF
    
    bil_eau_s(:) = (precip_rain(:) * dtime) - (evap(:) - snow_evap(:)) * dtime
#ifdef ISO
    snow_evap_diag(:) = snow_evap(:)
    coeff_rel_diag    = coeff_rel
#endif



!****************************************************************************************
! - Calculate melting snow
! - Calculate calving and decrement snow, if there are to much snow
! - Update temperature at surface

!****************************************************************************************

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

    DO i = 1, knon
       ! Y'a-t-il fonte de neige?
       neige_fond = (snow(i)>epsfra .OR. nisurf==is_sic .OR. nisurf==is_lic) .AND. tsurf_new(i)>=RTT
       IF (neige_fond) THEN
          fq_fonte     = MIN( MAX((tsurf_new(i)-RTT )/chasno,0.0),snow(i))
          ffonte(i)    = fq_fonte * RLMLT/dtime
          fqfonte(i)   = fq_fonte/dtime
          snow(i)      = MAX(0., snow(i) - fq_fonte)
          bil_eau_s(i) = bil_eau_s(i) + fq_fonte 
          tsurf_new(i) = tsurf_new(i) - fq_fonte * chasno  
#ifdef ISO
          fq_fonte_diag(i) = fq_fonte
#endif


!IM cf JLD OK     
!IM cf JLD/ GKtest fonte aussi pour la glace
          IF (nisurf == is_sic .OR. nisurf == is_lic ) THEN
             fq_fonte = MAX((tsurf_new(i)-RTT )/chaice,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
             ENDIF
             tsurf_new(i) = RTT
          ENDIF
          d_ts(i) = tsurf_new(i) - tsurf(i)
       ENDIF

       ! s'il y a une hauteur trop importante de neige, elle est ecretee
       fqcalving(i) = MAX(0., snow(i) - snow_max)/dtime
       snow(i)=MIN(snow(i),snow_max)
    ENDDO
#ifdef ISO
    DO i = 1, knon
       fqcalving_diag(i) = fqcalving(i)
       fqfonte_diag(i)   = fqfonte(i)
    ENDDO !DO i = 1, knon
#endif


    IF (nisurf == is_ter) THEN
       DO i = 1, knon
          qsol(i) = qsol(i) + bil_eau_s(i)
          run_off_ter(i) = run_off_ter(i) + 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) 
       ENDDO
    ELSE IF (nisurf == is_lic) THEN
       DO i = 1, knon
          j = knindex(i)
          !--temporal filtering
          run_off_lic(i)   = coeff_rel*fqcalving(i) + (1.-coeff_rel)*run_off_lic_0(j)
          run_off_lic_0(j) = run_off_lic(i)
          !--add melting snow and liquid precip to runoff of ice cap
          run_off_lic(i)   = run_off_lic(i) + fqfonte(i) + precip_rain(i)
       ENDDO
    ENDIF

#ifdef ISO
    DO i = 1, klon    
      run_off_lic_diag(i) = run_off_lic(i)
    ENDDO ! DO i = 1, knon    
#endif
    
!****************************************************************************************
! Save ffonte, fqfonte and fqcalving in global arrays for each 
! sub-surface separately

!****************************************************************************************
    DO i = 1, knon
       ffonte_global(knindex(i),nisurf)    = ffonte(i)
       fqfonte_global(knindex(i),nisurf)   = fqfonte(i)
       fqcalving_global(knindex(i),nisurf) = fqcalving(i)
    ENDDO

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

  END SUBROUTINE fonte_neige

!****************************************************************************************

  SUBROUTINE fonte_neige_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) &
                                ,'fonte_neige 413') &
        == 1) then
          WRITE(*,*) 'i=',i
          STOP
        ENDIF
      ENDDO !DO i=1,klon
    ENDIF !IF (iso_eau > 0) then 
#endif    
#endif



! Deallocation of all varaibles in the module
!   DEALLOCATE(run_off_lic_0, run_off_ter, run_off_lic, ffonte_global, &
!        fqfonte_global, fqcalving_global)

    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 fonte_neige_final

!****************************************************************************************

  SUBROUTINE fonte_neige_get_vars(pctsrf, fqcalving_out, &
              fqfonte_out, ffonte_out, run_off_lic_out &
#ifdef ISO     
         ,fxtcalving_out, fxtfonte_out,xtrun_off_lic_out &
#endif     
         )


! Cumulate ffonte, fqfonte and fqcalving respectively for
! all type of surfaces according to their fraction.

! This routine is called from physiq.F before histwrite.
!****************************************************************************************

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

    REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: pctsrf

    REAL, DIMENSION(klon), INTENT(OUT)      :: fqcalving_out
    REAL, DIMENSION(klon), INTENT(OUT)      :: fqfonte_out
    REAL, DIMENSION(klon), INTENT(OUT)      :: ffonte_out
    REAL, DIMENSION(klon), INTENT(OUT)      :: run_off_lic_out

#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
 
    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)
    ENDDO

    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)
        ENDDO !DO ixt=1,niso
      ENDDO !DO i=1,klon
    ENDDO !DO nisurf = 1, nbsrf
    xtrun_off_lic_out(:,:) = xtrunofflic_global(:,:)
#endif

  END SUBROUTINE fonte_neige_get_vars

!****************************************************************************************

!#ifdef ISO
!  SUBROUTINE fonte_neige_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  fonte_neige_export_xtrun_off_lic_0
!#endif 

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

        ! dans cette routine, on a besoin des variables globales de
        ! fonte_neige_mod. Il faut donc la mettre dans fonte_neige_mod
        ! le reste de gestion_neige est dans isotopes_routines_mod car sinon pb
        ! de dépendance circulaire.

    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,klon), INTENT(IN) :: xtprecip_snow, xtprecip_rain
    REAL, DIMENSION(niso,klon), INTENT(IN)  :: fxtfonte_neige,fxtcalving
    INTEGER, INTENT(IN)                     :: nisurf
    INTEGER, DIMENSION(klon), 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_fonte_neige 625')
        ENDDO
      ENDIF
    ENDIF
#endif

! calcul de run_off_lic 

    IF (nisurf == is_lic) THEN
!         coeff_rel = dtime/(tau_calv * rday)

      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)
        ENDDO !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_fonte_neige 1201a', &
                    errmax,errmaxrel) == 1) THEN
               WRITE(*,*) 'i,j=',i,j   
               WRITE(*,*) 'coeff_rel_diag=',coeff_rel_diag
               STOP
            ENDIF
          ENDIF
#endif
      ENDDO
    ENDIF !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)
      ENDDO !do ixt=1,niso
    ENDDO   

    IF (nisurf == is_lic) THEN
      DO i = 1, knon    
        DO ixt = 1, niso   
        xtrunofflic_global(ixt,knindex(i)) = xtrun_off_lic(ixt,i)
        ENDDO ! DO ixt=1,niso   
      ENDDO
    ENDIF
       
  END SUBROUTINE gestion_neige_besoin_varglob_fonte_neige
#endif


END MODULE fonte_neige_mod
