diff --git a/DESCRIPTION b/DESCRIPTION
index ca35bec4c2ad471bdf18db50b37d20d89a9de1f4..50051874d2d3f569095d27768451e1bc75b84bb6 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.6.3.68
+Version: 1.6.3.69
 Date: 2020-11-23
 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 bc9380ed5c47fa4015bc044149168d8b4f441168..e53fabac733471ca8529d881f3af6e1dac871fe5 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -4,7 +4,7 @@
 
 
 
-### 1.6.3.68 Release Notes (2020-11-23)
+### 1.6.3.69 Release Notes (2020-11-23)
 
 #### New features
 
diff --git a/src/frun_ETP.f90 b/src/frun_ETP.f90
index c0c6f20ddf9dd5936066c43ab5b9ef34ae9e09d0..9d064c3bc20eb00697f150d12aa9eaa7e539318f 100644
--- a/src/frun_ETP.f90
+++ b/src/frun_ETP.f90
@@ -1,6 +1,45 @@
+!------------------------------------------------------------------------------
+!    Subroutines relative to the Oudin potential evapotranspiration (PE) formula
+!------------------------------------------------------------------------------
+! TITLE   : airGR
+! PROJECT : airGR
+! FILE    : frun_ETP.f90
+!------------------------------------------------------------------------------
+! AUTHORS
+! Original code: L. Oudin
+! Cleaning and formatting for airGR: Fr. Bourgin
+! Further cleaning: O. Delaigue, G. Thirel
+!------------------------------------------------------------------------------
+! Creation date: 2004
+! Last modified: 20/10/2020
+!------------------------------------------------------------------------------
+! REFERENCES
+! Oudin, L., Hervieu, F., Michel, C., Perrin, C., Andréassian, V.,
+! Anctil, F. and Loumagne, C., 2005. Which potential evapotranspiration
+! input for a rainfall-runoff model? Part 2 - Towards a simple and
+! efficient PE model for rainfall-runoff modelling. Journal of Hydrology
+! 303(1-4), 290-306.
+!------------------------------------------------------------------------------
+! Quick description of public procedures:
+!         1. frun_etp_oudin
+!         2. PE_OUDIN
+!------------------------------------------------------------------------------
+
+
 !*******************************************************************************
-      SUBROUTINE frun_etp_oudin(LInputs,InputsLAT,InputsTT,InputsJJ,OutputsETP)
+      SUBROUTINE frun_etp_oudin(LInputs,InputsLAT,InputsTemp,InputsJJ,OutputsPE)
 !*******************************************************************************
+! Subroutine that performs the call to the PE_OUDIN subroutine at each time step, 
+! and stores the final values
+! Inputs
+!       LInputs      ! Integer, length of input and output series
+!       InputsLAT    ! Vector of real, input series of latitude [rad]
+!       InputsTemp   ! Vector of real, input series of air mean temperature [degC]
+!       InputsJJ     ! Vector of real, input series of Julian day [-]
+! Outputs
+!       OutputsPE    ! Vector of real, output series of potential evapotranspiration (PE) [mm/time step]
+
+
 
       !DEC$ ATTRIBUTES DLLEXPORT :: frun_etp_oudin
 
@@ -10,27 +49,27 @@
       ! in
       integer, intent(in) :: LInputs
       doubleprecision, dimension(LInputs), intent(in) :: InputsLAT
-      doubleprecision, dimension(LInputs), intent(in) :: InputsTT
+      doubleprecision, dimension(LInputs), intent(in) :: InputsTemp
       doubleprecision, dimension(LInputs), intent(in) :: InputsJJ
       
       ! out
-      doubleprecision, dimension(LInputs), intent(out) :: OutputsETP
+      doubleprecision, dimension(LInputs), intent(out) :: OutputsPE
     
       !! locals
-      integer :: k
-      real :: FI, tt, jj, ETPoud
+      integer         :: k
+      doubleprecision :: FI, tt, jj, PEoud
       
       !--------------------------------------------------------------
       ! Time loop
       !--------------------------------------------------------------
-      DO k=1,LInputs
-        tt = InputsTT(k)
+      DO k = 1, LInputs
+        tt = InputsTemp(k)
         jj = InputsJJ(k)
-        FI = InputsLAT(k) / 57.296
+        FI = InputsLAT(k)! 
         !model run on one time step
-        CALL ETP_OUDIN(FI,tt,jj,ETPoud)
+        CALL PE_OUDIN(FI, tt, jj, PEoud)
         !storage of outputs
-        OutputsETP(k) = ETPoud
+        OutputsPE(k) = PEoud
       ENDDO
 
       RETURN
@@ -38,14 +77,22 @@
       ENDSUBROUTINE
 
 
+
+
+
+!################################################################################################################################
+
+
+
+
 !*******************************************************************************
