diff --git a/DESCRIPTION b/DESCRIPTION
index 119690e713f663dbbce86c12b83bc02c66f53337..8b40fa90aa3ac4f463dee0a4998f80071ebe7fe1 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,8 +1,8 @@
 Package: airGR
 Type: Package
 Title: Suite of GR Hydrological Models for Precipitation-Runoff Modelling
-Version: 1.3.2.63
-Date: 2019-11-20
+Version: 1.3.2.64
+Date: 2019-11-21
 Authors@R: c(
   person("Laurent", "Coron", role = c("aut", "trl"), comment = c(ORCID = "0000-0002-1503-6204")),
   person("Olivier", "Delaigue", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-7668-8468"), email = "airGR@irstea.fr"),
diff --git a/NEWS.md b/NEWS.md
index c365ef72eae9d5c8c97f5f6789be0249f848471a..4f4bc8b15c4b428877ea97c231624c2ed7f7ff1e 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -2,7 +2,7 @@
 
 
 
-### 1.3.2.63 Release Notes (2019-11-20)
+### 1.3.2.64 Release Notes (2019-11-21)
 
 
 #### New features
diff --git a/src/frun_CEMANEIGE.f b/src/frun_CEMANEIGE.f
index c4ca0aeaefb272a33635313f20d783d046c19ff5..1ff16a83b29ff01aaf956d647e7d42f6e0af9288 100644
--- a/src/frun_CEMANEIGE.f
+++ b/src/frun_CEMANEIGE.f
@@ -1,31 +1,71 @@
+!------------------------------------------------------------------------------
+!    Subroutines relative to the annual GR4J model
+!------------------------------------------------------------------------------
+! TITLE   : airGR
+! PROJECT : airGR
+! FILE    : frun_CEMANEIGE.f
+!------------------------------------------------------------------------------
+! AUTHORS
+! Original code: A. Valéry, P. Riboust
+! Cleaning and formatting for airGR: L. Coron
+! Further cleaning: G. Thirel
+!------------------------------------------------------------------------------
+! Creation date: 2011
+! Last modified: 22/11/2019
+!------------------------------------------------------------------------------
+! REFERENCES
+! Riboust, P., G. Thirel, N. Le Moine and P. Ribstein (2019), Revisiting a 
+! simple degree-day model for integrating satellite data: implementation of 
+! SWE-SCA hystereses. Journal of Hydrology and Hydromechanics, 
+! doi:10.2478/johh-2018-0004, 67, 1, 70–81. 
+!
+! Valéry, A., V. Andréassian and C. Perrin (2014), "As simple as possible but 
+! not simpler": what is useful in a temperature-based snow-accounting routine? 
+! Part 1 - Comparison of six snow accounting routines on 380 catchments, 
+! Journal of Hydrology, doi:10.1016/j.jhydrol.2014.04.059. 
+!
+! Valéry, A., V. Andréassian and C. Perrin (2014), "As simple as possible but 
+! not simpler": What is useful in a temperature-based snow-accounting routine? 
+! Part 2 - Sensitivity analysis of the Cemaneige snow accounting routine on 
+! 380 catchments, Journal of Hydrology, doi:10.1016/j.jhydrol.2014.04.058.
+!------------------------------------------------------------------------------
+! Quick description of public procedures:
+!         1. frun_cemaneige
+!------------------------------------------------------------------------------
+
+
+      SUBROUTINE frun_cemaneige(LInputs,InputsPrecip,
+     &                          InputsFracSolidPrecip,InputsTemp,
+     &                          MeanAnSolidPrecip,NParam,Param,NStates,
+     &                          StateStart,IsHyst,NOutputs,IndOutputs,
+     &                          Outputs,StateEnd)
+! Subroutine that runs the CemaNeige model 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/time step]
+!       InputsFracSolidPrecip ! Vector of real, input series of fraction of solid precipitation [0-1]
+!       InputsTemp   ! Vector of real, input series of air mean temperature [degC]
+!       MeanAnSolidPrecip     ! Real, value of annual mean solid precip [mm/y]
+!       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] and [-] and thresholds [mm])
+!       IsHyst       ! integer, whether we should use the linear hysteresis [1] or not [0]
+!       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] and [-] and thresholds [mm])
 
 
 
-
-      SUBROUTINE frun_cemaneige(
-                                 ! inputs
-     &                             LInputs              , ! [integer] length of input and output series
-     &                             InputsPrecip         , ! [double]  input series of total precipitation [mm/time step]
-     &                             InputsFracSolidPrecip, ! [double]  input series of fraction of solid precipitation [0-1]
-     &                             InputsTemp           , ! [double]  input series of air mean temperature [degC]
-     &                             MeanAnSolidPrecip    , ! [double]  value of annual mean solid precip [mm/y]
-     &                             NParam               , ! [integer] number of model parameter
-     &                             Param                , ! [double]  parameter set
-     &                             NStates              , ! [integer] number of state variables used for model initialization = 4
-     &                             StateStart           , ! [double]  state variables used when the model run starts
-     &                             IsHyst               , ! [integer] whether we should use the linear hysteresis or not
-     &                             NOutputs             , ! [integer] number of output series
-     &                             IndOutputs           , ! [integer] indices of output series
-                                 ! outputs
-     &                             Outputs              , ! [double]  output series
-     &                             StateEnd             ) ! [double]  state variables at the end of the model run
-
-
       !DEC$ ATTRIBUTES DLLEXPORT :: frun_cemaneige
 
 
       Implicit None
-      ! input and output variables
+
+      !! dummies
+      ! in
       integer, intent(in) :: LInputs,NParam,NStates,NOutputs
       doubleprecision, intent(in) :: MeanAnSolidPrecip
       doubleprecision, intent(in), dimension(LInputs) :: InputsPrecip
@@ -34,19 +74,20 @@
       doubleprecision, intent(in), dimension(LInputs) :: InputsTemp
       doubleprecision, intent(in), dimension(NParam)  :: Param
       doubleprecision, intent(in), dimension(NStates) :: StateStart
-      integer, intent(in) :: IsHyst                       ! 1 if linear hysteresis is used, 0 otherwise
-      logical :: IsHystBool                               ! TRUE if linear hysteresis is used, FALSE otherwise
-      doubleprecision, intent(out), dimension(NStates) :: StateEnd
+      integer, intent(in) :: IsHyst ! 1 if linear hysteresis is used, 0 otherwise
       integer, intent(in), dimension(NOutputs) :: IndOutputs
+      ! out
+      doubleprecision, intent(out), dimension(NStates) :: StateEnd
       doubleprecision, intent(out), dimension(LInputs,NOutputs) :: 
      & Outputs
 
-      ! parameters, internal states and variables
-      doubleprecision CTG,Kf
-      doubleprecision G,eTG,PliqAndMelt,dG,prct
-      doubleprecision Tmelt,Gthreshold,MinSpeed,Gacc,Glocalmax
-      doubleprecision Pliq,Psol,Gratio,PotMelt,Melt,Ginit
-      integer I,K
+      !! locals
+      logical :: IsHystBool         ! TRUE if linear hysteresis is used, FALSE otherwise
+      doubleprecision :: CTG,Kf
+      doubleprecision :: G,eTG,PliqAndMelt,dG,prct
+      doubleprecision :: Tmelt,Gthreshold,MinSpeed,Gacc,Glocalmax
+      doubleprecision :: Pliq,Psol,Gratio,PotMelt,Melt,Ginit
+      integer :: I,K
 
       IF (IsHyst .EQ. 1) IsHystBool = .TRUE.
       IF (IsHyst .EQ. 0) IsHystBool = .FALSE.
diff --git a/src/frun_GR1A.f b/src/frun_GR1A.f
index 00df0efe7400e8132b73eeab261a7c814a418475..2789e0b9cd0304ea267842c1fac915dac4b88af7 100644
--- a/src/frun_GR1A.f
+++ b/src/frun_GR1A.f
@@ -1,41 +1,72 @@
-
-
-      SUBROUTINE frun_gr1a(
-                                 ! inputs
-     &                             LInputs      , ! [integer] length of input and output series
-     &                             InputsPrecip , ! [double] input series of total precipitation [mm/year]
-     &                             InputsPE     , ! [double] input series of potential evapotranspiration (PE) [mm/year]
-     &                             NParam       , ! [integer] number of model parameters
-     &                             Param        , ! [double] parameter set
-     &                             NStates      , ! [integer] number of state variables
-     &                             StateStart   , ! [double] state variables used when the model run starts (none here)
-     &                             NOutputs     , ! [integer] number of output series
-     &                             IndOutputs   , ! [integer] indices of output series
-                                 ! outputs
-     &                             Outputs      , ! [double] output series
-     &                             StateEnd     ) ! [double] state variables at the end of the model run (none here)
+!------------------------------------------------------------------------------
+!    Subroutines relative to the annual GR1A model
+!------------------------------------------------------------------------------
+! TITLE   : airGR
+! PROJECT : airGR
+! FILE    : frun_GR1A.f
+!------------------------------------------------------------------------------
+! AUTHORS
+! Original code: S. Mouelhi
+! Cleaning and formatting for airGR: L. Coron
+! Further cleaning: G. Thirel
+!------------------------------------------------------------------------------
+! Creation date: 2003
+! Last modified: 21/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.
+!------------------------------------------------------------------------------
+! Quick description of public procedures:
+!         1. frun_gr1a
+!         2. MOD_GR1A
+!------------------------------------------------------------------------------
+
+
+      SUBROUTINE frun_gr1a(LInputs,InputsPrecip,InputsPE,NParam,Param,
+     &                     NStates,StateStart,NOutputs,IndOutputs,
+     &                     Outputs,StateEnd)
+! Subroutine that initializes GR1A, get its parameters, performs the call 
+! to the MOD_GR1A 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/year]
+!       InputsPE     ! Vector of real, input series of potential evapotranspiration (PE) [mm/year]
+!       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 (none here)
+!       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 (none here)
 
 
       !DEC$ ATTRIBUTES DLLEXPORT :: frun_gr1a
 
 
       Implicit None
-      ! input and output variables
+
+      !! dummies
+      ! in
       integer, intent(in) :: LInputs,NParam,NStates,NOutputs
-      doubleprecision, dimension(LInputs)  :: InputsPrecip
-      doubleprecision, dimension(LInputs)  :: InputsPE
-      doubleprecision, dimension(NParam)   :: Param
-      doubleprecision, dimension(NStates)  :: StateStart
-      doubleprecision, dimension(NStates)  :: StateEnd
-      integer, dimension(NOutputs) :: IndOutputs
-      doubleprecision, dimension(LInputs,NOutputs) :: Outputs
-
-      ! parameters, internal states and variables
-      integer NMISC
-      parameter (NMISC=3)
-      doubleprecision MISC(NMISC)
+      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=3
+      doubleprecision, dimension(NMISC) :: MISC
       doubleprecision P0,P1,E1,Q
-      integer I,K
 
       !--------------------------------------------------------------
       ! Initializations
@@ -49,8 +80,8 @@
       MISC = -999.999
 !      Outputs = -999.999  ! initialization made in R
 
-      StateStart = -999.999  ! CRAN-compatibility updates
-      StateEnd = -999.999  ! CRAN-compatibility updates
+!      StateStart = -999.999  ! CRAN-compatibility updates
+!      StateEnd = -999.999  ! CRAN-compatibility updates
 
 
       !--------------------------------------------------------------
@@ -85,25 +116,31 @@
 
 !**********************************************************************
       SUBROUTINE MOD_GR1A(Param,P0,P1,E1,Q,MISC)
-! Calculation of streamflow on a single time step with the GR1A model
-! Inputs:
-!       Param  Vector of model parameters (Param(1) [-])
-!       P0     Value of rainfall during the previous time step [mm/year]
-!       P1     Value of rainfall during the current time step [mm/year]
-!       E1     Value of potential evapotranspiration during the current time step [mm/year]
-! Outputs:
-!       Q      Value of simulated flow at the catchment outlet for the current time step [mm/year]
-!       MISC   Vector of model outputs for the time step [mm/year]
+! Calculation of streamflow on a single time step (year) with the GR1A model
+! Inputs
+!       Param ! Vector of real, model parameters (Param(1) [-])
+!       P0    ! Real, value of rainfall during the previous time step [mm/year]
+!       P1    ! Real, value of rainfall during the current time step [mm/year]
+!       E1    ! Real, value of potential evapotranspiration during the current time step [mm/year]
+! Outputs
+!       Q     ! Real, value of simulated flow at the catchment outlet for the current time step [mm/year]
+!       MISC  ! Vector of real, model outputs for the time step [mm/year]
 !**********************************************************************
       Implicit None
-      INTEGER NMISC,NParam
-      PARAMETER (NMISC=3)
-      PARAMETER (NParam=1)
-      DOUBLEPRECISION Param(NParam)
-      DOUBLEPRECISION MISC(NMISC)
-      DOUBLEPRECISION P0,P1,E1,Q
-
-      DOUBLEPRECISION tt ! speed-up
+
+      !! locals
+      integer, parameter :: NMISC=3
+      integer, parameter :: NParam=1
+      doubleprecision :: tt ! speed-up
+      
+      !! dummies
+      ! in
+      doubleprecision, dimension(NParam), intent(in) :: Param
+      doubleprecision, intent(in) :: P0,P1,E1
+      ! out
+      doubleprecision, dimension(NMISC), intent(out) :: MISC
+      doubleprecision, intent(out) :: Q
+      
 
 ! Runoff
       ! speed-up
diff --git a/src/frun_GR2M.f b/src/frun_GR2M.f
index 904d9befe5003a71a3fa6535d2b4ee83470390c3..d23dc3b9590ee2b8765cd1fddd09132ca283ffc7 100644
--- a/src/frun_GR2M.f
+++ b/src/frun_GR2M.f
@@ -1,42 +1,77 @@
-
-
-      SUBROUTINE frun_gr2m(
-                                 ! inputs
-     &                             LInputs      , ! [integer] length of input and output series
-     &                             InputsPrecip , ! [double]  input series of total precipitation [mm/month]
-     &                             InputsPE     , ! [double]  input series of potential evapotranspiration (PE) [mm/month]
-     &                             NParam       , ! [integer] number of model parameters
-     &                             Param        , ! [double]  parameter set
-     &                             NStates      , ! [integer] number of state variables
-     &                             StateStart   , ! [double]  state variables used when the model run starts (store levels [mm])
-     &                             NOutputs     , ! [integer] number of output series
-     &                             IndOutputs   , ! [integer] indices of output series
-                                 ! outputs
-     &                             Outputs      , ! [double]  output series
-     &                             StateEnd     ) ! [double]  state variables at the end of the model run (store levels [mm])
+!------------------------------------------------------------------------------
+!    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: 21/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
-      ! input and output variables
+
+      !! dummies
+      ! in
       integer, intent(in) :: LInputs,NParam,NStates,NOutputs
-      doubleprecision, dimension(LInputs)  :: InputsPrecip
-      doubleprecision, dimension(LInputs)  :: InputsPE
-      doubleprecision, dimension(NParam)   :: Param
-      doubleprecision, dimension(NStates)  :: StateStart
-      doubleprecision, dimension(NStates)  :: StateEnd
-      integer, dimension(NOutputs) :: IndOutputs
-      doubleprecision, dimension(LInputs,NOutputs) :: Outputs
-
-      ! parameters, internal states and variables
-      integer NMISC
-      parameter (NMISC=30)
-      doubleprecision St(2)
-      doubleprecision MISC(NMISC)
-      doubleprecision P,E,Q
-      integer I,K
+      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
@@ -97,28 +132,33 @@
       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 model states at the beginning of the time step [mm]
-!       Param  Vector of model parameters (Param(1) [mm]; Param(2) [-])
-!       P      Value of rainfall during the time step [mm/month]
-!       E      Value of potential evapotranspiration during the time step [mm/month]
+!       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 model states at the end of the time step [mm]
-!       Q      Value of simulated flow at the catchment outlet for the time step [mm/month]
-!       MISC   Vector of model outputs for the time step [mm/month]
+!       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
-      INTEGER NMISC,NParam
-      PARAMETER (NMISC=30)
-      PARAMETER (NParam=2)
-      DOUBLEPRECISION St(2)
-      DOUBLEPRECISION Param(NParam)
-      DOUBLEPRECISION MISC(NMISC)
-      DOUBLEPRECISION P,E,Q
-      DOUBLEPRECISION WS,tanHyp,S1,S2
-      DOUBLEPRECISION P1,P2,P3,R1,R2,AE,EXCH
-
-      DOUBLEPRECISION TWS, Sr ! speed-up
 
