diff --git a/DESCRIPTION b/DESCRIPTION
index b5589b4e2bc17a542cefbad01ef9a6dbde37e7a1..6ce999db7bcb25a833a06af19c05391c444a5676 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: airGR
 Type: Package
 Title: Suite of GR Hydrological Models for Precipitation-Runoff Modelling
-Version: 1.3.2.67
+Version: 1.3.2.68
 Date: 2019-12-02
 Authors@R: c(
   person("Laurent", "Coron", role = c("aut", "trl"), comment = c(ORCID = "0000-0002-1503-6204")),
diff --git a/NEWS.md b/NEWS.md
index e95757284c35656c6246f36c5ae3f10ca43c5935..85fefe84a68a623d94eac4cba9a1c518c4422d79 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -2,7 +2,7 @@
 
 
 
-### 7 Release Notes (2019-12-02)
+### 1.3.2.68 Release Notes (2019-12-02)
 
 
 #### New features
@@ -21,6 +21,7 @@
 #### Minor user-visible changes
 
 - Added the diagram of GR2M in the <code>RunModel_GR2M ()</code> documentation.
+- Fotran code cleaned and translated from F77 to F90.
 
 
 #### CRAN-compatibility updates
diff --git a/src/frun_CEMANEIGE.f b/src/frun_CEMANEIGE.f90
similarity index 95%
rename from src/frun_CEMANEIGE.f
rename to src/frun_CEMANEIGE.f90
index 1ff16a83b29ff01aaf956d647e7d42f6e0af9288..b19e1aae4c9598601a25e012454d500bacd9b58b 100644
--- a/src/frun_CEMANEIGE.f
+++ b/src/frun_CEMANEIGE.f90
@@ -34,11 +34,11 @@
 !------------------------------------------------------------------------------
 
 
-      SUBROUTINE frun_cemaneige(LInputs,InputsPrecip,
-     &                          InputsFracSolidPrecip,InputsTemp,
-     &                          MeanAnSolidPrecip,NParam,Param,NStates,
-     &                          StateStart,IsHyst,NOutputs,IndOutputs,
-     &                          Outputs,StateEnd)
+      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
@@ -69,8 +69,7 @@
       integer, intent(in) :: LInputs,NParam,NStates,NOutputs
       doubleprecision, intent(in) :: MeanAnSolidPrecip
       doubleprecision, intent(in), dimension(LInputs) :: InputsPrecip
-      doubleprecision, intent(in), dimension(LInputs) :: 
-     & InputsFracSolidPrecip
+      doubleprecision, intent(in), dimension(LInputs) :: InputsFracSolidPrecip
       doubleprecision, intent(in), dimension(LInputs) :: InputsTemp
       doubleprecision, intent(in), dimension(NParam)  :: Param
       doubleprecision, intent(in), dimension(NStates) :: StateStart
@@ -78,8 +77,7 @@
       integer, intent(in), dimension(NOutputs) :: IndOutputs
       ! out
       doubleprecision, intent(out), dimension(NStates) :: StateEnd
-      doubleprecision, intent(out), dimension(LInputs,NOutputs) :: 
-     & Outputs
+      doubleprecision, intent(out), dimension(LInputs,NOutputs) :: Outputs
 
       !! locals
       logical :: IsHystBool         ! TRUE if linear hysteresis is used, FALSE otherwise
@@ -221,5 +219,5 @@
 
       RETURN
 
-      ENDSUBROUTINE
+      END SUBROUTINE
 
diff --git a/src/frun_GR1A.f b/src/frun_GR1A.f90
similarity index 96%
rename from src/frun_GR1A.f
rename to src/frun_GR1A.f90
index e35857c6b043cb19a84e4fa47d7891772824217a..a6ff7cf89d0b13b7bbb77597a225351394ac68fc 100644
--- a/src/frun_GR1A.f
+++ b/src/frun_GR1A.f90
@@ -24,9 +24,9 @@
 !------------------------------------------------------------------------------
 
 
-      SUBROUTINE frun_gr1a(LInputs,InputsPrecip,InputsPE,NParam,Param,
-     &                     NStates,StateStart,NOutputs,IndOutputs,
-     &                     Outputs,StateEnd)
+      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
@@ -59,8 +59,7 @@
       integer, dimension(NOutputs),        intent(in) :: IndOutputs
       ! out
       doubleprecision, dimension(NStates), intent(out) :: StateEnd
-      doubleprecision, dimension(LInputs,NOutputs), 
-     & intent(out) :: Outputs
+      doubleprecision, dimension(LInputs,NOutputs), intent(out) :: Outputs
       
       !! locals
       integer :: I,K
@@ -155,6 +154,6 @@
       MISC( 3)=Q             ! Qsim   ! [numeric] simulated outflow at catchment outlet [mm/year]
 
 
-      ENDSUBROUTINE
+      END SUBROUTINE
 
 
diff --git a/src/frun_GR2M.f b/src/frun_GR2M.f90
similarity index 97%
rename from src/frun_GR2M.f
rename to src/frun_GR2M.f90
index e2317c9f1ff14f0cbd5eb1964c77f5c5e7a353bf..9be6de91b8e8522b920a7d9a7f1158ac06dc0baa 100644
--- a/src/frun_GR2M.f
+++ b/src/frun_GR2M.f90
@@ -28,9 +28,9 @@
 !------------------------------------------------------------------------------
 
 
-      SUBROUTINE frun_gr2m(LInputs,InputsPrecip,InputsPE,NParam,Param,
-     &                     NStates,StateStart,NOutputs,IndOutputs,
-     &                     Outputs,StateEnd)
+      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
@@ -63,8 +63,7 @@
       integer, dimension(NOutputs),        intent(in) :: IndOutputs
       ! out
       doubleprecision, dimension(NStates), intent(out) :: StateEnd
-      doubleprecision, dimension(LInputs,NOutputs), 
-     & intent(out) :: Outputs
+      doubleprecision, dimension(LInputs,NOutputs), intent(out) :: Outputs
       
       !! locals
       integer :: I,K
@@ -220,6 +219,6 @@
       MISC(10)=Q             ! Qsim   ! [numeric] simulated outflow at catchment outlet [mm/month]
 
 
-      ENDSUBROUTINE
+      END SUBROUTINE
 
 
diff --git a/src/frun_GR4H.f b/src/frun_GR4H.f90
similarity index 97%
rename from src/frun_GR4H.f
rename to src/frun_GR4H.f90
index e4cf854cecb3f46ea1237ee26d5e30289fe6d55a..ec12cc787151ab74822469c44407ca9e945a16cd 100644
--- a/src/frun_GR4H.f
+++ b/src/frun_GR4H.f90
@@ -24,9 +24,9 @@
 !------------------------------------------------------------------------------
 
 
-      SUBROUTINE frun_gr4h(LInputs,InputsPrecip,InputsPE,NParam,Param,
-     &                     NStates,StateStart,NOutputs,IndOutputs,
-     &                     Outputs,StateEnd)
+      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
@@ -59,8 +59,7 @@
       integer, dimension(NOutputs),        intent(in) :: IndOutputs
       ! out
       doubleprecision, dimension(NStates), intent(out) :: StateEnd
-      doubleprecision, dimension(LInputs,NOutputs), 
-     & intent(out) :: Outputs
+      doubleprecision, dimension(LInputs,NOutputs), intent(out) :: Outputs
       
       !! locals
       integer :: I,K
@@ -151,8 +150,7 @@
 
 
 !**********************************************************************
-      SUBROUTINE MOD_GR4H(St,StUH1,StUH2,OrdUH1,OrdUH2,Param,P1,E,Q,
-     &MISC)
+      SUBROUTINE MOD_GR4H(St,StUH1,StUH2,OrdUH1,OrdUH2,Param,P1,E,Q,MISC)
 ! Calculation of streamflow on a single time step (hour) with the GR4H model
 ! Inputs:
 !       St     Vector of real, model states in stores at the beginning of the time step [mm]
@@ -319,6 +317,6 @@
 
 
 
-      ENDSUBROUTINE
+      END SUBROUTINE
 
 
diff --git a/src/frun_GR4J.f b/src/frun_GR4J.f90
similarity index 97%
rename from src/frun_GR4J.f
rename to src/frun_GR4J.f90
index f4484dba51c5bf6c0f09d87e0eeef31ecebf8863..fa486f31dd290c4501ca8f55f5251c88d2eca985 100644
--- a/src/frun_GR4J.f
+++ b/src/frun_GR4J.f90
@@ -24,9 +24,9 @@
 !------------------------------------------------------------------------------
 
 
-      SUBROUTINE frun_gr4j(LInputs,InputsPrecip,InputsPE,NParam,Param,
-     &                     NStates,StateStart,NOutputs,IndOutputs,
-     &                     Outputs,StateEnd)
+      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
@@ -59,8 +59,7 @@
       integer, dimension(NOutputs),        intent(in) :: IndOutputs
       ! out
       doubleprecision, dimension(NStates), intent(out) :: StateEnd
-      doubleprecision, dimension(LInputs,NOutputs), 
-     & intent(out) :: Outputs
+      doubleprecision, dimension(LInputs,NOutputs), intent(out) :: Outputs
       
       !! locals
       integer :: I,K
@@ -151,8 +150,7 @@
 
 
 !**********************************************************************
-      SUBROUTINE MOD_GR4J(St,StUH1,StUH2,OrdUH1,OrdUH2,Param,P1,E,Q,
-     &MISC)
+      SUBROUTINE MOD_GR4J(St,StUH1,StUH2,OrdUH1,OrdUH2,Param,P1,E,Q,MISC)
 ! Calculation of streamflow on a single time step (day) with the GR4J model
 ! Inputs:
 !       St     Vector of real, model states in stores at the beginning of the time step [mm]
@@ -316,6 +314,6 @@
 
 
 
-      ENDSUBROUTINE
+      END SUBROUTINE
 
 
diff --git a/src/frun_GR5J.f b/src/frun_GR5J.f90
similarity index 97%
rename from src/frun_GR5J.f
rename to src/frun_GR5J.f90
index ecb8f7fd4a716ae15b74b9018d9f1446df11ca1b..d7b12ba03a16c445711b34fbaed6bb3261bcd266 100644
--- a/src/frun_GR5J.f
+++ b/src/frun_GR5J.f90
@@ -29,9 +29,9 @@
 !------------------------------------------------------------------------------
 
 
-      SUBROUTINE frun_gr5j(LInputs,InputsPrecip,InputsPE,NParam,Param,
-     &                     NStates,StateStart,NOutputs,IndOutputs,
-     &                     Outputs,StateEnd)
+      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
@@ -64,8 +64,7 @@
       integer, dimension(NOutputs),        intent(in) :: IndOutputs
       ! out
       doubleprecision, dimension(NStates), intent(out) :: StateEnd
-      doubleprecision, dimension(LInputs,NOutputs), 
-     & intent(out) :: Outputs
+      doubleprecision, dimension(LInputs,NOutputs), intent(out) :: Outputs
       
       !! locals
       integer :: I,K
@@ -148,8 +147,7 @@
 
 
 !**********************************************************************
-      SUBROUTINE MOD_GR5J(St,StUH2,OrdUH2,Param,P1,E,Q,
-     &MISC)
+      SUBROUTINE MOD_GR5J(St,StUH2,OrdUH2,Param,P1,E,Q,MISC)
 ! Calculation of streamflow on a single time step (day) with the GR5J model
 ! Inputs:
 !       St     Vector of real, model states in stores at the beginning of the time step [mm]
@@ -301,6 +299,6 @@
       MISC(18)=Q             ! Qsim   ! simulated outflow at catchment outlet [mm/day]
 
 
-      ENDSUBROUTINE
+      END SUBROUTINE
 
 
diff --git a/src/frun_GR6J.f b/src/frun_GR6J.f90
similarity index 98%
rename from src/frun_GR6J.f
rename to src/frun_GR6J.f90
index 5d3600ac29794068b7229546a86f9712c590220a..751d9114b19b7d2342bb84754327db83bbb53808 100644
--- a/src/frun_GR6J.f
+++ b/src/frun_GR6J.f90
@@ -25,9 +25,9 @@
 !------------------------------------------------------------------------------
 
 
-      SUBROUTINE frun_gr6j(LInputs,InputsPrecip,InputsPE,NParam,Param,
-     &                     NStates,StateStart,NOutputs,IndOutputs,
-     &                     Outputs,StateEnd)
+      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
@@ -60,8 +60,7 @@
       integer, dimension(NOutputs),        intent(in) :: IndOutputs
       ! out
       doubleprecision, dimension(NStates), intent(out) :: StateEnd
-      doubleprecision, dimension(LInputs,NOutputs), 
-     & intent(out) :: Outputs
+      doubleprecision, dimension(LInputs,NOutputs), intent(out) :: Outputs
       
       !! locals
       integer :: I,K
@@ -156,8 +155,7 @@
 
 
 !**********************************************************************
-      SUBROUTINE MOD_GR6J(St,StUH1,StUH2,OrdUH1,OrdUH2,Param,P1,E,Q,
-     &MISC)
+      SUBROUTINE MOD_GR6J(St,StUH1,StUH2,OrdUH1,OrdUH2,Param,P1,E,Q,MISC)
 ! Calculation of streamflow on a single time step (day) with the GR6J model
 ! Inputs:
 !       St     Vector of real, model states in stores at the beginning of the time step [mm]
@@ -345,6 +343,6 @@
       MISC(20)=Q                  ! Qsim   ! simulated outflow at catchment outlet [mm/day]
 
 
-      ENDSUBROUTINE
+      END SUBROUTINE
 
 
diff --git a/src/utils_D.f b/src/utils_D.f90
similarity index 100%
rename from src/utils_D.f
rename to src/utils_D.f90
diff --git a/src/utils_H.f b/src/utils_H.f90
similarity index 100%
rename from src/utils_H.f
rename to src/utils_H.f90