Reusing old transformation of X5

R/CreateCalibOptions.R

@@ -116,9 +116,9 @@ CreateCalibOptions <- function(FUN_MOD,FUN_CALIB=Calibration_Michel,FUN_TRANSFO=
                              +5.55, -0.46, +3.75, -9.09, -4.69,
                              +6.10, -0.11, +4.43, -8.60, -0.66),ncol=NParam,byrow=TRUE); }
       if("GR6J"%in% ObjectClass){ 
-        ParamT <- matrix( c( +4.41, +0.41, +2.88, -9.10, -0.13, +0.81,
-                             +5.02, +0.61, +3.45, -8.68, +1.95, +2.27,
-                             +5.58, +0.78, +4.18, -8.12, +3.59, +3.56),ncol=NParam,byrow=TRUE); }
+        ParamT <- matrix( c( +3.60, -1.00, +3.30, -9.10, -0.90, +3.00,
+                             +3.90, -0.50, +4.10, -8.70, +0.10, +4.00,
+                             +4.50, +0.50, +5.00, -8.10, +1.10, +5.00),ncol=NParam,byrow=TRUE); }
       if("GR2M"%in% ObjectClass){ 
         ParamT <- matrix( c( +5.03, -7.15,
                              +5.22, -6.74,
@@ -140,9 +140,9 @@ CreateCalibOptions <- function(FUN_MOD,FUN_CALIB=Calibration_Michel,FUN_TRANSFO=
                              +5.55, -0.46, +3.75, -9.09, -4.69, -9.14, +6.90,
                              +6.10, -0.11, +4.43, -8.60, -0.66, +4.10, +7.21),ncol=NParam,byrow=TRUE); }
       if("CemaNeigeGR6J"%in% ObjectClass){ 
-        ParamT <- matrix( c( +4.41, +0.41, +2.88, -9.10, -0.13, +0.81, -9.96, +6.63,
-                             +5.02, +0.61, +3.45, -8.68, +1.95, +2.27, -9.14, +6.90,
-                             +5.58, +0.78, +4.18, -8.12, +3.59, +3.56, +4.10, +7.21),ncol=NParam,byrow=TRUE); }
+        ParamT <- matrix( c( +3.60, -1.00, +3.30, -9.10, -0.90, +3.00, -6.26, +0.55,
+                             +3.90, -0.50, +4.10, -8.70, +0.10, +4.00, -2.13, +0.92,
+                             +4.50, +0.50, +5.00, -8.10, +1.10, +5.00, +4.86, +1.40),ncol=NParam,byrow=TRUE); }
 
       StartParamList    <- NULL;
       StartParamDistrib <- TransfoParam(ParamIn=ParamT,Direction="TR",FUN_TRANSFO=FUN_TRANSFO);

R/TransfoParam_GR6J.R