+      !! 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.
diff --git a/src/frun_GR4H.f b/src/frun_GR4H.f
index 5ffe19e4cf3e35a267c70b8b73ec38bb2d77fd59..03cf31cce165877c9217631c0ed4083f50796e17 100644
--- a/src/frun_GR4H.f
+++ b/src/frun_GR4H.f
@@ -1,44 +1,75 @@
-
-
-      SUBROUTINE frun_gr4h(
-                                 ! inputs
-     &                             LInputs      , ! [integer] length of input and output series
-     &                             InputsPrecip , ! [double]  input series of total precipitation [mm/hour]
-     &                             InputsPE     , ! [double]  input series of potential evapotranspiration (PE) [mm/hour]
-     &                             NParam       , ! [integer] number of model parameters
-     &                             Param        , ! [double]  parameter set
-     &                             NStates      , ! [integer] number of state variables
-     &                             StateStart   , ! [double]  state variables used when the model run starts (store levels [mm] and Unit Hydrograph (UH) storages [mm])
-     &                             NOutputs     , ! [integer] number of output series
-     &                             IndOutputs   , ! [integer] indices of output series
-                                 ! outputs
-     &                             Outputs      , ! [double]  output series
-     &                             StateEnd     ) ! [double]  state variables at the end of the model run (store levels [mm] and Unit Hydrograph (UH) storages [mm])
+!------------------------------------------------------------------------------
+!    Subroutines relative to the annual GR4H model
+!------------------------------------------------------------------------------
+! TITLE   : airGR
+! PROJECT : airGR
+! FILE    : frun_GR4H.f
+!------------------------------------------------------------------------------
+! AUTHORS
+! Original code: C. Perrin
+! Cleaning and formatting for airGR: L. Coron
+! Further cleaning: G. Thirel
+!------------------------------------------------------------------------------
+! Creation date: 2003
+! Last modified: 21/11/2019
+!------------------------------------------------------------------------------
+! REFERENCES
+! Perrin, C., C. Michel and V. Andréassian (2003). Improvement of a 
+! parsimonious model for streamflow simulation. Journal of Hydrology, 
+! 279(1-4), 275-289, doi:10.1016/S0022-1694(03)00225-7. 
+!------------------------------------------------------------------------------
+! Quick description of public procedures:
+!         1. frun_gr4h
+!         2. MOD_GR4H
+!------------------------------------------------------------------------------
+
+
+      SUBROUTINE frun_gr4h(LInputs,InputsPrecip,InputsPE,NParam,Param,
+     &                     NStates,StateStart,NOutputs,IndOutputs,
+     &                     Outputs,StateEnd)
+! Subroutine that initializes GR4H, get its parameters, performs the call 
+! to the MOD_GR4H 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/hour]
+!       InputsPE     ! Vector of real, input series of potential evapotranspiration (PE) [mm/hour]
+!       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] and Unit Hydrograph (UH) storages [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] and Unit Hydrograph (UH) storages [mm])
 
 
       !DEC$ ATTRIBUTES DLLEXPORT :: frun_gr4h
 
 
       Implicit None
-      ! input and output variables
+
+      !! dummies
+      ! in
       integer, intent(in) :: LInputs,NParam,NStates,NOutputs
-      doubleprecision, dimension(LInputs)  :: InputsPrecip
-      doubleprecision, dimension(LInputs)  :: InputsPE
-      doubleprecision, dimension(NParam)   :: Param
-      doubleprecision, dimension(NStates)  :: StateStart
-      doubleprecision, dimension(NStates)  :: StateEnd
-      integer, dimension(NOutputs) :: IndOutputs
-      doubleprecision, dimension(LInputs,NOutputs) :: Outputs
-
-      ! parameters, internal states and variables
-      integer NH,NMISC
-      parameter (NH=480,NMISC=30)
-      doubleprecision St(2), StUH1(NH), StUH2(2*NH)
-      doubleprecision OrdUH1(NH), OrdUH2(2*NH)
-      doubleprecision MISC(NMISC)
-      doubleprecision D
-      doubleprecision P1,E,Q
-      integer I,K
+      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 :: NH=480,NMISC=30
+      doubleprecision, dimension(2)    :: St
+      doubleprecision, dimension(NH)   :: StUH1, OrdUH1
+      doubleprecision, dimension(2*NH) :: StUH2, OrdUH2
+      doubleprecision, dimension(NMISC) :: MISC
+      doubleprecision :: D,P1,E,Q
 
       !--------------------------------------------------------------
       ! Initializations
@@ -122,41 +153,48 @@
 !**********************************************************************
       SUBROUTINE MOD_GR4H(St,StUH1,StUH2,OrdUH1,OrdUH2,Param,P1,E,Q,
      &MISC)
-! Run on a single time step with the GR4H model
+! Calculation of streamflow on a single time step (hour) with the GR4H model
 ! Inputs:
-!       St     Vector of model states in stores at the beginning of the time step [mm]
-!       StUH1  Vector of model states in Unit Hydrograph 1 at the beginning of the time step [mm]
-!       StUH2  Vector of model states in Unit Hydrograph 2 at the beginning of the time step [mm]
-!       OrdUH1 Vector of ordinates in UH1 [-]
-!       OrdUH2 Vector of ordinates in UH2 [-]
-!       Param  Vector of model parameters [various units]
-!       P1     Value of rainfall during the time step [mm]
-!       E      Value of potential evapotranspiration during the time step [mm]
+!       St     Vector of real, model states in stores at the beginning of the time step [mm]
+!       StUH1  Vector of real, model states in Unit Hydrograph 1 at the beginning of the time step [mm]
+!       StUH2  Vector of real, model states in Unit Hydrograph 2 at the beginning of the time step [mm]
+!       OrdUH1 Vector of real, ordinates in UH1 [-]
+!       OrdUH2 Vector of real, ordinates in UH2 [-]
+!       Param  Vector of real, model parameters [various units]
+!       P1     Real, value of rainfall during the time step [mm]
+!       E      Real, value of potential evapotranspiration during the time step [mm]
 ! Outputs:
-!       St     Vector of model states in stores at the beginning of the time step [mm]
-!       StUH1  Vector of model states in Unit Hydrograph 1 at the beginning of the time step [mm]
-!       StUH2  Vector of model states in Unit Hydrograph 2 at the beginning of the time step [mm]
-!       Q      Value of simulated flow at the catchment outlet for the time step [mm/hour]
-!       MISC   Vector of model outputs for the time step [mm/hour]
+!       St     Vector of real, model states in stores at the beginning of the time step [mm]
+!       StUH1  Vector of real, model states in Unit Hydrograph 1 at the beginning of the time step [mm]
+!       StUH2  Vector of real, model states in Unit Hydrograph 2 at the beginning of the time step [mm]
+!       Q      Real, value of simulated flow at the catchment outlet for the time step [mm/hour]
+!       MISC   Vector of real, model outputs for the time step [mm/hour]
 !**********************************************************************
       Implicit None
-      INTEGER NH,NMISC,NParam
-      PARAMETER (NH=480,NMISC=30)
-      PARAMETER (NParam=4)
-      DOUBLEPRECISION St(2),StUH1(NH)
-      DOUBLEPRECISION StUH2(2*NH),OrdUH1(NH),OrdUH2(2*NH)
-      DOUBLEPRECISION Param(NParam)
-      DOUBLEPRECISION MISC(NMISC)
-      DOUBLEPRECISION P1,E,Q
-      DOUBLEPRECISION A,B,EN,ER,PN,PR,PS,WS,tanHyp
-      DOUBLEPRECISION PERC,PRHU1,PRHU2,EXCH,QR,QD
-      DOUBLEPRECISION AE,AEXCH1,AEXCH2
-      INTEGER K
-
-      DATA B/0.9/
-
-      DOUBLEPRECISION TWS, Sr, Rr ! speed-up
 
