Skip to content
GitLab
Commit 688c84b0 authored by Delaigue Olivier's avatar Delaigue Olivier
Browse files

v1.2.0.1 NEW: add CemaNeige hysteresis Fortran file

parent 65b10d1a
dev 115-add-pcict-date-format-managment 118-plot-outputsmodel-remove-case-sensitivity-to-the-which-argument 136-decreased-performance-of-calibration-execution-time 60-migrate-transfoparam-functions-and-createcaliboptions-to-s3-methods 96-parallelize-the-grid-screening-step-during-calibration airGRplus cleanAirGRplus master CRAN_v1.7.0 CRAN_v1.6.12 CRAN_v1.6.11 CRAN_v1.6.10.4 CRAN_v1.6.9.27 CRAN_v1.6.9.24 CRAN_v1.6.9.23 CRAN_v1.6.9.21 CRAN_v1.4.3.65 CRAN_v1.4.3.60 CRAN_v1.4.3.52 CRAN_v1.3.2.42 CRAN_v1.3.2.41 CRAN_v1.3.2.40 CRAN_v1.3.2.39 CRAN_v1.3.2.23 CRAN_v1.3.2.17 CRAN_v1.2.13.16 CRAN_v1.2.13.13
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 77 additions and 35 deletions
+77 -35
src/frun_CEMANEIGE.f
+ 77
− 35
View file @ 688c84b0
SUBROUTINE frun_CEMANEIGE( SUBROUTINE frun_CEMANEIGE(
!inputs !inputs
& LInputs , ! [integer] length of input and output series & LInputs , ! [integer] length of input and output series
...@@ -11,6 +13,7 @@ ...@@ -11,6 +13,7 @@
& Param , ! [double] parameter set & Param , ! [double] parameter set
& NStates , ! [integer] number of state variables used for model initialising = 2 & NStates , ! [integer] number of state variables used for model initialising = 2
& StateStart , ! [double] state variables used when the model run starts & 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 & NOutputs , ! [integer] number of output series
& IndOutputs , ! [integer] indices of output series & IndOutputs , ! [integer] indices of output series
!outputs !outputs
...@@ -25,39 +28,53 @@ ...@@ -25,39 +28,53 @@
!### input and output variables !### input and output variables
integer, intent(in) :: LInputs,NParam,NStates,NOutputs integer, intent(in) :: LInputs,NParam,NStates,NOutputs
doubleprecision, intent(in) :: MeanAnSolidPrecip doubleprecision, intent(in) :: MeanAnSolidPrecip
doubleprecision, dimension(LInputs) :: InputsPrecip doubleprecision, intent(in), dimension(LInputs) :: InputsPrecip
doubleprecision, dimension(LInputs) :: InputsFracSolidPrecip doubleprecision, intent(in), dimension(LInputs) ::
doubleprecision, dimension(LInputs) :: InputsTemp & InputsFracSolidPrecip
doubleprecision, dimension(NParam) :: Param doubleprecision, intent(in), dimension(LInputs) :: InputsTemp
doubleprecision, dimension(NStates) :: StateStart doubleprecision, intent(in), dimension(NParam) :: Param
doubleprecision, dimension(NStates) :: StateEnd doubleprecision, intent(in), dimension(NStates) :: StateStart
integer, dimension(NOutputs) :: IndOutputs logical, intent(in) :: IsHyst ! TRUE if linear hysteresis is used, FALSE otherwise
doubleprecision, dimension(LInputs,NOutputs) :: Outputs doubleprecision, intent(out), dimension(NStates) :: StateEnd
integer, intent(in), dimension(NOutputs) :: IndOutputs
doubleprecision, intent(out), dimension(LInputs,NOutputs) ::
& Outputs
!parameters, internal states and variables !parameters, internal states and variables
doubleprecision CTG,Kf doubleprecision CTG,Kf
doubleprecision G,eTG,PliqAndMelt doubleprecision G,eTG,PliqAndMelt,dG,prct
doubleprecision Tmelt,Gthreshold,MinSpeed doubleprecision Tmelt,Gthreshold,MinSpeed,Gacc,Glocalmax
doubleprecision Pliq,Psol,Gratio,PotMelt,Melt doubleprecision Pliq,Psol,Gratio,PotMelt,Melt,Ginit
integer I,K integer I,K
!-------------------------------------------------------------- !--------------------------------------------------------------
!Initialisations !Initialisations
!-------------------------------------------------------------- !--------------------------------------------------------------
!initilisation des constantes !initialisation of constants
Tmelt=0. Tmelt=0.
Gthreshold=0.9*MeanAnSolidPrecip
MinSpeed=0.1 MinSpeed=0.1
!initilisation of model states using StateStart !initialisation of model states using StateStart
G=StateStart(1) G=StateStart(1)
eTG=StateStart(2) eTG=StateStart(2)
Gthreshold=StateStart(3)
Glocalmax=StateStart(4)
Gratio=0.
PliqAndMelt=0. PliqAndMelt=0.
!setting parameter values !setting parameter values
CTG=Param(1) CTG=Param(1)
Kf=Param(2) 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 !initialisation of model outputs
c StateEnd = -999.999 !initialisation made in R c StateEnd = -999.999 !initialisation made in R
...@@ -71,14 +88,16 @@ c Outputs = -999.999 !initialisation made in R ...@@ -71,14 +88,16 @@ c Outputs = -999.999 !initialisation made in R
DO k=1,LInputs DO k=1,LInputs
!SolidPrecip and LiquidPrecip !SolidPrecip and LiquidPrecip
Pliq=(1-InputsFracSolidPrecip(k))*InputsPrecip(k) Pliq=(1.-InputsFracSolidPrecip(k))*InputsPrecip(k)
Psol=InputsFracSolidPrecip(k)*InputsPrecip(k) Psol=InputsFracSolidPrecip(k)*InputsPrecip(k)
!Snow pack volume before melt !Snow pack volume before melt
IF (IsHyst) Ginit=G
G=G+Psol G=G+Psol
!Snow pack thermal state before melt !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. IF(eTG.GT.0.) eTG=0.
!Potential melt !Potential melt
...@@ -89,46 +108,69 @@ c Outputs = -999.999 !initialisation made in R ...@@ -89,46 +108,69 @@ c Outputs = -999.999 !initialisation made in R
PotMelt=0. PotMelt=0.
ENDIF ENDIF
!Gratio
IF(G.LT.Gthreshold) THEN IF (IsHyst) THEN
Gratio=G/Gthreshold 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 ELSE
Gratio=1. IF(G.LT.Gthreshold) THEN
Gratio=G/Gthreshold
ELSE
Gratio=1.
ENDIF
ENDIF ENDIF
!Actual melt !Actual melt
Melt=((1-MinSpeed)*Gratio+MinSpeed)*PotMelt Melt=((1.-MinSpeed)*Gratio+MinSpeed)*PotMelt
!Update of snow pack volume !Update of snow pack volume
G=G-Melt G=G-Melt
IF (IsHyst) THEN
!Update of Gratio dG=G-Ginit !Melt in case of negative dG, accumulation otherwise
IF(G.LT.Gthreshold) THEN
Gratio=G/Gthreshold
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 ELSE
Gratio=1. IF(G.LT.Gthreshold) THEN
Gratio=G/Gthreshold
ELSE
Gratio=1.
ENDIF
ENDIF ENDIF
!Water volume to pass to the hydrological model !Water volume to pass to the hydrological model
PliqAndMelt=Pliq+Melt PliqAndMelt=Pliq+Melt
!Storage of outputs !Storage of outputs
DO I=1,NOutputs DO I=1,NOutputs
IF(IndOutputs(I).EQ.1) Outputs(k,I)=Pliq ! Pliq ! observed liquid precipitation [mm/day] 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.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.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.4) Outputs(k,I)=eTG ! ThermalState ! thermal state [°C]
IF(IndOutputs(I).EQ.5) Outputs(k,I)=Gratio ! Gratio ! Gratio [-] 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.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.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.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.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
ENDDO ENDDO
StateEnd(1)=G StateEnd(1)=G
StateEnd(2)=eTG StateEnd(2)=eTG
StateEnd(3)=Gthreshold
StateEnd(4)=Glocalmax
RETURN RETURN
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment