frun_GR2M.f 7.92 KiB
!------------------------------------------------------------------------------
!    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: 25/11/2019
!------------------------------------------------------------------------------
! 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
7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
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]
141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210
! 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,tanHyp,S1,S2 doubleprecision :: P1,P2,P3,R1,R2,AE,EXCH doubleprecision :: 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 TWS = tanHyp(WS) 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 WS=E/Param(1) IF(WS.GT.13.) WS=13. ! speed-up TWS = tanHyp(WS) 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]
211212213214215216217218219220221222223224
MISC( 3)=St(1) ! Prod ! [numeric] production store level (St(1)) [mm] MISC( 4)=P1 ! Pn ! [numeric] net rainfall (P1) [mm/month] MISC( 5)=AE ! AE ! [numeric] actual evapotranspiration [mm/month] MISC( 6)=P2 ! Perc ! [numeric] percolation (P2) [mm/month] MISC( 7)=P3 ! PR ! [numeric] P3=P1+P2 [mm/month] MISC( 8)=St(2) ! Rout ! [numeric] routing store level (St(2)) [mm] MISC( 9)=EXCH ! EXCH ! [numeric] groundwater exchange (EXCH) [mm/month] MISC(10)=Q ! Qsim ! [numeric] simulated outflow at catchment outlet [mm/month] ENDSUBROUTINE