+      !! locals
+      integer, parameter :: NParam=4,NMISC=30,NH=480
+      doubleprecision :: A,B,EN,ER,PN,PR,PS,WS,tanHyp
+      doubleprecision :: PERC,PRHU1,PRHU2,EXCH,QR,QD
+      doubleprecision :: AE,AEXCH1,AEXCH2
+      integer :: K
+      data B/0.9/
+      doubleprecision :: TWS, Sr, Rr ! speed-up
+
+      !! dummies
+      ! in
+      doubleprecision, dimension(NParam), intent(in)  :: Param
+      doubleprecision, intent(in) :: P1,E
+      doubleprecision, dimension(NH), intent(inout)   :: OrdUH1
+      doubleprecision, dimension(2*NH), intent(inout) :: OrdUH2
+      ! inout
+      doubleprecision, dimension(2), intent(inout)    :: St
+      doubleprecision, dimension(NH), intent(inout)   :: StUH1
+      doubleprecision, dimension(2*NH), intent(inout) :: StUH2
+      ! out
+      doubleprecision, intent(out) :: Q
+      doubleprecision, dimension(NMISC), intent(out)  :: MISC
+      
       A=Param(1)
 
 
diff --git a/src/frun_GR4J.f b/src/frun_GR4J.f
index 14fb0c1d41405e61880b8f6f38329dc015e32ff3..c3a88d4da3b906600269d5b52c1a34d28d68b08b 100644
--- a/src/frun_GR4J.f
+++ b/src/frun_GR4J.f
@@ -1,44 +1,75 @@
-
-
-      SUBROUTINE frun_gr4j(
-                                 ! inputs
-     &                             LInputs      , ! [integer] length of input and output series
-     &                             InputsPrecip , ! [double]  input series of total precipitation [mm/day]
-     &                             InputsPE     , ! [double]  input series of potential evapotranspiration (PE) [mm/day]
-     &                             NParam       , ! [integer] number of model parameters
-     &                             Param        , ! [double]  parameter set
-     &                             NStates      , ! [integer] number of state variables
-     &                             StateStart   , ! [double]  state variables used when the model run starts (store levels [mm] and Unit Hydrograph (UH) storages [mm])
-     &                             NOutputs     , ! [integer] number of output series
-     &                             IndOutputs   , ! [integer] indices of output series
-                                 ! outputs
-     &                             Outputs      , ! [double]  output series
-     &                             StateEnd     ) ! [double]  state variables at the end of the model run (store levels [mm] and Unit Hydrograph (UH) storages [mm])
+!------------------------------------------------------------------------------
+!    Subroutines relative to the annual GR4J model
+!------------------------------------------------------------------------------
+! TITLE   : airGR
+! PROJECT : airGR
+! FILE    : frun_GR4J.f
+!------------------------------------------------------------------------------
+! AUTHORS
+! Original code: C. Perrin
+! Cleaning and formatting for airGR: L. Coron
+! Further cleaning: G. Thirel
+!------------------------------------------------------------------------------
+! Creation date: 2000
+! Last modified: 21/11/2019
+!------------------------------------------------------------------------------
+! REFERENCES
+! Perrin, C., C. Michel and V. Andréassian (2003). Improvement of a 
+! parsimonious model for streamflow simulation. Journal of Hydrology, 
+! 279(1-4), 275-289. doi:10.1016/S0022-1694(03)00225-7.
+!------------------------------------------------------------------------------
+! Quick description of public procedures:
+!         1. frun_gr4j
+!         2. MOD_GR4J
+!------------------------------------------------------------------------------
+
+
+      SUBROUTINE frun_gr4j(LInputs,InputsPrecip,InputsPE,NParam,Param,
+     &                     NStates,StateStart,NOutputs,IndOutputs,
+     &                     Outputs,StateEnd)
+! Subroutine that initializes GR4J, get its parameters, performs the call 
+! to the MOD_GR4J 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/day]
+!       InputsPE     ! Vector of real, input series of potential evapotranspiration (PE) [mm/day]
+!       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] and Unit Hydrograph (UH) storages [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] and Unit Hydrograph (UH) storages [mm])
 
 
       !DEC$ ATTRIBUTES DLLEXPORT :: frun_gr4j
 
 
       Implicit None
-      ! input and output variables
+
+      !! dummies
+      ! in
       integer, intent(in) :: LInputs,NParam,NStates,NOutputs
-      doubleprecision, dimension(LInputs)  :: InputsPrecip
-      doubleprecision, dimension(LInputs)  :: InputsPE
-      doubleprecision, dimension(NParam)   :: Param
-      doubleprecision, dimension(NStates)  :: StateStart
-      doubleprecision, dimension(NStates)  :: StateEnd
-      integer, dimension(NOutputs) :: IndOutputs
-      doubleprecision, dimension(LInputs,NOutputs) :: Outputs
-
-      ! parameters, internal states and variables
-      integer NH,NMISC
-      parameter (NH=20,NMISC=30)
-      doubleprecision St(2), StUH1(NH), StUH2(2*NH)
-      doubleprecision OrdUH1(NH), OrdUH2(2*NH)
-      doubleprecision MISC(NMISC)
-      doubleprecision D
-      doubleprecision P1,E,Q
-      integer I,K
+      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 :: NH=20,NMISC=30
+      doubleprecision, dimension(2)     :: St
+      doubleprecision, dimension(NH)    :: StUH1, OrdUH1
+      doubleprecision, dimension(2*NH)  :: StUH2, OrdUH2
+      doubleprecision, dimension(NMISC) :: MISC
+      doubleprecision :: D,P1,E,Q
 
       !--------------------------------------------------------------
       ! Initializations
@@ -122,39 +153,47 @@
 !**********************************************************************
       SUBROUTINE MOD_GR4J(St,StUH1,StUH2,OrdUH1,OrdUH2,Param,P1,E,Q,
      &MISC)
-! Run on a single time step with the GR4J model
+! Calculation of streamflow on a single time step (day) with the GR4J model
 ! Inputs:
-!       St     Vector of model states in stores at the beginning of the time step [mm]
-!       StUH1  Vector of model states in Unit Hydrograph 1 at the beginning of the time step [mm]
-!       StUH2  Vector of model states in Unit Hydrograph 2 at the beginning of the time step [mm]
-!       OrdUH1 Vector of ordinates in UH1 [-]
-!       OrdUH2 Vector of ordinates in UH2 [-]
-!       Param  Vector of model parameters [various units]
-!       P1     Value of rainfall during the time step [mm/day]
-!       E      Value of potential evapotranspiration during the time step [mm/day]
+!       St     Vector of real, model states in stores at the beginning of the time step [mm]
+!       StUH1  Vector of real, model states in Unit Hydrograph 1 at the beginning of the time step [mm]
+!       StUH2  Vector of real, model states in Unit Hydrograph 2 at the beginning of the time step [mm]
+!       OrdUH1 Vector of real, ordinates in UH1 [-]
+!       OrdUH2 Vector of real, ordinates in UH2 [-]
+!       Param  Vector of real, model parameters [various units]
+!       P1     Real, value of rainfall during the time step [mm/day]
+!       E      Real, value of potential evapotranspiration during the time step [mm/day]
 ! Outputs:
-!       St     Vector of model states in stores at the end of the time step [mm]
-!       StUH1  Vector of model states in Unit Hydrograph 1 at the end of the time step [mm]
-!       StUH2  Vector of model states in Unit Hydrograph 2 at the end of the time step [mm]
-!       Q      Value of simulated flow at the catchment outlet for the time step [mm/day]
-!       MISC   Vector of model outputs for the time step [mm/day]
+!       St     Vector of real, model states in stores at the end of the time step [mm]
+!       StUH1  Vector of real, model states in Unit Hydrograph 1 at the end of the time step [mm]
+!       StUH2  Vector of real, model states in Unit Hydrograph 2 at the end of the time step [mm]
+!       Q      Real, value of simulated flow at the catchment outlet for the time step [mm/day]
+!       MISC   Vector of real, model outputs for the time step [mm/day]
 !**********************************************************************
       Implicit None
