!------------------------------------------------------------------------------ ! Subroutines relative to the annual GR2M model !------------------------------------------------------------------------------ ! TITLE : airGR ! PROJECT : airGR ! FILE : frun_GR2M.f !------------------------------------------------------------------------------ ! AUTHORS ! Original code: S. Mouelhi ! Cleaning and formatting for airGR: L. Coron ! Further cleaning: G. Thirel !------------------------------------------------------------------------------ ! Creation date: 2003 ! Last modified: 16/04/2020 !------------------------------------------------------------------------------ ! REFERENCES ! Mouelhi S. (2003). Vers une chaîne cohérente de modèles pluie-débit ! conceptuels globaux aux pas de temps pluriannuel, annuel, mensuel et ! journalier. PhD thesis (in French), ENGREF, Cemagref Antony, France. ! ! Mouelhi, S., C. Michel, C. Perrin and V. Andréassian (2006). Stepwise ! development of a two-parameter monthly water balance model. Journal of ! Hydrology, 318(1-4), 200-214. doi:10.1016/j.jhydrol.2005.06.014. !------------------------------------------------------------------------------ ! Quick description of public procedures: ! 1. frun_gr2m ! 2. MOD_GR2M !------------------------------------------------------------------------------ SUBROUTINE frun_gr2m(LInputs,InputsPrecip,InputsPE,NParam,Param, & NStates,StateStart,NOutputs,IndOutputs, & Outputs,StateEnd) ! Subroutine that initializes GR2M, get its parameters, performs the call ! to the MOD_GR2M subroutine at each time step, and stores the final states ! Inputs ! LInputs ! Integer, length of input and output series ! InputsPrecip ! Vector of real, input series of total precipitation [mm/month] ! InputsPE ! Vector of real, input series of potential evapotranspiration (PE) [mm/month] ! NParam ! Integer, number of model parameters ! Param ! Vector of real, parameter set ! NStates ! Integer, number of state variables ! StateStart ! Vector of real, state variables used when the model run starts (store levels [mm]) ! NOutputs ! Integer, number of output series ! IndOutputs ! Vector of integer, indices of output series ! Outputs ! Outputs ! Vector of real, output series ! StateEnd ! Vector of real, state variables at the end of the model run (store levels [mm]) !DEC$ ATTRIBUTES DLLEXPORT :: frun_gr2m Implicit None !! dummies ! in integer, intent(in) :: LInputs,NParam,NStates,NOutputs doubleprecision, dimension(LInputs), intent(in) :: InputsPrecip doubleprecision, dimension(LInputs), intent(in) :: InputsPE doubleprecision, dimension(NParam), intent(in) :: Param doubleprecision, dimension(NStates), intent(in) :: StateStart integer, dimension(NOutputs), intent(in) :: IndOutputs ! out doubleprecision, dimension(NStates), intent(out) :: StateEnd doubleprecision, dimension(LInputs,NOutputs), intent(out) :: Outputs !! locals integer :: I,K integer, parameter :: NMISC=30 doubleprecision, dimension(2) :: St doubleprecision, dimension(NMISC) :: MISC doubleprecision :: P,E,Q !-------------------------------------------------------------- ! Initializations !-------------------------------------------------------------- ! initialization of model states to zero St=0. ! initialization of model states using StateStart St(1)=StateStart(1) St(2)=StateStart(2) ! parameter values ! Param(1) : production store capacity [mm] ! Param(2) : groundwater exchange coefficient [-] ! initialization of model outputs Q = -999.999 MISC = -999.999 ! StateEnd = -999.999 ! initialization made in R ! Outputs = -999.999 ! initialization made in R !-------------------------------------------------------------- ! Time loop !-------------------------------------------------------------- DO k=1,LInputs P =InputsPrecip(k) E =InputsPE(k) ! Q = -999.999 ! MISC = -999.999 !model run on one time step CALL MOD_GR2M(St,Param,P,E,Q,MISC) !storage of outputs DO I=1,NOutputs Outputs(k,I)=MISC(IndOutputs(I)) ENDDO ENDDO ! model states at the end of the run StateEnd(1)=St(1) StateEnd(2)=St(2) RETURN ENDSUBROUTINE !################################################################################################################################ !********************************************************************** SUBROUTINE MOD_GR2M(St,Param,P,E,Q,MISC) ! Calculation of streamflow on a single time step (month) with the GR2M model ! Inputs: ! St Vector of real, model states at the beginning of the time step [mm] ! Param Vector of real, model parameters (Param(1) [mm]; Param(2) [-]) ! P Real, value of rainfall during the time step [mm/month] ! E Real, value of potential evapotranspiration during the time step [mm/month] ! Outputs: ! St Vector of real, model states at the end of the time step [mm] ! Q Real, value of simulated flow at the catchment outlet for the time step [mm/month] ! MISC Vector of real, model outputs for the time step [mm/month] !********************************************************************** Implicit None !! locals integer, parameter :: NParam=2,NMISC=30 doubleprecision :: WS,S1,S2 doubleprecision :: P1,P2,P3,R1,R2,AE,EXCH,PS doubleprecision :: expWS, TWS, Sr ! speed-up !! dummies ! in doubleprecision, dimension(NParam), intent(in) :: Param doubleprecision, intent(in) :: P,E ! inout doubleprecision, dimension(2), intent(inout) :: St ! out doubleprecision, intent(out) :: Q doubleprecision, dimension(NMISC), intent(out) :: MISC ! Production store WS=P/Param(1) IF(WS.GT.13.) WS=13. ! speed-up expWS = exp(2.*WS) TWS = (expWS - 1.)/(expWS + 1.) S1=(St(1)+Param(1)*TWS)/(1.+St(1)/Param(1)*TWS) ! S1=(X(1)+Param(1)*tanHyp(WS))/(1.+X(1)/Param(1)*tanHyp(WS)) ! end speed-up P1=P+St(1)-S1 PS = P - P1 WS=E/Param(1) IF(WS.GT.13.) WS=13. ! speed-up expWS = exp(2.*WS) TWS = (expWS - 1.)/(expWS + 1.) S2=S1*(1.-TWS)/(1.+(1.-S1/Param(1))*TWS) ! S2=S1*(1.-tanHyp(WS))/(1.+(1.-S1/Param(1))*tanHyp(WS)) ! end speed-up AE = S1 - S2 ! Percolation ! speed-up Sr = S2/Param(1) Sr = Sr * Sr * Sr + 1. St(1)=S2/Sr**(1./3.) ! X(1)=S2/(1+(S2/Param(1))**3.)**(1./3.) ! end speed-up P2=S2-St(1) P3=P1+P2 ! QR calculation (routing store) R1=St(2)+P3 ! Water exchange R2=Param(2)*R1 EXCH = R2 - R1 ! Total runoff Q=R2*R2/(R2+60.) ! Updating store level St(2)=R2-Q ! Variables storage MISC( 1)=E ! PE ! [numeric] observed potential evapotranspiration [mm/month] MISC( 2)=P ! Precip ! [numeric] observed total precipitation [mm/month] MISC( 3)=St(1) ! Prod ! [numeric] production store level (St(1)) [mm] MISC( 4)=P1 ! Pn ! [numeric] net rainfall (P1) [mm/month] MISC( 5)=PS ! Ps ! [numeric] part of P filling the production store [mm/month] MISC( 6)=AE ! AE ! [numeric] actual evapotranspiration [mm/month] MISC( 7)=P2 ! Perc ! [numeric] percolation (P2) [mm/month] MISC( 8)=P3 ! PR ! [numeric] P3=P1+P2 [mm/month] MISC( 9)=St(2) ! Rout ! [numeric] routing store level (St(2)) [mm] MISC(10)=EXCH ! EXCH ! [numeric] groundwater exchange (EXCH) [mm/month] MISC(11)=Q ! Qsim ! [numeric] simulated outflow at catchment outlet [mm/month] END SUBROUTINE