diff --git a/src/frun_CEMANEIGE.f b/src/frun_CEMANEIGE.f
index 3e66f254f3ba47c76e7ec27b469a0d78065cd39f..13dd51578177b31802da64b7e6e3e51557a63ba0 100644
--- a/src/frun_CEMANEIGE.f
+++ b/src/frun_CEMANEIGE.f
@@ -1,5 +1,7 @@
+
+
SUBROUTINE frun_CEMANEIGE(
!inputs
& LInputs , ! [integer] length of input and output series
@@ -11,6 +13,7 @@
& Param , ! [double] parameter set
& NStates , ! [integer] number of state variables used for model initialising = 2
& StateStart , ! [double] state variables used when the model run starts
+ & IsHyst , ! [logical] whether we should use the linear hysteresis or not
& NOutputs , ! [integer] number of output series
& IndOutputs , ! [integer] indices of output series
!outputs
@@ -25,39 +28,53 @@
!### input and output variables
integer, intent(in) :: LInputs,NParam,NStates,NOutputs
doubleprecision, intent(in) :: MeanAnSolidPrecip
- doubleprecision, dimension(LInputs) :: InputsPrecip
- doubleprecision, dimension(LInputs) :: InputsFracSolidPrecip
- doubleprecision, dimension(LInputs) :: InputsTemp
- doubleprecision, dimension(NParam) :: Param
- doubleprecision, dimension(NStates) :: StateStart
- doubleprecision, dimension(NStates) :: StateEnd
- integer, dimension(NOutputs) :: IndOutputs
- doubleprecision, dimension(LInputs,NOutputs) :: Outputs
+ doubleprecision, intent(in), dimension(LInputs) :: InputsPrecip
+ doubleprecision, intent(in), dimension(LInputs) ::
+ & InputsFracSolidPrecip
+ doubleprecision, intent(in), dimension(LInputs) :: InputsTemp
+ doubleprecision, intent(in), dimension(NParam) :: Param
+ doubleprecision, intent(in), dimension(NStates) :: StateStart
+ logical, intent(in) :: IsHyst ! TRUE if linear hysteresis is used, FALSE otherwise
+ doubleprecision, intent(out), dimension(NStates) :: StateEnd
+ integer, intent(in), dimension(NOutputs) :: IndOutputs
+ doubleprecision, intent(out), dimension(LInputs,NOutputs) ::
+ & Outputs
!parameters, internal states and variables
doubleprecision CTG,Kf
- doubleprecision G,eTG,PliqAndMelt
- doubleprecision Tmelt,Gthreshold,MinSpeed
- doubleprecision Pliq,Psol,Gratio,PotMelt,Melt
+ doubleprecision G,eTG,PliqAndMelt,dG,prct
+ doubleprecision Tmelt,Gthreshold,MinSpeed,Gacc,Glocalmax
+ doubleprecision Pliq,Psol,Gratio,PotMelt,Melt,Ginit
integer I,K
!--------------------------------------------------------------
!Initialisations
!--------------------------------------------------------------
- !initilisation des constantes
+ !initialisation of constants
Tmelt=0.
- Gthreshold=0.9*MeanAnSolidPrecip
MinSpeed=0.1
- !initilisation of model states using StateStart
+ !initialisation of model states using StateStart
G=StateStart(1)
eTG=StateStart(2)
+ Gthreshold=StateStart(3)
+ Glocalmax=StateStart(4)
+ Gratio=0.
PliqAndMelt=0.
!setting parameter values
CTG=Param(1)
Kf=Param(2)
+ IF (IsHyst) THEN
+ Gacc=Param(3)
+ prct=Param(4)
+
+ Gthreshold=prct*MeanAnSolidPrecip
+ Glocalmax=Gthreshold
+ ELSE
+ Gthreshold=0.9*MeanAnSolidPrecip
+ ENDIF
!initialisation of model outputs
c StateEnd = -999.999 !initialisation made in R
@@ -71,14 +88,16 @@ c Outputs = -999.999 !initialisation made in R
DO k=1,LInputs
!SolidPrecip and LiquidPrecip
- Pliq=(1-InputsFracSolidPrecip(k))*InputsPrecip(k)
+ Pliq=(1.-InputsFracSolidPrecip(k))*InputsPrecip(k)
Psol=InputsFracSolidPrecip(k)*InputsPrecip(k)
!Snow pack volume before melt
+ IF (IsHyst) Ginit=G
G=G+Psol
+
!Snow pack thermal state before melt
- eTG=CTG*eTG + (1-CTG)*InputsTemp(k)
+ eTG=CTG*eTG + (1.-CTG)*InputsTemp(k)
IF(eTG.GT.0.) eTG=0.
!Potential melt
@@ -89,46 +108,69 @@ c Outputs = -999.999 !initialisation made in R
PotMelt=0.
ENDIF
- !Gratio
- IF(G.LT.Gthreshold) THEN
- Gratio=G/Gthreshold
+
+ IF (IsHyst) THEN
+ IF (Potmelt.GT.0.) THEN
+ IF (G.LT.Glocalmax.AND.Gratio.EQ.1.) Glocalmax=G !Update in case of potential melt and G lower than Gseuil
+ Gratio=MIN(G/Glocalmax,1.)
+ ENDIF
ELSE
- Gratio=1.
+ IF(G.LT.Gthreshold) THEN
+ Gratio=G/Gthreshold
+ ELSE
+ Gratio=1.
+ ENDIF
ENDIF
!Actual melt
- Melt=((1-MinSpeed)*Gratio+MinSpeed)*PotMelt
+ Melt=((1.-MinSpeed)*Gratio+MinSpeed)*PotMelt
!Update of snow pack volume
G=G-Melt
-
- !Update of Gratio
- IF(G.LT.Gthreshold) THEN
- Gratio=G/Gthreshold
+ IF (IsHyst) THEN
+ dG=G-Ginit !Melt in case of negative dG, accumulation otherwise
+
+
+ IF (dG.GT.0.) THEN
+ Gratio = MIN(Gratio+(Psol-Melt)/Gacc,1.) !Psol - Melt = dG
+ IF (Gratio.EQ.1.) Glocalmax = Gthreshold
+ ENDIF
+ IF (dG.LT.0.) THEN
+ Gratio=MIN(G/Glocalmax,1.)
+ ENDIF
ELSE
- Gratio=1.
+ IF(G.LT.Gthreshold) THEN
+ Gratio=G/Gthreshold
+ ELSE
+ Gratio=1.
+ ENDIF
ENDIF
+
!Water volume to pass to the hydrological model
PliqAndMelt=Pliq+Melt
!Storage of outputs
DO I=1,NOutputs
- IF(IndOutputs(I).EQ.1) Outputs(k,I)=Pliq ! Pliq ! observed liquid precipitation [mm/day]
- IF(IndOutputs(I).EQ.2) Outputs(k,I)=Psol ! Psol ! observed solid precipitation [mm/day]
- IF(IndOutputs(I).EQ.3) Outputs(k,I)=G ! SnowPack ! snow pack [mm]
- IF(IndOutputs(I).EQ.4) Outputs(k,I)=eTG ! ThermalState ! thermal state [°C]
- IF(IndOutputs(I).EQ.5) Outputs(k,I)=Gratio ! Gratio ! Gratio [-]
- IF(IndOutputs(I).EQ.6) Outputs(k,I)=PotMelt ! PotMelt ! potential snow melt [mm/day]
- IF(IndOutputs(I).EQ.7) Outputs(k,I)=Melt ! Melt ! melt [mm/day]
- IF(IndOutputs(I).EQ.8) Outputs(k,I)=PliqAndMelt ! PliqAndMelt ! liquid precipitation + melt [mm/day]
- IF(IndOutputs(I).EQ.9) Outputs(k,I)=InputsTemp(k) ! Temp ! air temperature [°C]
+ IF(IndOutputs(I).EQ.1) Outputs(k,I)=Pliq ! Pliq ! observed liquid precipitation [mm/day]
+ IF(IndOutputs(I).EQ.2) Outputs(k,I)=Psol ! Psol ! observed solid precipitation [mm/day]
+ IF(IndOutputs(I).EQ.3) Outputs(k,I)=G ! SnowPack ! snow pack [mm]
+ IF(IndOutputs(I).EQ.4) Outputs(k,I)=eTG ! ThermalState ! thermal state [°C]
+ IF(IndOutputs(I).EQ.5) Outputs(k,I)=Gratio ! Gratio ! Gratio [-]
+ IF(IndOutputs(I).EQ.6) Outputs(k,I)=PotMelt ! PotMelt ! potential snow melt [mm/day]
+ IF(IndOutputs(I).EQ.7) Outputs(k,I)=Melt ! Melt ! melt [mm/day]
+ IF(IndOutputs(I).EQ.8) Outputs(k,I)=PliqAndMelt ! PliqAndMelt ! liquid precipitation + melt [mm/day]
+ IF(IndOutputs(I).EQ.9) Outputs(k,I)=InputsTemp(k) ! Temp ! air temperature [°C]
+ IF(IndOutputs(I).EQ.10) Outputs(k,I)=Gthreshold ! Gthreshold ! melt threshold [mm]
+ IF(IndOutputs(I).EQ.11) Outputs(k,I)=Glocalmax ! Glocalmax ! local melt threshold for hysteresis [mm]
ENDDO
ENDDO
StateEnd(1)=G
StateEnd(2)=eTG
+ StateEnd(3)=Gthreshold
+ StateEnd(4)=Glocalmax
RETURN