-      INTEGER NH,NMISC,NParam
-      PARAMETER (NH=20,NMISC=30)
-      PARAMETER (NParam=4)
-      DOUBLEPRECISION St(2),StUH1(NH),StUH2(2*NH)
-      DOUBLEPRECISION OrdUH1(NH),OrdUH2(2*NH)
-      DOUBLEPRECISION Param(NParam)
-      DOUBLEPRECISION MISC(NMISC)
-      DOUBLEPRECISION P1,E,Q
-      DOUBLEPRECISION A,B,EN,ER,PN,PR,PS,WS,tanHyp
-      DOUBLEPRECISION PERC,PRHU1,PRHU2,EXCH,QR,QD
-      DOUBLEPRECISION AE,AEXCH1,AEXCH2
-      INTEGER K
-
-      DATA B/0.9/
-      DOUBLEPRECISION TWS, Sr, Rr ! speed-up
+
+      !! locals
+      integer, parameter :: NParam=4,NMISC=30,NH=20
+      doubleprecision :: A,B,EN,ER,PN,PR,PS,WS,tanHyp
+      doubleprecision :: PERC,PRHU1,PRHU2,EXCH,QR,QD
+      doubleprecision :: AE,AEXCH1,AEXCH2
+      integer :: K
+      data B/0.9/
+      doubleprecision :: TWS, Sr, Rr ! speed-up
+
+      !! dummies
+      ! in
+      doubleprecision, dimension(NParam), intent(in)  :: Param
+      doubleprecision, intent(in) :: P1,E
+      doubleprecision, dimension(NH), intent(inout)   :: OrdUH1
+      doubleprecision, dimension(2*NH), intent(inout) :: OrdUH2
+      ! inout
+      doubleprecision, dimension(2), intent(inout)    :: St
+      doubleprecision, dimension(NH), intent(inout)   :: StUH1
+      doubleprecision, dimension(2*NH), intent(inout) :: StUH2
+      ! out
+      doubleprecision, intent(out) :: Q
+      doubleprecision, dimension(NMISC), intent(out)  :: MISC
 
       A=Param(1)
 
diff --git a/src/frun_GR5J.f b/src/frun_GR5J.f
index c9be386bafb934b81342e8c5b427a3708c37adc8..a5baab6a2c279a63ef7536efab565ee4d73fc2a8 100644
--- a/src/frun_GR5J.f
+++ b/src/frun_GR5J.f
@@ -1,44 +1,79 @@
-
-
-      SUBROUTINE frun_gr5j(
-                                 ! inputs
-     &                             LInputs      , ! [integer] length of input and output series
-     &                             InputsPrecip , ! [double]  input series of total precipitation [mm/day]
-     &                             InputsPE     , ! [double]  input series of potential evapotranspiration (PE) [mm/day]
-     &                             NParam       , ! [integer] number of model parameters
-     &                             Param        , ! [double]  parameter set
-     &                             NStates      , ! [integer] number of state variables
-     &                             StateStart   , ! [double]  state variables used when the model run starts (store levels [mm] and Unit Hydrograph (UH) storages [mm])
-     &                             NOutputs     , ! [integer] number of output series
-     &                             IndOutputs   , ! [integer] indices of output series
-                                 ! outputs
-     &                             Outputs      , ! [double]  output series
-     &                             StateEnd     ) ! [double]  state variables at the end of the model run (store levels [mm] and Unit Hydrograph (UH) storages [mm])
+!------------------------------------------------------------------------------
+!    Subroutines relative to the annual GR5J model
+!------------------------------------------------------------------------------
+! TITLE   : airGR
+! PROJECT : airGR
+! FILE    : frun_GR5J.f
+!------------------------------------------------------------------------------
+! AUTHORS
+! Original code: N. Le Moine
+! Cleaning and formatting for airGR: L. Coron
+! Further cleaning: G. Thirel
+!------------------------------------------------------------------------------
+! Creation date: 2006
+! Last modified: 21/11/2019
+!------------------------------------------------------------------------------
+! REFERENCES
+! Le Moine, N. (2008). Le bassin versant de surface vu par le souterrain : une 
+! voie d'amélioration des performances et du réalisme des modèles 
+! pluie-débit ? PhD thesis (French), UPMC, Paris, France. 
+!
+! Pushpalatha, R., C. Perrin, N. Le Moine, T. Mathevet, and V. Andréassian 
+! (2011). A downward structural sensitivity analysis of hydrological models to 
+! improve low-flow simulation. Journal of Hydrology, 411(1-2), 66-76. 
+! doi:10.1016/j.jhydrol.2011.09.034.
+!------------------------------------------------------------------------------
+! Quick description of public procedures:
+!         1. frun_gr5j
+!         2. MOD_GR5J
+!------------------------------------------------------------------------------
+
+
+      SUBROUTINE frun_gr5j(LInputs,InputsPrecip,InputsPE,NParam,Param,
+     &                     NStates,StateStart,NOutputs,IndOutputs,
+     &                     Outputs,StateEnd)
+! Subroutine that initializes GR5J, get its parameters, performs the call 
+! to the MOD_GR5J 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/day]
+!       InputsPE     ! Vector of real, input series of potential evapotranspiration (PE) [mm/day]
+!       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] and Unit Hydrograph (UH) storages [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] and Unit Hydrograph (UH) storages [mm])
 
 
       !DEC$ ATTRIBUTES DLLEXPORT :: frun_gr5j
 
 
       Implicit None
-      ! input and output variables
+
+      !! dummies
+      ! in
       integer, intent(in) :: LInputs,NParam,NStates,NOutputs
-      doubleprecision, dimension(LInputs)  :: InputsPrecip
-      doubleprecision, dimension(LInputs)  :: InputsPE
-      doubleprecision, dimension(NParam)   :: Param
-      doubleprecision, dimension(NStates)  :: StateStart
-      doubleprecision, dimension(NStates)  :: StateEnd
-      integer, dimension(NOutputs) :: IndOutputs
-      doubleprecision, dimension(LInputs,NOutputs) :: Outputs
-
-      ! parameters, internal states and variables
-      integer NH,NMISC
-      parameter (NH=20,NMISC=30)
-      doubleprecision St(2), StUH2(2*NH)
-      doubleprecision OrdUH2(2*NH)
-      doubleprecision MISC(NMISC)
-      doubleprecision D
-      doubleprecision P1,E,Q
-      integer I,K
+      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 :: NH=20,NMISC=30
+      doubleprecision, dimension(2)     :: St
+      doubleprecision, dimension(2*NH)  :: StUH2, OrdUH2
+      doubleprecision, dimension(NMISC) :: MISC
+      doubleprecision :: D,P1,E,Q
 
       !--------------------------------------------------------------
       ! Initializations
@@ -115,36 +150,42 @@
 !**********************************************************************
       SUBROUTINE MOD_GR5J(St,StUH2,OrdUH2,Param,P1,E,Q,
      &MISC)
-! Run on a single time step with the GR5J model
+! Calculation of streamflow on a single time step (day) with the GR5J model
 ! Inputs:
-!       St     Vector of model states in stores at the beginning of the time step [mm]
-!       StUH2  Vector of model states in Unit Hydrograph 2 at the beginning of the time step [mm]
-!       OrdUH2 Vector of ordinates in UH2 [-]
-!       Param  Vector of model parameters [various units]
-!       P1     Value of rainfall during the time step [mm]
-!       E      Value of potential evapotranspiration during the time step [mm]
+!       St     Vector of real, model states in stores at the beginning of the time step [mm]
+!       StUH2  Vector of real, model states in Unit Hydrograph 2 at the beginning of the time step [mm]
+!       OrdUH2 Vector of real, ordinates in UH2 [-]
+!       Param  Vector of real, model parameters [various units]
+!       P1     Real, value of rainfall during the time step [mm]
+!       E      Real, value of potential evapotranspiration during the time step [mm]
 ! Outputs:
-!       St     Vector of model states in stores at the end of the time step [mm]
-!       StUH2  Vector of model states in Unit Hydrograph 2 at the end of the time step [mm]
-!       Q      Value of simulated flow at the catchment outlet for the time step [mm/day]
-!       MISC   Vector of model outputs for the time step [mm or mm/day]
+!       St     Vector of real, model states in stores at the end of the time step [mm]
+!       StUH2  Vector of real, model states in Unit Hydrograph 2 at the end of the time step [mm]
+!       Q      Real, value of simulated flow at the catchment outlet for the time step [mm/day]
+!       MISC   Vector of real, model outputs for the time step [mm or mm/day]
 !**********************************************************************
       Implicit None