-      SUBROUTINE ETP_OUDIN(FI,DT,JD,DPE)
+      SUBROUTINE PE_OUDIN(FI,DT,JD,DPE)
 !*******************************************************************************
-!     This subroutine calculates daily potential evapotranspiration (DPE)
-!     using daily temperature and daily extra-atmospheric global radiation
+!     Calculation of potential evapotranspiration (DPE) on a single time step
+!     using air temperature and daily extra-atmospheric global radiation
 !     (that depends only on Julian day)
 !
-!     The PE formula is is that described in:
+!     The PE formula is described in:
 !     Oudin, L., Hervieu, F., Michel, C., Perrin, C., Andréassian, V.,
 !     Anctil, F. and Loumagne, C., 2005. Which potential evapotranspiration
 !     input for a rainfall-runoff model? Part 2 - Towards a simple and
@@ -58,67 +105,68 @@
 !     of hydrology. Journal of Hydrology 66 (1/4), 1-76.
 !
 !***************************************************************
-!     Inputs:
-!     xLAT: Latitude in decimal degrees
-!     DT: Temperature in degree C
-!     JD: Julian day
+! Inputs
+!       FI           ! Latitude [rad]
+!       DT           ! Air Temperature [degC]
+!       JD           ! Julian day [-]
 !
-!     Output:
-!     DPE: Daily potential evapotranspiration in mm
+! Outputs
+!       DPE          ! Potential evapotranspiration [mm/time step]
 !***************************************************************
       IMPLICIT NONE
 
-      REAL :: xLAT, FI, COSFI, TETA, COSTETA, COSGZ, GZ, COSGZ2
-      REAL :: SINGZ, COSOM, COSOM2, SINOM, COSPZ, OM, GE
-      REAL ::  ETA, DPE, DT, JD, RD
-
-!      DATA RD/57.296/
+      !! dummies
+      ! in
+      doubleprecision, intent(in) :: FI, DT, JD
+      ! out
+      doubleprecision, intent(out) :: DPE
+      !! locals
+      doubleprecision :: COSFI, TETA, COSTETA, COSGZ, GZ, COSGZ2
+      doubleprecision :: SINGZ, COSOM, COSOM2, SINOM, COSPZ, OM, GE
+      doubleprecision :: ETA
 
 ! Calculation of extra-atmospheric global radiation (Appendix C in Morton
 ! (1983), Eq. C-6 to C-11, p.60-61)
-! Converts latitude in radians
-!      FI=xLAT/RD
       COSFI=COS(FI)
-!      AFI=ABS(xLAT/42.)
 
 ! TETA: Declination of the sun in radians
       TETA=0.4093*SIN(JD/58.1-1.405)
       COSTETA=COS(TETA)
-      COSGZ=MAX(0.001,COS(FI-TETA))
+      COSGZ=MAX(0.001d0,COS(FI-TETA))
 
 ! GZ: Noon angular zenith distance of the sun
       GZ=ACOS(COSGZ)
       COSGZ2=COSGZ*COSGZ
 
-      IF(COSGZ2.GE.1.)THEN
-      SINGZ=0.
+      IF(COSGZ2.GE.1.) THEN
+        SINGZ=0.
       ELSE
-      SINGZ=SQRT(1.-COSGZ2)
+        SINGZ=SQRT(1.-COSGZ2)
       ENDIF
 
       COSOM=1.-COSGZ/COSFI/COSTETA
-      IF(COSOM.LT.-1.)COSOM=-1.
-      IF(COSOM.GT.1.)COSOM=1.
+      IF(COSOM.LT.-1.) COSOM=-1.
+      IF(COSOM.GT.1.) COSOM=1.
       COSOM2=COSOM*COSOM
-      IF(COSOM2.GE.1.)THEN
-      SINOM=0.
+      IF(COSOM2.GE.1.) THEN
+        SINOM=0.
       ELSE
-      SINOM=SQRT(1.-COSOM2)
+        SINOM=SQRT(1.-COSOM2)
       ENDIF
 
       OM=ACOS(COSOM)
 ! PZ: Average angular zenith distance of the sun
       COSPZ=COSGZ+COSFI*COSTETA*(SINOM/OM-1.)
-      IF(COSPZ.LT.0.001)COSPZ=0.001
+      IF(COSPZ.LT.0.001) COSPZ=0.001
 ! ETA: Radius vector of the sun
       ETA=1.+COS(JD/58.1)/30.
 ! GE: extra-atmospheric global radiation
       GE=446.*OM*COSPZ*ETA
 
 ! Daily PE by Oudin et al. (2006) formula:
-      DPE=MAX(0.,GE*(DT+5.)/100./28.5)
+      DPE=MAX(0.d0,GE*(DT+5.)/100./28.5)
 
       RETURN
 
-	END SUBROUTINE ETP_OUDIN
+      END SUBROUTINE PE_OUDIN
 !*******************************************************************************