@@ -11,7 +11,7 @@ TransfoParam_GR6J <- function(ParamIn,Direction){
     ParamOut[,2] <- sinh(ParamIn[,2]);                   ### GR6J X2 (groundwater exchange coefficient 1)
     ParamOut[,3] <- exp(ParamIn[,3]);                    ### GR6J X3 (routing store capacity)
     ParamOut[,4] <- 20+19.5*(ParamIn[,4]-9.99)/19.98;    ### GR6J X4 (unit hydrograph time constant)
-    ParamOut[,5] <- (ParamIn[,5] + 9.99) / 19.98;        ### GR5J X5 (groundwater exchange coefficient 2)
+    ParamOut[,5] <- ParamIn[,5]/5.;                      ### GR5J X5 (groundwater exchange coefficient 2)
     ParamOut[,6] <- exp(ParamIn[,6]);                    ### GR6J X6 (coefficient for emptying exponential store)
   }	
   if(Direction=="RT"){
@@ -20,7 +20,7 @@ TransfoParam_GR6J <- function(ParamIn,Direction){
     ParamOut[,2] <- asinh(ParamIn[,2]);                  ### GR6J X2 (groundwater exchange coefficient 1)
     ParamOut[,3] <- log(ParamIn[,3]);                    ### GR6J X3 (routing store capacity)
     ParamOut[,4] <- 9.99+19.98*(ParamIn[,4]-20)/19.5;    ### GR6J X4 (unit hydrograph time constant)
-    ParamOut[,5] <- ParamIn[,5] * 19.98 - 9.99;          ### GR5J X5 (groundwater exchange coefficient 2)
+    ParamOut[,5] <- ParamIn[,5] * 5.;                    ### GR5J X5 (groundwater exchange coefficient 2)
     ParamOut[,6] <- log(ParamIn[,6]);                    ### GR6J X6 (coefficient for emptying exponential store)
   }	
 

src/frun_GR6J.f

@@ -137,23 +137,23 @@ C       Param  Vector of model parameters [various units]
 C       P1     Value of rainfall during the time step [mm]
 C       E      Value of potential evapotranspiration during the time step [mm]
 C Outputs:
-C       St     Vector of model states in stores at the beginning of the time step [mm]
-C       StUH1  Vector of model states in Unit Hydrograph 1 at the beginning of the time step [mm]
-C       StUH2  Vector of model states in Unit Hydrograph 2 at the beginning of the time step [mm]
-C       Q      Value of simulated flow at the catchment outlet for the time step [mm]
-C       MISC   Vector of model outputs for the time step [mm]
+C       St     Vector of model states in stores at the end of the time step [mm]
+C       StUH1  Vector of model states in Unit Hydrograph 1 at the end of the time step [mm]
+C       StUH2  Vector of model states in Unit Hydrograph 2 at the end of the time step [mm]
+C       Q      Value of simulated flow at the catchment outlet for the time step [mm/day]
+C       MISC   Vector of model outputs for the time step [mm or mm/day]
 C**********************************************************************
       Implicit None
       INTEGER NH,NMISC,NParam
       PARAMETER (NH=20,NMISC=30)
       PARAMETER (NParam=6)
-      DOUBLEPRECISION St(3),StUH1(NH)
-      DOUBLEPRECISION StUH2(2*NH),OrdUH1(NH),OrdUH2(2*NH)
+      DOUBLEPRECISION St(3)
+	  DOUBLEPRECISION StUH1(NH),StUH2(2*NH),OrdUH1(NH),OrdUH2(2*NH)
       DOUBLEPRECISION Param(NParam)
       DOUBLEPRECISION MISC(NMISC)
       DOUBLEPRECISION P1,E,Q
       DOUBLEPRECISION A,B,C,EN,ER,PN,PR,PS,WS,tanHyp,AR
-      DOUBLEPRECISION PERC,PRHU1,PRHU2,EXCH,QR,QD,QR1
+      DOUBLEPRECISION PERC,PRUH1,PRUH2,EXCH,QR,QD,QR1
       DOUBLEPRECISION AE,AEXCH1,AEXCH2
       INTEGER K
 
@@ -215,27 +215,27 @@ C Percolation from production store
       PR=PR+PERC
 
 C Split of effective rainfall into the two routing components
-      PRHU1=PR*B
-      PRHU2=PR*(1.-B)
+      PRUH1=PR*B
+      PRUH2=PR*(1.-B)
 
 C Convolution of unit hydrograph UH1
       DO K=1,MAX(1,MIN(NH-1,INT(Param(4)+1.)))
-      StUH1(K)=StUH1(K+1)+OrdUH1(K)*PRHU1
+      StUH1(K)=StUH1(K+1)+OrdUH1(K)*PRUH1
       ENDDO
-      StUH1(NH)=OrdUH1(NH)*PRHU1
+      StUH1(NH)=OrdUH1(NH)*PRUH1
 
 C Convolution of unit hydrograph UH2
       DO K=1,MAX(1,MIN(2*NH-1,2*INT(Param(4)+1.)))
-      StUH2(K)=StUH2(K+1)+OrdUH2(K)*PRHU2
+      StUH2(K)=StUH2(K+1)+OrdUH2(K)*PRUH2
       ENDDO
-      StUH2(2*NH)=OrdUH2(2*NH)*PRHU2
+      StUH2(2*NH)=OrdUH2(2*NH)*PRUH2
 
 C Potential intercatchment semi-exchange
       EXCH=Param(2)*(St(2)/Param(3)-Param(5))
 
 C Routing store
       AEXCH1=EXCH
-      IF((St(2)+StUH1(1)+EXCH).LT.0) AEXCH1=-St(2)-StUH1(1)
+      IF((St(2)+(1-C)*StUH1(1)+EXCH).LT.0) AEXCH1=-St(2)-(1-C)*StUH1(1)
       St(2)=St(2)+(1-C)*StUH1(1)+EXCH
       IF(St(2).LT.0.)St(2)=0.
 	  
@@ -289,8 +289,8 @@ C Variables storage
       MISC( 7)=StUH1(1)      ! Q9     ! outflow from UH1 (Q9) [mm/day]
       MISC( 8)=StUH2(1)      ! Q1     ! outflow from UH2 (Q1) [mm/day]
       MISC( 9)=St(2)         ! Rout   ! routing store level (St(2)) [mm]
-      MISC(10)=EXCH          ! Exch   ! potential semi-exchange between catchments (EXCH) [mm/day]
-      MISC(11)=AEXCH1+AEXCH2 ! AExch  ! actual total exchange between catchments (AEXCH1+AEXCH2) [mm/day]
+      MISC(10)=EXCH          ! Exch   ! potential third-exchange between catchments (EXCH) [mm/day]
+      MISC(11)=AEXCH1+AEXCH2+EXCH ! AExch  ! actual total exchange between catchments (AEXCH1+AEXCH2+EXCH) [mm/day]
       MISC(12)=QR            ! QR     ! outflow from routing store (QR) [mm/day]
       MISC(13)=QR1           ! QR1    ! outflow from exponential store (QR1) [mm/day]
       MISC(14)=St(3)         ! Exp    ! exponential store level (St(3)) (negative) [mm]