-      INTEGER NH,NMISC,NParam
-      PARAMETER (NH=20,NMISC=30)
-      PARAMETER (NParam=5)
-      DOUBLEPRECISION St(2),StUH2(2*NH)
-      DOUBLEPRECISION OrdUH2(2*NH)
-      DOUBLEPRECISION Param(NParam)
-      DOUBLEPRECISION MISC(NMISC)
-      DOUBLEPRECISION P1,E,Q
-      DOUBLEPRECISION A,B,EN,ER,PN,PR,PS,WS,tanHyp
-      DOUBLEPRECISION PERC,EXCH,QR,QD,Q1,Q9
-      DOUBLEPRECISION AE,AEXCH1,AEXCH2
-      INTEGER K
-
-      DATA B/0.9/
-      DOUBLEPRECISION TWS, Sr, Rr ! speed-up
+
+      !! locals
+      integer, parameter :: NParam=5,NMISC=30,NH=20
+      doubleprecision :: A,B,EN,ER,PN,PR,PS,WS,tanHyp
+      doubleprecision :: PERC,EXCH,QR,QD,Q1,Q9
+      doubleprecision :: AE,AEXCH1,AEXCH2
+      integer :: K
+      data B/0.9/
+      doubleprecision :: TWS, Sr, Rr ! speed-up
+
+      !! dummies
+      ! in
+      doubleprecision, dimension(NParam), intent(in)  :: Param
+      doubleprecision, intent(in) :: P1,E
+      doubleprecision, dimension(2*NH), intent(inout) :: OrdUH2
+      ! inout
+      doubleprecision, dimension(2), intent(inout)    :: St
+      doubleprecision, dimension(2*NH), intent(inout) :: StUH2
+      ! out
+      doubleprecision, intent(out) :: Q
+      doubleprecision, dimension(NMISC), intent(out)  :: MISC
 
       A=Param(1)
 
diff --git a/src/frun_GR6J.f b/src/frun_GR6J.f
index 17927fd4dd436a89c847d1dc55311b7a0f40d935..d2c7d2cdf8b8aa94d2b4fcb6e0e9071212719dd6 100644
--- a/src/frun_GR6J.f
+++ b/src/frun_GR6J.f
@@ -1,44 +1,76 @@
-
-
-      SUBROUTINE frun_gr6j(
-                                 ! inputs
-     &                             LInputs      , ! [integer] length of input and output series
-     &                             InputsPrecip , ! [double]  input series of total precipitation [mm/day]
-     &                             InputsPE     , ! [double]  input series of potential evapotranspiration (PE) [mm/day]
-     &                             NParam       , ! [integer] number of model parameters
-     &                             Param        , ! [double]  parameter set
-     &                             NStates      , ! [integer] number of state variables
-     &                             StateStart   , ! [double]  state variables used when the model run starts (store levels [mm] and Unit Hydrograph (UH) storages [mm])
-     &                             NOutputs     , ! [integer] number of output series
-     &                             IndOutputs   , ! [integer] indices of output series
-                                 ! outputs
-     &                             Outputs      , ! [double]  output series
-     &                             StateEnd     ) ! [double]  state variables at the end of the model run (store levels [mm] and Unit Hydrograph (UH) storages [mm])
+!------------------------------------------------------------------------------
+!    Subroutines relative to the annual GR6J model
+!------------------------------------------------------------------------------
+! TITLE   : airGR
+! PROJECT : airGR
+! FILE    : frun_GR6J.f
+!------------------------------------------------------------------------------
+! AUTHORS
+! Original code: R. Pushpalatha
+! Cleaning and formatting for airGR: L. Coron
+! Further cleaning: G. Thirel
+!------------------------------------------------------------------------------
+! Creation date: 2010
+! Last modified: 21/11/2019
+!------------------------------------------------------------------------------
+! REFERENCES
+! Pushpalatha, R., C. Perrin, N. Le Moine, T. Mathevet and V. Andréassian! 
+! (2011). A downward structural sensitivity analysis of hydrological models to 
+! improve low-flow simulation. Journal of Hydrology, 411(1-2), 66-76. 
+! doi:10.1016/j.jhydrol.2011.09.034.
+!------------------------------------------------------------------------------
+! Quick description of public procedures:
+!         1. frun_gr6j
+!         2. MOD_GR6J
+!------------------------------------------------------------------------------
+
+
+      SUBROUTINE frun_gr6j(LInputs,InputsPrecip,InputsPE,NParam,Param,
+     &                     NStates,StateStart,NOutputs,IndOutputs,
+     &                     Outputs,StateEnd)
+! Subroutine that initializes GR6J, get its parameters, performs the call 
+! to the MOD_GR6J 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/day]
+!       InputsPE     ! Vector of real, input series of potential evapotranspiration (PE) [mm/day]
+!       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] and Unit Hydrograph (UH) storages [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] and Unit Hydrograph (UH) storages [mm])
 
 
       !DEC$ ATTRIBUTES DLLEXPORT :: frun_gr6j
 
 
       Implicit None
-      ! input and output variables
+
+      !! dummies
+      ! in
       integer, intent(in) :: LInputs,NParam,NStates,NOutputs
-      doubleprecision, dimension(LInputs)  :: InputsPrecip
-      doubleprecision, dimension(LInputs)  :: InputsPE
-      doubleprecision, dimension(NParam)   :: Param
-      doubleprecision, dimension(NStates)  :: StateStart
-      doubleprecision, dimension(NStates)  :: StateEnd
-      integer, dimension(NOutputs) :: IndOutputs
-      doubleprecision, dimension(LInputs,NOutputs) :: Outputs
-
-      ! parameters, internal states and variables
-      integer NH,NMISC
-      parameter (NH=20,NMISC=30)
-      doubleprecision St(3), StUH1(NH), StUH2(2*NH)
-      doubleprecision OrdUH1(NH), OrdUH2(2*NH)
-      doubleprecision MISC(NMISC)
-      doubleprecision D
-      doubleprecision P1,E,Q
-      integer I,K
+      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 :: NH=20,NMISC=30
+      doubleprecision, dimension(3)     :: St
+      doubleprecision, dimension(NH)    :: StUH1, OrdUH1
+      doubleprecision, dimension(2*NH)  :: StUH2, OrdUH2
+      doubleprecision, dimension(NMISC) :: MISC
+      doubleprecision :: D,P1,E,Q
 
       !--------------------------------------------------------------
       ! Initializations
@@ -126,40 +158,48 @@
 !**********************************************************************
       SUBROUTINE MOD_GR6J(St,StUH1,StUH2,OrdUH1,OrdUH2,Param,P1,E,Q,
      &MISC)
-! Run on a single time step with the GR6J model
+! Calculation of streamflow on a single time step (day) with the GR6J model
 ! Inputs:
-!       St     Vector of model states in stores at the beginning of the time step [mm]
-!       StUH1  Vector of model states in Unit Hydrograph 1 at the beginning of the time step [mm]
-!       StUH2  Vector of model states in Unit Hydrograph 2 at the beginning of the time step [mm]
-!       OrdUH1 Vector of ordinates in UH1 [-]
-!       OrdUH2 Vector of ordinates in UH2 [-]
-!       Param  Vector of model parameters [various units]
-!       P1     Value of rainfall during the time step [mm]
-!       E      Value of potential evapotranspiration during the time step [mm]
+!       St     Vector of real, model states in stores at the beginning of the time step [mm]
+!       StUH1  Vector of real, model states in Unit Hydrograph 1 at the beginning of the time step [mm]
+!       StUH2  Vector of real, model states in Unit Hydrograph 2 at the beginning of the time step [mm]
+!       OrdUH1 Vector of real, ordinates in UH1 [-]
+!       OrdUH2 Vector of real, ordinates in UH2 [-]
+!       Param  Vector of real, model parameters [various units]
+!       P1     Real, value of rainfall during the time step [mm]
+!       E      Real, value of potential evapotranspiration during the time step [mm]
 ! Outputs:
-!       St     Vector of model states in stores at the end of the time step [mm]
-!       StUH1  Vector of model states in Unit Hydrograph 1 at the end of the time step [mm]
-!       StUH2  Vector of model states in Unit Hydrograph 2 at the end of the time step [mm]
-!       Q      Value of simulated flow at the catchment outlet for the time step [mm/day]
-!       MISC   Vector of model outputs for the time step [mm or mm/day]
+!       St     Vector of real, model states in stores at the end of the time step [mm]
+!       StUH1  Vector of real, model states in Unit Hydrograph 1 at the end of the time step [mm]
+!       StUH2  Vector of real, model states in Unit Hydrograph 2 at the end of the time step [mm]
+!       Q      Real, value of simulated flow at the catchment outlet for the time step [mm/day]
+!       MISC   Vector of real, model outputs for the time step [mm or mm/day]
 !**********************************************************************
       Implicit None
-      INTEGER NH,NMISC,NParam
-      PARAMETER (NH=20,NMISC=30)
-      PARAMETER (NParam=6)
-      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,PRUH1,PRUH2,EXCH,QR,QD,QRExp
-      DOUBLEPRECISION AE,AEXCH1,AEXCH2
-      INTEGER K
-
-      DATA B/0.9/
-      DATA C/0.4/
-      DOUBLEPRECISION TWS, Sr, Rr ! speed-up
+
+      !! locals
+      integer, parameter :: NParam=6,NMISC=30,NH=20
+      doubleprecision :: A,B,C,EN,ER,PN,PR,PS,WS,tanHyp,AR
+      doubleprecision :: PERC,PRUH1,PRUH2,EXCH,QR,QD,QRExp
+      doubleprecision :: AE,AEXCH1,AEXCH2
+      integer :: K
+      data B/0.9/
+      data C/0.4/
+      doubleprecision :: TWS, Sr, Rr ! speed-up
+
+      !! dummies
+      ! in
+      doubleprecision, dimension(NParam), intent(in)  :: Param
+      doubleprecision, intent(in) :: P1,E
+      doubleprecision, dimension(NH), intent(inout)   :: OrdUH1
+      doubleprecision, dimension(2*NH), intent(inout) :: OrdUH2
+      ! inout
+      doubleprecision, dimension(3), intent(inout)    :: St
+      doubleprecision, dimension(NH), intent(inout)   :: StUH1
+      doubleprecision, dimension(2*NH), intent(inout) :: StUH2
+      ! out
+      doubleprecision, intent(out) :: Q
+      doubleprecision, dimension(NMISC), intent(out)  :: MISC
 
       A=Param(1)
 
diff --git a/src/utils_D.f b/src/utils_D.f
index eaf0d6ff9cd49091dfc13bd34b4820721fa29d67..d12b96d3cc119e36bad19395595f86e68544a03e 100644
--- a/src/utils_D.f
+++ b/src/utils_D.f
@@ -1,69 +1,123 @@
+!------------------------------------------------------------------------------
+!    Subroutines relative to the daily unit hydrographs and hyperbolic
+!    tangent calculation
+!------------------------------------------------------------------------------
+! TITLE   : airGR
+! PROJECT : airGR
+! FILE    : utils_D.f
+!------------------------------------------------------------------------------
+! AUTHORS
+! Original code: Unknown soldier
+! Cleaning and formatting for airGR: L. Coron
+! Further cleaning: G. Thirel
+!------------------------------------------------------------------------------
+! Creation date: 2000
+! Last modified: 21/11/2019
+!------------------------------------------------------------------------------
+! REFERENCES
+!
+!------------------------------------------------------------------------------
+! Quick description of public procedures:
+!         1. UH1
+!         2. UH2
+!         3. SS1
+!         4. SS2
+!         5. tanHyp
+!------------------------------------------------------------------------------
+
+
 !**********************************************************************
       SUBROUTINE UH1(OrdUH1,C,D)
-! Computation of ordinates of GR unit hydrograph UH1 using successive differences on the S curve SS1
-! Inputs:
-!    C: time constant
-!    D: exponent
-! Outputs:
-!    OrdUH1: NH ordinates of discrete hydrograph
+! Subroutine that computes the ordinates of the daily GR unit hydrograph 
+! UH1 using successive differences on the S curve SS1
+! Inputs
+!    C ! Real, time constant
+!    D ! Real, exponent
+! Outputs
+!    OrdUH1: Vector of real, NH ordinates of the discrete hydrograph
 !**********************************************************************
+
       Implicit None
-      INTEGER NH
-      PARAMETER (NH=20)
-      DOUBLEPRECISION OrdUH1(NH)
-      DOUBLEPRECISION C,D,SS1
-      INTEGER I
 
+      !! locals
+      integer, parameter :: NH=20
+      doubleprecision :: SS1
+      integer :: I
+
+      !! dummies
+      ! in
+      double precision, intent(in) :: C,D
+      ! out
+      double precision, dimension (NH), intent(out) :: OrdUH1
+      
       DO I=1,NH
-      OrdUH1(I)=SS1(I,C,D)-SS1(I-1,C,D)
+        OrdUH1(I)=SS1(I,C,D)-SS1(I-1,C,D)
       ENDDO
       ENDSUBROUTINE
 
 
 !**********************************************************************
       SUBROUTINE UH2(OrdUH2,C,D)
-! Computation of ordinates of GR unit hydrograph HU2 using successive differences on the S curve SS2
-! Inputs:
-!    C: time constant
-!    D: exponent
-! Outputs:
-!    OrdUH2: 2*NH ordinates of discrete hydrograph
+! Subroutine that computes the ordinates of the daily GR unit hydrograph 
+! UH2 using successive differences on the S curve SS2
+! Inputs
+!    C ! Real, time constant
+!    D ! Real, exponent
+! Outputs
+!    OrdUH2: Vector of real, 2*NH ordinates of the discrete hydrograph
 !**********************************************************************
+
       Implicit None
-      INTEGER NH
-      PARAMETER (NH=20)
-      DOUBLEPRECISION OrdUH2(2*NH)
-      DOUBLEPRECISION C,D,SS2
-      INTEGER I
 
+      !! locals
+      integer, parameter :: NH=20
+      doubleprecision :: SS2
+      integer :: I
+
+      !! dummies
+      ! in
+      double precision, intent(in) :: C,D
+      ! out
+      double precision, dimension (2*NH), intent(out) :: OrdUH2
+      
       DO I =1,2*NH
-      OrdUH2(I)=SS2(I,C,D)-SS2(I-1,C,D)
+        OrdUH2(I)=SS2(I,C,D)-SS2(I-1,C,D)
       ENDDO
       ENDSUBROUTINE
 
 
 !**********************************************************************
       FUNCTION SS1(I,C,D)
-! Values of the S curve (cumulative HU curve) of GR unit hydrograph UH1
-! Inputs:
-!    C: time constant
-!    D: exponent
-!    I: time-step
-! Outputs:
-!    SS1: Values of the S curve for I
+! Function that computes the values of the S curve (cumulative UH 
+! curve) of GR unit hydrograph UH1
+! Inputs
+!    C ! Real, time constant
+!    D ! Real, exponent
+!    I ! Integer, time-step
+! Outputs
+!    SS1 ! Real, value of the S curve for I
 !**********************************************************************
+
       Implicit None
-      DOUBLEPRECISION C,D,SS1
-      INTEGER I,FI
+
+      !! dummies
+      ! in
+      doubleprecision, intent(in) :: C,D
+      integer, intent(in) :: I
+      ! out
+      doubleprecision ::SS1
+      
+      !! locals
+      integer :: FI
 
       FI=I
       IF(FI.LE.0) THEN
-      SS1=0.
-      RETURN
+        SS1=0.
+        RETURN
       ENDIF
       IF(FI.LT.C) THEN
-      SS1=(FI/C)**D
-      RETURN
+        SS1=(FI/C)**D
+        RETURN
       ENDIF
       SS1=1.
       ENDFUNCTION
@@ -71,30 +125,40 @@
 
 !**********************************************************************
       FUNCTION SS2(I,C,D)
-! Values of the S curve (cumulative HU curve) of GR unit hydrograph UH2
-! Inputs:
-!    C: time constant
-!    D: exponent
-!    I: time-step
-! Outputs:
-!    SS2: Values of the S curve for I
+! Function that computes the values of the S curve (cumulative UH 
+! curve) of GR unit hydrograph UH2
+! Inputs
+!    C ! Real, time constant
+!    D ! Real, exponent
+!    I ! Integer, time-step
+! Outputs
+!    SS2 ! Real, value of the S curve for I
 !**********************************************************************
+
       Implicit None
-      DOUBLEPRECISION C,D,SS2
-      INTEGER I,FI
+
+      !! dummies
+      ! in
+      doubleprecision, intent(in) :: C,D
+      integer, intent(in) :: I
+      ! out
+      doubleprecision ::SS2
+      
+      !! locals
+      integer :: FI
 
       FI=I
       IF(FI.LE.0) THEN
-      SS2=0.
-      RETURN
+        SS2=0.
+        RETURN
       ENDIF
       IF(FI.LE.C) THEN
-      SS2=0.5*(FI/C)**D
-      RETURN
+        SS2=0.5*(FI/C)**D
+        RETURN
       ENDIF
       IF(FI.LT.2.*C) THEN
-      SS2=1.-0.5*(2.-FI/C)**D
-      RETURN
+        SS2=1.-0.5*(2.-FI/C)**D
+        RETURN
       ENDIF
       SS2=1.
       ENDFUNCTION
@@ -102,10 +166,23 @@
 
 !**********************************************************************
       FUNCTION tanHyp(Val)
-! Computation of hyperbolic tangent
+! Function that calculates the hyperbolic tangent
 !**********************************************************************
+! Inputs
+!      Val ! Real, value
+! Outputs
+!      tanHyp ! Real, value of the hyperbolic tangent
+
       Implicit None
-      DOUBLEPRECISION Val,ValExp,tanHyp
+
+      !! dummies
+      ! in
+      doubleprecision, intent(in) :: Val
+      ! out
+      doubleprecision :: tanHyp
+      
+      !! locals
+      doubleprecision :: ValExp
 
       ValExp=EXP(Val)
       tanHyp=(ValExp - 1./ValExp)/(ValExp + 1./ValExp)
diff --git a/src/utils_H.f b/src/utils_H.f
index e3c2ada84d03cb62a29e943aa2b5b94bbf00e256..6d9021392dc6c6fbd8119a0532539036d75ab111 100644
--- a/src/utils_H.f
+++ b/src/utils_H.f
@@ -1,71 +1,121 @@
+!------------------------------------------------------------------------------
+!    Subroutines relative to the hourly unit hydrographs
+!------------------------------------------------------------------------------
+! TITLE   : airGR
+! PROJECT : airGR
+! FILE    : utils_H.f
+!------------------------------------------------------------------------------
+! AUTHORS
+! Original code: Unknown soldier
+! Cleaning and formatting for airGR: L. Coron
+! Further cleaning: G. Thirel
+!------------------------------------------------------------------------------
+! Creation date: 2000
+! Last modified: 22/11/2019
+!------------------------------------------------------------------------------
+! REFERENCES
+!
+!------------------------------------------------------------------------------
+! Quick description of public procedures:
+!         1. UH1_H
+!         2. UH2_H
+!         3. SS1_H
+!         4. SS2_H
+!------------------------------------------------------------------------------
 
 
 !**********************************************************************
       SUBROUTINE UH1_H(OrdUH1,C,D)
-! Computation of ordinates of GR unit hydrograph UH2 using successive differences on the S curve SS1
-! Inputs:
-!    C: time constant
-!    D: exponent
-! Outputs:
-!    OrdUH1: NH ordinates of discrete hydrograph
+! Subroutine that computes the ordinates of the hourly GR unit hydrograph 
+! UH1 using successive differences on the S curve SS1
+! Inputs
+!    C ! Real, time constant
+!    D ! Real, exponent
+! Outputs
+!    OrdUH1: Vector of real, NH ordinates of the discrete hydrograph
 !**********************************************************************
+
       Implicit None
-      INTEGER NH
-      PARAMETER (NH=480)
-      DOUBLEPRECISION OrdUH1(NH)
-      DOUBLEPRECISION C,D,SS1_H
-      INTEGER I
 
+      !! locals
+      integer, parameter :: NH=480
+      doubleprecision :: SS1_H
+      integer :: I
+
+      !! dummies
+      ! in
+      double precision, intent(in) :: C,D
+      ! out
+      double precision, dimension (NH), intent(out) :: OrdUH1
+      
       DO I=1,NH
-      OrdUH1(I)=SS1_H(I,C,D)-SS1_H(I-1,C,D)
+        OrdUH1(I)=SS1_H(I,C,D)-SS1_H(I-1,C,D)
       ENDDO
       ENDSUBROUTINE
 
 
 !**********************************************************************
       SUBROUTINE UH2_H(OrdUH2,C,D)
-! Computation of ordinates of GR unit hydrograph UH2 using successive differences on the S curve SS2
-! Inputs:
-!    C: time constant
-!    D: exponent
-! Outputs:
-!    OrdUH2: 2*NH ordinates of discrete hydrograph
+! Subroutine that computes the ordinates of the hourly GR unit hydrograph 
+! UH2 using successive differences on the S curve SS2
+! Inputs
+!    C ! Real, time constant
+!    D ! Real, exponent
+! Outputs
+!    OrdUH2: Vector of real, 2*NH ordinates of the discrete hydrograph
 !**********************************************************************
+
       Implicit None
-      INTEGER NH
-      PARAMETER (NH=480)
-      DOUBLEPRECISION OrdUH2(2*NH)
-      DOUBLEPRECISION C,D,SS2_H
-      INTEGER I
+
+      !! locals
+      integer, parameter :: NH=480
+      doubleprecision :: SS2_H
+      integer :: I
+
+      !! dummies
+      ! in
+      double precision, intent(in) :: C,D
+      ! out
+      double precision, dimension (2*NH), intent(out) :: OrdUH2
 
       DO I =1,2*NH
-      OrdUH2(I)=SS2_H(I,C,D)-SS2_H(I-1,C,D)
+        OrdUH2(I)=SS2_H(I,C,D)-SS2_H(I-1,C,D)
       ENDDO
       ENDSUBROUTINE
 
 
 !**********************************************************************
       FUNCTION SS1_H(I,C,D)
-! Values of the S curve (cumulative HU curve) of GR unit hydrograph HU1
-! Inputs:
-!    C: time constant
-!    D: exponent
-!    I: time-step
-! Outputs:
-!    SS1: Values of the S curve for I
+! Function that computes the values of the S curve (cumulative UH 
+! curve) of GR unit hydrograph UH1
+! Inputs
+!    C ! Real, time constant
+!    D ! Real, exponent
+!    I ! Integer, time-step
+! Outputs
+!    SS1_H ! Real, value of the S curve for I
 !**********************************************************************
+ 
       Implicit None
-      DOUBLEPRECISION C,D,SS1_H
-      INTEGER I,FI
+
+      !! dummies
+      ! in
+      doubleprecision, intent(in) :: C,D
+      integer, intent(in) :: I
+      ! out
+      doubleprecision ::SS1_H
+      
+      !! locals
+      integer :: FI
 
       FI=I
       IF(FI.LE.0) THEN
-      SS1_H=0.
-      RETURN
+        SS1_H=0.
+        RETURN
       ENDIF
       IF(FI.LT.C) THEN
-      SS1_H=(FI/C)**D
-      RETURN
+        SS1_H=(FI/C)**D
+        RETURN
       ENDIF
       SS1_H=1.
       ENDFUNCTION
@@ -73,17 +123,27 @@
 
 !**********************************************************************
       FUNCTION SS2_H(I,C,D)
-! Values of the S curve (cumulative HU curve) of GR unit hydrograph HU2
-! Inputs:
-!    C: time constant
-!    D: exponent
-!    I: time-step
-! Outputs:
-!    SS2: Values of the S curve for I
+! Function that computes the values of the S curve (cumulative UH 
+! curve) of GR unit hydrograph UH2
+! Inputs
+!    C ! Real, time constant
+!    D ! Real, exponent
+!    I ! Integer, time-step
+! Outputs
+!    SS2_H ! Real, value of the S curve for I
 !**********************************************************************
+
       Implicit None
-      DOUBLEPRECISION C,D,SS2_H
-      INTEGER I,FI
+
+      !! dummies
+      ! in
+      doubleprecision, intent(in) :: C,D
+      integer, intent(in) :: I
+      ! out
+      doubleprecision ::SS2_H
+      
+      !! locals
+      integer :: FI
 
       FI=I
       IF(FI.LE.0) THEN
@@ -91,12 +151,12 @@
       RETURN
       ENDIF
       IF(FI.LE.C) THEN
-      SS2_H=0.5*(FI/C)**D
-      RETURN
+        SS2_H=0.5*(FI/C)**D
+        RETURN
       ENDIF
       IF(FI.LT.2.*C) THEN
-      SS2_H=1.-0.5*(2.-FI/C)**D
-      RETURN
+        SS2_H=1.-0.5*(2.-FI/C)**D
+        RETURN
       ENDIF
       SS2_H=1.
       ENDFUNCTION