diff --git a/src/MNH/default_desfmn.f90 b/src/MNH/default_desfmn.f90
old mode 100644
new mode 100755
index e460a683600a397cb7a0eec78a68182a9c6f07a6..6dd86e70f07ffcf953be0858722ff751cf83cd3d
--- a/src/MNH/default_desfmn.f90
+++ b/src/MNH/default_desfmn.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 1994-2020 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
!MNH_LIC for details. version 1.
@@ -253,6 +253,10 @@ USE MODD_CONDSAMP
USE MODD_MEAN_FIELD
USE MODD_DRAGTREE_n
USE MODD_DRAGBLDG_n
+USE MODD_EOL_MAIN
+USE MODD_EOL_ADNR
+USE MODD_EOL_ALM
+USE MODD_EOL_SHARED_IO
!
!
USE MODD_PARAM_LIMA, ONLY : LCOLD, LNUCL, LSEDI, LHHONI, LSNOW, LHAIL, LMEYERS,&
@@ -531,6 +535,32 @@ XVDEPOTREE = 0.02 ! 2 cm/s
!
LDRAGBLDG = .FALSE.
!
+!* 10d. SET DEFAULT VALUES FOR MODD_EOL* :
+! ----------------------------------
+!
+! 10d.i) MODD_EOL_MAIN
+!
+LMAIN_EOL = .FALSE.
+CMETH_EOL = 'ADNR'
+CSMEAR = '3LIN'
+NMODEL_EOL = 1
+!
+! 10d.ii) MODD_EOL_SHARED_IO
+!
+CFARM_CSVDATA = 'data_farm.csv'
+CTURBINE_CSVDATA = 'data_turbine.csv'
+CBLADE_CSVDATA = 'data_blade.csv'
+CAIRFOIL_CSVDATA = 'data_airfoil.csv'
+!
+CINTERP = 'CLS'
+!
+! 10d.iii) MODD_EOL_ALM
+!
+NNB_BLAELT = 42
+LTIMESPLIT = .FALSE.
+LTIPLOSSG = .TRUE.
+LTECOUTPTS = .FALSE.
+!
!-------------------------------------------------------------------------------
!
!* 11. SET DEFAULT VALUES FOR MODD_BUDGET :
diff --git a/src/MNH/diag.f90 b/src/MNH/diag.f90
old mode 100644
new mode 100755
index 24d921c97a8ed136ccbfe6acd0e29110e9b8b4a0..27d3244abd5f99a5aa06b83c1d033693ec57403b
--- a/src/MNH/diag.f90
+++ b/src/MNH/diag.f90
@@ -190,7 +190,7 @@ CHARACTER (LEN=4) :: YTURB ! initial flag to call to turbulence schemes
! CHARACTER (LEN=40) :: YFMT,YFMT2! format for cpu analysis printing
INTEGER :: ILUOUT0 ! Logical unit number for the output listing
REAL(kind=MNHTIME), DIMENSION(2) :: ZTIME0, ZTIME1, ZTIME2, ZRAD, ZDCONV, ZSHADOWS, ZGROUND, &
- ZTRACER, ZDRAG, ZTURB, ZMAFL, ZCHEM, ZTIME_BU ! CPU times
+ ZTRACER, ZDRAG, ZTURB, ZMAFL, ZCHEM, ZTIME_BU, ZEOL ! CPU times
REAL(kind=MNHTIME), DIMENSION(2) :: ZSTART, ZINIT, ZWRIT, ZBALL, ZPHYS, ZSURF, ZWRITS, ZTRAJ ! storing variables
INTEGER(KIND=LFIINT) :: INPRAR ! number of articles predicted in the LFIFM file
INTEGER :: ISTEPBAL ! loop indice for balloons and aircraft
@@ -691,11 +691,12 @@ ZTURB = 0.0_MNHTIME
ZDRAG = 0.0_MNHTIME
ZMAFL = 0.0_MNHTIME
ZCHEM = 0.0_MNHTIME
+ZEOL = 0.0_MNHTIME
XTIME_LES = 0.0_MNHTIME
XTIME_LES_BU_PROCESS = 0.0_MNHTIME
XTIME_BU_PROCESS = 0.0_MNHTIME
CALL PHYS_PARAM_n( 1, TOUTDATAFILE, &
- ZRAD, ZSHADOWS, ZDCONV, ZGROUND, ZMAFL, ZDRAG, &
+ ZRAD, ZSHADOWS, ZDCONV, ZGROUND, ZMAFL, ZDRAG,ZEOL, &
ZTURB, ZTRACER, ZTIME_BU, ZWETDEPAER, GMASKkids, GCLOUD_ONLY )
WRITE(ILUOUT0,*) 'DIAG AFTER PHYS_PARAM1'
IF (LCHEMDIAG) THEN
diff --git a/src/MNH/eol_adnr.f90 b/src/MNH/eol_adnr.f90
new file mode 100755
index 0000000000000000000000000000000000000000..79ae8a7ad9441ba200cf9036d2bc27accbe2841a
--- /dev/null
+++ b/src/MNH/eol_adnr.f90
@@ -0,0 +1,249 @@
+!MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
+!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
+!MNH_LIC for details. version 1.
+! #######################
+ MODULE MODI_EOL_ADNR
+! #######################
+!
+INTERFACE
+!
+SUBROUTINE EOL_ADNR(PDXX, PDYY, PDZZ, &
+ PRHO_M, PUT_M, &
+ PFX_RG )
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX,PDYY,PDZZ ! mesh size
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHO_M ! dry Density
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PUT_M ! Wind speed at mass point
+!
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PFX_RG ! Aerodynamic force (cartesian mesh..
+ ! .. x axis, global
+ ! ..frame)
+!
+!
+END SUBROUTINE EOL_ADNR
+!
+END INTERFACE
+!
+END MODULE MODI_EOL_ADNR
+!
+! ###################################################################
+ SUBROUTINE EOL_ADNR(PDXX, PDYY, PDZZ, &
+ PRHO_M, PUT_M, &
+ PFX_RG )
+! ###################################################################
+!
+!!**** *MODI_EOL_ADNR* -
+!!
+!! PURPOSE
+!! -------
+!! It is possible to include wind turbines parameterization in Meso-NH,
+!! and several methods are available. One of the models is the Non-
+!! Rotating Actuator Disk Non-Rotating model (ADNR). It allows to
+!! compute aerodynamic forces according the wind speed and the
+!! caracteristics of the wind turbine.
+!!
+!!** METHOD
+!! ------
+!! The actuator disc flow model, in this routine, is computed without
+!! rotation. It consists in applying a thrust force over the disc drawn
+!! by the blades. This aerodynamic force acts against the wind to disturb
+!! the flow.
+!!
+!! REFERENCE
+!! ---------
+!! PA. Joulin PhD Thesis. 2020.
+!!
+!! AUTHOR
+!! ------
+!! PA. Joulin *CNRM & IFPEN*
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 01/2017
+!! Modification 19/10/20 (PA. Joulin) Updated for a main version
+!!
+!!---------------------------------------------------------------
+!
+!
+!* 0. DECLARATIONS
+! ------------
+!
+! To work with ADNR
+USE MODD_EOL_ADNR
+!
+USE MODD_EOL_SHARED_IO, ONLY: CINTERP
+USE MODD_EOL_SHARED_IO, ONLY: XTHRUT
+USE MODI_EOL_MATHS, ONLY: INTERP_LIN8NB
+! To know the grid
+USE MODD_GRID_n, ONLY: XXHAT,XYHAT,XZS,XZZ
+USE MODE_ll, ONLY: GET_INDICE_ll
+USE MODD_PARAMETERS, ONLY: JPVEXT
+! To play with MPI
+USE MODD_VAR_ll, ONLY: NMNH_COMM_WORLD, IP
+use MODD_PRECISION, only: MNHREAL_MPI
+USE MODD_MPIF, ONLY: MPI_SUM
+!
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of dummy arguments :
+! Meso-NH
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX,PDYY,PDZZ ! Mesh size
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHO_M ! Dry Density * Jacobian
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PUT_M ! Wind speed at mass point
+! Wind turbine aerodynamic
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PFX_RG ! Aerodynamic force (cartesian mesh, x axis, RG frame) [N]
+!
+!* 0.2 Declarations of local variables :
+!
+! Indicies
+INTEGER :: IIB,IJB,IKB ! Begin of a CPU domain
+INTEGER :: IIE,IJE,IKE ! End of a CPU domain
+INTEGER :: IKU ! Vertical size of the domain
+INTEGER :: JI, JJ, JK ! Loop index
+INTEGER :: JROT ! Wind turbine index
+!
+! Wind
+REAL :: ZRHO_I ! Interpolated density [kg/m3]
+REAL :: ZUT_I ! Interpolated wind speed U (RG) [m/s]
+REAL, DIMENSION(SIZE(PUT_M,1),SIZE(PUT_M,2),SIZE(PUT_M,3)) :: ZZH ! True heigth to interpolate 8N
+!
+! Wind turbine
+REAL, DIMENSION(TFARM%NNB_TURBINES) :: ZTHRUT ! Thrust [N]
+! Geometry
+REAL :: ZY_DIST ! Distance Hub - Cell on Y [m]
+REAL :: ZZ_DIST ! Distance Hub - Cell on Z [m]
+REAL :: ZR_DIST ! Radial distance Hub - Cell [m]
+REAL, DIMENSION(3) :: ZPOS ! Element position [m]
+!
+!Numerical
+INTEGER :: IINFO ! code info return
+!
+!* 0.3 Implicit arguments
+!
+! A. From MODD_EOL_ADNR:
+!TYPE(FARM) :: TFARM
+!TYPE(TURBINE) :: TTURBINE
+!REAL, DIMENSION(:), ALLOCATABLE :: XA_INDU ! Induction factor [-]
+!REAL, DIMENSION(:), ALLOCATABLE :: XCT_D ! Adapted thrust coef (for U_d) [-]
+!
+!
+! B. From MODD_EOL_SHARED_IO:
+! for NAM_EOL_ADNR:
+!CHARACTER(LEN=100) :: CFARM_CSVDATA ! File to read, with farm data
+!CHARACTER(LEN=100) :: CTURBINE_CSVDATA ! File to read, turbine data
+!CHARACTER(LEN=3) :: CINTERP ! Interpolation method for wind speed
+! for outputs
+!REAL, DIMENSION(:), ALLOCATABLE :: XTHRUT ! Thrust [N]
+!
+!
+!-------------------------------------------------------------------------------
+!
+!* 1. INITIALIZATIONS
+! ---------------
+!
+!* 1.1 Subdomain (CPU) indices
+!
+CALL GET_INDICE_ll(IIB,IJB,IIE,IJE) ! Get begin and end domain index (CPU)
+IKU = SIZE(PUT_M,3) ! Top of the domain end index
+IKB=1+JPVEXT ! Vertical begin index
+IKE=IKU-JPVEXT ! Vertical end index
+!
+!* 1.2 Induction factor and adapted thrust coef (that will use U_disc)
+!
+DO JROT=1,TFARM%NNB_TURBINES
+ XA_INDU(JROT) = 0.5*(1-(1-TFARM%XCT_INF(JROT))**0.5)
+ XCT_D(JROT) = 4*XA_INDU(JROT)/(1-XA_INDU(JROT))
+END DO
+!
+CALL PRINTMER_CPU1('ADNR : At ', 1.3)
+!* 1.3 Inits
+!
+ZTHRUT(:) = 0.
+XTHRUT(:) = 0.
+!
+CALL PRINTMER_CPU1('ADNR : At ', 1.4)
+!* 1.4 Vertical coordinate in case of interpolation
+!
+IF (CINTERP=='8NB') THEN
+ DO JK=1,IKU-1
+ ZZH(:,:,JK) = (0.5*(XZZ(:,:,JK)+XZZ(:,:,JK+1))-XZS(:,:))
+ END DO
+ ZZH(:,:,IKU) = 2*ZZH(:,:,IKU-1) - ZZH(:,:,IKU-2)
+END IF
+!
+!-------------------------------------------------------------------------------
+!
+!* 2. COMPUTES THRUST FORCE THAT ACTS ON THE ROTOR DUE TO THE WIND
+! ------------------------------------------------------------
+!
+!* 2.1 Finding the position of wind turbines
+!
+! Loop over domain
+DO JK=IKB,IKE
+ DO JJ=IJB,IJE
+ DO JI=IIB,IIE
+ ! Loop over wind turbines
+ DO JROT=1,TFARM%NNB_TURBINES
+ ! X axis position test :
+ IF (TFARM%XPOS_X(JROT) >= XXHAT(JI) .AND. &
+ TFARM%XPOS_X(JROT) < XXHAT(JI) + PDXX(JI,JJ,JK)) THEN
+ ! YZ plane distances calculations
+ ZY_DIST = TFARM%XPOS_Y(JROT)-XYHAT(JJ)
+ ZZ_DIST = TTURBINE%XH_HEIGHT-(XZZ(JI,JJ,JK)-XZS(JI,JJ))
+ ZR_DIST = (ZY_DIST**2 + ZZ_DIST**2 )**(1./2.)
+ !
+ ! Disc position test
+ IF (ZR_DIST <= TTURBINE%XR_MAX) THEN
+!
+!* 2.2 Interpolating the wind
+!
+ ZPOS(1) = XXHAT(JI)
+ ZPOS(2) = XYHAT(JJ)
+ ZPOS(3) = XZZ(JI,JJ,JK)-XZS(JI,JJ)
+ SELECT CASE(CINTERP)
+ CASE('CLS')
+ ZUT_I = PUT_M(JI,JJ,JK)
+ ZRHO_I = PRHO_M(JI,JJ,JK)
+ CASE('8NB')
+ ZUT_I = INTERP_LIN8NB(ZPOS(:),JI,JJ,JK,PUT_M,ZZH)
+ ZRHO_I = INTERP_LIN8NB(ZPOS(:),JI,JJ,JK,PRHO_M,ZZH)
+ END SELECT
+!
+!* 2.3 Calculating the thrust of a cell wind->rotor
+!
+ PFX_RG(JI,JJ,JK) = PFX_RG(JI,JJ,JK) &
+ + 0.5*ZRHO_I*XCT_D(JROT) &
+ *PDYY(JI,JJ,JK)*PDZZ(JI,JJ,JK) &
+ *ZUT_I**2
+!
+!* 2.4 Calculating the thrust of the rotor wind->rotor
+! in a pseudo hub coordinate system (-)
+!
+ ZTHRUT(JROT) = ZTHRUT(JROT) &
+ - 0.5*ZRHO_I*XCT_D(JROT) &
+ *PDYY(JI,JJ,JK)*PDZZ(JI,JJ,JK) &
+ *ZUT_I**2
+!
+ END IF ! Disc position test
+ END IF ! X axis position test
+ END DO ! WT loop
+ END DO ! X domain loop
+ END DO ! Y domain loop
+END DO ! Z domain loop
+!
+!* 2.4 Bottom and top boundaries
+!
+PFX_RG(:,:,IKB-1) = PFX_RG(:,:,IKB)
+PFX_RG(:,:,IKE+1) = PFX_RG(:,:,IKE)
+!
+!* 3. SHARING THE DATAS OVER THE CPUS
+! -------------------------------
+!
+CALL MPI_ALLREDUCE(ZTHRUT, XTHRUT, SIZE(XTHRUT), &
+ MNHREAL_MPI,MPI_SUM,NMNH_COMM_WORLD,IINFO)
+!
+!
+!
+END SUBROUTINE EOL_ADNR
diff --git a/src/MNH/eol_alm.f90 b/src/MNH/eol_alm.f90
new file mode 100755
index 0000000000000000000000000000000000000000..585dd2d5c5ba8bb8e886aa9fc6b6a3afb017cd7b
--- /dev/null
+++ b/src/MNH/eol_alm.f90
@@ -0,0 +1,504 @@
+!MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
+!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
+!MNH_LIC for details. version 1.
+! #######################
+ MODULE MODI_EOL_ALM
+! #######################
+!
+INTERFACE
+!
+SUBROUTINE EOL_ALM(KTCOUNT, PTSTEP, &
+ PDXX, PDYY, PDZZ, &
+ PRHO_M, &
+ PUT_M, PVT_M, PWT_M, &
+ PFX_RG, PFY_RG, PFZ_RG )
+
+!
+INTEGER, INTENT(IN) :: KTCOUNT ! iteration count
+REAL, INTENT(IN) :: PTSTEP ! timestep except
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX,PDYY,PDZZ ! mesh size
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHO_M ! dry Density
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PUT_M,PVT_M,PWT_M ! Wind speed at mass point
+!
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PFX_RG ! Aerodynamic force ..
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PFY_RG ! .. cartesian mesh ..
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PFZ_RG ! .. global frame)
+!
+!
+END SUBROUTINE EOL_ALM
+!
+END INTERFACE
+!
+END MODULE MODI_EOL_ALM
+!
+! ###################################################################
+ SUBROUTINE EOL_ALM(KTCOUNT, PTSTEP, &
+ PDXX, PDYY, PDZZ, &
+ PRHO_M, &
+ PUT_M, PVT_M, PWT_M, &
+ PFX_RG, PFY_RG, PFZ_RG )
+! ###################################################################
+!
+!!**** *MODI_EOL_ALM* -
+!!
+!! PURPOSE
+!! -------
+!! It is possible to include wind turbines parameterization in Meso-NH,
+!! and several methods are available. One of the models is the Actuator
+!! Line Method (ALM). It allows to compute aerodynamic forces according
+!! the wind speed and the caracteristics of the wind turbine.
+!!
+!!** METHOD
+!! ------
+!! The ALM consists in modeling each blade by one line divided into blade
+!! element points (Sørensen and Shen, 2002). These points are applying
+!! aerodynamic forces into the flow.
+!! Each point carries a two-dimensional (2D) airfoil, and its characteris-
+!! tics, as lift and drag coefficients. Knowing these coefficients, and
+!! the angle of attack, the lift and drag forces can be evaluated.
+!!
+!! REFERENCE
+!! ---------
+!! PA. Joulin PhD Thesis. 2020.
+!!
+!!
+!! AUTHOR
+!! ------
+!! PA. Joulin *CNRM & IFPEN*
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 24/01/17
+!! Modification 14/10/20 (PA. Joulin) Updated for a main version
+!!
+!!---------------------------------------------------------------
+!
+!
+!* 0. DECLARATIONS
+! ------------
+!
+! To work with wind turbines
+USE MODD_EOL_ALM
+USE MODD_EOL_KINE_ALM
+!
+USE MODD_EOL_SHARED_IO, ONLY: CINTERP
+USE MODD_EOL_SHARED_IO, ONLY: XTHRUT, XTORQT, XPOWT
+!
+USE MODI_EOL_MATHS
+USE MODI_EOL_READER, ONLY: GET_AIRFOIL_ID
+USE MODI_EOL_PRINTER, ONLY: PRINT_TSPLIT
+USE MODI_EOL_ERROR, ONLY: EOL_WTCFL_ERROR
+! Math
+USE MODD_CST, ONLY: XPI
+! To know the grid
+USE MODD_GRID_n, ONLY: XXHAT,XYHAT,XZS,XZZ
+USE MODE_ll, ONLY: GET_INDICE_ll
+USE MODD_PARAMETERS, ONLY: JPVEXT
+! MPI stuffs
+USE MODD_VAR_ll, ONLY: NMNH_COMM_WORLD
+USE MODD_PRECISION, ONLY: MNHREAL_MPI
+USE MODD_MPIF, ONLY: MPI_SUM
+USE MODE_SUM_ll, ONLY: MIN_ll
+USE MODD_VAR_ll, ONLY: IP
+!
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of dummy arguments :
+!
+INTEGER, INTENT(IN) :: KTCOUNT ! iteration count
+REAL, INTENT(IN) :: PTSTEP ! timestep except
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX,PDYY,PDZZ ! mesh size
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHO_M ! dry Density
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PUT_M,PVT_M,PWT_M ! Wind speed at mass point
+!
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PFX_RG ! Aerodynamic force ..
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PFY_RG ! .. cartesian mesh ..
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PFZ_RG ! .. global frame)
+!
+!
+!* 0.2 Declarations of local variables :
+!
+! Indicies Compteurs
+INTEGER :: IIB,IJB,IKB ! Begin of a CPU domain
+INTEGER :: IIE,IJE,IKE ! End of a CPU domain
+INTEGER :: IKU ! Vertical size of the domain
+INTEGER :: JI, JJ, JK ! Loop index
+INTEGER :: JROT, JBLA, JBELT ! Rotor, blade, and blade element indicies
+!
+! Averages variables over all sub-timestep (if Time splitting)
+REAL, DIMENSION(TFARM%NNB_TURBINES,TTURBINE%NNB_BLADES,TBLADE%NNB_BLAELT) :: ZAOA_ATS ! Angle of attack of an element, hub frame [rad]
+REAL, DIMENSION(TFARM%NNB_TURBINES,TTURBINE%NNB_BLADES,TBLADE%NNB_BLAELT) :: ZFLIFT_ATS ! Aerodynamic lift force, parallel to Urel [N]
+REAL, DIMENSION(TFARM%NNB_TURBINES,TTURBINE%NNB_BLADES,TBLADE%NNB_BLAELT) :: ZFDRAG_ATS ! Aerodynamic drag force, perpendicular to Urel [N]
+REAL, DIMENSION(TFARM%NNB_TURBINES,TTURBINE%NNB_BLADES,TBLADE%NNB_BLAELT,3) :: ZFAERO_RE_ATS ! Aerodynamic force (lift+drag) in RE [N]
+REAL, DIMENSION(TFARM%NNB_TURBINES,TTURBINE%NNB_BLADES,TBLADE%NNB_BLAELT,3) :: ZFAERO_RG_ATS ! Aerodynamic force (lift+drag) in RG [N]
+
+
+! -- Wind --
+REAL :: ZRHO_I ! Interpolated density [kg/m3]
+REAL :: ZUT_I ! Interpolated wind speed U (RG) [m/s]
+REAL :: ZVT_I ! Interpolated wind speed V (RG) [m/s]
+REAL :: ZWT_I ! Interpolated wind speed W (RG) [m/s]
+REAL, DIMENSION(3) :: ZWIND_VEL_RG ! Wind velocity in RG frame [m/s]
+REAL, DIMENSION(3) :: ZWIND_VEL_RE ! Wind velocity in RE frame [m/s]
+REAL, DIMENSION(3) :: ZWINDREL_VEL_RE ! Relative wind velocity in RE frame [m/s]
+REAL :: ZWINDREL_VEL ! Norm of the relative wind velocity [m/s]
+REAL, DIMENSION(SIZE(PUT_M,1),SIZE(PUT_M,2),SIZE(PUT_M,3)) :: ZZH ! True heigth to interpolate 8NB
+!
+! -- Wind turbine --
+INTEGER :: INB_WT, INB_B, INB_BELT ! Total numbers
+REAL :: ZRAD ! Blade radius [m]
+INTEGER :: IAID ! Airfoil index [-]
+!
+! -- Aero --
+REAL :: ZAOA ! Attack angle of an element [rad]
+REAL :: ZCDRAG ! Drag coefficient of an element []
+REAL :: ZCLIFT ! Lift coefficient of an element []
+REAL :: ZFDRAG ! Drag force of an element, parallel to Urel [N]
+REAL :: ZFLIFT ! Lift force of an element, perpendicular to Urel [N]
+REAL, DIMENSION(3) :: ZFAERO_RE ! Aerodynamic force (lift+drag) in RE [N]
+REAL, DIMENSION(3) :: ZFAERO_RG ! Aerodynamic force (lift+drag) in RG [N]
+! Tip loss
+REAL :: ZFTIPL ! tip loss function
+REAL :: ZPHI ! angle twist+pitch+aa
+!
+! Thrust, Torque and Power
+REAL, DIMENSION(3) :: ZFAERO_RH ! Aerodynamic force (lift+drag) in RH [N] (thrust/torque)
+REAL, DIMENSION(3) :: ZDIST_HBELT_RH ! Distance between blade element and hub, in RH [m]
+REAL, DIMENSION(3) :: ZDIST_HBELT_RG ! Distances between blade element and hub, in RG [m]
+REAL, DIMENSION(3) :: Z3D_TORQT ! Full torque force (3D) of the wind turbine [N]
+!
+! -- Time spliting --
+INTEGER :: KTSUBCOUNT, INBSUBCOUNT ! sub iteration count
+REAL :: ZTSUBSTEP ! sub timestep
+REAL :: ZMAXTSTEP ! Max value for timestep to respect WTCFL criteria
+!
+! -- Numerical --
+INTEGER :: IINFO ! code info return
+!
+!
+!* 0.3 Implicit arguments
+!
+! A. From MODD_EOL_ALM
+!TYPE(FARM) :: TFARM
+!TYPE(TURBINE) :: TTURBINE
+!TYPE(BLADE) :: TBLADE
+!TYPE(AIRFOIL), DIMENSION(:), ALLOCATABLE :: TAIRFOIL
+!
+!REAL, DIMENSION(:,:,:), ALLOCATABLE :: XELT_RAD ! Blade elements radius [m]
+!REAL, DIMENSION(:,:,:), ALLOCATABLE :: XAOA_GLB ! Angle of attack of an element [rad]
+!REAL, DIMENSION(:,:,:), ALLOCATABLE :: XFLIFT_GLB ! Lift force, parallel to Urel [N]
+!REAL, DIMENSION(:,:,:), ALLOCATABLE :: XFDRAG_GLB ! Drag force, perpendicular to Urel [N]
+!REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: XFAERO_RE_GLB ! Aerodyn. force (lift+drag) in RE [N]
+!REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: XFAERO_RG_GLB ! Aerodyn. force (lift+drag) in RG [N]
+!
+!INTEGER :: NNB_BLAELT ! Number of blade elements
+!LOGICAL :: LTIMESPLIT ! Flag to apply Time splitting method
+!LOGICAL :: LTIPLOSSG ! Flag to apply Glauert's tip loss correction
+!LOGICAL :: LTECOUTPTS ! Flag to get Tecplot file output of element points
+!
+! B. From MODD_EOL_SHARED_IO:
+! for namelist NAM_EOL_ALM
+!CHARACTER(LEN=100) :: CFARM_CSVDATA ! Farm file to read
+!CHARACTER(LEN=100) :: CTURBINE_CSVDATA ! Turbine file to read
+!CHARACTER(LEN=100) :: CBLADE_CSVDATA ! Blade file to read
+!CHARACTER(LEN=100) :: CAIRFOIL_CSVDATA ! Airfoil file to read
+!CHARACTER(LEN=3) :: CINTERP ! Interpolation method for wind speed
+! for output
+!REAL, DIMENSION(:), ALLOCATABLE :: XTHRUT ! Thrust [N]
+!REAL, DIMENSION(:), ALLOCATABLE :: XTORQT ! Torque [Nm]
+!REAL, DIMENSION(:), ALLOCATABLE :: XPOWT ! Power [W]
+!
+!
+!-------------------------------------------------------------------------------
+!
+!
+!* 1. INITIALIZATIONS
+! ---------------
+!
+!* 1.1 Subdomain (CPU) indices
+!
+CALL GET_INDICE_ll(IIB,IJB,IIE,IJE) ! Get begin and end domain index (CPU)
+IKU = SIZE(PUT_M,3) ! Top of the domain end index
+IKB=1+JPVEXT ! Vertical begin index
+IKE=IKU-JPVEXT ! Vertical end index
+!
+!* 1.2 Some usefull integers
+!
+INB_WT = TFARM%NNB_TURBINES
+INB_B = TTURBINE%NNB_BLADES
+INB_BELT = TBLADE%NNB_BLAELT
+!
+!* 1.3 Vertical coordinate in case of interpolation
+!
+IF (CINTERP=='8NB') THEN
+ DO JK=1,IKU-1
+ ZZH(:,:,JK) = (0.5*(XZZ(:,:,JK)+XZZ(:,:,JK+1))-XZS(:,:))
+ END DO
+ ZZH(:,:,IKU) = 2*ZZH(:,:,IKU-1) - ZZH(:,:,IKU-2)
+END IF
+!
+!* 1.4 Set to zeros at each MNH time steps
+!
+! Averaged variables (over time splitting)
+ZAOA_ATS(:,:,:) = 0.
+ZFLIFT_ATS(:,:,:) = 0.
+ZFDRAG_ATS(:,:,:) = 0.
+ZFAERO_RE_ATS(:,:,:,:) = 0.
+ZFAERO_RG_ATS(:,:,:,:) = 0.
+!
+! Global variables (seen by all CPU)
+XAOA_GLB(:,:,:) = 0.
+XFLIFT_GLB(:,:,:) = 0.
+XFDRAG_GLB(:,:,:) = 0.
+XFAERO_RE_GLB(:,:,:,:) = 0.
+XFAERO_RG_GLB(:,:,:,:) = 0.
+!
+XTHRUT(:) = 0.
+XTORQT(:) = 0.
+!
+!
+!-------------------------------------------------------------------------------
+!
+!* 2. COMPUTES WTCFL CRITERIA
+! -----------------------
+!
+!* 2.1 Computing the highest timestep acceptable
+ZMAXTSTEP = ABS( MIN(MIN_ll(PDXX(:,:,:),IINFO),&
+ MIN_ll(PDYY(:,:,:),IINFO),&
+ MIN_ll(PDZZ(:,:,:),IINFO))&
+ /(MAXVAL(TFARM%XOMEGA(:))*TTURBINE%XR_MAX))
+!
+IF (.NOT.LTIMESPLIT) THEN
+!* 2.2 Checking conditions
+! If time step too high : abort
+ IF (PTSTEP > ZMAXTSTEP) THEN
+ CALL EOL_WTCFL_ERROR(ZMAXTSTEP)
+ STOP
+! If time step ok, continue
+ ELSE
+ INBSUBCOUNT = 1
+ ZTSUBSTEP = PTSTEP/INBSUBCOUNT
+ END IF
+ELSE
+!* 2.3 Timesplitting : new sub-timestep
+ INBSUBCOUNT = INT(PTSTEP/ZMAXTSTEP) + 1
+ ZTSUBSTEP = PTSTEP/INBSUBCOUNT
+ CALL PRINT_TSPLIT(INBSUBCOUNT, ZTSUBSTEP)
+END IF
+!
+!* 2.4 Start looping over sub-timesteps
+DO KTSUBCOUNT=1,INBSUBCOUNT
+!
+!
+!-------------------------------------------------------------------------------
+!
+!* 3. KINEMATICS COMPUTATIONS
+! -----------------------
+!
+ CALL EOL_KINE_ALM(KTCOUNT, KTSUBCOUNT, ZTSUBSTEP, PTSTEP)
+!
+!
+!-------------------------------------------------------------------------------
+!
+!* 4. COMPUTES AERODYNAMIC FORCES THAT ACTS ON THE BLADES DUE TO THE WIND
+! --------------------------------------------------------------
+!
+!* 4.1 Finding the position of wind turbines
+!
+! Loop over domain
+ DO JK=IKB,IKE
+ DO JJ=IJB,IJE
+ DO JI=IIB,IIE
+ ! Loop over wind turbines
+ DO JROT=1, INB_WT
+ DO JBLA=1, INB_B
+ DO JBELT=1, INB_BELT
+ ! Position test
+ IF (XPOS_ELT_RG(JROT,JBLA,JBELT,1) >= XXHAT(JI) .AND. &
+ XPOS_ELT_RG(JROT,JBLA,JBELT,1) < XXHAT(JI) + PDXX(JI,JJ,JK)) THEN
+!
+ IF (XPOS_ELT_RG(JROT,JBLA,JBELT,2) >= XYHAT(JJ) .AND. &
+ XPOS_ELT_RG(JROT,JBLA,JBELT,2) < XYHAT(JJ) + PDYY(JI,JJ,JK)) THEN
+!
+ IF (XPOS_ELT_RG(JROT,JBLA,JBELT,3) >= XZZ(JI,JJ,JK) .AND. &
+ XPOS_ELT_RG(JROT,JBLA,JBELT,3) < XZZ(JI,JJ,JK) + PDZZ(JI,JJ,JK)) THEN
+!
+!* 4.2 Extracting the wind
+!
+ SELECT CASE(CINTERP)
+ CASE('CLS')
+ ZUT_I = PUT_M(JI,JJ,JK)
+ ZVT_I = PVT_M(JI,JJ,JK)
+ ZWT_I = PWT_M(JI,JJ,JK)
+ ZRHO_I = PRHO_M(JI,JJ,JK)
+ CASE('8NB')
+ ZUT_I = INTERP_LIN8NB(XPOS_ELT_RG(JROT,JBLA,JBELT,:),&
+ JI,JJ,JK,PUT_M,ZZH)
+ ZVT_I = INTERP_LIN8NB(XPOS_ELT_RG(JROT,JBLA,JBELT,:),&
+ JI,JJ,JK,PVT_M,ZZH)
+ ZWT_I = INTERP_LIN8NB(XPOS_ELT_RG(JROT,JBLA,JBELT,:),&
+ JI,JJ,JK,PWT_M,ZZH)
+ ZRHO_I = INTERP_LIN8NB(XPOS_ELT_RG(JROT,JBLA,JBELT,:),&
+ JI,JJ,JK,PRHO_M,ZZH)
+ END SELECT
+ ZWIND_VEL_RG(1) = ZUT_I
+ ZWIND_VEL_RG(2) = ZVT_I
+ ZWIND_VEL_RG(3) = ZWT_I
+!
+!* 4.3 Calculating the wind in RE frame
+!
+ ZWIND_VEL_RE(:) = MATMUL(XMAT_RE_RG(JROT,JBLA,JBELT,:,:), ZWIND_VEL_RG(:))
+!
+!* 4.4 Calculating the relative wind speed in RE frame + norm
+!
+ ZWINDREL_VEL_RE(:) = ZWIND_VEL_RE(:) - XTVEL_ELT_RE(JROT,JBLA,JBELT,:)
+ ZWINDREL_VEL = NORM(ZWINDREL_VEL_RE)
+!
+!* 4.5 Calculating the angle of attack
+!
+ ZAOA = ATAN2(ZWINDREL_VEL_RE(1), ZWINDREL_VEL_RE(2))
+!
+!* 4.6 Getting aerodynamic coefficients from tabulated data
+!
+ ZRAD = XELT_RAD(JROT,JBLA,JBELT) ! Radius of the element
+ IAID = GET_AIRFOIL_ID(TTURBINE,TBLADE,TAIRFOIL,ZRAD) ! ID of the airfoil
+ ZCLIFT = INTERP_SPLCUB(ZAOA*180/XPI, &
+ TAIRFOIL(IAID)%XAA,&
+ TAIRFOIL(IAID)%XCL)
+ ZCDRAG = INTERP_SPLCUB(ZAOA*180/XPI, &
+ TAIRFOIL(IAID)%XAA,&
+ TAIRFOIL(IAID)%XCD)
+!
+!* 4.7 Tip loss correction (Glauert)
+!
+ IF (LTIPLOSSG) THEN
+ ZPHI = + ZAOA &
+ + TFARM%XBLA_PITCH(JROT) &
+ + XTWIST_ELT(JROT,JBLA,JBELT)
+ IF (ZPHI > 0.0) THEN
+ ZFTIPL = (2.0/XPI)*ACOS(MIN( &
+ 1.0, EXP(-(TTURBINE%NNB_BLADES/2.0) &
+ *(TTURBINE%XR_MAX-ZRAD)/(ZRAD*SIN(ZPHI)))))
+ ELSE
+ ZFTIPL = 1.0
+ END IF
+ ZCLIFT = ZFTIPL*ZCLIFT
+ ZCDRAG = ZFTIPL*ZCDRAG
+ END IF
+!
+!* 4.8 Computing aerodynamic forces in relative frame
+! that act on blades (wind->blade)
+ ZFLIFT = 0.5*ZRHO_I*XSURF_ELT(JROT,JBLA,JBELT)*ZCLIFT*ZWINDREL_VEL**2
+ ZFDRAG = 0.5*ZRHO_I*XSURF_ELT(JROT,JBLA,JBELT)*ZCDRAG*ZWINDREL_VEL**2
+!
+!* 4.9 Evaluating the aerodynamiques forces in RE frame
+! that act on blades (wind->blade)
+ ZFAERO_RE(1) = SIN(ZAOA)*ZFDRAG + COS(ZAOA)*ZFLIFT
+ ZFAERO_RE(2) = COS(ZAOA)*ZFDRAG - SIN(ZAOA)*ZFLIFT
+ ZFAERO_RE(3) = .0 ! 2D flow around arifoil assumption
+!
+!* 4.10 Evaluating the aerodynamiques forces in RG frame
+! that act on blades (wind->blade)
+ ZFAERO_RG(:) = MATMUL(XMAT_RG_RE(JROT,JBLA,JBELT,:,:), ZFAERO_RE(:))
+!
+!* 4.11 Adding it to the cell of Meso-NH
+ PFX_RG(JI,JJ,JK) = PFX_RG(JI,JJ,JK) + ZFAERO_RG(1) / FLOAT(INBSUBCOUNT)
+ PFY_RG(JI,JJ,JK) = PFY_RG(JI,JJ,JK) + ZFAERO_RG(2) / FLOAT(INBSUBCOUNT)
+ PFZ_RG(JI,JJ,JK) = PFZ_RG(JI,JJ,JK) + ZFAERO_RG(3) / FLOAT(INBSUBCOUNT)
+!
+!* 4.12 Storing mean values over one full MNH timestep
+! (all the sub-timesteps values are averaged)
+ ZAOA_ATS(JROT,JBLA,JBELT) = ZAOA_ATS(JROT,JBLA,JBELT) &
+ + ZAOA / FLOAT(INBSUBCOUNT)
+ ZFLIFT_ATS(JROT,JBLA,JBELT) = ZFLIFT_ATS(JROT,JBLA,JBELT) &
+ + ZFLIFT / FLOAT(INBSUBCOUNT)
+ ZFDRAG_ATS(JROT,JBLA,JBELT) = ZFDRAG_ATS(JROT,JBLA,JBELT) &
+ + ZFDRAG / FLOAT(INBSUBCOUNT)
+ ZFAERO_RE_ATS(JROT,JBLA,JBELT,:)= ZFAERO_RE_ATS(JROT,JBLA,JBELT,:) &
+ + ZFAERO_RE(:) / FLOAT(INBSUBCOUNT)
+ ZFAERO_RG_ATS(JROT,JBLA,JBELT,:)= ZFAERO_RG_ATS(JROT,JBLA,JBELT,:) &
+ + ZFAERO_RG(:) / FLOAT(INBSUBCOUNT)
+!
+ ! End of position tests
+ END IF
+ END IF
+ END IF
+ ! End of wind turbine loops
+ END DO
+ END DO
+ END DO
+ ! End of domain loops
+ END DO
+ END DO
+ END DO
+! End of sub-time loop
+END DO
+!
+!
+!* 4.13 Top and bottom conditions
+PFX_RG(:,:,IKB-1) = PFX_RG(:,:,IKB)
+PFX_RG(:,:,IKE+1) = PFX_RG(:,:,IKE)
+!
+PFY_RG(:,:,IKB-1) = PFY_RG(:,:,IKB)
+PFY_RG(:,:,IKE+1) = PFY_RG(:,:,IKE)
+!
+PFZ_RG(:,:,IKB-1) = PFZ_RG(:,:,IKB)
+PFZ_RG(:,:,IKE+1) = PFZ_RG(:,:,IKE)
+!
+!
+!-------------------------------------------------------------------------------
+!
+!* 5. SHARING THE DATAS OVER THE CPUS
+! -------------------------------
+CALL MPI_ALLREDUCE(ZAOA_ATS, XAOA_GLB, SIZE(XAOA_GLB), &
+ MNHREAL_MPI,MPI_SUM,NMNH_COMM_WORLD,IINFO)
+CALL MPI_ALLREDUCE(ZFLIFT_ATS, XFLIFT_GLB, SIZE(XFLIFT_GLB), &
+ MNHREAL_MPI,MPI_SUM,NMNH_COMM_WORLD,IINFO)
+CALL MPI_ALLREDUCE(ZFDRAG_ATS, XFDRAG_GLB, SIZE(XFDRAG_GLB), &
+ MNHREAL_MPI,MPI_SUM,NMNH_COMM_WORLD,IINFO)
+CALL MPI_ALLREDUCE(ZFAERO_RE_ATS, XFAERO_RE_GLB, SIZE(XFAERO_RE_GLB),&
+ MNHREAL_MPI,MPI_SUM,NMNH_COMM_WORLD,IINFO)
+CALL MPI_ALLREDUCE(ZFAERO_RG_ATS, XFAERO_RG_GLB, SIZE(XFAERO_RG_GLB),&
+ MNHREAL_MPI,MPI_SUM,NMNH_COMM_WORLD,IINFO)
+!
+!
+!-------------------------------------------------------------------------------
+!
+!* 6. COMPUTING THRUST, TORQUE AND POWER
+! ---------------------------------
+!
+IF(IP == 1) THEN
+ DO JROT=1,TFARM%NNB_TURBINES
+ DO JBLA=1, TTURBINE%NNB_BLADES
+ DO JBELT=1, TBLADE%NNB_BLAELT
+!
+!* 6.1 Preliminaries
+! Aerodynamic load (wind->blade) in RH
+ ZFAERO_RH(:) = MATMUL(XMAT_RH_RG(JROT,:,:), &
+ XFAERO_RG_GLB(JROT,JBLA,JBELT,:))
+! Distance between element and hub in RG
+ ZDIST_HBELT_RG(:) = XPOS_ELT_RG(JROT,JBLA,JBELT,:) - XPOS_HUB_RG(JROT,:)
+! Distance between element and hub in RH
+ ZDIST_HBELT_RH(:) = MATMUL(XMAT_RH_RG(JROT,:,:),ZDIST_HBELT_RG(:))
+!
+!* 6.2 Thrust (wind->rotor): in RH
+ XTHRUT(JROT) = XTHRUT(JROT) + ZFAERO_RH(3) ! Only Z component
+!* 6.3 Torque (wind->rotor) in RH
+ Z3D_TORQT = CROSS(ZDIST_HBELT_RH(:),ZFAERO_RH(:))
+ XTORQT(JROT) = XTORQT(JROT) + Z3D_TORQT(3) ! Only Z component
+ END DO
+ END DO
+!
+!* 6.4 Power (wind->rotor)
+ XPOWT(JROT) = XTORQT(JROT) * TFARM%XOMEGA(JROT)
+ END DO
+END IF
+!
+!
+END SUBROUTINE EOL_ALM
diff --git a/src/MNH/eol_debugger.f90 b/src/MNH/eol_debugger.f90
new file mode 100755
index 0000000000000000000000000000000000000000..0f9f21802f65b9737e5dc6c8d8286a17fb4c0424
--- /dev/null
+++ b/src/MNH/eol_debugger.f90
@@ -0,0 +1,239 @@
+!MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
+!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
+!MNH_LIC for details. version 1.
+! #######################
+ MODULE MODI_EOL_DEBUGGER
+! #######################
+!
+INTERFACE
+!
+!
+SUBROUTINE PRINTMER_ll(HNAME,PVAR)
+ CHARACTER(LEN=*), INTENT(IN) :: HNAME ! Name of the variable,
+ REAL, INTENT(IN) :: PVAR ! Variable,
+END SUBROUTINE PRINTMER_ll
+!
+SUBROUTINE PRINTMEI_ll(HNAME,KVAR)
+ CHARACTER(LEN=*), INTENT(IN) :: HNAME ! Name of the variable,
+ INTEGER, INTENT(IN) :: KVAR ! Variable,
+END SUBROUTINE PRINTMEI_ll
+!
+SUBROUTINE PRINTMEC_ll(HNAME,CVAR)
+ CHARACTER(LEN=*), INTENT(IN) :: HNAME ! Name of the variable,
+ CHARACTER(LEN=*), INTENT(IN) :: CVAR ! Variable,
+END SUBROUTINE PRINTMEC_ll
+!
+SUBROUTINE PRINTMER_CPU1(HNAME,PVAR)
+ CHARACTER(LEN=*), INTENT(IN) :: HNAME ! Name of the variable,
+ REAL, INTENT(IN) :: PVAR ! Variable,
+END SUBROUTINE PRINTMER_CPU1
+!
+SUBROUTINE PRINTMEI_CPU1(HNAME,KVAR)
+ CHARACTER(LEN=*), INTENT(IN) :: HNAME ! Name of the variable,
+ INTEGER, INTENT(IN) :: KVAR ! Variable,
+END SUBROUTINE PRINTMEI_CPU1
+!
+SUBROUTINE PRINTMEC_CPU1(HNAME,CVAR)
+ CHARACTER(LEN=*), INTENT(IN) :: HNAME ! Name of the variable,
+ CHARACTER(LEN=*), INTENT(IN) :: CVAR ! Variable,
+END SUBROUTINE PRINTMEC_CPU1
+!
+SUBROUTINE PRINTMER_ELT1(KROT,KBLA,KELT,HNAME,PVAR)
+ INTEGER, INTENT(IN) :: KROT, KBLA, KELT ! Loop index,
+ CHARACTER(LEN=*), INTENT(IN) :: HNAME ! Name of the variable,
+ REAL, INTENT(IN) :: PVAR ! Variable,
+END SUBROUTINE PRINTMER_ELT1
+!
+SUBROUTINE PRINTMER_ELT42(KROT,KBLA,KELT,HNAME,PVAR)
+ INTEGER, INTENT(IN) :: KROT, KBLA, KELT ! Loop index,
+ CHARACTER(LEN=*), INTENT(IN) :: HNAME ! Name of the variable,
+ REAL, INTENT(IN) :: PVAR ! Variable,
+END SUBROUTINE PRINTMER_ELT42
+!
+SUBROUTINE PRINTMER_3BELT42(KROT,KBLA,KELT,HNAME,PVAR)
+ INTEGER, INTENT(IN) :: KROT, KBLA, KELT ! Loop index,
+ CHARACTER(LEN=*), INTENT(IN) :: HNAME ! Name of the variable,
+ REAL, INTENT(IN) :: PVAR ! Variable,
+END SUBROUTINE PRINTMER_3BELT42
+!
+END INTERFACE
+!
+END MODULE MODI_EOL_DEBUGGER
+!-------------------------------------------------------------------
+!
+!!**** *EOL_PRINTER* -
+!!
+!! PURPOSE
+!! -------
+!! Some usefull toold to debbug my code
+!!
+!! AUTHOR
+!! ------
+!! PA. Joulin *CNRM & IFPEN*
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 26/10/2020
+!!
+!!---------------------------------------------------------------
+!
+!#########################################################
+SUBROUTINE PRINTMER_ll(HNAME,PVAR)
+!
+USE MODD_VAR_ll, ONLY: IP
+!
+IMPLICIT NONE
+!
+CHARACTER(LEN=*), INTENT(IN) :: HNAME ! Name of the variable,
+REAL, INTENT(IN) :: PVAR ! Variable,
+!
+!
+PRINT*, 'CPU n. ', IP, ' : ', HNAME, ' = ', PVAR
+!
+END SUBROUTINE PRINTMER_ll
+!#########################################################
+!
+!#########################################################
+SUBROUTINE PRINTMEI_ll(HNAME,KVAR)
+!
+USE MODD_VAR_ll, ONLY: IP
+!
+IMPLICIT NONE
+!
+CHARACTER(LEN=*), INTENT(IN) :: HNAME ! Name of the variable,
+INTEGER, INTENT(IN) :: KVAR ! Variable,
+!
+!
+PRINT*, 'CPU n. ', IP, ' : ', HNAME, ' = ', KVAR
+!
+END SUBROUTINE PRINTMEI_ll
+!#########################################################
+!
+!#########################################################
+SUBROUTINE PRINTMEC_ll(HNAME,CVAR)
+!
+USE MODD_VAR_ll, ONLY: IP
+!
+IMPLICIT NONE
+!
+CHARACTER(LEN=*), INTENT(IN) :: HNAME ! Name of the variable,
+CHARACTER(LEN=*), INTENT(IN) :: CVAR ! Variable,
+!
+!
+PRINT*, 'CPU n. ', IP, ' : ', HNAME, ' = ', CVAR
+!
+END SUBROUTINE PRINTMEC_ll
+!#########################################################
+!
+!#########################################################
+SUBROUTINE PRINTMER_CPU1(HNAME, PVAR)
+!
+USE MODD_VAR_ll, ONLY: IP
+!
+IMPLICIT NONE
+!
+CHARACTER(LEN=*), INTENT(IN) :: HNAME ! Name of the variable,
+REAL, INTENT(IN) :: PVAR ! Variable,
+!
+IF (IP==1) THEN
+ PRINT*, HNAME, ' = ', PVAR
+END IF
+!
+END SUBROUTINE PRINTMER_CPU1
+!#########################################################
+!
+!#########################################################
+SUBROUTINE PRINTMEI_CPU1(HNAME, KVAR)
+!
+USE MODD_VAR_ll, ONLY: IP
+!
+IMPLICIT NONE
+!
+CHARACTER(LEN=*), INTENT(IN) :: HNAME ! Name of the variable,
+INTEGER, INTENT(IN) :: KVAR ! Variable,
+!
+IF (IP==1) THEN
+ PRINT*, HNAME, ' = ', KVAR
+END IF
+!
+END SUBROUTINE PRINTMEI_CPU1
+!#########################################################
+!
+!#########################################################
+SUBROUTINE PRINTMEC_CPU1(HNAME, CVAR)
+!
+USE MODD_VAR_ll, ONLY: IP
+!
+IMPLICIT NONE
+!
+CHARACTER(LEN=*), INTENT(IN) :: HNAME ! Name of the variable,
+CHARACTER(LEN=*), INTENT(IN) :: CVAR ! Variable,
+!
+IF (IP==1) THEN
+ PRINT*, HNAME, ' = ', CVAR
+END IF
+!
+END SUBROUTINE PRINTMEC_CPU1
+!#########################################################
+!
+!#########################################################
+SUBROUTINE PRINTMER_ELT1(KROT, KBLA, KELT, HNAME, PVAR)
+!
+USE MODD_CST, ONLY: XPI
+!
+IMPLICIT NONE
+!
+INTEGER, INTENT(IN) :: KROT, KBLA, KELT ! Loop index,
+CHARACTER(LEN=*), INTENT(IN) :: HNAME ! Name of the variable,
+REAL, INTENT(IN) :: PVAR ! Variable,
+!
+IF ((KROT==1).AND.(KBLA==1).AND.(KELT==1)) THEN
+ PRINT*, HNAME, PVAR
+END IF
+!
+END SUBROUTINE PRINTMER_ELT1
+!#########################################################
+!
+!#########################################################
+SUBROUTINE PRINTMER_ELT42(KROT, KBLA, KELT, HNAME, PVAR)
+!
+USE MODD_CST, ONLY: XPI
+!
+IMPLICIT NONE
+!
+INTEGER, INTENT(IN) :: KROT, KBLA, KELT ! Loop index,
+CHARACTER(LEN=*), INTENT(IN) :: HNAME ! Name of the variable,
+REAL, INTENT(IN) :: PVAR ! Variable,
+!
+IF ((KROT==1).AND.(KBLA==1).AND.(KELT==42)) THEN
+ PRINT*, HNAME, PVAR
+END IF
+!
+END SUBROUTINE PRINTMER_ELT42
+!#########################################################
+!
+!#########################################################
+SUBROUTINE PRINTMER_3BELT42(KROT, KBLA, KELT, HNAME, PVAR)
+!
+USE MODD_CST, ONLY: XPI
+!
+IMPLICIT NONE
+!
+INTEGER, INTENT(IN) :: KROT, KBLA, KELT ! Loop index,
+CHARACTER(LEN=*), INTENT(IN) :: HNAME ! Name of the variable,
+REAL, INTENT(IN) :: PVAR ! Variable,
+!
+IF ((KROT==1).AND.(KBLA==1).AND.(KELT==42)) THEN
+ PRINT*, HNAME, 'B1 = ', PVAR*180/XPI
+END IF
+IF ((KROT==1).AND.(KBLA==2).AND.(KELT==42)) THEN
+ PRINT*, HNAME, 'B2 = ', PVAR*180/XPI
+END IF
+IF ((KROT==1).AND.(KBLA==3).AND.(KELT==42)) THEN
+ PRINT*, HNAME, 'B3 = ', PVAR*180/XPI
+END IF
+!
+END SUBROUTINE PRINTMER_3BELT42
+!#########################################################
+!
diff --git a/src/MNH/eol_error.f90 b/src/MNH/eol_error.f90
new file mode 100755
index 0000000000000000000000000000000000000000..06f8f6fae83d2b3510d2da94242caabf2d009622
--- /dev/null
+++ b/src/MNH/eol_error.f90
@@ -0,0 +1,192 @@
+!MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
+!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
+!MNH_LIC for details. version 1.
+! #######################
+ MODULE MODI_EOL_ERROR
+! #######################
+!
+INTERFACE
+!
+! **********
+! EOL_READER
+! **********
+!
+SUBROUTINE EOL_CSVNOTFOUND_ERROR(HFILE)
+ CHARACTER(LEN=*), INTENT(IN) :: HFILE ! file read
+END SUBROUTINE EOL_CSVNOTFOUND_ERROR
+!
+SUBROUTINE EOL_CSVEMPTY_ERROR(HFILE,KNBLINE)
+ CHARACTER(LEN=*), INTENT(IN) :: HFILE ! file read
+ INTEGER, INTENT(IN) :: KNBLINE ! number of lines
+END SUBROUTINE EOL_CSVEMPTY_ERROR
+!
+!
+! ***
+! ALM
+! ***
+!
+SUBROUTINE EOL_AIRFOILNOTFOUND_ERROR(HFILE,HVAR)
+ CHARACTER(LEN=*), INTENT(IN) :: HFILE ! file read
+ CHARACTER(LEN=*), INTENT(IN) :: HVAR ! missing data
+END SUBROUTINE EOL_AIRFOILNOTFOUND_ERROR
+!
+SUBROUTINE EOL_WTCFL_ERROR(PMAXTSTEP)
+ REAL, INTENT(IN) :: PMAXTSTEP ! maximum acceptable time-step
+END SUBROUTINE EOL_WTCFL_ERROR
+!
+SUBROUTINE EOL_BLADEDATA_ERROR(PDELTARAD)
+ REAL, INTENT(IN) :: PDELTARAD ! Span lenght of an element
+END SUBROUTINE EOL_BLADEDATA_ERROR
+!
+END INTERFACE
+!
+END MODULE MODI_EOL_ERROR
+!-------------------------------------------------------------------
+!
+!!**** *EOL_ERROR* -
+!!
+!! PURPOSE
+!! -------
+!! Some usefull subs to manage errors linked to wind turbines
+!!
+!! AUTHOR
+!! ------
+!! PA. Joulin *CNRM & IFPEN*
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 15/12/2020
+!!
+!!---------------------------------------------------------------
+!
+!#########################################################
+SUBROUTINE EOL_CSVNOTFOUND_ERROR(HFILE)
+!
+USE MODD_LUNIT_n, ONLY: TLUOUT
+USE MODE_IO_FILE, ONLY: IO_File_close
+USE MODE_MSG
+!
+CHARACTER(LEN=*), INTENT(IN) :: HFILE ! file read
+!
+WRITE(TLUOUT%NLU,*) ''
+WRITE(TLUOUT%NLU,'(A)') &
+ ' EOL Initialization error. Program aborted.'
+WRITE(TLUOUT%NLU,'(A,A,A)') &
+ ' File "', TRIM(HFILE), '" not found in directory.'
+WRITE(TLUOUT%NLU,*) ''
+!
+CALL IO_File_close(TLUOUT)
+CALL PRINT_MSG( NVERB_FATAL, 'GEN', 'EOL_CSVNOTFOUND_ERROR', 'CSV file for wind turbine missing' )
+!
+END SUBROUTINE EOL_CSVNOTFOUND_ERROR
+!#########################################################
+
+!
+!#########################################################
+SUBROUTINE EOL_CSVEMPTY_ERROR(HFILE,KNBLINE)
+!
+USE MODD_LUNIT_n, ONLY: TLUOUT
+USE MODE_IO_FILE, ONLY: IO_File_close
+USE MODE_MSG
+!
+CHARACTER(LEN=*), INTENT(IN) :: HFILE ! file read
+INTEGER, INTENT(IN) :: KNBLINE ! number of lines
+!
+WRITE(TLUOUT%NLU,*) ''
+WRITE(TLUOUT%NLU,'(A)') &
+ ' EOL Initialization error. Program aborted.'
+WRITE(TLUOUT%NLU,'(I2,A,A,A)') &
+ KNBLINE, ' line have been read in file: ', TRIM(HFILE), '.'
+WRITE(TLUOUT%NLU,'(A)') &
+ ' At least 2 should be there: header + data.'
+WRITE(TLUOUT%NLU,*) ''
+!
+CALL IO_File_close(TLUOUT)
+CALL PRINT_MSG( NVERB_FATAL, 'GEN', 'EOL_CSVEMPTY_ERROR', 'CSV file for wind turbine missing data' )
+!
+END SUBROUTINE EOL_CSVEMPTY_ERROR
+!#########################################################
+!
+!#########################################################
+SUBROUTINE EOL_AIRFOILNOTFOUND_ERROR(HFILE,HVAR)
+!
+USE MODD_LUNIT_n, ONLY: TLUOUT
+USE MODE_IO_FILE, ONLY: IO_File_close
+USE MODE_MSG
+!
+CHARACTER(LEN=*), INTENT(IN) :: HFILE ! file read
+CHARACTER(LEN=*), INTENT(IN) :: HVAR ! missing data
+!
+WRITE(TLUOUT%NLU,*) ''
+WRITE(TLUOUT%NLU,'(A)') &
+ ' EOL Initialization error. Program aborted.'
+WRITE(TLUOUT%NLU,'(A,A,A)') &
+ ' I am looking for the characteristics of ', TRIM(HVAR), '.'
+WRITE(TLUOUT%NLU,'(A,A)') &
+ ' But I cannot find them in file: ', TRIM(HFILE)
+WRITE(TLUOUT%NLU,*) ''
+!
+CALL IO_File_close(TLUOUT)
+CALL PRINT_MSG( NVERB_FATAL, 'GEN', 'EOL_AIRFOILNOTFOUND_ERROR', 'File for Airfol missing' )
+!
+END SUBROUTINE EOL_AIRFOILNOTFOUND_ERROR
+!#########################################################
+!
+!#########################################################
+SUBROUTINE EOL_WTCFL_ERROR(PMAXTSTEP)
+!
+USE MODD_LUNIT_n, ONLY: TLUOUT
+USE MODE_IO_FILE, ONLY: IO_File_close
+USE MODE_MSG
+!
+REAL, INTENT(IN) :: PMAXTSTEP ! maximum acceptable time-step
+!
+WRITE(TLUOUT%NLU,*) ''
+WRITE(TLUOUT%NLU,'(A)') &
+ 'Sorry but I had to stop the simulation. '
+WRITE(TLUOUT%NLU,'(A)') &
+ 'The time step XTSTEP is too large: '
+WRITE(TLUOUT%NLU,'(A)') &
+ 'the blades can jump over one or several cells. '
+WRITE(TLUOUT%NLU,'(A)') &
+ 'Please, turn on the time-splitting method (LTIMESPLIT=.TRUE.), '
+WRITE(TLUOUT%NLU,'(A,F10.8,A)') &
+ 'or decrease XTSTEP to a value lower than ', PMAXTSTEP, ' sec.'
+WRITE(TLUOUT%NLU,*) ''
+!
+CALL IO_File_close(TLUOUT)
+CALL PRINT_MSG( NVERB_FATAL, 'GEN', 'EOL_WTCFL_ERROR', 'WRONG TIME-STEP with wind turbine' )
+!
+END SUBROUTINE EOL_WTCFL_ERROR
+!#########################################################
+!
+!#########################################################
+SUBROUTINE EOL_BLADEDATA_ERROR(PDELTARAD)
+!
+USE MODD_LUNIT_n, ONLY: TLUOUT
+USE MODE_IO_FILE, ONLY: IO_File_close
+USE MODE_MSG
+USE MODD_EOL_SHARED_IO, ONLY: CBLADE_CSVDATA
+!
+REAL, INTENT(IN) :: PDELTARAD ! hals section width
+!
+WRITE(TLUOUT%NLU,*) ''
+WRITE(TLUOUT%NLU,'(A)') &
+ ' EOL Initialization error. Program aborted.'
+WRITE(TLUOUT%NLU,'(A,F4.2,A)') &
+ 'A blade element center position is set to ', PDELTARAD, '.'
+WRITE(TLUOUT%NLU,'(A)') &
+ 'As a blade element center, it has to be set in ]0%;100%[. '
+WRITE(TLUOUT%NLU,'(A,A,A)') &
+ 'Please, check your blade data in ', TRIM(CBLADE_CSVDATA), ','
+WRITE(TLUOUT%NLU,'(A)') &
+'and make sure it is element centers (not nodes) along the blade.'
+WRITE(TLUOUT%NLU,*) ''
+!
+CALL IO_File_close(TLUOUT)
+CALL PRINT_MSG( NVERB_FATAL, 'GEN', 'EOL_BLADEDATA_ERROR', 'ERROR IN BLADE DATA' )
+!
+END SUBROUTINE EOL_BLADEDATA_ERROR
+!#########################################################
+!
diff --git a/src/MNH/eol_kine_alm.f90 b/src/MNH/eol_kine_alm.f90
new file mode 100755
index 0000000000000000000000000000000000000000..fd257b3f6f5235c4b38e600cbc08014ffa9c9871
--- /dev/null
+++ b/src/MNH/eol_kine_alm.f90
@@ -0,0 +1,334 @@
+!MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
+!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
+!MNH_LIC for details. version 1.
+!
+!#######################
+MODULE MODI_EOL_KINE_ALM
+!
+INTERFACE !----------------------------------------
+!
+SUBROUTINE EOL_KINE_ALM(KTCOUNT,KTSUBCOUNT,PTSUBSTEP,PTSTEP)
+!
+ INTEGER, INTENT(IN) :: KTCOUNT ! iteration count
+ INTEGER, INTENT(IN) :: KTSUBCOUNT ! sub iteration count
+ REAL, INTENT(IN) :: PTSUBSTEP ! sub timestep
+ REAL, INTENT(IN) :: PTSTEP ! timestep
+!
+END SUBROUTINE EOL_KINE_ALM
+!
+END INTERFACE !------------------------------------
+!
+END MODULE MODI_EOL_KINE_ALM
+!#######################
+!
+!
+!
+!
+!###################################################################
+SUBROUTINE EOL_KINE_ALM(KTCOUNT,KTSUBCOUNT,PTSUBSTEP,PTSTEP)
+!
+!!**** *EOL_KINEMATICS * -
+!!
+!! PURPOSE
+!! -------
+!! Compute positions, oritentations and velocities of all the
+!! elements of the wind turbine
+!!
+!!** METHOD
+!! ------
+!!
+!! REFERENCE
+!! ---------
+!!
+!!
+!! AUTHOR
+!! ------
+!! PA. Joulin *CNRM & IFPEN*
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 04/2017
+!! Modification 10/11/20 (PA. Joulin) Updated for a main version
+!!
+!!---------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+!* 0.1 Modules
+USE MODD_EOL_KINE_ALM
+USE MODD_EOL_ALM
+USE MODI_EOL_MATHS
+USE MODD_TIME_n, ONLY : TDTCUR
+USE MODD_CST, ONLY : XPI
+!
+IMPLICIT NONE
+!
+!* 0.2 Declarations of dummy arguments :
+!
+INTEGER, INTENT(IN) :: KTCOUNT ! iteration count
+INTEGER, INTENT(IN) :: KTSUBCOUNT ! sub iteration count
+REAL, INTENT(IN) :: PTSUBSTEP ! sub timestep
+REAL, INTENT(IN) :: PTSTEP ! timestep
+!
+!* 0.3 Local variables
+DOUBLE PRECISION, DIMENSION(3,3) :: ZORI_MAT_X, ZORI_MAT_Y, ZORI_MAT_Z
+DOUBLE PRECISION, DIMENSION(3) :: ZADD_TO_POS
+!
+DOUBLE PRECISION, DIMENSION(3) :: ZDIST_TOWO_TELT_RG ! Distance between tower elmt and tower base
+DOUBLE PRECISION, DIMENSION(3) :: ZDIST_TOWO_NELT_RG ! Distance between nacelle and base of tower
+DOUBLE PRECISION, DIMENSION(3) :: ZDIST_NAC_HUB_RG ! Distance between hub and base of nacelle
+DOUBLE PRECISION, DIMENSION(3) :: ZDIST_HUB_BLA_RG ! Distance between blade and base of hub
+DOUBLE PRECISION, DIMENSION(3) :: ZDIST_BLA_ELT_RG ! Distance between blade and elements
+!
+DOUBLE PRECISION, DIMENSION(3) :: ZPOS_ELTLE_RE ! Leading Edge (LE) position, in RE
+DOUBLE PRECISION, DIMENSION(3) :: ZPOS_ELTLE_RG ! Leading Edge (LE) position, in RG
+DOUBLE PRECISION, DIMENSION(3) :: ZPOS_ELTTE_RE ! Trailing Edge (TE) position, in RE
+DOUBLE PRECISION, DIMENSION(3) :: ZPOS_ELTTE_RG ! Trailing Edge (TE) position, in RG
+!
+REAL :: ZTIME ! TIME
+INTEGER :: JROT, JBLA, JTELT, JNELT, JBELT ! Loop control
+INTEGER :: INB_WT, INB_B, INB_BELT ! Total numbers
+INTEGER :: INB_TELT, INB_NELT ! Total numbers
+!
+!
+!---------------------------------------------------------------
+!
+!* 1. PRELIMINARIES
+! -------------
+!
+!* 1.1 Some usefull integers
+INB_WT = TFARM%NNB_TURBINES
+INB_B = TTURBINE%NNB_BLADES
+INB_TELT = 2
+INB_NELT = 2
+INB_BELT = TBLADE%NNB_BLAELT
+!
+!* 1.2 Sub-time computation
+ZTIME = TDTCUR%xtime+(KTSUBCOUNT)*PTSUBSTEP
+!
+!* 1.3 Tecplotfile : opening + headers
+IF (LTECOUTPTS) THEN
+ CALL OPEN_TECOUT(45, KTCOUNT, KTSUBCOUNT)
+END IF
+!
+!
+!---------------------------------------------------------------
+!
+!* 2. COMPUTATIONS
+! ------------
+!
+DO JROT=1, INB_WT
+!
+! ---- TOWER ----
+!
+!* T.0 Update origin positions in RG (Base position + floating)
+ XPOS_TOWO_RG(JROT,:) = XPOS_REF(:) + XPOSINI_TOWO_RG(JROT,:) &
+ + XTVEL_TOWO_RG(JROT,:)*ZTIME
+!
+!* T.1 Update orientation
+ CALL GET_ORI_MAT_X(XANGINI_TOW_RG(JROT,1) + XRVEL_RT_RG(JROT,1)*ZTIME, ZORI_MAT_X)
+ CALL GET_ORI_MAT_Y(XANGINI_TOW_RG(JROT,2) + XRVEL_RT_RG(JROT,2)*ZTIME, ZORI_MAT_Y)
+ CALL GET_ORI_MAT_Z(XANGINI_TOW_RG(JROT,3) + XRVEL_RT_RG(JROT,3)*ZTIME, ZORI_MAT_Z)
+! Compute orientation matrix
+ XMAT_RG_RT(JROT,:,:) = MATMUL(ZORI_MAT_X, MATMUL(ZORI_MAT_Y, ZORI_MAT_Z))
+!
+!* T.2 Update positions in RG
+ DO JTELT=1, INB_TELT
+ XPOS_TELT_RG(JROT,JTELT,:) = XPOS_TOWO_RG(JROT,:) &
+ + MATMUL(XMAT_RG_RT(JROT,:,:),XPOS_TELT_RT(JROT,JTELT,:))
+ END DO
+!
+!* T.3 Update structural velocities
+ DO JTELT=1, INB_TELT
+ ! Rotation of tower already in RG
+ ! Translation of elements
+ ZDIST_TOWO_TELT_RG(:) = XPOS_TELT_RG(JROT,JTELT,:) - XPOS_TOWO_RG(JROT,:)
+ XTVEL_TELT_RG(JROT,JTELT,:) = XTVEL_TOWO_RG(JROT,:) &
+ + CROSS(XRVEL_RT_RG(JROT,:),ZDIST_TOWO_TELT_RG(:))
+ ENDDO
+!
+!* T.4 Print in tecplot file
+ IF (LTECOUTPTS) THEN
+ DO JTELT=1, INB_TELT
+ CALL PRINT_TECOUT(45, XPOS_TELT_RG(JROT,JTELT,:))
+ END DO
+ END IF
+!
+!
+! ---- NACELLE ----
+!
+!* N.0 Update origin positions in RG
+ XPOS_NACO_RG(JROT,:) = XPOS_TELT_RG(JROT,INB_TELT,:) &
+ + MATMUL(XMAT_RG_RT(JROT,:,:),XPOSINI_NACO_RT(JROT,:))
+!
+!* N.1 Update orientation
+ CALL GET_ORI_MAT_X(XANGINI_NAC_RT(JROT,1) + XRVEL_RN_RT(JROT,1)*ZTIME, ZORI_MAT_X)
+ CALL GET_ORI_MAT_Y(XANGINI_NAC_RT(JROT,2) + XRVEL_RN_RT(JROT,2)*ZTIME, ZORI_MAT_Y)
+ CALL GET_ORI_MAT_Z(XANGINI_NAC_RT(JROT,3) + XRVEL_RN_RT(JROT,3)*ZTIME, ZORI_MAT_Z)
+!
+! Orientation matrix
+ XMAT_RT_RN(JROT,:,:) = MATMUL(ZORI_MAT_X, MATMUL(ZORI_MAT_Y, ZORI_MAT_Z))
+ XMAT_RG_RN(JROT,:,:) = MATMUL(XMAT_RG_RT(JROT,:,:), XMAT_RT_RN(JROT,:,:))
+!
+!* N.2 Update positions in RG
+ DO JNELT=1, INB_NELT
+ XPOS_NELT_RG(JROT,JNELT,:) = XPOS_NACO_RG(JROT,:) &
+ + MATMUL(XMAT_RG_RN(JROT,:,:),XPOS_NELT_RN(JROT,JNELT,:))
+ END DO
+!
+!* N.3 Update structural velocities
+ ! Rotation of nacelle in RG
+ XRVEL_RN_RG(JROT,:) = MATMUL(XMAT_RG_RT(JROT,:,:),XRVEL_RN_RT(JROT,:)) &
+ + XRVEL_RT_RG(JROT,:)
+ DO JNELT=1, INB_NELT
+ ! Translation of elements in RG
+ ZDIST_TOWO_NELT_RG(:) = XPOS_NELT_RG(JROT,JNELT,:) - XPOS_TOWO_RG(JROT,:)
+ XTVEL_NELT_RG(JROT,JNELT,:) = XTVEL_TOWO_RG(JROT,:) &
+ + CROSS(XRVEL_RN_RG(JROT,:),ZDIST_TOWO_NELT_RG(:))
+ END DO
+!
+!* N.4 Print in tecplot file
+ IF (LTECOUTPTS) THEN
+ DO JNELT=1, INB_NELT
+ CALL PRINT_TECOUT(45, XPOS_NELT_RG(JROT,JNELT,:))
+ END DO
+ END IF
+!
+!
+! ---- HUB ----
+!
+!* H.1 Update positions
+ XPOS_HUB_RG(JROT,:) = XPOS_NELT_RG(JROT,INB_NELT,:) &
+ + MATMUL(XMAT_RG_RN(JROT,:,:),XPOSINI_HUB_RN(JROT,:))
+!
+!* H.2 Update orientation
+ CALL GET_ORI_MAT_X(XANGINI_HUB_RN(JROT,1) + XRVEL_RH_RN(JROT,1)*ZTIME, ZORI_MAT_X)
+ CALL GET_ORI_MAT_Y(XANGINI_HUB_RN(JROT,2) + XRVEL_RH_RN(JROT,2)*ZTIME, ZORI_MAT_Y)
+ CALL GET_ORI_MAT_Z(XANGINI_HUB_RN(JROT,3) + XRVEL_RH_RN(JROT,3)*ZTIME, ZORI_MAT_Z)
+! Orientation matrix
+ XMAT_RN_RH(JROT,:,:) = MATMUL(ZORI_MAT_X, MATMUL(ZORI_MAT_Y, ZORI_MAT_Z))
+ XMAT_RG_RH(JROT,:,:) = MATMUL(XMAT_RG_RN(JROT,:,:), XMAT_RN_RH(JROT,:,:))
+ XMAT_RH_RG(JROT,:,:) = TRANSPOSE(XMAT_RG_RH(JROT,:,:))
+!
+!* H.3 Update structural velocities
+! Rotation of hub in RG
+ XRVEL_RH_RG(JROT,:) = MATMUL(XMAT_RG_RH(JROT,:,:),XRVEL_RH_RN(JROT,:)) &
+ + XRVEL_RN_RG(JROT,:)
+! Translation of hub in RG
+ ZDIST_NAC_HUB_RG(:) = XPOS_HUB_RG(JROT,:) - XPOS_NELT_RG(JROT,INB_NELT,:)
+ XTVEL_HUB_RG(JROT,:) = XTVEL_NELT_RG(JROT,INB_NELT,:) + CROSS(XRVEL_RH_RG(JROT,:),ZDIST_NAC_HUB_RG(:))
+!
+!* H.4 Print in tecplot file
+ IF (LTECOUTPTS) THEN
+ CALL PRINT_TECOUT(45, XPOS_HUB_RG(JROT,:))
+ END IF
+!
+!
+! ---- BLADES ----
+!
+ DO JBLA=1, INB_B
+!* B.1 Update positions
+ XPOS_BLA_RG(JROT,JBLA,:) = XPOS_HUB_RG(JROT,:) &
+ + MATMUL(XMAT_RG_RH(JROT,:,:),XPOSINI_BLA_RH(JROT,JBLA,:))
+!
+!* B.2 Update orientation
+ CALL GET_ORI_MAT_X(XANGINI_BLA_RH(JROT,JBLA,1) + XRVEL_RB_RH(JROT,JBLA,1)*ZTIME, ZORI_MAT_X)
+ CALL GET_ORI_MAT_Y(XANGINI_BLA_RH(JROT,JBLA,2) + XRVEL_RB_RH(JROT,JBLA,2)*ZTIME, ZORI_MAT_Y)
+ CALL GET_ORI_MAT_Z(XANGINI_BLA_RH(JROT,JBLA,3) + XRVEL_RB_RH(JROT,JBLA,3)*ZTIME, ZORI_MAT_Z)
+! Orientation matrix
+ XMAT_RH_RB(JROT,JBLA,:,:) = MATMUL(ZORI_MAT_X, MATMUL(ZORI_MAT_Y, ZORI_MAT_Z))
+ XMAT_RG_RB(JROT,JBLA,:,:) = MATMUL(XMAT_RG_RH(JROT,:,:), XMAT_RH_RB(JROT,JBLA,:,:))
+!
+!* B.3 Update structural velocities
+! Rotation of blade in RG
+ XRVEL_RB_RG(JROT,JBLA,:) = XRVEL_RH_RG(JROT,:) &
+ + MATMUL(XMAT_RG_RB(JROT,JBLA,:,:),XRVEL_RB_RH(JROT,JBLA,:))
+! Translation of blade in RG
+ ZDIST_HUB_BLA_RG(:) = XPOS_BLA_RG(JROT,JBLA,:) - XPOS_HUB_RG(JROT,:)
+ XTVEL_BLA_RG(JROT,JBLA,:) = XTVEL_HUB_RG(JROT,:) &
+ + CROSS(XRVEL_RB_RG(JROT,JBLA,:),ZDIST_HUB_BLA_RG(:))
+!
+!* B.4 Print in tecplot file
+ IF (LTECOUTPTS) THEN
+ CALL PRINT_TECOUT(45, XPOS_BLA_RG(JROT,JBLA,:))
+ END IF
+!
+!
+! ---- ELEMENTS ----
+!* E.0 Positioning sections (cuts) in RG
+ DO JBELT=1, INB_BELT+1
+ XPOS_SEC_RG(JROT,JBLA,JBELT,:) = XPOS_BLA_RG(JROT,JBLA,:) &
+ + MATMUL(XMAT_RG_RB(JROT,JBLA,:,:),XPOS_SEC_RB(JROT,JBLA,JBELT,:))
+ ENDDO
+!
+!
+ DO JBELT=1, INB_BELT
+!* E.1 Positioning sections centers (application points) in RG
+ XPOS_ELT_RG(JROT,JBLA,JBELT,:) = XPOS_BLA_RG(JROT,JBLA,:) &
+ + MATMUL(XMAT_RG_RB(JROT,JBLA,:,:),XPOS_ELT_RB(JROT,JBLA,JBELT,:))
+!* E.2 Update orientation
+ CALL GET_ORI_MAT_X(XANGINI_ELT_RB(JROT,JBLA,JBELT,1) &
+ + XRVEL_RE_RB(JROT,JBLA,JBELT,1)*ZTIME, ZORI_MAT_X)
+ CALL GET_ORI_MAT_Y(XANGINI_ELT_RB(JROT,JBLA,JBELT,2) &
+ + XRVEL_RE_RB(JROT,JBLA,JBELT,2)*ZTIME, ZORI_MAT_Y)
+ CALL GET_ORI_MAT_Z(XANGINI_ELT_RB(JROT,JBLA,JBELT,3) &
+ + XRVEL_RE_RB(JROT,JBLA,JBELT,3)*ZTIME, ZORI_MAT_Z)
+! Orientation matrix
+ XMAT_RB_RE(JROT,JBLA,JBELT,:,:) = MATMUL(ZORI_MAT_X, MATMUL(ZORI_MAT_Y,ZORI_MAT_Z))
+ XMAT_RG_RE(JROT,JBLA,JBELT,:,:) = MATMUL(XMAT_RG_RB(JROT,JBLA,:,:), XMAT_RB_RE(JROT,JBLA,JBELT,:,:))
+ XMAT_RE_RG(JROT,JBLA,JBELT,:,:) = TRANSPOSE(XMAT_RG_RE(JROT,JBLA,JBELT,:,:))
+!
+!* E.3 Update structural velocities
+! Rotation of elements in RG
+ XRVEL_RE_RG(JROT,JBLA,JBELT,:) = XRVEL_RB_RG(JROT,JBLA,:) &
+ + MATMUL(XMAT_RG_RE(JROT,JBLA,JBELT,:,:),&
+ XRVEL_RE_RB(JROT,JBLA,JBELT,:))
+! Translation of elements in RG
+ ZDIST_BLA_ELT_RG(:) = XPOS_ELT_RG(JROT,JBLA,JBELT,:) - XPOS_BLA_RG(JROT,JBLA,:)
+ XTVEL_ELT_RG(JROT,JBLA,JBELT,:) = XTVEL_BLA_RG(JROT,JBLA,:) &
+ + CROSS(XRVEL_RE_RG(JROT,JBLA,JBELT,:),ZDIST_BLA_ELT_RG(:))
+ XTVEL_ELT_RE(JROT,JBLA,JBELT,:) = MATMUL(XMAT_RE_RG(JROT,JBLA,JBELT,:,:),&
+ XTVEL_ELT_RG(JROT,JBLA,JBELT,:))
+!
+!* E.4 Print in tecplot file
+ IF (LTECOUTPTS) THEN
+ CALL PRINT_TECOUT(45, XPOS_ELT_RG(JROT,JBLA,JBELT,:))
+ END IF
+!
+! ---- Leading Edge and Trailing Edge ----
+! It is just to have a tecplot more beautiful
+! For the moment, they are useless for computations
+ IF (LTECOUTPTS) THEN
+! LE.1 Update Positions
+ ZPOS_ELTLE_RE(1) = 0d0
+ ZPOS_ELTLE_RE(2) = - 1d0/4d0 * XCHORD_ELT(JROT,JBLA,JBELT)
+ ZPOS_ELTLE_RE(3) = 0d0
+ ZPOS_ELTTE_RE(1) = 0d0
+ ZPOS_ELTTE_RE(2) = + 3d0/4d0 * XCHORD_ELT(JROT,JBLA,JBELT)
+ ZPOS_ELTTE_RE(3) = 0d0
+!
+ ZPOS_ELTLE_RG(:) = XPOS_ELT_RG(JROT,JBLA,JBELT,:) &
+ + MATMUL(XMAT_RG_RE(JROT,JBLA,JBELT,:,:), &
+ ZPOS_ELTLE_RE(:))
+ ZPOS_ELTTE_RG(:) = XPOS_ELT_RG(JROT,JBLA,JBELT,:) &
+ + MATMUL(XMAT_RG_RE(JROT,JBLA,JBELT,:,:), &
+ ZPOS_ELTTE_RE(:))
+
+!* LE.2 Print in tecplot file
+ CALL PRINT_TECOUT(45, ZPOS_ELTLE_RG(:))
+ CALL PRINT_TECOUT(45, ZPOS_ELTTE_RG(:))
+ END IF
+!
+!
+ END DO ! Blade element loop
+ END DO ! Blade loop
+END DO ! Rotor loop
+!
+! Closing tec file
+IF (LTECOUTPTS) THEN
+ CLOSE(45)
+END IF
+!
+END SUBROUTINE EOL_KINE_ALM
diff --git a/src/MNH/eol_main.f90 b/src/MNH/eol_main.f90
new file mode 100755
index 0000000000000000000000000000000000000000..5042b4418412b0babe3971adaa057727aac35d1b
--- /dev/null
+++ b/src/MNH/eol_main.f90
@@ -0,0 +1,283 @@
+!MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
+!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
+!MNH_LIC for details. version 1.
+! #######################
+ MODULE MODI_EOL_MAIN
+! #######################
+!
+INTERFACE
+!
+SUBROUTINE EOL_MAIN(KTCOUNT, PTSTEP, &
+ PDXX, PDYY, PDZZ, &
+ PRHODJ, PUT, PVT, PWT, &
+ PRUS, PRVS, PRWS )
+!
+INTEGER, INTENT(IN) :: KTCOUNT ! iteration count
+REAL, INTENT(IN) :: PTSTEP ! timestep except
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX,PDYY,PDZZ ! mesh size
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! dry Density * Jacobian
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PUT, PVT, PWT ! wind speed variables
+!
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRUS, PRVS, PRWS ! Sources of Momentum
+!
+END SUBROUTINE EOL_MAIN
+!
+END INTERFACE
+!
+END MODULE MODI_EOL_MAIN
+!
+! ###################################################################
+ SUBROUTINE EOL_MAIN(KTCOUNT, PTSTEP, &
+ PDXX, PDYY, PDZZ, &
+ PRHODJ, PUT, PVT, PWT, &
+ PRUS, PRVS, PRWS )
+! ###################################################################
+!
+!!**** *EOL_MAIN * -
+!!
+!! PURPOSE
+!! -------
+!! It is possible to include wind turbines parameterization in Meso-NH,
+!! and several models are available. EOL_MAIN is the main subroutine
+!! to compute the aerodynamics of the wind turbine, according to the
+!! model chosen.
+!!
+!!** METHOD
+!! ------
+!!
+!! REFERENCE
+!! ---------
+!!
+!!
+!! AUTHOR
+!! ------
+!! PA. Joulin *CNRM & IFPEN*
+!!
+!! MODIFICATIONS
+!! -------------
+!! 21/10/20 Original
+!!
+!!---------------------------------------------------------------
+!
+!
+!* 0. DECLARATIONS
+! ------------
+!
+! To work with wind turbines
+USE MODD_EOL_MAIN
+USE MODI_EOL_ADNR
+USE MODI_EOL_ALM
+USE MODI_EOL_SMEAR
+! To play with MPI
+USE MODD_ARGSLIST_ll, ONLY: LIST_ll
+USE MODE_ll , ONLY: ADD3DFIELD_ll
+USE MODE_ll , ONLY: UPDATE_HALO_ll
+USE MODE_ll , ONLY: CLEANLIST_ll
+! To use some toolkit
+USE MODI_SHUMAN , ONLY: MXF, MYF, MZF
+USE MODI_SHUMAN , ONLY: MXM, MYM, MZM
+!
+use modd_budget, only: lbudget_u, lbudget_v, lbudget_w, NBUDGET_U, NBUDGET_V, NBUDGET_W, tbudgets
+use mode_budget, only: budget_store_init, budget_store_end
+!
+!
+IMPLICIT NONE
+!
+!* 0.1 Declarations of dummy arguments :
+!
+INTEGER, INTENT(IN) :: KTCOUNT ! iteration count
+REAL, INTENT(IN) :: PTSTEP ! timestep except
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX,PDYY,PDZZ ! mesh size
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! dry Density * Jacobian
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PUT, PVT, PWT ! Wind speed
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRUS, PRVS, PRWS ! Sources of Momentum
+!
+!
+!* 0.2 Declarations of local variables :
+!
+! Pointeurs and exchanges
+TYPE(LIST_ll), POINTER :: TZFIELDS_W_ll ! Field list of Wind for exchange
+TYPE(LIST_ll), POINTER :: TZFIELDS_F_ll ! Field list of aero Forces for exchange
+TYPE(LIST_ll), POINTER :: TZFIELDS_S_ll ! Field list of Smeared aero forces for exchange
+TYPE(LIST_ll), POINTER :: TZFIELDS_R_ll ! Field list of mnh foRces for exchange
+INTEGER :: IINFO ! Info integer
+INTEGER :: IKU ! Vertical size of the domain
+!
+! ABL
+REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: ZUT_M, ZVT_M, ZWT_M ! Wind speed at mass point
+REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: ZRHO_M ! Air density at mass point
+!
+!
+!* 0.3 Implicit arguments
+!* From MODD_EOL_MAIN
+! Aerodynamic forces in cartesian mesh
+!REAL, DIMENSION(:,:,:), ALLOCATABLE :: XFX_RG ! Along X in RG frame [F]
+!REAL, DIMENSION(:,:,:), ALLOCATABLE :: XFY_RG ! Along Y in RG frame [F]
+!REAL, DIMENSION(:,:,:), ALLOCATABLE :: XFZ_RG ! Along Z in RG frame [F]
+! Smeared forces
+!REAL, DIMENSION(:,:,:), ALLOCATABLE :: XFX_SMR_RG ! Along X in RG frame [F]
+!REAL, DIMENSION(:,:,:), ALLOCATABLE :: XFY_SMR_RG ! Along Y in RG frame [F]
+!REAL, DIMENSION(:,:,:), ALLOCATABLE :: XFZ_SMR_RG ! ALong Z in RG frame [F]
+!
+!* From NAM_EOL namelist
+!CHARACTER(LEN=4) :: CMETH_EOL ! Aerodynamic method
+!CHARACTER(LEN=4) :: CSMEAR ! Type of smearing
+!
+!
+!--------------------------------------------------------
+!
+!* 1. INITIALIZATIONS
+! ---------------
+!
+!* 1.1 Indices
+!
+IKU = SIZE(PUT,3)
+!
+!* 1.2 Pointers
+!
+NULLIFY(TZFIELDS_W_ll)
+NULLIFY(TZFIELDS_F_ll)
+NULLIFY(TZFIELDS_S_ll)
+NULLIFY(TZFIELDS_R_ll)
+!
+!* 1.3 Forces
+!
+XFX_RG(:,:,:) = 0.
+XFY_RG(:,:,:) = 0.
+XFZ_RG(:,:,:) = 0.
+XFX_SMR_RG(:,:,:) = 0.
+XFY_SMR_RG(:,:,:) = 0.
+XFZ_SMR_RG(:,:,:) = 0.
+!
+!
+!-----------------------------------------------------------------------
+!
+!* 2. COMPUTES VELOCITY COMPONENTS AND DENSITY AT MASS POINT
+! ------------------------------------------------------
+!
+!* 2.1 Sharing the input
+!
+CALL ADD3DFIELD_ll( TZFIELDS_W_ll,PUT, 'EOL_MAIN::PUT')
+CALL ADD3DFIELD_ll( TZFIELDS_W_ll,PWT, 'EOL_MAIN::PWT')
+CALL ADD3DFIELD_ll( TZFIELDS_W_ll,PVT, 'EOL_MAIN::PVT')
+CALL UPDATE_HALO_ll(TZFIELDS_W_ll,IINFO)
+CALL CLEANLIST_ll( TZFIELDS_W_ll)
+!
+!* 2.2 Masss point evaluation
+!
+ZUT_M(:,:,:) = MXF( PUT(:,:,:) )
+ZVT_M(:,:,:) = MYF( PVT(:,:,:) )
+ZWT_M(:,:,:) = MZF( PWT(:,:,:) )
+ZRHO_M(:,:,:) = PRHODJ(:,:,:)/(PDXX(:,:,:)*PDYY(:,:,:)*PDZZ(:,:,:))
+!
+!* 2.3 Sharing the new wind
+!
+CALL ADD3DFIELD_ll( TZFIELDS_W_ll,ZUT_M, 'EOL_MAIN::ZUT_M')
+CALL ADD3DFIELD_ll( TZFIELDS_W_ll,ZWT_M, 'EOL_MAIN::ZWT_M')
+CALL ADD3DFIELD_ll( TZFIELDS_W_ll,ZVT_M, 'EOL_MAIN::ZVT_M')
+CALL UPDATE_HALO_ll(TZFIELDS_W_ll,IINFO)
+CALL CLEANLIST_ll( TZFIELDS_W_ll)
+!
+!
+!--------------------------------------------------------
+!
+!* 3. COMPUTES AERODYNAMICS FORCES
+! ----------------------------
+!
+!* 3.1 Model selection
+!
+!
+SELECT CASE(CMETH_EOL)
+!
+ CASE('ADNR') ! Actuator Disc Non-Rotating
+ CALL EOL_ADNR(PDXX, PDYY, PDZZ, &
+ ZRHO_M, &
+ ZUT_M, &
+ XFX_RG )
+!
+ CASE('ALM') ! Actuator Line Method
+ CALL EOL_ALM(KTCOUNT, PTSTEP, &
+ PDXX, PDYY, PDZZ, &
+ ZRHO_M, &
+ ZUT_M, ZVT_M, ZWT_M, &
+ XFX_RG, XFY_RG, XFZ_RG )
+!
+END SELECT
+!
+!* 3.2 Sharing 3D field
+!
+CALL ADD3DFIELD_ll( TZFIELDS_F_ll,XFX_RG, 'EOL_MAIN::XFX_RG' )
+CALL ADD3DFIELD_ll( TZFIELDS_F_ll,XFY_RG, 'EOL_MAIN::XFY_RG' )
+CALL ADD3DFIELD_ll( TZFIELDS_F_ll,XFZ_RG, 'EOL_MAIN::XFZ_RG' )
+CALL UPDATE_HALO_ll(TZFIELDS_F_ll,IINFO)
+CALL CLEANLIST_ll( TZFIELDS_F_ll)
+!
+!
+!--------------------------------------------------------
+!
+!* 4. SMEARING THE FORCES
+! -------------------
+!
+!* 4.1 Smearing technique selection
+!
+SELECT CASE (CSMEAR)
+!
+ CASE( 'NULL' ) ! No smearing
+ XFX_SMR_RG(:,:,:) = XFX_RG(:,:,:)
+ XFY_SMR_RG(:,:,:) = XFY_RG(:,:,:)
+ XFZ_SMR_RG(:,:,:) = XFZ_RG(:,:,:)
+!
+ CASE( '1LIN' ) ! Linear smearing
+ CALL SMEAR_1LIN(XFX_RG, &
+ XFX_SMR_RG)
+!
+ CASE( '3LIN' ) ! Linear smearing
+ CALL SMEAR_3LIN(XFX_RG, &
+ XFY_RG, &
+ XFZ_RG, &
+ XFX_SMR_RG,&
+ XFY_SMR_RG,&
+ XFZ_SMR_RG)
+!
+END SELECT
+!
+!* 4.2 Sharing 3D field
+!
+CALL ADD3DFIELD_ll( TZFIELDS_S_ll,XFX_SMR_RG, 'EOL_MAIN::XFX_SMR_RG' )
+CALL ADD3DFIELD_ll( TZFIELDS_S_ll,XFY_SMR_RG, 'EOL_MAIN::XFY_SMR_RG' )
+CALL ADD3DFIELD_ll( TZFIELDS_S_ll,XFZ_SMR_RG, 'EOL_MAIN::XFZ_SMR_RG' )
+CALL UPDATE_HALO_ll(TZFIELDS_S_ll,IINFO)
+CALL CLEANLIST_ll( TZFIELDS_S_ll)
+!
+!
+!-------------------------------------------------------------------------------
+!
+!* 5. ADDING THE FORCES TO THE FIELD
+! ------------------------------
+!
+!* 5.1 Adding them to flux points, rotor->wind
+!
+if (lbudget_u) call Budget_store_init( tbudgets(NBUDGET_U), 'EOL', prus(:,:,:) )
+if (lbudget_v) call Budget_store_init( tbudgets(NBUDGET_V), 'EOL', prvs(:,:,:) )
+if (lbudget_w) call Budget_store_init( tbudgets(NBUDGET_W), 'EOL', prws(:,:,:) )
+!
+PRUS(:,:,:)=PRUS(:,:,:)-MXM(XFX_SMR_RG(:,:,:))
+PRVS(:,:,:)=PRVS(:,:,:)-MYM(XFY_SMR_RG(:,:,:))
+PRWS(:,:,:)=PRWS(:,:,:)-MZM(XFZ_SMR_RG(:,:,:))
+!
+if (lbudget_u) call Budget_store_end( tbudgets(NBUDGET_U), 'EOL', prus(:,:,:) )
+if (lbudget_v) call Budget_store_end( tbudgets(NBUDGET_V), 'EOL', prvs(:,:,:) )
+if (lbudget_w) call Budget_store_end( tbudgets(NBUDGET_W), 'EOL', prws(:,:,:) )
+!
+!
+!* 5.2 Sharing the field
+!
+CALL ADD3DFIELD_ll( TZFIELDS_R_ll,PRUS,'EOL_MAIN::PRUS' )
+CALL ADD3DFIELD_ll( TZFIELDS_R_ll,PRVS,'EOL_MAIN::PRVS' )
+CALL ADD3DFIELD_ll( TZFIELDS_R_ll,PRWS,'EOL_MAIN::PRWS' )
+CALL UPDATE_HALO_ll(TZFIELDS_R_ll,IINFO)
+CALL CLEANLIST_ll( TZFIELDS_R_ll)
+!
+END SUBROUTINE EOL_MAIN
diff --git a/src/MNH/eol_maths.f90 b/src/MNH/eol_maths.f90
new file mode 100755
index 0000000000000000000000000000000000000000..6d2d659cca7b35c038041916281939854b3179ff
--- /dev/null
+++ b/src/MNH/eol_maths.f90
@@ -0,0 +1,375 @@
+!MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
+!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
+!MNH_LIC for details. version 1.
+! #######################
+ MODULE MODI_EOL_MATHS
+! #######################
+!
+INTERFACE
+!
+FUNCTION CROSS(PA, PB)
+ DOUBLE PRECISION, DIMENSION(3) :: CROSS
+ DOUBLE PRECISION, DIMENSION(3), INTENT(IN) :: PA, PB
+END FUNCTION CROSS
+!
+FUNCTION NORM(PA)
+ DOUBLE PRECISION :: NORM
+ DOUBLE PRECISION, DIMENSION(3), INTENT(IN) :: PA
+END FUNCTION NORM
+!
+SUBROUTINE GET_ORI_MAT_X(PTHETA, PORI_MAT_X)
+ DOUBLE PRECISION, INTENT(IN) :: PTHETA ! Angle
+ DOUBLE PRECISION, DIMENSION(3,3), INTENT(OUT) :: PORI_MAT_X ! Matrix
+END SUBROUTINE GET_ORI_MAT_X
+!
+SUBROUTINE GET_ORI_MAT_Y(PTHETA, PORI_MAT_Y)
+ DOUBLE PRECISION, INTENT(IN) :: PTHETA ! Angle
+ DOUBLE PRECISION, DIMENSION(3,3), INTENT(OUT) :: PORI_MAT_Y ! Matrix
+END SUBROUTINE GET_ORI_MAT_Y
+!
+SUBROUTINE GET_ORI_MAT_Z(PTHETA, PORI_MAT_Z)
+ DOUBLE PRECISION, INTENT(IN) :: PTHETA ! Angle
+ DOUBLE PRECISION, DIMENSION(3,3), INTENT(OUT) :: PORI_MAT_Z ! Matrix
+END SUBROUTINE GET_ORI_MAT_Z
+!
+FUNCTION INTERP_SPLCUB(PAV, PX, PY)
+ REAL :: INTERP_SPLCUB ! interface
+ REAL, INTENT(IN) :: PAV ! Abscissa where spline is to be evaluate
+ REAL, DIMENSION(:), INTENT(IN) :: PX, PY
+END FUNCTION INTERP_SPLCUB
+!
+FUNCTION INTERP_LIN8NB(PPOS, KI, KJ, KK, PVAR, PZH)
+ REAL :: INTERP_LIN8NB ! interface
+ REAL, DIMENSION(3), INTENT(IN) :: PPOS ! Position where we want to evaluate
+ INTEGER, INTENT(IN) :: KI, KJ, KK ! Meso-NH cell index
+ REAL, DIMENSION(:,:,:), INTENT(IN) :: PVAR,PZH ! Variable to interpolate
+END FUNCTION INTERP_LIN8NB
+!
+END INTERFACE
+!
+END MODULE MODI_EOL_MATHS
+!-------------------------------------------------------------------
+!
+!!**** *EOL_MATHS* -
+!!
+!! PURPOSE
+!! -------
+!! Some usefull tools for wind turbine study
+!!
+!! AUTHOR
+!! ------
+!! PA. Joulin *CNRM & IFPEN*
+!!
+!! MODIFICATIONS
+!! -------------
+!! 04/2018 Original
+!!
+!!---------------------------------------------------------------
+!#########################################################
+FUNCTION CROSS(PA, PB)
+! Vectorial product 3D : PA * PB
+!
+ DOUBLE PRECISION, DIMENSION(3) :: CROSS
+ DOUBLE PRECISION, DIMENSION(3), INTENT(IN) :: PA, PB
+!
+ CROSS(1) = PA(2) * PB(3) - PA(3) * PB(2)
+ CROSS(2) = PA(3) * PB(1) - PA(1) * PB(3)
+ CROSS(3) = PA(1) * PB(2) - PA(2) * PB(1)
+!
+END FUNCTION CROSS
+!#########################################################
+!
+!#########################################################
+FUNCTION NORM(PA)
+! Eulerian norm of 3D vector :
+!
+ DOUBLE PRECISION :: NORM
+ DOUBLE PRECISION, DIMENSION(3), INTENT(IN) :: PA
+!
+ NORM = SQRT( PA(1)**2 + PA(2)**2 + PA(3)**2 )
+!
+END FUNCTION NORM
+!
+!
+!#########################################################
+SUBROUTINE GET_ORI_MAT_X(PTHETA, PORI_MAT_X)
+! Rotation matrix of PTHETA angle around X
+!
+ DOUBLE PRECISION, INTENT(IN) :: PTHETA ! Angle
+ DOUBLE PRECISION, DIMENSION(3,3), INTENT(OUT) :: PORI_MAT_X ! Matrix
+!
+ PORI_MAT_X (1,1) = 1d0
+ PORI_MAT_X (1,2) = 0d0
+ PORI_MAT_X (1,3) = 0d0
+ PORI_MAT_X (2,1) = 0d0
+ PORI_MAT_X (2,2) = +COS(PTHETA)
+ PORI_MAT_X (2,3) = -SIN(PTHETA)
+ PORI_MAT_X (3,1) = 0d0
+ PORI_MAT_X (3,2) = +SIN(PTHETA)
+ PORI_MAT_X (3,3) = +COS(PTHETA)
+!
+END SUBROUTINE GET_ORI_MAT_X
+!#########################################################
+!
+!#########################################################
+SUBROUTINE GET_ORI_MAT_Y(PTHETA, PORI_MAT_Y)
+! Rotation matrix of PTHETA angle around Y
+!
+ DOUBLE PRECISION, INTENT(IN) :: PTHETA ! Angle
+ DOUBLE PRECISION, DIMENSION(3,3), INTENT(OUT) :: PORI_MAT_Y ! Matrix
+!
+ PORI_MAT_Y (1,1) = +COS(PTHETA)
+ PORI_MAT_Y (1,2) = 0d0
+ PORI_MAT_Y (1,3) = +SIN(PTHETA)
+ PORI_MAT_Y (2,1) = 0d0
+ PORI_MAT_Y (2,2) = 1d0
+ PORI_MAT_Y (2,3) = 0d0
+ PORI_MAT_Y (3,1) = -SIN(PTHETA)
+ PORI_MAT_Y (3,2) = 0d0
+ PORI_MAT_Y (3,3) = +COS(PTHETA)
+!
+END SUBROUTINE GET_ORI_MAT_Y
+!#########################################################
+!
+!#########################################################
+SUBROUTINE GET_ORI_MAT_Z(PTHETA, PORI_MAT_Z)
+! Rotation matrix of PTHETA angle around Z
+!
+ DOUBLE PRECISION, INTENT(IN) :: PTHETA ! Angle
+ DOUBLE PRECISION, DIMENSION(3,3), INTENT(OUT) :: PORI_MAT_Z ! Matrix
+!
+ PORI_MAT_Z (1,1) = +COS(PTHETA)
+ PORI_MAT_Z (1,2) = -SIN(PTHETA)
+ PORI_MAT_Z (1,3) = 0d0
+ PORI_MAT_Z (2,1) = +SIN(PTHETA)
+ PORI_MAT_Z (2,2) = +COS(PTHETA)
+ PORI_MAT_Z (2,3) = 0d0
+ PORI_MAT_Z (3,1) = 0d0
+ PORI_MAT_Z (3,2) = 0d0
+ PORI_MAT_Z (3,3) = 1d0
+!
+END SUBROUTINE GET_ORI_MAT_Z
+!#########################################################
+!
+!#########################################################
+FUNCTION INTERP_SPLCUB(PAV, PX, PY)
+! adapted from https://ww2.odu.edu/~agodunov/computing/programs/book2/Ch01/spline.f90
+!
+IMPLICIT NONE
+REAL, INTENT(IN) :: PAV ! Abscissa where spline is to be evaluate
+REAL, DIMENSION(:), INTENT(IN) :: PX, PY
+!
+INTEGER :: INBVAL ! Nb points of data
+REAL, ALLOCATABLE :: ZCOEF1(:),ZCOEF2(:),ZCOEF3(:) ! Coefficients
+!
+INTEGER :: II, IJ, IBOT, ITOP, IMID
+REAL :: ZH
+REAL :: PDX
+REAL :: INTERP_SPLCUB ! function
+!
+! --------- Intialisations ---------
+INBVAL = SIZE(PX)
+ALLOCATE(ZCOEF1(INBVAL))
+ALLOCATE(ZCOEF2(INBVAL))
+ALLOCATE(ZCOEF3(INBVAL))
+!
+! --------- Calculs des coefficients ---------
+!
+! - Check size of input data
+IF (INBVAL < 2 ) THEN
+ RETURN
+END IF
+IF (INBVAL < 3 ) THEN
+ ZCOEF1(1) = (PY(2)-PY(1))/(PX(2)-PX(1))
+ ZCOEF2(1) = 0.
+ ZCOEF3(1) = 0.
+ ZCOEF1(2) = ZCOEF1(1)
+ ZCOEF2(2) = 0.
+ ZCOEF3(2) = 0.
+ RETURN
+END IF
+!
+! - Preliminaries
+ZCOEF3(1) = PX(2) - PX(1)
+ZCOEF2(2) = (PY(2) - PY(1))/ZCOEF3(1)
+DO II = 2, INBVAL-1
+ ZCOEF3(II) = PX(II+1) - PX(II)
+ ZCOEF1(II) = 2.0 * (ZCOEF3(II-1) + ZCOEF3(II))
+ ZCOEF2(II+1) = (PY(II+1) - PY(II))/ZCOEF3(II)
+ ZCOEF2(II) = ZCOEF2(II+1) - ZCOEF2(II)
+END DO
+!
+! - Boundaries
+ZCOEF1(1) = - ZCOEF3(1)
+ZCOEF1(INBVAL) = - ZCOEF3(INBVAL-1)
+ZCOEF2(1) = 0.0
+ZCOEF2(INBVAL) = 0.0
+IF (INBVAL /= 3) THEN
+ ZCOEF2(1) = ZCOEF2(3)/(PX(4)-PX(2)) - ZCOEF2(2)/(PX(3)-PX(1))
+ ZCOEF2(INBVAL) = ZCOEF2(INBVAL-1)/(PX(INBVAL)-PX(INBVAL-2)) &
+ -ZCOEF2(INBVAL-2)/(PX(INBVAL-1)-PX(INBVAL-3))
+ ZCOEF2(1) = ZCOEF2(1)*ZCOEF3(1)**2 / (PX(4)-PX(1))
+ ZCOEF2(INBVAL) =-ZCOEF2(INBVAL)*ZCOEF3(INBVAL-1)**2/(PX(INBVAL)-PX(INBVAL-3))
+END IF
+!
+! - Forward elemination
+DO II = 2, INBVAL
+ ZH = ZCOEF3(II-1)/ZCOEF1(II-1)
+ ZCOEF1(II) = ZCOEF1(II) - ZH*ZCOEF3(II-1)
+ ZCOEF2(II) = ZCOEF2(II) - ZH*ZCOEF2(II-1)
+END DO
+!
+! - Back substitution
+ZCOEF2(INBVAL) = ZCOEF2(INBVAL)/ZCOEF1(INBVAL)
+DO IJ = 1, INBVAL-1
+ II = INBVAL-IJ
+ ZCOEF2(II) = (ZCOEF2(II) - ZCOEF3(II)*ZCOEF2(II+1))/ZCOEF1(II)
+END DO
+!
+! - Spline coefficient calculations
+ZCOEF1(INBVAL) = (PY(INBVAL) - PY(INBVAL-1))/ZCOEF3(INBVAL-1) &
+ + ZCOEF3(INBVAL-1)*(ZCOEF2(INBVAL-1) + 2.0*ZCOEF2(INBVAL))
+DO II = 1, INBVAL-1
+ ZCOEF1(II) = (PY(II+1) - PY(II))/ZCOEF3(II) &
+ - ZCOEF3(II)*(ZCOEF2(II+1) + 2.0*ZCOEF2(II))
+ ZCOEF3(II) = (ZCOEF2(II+1) - ZCOEF2(II))/ZCOEF3(II)
+ ZCOEF2(II) = 3.0*ZCOEF2(II)
+END DO
+ZCOEF2(INBVAL) = 3.0*ZCOEF2(INBVAL)
+ZCOEF3(INBVAL) = ZCOEF3(INBVAL-1)
+
+! --------- Spline cubic interpolation ---------
+
+! If the absciss PAV is out of range
+! The ordinate will be the limit value (left or right)
+IF (PAV <= PX(1)) THEN
+ INTERP_SPLCUB = PY(1)
+ RETURN
+END IF
+IF (PAV >= PX(INBVAL)) THEN
+ INTERP_SPLCUB = PY(INBVAL)
+ RETURN
+END IF
+
+! Dichotomie research for IBOT, tq : PX(IBOT) <= PAV <= PX(IBOT+1)
+IBOT = 1
+ITOP = INBVAL +1
+DO WHILE (ITOP > IBOT+1)
+ IMID = (IBOT + ITOP)/2
+ IF (PAV < PX(IMID)) THEN
+ ITOP = IMID
+ ELSE
+ IBOT = IMID
+ END IF
+END DO
+
+! Evaluation of spline interpolation
+PDX = PAV - PX(IBOT)
+INTERP_SPLCUB = PY(IBOT)+PDX*(ZCOEF1(IBOT)+PDX*(ZCOEF2(IBOT)+PDX*ZCOEF3(IBOT)))
+
+! Endings
+DEALLOCATE(ZCOEF1)
+DEALLOCATE(ZCOEF2)
+DEALLOCATE(ZCOEF3)
+
+END FUNCTION INTERP_SPLCUB
+!#########################################################
+!
+!#########################################################
+FUNCTION INTERP_LIN8NB(PPOS, KI, KJ, KK, PVAR, PZH)
+!
+USE MODD_GRID_n, ONLY: XXHAT,XYHAT
+!
+REAL :: INTERP_LIN8NB ! Return
+REAL, DIMENSION(3), INTENT(IN) :: PPOS ! Position where we want to evaluate
+INTEGER, INTENT(IN) :: KI, KJ, KK ! Meso-NH cell index
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PVAR ! Variable to interpolate
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PZH ! Vertical height to interpolate
+!
+INTEGER :: IIP, IJP, IKP ! Previous cell index : P = i + 1
+INTEGER :: IIN, IJN, IKN ! Next cell index : N = i - 1
+!
+REAL :: ZUXNN, ZUXNP, ZUXPP, ZUXPN ! Interpolated variables (VAR) in X plane (VAR = A*POS + B)
+!
+REAL :: ZHXNN, ZHXNP, ZHXPP, ZHXPN ! Interpotaled variables (VAR) in X plane (VAR = A*POS + B)
+!
+REAL :: ZUXN, ZUXP ! Interpolated variables (VAR) in Y plance (VAR = A*POS + B)
+!
+!
+REAL :: ZALPHAX, ZALPHAY, ZALPHAZ ! Interpolated variables (VAR) in Z plane (VAR = A*POS + B)
+
+REAL :: ZUX ! Interpolated variable (VAR) in Z plane (VAR = A*POS + B)
+!
+! -----------------------------------------------
+!
+! FINDING 8 NEIGHBOORS
+! -- X axis
+IF (PPOS(1) <= 0.5*(XXHAT(KI)+XXHAT(KI+1))) THEN
+ IIP = KI - 1
+ IIN = KI
+ELSE
+ IIP = KI
+ IIN = KI + 1
+END IF
+! -- Y axis
+IF (PPOS(2) <= 0.5*((XYHAT(KJ)+XYHAT(KJ+1)))) THEN
+ IJP = KJ - 1
+ IJN = KJ
+ELSE
+ IJP = KJ
+ IJN = KJ + 1
+END IF
+! -- Z axis
+IF (PPOS(3) <= PZH(KI,KJ,KK)) THEN
+ IKP = KK - 1
+ IKN = KK
+ELSE
+ IKP = KK
+ IKN = KK + 1
+END IF
+!
+! INTERPOLATION
+! -- Along X
+! -- -- Alpha
+ZALPHAX = (PPOS(1) - 0.5*(XXHAT(IIP)+XXHAT(IIN))) / (XXHAT(IIN) - XXHAT(IIP))
+!!PRINT*, "ZALPHAX = ", ZALPHAX
+! -- -- -- Wind
+! -- -- Interpolated variable in temporary plane X
+ZUXNN = (1-ZALPHAX)*PVAR(IIP,IJN,IKN) + ZALPHAX*PVAR(IIN,IJN,IKN)
+ZUXNP = (1-ZALPHAX)*PVAR(IIP,IJN,IKP) + ZALPHAX*PVAR(IIN,IJN,IKP)
+ZUXPP = (1-ZALPHAX)*PVAR(IIP,IJP,IKP) + ZALPHAX*PVAR(IIN,IJP,IKP)
+ZUXPN = (1-ZALPHAX)*PVAR(IIP,IJP,IKN) + ZALPHAX*PVAR(IIN,IJP,IKN)
+! -- -- -- Height
+ZHXNN = (1-ZALPHAX)*PZH(IIP,IJN,IKN) + ZALPHAX*PZH(IIN,IJN,IKN)
+ZHXNP = (1-ZALPHAX)*PZH(IIP,IJN,IKP) + ZALPHAX*PZH(IIN,IJN,IKP)
+ZHXPP = (1-ZALPHAX)*PZH(IIP,IJP,IKP) + ZALPHAX*PZH(IIN,IJP,IKP)
+ZHXPN = (1-ZALPHAX)*PZH(IIP,IJP,IKN) + ZALPHAX*PZH(IIN,IJP,IKN)
+!
+!
+! -- Along Y
+! -- -- Alpha
+ZALPHAY = (PPOS(2) - 0.5*(XYHAT(IJP)+XYHAT(IJN))) / (XYHAT(IJN) - XYHAT(IJP))
+!PRINT*, "ZALPHAY = ", ZALPHAY
+! -- -- Interpolated variable in temporary plane Y
+! -- -- -- Wind
+ZUXN = (1-ZALPHAY)*ZUXPN + ZALPHAY*ZUXNN
+ZUXP = (1-ZALPHAY)*ZUXPP + ZALPHAY*ZUXNP
+! -- -- -- Height
+ZHXN = (1-ZALPHAY)*ZHXPN + ZALPHAY*ZHXNN
+ZHXP = (1-ZALPHAY)*ZHXPP + ZALPHAY*ZHXNP
+!
+!
+! -- Along Z
+! -- -- Alpha Z
+ZALPHAZ = (PPOS(3) - ZHXP) / (ZHXN - ZHXP)
+!PRINT*, "ZALPHAZ = ", ZALPHAZ
+ZUX = (1 - ZALPHAZ)*ZUXP + ZALPHAZ*ZUXN
+!
+!
+INTERP_LIN8NB = ZUX
+!
+!
+!
+END FUNCTION INTERP_LIN8NB
+!#########################################################
diff --git a/src/MNH/eol_printer.f90 b/src/MNH/eol_printer.f90
new file mode 100755
index 0000000000000000000000000000000000000000..c8b33e32ce1a107dfcb8eb03ce00342fa2425669
--- /dev/null
+++ b/src/MNH/eol_printer.f90
@@ -0,0 +1,356 @@
+!MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
+!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
+!MNH_LIC for details. version 1.
+! #########################
+ MODULE MODI_EOL_PRINTER
+! #########################
+!
+INTERFACE
+!
+! ****
+! ADNR
+! ****
+!
+SUBROUTINE PRINT_DATA_FARM_ADNR(KFILE,TPFARM)
+ USE MODD_EOL_ADNR, ONLY: FARM
+ INTEGER, INTENT(IN) :: KFILE ! output file
+ TYPE(FARM), INTENT(IN) :: TPFARM ! stored farm data
+END SUBROUTINE PRINT_DATA_FARM_ADNR
+!
+SUBROUTINE PRINT_DATA_TURBINE_ADNR(KFILE,TPTURBINE)
+ USE MODD_EOL_ADNR, ONLY : TURBINE
+ INTEGER, INTENT(IN) :: KFILE ! output file
+ TYPE(TURBINE), INTENT(IN) :: TPTURBINE ! stored turbine data
+END SUBROUTINE PRINT_DATA_TURBINE_ADNR
+!
+! ***
+! ALM
+! ***
+!
+SUBROUTINE PRINT_DATA_FARM_ALM(KFILE,TPFARM)
+ USE MODD_EOL_ALM, ONLY: FARM
+ INTEGER, INTENT(IN) :: KFILE ! output file
+ TYPE(FARM), INTENT(IN) :: TPFARM ! stored farm data
+END SUBROUTINE PRINT_DATA_FARM_ALM
+!
+SUBROUTINE PRINT_DATA_TURBINE_ALM(KFILE,TPTURBINE)
+ USE MODD_EOL_ALM, ONLY : TURBINE
+ INTEGER, INTENT(IN) :: KFILE ! output file
+ TYPE(TURBINE), INTENT(IN) :: TPTURBINE ! stored turbine data
+END SUBROUTINE PRINT_DATA_TURBINE_ALM
+!
+SUBROUTINE PRINT_DATA_BLADE_ALM(KFILE,TPBLADE)
+ USE MODD_EOL_ALM, ONLY : BLADE
+ INTEGER, INTENT(IN) :: KFILE ! output file
+ TYPE(BLADE), INTENT(IN) :: TPBLADE ! stored blade data
+END SUBROUTINE PRINT_DATA_BLADE_ALM
+!
+SUBROUTINE PRINT_DATA_AIRFOIL_ALM(KFILE,TPAIRFOIL)
+ USE MODD_EOL_ALM, ONLY : AIRFOIL
+ INTEGER, INTENT(IN) :: KFILE ! output file
+ TYPE(AIRFOIL), DIMENSION(:), INTENT(IN) :: TPAIRFOIL ! stored airfoil data
+END SUBROUTINE PRINT_DATA_AIRFOIL_ALM
+!
+SUBROUTINE OPEN_TECOUT(KFILE, KTCOUNT, KTSUBCOUNT)
+ INTEGER, INTENT(IN) :: KFILE ! File index
+ INTEGER, INTENT(IN) :: KTCOUNT ! Time step index
+ INTEGER, INTENT(IN) :: KTSUBCOUNT ! Subtime step index
+END SUBROUTINE OPEN_TECOUT
+!
+SUBROUTINE PRINT_TECOUT(KFILE,PVAR)
+ INTEGER, INTENT(IN) :: KFILE ! File index
+ REAL, DIMENSION(3), INTENT(IN) :: PVAR ! Vector to plot
+END SUBROUTINE PRINT_TECOUT
+!
+SUBROUTINE PRINT_TSPLIT(KNBSUBCOUNT,PTSUBSTEP)
+ INTEGER, INTENT(IN) :: KNBSUBCOUNT ! splitting value
+ REAL, INTENT(IN) :: PTSUBSTEP ! sub timestep
+END SUBROUTINE PRINT_TSPLIT
+!
+!
+END INTERFACE
+!
+END MODULE MODI_EOL_PRINTER
+!-------------------------------------------------------------------
+!
+!!**** *EOL_PRINT* -
+!!
+!! PURPOSE
+!! -------
+!! Some usefull subs to print wind turbine's datas
+!!
+!! AUTHOR
+!! ------
+!! PA. Joulin *CNRM & IFPEN*
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 26/10/2020
+!!
+!!---------------------------------------------------------------
+!
+!#########################################################
+SUBROUTINE PRINT_DATA_FARM_ADNR(KFILE,TPFARM)
+!
+USE MODD_EOL_ADNR, ONLY : FARM
+USE MODD_EOL_SHARED_IO, ONLY : CFARM_CSVDATA
+USE MODD_VAR_ll, ONLY : IP ! only master cpu
+!
+IMPLICIT NONE
+!
+INTEGER, INTENT(IN) :: KFILE ! File index
+TYPE(FARM), INTENT(IN) :: TPFARM ! dummy stored farm data
+!
+INTEGER :: JROT ! Loop index
+!
+IF (IP==1) THEN
+ WRITE(KFILE,*) ''
+ WRITE(KFILE,*) '======================== WIND TURBINE DATA ========================'
+ WRITE(KFILE,*) ''
+ WRITE(KFILE,*) '---- Farm ----'
+ WRITE(KFILE,*) 'Data from file : ', TRIM(CFARM_CSVDATA)
+ WRITE(KFILE,*) 'Number of turbines : ', TPFARM%NNB_TURBINES
+ WRITE(KFILE,*) 'Positions [m] and thrust coef [-] : '
+ DO JROT=1, TPFARM%NNB_TURBINES
+ WRITE(KFILE,'(1X,A,I3,A,F10.1,A,F10.1,A,F10.3)') 'n.', JROT,&
+ ' : X = ', TPFARM%XPOS_X(JROT),&
+ ' ; Y = ', TPFARM%XPOS_Y(JROT),&
+ ' ; CT_inf = ', TPFARM%XCT_INF(JROT)
+ END DO
+ WRITE(KFILE,*) ''
+END IF
+!
+END SUBROUTINE PRINT_DATA_FARM_ADNR
+!#########################################################
+!
+!#########################################################
+SUBROUTINE PRINT_DATA_TURBINE_ADNR(KFILE,TPTURBINE)
+!
+USE MODD_EOL_ADNR, ONLY : TURBINE
+USE MODD_EOL_SHARED_IO, ONLY : CTURBINE_CSVDATA
+USE MODD_VAR_ll, ONLY : IP ! only master cpu
+!
+IMPLICIT NONE
+!
+INTEGER, INTENT(IN) :: KFILE ! File index
+TYPE(TURBINE), INTENT(IN) :: TPTURBINE ! dummy stored turbine data
+!
+IF (IP==1) THEN
+ WRITE(KFILE,*) '---- Turbine ----'
+ WRITE(KFILE,* ) 'Data from file : ', TRIM(CTURBINE_CSVDATA)
+ WRITE(KFILE,'(1X,A,A10)' ) 'Wind turbine : ', TPTURBINE%CNAME
+ WRITE(KFILE,'(1X,A,F10.1)') 'Hub height [m] : ', TPTURBINE%XH_HEIGHT
+ WRITE(KFILE,'(1X,A,F10.3)') 'Blade radius [m] : ', TPTURBINE%XR_MAX
+ WRITE(KFILE,*) ''
+ WRITE(KFILE,*) '==================================================================='
+ WRITE(KFILE,*) ''
+END IF
+!
+END SUBROUTINE PRINT_DATA_TURBINE_ADNR
+!#########################################################
+!
+!#########################################################
+SUBROUTINE PRINT_DATA_FARM_ALM(KFILE,TPFARM)
+!
+USE MODD_EOL_ALM, ONLY : FARM
+USE MODD_EOL_SHARED_IO, ONLY : CFARM_CSVDATA
+USE MODD_VAR_ll, ONLY : IP ! only master cpu
+!
+IMPLICIT NONE
+!
+INTEGER, INTENT(IN) :: KFILE ! File index
+TYPE(FARM), INTENT(IN) :: TPFARM ! dummy stored farm data
+!
+INTEGER :: JROT ! Loop index
+!
+IF (IP==1) THEN
+ WRITE(KFILE,*) ''
+ WRITE(KFILE,*) '======================== WIND TURBINE DATA ========================'
+ WRITE(KFILE,*) ''
+ WRITE(KFILE,*) '---- Farm ----'
+ WRITE(KFILE,*) 'Data from file : ', TRIM(CFARM_CSVDATA)
+ WRITE(KFILE,*) 'Number of turbines : ', TPFARM%NNB_TURBINES
+ WRITE(KFILE,*) 'Tower base positions (X,Y) [m] : '
+ DO JROT=1, TPFARM%NNB_TURBINES
+ WRITE(KFILE, '(1X,A,I3,A,F10.1,A,F10.1,A)') 'n.', JROT,&
+ ': (', TPFARM%XPOS_X(JROT),',',TPFARM%XPOS_Y(JROT),')'
+ END DO
+ WRITE(KFILE,*) 'Working state (rad/s,rad,rad) : '
+ DO JROT=1, TPFARM%NNB_TURBINES
+ WRITE(KFILE, '(1X,A,I3,A,F10.5,A,F10.5,A,F10.5)') 'n.', JROT,&
+ ': Omega = ', TPFARM%XOMEGA(JROT), &
+ ' ; Yaw = ', TPFARM%XNAC_YAW(JROT),&
+ ' ; Pitch = ', TPFARM%XBLA_PITCH(JROT)
+ END DO
+ WRITE(KFILE,*) ''
+END IF
+!
+END SUBROUTINE PRINT_DATA_FARM_ALM
+!#########################################################
+!
+!#########################################################
+SUBROUTINE PRINT_DATA_TURBINE_ALM(KFILE,TPTURBINE)
+!
+USE MODD_EOL_ALM, ONLY : TURBINE
+USE MODD_EOL_SHARED_IO, ONLY : CTURBINE_CSVDATA
+USE MODD_VAR_ll, ONLY : IP ! only master cpu
+!
+IMPLICIT NONE
+!
+INTEGER, INTENT(IN) :: KFILE ! File index
+TYPE(TURBINE), INTENT(IN) :: TPTURBINE ! dummy stored turbine data
+!
+IF (IP==1) THEN
+ WRITE(KFILE,*) '---- Turbine ----'
+ WRITE(KFILE,* ) 'Data from file : ', TRIM(CTURBINE_CSVDATA)
+ WRITE(KFILE,'(1X,A,A10)' ) 'Wind turbine : ', TPTURBINE%CNAME
+ WRITE(KFILE,'(1X,A,I10)' ) 'Number of blades : ', TPTURBINE%NNB_BLADES
+ WRITE(KFILE,'(1X,A,F10.1)') 'Hub height [m] : ', TPTURBINE%XH_HEIGHT
+ WRITE(KFILE,'(1X,A,F10.3)') 'Blade min radius [m] : ', TPTURBINE%XR_MIN
+ WRITE(KFILE,'(1X,A,F10.3)') 'Blade max radius [m] : ', TPTURBINE%XR_MAX
+ WRITE(KFILE,'(1X,A,F10.3)') 'Nacelle tilt [rad] : ', TPTURBINE%XNAC_TILT
+ WRITE(KFILE,'(1X,A,F10.3)') 'Hub deport [m] : ', TPTURBINE%XH_DEPORT
+ WRITE(KFILE,*) ''
+END IF
+!
+END SUBROUTINE PRINT_DATA_TURBINE_ALM
+!#########################################################
+!
+!#########################################################
+SUBROUTINE PRINT_DATA_BLADE_ALM(KFILE,TPBLADE)
+!
+USE MODD_EOL_ALM, ONLY : BLADE
+USE MODD_EOL_SHARED_IO, ONLY : CBLADE_CSVDATA
+USE MODD_VAR_ll, ONLY : IP ! only master cpu
+!
+IMPLICIT NONE
+!
+INTEGER, INTENT(IN) :: KFILE ! File index
+TYPE(BLADE), INTENT(IN) :: TPBLADE ! dummy stored blade data
+!
+IF (IP==1) THEN
+ WRITE(KFILE,*) '---- Blade ----'
+ WRITE(KFILE,* ) 'Data from file : ', TRIM(CBLADE_CSVDATA)
+ WRITE(KFILE,'(1X,A,I10)' ) 'Nb of data (from data file) : ', TPBLADE%NNB_BLADAT
+ WRITE(KFILE,'(1X,A,F10.1)') 'First node radius [m] : ', TPBLADE%XRAD(1)
+ WRITE(KFILE,'(1X,A,F10.1)') 'Last node radius [m] : ', TPBLADE%XRAD(TPBLADE%NNB_BLADAT)
+ WRITE(KFILE,'(1X,A,F10.1)') 'Chord max. [m] : ', MAXVAL(TPBLADE%XCHORD(:))
+ WRITE(KFILE,'(1X,A,I10)' ) 'Nb of blade element (from nam) : ', TPBLADE%NNB_BLAELT
+ WRITE(KFILE,*) ''
+END IF
+!
+END SUBROUTINE PRINT_DATA_BLADE_ALM
+!#########################################################
+!
+!#########################################################
+SUBROUTINE PRINT_DATA_AIRFOIL_ALM(KFILE,TPAIRFOIL)
+!
+USE MODD_EOL_ALM, ONLY : AIRFOIL
+USE MODD_EOL_SHARED_IO, ONLY : CAIRFOIL_CSVDATA
+USE MODD_VAR_ll, ONLY : IP ! only master cpu
+!
+IMPLICIT NONE
+!
+INTEGER, INTENT(IN) :: KFILE ! File index
+TYPE(AIRFOIL), DIMENSION(:), INTENT(IN) :: TPAIRFOIL ! dummy stored airfoil data
+!
+INTEGER :: JA
+!
+IF (IP==1) THEN
+ WRITE(KFILE,*) '---- Airfoils ----'
+ WRITE(KFILE,* ) 'Data from file : ', TRIM(CAIRFOIL_CSVDATA)
+ WRITE(KFILE,'(1X,A,I10)' ) 'Nb of airfoils (from data file) : ', SIZE(TPAIRFOIL)
+ WRITE(KFILE,'(1X,A)' ) 'Different airfoils : '
+ DO JA=1,SIZE(TPAIRFOIL)
+ WRITE(KFILE,'(1X,A,I3,A,A)') 'Airfoil n.', JA,&
+ ': ', TPAIRFOIL(JA)%CNAME
+ END DO
+ WRITE(KFILE,*) ''
+ WRITE(KFILE,*) '==================================================================='
+ WRITE(KFILE,*) ''
+END IF
+!
+END SUBROUTINE PRINT_DATA_AIRFOIL_ALM
+!#########################################################
+!
+!
+!#########################################################
+SUBROUTINE OPEN_TECOUT(KFILE, KTCOUNT, KTSUBCOUNT)
+!
+USE MODD_EOL_ALM, ONLY:TFARM,TTURBINE,TBLADE
+!
+IMPLICIT NONE
+!
+INTEGER, INTENT(IN) :: KFILE ! File index
+INTEGER, INTENT(IN) :: KTCOUNT ! Time step index
+INTEGER, INTENT(IN) :: KTSUBCOUNT ! Subtime step index
+!
+INTEGER :: INB_WT, INB_B, INB_BELT ! Total numbers of wind turbines, blades, and blade elt
+INTEGER :: INB_TELT, INB_NELT ! Total numbers of tower elt, and nacelle elt
+INTEGER :: ITOTELT ! Total number of points
+!
+CHARACTER(LEN=1024) :: HFILE ! File name
+!
+INB_WT = TFARM%NNB_TURBINES
+INB_B = TTURBINE%NNB_BLADES
+INB_BELT = TBLADE%NNB_BLAELT
+! Hard coded variables, but they will be usefull in next updates
+INB_TELT = 2
+INB_NELT = 2
+!
+ITOTELT = INB_WT*(INB_TELT+INB_NELT+INB_B*(1+INB_BELT*3))
+!
+! File name and opening
+WRITE(HFILE, "(A18,I4.4,I2.2,A3)") "Tecplot2.0_Output_", KTCOUNT, KTSUBCOUNT,".tp"
+OPEN( KFILE, file=HFILE, form="FORMATTED")
+!
+! Tecplot Header
+WRITE(KFILE,*) 'TITLE="Wind Turbines Points"'
+WRITE(KFILE,*) 'VARIABLES="X" "Y" "Z"'
+WRITE(KFILE,*) 'ZONE I=',ITOTELT,' J=3 K=1 DATAPACKING=POINT'
+!
+END SUBROUTINE OPEN_TECOUT
+!#########################################################
+!
+!#########################################################
+SUBROUTINE PRINT_TECOUT(KFILE,PVAR)
+IMPLICIT NONE
+INTEGER, INTENT(IN) :: KFILE ! File index
+REAL, DIMENSION(3), INTENT(IN) :: PVAR ! Vector to plot
+!
+! It plots two points, slightly different, to get a thickness
+!
+WRITE(KFILE,*) PVAR(1), &
+ PVAR(2), &
+ PVAR(3)
+!
+END SUBROUTINE PRINT_TECOUT
+!#########################################################
+!
+!#########################################################
+SUBROUTINE PRINT_TSPLIT(KNBSUBCOUNT,PTSUBSTEP)
+USE MODD_LUNIT_n , ONLY: TLUOUT
+!
+INTEGER, INTENT(IN) :: KNBSUBCOUNT ! splitting value
+REAL, INTENT(IN) :: PTSUBSTEP ! sub-time step
+!
+WRITE(TLUOUT%NLU,'(A)') 'From EOL - Actuator Line Model. Time-splitting is activated:'
+WRITE(TLUOUT%NLU,'(1X,A,I2)') 'Number of splitting: ', KNBSUBCOUNT
+WRITE(TLUOUT%NLU,'(1X,A,F6.4)') 'Sub-time step value: ', PTSUBSTEP
+!
+END SUBROUTINE PRINT_TSPLIT
+!#########################################################
+!
+!#########################################################
+SUBROUTINE PRINT_ERROR_WTCFL(KFILE,PMAXTSTEP)
+INTEGER, INTENT(IN) :: KFILE ! output file
+REAL, INTENT(IN) :: PMAXTSTEP ! maximum acceptable time-step
+WRITE(KFILE,'(A,A,A,A,A,F10.8,A)') &
+ 'Sorry but I had to stop the simulation. ' ,&
+ 'The time step XTSTEP is too high: ' ,&
+ 'the blades can jump over one or several cells. ' ,&
+ 'Please, turn on the time-splitting method (LTIMESPLIT=.TRUE.), ',&
+ 'or decrease XTSTEP to a value lower than ', PMAXTSTEP, ' sec.'
+END SUBROUTINE PRINT_ERROR_WTCFL
+!
diff --git a/src/MNH/eol_reader.f90 b/src/MNH/eol_reader.f90
new file mode 100755
index 0000000000000000000000000000000000000000..9dfa91fbf2531ad8cc505eca26da39f0c923cb77
--- /dev/null
+++ b/src/MNH/eol_reader.f90
@@ -0,0 +1,676 @@
+!MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
+!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
+!MNH_LIC for details. version 1.
+! #######################
+ MODULE MODI_EOL_READER
+! #######################
+!
+INTERFACE
+!
+! ADNR
+!
+SUBROUTINE READ_CSVDATA_FARM_ADNR(KLUNAM,HFILE,TPFARM)
+ USE MODD_EOL_ADNR, ONLY: FARM
+ INTEGER, INTENT(IN) :: KLUNAM ! logical unit of the file
+ CHARACTER(LEN=*), INTENT(IN) :: HFILE ! file to read
+ TYPE(FARM), INTENT(OUT) :: TPFARM ! stored farm data
+END SUBROUTINE READ_CSVDATA_FARM_ADNR
+!
+SUBROUTINE READ_CSVDATA_TURBINE_ADNR(KLUNAM,HFILE,TPTURBINE)
+ USE MODD_EOL_ADNR, ONLY : TURBINE
+ INTEGER, INTENT(IN) :: KLUNAM ! logical unit of the file
+ CHARACTER(LEN=*), INTENT(IN) :: HFILE ! file to read
+ TYPE(TURBINE), INTENT(OUT) :: TPTURBINE ! stored turbine data
+END SUBROUTINE READ_CSVDATA_TURBINE_ADNR
+!
+! ALM
+!
+SUBROUTINE READ_CSVDATA_FARM_ALM(KLUNAM,HFILE,TPFARM)
+ USE MODD_EOL_ALM, ONLY: FARM
+ INTEGER, INTENT(IN) :: KLUNAM ! logical unit of the file
+ CHARACTER(LEN=*), INTENT(IN) :: HFILE ! file to read
+ TYPE(FARM), INTENT(OUT) :: TPFARM ! stored farm data
+END SUBROUTINE READ_CSVDATA_FARM_ALM
+!
+SUBROUTINE READ_CSVDATA_TURBINE_ALM(KLUNAM,HFILE,TPTURBINE)
+ USE MODD_EOL_ALM, ONLY : TURBINE
+ INTEGER, INTENT(IN) :: KLUNAM ! logical unit of the file
+ CHARACTER(LEN=*), INTENT(IN) :: HFILE ! file to read
+ TYPE(TURBINE), INTENT(OUT) :: TPTURBINE ! stored turbine data
+END SUBROUTINE READ_CSVDATA_TURBINE_ALM
+
+SUBROUTINE READ_CSVDATA_BLADE_ALM(KLUNAM,HFILE,TPTURBINE,TPBLADE)
+ USE MODD_EOL_ALM, ONLY : TURBINE, BLADE
+ INTEGER, INTENT(IN) :: KLUNAM ! logical unit of the file
+ CHARACTER(LEN=*), INTENT(IN) :: HFILE ! file to read
+ TYPE(TURBINE), INTENT(IN) :: TPTURBINE ! stored turbine data
+ TYPE(BLADE), INTENT(OUT) :: TPBLADE ! stored blade data
+END SUBROUTINE READ_CSVDATA_BLADE_ALM
+!
+SUBROUTINE READ_CSVDATA_AIRFOIL_ALM(KLUNAM,HFILE,TPBLADE,TPAIRFOIL)
+ USE MODD_EOL_ALM, ONLY : BLADE, AIRFOIL
+ INTEGER, INTENT(IN) :: KLUNAM ! logical unit of the file
+ CHARACTER(LEN=*), INTENT(IN) :: HFILE ! file to read
+ TYPE(BLADE), INTENT(IN) :: TPBLADE ! stored blade data (to select airfoils)
+ TYPE(AIRFOIL), DIMENSION(:), ALLOCATABLE, INTENT(OUT) :: TPAIRFOIL ! stored airfoil data
+END SUBROUTINE READ_CSVDATA_AIRFOIL_ALM
+!
+SUBROUTINE HOW_MANY_LINES_OF(KLUNAM,HFILE,HNAME,KLINE)
+ INTEGER, INTENT(IN) :: KLUNAM ! logical unit of the file
+ CHARACTER(LEN=*), INTENT(IN) :: HFILE ! file to read
+ CHARACTER(LEN=*), INTENT(IN) :: HNAME ! turbine's name
+ INTEGER, INTENT(OUT) :: KLINE
+END SUBROUTINE HOW_MANY_LINES_OF
+!
+FUNCTION GET_AIRFOIL_ID(TPTURBINE,TPBLADE,TPAIRFOIL,PRADIUS)
+ USE MODD_EOL_ALM, ONLY : TURBINE, BLADE, AIRFOIL
+ IMPLICIT NONE
+ INTEGER :: GET_AIRFOIL_ID
+ TYPE(TURBINE), INTENT(IN) :: TPTURBINE ! stored turbine data
+ TYPE(BLADE), INTENT(IN) :: TPBLADE ! stored blade data
+ TYPE(AIRFOIL), DIMENSION(:), INTENT(IN) :: TPAIRFOIL ! stored airfoil data
+ REAL, INTENT(IN) :: PRADIUS ! Radius position studied
+END FUNCTION GET_AIRFOIL_ID
+!
+END INTERFACE
+!
+END MODULE MODI_EOL_READER
+!-------------------------------------------------------------------
+!
+!!**** *EOL_READER* -
+!!
+!! PURPOSE
+!! -------
+!! Some usefull subs to read wind turbine's datas
+!!
+!! AUTHOR
+!! ------
+!! PA. Joulin *CNRM & IFPEN*
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 05/2018
+!! Modification 21/10/20 (PA. Joulin) Updated for a main version
+!!
+!!---------------------------------------------------------------
+!
+!#########################################################
+SUBROUTINE READ_CSVDATA_FARM_ADNR(KLUNAM,HFILE,TPFARM)
+!
+USE MODD_EOL_ADNR, ONLY: FARM
+USE MODI_EOL_ERROR, ONLY: EOL_CSVNOTFOUND_ERROR, EOL_CSVEMPTY_ERROR
+!
+IMPLICIT NONE
+!
+INTEGER, INTENT(IN) :: KLUNAM ! logical unit of the file
+CHARACTER(LEN=*), INTENT(IN) :: HFILE ! file to read
+TYPE(FARM), INTENT(OUT) :: TPFARM ! dummy stored data blade
+!
+LOGICAL :: GEXIST ! Existence of file
+!
+INTEGER :: INBLINE ! Nb of line in csv file
+!
+CHARACTER(LEN=400) :: YSTRING
+!
+! Read data
+REAL :: ZPOS_X
+REAL :: ZPOS_Y
+REAL :: ZCT_INF
+!
+! Checking file existence
+INQUIRE(FILE=HFILE, EXIST=GEXIST)
+IF (.NOT.GEXIST) THEN
+ CALL EOL_CSVNOTFOUND_ERROR(HFILE)
+END IF
+!
+! Opening the file
+OPEN(UNIT=KLUNAM,FILE=HFILE, FORM='formatted', STATUS='OLD')
+! Counting number of line
+REWIND(KLUNAM)
+INBLINE=0
+DO
+ READ(KLUNAM,END=101,FMT='(A400)') YSTRING
+ IF (LEN_TRIM(YSTRING) > 0) THEN
+ INBLINE = INBLINE + 1
+ END IF
+END DO
+!
+101 CONTINUE
+IF (INBLINE < 2) THEN
+ CALL EOL_CSVEMPTY_ERROR(HFILE,INBLINE)
+ STOP
+ELSE
+ ! Saving number of wind turbine
+ TPFARM%NNB_TURBINES = INBLINE - 1
+ ! Allocations
+ ALLOCATE(TPFARM%XPOS_X(TPFARM%NNB_TURBINES))
+ ALLOCATE(TPFARM%XPOS_Y(TPFARM%NNB_TURBINES))
+ ALLOCATE(TPFARM%XCT_INF(TPFARM%NNB_TURBINES))
+ !
+ ! New read
+ REWIND(KLUNAM)
+ READ(KLUNAM,FMT='(A400)') YSTRING ! Header reading
+ !
+ ! Saving data
+ DO INBLINE=1, TPFARM%NNB_TURBINES
+ READ(KLUNAM,FMT='(A400)') YSTRING
+ READ(YSTRING,*) ZPOS_X, ZPOS_Y, ZCT_INF
+ TPFARM%XPOS_X(INBLINE) = ZPOS_X
+ TPFARM%XPOS_Y(INBLINE) = ZPOS_Y
+ TPFARM%XCT_INF(INBLINE) = ZCT_INF
+ END DO
+ CLOSE(KLUNAM)
+ RETURN
+END IF
+!
+END SUBROUTINE READ_CSVDATA_FARM_ADNR
+!#########################################################
+!
+!#########################################################
+SUBROUTINE READ_CSVDATA_TURBINE_ADNR(KLUNAM,HFILE,TPTURBINE)
+!
+USE MODD_EOL_ADNR, ONLY: TURBINE
+USE MODI_EOL_ERROR, ONLY: EOL_CSVNOTFOUND_ERROR, EOL_CSVEMPTY_ERROR
+!
+IMPLICIT NONE
+!
+INTEGER, INTENT(IN) :: KLUNAM ! logical unit of the file
+CHARACTER(LEN=*), INTENT(IN) :: HFILE ! file to read
+TYPE(TURBINE), INTENT(OUT) :: TPTURBINE ! dummy stored data turbine
+!
+LOGICAL :: GEXIST ! Existence of file
+!
+INTEGER :: INBLINE ! Nb of line in csv file
+!
+CHARACTER(LEN=400) :: YSTRING
+!
+CHARACTER(LEN=80) :: YWT_NAME
+REAL :: ZH_HEIGHT
+REAL :: ZR_MAX
+!
+! Checking file existence
+INQUIRE(FILE=HFILE, EXIST=GEXIST)
+IF (.NOT.GEXIST) THEN
+ CALL EOL_CSVNOTFOUND_ERROR(HFILE)
+END IF
+!
+! Opening
+OPEN(UNIT=KLUNAM,FILE=HFILE, FORM='formatted', STATUS='OLD')
+!
+! Counting number of line
+REWIND(KLUNAM)
+INBLINE=0
+DO
+ READ(KLUNAM,END=101,FMT='(A400)') YSTRING
+ IF (LEN_TRIM(YSTRING) > 0) THEN
+ INBLINE = INBLINE + 1
+ END IF
+END DO
+!
+101 CONTINUE
+IF (INBLINE /= 2) THEN
+ CALL EOL_CSVEMPTY_ERROR(HFILE,INBLINE)
+ STOP
+ELSE
+ REWIND(KLUNAM)
+ READ(KLUNAM,FMT='(A400)') YSTRING ! Header reading
+ READ(KLUNAM,FMT='(A400)') YSTRING ! Reading next line
+ ! Read data
+ READ(YSTRING,*) YWT_NAME, ZH_HEIGHT, ZR_MAX ! reading data
+ TPTURBINE%CNAME = YWT_NAME ! Saving them
+ TPTURBINE%XH_HEIGHT = ZH_HEIGHT
+ TPTURBINE%XR_MAX = ZR_MAX
+ REWIND(KLUNAM) ! Rembobinage, plutôt 2 fois qu'1 !
+ RETURN
+ CLOSE(KLUNAM)
+END IF
+!
+END SUBROUTINE READ_CSVDATA_TURBINE_ADNR
+!#########################################################
+!
+!#########################################################
+SUBROUTINE READ_CSVDATA_FARM_ALM(KLUNAM,HFILE,TPFARM)
+!
+USE MODD_EOL_ALM, ONLY: FARM
+USE MODI_EOL_ERROR, ONLY: EOL_CSVNOTFOUND_ERROR, EOL_CSVEMPTY_ERROR
+!
+IMPLICIT NONE
+!
+INTEGER, INTENT(IN) :: KLUNAM ! logical unit of the file
+CHARACTER(LEN=*), INTENT(IN) :: HFILE ! file to read
+TYPE(FARM), INTENT(OUT) :: TPFARM ! dummy stored data blade
+!
+LOGICAL :: GEXIST ! Existence of file
+!
+INTEGER :: INBLINE ! Nb of line in csv file
+!
+CHARACTER(LEN=400) :: YSTRING
+!
+! Read data
+REAL :: ZPOS_X
+REAL :: ZPOS_Y
+REAL :: ZOMEGA
+REAL :: ZYAW
+REAL :: ZPITCH
+!
+! Checking file existence
+INQUIRE(FILE=HFILE, EXIST=GEXIST)
+IF (.NOT.GEXIST) THEN
+ CALL EOL_CSVNOTFOUND_ERROR(HFILE)
+END IF
+!
+! Opening the file
+OPEN(UNIT=KLUNAM,FILE=HFILE, FORM='formatted', STATUS='OLD')
+! Counting number of line
+REWIND(KLUNAM)
+INBLINE=0
+DO
+ READ(KLUNAM,END=101,FMT='(A400)') YSTRING
+ IF (LEN_TRIM(YSTRING) > 0) THEN
+ INBLINE = INBLINE + 1
+ END IF
+END DO
+!
+101 CONTINUE
+IF (INBLINE < 2) THEN
+ CALL EOL_CSVEMPTY_ERROR(HFILE,INBLINE)
+ STOP
+ELSE
+ ! Saving number of wind turbine
+ TPFARM%NNB_TURBINES = INBLINE - 1
+ ! Allocations
+ ALLOCATE(TPFARM%XPOS_X (TPFARM%NNB_TURBINES))
+ ALLOCATE(TPFARM%XPOS_Y (TPFARM%NNB_TURBINES))
+ ALLOCATE(TPFARM%XOMEGA (TPFARM%NNB_TURBINES))
+ ALLOCATE(TPFARM%XNAC_YAW (TPFARM%NNB_TURBINES))
+ ALLOCATE(TPFARM%XBLA_PITCH(TPFARM%NNB_TURBINES))
+ !
+ ! New read
+ REWIND(KLUNAM)
+ READ(KLUNAM,FMT='(A400)') YSTRING ! Header reading
+ !
+ ! Saving data
+ DO INBLINE=1, TPFARM%NNB_TURBINES
+ READ(KLUNAM,FMT='(A400)') YSTRING
+ READ(YSTRING,*) ZPOS_X, ZPOS_Y, ZOMEGA, ZYAW, ZPITCH
+ TPFARM%XPOS_X(INBLINE) = ZPOS_X
+ TPFARM%XPOS_Y(INBLINE) = ZPOS_Y
+ TPFARM%XOMEGA(INBLINE) = ZOMEGA
+ TPFARM%XNAC_YAW(INBLINE) = ZYAW
+ TPFARM%XBLA_PITCH(INBLINE) = ZPITCH
+ END DO
+ CLOSE(KLUNAM)
+ RETURN
+END IF
+!
+END SUBROUTINE READ_CSVDATA_FARM_ALM
+!#########################################################
+!
+!#########################################################
+SUBROUTINE READ_CSVDATA_TURBINE_ALM(KLUNAM,HFILE,TPTURBINE)
+!
+USE MODD_EOL_ALM, ONLY: TURBINE
+USE MODI_EOL_ERROR, ONLY: EOL_CSVNOTFOUND_ERROR, EOL_CSVEMPTY_ERROR
+!
+IMPLICIT NONE
+!
+INTEGER, INTENT(IN) :: KLUNAM ! logical unit of the file
+CHARACTER(LEN=*), INTENT(IN) :: HFILE ! file to read
+TYPE(TURBINE), INTENT(OUT) :: TPTURBINE ! dummy stored data turbine
+!
+LOGICAL :: GEXIST ! Existence of file
+!
+INTEGER :: INBLINE ! Nb of line in csv file
+!
+CHARACTER(LEN=400) :: YSTRING
+!
+CHARACTER(LEN=80) :: YWT_NAME
+INTEGER :: INB_BLADE
+REAL :: ZH_HEIGHT
+REAL :: ZR_MIN
+REAL :: ZR_MAX
+REAL :: ZNAC_TILT
+REAL :: ZHUB_DEPORT
+!
+! Checking file existence
+INQUIRE(FILE=HFILE, EXIST=GEXIST)
+IF (.NOT.GEXIST) THEN
+ CALL EOL_CSVNOTFOUND_ERROR(HFILE)
+END IF
+!
+! Opening
+OPEN(UNIT=KLUNAM,FILE=HFILE, FORM='formatted', STATUS='OLD')
+!
+! Counting number of line
+REWIND(KLUNAM)
+INBLINE=0
+DO
+ READ(KLUNAM,END=101,FMT='(A400)') YSTRING
+ IF (LEN_TRIM(YSTRING) > 0) THEN
+ INBLINE = INBLINE + 1
+ END IF
+END DO
+!
+101 CONTINUE
+IF (INBLINE /= 2) THEN
+ CALL EOL_CSVEMPTY_ERROR(HFILE,INBLINE)
+ STOP
+ELSE
+ REWIND(KLUNAM)
+ READ(KLUNAM,FMT='(A400)') YSTRING ! Header reading
+ READ(KLUNAM,FMT='(A400)') YSTRING ! Reading next line
+ ! Read data
+ READ(YSTRING,*) YWT_NAME, INB_BLADE, ZH_HEIGHT,& ! reading data
+ ZR_MIN, ZR_MAX, ZNAC_TILT, &
+ ZHUB_DEPORT
+ TPTURBINE%CNAME = YWT_NAME ! Saving them
+ TPTURBINE%NNB_BLADES = INB_BLADE
+ TPTURBINE%XH_HEIGHT = ZH_HEIGHT
+ TPTURBINE%XR_MIN = ZR_MIN
+ TPTURBINE%XR_MAX = ZR_MAX
+ TPTURBINE%XNAC_TILT = ZNAC_TILT
+ TPTURBINE%XH_DEPORT = ZHUB_DEPORT
+ REWIND(KLUNAM)
+ RETURN
+ CLOSE(KLUNAM)
+END IF
+!
+END SUBROUTINE READ_CSVDATA_TURBINE_ALM
+!#########################################################
+!
+!#########################################################
+SUBROUTINE READ_CSVDATA_BLADE_ALM(KLUNAM,HFILE,TPTURBINE,TPBLADE)
+!
+USE MODD_EOL_ALM, ONLY: TURBINE, BLADE, NNB_BLAELT
+USE MODI_EOL_ERROR, ONLY: EOL_CSVNOTFOUND_ERROR, EOL_CSVEMPTY_ERROR
+USE MODI_EOL_ERROR, ONLY: EOL_BLADEDATA_ERROR
+!
+INTEGER, INTENT(IN) :: KLUNAM ! logical unit of the file
+CHARACTER(LEN=*), INTENT(IN) :: HFILE ! file to read
+TYPE(TURBINE), INTENT(IN) :: TPTURBINE ! stored turbine data
+TYPE(BLADE), INTENT(OUT) :: TPBLADE ! dummy stored data blade
+!
+LOGICAL :: GEXIST ! Existence of file
+!
+INTEGER :: INBLINE ! Nb of line in csv file
+INTEGER :: INBDATA ! Nb of data (line/section) of blade
+!
+CHARACTER(LEN=400) :: YSTRING
+!
+REAL :: ZCENTER ! Center pos. of elmt [m]
+REAL :: ZCHORD ! Blade chord of elmt [m]
+REAL :: ZTWIST ! Twist of elmt [rad]
+CHARACTER(LEN=20) :: YAIRFOIL ! Airfoil name [-]
+!
+! Checking file existence
+INQUIRE(FILE=HFILE, EXIST=GEXIST)
+IF (.NOT.GEXIST) THEN
+ CALL EOL_CSVNOTFOUND_ERROR(HFILE)
+END IF
+!
+! Ouverture
+OPEN(UNIT=KLUNAM,FILE=HFILE, FORM='formatted', STATUS='OLD')
+! Counting number of line
+REWIND(KLUNAM)
+INBLINE=0
+DO
+ READ(KLUNAM,END=101,FMT='(A400)') YSTRING
+ IF (LEN_TRIM(YSTRING) > 0) THEN
+ INBLINE = INBLINE + 1
+ END IF
+END DO
+!
+101 CONTINUE
+IF (INBLINE < 2) THEN
+ CALL EOL_CSVEMPTY_ERROR(HFILE,INBLINE)
+ STOP
+ELSE
+ TPBLADE%NNB_BLAELT = NNB_BLAELT
+ ! Saving number of data
+ TPBLADE%NNB_BLADAT = INBLINE - 1
+ ALLOCATE(TPBLADE%XRAD(TPBLADE%NNB_BLADAT))
+ ALLOCATE(TPBLADE%XCHORD(TPBLADE%NNB_BLADAT))
+ ALLOCATE(TPBLADE%XTWIST(TPBLADE%NNB_BLADAT))
+ ALLOCATE(TPBLADE%CAIRFOIL(TPBLADE%NNB_BLADAT))
+ !
+ ! New read
+ REWIND(KLUNAM)
+ READ(KLUNAM,FMT='(A400)') YSTRING ! Header reading
+ !
+ ! Saving data
+ DO INBLINE=1, TPBLADE%NNB_BLADAT
+ READ(KLUNAM,FMT='(A400)') YSTRING
+ READ(YSTRING,*) ZCENTER, ZCHORD, ZTWIST, YAIRFOIL ! Reading data
+ IF ((ZCENTER<=0.0) .OR. (ZCENTER>= 1.0)) THEN
+ ! Checking data
+ CALL EOL_BLADEDATA_ERROR(ZCENTER)
+ ELSE
+ ! Storing them
+ TPBLADE%XRAD(INBLINE) = ZCENTER*(TPTURBINE%XR_MAX-TPTURBINE%XR_MIN) & ! Data in %
+ + TPTURBINE%XR_MIN
+ TPBLADE%XCHORD(INBLINE) = ZCHORD
+ TPBLADE%XTWIST(INBLINE) = ZTWIST
+ TPBLADE%CAIRFOIL(INBLINE) = YAIRFOIL
+ END IF
+ END DO
+ CLOSE(KLUNAM)
+ RETURN
+END IF
+!
+END SUBROUTINE READ_CSVDATA_BLADE_ALM
+!#########################################################
+!
+!#########################################################
+SUBROUTINE READ_CSVDATA_AIRFOIL_ALM(KLUNAM,HFILE,TPBLADE,TPAIRFOIL)
+!
+USE MODD_EOL_ALM, ONLY: BLADE, AIRFOIL
+USE MODI_EOL_ERROR, ONLY: EOL_CSVNOTFOUND_ERROR, EOL_CSVEMPTY_ERROR
+USE MODI_EOL_ERROR, ONLY: EOL_AIRFOILNOTFOUND_ERROR
+!
+INTEGER, INTENT(IN) :: KLUNAM ! logical unit of the file
+CHARACTER(LEN=*), INTENT(IN) :: HFILE ! file to read
+TYPE(BLADE), INTENT(IN) :: TPBLADE ! stored blade data (to select airfoils)
+TYPE(AIRFOIL), DIMENSION(:), ALLOCATABLE, INTENT(OUT) :: TPAIRFOIL ! dummy stored data blade
+!
+LOGICAL :: GEXIST ! Existence of file
+!
+INTEGER :: INBDATA ! Nb of data (line/section) per airfoil
+INTEGER :: INBLINE ! Nb of line in csv file
+LOGICAL :: GAIRFLAG ! Flag for airfoil counting
+!
+INTEGER :: JI, JJ, IA ! loop control
+INTEGER :: INBAIRFOIL ! Nb of differents airfoils on one blade
+!
+CHARACTER(LEN=400) :: YSTRING
+!
+CHARACTER(LEN=15), DIMENSION(:), ALLOCATABLE :: YAIRFOIL
+!
+CHARACTER(LEN=15) :: YREAD_NAME
+REAL :: ZAA ! Attack Angle [rad]
+REAL :: ZRE ! Reynolds Number [-]
+REAL :: ZCL ! Lift Coef [-]
+REAL :: ZCD ! Drag Coef [-]
+REAL :: ZCM ! Moment Coef [-]
+!
+! Checking file existence
+INQUIRE(FILE=HFILE, EXIST=GEXIST)
+IF (.NOT.GEXIST) THEN
+ CALL EOL_CSVNOTFOUND_ERROR(HFILE)
+END IF
+!
+! Ouverture
+OPEN(UNIT=KLUNAM,FILE=HFILE, FORM='formatted', STATUS='OLD')
+!
+! 1. Counting number of differents airfoils along the blade and selection :
+!
+! Allcation of local airfoil array
+ALLOCATE(YAIRFOIL(SIZE(TPBLADE%CAIRFOIL)))
+!
+YAIRFOIL(:) = ''
+INBAIRFOIL = 1
+YAIRFOIL(1) = TRIM(TPBLADE%CAIRFOIL(1))
+!
+DO JI=1, SIZE(TPBLADE%CAIRFOIL)
+ GAIRFLAG = .FALSE.
+ DO JJ=1, INBAIRFOIL
+ IF (TRIM(TPBLADE%CAIRFOIL(JI)) == TRIM(YAIRFOIL(JJ))) THEN
+ GAIRFLAG = .TRUE.
+ END IF
+ END DO
+ IF (GAIRFLAG .EQV. .FALSE.) THEN
+ INBAIRFOIL = INBAIRFOIL + 1
+ YAIRFOIL(INBAIRFOIL) = TRIM(TPBLADE%CAIRFOIL(JI))
+ END IF
+END DO
+ALLOCATE(TPAIRFOIL(INBAIRFOIL))
+!
+! 2. Reading and storing data :
+!
+DO IA = 1, INBAIRFOIL
+ ! Array allocation
+ CALL HOW_MANY_LINES_OF(KLUNAM,HFILE,YAIRFOIL(IA),INBDATA)
+ ALLOCATE(TPAIRFOIL(IA)%XAA(INBDATA))
+ ALLOCATE(TPAIRFOIL(IA)%XRE(INBDATA))
+ ALLOCATE(TPAIRFOIL(IA)%XCL(INBDATA))
+ ALLOCATE(TPAIRFOIL(IA)%XCD(INBDATA))
+ ALLOCATE(TPAIRFOIL(IA)%XCM(INBDATA))
+ !
+ REWIND(KLUNAM)
+ INBLINE = 0
+ DO
+ READ(KLUNAM,END=101,FMT='(A400)') YSTRING ! Header
+ !* reads the string
+ IF (LEN_TRIM(YSTRING)>0) THEN
+ READ(YSTRING,FMT=*) YREAD_NAME
+ IF (TRIM(YREAD_NAME)==TRIM(YAIRFOIL(IA))) THEN ! Read data
+ INBLINE = INBLINE + 1
+ READ(YSTRING,*) YREAD_NAME, ZAA, ZRE, &
+ ZCL, ZCD, ZCM
+ ! Storing data
+ TPAIRFOIL(IA)%CNAME = YREAD_NAME
+ TPAIRFOIL(IA)%XAA(INBLINE) = ZAA
+ TPAIRFOIL(IA)%XRE(INBLINE) = ZRE
+ TPAIRFOIL(IA)%XCL(INBLINE) = ZCL
+ TPAIRFOIL(IA)%XCD(INBLINE) = ZCD
+ TPAIRFOIL(IA)%XCM(INBLINE) = ZCM
+ ELSE ! The name doesnt appear during a new read..
+ IF (INBLINE > 0) THEN ! .. but it has already been found, ..
+ REWIND(KLUNAM) ! .. so it is the end of the data ..
+ EXIT ! .. and we can exit :)
+ END IF
+ END IF
+ END IF
+ END DO
+END DO
+!
+CLOSE(KLUNAM)
+101 CONTINUE
+ IF (INBLINE == 0) THEN
+ CALL EOL_AIRFOILNOTFOUND_ERROR(HFILE,YAIRFOIL(IA))
+ STOP
+ END IF
+END SUBROUTINE READ_CSVDATA_AIRFOIL_ALM
+!#########################################################
+!
+!#########################################################
+SUBROUTINE HOW_MANY_LINES_OF(KLUNAM,HFILE,HNAME,KLINE)
+!
+USE MODI_EOL_ERROR, ONLY: EOL_AIRFOILNOTFOUND_ERROR
+!
+IMPLICIT NONE
+!
+INTEGER, INTENT(IN) :: KLUNAM ! logical unit of the file
+CHARACTER(LEN=*), INTENT(IN) :: HFILE ! file to read
+CHARACTER(LEN=*), INTENT(IN) :: HNAME ! turbine's name
+INTEGER, INTENT(OUT) :: KLINE
+!
+!
+CHARACTER(LEN=400) :: YSTRING
+CHARACTER(LEN=80) :: HREAD_NAME
+!
+REWIND(KLUNAM)
+KLINE=0
+DO
+ READ(KLUNAM,END=101,FMT='(A400)') YSTRING
+!* reads the string
+ IF (LEN_TRIM(YSTRING) > 0) THEN
+ READ(YSTRING,FMT=*) HREAD_NAME
+ IF (TRIM(HREAD_NAME)==TRIM(HNAME)) THEN
+ KLINE = KLINE + 1
+ ELSE ! The name doesnt appear during a new read..
+ IF (KLINE > 0) THEN ! .. but it has already been found, ..
+ REWIND(KLUNAM) ! .. so it is the end of the data ..
+ RETURN ! .. and we can return :)
+ END IF
+ END IF
+ END IF
+END DO
+101 CONTINUE
+ IF (KLINE == 0) THEN
+ CALL EOL_AIRFOILNOTFOUND_ERROR(HFILE,HNAME)
+ STOP
+ END IF
+END SUBROUTINE HOW_MANY_LINES_OF
+!#########################################################
+!
+!#########################################################
+FUNCTION GET_AIRFOIL_ID(TPTURBINE,TPBLADE,TPAIRFOIL,PRADIUS)
+! Allows to link an airfoil from a TPBLADE, at a specific radius (PRADIUS)
+! to the airfoils characteristics (TPAIRFOIL).
+! The result is an integer (IAID) that should be used like that :
+! IAID = GET_GET_AIRFOIL_ID(TPTURBINE,TPBLADE,TPAIRFOIL,PRADIUS)
+! TPAIRFOIL(IAID)%MEMBER_OF_TPAIRFOIL
+!
+USE MODD_EOL_ALM, ONLY : TURBINE, BLADE, AIRFOIL
+!
+IMPLICIT NONE
+!
+TYPE(TURBINE), INTENT(IN) :: TPTURBINE ! stored turbine data
+TYPE(BLADE), INTENT(IN) :: TPBLADE ! stored blade data
+TYPE(AIRFOIL), DIMENSION(:), INTENT(IN) :: TPAIRFOIL ! stored arifoil data
+REAL, INTENT(IN) :: PRADIUS ! Radius position studied
+INTEGER :: GET_AIRFOIL_ID
+!
+INTEGER :: INB_BDATA ! Total number of blade data
+INTEGER :: JBDATA ! Index over blade's data
+INTEGER :: JA ! Index over diffetents airfoils
+REAL, DIMENSION(SIZE(TPBLADE%XRAD)) :: ZDELTARAD ! 2*ZDELTARAD = section lenght
+!
+! Checking data
+IF ((PRADIUS < TPTURBINE%XR_MIN) .OR. (PRADIUS > TPTURBINE%XR_MAX)) THEN
+ PRINT*, 'The studied radius R = ', PRADIUS, ' is out of blade range : [', &
+ TPTURBINE%XR_MIN, ';', TPTURBINE%XR_MAX, ']'
+ RETURN
+END IF
+!
+! Preliminaires
+INB_BDATA = SIZE(TPBLADE%XRAD)
+!
+! Computes half length of sections
+ZDELTARAD(1) = TPBLADE%XRAD(1) - TPTURBINE%XR_MIN
+DO JBDATA=2,INB_BDATA-1
+ ZDELTARAD(JBDATA) = TPBLADE%XRAD(JBDATA) &
+ - TPBLADE%XRAD(JBDATA-1) &
+ - ZDELTARAD(JBDATA-1)
+END DO
+ZDELTARAD(INB_BDATA) = TPTURBINE%XR_MAX - TPBLADE%XRAD(INB_BDATA)
+!
+! Looking for the section at r=PRADIUS
+DO JBDATA=1,INB_BDATA
+ IF ((PRADIUS >= TPBLADE%XRAD(JBDATA)-ZDELTARAD(JBDATA)) .AND. &
+ (PRADIUS < TPBLADE%XRAD(JBDATA)+ZDELTARAD(JBDATA))) THEN
+! Looking for the ID of the airfoil of this section
+ DO JA=1,SIZE(TPAIRFOIL)
+ IF (TRIM(TPBLADE%CAIRFOIL(JBDATA)) == TRIM(TPAIRFOIL(JA)%CNAME)) THEN
+ GET_AIRFOIL_ID = JA
+ EXIT
+ END IF
+ END DO
+!
+ EXIT
+ END IF
+END DO
+!
+END FUNCTION GET_AIRFOIL_ID
+!#########################################################
diff --git a/src/MNH/eol_smear.f90 b/src/MNH/eol_smear.f90
new file mode 100755
index 0000000000000000000000000000000000000000..00ca039d0e365ab854320d0ec29f72310c7f5336
--- /dev/null
+++ b/src/MNH/eol_smear.f90
@@ -0,0 +1,120 @@
+!MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
+!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
+!MNH_LIC for details. version 1.
+! #######################
+ MODULE MODI_EOL_SMEAR
+! #######################
+!
+INTERFACE
+!
+SUBROUTINE SMEAR_1LIN(PFX, PFX_SMR)
+ REAL, DIMENSION(:,:,:), INTENT(IN) :: PFX ! Force to smear
+ REAL, DIMENSION(:,:,:), INTENT(OUT) :: PFX_SMR ! Smeared force
+END SUBROUTINE SMEAR_1LIN
+!
+SUBROUTINE SMEAR_3LIN(PFX, PFY, PFZ, PFX_SMR, PFY_SMR, PFZ_SMR)
+ REAL, DIMENSION(:,:,:), INTENT(IN) :: PFX, PFY, PFZ ! Force to smear
+ REAL, DIMENSION(:,:,:), INTENT(OUT) :: PFX_SMR, PFY_SMR, PFZ_SMR ! Smeared force
+END SUBROUTINE SMEAR_3LIN
+!
+END INTERFACE
+!
+END MODULE MODI_EOL_SMEAR
+!-------------------------------------------------------------------
+!
+!!**** *MODI_EOL_SMEAR* -
+!!
+!! PURPOSE
+!! -------
+!! Smear the forces of wind turbine to avoid numerical instabilities
+!!
+!! AUTHOR
+!! ------
+!! PA. Joulin *CNRM & IFPEN*
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 04/2018
+!! Modification 21/10/20 (PA. Joulin) Updated for a main version
+!!
+!!---------------------------------------------------------------
+!
+!#########################################################
+SUBROUTINE SMEAR_1LIN(PFX, PFX_SMR)
+!!
+!! METHOD
+!! ------
+!! A linear smearing is done is this subroutine.
+!!
+!---------------------------------------------------------
+!
+!* 0. DECLARATIONS
+!
+ USE MODI_SHUMAN, ONLY: MXF, MXM
+!
+!* 0.1 declarations of arguments
+!
+ REAL, DIMENSION(:,:,:), INTENT(IN) :: PFX ! Force to smear
+ REAL, DIMENSION(:,:,:), INTENT(OUT) :: PFX_SMR ! Smeared force
+!
+!* 1. SMEARING
+!
+ PFX_SMR(:,:,:) = MXF(MXM(PFX(:,:,:)))
+!
+END SUBROUTINE SMEAR_1LIN
+!#########################################################
+!
+!#########################################################
+SUBROUTINE SMEAR_3LIN(PFX, PFY, PFZ, PFX_SMR, PFY_SMR, PFZ_SMR)
+!!
+!! METHOD
+!! ------
+!! A linear smearing is done is this subroutine.
+!!
+!---------------------------------------------------------
+!
+!* 0. DECLARATIONS
+!
+USE MODD_PARAMETERS, ONLY: JPVEXT
+USE MODI_SHUMAN, ONLY: MXF, MYF, MZF
+USE MODI_SHUMAN, ONLY: MXM, MYM, MZM
+!
+!* 0.1 declarations of arguments
+!
+ REAL, DIMENSION(:,:,:), INTENT(IN) :: PFX, PFY, PFZ ! Force to smear
+ REAL, DIMENSION(:,:,:), INTENT(OUT) :: PFX_SMR, PFY_SMR, PFZ_SMR ! Smeared force
+!
+!* 0.2 declarations of local variables
+!
+ INTEGER :: IKU,IKB,IKE ! Vertical indices
+!
+!* 1. INITIALIZATIONS
+!
+ IKU = SIZE(PFX,3) ! Top of the domain end index
+ IKB=1+JPVEXT ! Vertical begin index
+ IKE=IKU-JPVEXT ! Vertical end index
+!
+!* 2. SMEARING
+!
+ PFX_SMR(:,:,:) = MXF(MXM(MYF(MYM(PFX(:,:,:)))))
+ PFY_SMR(:,:,:) = MXF(MXM(MYF(MYM(PFY(:,:,:)))))
+ PFZ_SMR(:,:,:) = MXF(MXM(MYF(MYM(PFZ(:,:,:)))))
+!
+ PFX_SMR(:,:,:) = MZF(MZM(PFX_SMR(:,:,:)))
+ PFY_SMR(:,:,:) = MZF(MZM(PFY_SMR(:,:,:)))
+ PFZ_SMR(:,:,:) = MZF(MZM(PFZ_SMR(:,:,:)))
+!
+!* 3. BOUNDARY VALUES
+!
+ PFX_SMR(:,:,IKB-1) = PFX_SMR(:,:,IKB)
+ PFX_SMR(:,:,IKE+1) = PFX_SMR(:,:,IKE)
+!
+ PFY_SMR(:,:,IKB-1) = PFY_SMR(:,:,IKB)
+ PFY_SMR(:,:,IKE+1) = PFY_SMR(:,:,IKE)
+!
+ PFZ_SMR(:,:,IKB-1) = PFZ_SMR(:,:,IKB)
+ PFZ_SMR(:,:,IKE+1) = PFZ_SMR(:,:,IKE)
+!
+END SUBROUTINE SMEAR_3LIN
+!#########################################################
diff --git a/src/MNH/ini_budget.f90 b/src/MNH/ini_budget.f90
old mode 100644
new mode 100755
index 30d77166c7bfb08f66e72fcf5f4f9e94e0233d4b..fcbd15b4b9d855515a5f4f2e776580c8c589cb64
--- a/src/MNH/ini_budget.f90
+++ b/src/MNH/ini_budget.f90
@@ -106,7 +106,7 @@ end subroutine Budget_preallocate
OHORELAX_UVWTH,OHORELAX_RV,OHORELAX_RC,OHORELAX_RR, &
OHORELAX_RI,OHORELAX_RS, OHORELAX_RG, OHORELAX_RH,OHORELAX_TKE, &
OHORELAX_SV, OVE_RELAX, ove_relax_grd, OCHTRANS, &
- ONUDGING,ODRAGTREE,ODEPOTREE, &
+ ONUDGING,ODRAGTREE,ODEPOTREE, OAERO_EOL, &
HRAD,HDCONV,HSCONV,HTURB,HTURBDIM,HCLOUD )
! #################################################################
!
@@ -297,6 +297,7 @@ LOGICAL, INTENT(IN) :: OCHTRANS ! switch to activate convective
LOGICAL, INTENT(IN) :: ONUDGING ! switch to activate nudging
LOGICAL, INTENT(IN) :: ODRAGTREE ! switch to activate vegetation drag
LOGICAL, INTENT(IN) :: ODEPOTREE ! switch to activate droplet deposition on tree
+LOGICAL, INTENT(IN) :: OAERO_EOL ! switch to activate wind turbine wake
CHARACTER (LEN=*), INTENT(IN) :: HRAD ! type of the radiation scheme
CHARACTER (LEN=*), INTENT(IN) :: HDCONV ! type of the deep convection scheme
CHARACTER (LEN=*), INTENT(IN) :: HSCONV ! type of the shallow convection scheme
@@ -615,10 +616,15 @@ if ( lbu_ru ) then
call Budget_source_add( tbudgets(NBUDGET_U), tzsource )
tzsource%cmnhname = 'DRAG'
- tzsource%clongname = 'drag force'
+ tzsource%clongname = 'drag force due to trees'
tzsource%lavailable = odragtree
call Budget_source_add( tbudgets(NBUDGET_U), tzsource )
+ tzsource%cmnhname = 'DRAGEOL'
+ tzsource%clongname = 'drag force due to wind turbine'
+ tzsource%lavailable = OAERO_EOL
+ call Budget_source_add( tbudgets(NBUDGET_U), tzsource )
+
tzsource%cmnhname = 'DRAGB'
tzsource%clongname = 'drag force due to buildings'
tzsource%lavailable = ldragbldg
@@ -750,10 +756,15 @@ if ( lbu_rv ) then
call Budget_source_add( tbudgets(NBUDGET_V), tzsource )
tzsource%cmnhname = 'DRAG'
- tzsource%clongname = 'drag force'
+ tzsource%clongname = 'drag force due to trees'
tzsource%lavailable = odragtree
call Budget_source_add( tbudgets(NBUDGET_V), tzsource )
+ tzsource%cmnhname = 'DRAGEOL'
+ tzsource%clongname = 'drag force due to wind turbine'
+ tzsource%lavailable = OAERO_EOL
+ call Budget_source_add( tbudgets(NBUDGET_V), tzsource )
+
tzsource%cmnhname = 'DRAGB'
tzsource%clongname = 'drag force due to buildings'
tzsource%lavailable = ldragbldg
diff --git a/src/MNH/ini_eol_adnr.f90 b/src/MNH/ini_eol_adnr.f90
new file mode 100755
index 0000000000000000000000000000000000000000..9dad33be28eb442336d8517fcd0cb08e2b08d04a
--- /dev/null
+++ b/src/MNH/ini_eol_adnr.f90
@@ -0,0 +1,135 @@
+!MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
+!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
+!MNH_LIC for details. version 1.
+!-----------------------------------------------------------------
+! ##########################
+ MODULE MODI_INI_EOL_ADNR
+! ##########################
+!
+INTERFACE
+!
+SUBROUTINE INI_EOL_ADNR
+!
+END SUBROUTINE INI_EOL_ADNR
+!
+END INTERFACE
+!
+END MODULE MODI_INI_EOL_ADNR
+!
+! ############################################################
+ SUBROUTINE INI_EOL_ADNR
+! ############################################################
+!
+!!**** *INI_EOL_ADNR*
+!!
+!! PURPOSE
+!! -------
+!! Routine to initialized the ADNR (wind turbine model) variables
+!!
+!!** METHOD
+!! ------
+!!
+!! EXTERNAL
+!! --------
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!!
+!! **MODD_EOL_SHARED_IO:
+!! for namelist NAM_EOL_ADNR (INPUT) :
+!! CHARACTER(LEN=100) :: CFARM_CSVDATA ! File to read, with farm data
+!! CHARACTER(LEN=100) :: CTURBINE_CSVDATA ! File to read, turbine data
+!! for ouputs :
+!! REAL, DIMENSION(:), ALLOCATABLE :: XTHRUT ! Thrust [N]
+!! REAL, DIMENSION(:), ALLOCATABLE :: XTORQT ! Torque [Nm]
+!! REAL, DIMENSION(:), ALLOCATABLE :: XPOWT ! Power [W]
+!! REAL, DIMENSION(:), ALLOCATABLE :: XTHRU_SUM ! Sum of thrust (N)
+!! REAL, DIMENSION(:), ALLOCATABLE :: XTORQ_SUM ! Sum of torque (Nm)
+!! REAL, DIMENSION(:), ALLOCATABLE :: XPOW_SUM ! Sum of power (W)
+!!
+!! **MODD_EOL_ADNR (OUTPUT):
+!! TYPE(FARM) :: TFARM
+!! TYPE(TURBINE) :: TTURBINE
+!! REAL, DIMENSION(:), ALLOCATABLE :: XA_INDU ! Induction factor(NumEol,données)
+!! REAL, DIMENSION(:), ALLOCATABLE :: XCT_D ! Adapted thrust coef (for U_d) [-]
+!!
+!! REFERENCE
+!! ---------
+!!
+!!
+!! AUTHOR
+!! ------
+!! PA. Joulin * Meteo France & IFPEN *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 31/05/18
+!! Modification 14/10/20 (PA. Joulin) Updated for a main version
+!!
+!--------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+!
+USE MODD_EOL_ADNR
+USE MODD_EOL_SHARED_IO, ONLY: CFARM_CSVDATA, CTURBINE_CSVDATA
+USE MODD_EOL_SHARED_IO, ONLY: XTHRUT, XTHRU_SUM
+USE MODI_EOL_READER, ONLY: READ_CSVDATA_FARM_ADNR
+USE MODI_EOL_READER, ONLY: READ_CSVDATA_TURBINE_ADNR
+USE MODI_EOL_PRINTER, ONLY: PRINT_DATA_FARM_ADNR
+USE MODI_EOL_PRINTER, ONLY: PRINT_DATA_TURBINE_ADNR
+! To print in output listing
+USE MODD_LUNIT_n, ONLY: TLUOUT
+!
+IMPLICIT NONE
+!
+! Integers
+INTEGER :: ILUOUT ! Output listing file
+!
+!
+!-------------------------------------------------------------------
+!
+!* 1. READING AND ALLOCATING DATA
+! ---------------------------
+! Reading in csv files
+! Allocation of TFARM and TTURBINE inside the function
+!
+!* 1.1 Wind farm data
+!
+CALL READ_CSVDATA_FARM_ADNR(40,TRIM(CFARM_CSVDATA),TFARM)
+!
+!* 1.2 Wind turbine data
+!
+CALL READ_CSVDATA_TURBINE_ADNR(41,TRIM(CTURBINE_CSVDATA),TTURBINE)
+!
+!
+!-------------------------------------------------------------------
+!
+!* 2. PRINTING DATA
+! -------------
+!
+!* 2.0 Output listing index
+ILUOUT= TLUOUT%NLU
+!
+!* 2.1 Wind farm data
+!
+CALL PRINT_DATA_FARM_ADNR(ILUOUT,TFARM)
+!
+!* 2.2 Wind turbine data
+!
+CALL PRINT_DATA_TURBINE_ADNR(ILUOUT,TTURBINE)
+!
+!
+!-------------------------------------------------------------------
+!
+!* 3. ALLOCATING VARIABLES
+! --------------------
+!
+ALLOCATE(XA_INDU (TFARM%NNB_TURBINES))
+ALLOCATE(XCT_D (TFARM%NNB_TURBINES))
+ALLOCATE(XTHRUT (TFARM%NNB_TURBINES))
+ALLOCATE(XTHRU_SUM(TFARM%NNB_TURBINES))
+!
+END SUBROUTINE INI_EOL_ADNR
diff --git a/src/MNH/ini_eol_alm.f90 b/src/MNH/ini_eol_alm.f90
new file mode 100755
index 0000000000000000000000000000000000000000..73b4001fde46ed9866b31543af58906e6f7091c2
--- /dev/null
+++ b/src/MNH/ini_eol_alm.f90
@@ -0,0 +1,439 @@
+!MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
+!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
+!MNH_LIC for details. version 1.
+!-----------------------------------------------------------------
+! ##########################
+ MODULE MODI_INI_EOL_ALM
+! ##########################
+!
+INTERFACE
+!
+SUBROUTINE INI_EOL_ALM(PDXX,PDYY)
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX,PDYY ! mesh size
+!
+END SUBROUTINE INI_EOL_ALM
+!
+END INTERFACE
+!
+END MODULE MODI_INI_EOL_ALM
+!
+! ############################################################
+ SUBROUTINE INI_EOL_ALM(PDXX,PDYY)
+! ############################################################
+!
+!!**** *INI_EOL_ALM* - routine to initialize the Actuator Line Model
+!! to simulate wind turbines
+!!
+!! PURPOSE
+!! -------
+!! Routine to initialized the ALM (wind turbine model) variables
+!!
+!!** METHOD
+!! ------
+!!
+!! EXTERNAL
+!! --------
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!!
+!! *MODD_EOL_SHARED_IO:
+!! *Namelist NAM_EOL_ADNR (INPUT) :
+!! CHARACTER(LEN=100) :: CFARM_CSVDATA ! File to read, with farm data
+!! CHARACTER(LEN=100) :: CTURBINE_CSVDATA ! File to read, turbine data
+!! CHARACTER(LEN=100) :: CBLADE_CSVDATA ! Blade file to read
+!! CHARACTER(LEN=100) :: CAIRFOIL_CSVDATA ! Airfoil file to read
+!! *for ouputs :
+!! REAL, DIMENSION(:), ALLOCATABLE :: XTHRUT ! Thrust [N]
+!! REAL, DIMENSION(:), ALLOCATABLE :: XTORQT ! Torque [Nm]
+!! REAL, DIMENSION(:), ALLOCATABLE :: XPOWT ! Power [W]
+!! REAL, DIMENSION(:), ALLOCATABLE :: XTHRU_SUM ! Sum of thrust (N)
+!! REAL, DIMENSION(:), ALLOCATABLE :: XTORQ_SUM ! Sum of torque (Nm)
+!! REAL, DIMENSION(:), ALLOCATABLE :: XPOW_SUM ! Sum of power (W)
+!!
+!! INTEGER :: NNB_BLAELT ! Number of blade elements
+!!
+!! *MODD_EOL_ALM (OUTPUT):
+!! TYPE(FARM) :: TFARM
+!! TYPE(TURBINE) :: TTURBINE
+!! TYPE(BLADE) :: TBLADE
+!! TYPE(AIRFOIL), DIMENSION(:), ALLOCATABLE :: TAIRFOIL
+!! REAL, DIMENSION(:,:,:), ALLOCATABLE :: XELT_RAD ! Blade elements radius [m]
+!! REAL, DIMENSION(:,:,:), ALLOCATABLE :: XAOA_GLB ! Angle of attack of an element [rad]
+!! REAL, DIMENSION(:,:,:), ALLOCATABLE :: XFLIFT_GLB ! Lift force, parallel to Urel [N]
+!! REAL, DIMENSION(:,:,:), ALLOCATABLE :: XFDRAG_GLB ! Drag force, perpendicular to Urel [N]
+!! REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: XFAERO_RE_GLB ! Aerodyn. force (lift+drag) in RE [N]
+!! REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: XFAERO_RG_GLB ! Aerodyn. force (lift+drag) in RG [N]
+!! REAL, DIMENSION(:,:,:), ALLOCATABLE :: XAOA_SUM ! Sum of angle of attack [rad]
+!! REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: XFAERO_RE_SUM ! Sum of aerodyn. force (lift+drag) in RE [N]
+!
+!! *MODD_EOL_KINEMATICS
+!! Positions
+!! Orientations
+!! Velocities
+!!
+!! REFERENCE
+!! ---------
+!!
+!!
+!! AUTHOR
+!! ------
+!! PA. Joulin * Meteo France & IFPEN *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 31/05/18
+!! Modification 10/11/20 (PA. Joulin) Updated for a main version
+!!
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+!
+!* 0.1 Modules
+!
+USE MODD_EOL_ALM
+USE MODD_EOL_SHARED_IO, ONLY: CFARM_CSVDATA
+USE MODD_EOL_SHARED_IO, ONLY: CTURBINE_CSVDATA
+USE MODD_EOL_SHARED_IO, ONLY: CBLADE_CSVDATA
+USE MODD_EOL_SHARED_IO, ONLY: CAIRFOIL_CSVDATA
+USE MODD_EOL_SHARED_IO, ONLY: XTHRUT, XTORQT, XPOWT
+USE MODD_EOL_SHARED_IO, ONLY: XTHRU_SUM, XTORQ_SUM, XPOW_SUM
+USE MODI_EOL_READER, ONLY: READ_CSVDATA_FARM_ALM
+USE MODI_EOL_READER, ONLY: READ_CSVDATA_TURBINE_ALM
+USE MODI_EOL_READER, ONLY: READ_CSVDATA_BLADE_ALM
+USE MODI_EOL_READER, ONLY: READ_CSVDATA_AIRFOIL_ALM
+USE MODI_EOL_PRINTER, ONLY: PRINT_DATA_FARM_ALM
+USE MODI_EOL_PRINTER, ONLY: PRINT_DATA_TURBINE_ALM
+USE MODI_EOL_PRINTER, ONLY: PRINT_DATA_BLADE_ALM
+USE MODI_EOL_PRINTER, ONLY: PRINT_DATA_AIRFOIL_ALM
+USE MODD_EOL_KINE_ALM
+USE MODI_EOL_MATHS
+! To print in output listing
+USE MODD_LUNIT_n, ONLY: TLUOUT
+! Constant
+USE MODD_CST, ONLY: XPI
+! To know the grid
+USE MODD_GRID_n, ONLY: XXHAT,XYHAT,XZS
+USE MODE_ll, ONLY: GET_INDICE_ll
+USE MODD_PARAMETERS, ONLY: JPVEXT
+! MPI stuffs
+USE MODD_VAR_ll, ONLY: NMNH_COMM_WORLD
+USE MODD_PRECISION, ONLY: MNHREAL_MPI
+USE MODD_MPIF, ONLY: MPI_SUM
+!
+!* 0.2 Variables
+!
+IMPLICIT NONE
+! Interface
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX,PDYY ! mesh size
+!
+! Some loop controlers
+! .. for wind turbines
+INTEGER :: JROT ! Rotor index
+INTEGER :: JTELT ! Tower element index
+INTEGER :: JNELT ! Nacelle element index
+INTEGER :: JBLA ! Blade index
+INTEGER :: JBELT ! Balde element index
+! .. for domain
+INTEGER :: IIB,IJB,IKB ! Begin of a CPU domain
+INTEGER :: IIE,IJE ! End of a CPU domain
+INTEGER :: JI,JJ ! Domain index
+! Some variables to be coder-friendly
+INTEGER :: INB_WT, INB_B, INB_BELT ! Total numbers of wind turbines, blades, and blade elt
+INTEGER :: INB_TELT, INB_NELT ! Total numbers of tower elt, and nacelle elt
+REAL :: ZRAD ! Radius along the blade
+! Tower base folowing the terrain
+REAL,DIMENSION(:,:),ALLOCATABLE :: ZPOSINI_TOWO_RG ! Initial tower origin position
+INTEGER :: IINFO ! code info return
+!
+!-------------------------------------------------------------------
+!
+!* 1. READING AND ALLOCATING DATA
+! ---------------------------
+! Reading in csv files
+! Allocation of TFARM, TTURBINE, TBLADE and TAIRFOILS inside the function
+!
+!* 1.1 Wind farm data
+!
+CALL READ_CSVDATA_FARM_ALM(40,TRIM(CFARM_CSVDATA),TFARM)
+!
+!* 1.2 Wind turbine data
+!
+CALL READ_CSVDATA_TURBINE_ALM(41,TRIM(CTURBINE_CSVDATA),TTURBINE)
+!
+!* 1.3 Blade data
+!
+CALL READ_CSVDATA_BLADE_ALM(42,TRIM(CBLADE_CSVDATA),TTURBINE,TBLADE)
+!
+!* 1.4 Airfoil data
+!
+CALL READ_CSVDATA_AIRFOIL_ALM(43,TRIM(CAIRFOIL_CSVDATA),TBLADE,TAIRFOIL)
+!
+!
+!-------------------------------------------------------------------
+!
+!* 2. PRINTING DATA
+! -------------
+!
+!* 2.1 Wind farm data
+!
+CALL PRINT_DATA_FARM_ALM(TLUOUT%NLU,TFARM)
+!
+!* 2.2 Wind turbine data
+!
+CALL PRINT_DATA_TURBINE_ALM(TLUOUT%NLU,TTURBINE)
+!
+!* 2.3 Blade data
+!
+CALL PRINT_DATA_BLADE_ALM(TLUOUT%NLU,TBLADE)
+!
+!* 2.4 Airfoil data
+!
+CALL PRINT_DATA_AIRFOIL_ALM(TLUOUT%NLU,TAIRFOIL)
+!
+!
+!-------------------------------------------------------------------
+!
+!* 3. ALLOCATING VARIABLES
+! --------------------
+!
+!* 3.0 Preliminaries
+INB_WT = TFARM%NNB_TURBINES
+INB_B = TTURBINE%NNB_BLADES
+INB_BELT = TBLADE%NNB_BLAELT
+! Hard coded variables, but they will be usefull in next updates
+INB_TELT = 2
+INB_NELT = 2
+!
+!* 3.1 MODD_EOL_ALM variables
+! at t
+ALLOCATE(XELT_RAD (INB_WT,INB_B,INB_BELT ))
+ALLOCATE(XAOA_GLB (INB_WT,INB_B,INB_BELT ))
+ALLOCATE(XFDRAG_GLB (INB_WT,INB_B,INB_BELT ))
+ALLOCATE(XFLIFT_GLB (INB_WT,INB_B,INB_BELT ))
+ALLOCATE(XFAERO_RE_GLB(INB_WT,INB_B,INB_BELT,3))
+ALLOCATE(XFAERO_RG_GLB(INB_WT,INB_B,INB_BELT,3))
+ALLOCATE(XTHRUT (INB_WT ))
+ALLOCATE(XTORQT (INB_WT ))
+ALLOCATE(XPOWT (INB_WT ))
+! for mean values
+ALLOCATE(XAOA_SUM (INB_WT,INB_B,INB_BELT ))
+ALLOCATE(XFAERO_RE_SUM(INB_WT,INB_B,INB_BELT,3))
+ALLOCATE(XTHRU_SUM (INB_WT))
+ALLOCATE(XTORQ_SUM (INB_WT))
+ALLOCATE(XPOW_SUM (INB_WT))
+!
+!* 3.2 MODD_EOL_KINE_ALM variables
+!
+! - ORIENTATION MATRIX -
+ALLOCATE(XMAT_RG_RT(INB_WT,3,3))
+ALLOCATE(XMAT_RG_RN(INB_WT,3,3))
+ALLOCATE(XMAT_RT_RN(INB_WT,3,3))
+ALLOCATE(XMAT_RG_RH(INB_WT,3,3))
+ALLOCATE(XMAT_RH_RG(INB_WT,3,3))
+ALLOCATE(XMAT_RN_RH(INB_WT,3,3))
+ALLOCATE(XMAT_RG_RB(INB_WT,INB_B,3,3))
+ALLOCATE(XMAT_RH_RB(INB_WT,INB_B,3,3))
+ALLOCATE(XMAT_RG_RE(INB_WT,INB_B,INB_BELT,3,3))
+ALLOCATE(XMAT_RE_RG(INB_WT,INB_B,INB_BELT,3,3))
+ALLOCATE(XMAT_RB_RE(INB_WT,INB_B,INB_BELT,3,3))
+!
+! - POSITIONS & ORIENTATIONS -
+! Tower
+ALLOCATE(ZPOSINI_TOWO_RG (INB_WT,3) ) ! Initial tower origin pos. in RG
+ALLOCATE(XPOSINI_TOWO_RG (INB_WT,3) ) ! Initial tower origin pos. in RG
+ALLOCATE(XPOS_TOWO_RG (INB_WT,3) ) ! Current tower origin pos. in RG
+ALLOCATE(XPOS_TELT_RG (INB_WT,INB_TELT,3) ) ! Current tower element pos. in RG
+ALLOCATE(XPOS_TELT_RT (INB_WT,INB_TELT,3) ) ! Current tower element pos. in RT
+ALLOCATE(XANGINI_TOW_RG (INB_WT,3) ) ! Initial tower ori. in RG
+! Nacelle
+ALLOCATE(XPOSINI_NACO_RT (INB_WT,3) ) ! Initial nacelle origin pos. in RT
+ALLOCATE(XPOS_NACO_RG (INB_WT,3) ) ! Current nacelle origin pos. in RG
+ALLOCATE(XPOS_NELT_RG (INB_WT,INB_NELT,3) ) ! Current nacelle element pos. in RG
+ALLOCATE(XPOS_NELT_RN (INB_WT,INB_NELT,3) ) ! Current nacelle element pos. in RN
+ALLOCATE(XANGINI_NAC_RT (INB_WT,3) ) ! Initial nacelle ori. in RT
+! Hub
+ALLOCATE(XPOSINI_HUB_RN (INB_WT,3) ) ! Initial hub pos. in RN
+ALLOCATE(XPOS_HUB_RG (INB_WT,3) ) ! Current hub pos. in RG
+ALLOCATE(XANGINI_HUB_RN (INB_WT,3) ) ! Initial hub ori. in RN
+! Blade
+ALLOCATE(XPOSINI_BLA_RH (INB_WT,INB_B,3) ) ! Initial blade root pos. in RH
+ALLOCATE(XPOS_BLA_RG (INB_WT,INB_B,3) ) ! Current blade root pos. in RG
+ALLOCATE(XANGINI_BLA_RH (INB_WT,INB_B,3) ) ! Initial blade ori. RH
+! Element
+ALLOCATE(XPOS_ELT_RB (INB_WT,INB_B,INB_BELT,3) ) ! Element pos. in RB
+ALLOCATE(XPOS_ELT_RG (INB_WT,INB_B,INB_BELT,3) ) ! Element pos. in RG
+ALLOCATE(XPOS_SEC_RB (INB_WT,INB_B,INB_BELT+1,3)) ! Section pos. in RB
+ALLOCATE(XPOS_SEC_RG (INB_WT,INB_B,INB_BELT+1,3)) ! Section pos. in RG
+ALLOCATE(XANGINI_ELT_RB (INB_WT,INB_B,INB_BELT,3) ) ! Initial element ori. in RB
+ALLOCATE(XTWIST_ELT (INB_WT,INB_B,INB_BELT) ) ! Element twist in RB
+ALLOCATE(XCHORD_ELT (INB_WT,INB_B,INB_BELT) ) ! Element chord lenght
+ALLOCATE(XSURF_ELT (INB_WT,INB_B,INB_BELT) ) ! Element lift surface
+!
+! - STRUCTURAL VELOCITIES -
+! Tower
+ALLOCATE(XTVEL_TOWO_RG (INB_WT,3) ) ! Tower base trans. vel. in RG
+ALLOCATE(XTVEL_TELT_RG (INB_WT,INB_TELT,3) ) ! Tower element trans. vel. in RG
+ALLOCATE(XRVEL_RT_RG (INB_WT,3) ) ! RT/RG rot. vel.
+! Nacelle
+ALLOCATE(XTVEL_NACO_RT (INB_WT,3) ) ! Nacelle base trans. vel. in RT
+ALLOCATE(XTVEL_NELT_RG (INB_WT,INB_NELT,3) ) ! Nacelle element trans. vel. in RG
+ALLOCATE(XRVEL_RN_RT (INB_WT,3) ) ! RN/RT rot. vel.
+ALLOCATE(XRVEL_RN_RG (INB_WT,3) ) ! RN/RG rot. vel.
+! Hub
+ALLOCATE(XTVEL_HUB_RN (INB_WT,3) ) ! Hub base trans. vel. in RN
+ALLOCATE(XTVEL_HUB_RG (INB_WT,3) ) ! Hub base trans. vel. in RG
+ALLOCATE(XRVEL_RH_RN (INB_WT,3) ) ! RH/RT rot. vel.
+ALLOCATE(XRVEL_RH_RG (INB_WT,3) ) ! RH/RG rot. vel.
+! Blade
+ALLOCATE(XTVEL_BLA_RH (INB_WT,INB_B,3) ) ! Blade base trans. vel. in RH
+ALLOCATE(XTVEL_BLA_RG (INB_WT,INB_B,3) ) ! Blade base trans. vel. in RG
+ALLOCATE(XRVEL_RB_RH (INB_WT,INB_B,3) ) ! RB/RH rot. vel.
+ALLOCATE(XRVEL_RB_RG (INB_WT,INB_B,3) ) ! RB/RG rot. vel.
+! Elements
+ALLOCATE(XTVEL_ELT_RB (INB_WT,INB_B,INB_BELT,3) ) ! Element base trans. vel. in RB
+ALLOCATE(XTVEL_ELT_RG (INB_WT,INB_B,INB_BELT,3) ) ! Element base trans. vel. in RG
+ALLOCATE(XTVEL_ELT_RE (INB_WT,INB_B,INB_BELT,3) ) ! Element base trans. vel. in RE
+ALLOCATE(XRVEL_RE_RB (INB_WT,INB_B,INB_BELT,3) ) ! RE/RB rot. vel.
+ALLOCATE(XRVEL_RE_RG (INB_WT,INB_B,INB_BELT,3) ) ! RE/RG rot. vel.
+!
+!-------------------------------------------------------------------
+!
+!* 4. FIRST BUILDING OF WIND TURBINES
+! -------------------------------
+!
+!* 4.1 Preliminaries
+CALL GET_INDICE_ll(IIB,IJB,IIE,IJE) ! Get begin and end domain index (CPU)
+IKB=1+JPVEXT ! Vertical begin index
+!
+XPOS_REF(:) = 0d0 ! Global Origin
+!
+!
+!* 4.2 Tower
+DO JROT=1, INB_WT
+! Velocities (Rotor, (XYZ))
+ XTVEL_TOWO_RG(JROT,:) = 0d0
+ XRVEL_RT_RG(JROT,:) = 0d0
+! Positions
+ ! Init
+ ZPOSINI_TOWO_RG(JROT,1) = 0
+ ZPOSINI_TOWO_RG(JROT,2) = 0
+ ZPOSINI_TOWO_RG(JROT,3) = 0
+ ! Finding the elevation at the WT position
+ DO JJ=IJB,IJE
+ DO JI=IIB,IIE
+ IF (TFARM%XPOS_X(JROT) >= XXHAT(JI) .AND. &
+ TFARM%XPOS_X(JROT) < XXHAT(JI) + PDXX(JI,JJ,IKB)) THEN
+ IF (TFARM%XPOS_Y(JROT) >= XYHAT(JJ) .AND. &
+ TFARM%XPOS_Y(JROT) < XYHAT(JJ) + PDYY(JI,JJ,IKB)) THEN
+ ! Horizontal position
+ ZPOSINI_TOWO_RG(JROT,1) = TFARM%XPOS_X(JROT) ! Tower base position
+ ZPOSINI_TOWO_RG(JROT,2) = TFARM%XPOS_Y(JROT)
+ ! Finding the elevation at the WT position
+ ZPOSINI_TOWO_RG(JROT,3) = XZS(JI,JJ)
+ END IF
+ END IF
+ END DO
+ END DO
+ ! Sharing information
+ CALL MPI_ALLREDUCE(ZPOSINI_TOWO_RG, XPOSINI_TOWO_RG, SIZE(XPOSINI_TOWO_RG),&
+ MNHREAL_MPI,MPI_SUM,NMNH_COMM_WORLD,IINFO)
+ ! Tower elements
+ DO JTELT=1, INB_TELT
+ XPOS_TELT_RT(JROT,JTELT,1) = 0d0
+ XPOS_TELT_RT(JROT,JTELT,2) = 0d0
+ XPOS_TELT_RT(JROT,JTELT,3) = (JTELT-1)*(TTURBINE%XH_HEIGHT)/(INB_TELT-1)
+ END DO
+ ! Angles
+ XANGINI_TOW_RG(JROT,1) = 0d0
+ XANGINI_TOW_RG(JROT,2) = 0d0
+ XANGINI_TOW_RG(JROT,3) = 0d0
+!
+!
+!* 4.3 Nacelle
+! Velocities
+ XTVEL_NACO_RT(JROT,:) = 0d0
+ XRVEL_RN_RT(JROT,:) = 0d0
+! Positions (Rotor, (XYZ)) ! From last point of tower
+ ! Origin
+ XPOSINI_NACO_RT(JROT,:) = 0d0 ! Distance between nacelle base and tower top
+ ! Elements
+ DO JNELT=1, INB_NELT
+ XPOS_NELT_RN(JROT,JNELT,1) = (JNELT-1)*TTURBINE%XH_DEPORT/(INB_NELT-1)
+ XPOS_NELT_RN(JROT,JNELT,2) = 0d0
+ XPOS_NELT_RN(JROT,JNELT,3) = 0d0
+ END DO
+ ! Angles
+ XANGINI_NAC_RT(JROT,1) = 0d0
+ XANGINI_NAC_RT(JROT,2) = 0d0
+ XANGINI_NAC_RT(JROT,3) = XPI + TFARM%XNAC_YAW(JROT)
+!
+!
+!* 4.4 Hub
+! Velocities
+ XTVEL_HUB_RN(JROT,:) = 0d0
+ XRVEL_RH_RN(JROT,1) = 0d0
+ XRVEL_RH_RN(JROT,2) = 0d0
+ XRVEL_RH_RN(JROT,3) = TFARM%XOMEGA(JROT)
+! Position (Rotor, (XYZ)) ! From nacelle last point
+ XPOSINI_HUB_RN(JROT,1) = 0d0
+ XPOSINI_HUB_RN(JROT,2) = 0d0
+ XPOSINI_HUB_RN(JROT,3) = 0d0
+ XANGINI_HUB_RN(JROT,1) = 0d0
+ XANGINI_HUB_RN(JROT,2) = XPI/2d0 - TTURBINE%XNAC_TILT
+ XANGINI_HUB_RN(JROT,3) = 0d0
+!
+!
+!* 4.5 Blades
+ DO JBLA=1, INB_B
+! Velocities
+ XRVEL_RB_RH(JROT,JBLA,:) = 0d0
+ XTVEL_BLA_RH(JROT,JBLA,:) = 0d0
+! Position (Rotor, Blade, (XYZ)) ! From hub point
+ XPOSINI_BLA_RH(JROT,JBLA,:) = 0d0
+ XANGINI_BLA_RH(JROT,JBLA,1) = 0d0
+ XANGINI_BLA_RH(JROT,JBLA,2) = 0d0
+ XANGINI_BLA_RH(JROT,JBLA,3) = (JBLA-1)*2d0*XPI/INB_B
+!
+!
+!* 4.5 Elements
+! - Positioning of sections (cuts)
+ DO JBELT=1, INB_BELT+1
+
+ XPOS_SEC_RB(JROT,JBLA,JBELT,1) = TTURBINE%XR_MIN + (JBELT-1) &
+ * (TTURBINE%XR_MAX - TTURBINE%XR_MIN)/INB_BELT
+ XPOS_SEC_RB(JROT,JBLA,JBELT,2) = 0d0
+ XPOS_SEC_RB(JROT,JBLA,JBELT,3) = 0d0
+ ENDDO
+ DO JBELT=1, INB_BELT
+! - Positioning of centers (points of application)
+ XPOS_ELT_RB(JROT,JBLA,JBELT,1) = XPOS_SEC_RB(JROT,JBLA,JBELT,1) &
+ + (XPOS_SEC_RB(JROT,JBLA,JBELT+1,1) &
+ - XPOS_SEC_RB(JROT,JBLA,JBELT,1))/2d0
+ XPOS_ELT_RB(JROT,JBLA,JBELT,2) = 0d0
+ XPOS_ELT_RB(JROT,JBLA,JBELT,3) = 0d0
+ XELT_RAD(JROT,JBLA,JBELT) = XPOS_ELT_RB(JROT,JBLA,JBELT,1)
+! - Calculating chord and twist
+ ZRAD = XELT_RAD(JROT,JBLA,JBELT)
+ XCHORD_ELT(JROT,JBLA,JBELT) = INTERP_SPLCUB(ZRAD, TBLADE%XRAD, TBLADE%XCHORD)
+ XTWIST_ELT(JROT,JBLA,JBELT) = INTERP_SPLCUB(ZRAD, TBLADE%XRAD, TBLADE%XTWIST)
+! - Calculating lifting surface
+ XSURF_ELT(JROT,JBLA,JBELT) = XCHORD_ELT(JROT,JBLA,JBELT) &
+ * (XPOS_SEC_RB(JROT,JBLA,JBELT+1,1) &
+ - XPOS_SEC_RB(JROT,JBLA,JBELT,1))
+! Velocities
+ XRVEL_RE_RB(JROT,JBLA,JBELT,:) = 0d0
+ XTVEL_ELT_RB(JROT,JBLA,JBELT,:) = 0d0
+! Orientation
+ XANGINI_ELT_RB(JROT,JBLA,JBELT,1) = -XPI/2d0 + TFARM%XBLA_PITCH(JROT) &
+ + XTWIST_ELT(JROT,JBLA,JBELT)
+ XANGINI_ELT_RB(JROT,JBLA,JBELT,2) = XPI/2d0
+ XANGINI_ELT_RB(JROT,JBLA,JBELT,3) = XPI/2d0
+ END DO ! Loop element
+ END DO ! Loop blade
+END DO ! Loop turbine
+!
+END SUBROUTINE INI_EOL_ALM
diff --git a/src/MNH/ini_mean_field.f90 b/src/MNH/ini_mean_field.f90
old mode 100644
new mode 100755
index 36eafb4586599ef980f804c73ac1674cea7448c8..e9c5161a9ecec74a65183e26507495ea2f682a3f
--- a/src/MNH/ini_mean_field.f90
+++ b/src/MNH/ini_mean_field.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
!MNH_LIC for details. version 1.
@@ -57,11 +57,17 @@ END MODULE MODI_INI_MEAN_FIELD
USE MODD_MEAN_FIELD_n
USE MODD_MEAN_FIELD
USE MODD_PARAM_n
-
+USE MODD_EOL_MAIN, ONLY: LMAIN_EOL, CMETH_EOL, NMODEL_EOL
+USE MODD_EOL_SHARED_IO, ONLY: XTHRU_SUM, XTORQ_SUM, XPOW_SUM
+USE MODD_EOL_ALM
+USE MODE_MODELN_HANDLER
+!
IMPLICIT NONE
!
+INTEGER :: IMI !Current model index
+!
MEAN_COUNT = 0
-
+!
XUM_MEAN = 0.0
XVM_MEAN = 0.0
XWM_MEAN = 0.0
@@ -69,14 +75,29 @@ XTHM_MEAN = 0.0
XTEMPM_MEAN = 0.0
IF (CTURB /= 'NONE') XTKEM_MEAN = 0.0
XPABSM_MEAN = 0.0
-
+!
XU2_MEAN = 0.0
XV2_MEAN = 0.0
XW2_MEAN = 0.0
XTH2_MEAN = 0.0
XTEMP2_MEAN = 0.0
XPABS2_MEAN = 0.0
-
+!
+IMI = GET_CURRENT_MODEL_INDEX()
+IF (LMAIN_EOL .AND. IMI==NMODEL_EOL) THEN
+ SELECT CASE(CMETH_EOL)
+ CASE('ADNR') ! Actuator Disc Non-Rotating
+ XTHRU_SUM = 0.0
+ CASE('ALM') ! Actuator Line Method
+ XAOA_SUM = 0.0
+ XFAERO_RE_SUM = 0.0
+ XTHRU_SUM = 0.0
+ XTORQ_SUM = 0.0
+ XPOW_SUM = 0.0
+ END SELECT
+END IF
+!
+!
XUM_MAX = -1.E20
XVM_MAX = -1.E20
XWM_MAX = -1.E20
diff --git a/src/MNH/ini_modeln.f90 b/src/MNH/ini_modeln.f90
old mode 100644
new mode 100755
index dea7f4bcc0509d1db7b63b85daec8b7c3e4d1747..49a15b336a7c5c6b955031b060089f4a9c8af0f7
--- a/src/MNH/ini_modeln.f90
+++ b/src/MNH/ini_modeln.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 1994-2020 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
!MNH_LIC for details. version 1.
@@ -339,6 +339,7 @@ USE MODD_DYN_n
USE MODD_DYNZD
USE MODD_DYNZD_n
USE MODD_ELEC_n, only: XCION_POS_FW, XCION_NEG_FW
+USE MODD_EOL_MAIN
USE MODD_FIELD_n
#ifdef MNH_FOREFIRE
USE MODD_FOREFIRE
@@ -427,6 +428,8 @@ USE MODI_INI_DEEP_CONVECTION
USE MODI_INI_DRAG
USE MODI_INI_DYNAMICS
USE MODI_INI_ELEC_n
+USE MODI_INI_EOL_ADNR
+USE MODI_INI_EOL_ALM
USE MODI_INI_LES_N
USE MODI_INI_LG
USE MODI_INI_LW_SETUP
@@ -1671,7 +1674,7 @@ IF ( CBUTYPE /= "NONE" .AND. NBUMOD == KMI ) THEN
LHORELAX_UVWTH,LHORELAX_RV, LHORELAX_RC,LHORELAX_RR, &
LHORELAX_RI,LHORELAX_RS,LHORELAX_RG, LHORELAX_RH,LHORELAX_TKE, &
LHORELAX_SV, LVE_RELAX, LVE_RELAX_GRD, &
- LCHTRANS,LNUDGING,LDRAGTREE,LDEPOTREE, &
+ LCHTRANS,LNUDGING,LDRAGTREE,LDEPOTREE,LMAIN_EOL, &
CRAD,CDCONV,CSCONV,CTURB,CTURBDIM,CCLOUD )
END IF
!
@@ -2508,6 +2511,24 @@ IF (LCHEMDIAG) THEN
ELSE
ALLOCATE(XTCHEM(0))
END IF
-
+!-------------------------------------------------------------------------------
+!
+!* 32. Wind turbine
+!
+IF (LMAIN_EOL .AND. KMI == NMODEL_EOL) THEN
+ ALLOCATE(XFX_RG(IIU,IJU,IKU))
+ ALLOCATE(XFY_RG(IIU,IJU,IKU))
+ ALLOCATE(XFZ_RG(IIU,IJU,IKU))
+ ALLOCATE(XFX_SMR_RG(IIU,IJU,IKU))
+ ALLOCATE(XFY_SMR_RG(IIU,IJU,IKU))
+ ALLOCATE(XFZ_SMR_RG(IIU,IJU,IKU))
+ SELECT CASE(CMETH_EOL)
+ CASE('ADNR')
+ CALL INI_EOL_ADNR
+ CASE('ALM')
+ CALL INI_EOL_ALM(XDXX,XDYY)
+ END SELECT
+END IF
+!
END SUBROUTINE INI_MODEL_n
diff --git a/src/MNH/mean_field.f90 b/src/MNH/mean_field.f90
old mode 100644
new mode 100755
index b2fd6e083d992ea0e245bd333d2ba7402327896e..9ff392a18aa9940e7fcd1268cf923007366ae40a
--- a/src/MNH/mean_field.f90
+++ b/src/MNH/mean_field.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
!MNH_LIC for details. version 1.
@@ -10,12 +10,11 @@
!
INTERFACE
- SUBROUTINE MEAN_FIELD(PUT, PVT, PWT, PTHT, PTKET,PPABST)
+ SUBROUTINE MEAN_FIELD(PUT, PVT, PWT, PTHT, PTKET,PPABST)
REAL, DIMENSION(:,:,:), INTENT(IN) :: PUT, PVT, PWT ! variables
REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHT, PTKET ! variables
REAL, DIMENSION(:,:,:), INTENT(IN) :: PPABST ! variables
-
END SUBROUTINE MEAN_FIELD
END INTERFACE
@@ -46,6 +45,7 @@ END MODULE MODI_MEAN_FIELD
!! -------------
!! Original 07/2009
!! (C.Lac) 09/2016 Max values
+!! (PA.Joulin) 12/2020 Wind turbine variables
!!---------------------------------------------------------------
!
!
@@ -57,8 +57,14 @@ USE MODD_MEAN_FIELD_n
USE MODD_PARAM_n
USE MODD_MEAN_FIELD
USE MODD_CST
-
-!
+!
+USE MODD_EOL_MAIN, ONLY: LMAIN_EOL, CMETH_EOL, NMODEL_EOL
+USE MODD_EOL_SHARED_IO, ONLY: XTHRUT, XTORQT, XPOWT
+USE MODD_EOL_SHARED_IO, ONLY: XTHRU_SUM, XTORQ_SUM, XPOW_SUM
+USE MODD_EOL_ALM
+USE MODD_EOL_ADNR
+USE MODE_MODELN_HANDLER
+!
IMPLICIT NONE
!* 0.1 Declarations of dummy arguments :
@@ -66,12 +72,15 @@ IMPLICIT NONE
REAL, DIMENSION(:,:,:), INTENT(IN) :: PUT, PVT, PWT ! variables
REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHT, PTKET ! variables
REAL, DIMENSION(:,:,:), INTENT(IN) :: PPABST ! variables
-
!
!* 0.2 Declarations of local variables :
REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: ZTEMPT
INTEGER :: IIU,IJU,IKU,IIB,IJB,IKB,IIE,IJE,IKE ! Arrays bounds
INTEGER :: JI,JJ,JK ! Loop indexes
+!
+INTEGER :: IMI !Current model index
+!
+!
!-----------------------------------------------------------------------
!
!* 0. ARRAYS BOUNDS INITIALIZATION
@@ -103,6 +112,21 @@ IKE=IKU-JPVEXT
XTH2_MEAN = PTHT**2 + XTH2_MEAN
XTEMP2_MEAN = ZTEMPT**2 + XTEMP2_MEAN
XPABS2_MEAN = PPABST**2 + XPABS2_MEAN
+!
+! Wind turbine variables
+ IMI = GET_CURRENT_MODEL_INDEX()
+ IF (LMAIN_EOL .AND. IMI==NMODEL_EOL) THEN
+ SELECT CASE(CMETH_EOL)
+ CASE('ADNR') ! Actuator Disc Non-Rotating
+ XTHRU_SUM = XTHRUT + XTHRU_SUM
+ CASE('ALM') ! Actuator Line Method
+ XAOA_SUM = XAOA_GLB + XAOA_SUM
+ XFAERO_RE_SUM = XFAERO_RE_GLB + XFAERO_RE_SUM
+ XTHRU_SUM = XTHRUT + XTHRU_SUM
+ XTORQ_SUM = XTORQT + XTORQ_SUM
+ XPOW_SUM = XPOWT + XPOW_SUM
+ END SELECT
+ END IF
!
MEAN_COUNT = MEAN_COUNT + 1
!
diff --git a/src/MNH/modd_eol_adnr.f90 b/src/MNH/modd_eol_adnr.f90
new file mode 100755
index 0000000000000000000000000000000000000000..b2342570938d9508292e7a7003fa01bd90094cf4
--- /dev/null
+++ b/src/MNH/modd_eol_adnr.f90
@@ -0,0 +1,68 @@
+!MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
+!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
+!MNH_LIC for details. version 1.
+!!
+!! #####################
+ MODULE MODD_EOL_ADNR
+!! #####################
+!!
+!!*** *MODD_EOL_ADNR*
+!!
+!! PURPOSE
+!! -------
+!! It is possible to include wind turbines parameterization in Meso-NH,
+!! and several models are available. One of the models is the Non-Rotating
+!! Actuator Disk (ADNR). MODD_EOL_ADNR contains all the declarations for
+!! the ADNR model.
+!!
+!!** AUTHOR
+!! ------
+!! PA. Joulin *CNRM & IFPEN*
+!
+!! MODIFICATIONS
+!! -------------
+!! Original 04/16
+!! Modification 14/10/20 (PA. Joulin) Updated for a main version
+!!
+!-----------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! -----------------
+!
+IMPLICIT NONE
+!
+! ------ TYPES ------
+TYPE FARM
+ INTEGER :: NNB_TURBINES ! Number of wind turbines [-]
+ REAL, DIMENSION(:), ALLOCATABLE :: XPOS_X ! Tower base position, X coord [m]
+ REAL, DIMENSION(:), ALLOCATABLE :: XPOS_Y ! Tower base position, Y coord [m]
+ REAL, DIMENSION(:), ALLOCATABLE :: XCT_INF ! Thrust coefficient from U_infinite [-]
+END TYPE FARM
+!
+TYPE TURBINE
+ CHARACTER(LEN=10) :: CNAME ! Wind turbine name [-]
+ REAL :: XH_HEIGHT ! Hub height [m]
+ REAL :: XR_MAX ! Radius max of the blade [m]
+END TYPE TURBINE
+!
+! ------ VARIABLES ------
+! Farm data
+TYPE(FARM) :: TFARM
+TYPE(TURBINE) :: TTURBINE
+!
+! Global (CPU) variables
+REAL, DIMENSION(:), ALLOCATABLE :: XA_INDU ! Induction factor
+REAL, DIMENSION(:), ALLOCATABLE :: XCT_D ! Adapted thrust coef (for U_d) [-]
+!
+! Implicit from MODD_EOL_SHARED_IO:
+!REAL, DIMENSION(:), ALLOCATABLE :: XTHRUT ! Thrust [N]
+!REAL, DIMENSION(:), ALLOCATABLE :: XTHRU_SUM ! Sum of thrust (N)
+!
+! Namelist NAM_EOL_ADNR :
+! Implicit from MODD_EOL_SHARED_IO:
+!CHARACTER(LEN=100) :: CFARM_CSVDATA ! File to read, with farm data
+!CHARACTER(LEN=100) :: CTURBINE_CSVDATA ! File to read, turbine data
+!CHARACTER(LEN=3) :: CINTERP ! Interpolation method for wind speed
+!
+END MODULE MODD_EOL_ADNR
diff --git a/src/MNH/modd_eol_alm.f90 b/src/MNH/modd_eol_alm.f90
new file mode 100755
index 0000000000000000000000000000000000000000..187b4b88ab04dc43106261871fec6741d2f86753
--- /dev/null
+++ b/src/MNH/modd_eol_alm.f90
@@ -0,0 +1,117 @@
+!MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
+!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
+!MNH_LIC for details. version 1.
+!!
+!! #####################
+ MODULE MODD_EOL_ALM
+!! #####################
+!!
+!!*** *MODD_EOL_ALM*
+!!
+!! PURPOSE
+!! -------
+!! It is possible to include wind turbines parameterization in Meso-NH,
+!! and several models are available. One of the models is the Actuator
+!! Line Method (ALM). MODD_EOL_ALM contains all the declarations for
+!! the ALM model.
+!!
+!!
+!!** AUTHOR
+!! ------
+!! PA.Joulin *CNRM & IFPEN*
+!
+!! MODIFICATIONS
+!! -------------
+!! Original 04/01/17
+!! Modification 14/10/20 (PA. Joulin) Updated for a main version
+!!
+!-----------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! -----------------
+!
+IMPLICIT NONE
+!
+! ------ TYPES ------
+!
+TYPE FARM
+ INTEGER :: NNB_TURBINES ! Number of wind turbines [-]
+ REAL, DIMENSION(:), ALLOCATABLE :: XPOS_X ! Tower base position, X coord [m]
+ REAL, DIMENSION(:), ALLOCATABLE :: XPOS_Y ! Tower base position, Y coord [m]
+ REAL, DIMENSION(:), ALLOCATABLE :: XOMEGA ! Rotor rotation speed [rad/s]
+ REAL, DIMENSION(:), ALLOCATABLE :: XNAC_YAW ! Nacelle yaw angle [rad]
+ REAL, DIMENSION(:), ALLOCATABLE :: XBLA_PITCH ! Blade picth angle [rad]
+END TYPE FARM
+!
+TYPE TURBINE
+ CHARACTER(LEN=10) :: CNAME ! Nom de la turbine [-]
+ INTEGER :: NNB_BLADES ! Number of blades [-]
+ REAL :: XH_HEIGHT ! Hub height [m]
+ REAL :: XH_DEPORT ! Hub deport [m]
+ REAL :: XNAC_TILT ! Tilt of the nacelle [m]
+ REAL :: XR_MIN ! Minimum blade radius [m]
+ REAL :: XR_MAX ! Maximum blade radius [m]
+END TYPE TURBINE
+!
+TYPE BLADE
+ INTEGER :: NNB_BLAELT ! Number of blade element
+ INTEGER :: NNB_BLADAT ! Number of blade data
+ REAL, DIMENSION(:), ALLOCATABLE :: XRAD ! Data node radius [m]
+ REAL, DIMENSION(:), ALLOCATABLE :: XCHORD ! Element chord [m]
+ REAL, DIMENSION(:), ALLOCATABLE :: XTWIST ! Element twist [rad]
+ CHARACTER(LEN=20), DIMENSION(:), ALLOCATABLE :: CAIRFOIL ! Arifoil name [-]
+END TYPE BLADE
+!
+TYPE AIRFOIL
+ CHARACTER(LEN=15) :: CNAME ! Airfoil name [-]
+ REAL, DIMENSION(:), ALLOCATABLE :: XAA ! Attack Angle [rad]
+ REAL, DIMENSION(:), ALLOCATABLE :: XRE ! Reynolds Number [-]
+ REAL, DIMENSION(:), ALLOCATABLE :: XCL ! Lift coefficient [-]
+ REAL, DIMENSION(:), ALLOCATABLE :: XCD ! Drag coefficient [-]
+ REAL, DIMENSION(:), ALLOCATABLE :: XCM ! Moment coefficient [-]
+END TYPE AIRFOIL
+!
+! ------ VARIABLES ------
+! --- Farm data ---
+TYPE(FARM) :: TFARM
+TYPE(TURBINE) :: TTURBINE
+TYPE(BLADE) :: TBLADE
+TYPE(AIRFOIL), DIMENSION(:), ALLOCATABLE :: TAIRFOIL
+!
+! --- Global variables (Code & CPU) ---
+REAL, DIMENSION(:,:,:), ALLOCATABLE :: XELT_RAD ! Blade elements radius [m]
+REAL, DIMENSION(:,:,:), ALLOCATABLE :: XAOA_GLB ! Angle of attack of an element [rad]
+REAL, DIMENSION(:,:,:), ALLOCATABLE :: XFLIFT_GLB ! Lift force, parallel to Urel [N]
+REAL, DIMENSION(:,:,:), ALLOCATABLE :: XFDRAG_GLB ! Drag force, perpendicular to Urel [N]
+REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: XFAERO_RE_GLB ! Aerodyn. force (lift+drag) in RE [N]
+REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: XFAERO_RG_GLB ! Aerodyn. force (lift+drag) in RG [N]
+!
+! Mean values
+REAL, DIMENSION(:,:,:), ALLOCATABLE :: XAOA_SUM ! Sum of angle of attack [rad]
+REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: XFAERO_RE_SUM ! Sum of aerodyn. force (lift+drag) in RE [N]
+!
+! Implicit from MODD_EOL_SHARED_IO :
+!REAL, DIMENSION(:), ALLOCATABLE :: XTHRUT ! Thrust [N]
+!REAL, DIMENSION(:), ALLOCATABLE :: XTORQT ! Torque [Nm]
+!REAL, DIMENSION(:), ALLOCATABLE :: XPOWT ! Power [W]
+!REAL, DIMENSION(:), ALLOCATABLE :: XTHRU_SUM ! Sum of thrust (N)
+!REAL, DIMENSION(:), ALLOCATABLE :: XTORQ_SUM ! Sum of torque (Nm)
+!REAL, DIMENSION(:), ALLOCATABLE :: XPOW_SUM ! Sum of power (W)
+!
+!
+! --- Namelist NAM_EOL_ALM ---
+! Implicit from MODD_EOL_SHARED_IO :
+!CHARACTER(LEN=100) :: CFARM_CSVDATA ! Farm file to read
+!CHARACTER(LEN=100) :: CTURBINE_CSVDATA ! Turbine file to read
+!CHARACTER(LEN=100) :: CBLADE_CSVDATA ! Blade file to read
+!CHARACTER(LEN=100) :: CAIRFOIL_CSVDATA ! Airfoil file to read
+!CHARACTER(LEN=3) :: CINTERP ! Interpolation method for wind speed
+!
+INTEGER :: NNB_BLAELT ! Number of blade elements
+!
+LOGICAL :: LTIMESPLIT ! Flag to apply Time splitting method
+LOGICAL :: LTIPLOSSG ! Flag to apply Glauert's tip loss correction
+LOGICAL :: LTECOUTPTS ! Flag to get Tecplot file output of element points
+!
+END MODULE MODD_EOL_ALM
diff --git a/src/MNH/modd_eol_kine_alm.f90 b/src/MNH/modd_eol_kine_alm.f90
new file mode 100755
index 0000000000000000000000000000000000000000..a78e65673fcdecc58787112cd9f0ebc5b590638b
--- /dev/null
+++ b/src/MNH/modd_eol_kine_alm.f90
@@ -0,0 +1,107 @@
+!MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
+!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
+!MNH_LIC for details. version 1.
+!!
+!! #####################
+ MODULE MODD_EOL_KINE_ALM
+!! #####################
+!!
+!!*** *MODD_EOL_KINE_ALM*
+!!
+!! PURPOSE
+!! -------
+! Declaration to take into account wind turbine motion
+!
+!!
+!!** AUTHOR
+!! ------
+!! PA.Joulin *CNRM & IFPEN*
+!
+!! MODIFICATIONS
+!! -------------
+!! Original 04/18
+!!
+!-----------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! -----------------
+!
+!
+! - Matrix to move from one fram to an other -
+DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE :: XMAT_RG_RT ! RG = Repere GLOBAL, RT = Repere TOWER
+DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE :: XMAT_RG_RN, XMAT_RT_RN ! RN = Repere NACELLE
+DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE :: XMAT_RG_RH, XMAT_RH_RG, XMAT_RN_RH ! RH = Repere HUB
+DOUBLE PRECISION, DIMENSION(:,:,:,:), ALLOCATABLE :: XMAT_RG_RB, XMAT_RH_RB ! RB = Repere BLADE
+DOUBLE PRECISION, DIMENSION(:,:,:,:,:), ALLOCATABLE :: XMAT_RG_RE, XMAT_RE_RG, XMAT_RB_RE ! RE = Repere ELEMENT
+
+! - POSITIONS & ORIENTATIONS -
+DOUBLE PRECISION, DIMENSION(3) :: XPOS_REF ! Reference position
+
+! Tower
+DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: XPOSINI_TOWO_RG ! Initial tower origin position, in global reference frame
+DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: XPOS_TOWO_RG ! Current tower origin real position, in global reference frame
+DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE :: XPOS_TELT_RG ! Current tower element position, in global reference frame
+DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE :: XPOS_TELT_RT ! Current tower element position, in tower frame
+DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: XANGINI_TOW_RG ! Initial tower orientation in global ref frame
+
+! Nacelle
+DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: XPOSINI_NACO_RT ! Initial nacelle position, in tower reference frame
+DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: XPOS_NACO_RG ! Initial nacelle position, in global reference frame
+DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE :: XPOS_NELT_RG ! Initial nacelle position, in global reference frame
+DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE :: XPOS_NELT_RN ! Initial nacelle position, in nacelle reference frame
+DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: XANGINI_NAC_RT ! Initial nacelle orientation, in tower reference frame
+
+! Hub
+DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: XPOSINI_HUB_RN ! Initial hub position, in nacelle reference frame
+DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: XPOS_HUB_RG ! Current hub position, in global reference frame
+DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: XANGINI_HUB_RN ! Initial hub orientation, in nacelle reference frame
+
+! Blade
+DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE :: XPOSINI_BLA_RH ! Initial blade root position, in hub reference frame
+DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE :: XPOS_BLA_RG ! Current blade root position, in global reference frame
+DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE :: XANGINI_BLA_RH ! Initial blade orientation, in hub reference frame
+
+! Element
+DOUBLE PRECISION, DIMENSION(:,:,:,:), ALLOCATABLE :: XPOS_ELT_RB ! Element position, in blade reference frame
+DOUBLE PRECISION, DIMENSION(:,:,:,:), ALLOCATABLE :: XPOS_ELT_RG ! Element position, in global reference frame
+DOUBLE PRECISION, DIMENSION(:,:,:,:), ALLOCATABLE :: XPOS_SEC_RB ! Section position, in blade reference frame
+DOUBLE PRECISION, DIMENSION(:,:,:,:), ALLOCATABLE :: XPOS_SEC_RG ! Section position, in global reference frame
+DOUBLE PRECISION, DIMENSION(:,:,:,:), ALLOCATABLE :: XANGINI_ELT_RB ! Initial element orientation in blade reference frame
+DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE :: XTWIST_ELT ! Element twist, interpolated from data
+DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE :: XCHORD_ELT ! Element chord lenght, interpolated from data
+DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE :: XSURF_ELT ! Element lift surface
+!
+!
+! - STRUCTURAL VELOCITIES -
+! Tower
+DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: XTVEL_TOWO_RG ! Tower base translation velocity, in global reference frame
+DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE :: XTVEL_TELT_RG ! Tower element velocity, in global reference frame
+DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: XRVEL_RT_RG ! RT/RG rotational velocity
+
+! Nacelle
+DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: XTVEL_NACO_RT ! Nacelle base translation velocity, in tower reference frame
+DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE :: XTVEL_NELT_RG ! Nacelle element translation velocity, in global reference frame
+DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: XRVEL_RN_RT ! RN/RT rotational velocity
+DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: XRVEL_RN_RG ! RN/RG rotational velocity
+
+! Hub
+DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: XTVEL_HUB_RN ! Hub base translation velocity, in global reference frame
+DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: XTVEL_HUB_RG ! Hub base translation velocity, in global reference frame
+DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: XRVEL_RH_RN ! RH/RN rotational velocity
+DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: XRVEL_RH_RG ! RH/RG rotational velocity
+
+! Blade
+DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE :: XTVEL_BLA_RH ! Blade base translation velocity, in global reference frame
+DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE :: XTVEL_BLA_RG ! Blade base translation velocity, in global reference frame
+DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE :: XRVEL_RB_RH ! RB/RH rotational velocity
+DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE :: XRVEL_RB_RG ! RB/RG rotational velocity
+
+! Elements
+DOUBLE PRECISION, DIMENSION(:,:,:,:), ALLOCATABLE :: XTVEL_ELT_RB ! Element base translation velocity, in global reference frame
+DOUBLE PRECISION, DIMENSION(:,:,:,:), ALLOCATABLE :: XTVEL_ELT_RG ! Element base translation velocity, in global reference frame
+DOUBLE PRECISION, DIMENSION(:,:,:,:), ALLOCATABLE :: XTVEL_ELT_RE ! Element base translation velocity, in element reference frame
+DOUBLE PRECISION, DIMENSION(:,:,:,:), ALLOCATABLE :: XRVEL_RE_RB ! RE/RB rotational velocity
+DOUBLE PRECISION, DIMENSION(:,:,:,:), ALLOCATABLE :: XRVEL_RE_RG ! RE/RG rotational velocity
+
+END MODULE MODD_EOL_KINE_ALM
diff --git a/src/MNH/modd_eol_main.f90 b/src/MNH/modd_eol_main.f90
new file mode 100755
index 0000000000000000000000000000000000000000..69d7baa1e0b4c2a100abcaafabf10eb98e9c58b9
--- /dev/null
+++ b/src/MNH/modd_eol_main.f90
@@ -0,0 +1,54 @@
+!MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
+!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
+!MNH_LIC for details. version 1.
+!!
+!! #####################
+ MODULE MODD_EOL_MAIN
+!! #####################
+!!
+!!*** *MODD_EOL_MAIN*
+!!
+!! PURPOSE
+!! -------
+!! It is possible to include wind turbines parameterization in Meso-NH,
+!! and several methods are available. MODD_EOL_MAIN contains all the
+!! main declarations.
+!!
+!!** AUTHOR
+!! ------
+!! PA.Joulin *CNRM & IFPEN*
+!
+!! MODIFICATIONS
+!! -------------
+!! Original 24/01/17
+!! Modification 14/10/20 (PA. Joulin) Updated for a main version
+!!
+!-----------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! -----------------
+USE MODD_PARAMETERS
+!
+!
+! Necessary for each models
+IMPLICIT NONE
+!
+! ------ VARIABLES ------
+!
+! Aerodynamic forces in cartesian mesh
+REAL, DIMENSION(:,:,:), ALLOCATABLE :: XFX_RG ! Along X in RG frame [F]
+REAL, DIMENSION(:,:,:), ALLOCATABLE :: XFY_RG ! Along Y in RG frame [F]
+REAL, DIMENSION(:,:,:), ALLOCATABLE :: XFZ_RG ! Along Z in RG frame [F]
+! Smeared forces
+REAL, DIMENSION(:,:,:), ALLOCATABLE :: XFX_SMR_RG ! Along X in RG frame [F]
+REAL, DIMENSION(:,:,:), ALLOCATABLE :: XFY_SMR_RG ! Along Y in RG frame [F]
+REAL, DIMENSION(:,:,:), ALLOCATABLE :: XFZ_SMR_RG ! ALong Z in RG frame [F]
+!
+! Namelist NAM_EOL :
+LOGICAL :: LMAIN_EOL ! Flag to take into account wind turbine
+CHARACTER(LEN=4) :: CMETH_EOL ! Aerodynamic method
+CHARACTER(LEN=4) :: CSMEAR ! Type of smearing
+INTEGER :: NMODEL_EOL ! Son number, where the wind farm is
+!
+END MODULE MODD_EOL_MAIN
diff --git a/src/MNH/modd_eol_shared_io.f90 b/src/MNH/modd_eol_shared_io.f90
new file mode 100755
index 0000000000000000000000000000000000000000..2838aa92ce355b9cc1ba51ab40b8ccce9cd67ec0
--- /dev/null
+++ b/src/MNH/modd_eol_shared_io.f90
@@ -0,0 +1,57 @@
+!MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
+!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
+!MNH_LIC for details. version 1.
+!!
+!! #####################
+ MODULE MODD_EOL_SHARED_IO
+!! #####################
+!!
+!!*** *MODD_EOL_SHARED_IO*
+!!
+!! PURPOSE
+!! -------
+!! It is possible to include wind turbines parameterization in Meso-NH,
+!! and several models are available.
+!!
+!! MODD_EOL_SHARED_IO contains the declarations for the inputs/output
+!! shared by the differents models.
+!!
+!!
+!!** AUTHOR
+!! ------
+!! PA.Joulin *CNRM & IFPEN*
+!
+!! MODIFICATIONS
+!! -------------
+!! Original 17/11/20
+!!
+!-----------------------------------------------------------------------------
+!
+IMPLICIT NONE
+!
+!* 1. INPUTS VAR
+! ---------------
+!
+! --- Namelistis NAM_EOL_ADNR/NAM_EOL_ALM ---
+! * .csv files
+CHARACTER(LEN=100) :: CFARM_CSVDATA ! Farm file to read
+CHARACTER(LEN=100) :: CTURBINE_CSVDATA ! Turbine file to read
+CHARACTER(LEN=100) :: CBLADE_CSVDATA ! Blade file to read
+CHARACTER(LEN=100) :: CAIRFOIL_CSVDATA ! Airfoil file to read
+! * flags
+CHARACTER(LEN=3) :: CINTERP ! Interpolation method for wind speed
+!
+!
+!* 2. OUTPUTS VAR
+! -----------------
+!
+! --- Thruts torque and power ---
+REAL, DIMENSION(:), ALLOCATABLE :: XTHRUT ! Thrust [N]
+REAL, DIMENSION(:), ALLOCATABLE :: XTORQT ! Torque [Nm]
+REAL, DIMENSION(:), ALLOCATABLE :: XPOWT ! Power [W]
+REAL, DIMENSION(:), ALLOCATABLE :: XTHRU_SUM ! Sum of thrust (N)
+REAL, DIMENSION(:), ALLOCATABLE :: XTORQ_SUM ! Sum of torque (Nm)
+REAL, DIMENSION(:), ALLOCATABLE :: XPOW_SUM ! Sum of power (W)
+!
+END MODULE MODD_EOL_SHARED_IO
diff --git a/src/MNH/modd_sub_modeln.f90 b/src/MNH/modd_sub_modeln.f90
old mode 100644
new mode 100755
index 7bf9f5301aafec7fca42c78aaf5fe7b30c22d564..351bb95b074094b5dc3add37f8a18b587ca221de
--- a/src/MNH/modd_sub_modeln.f90
+++ b/src/MNH/modd_sub_modeln.f90
@@ -41,7 +41,7 @@ TYPE SUB_MODEL_t
integer :: nfile_output_current = 0 ! Number of the current output file
REAL(kind=MNHTIME), DIMENSION(2) :: XT_START
REAL(kind=MNHTIME), DIMENSION(2) :: XT_STORE, XT_BOUND, XT_GUESS
- REAL(kind=MNHTIME), DIMENSION(2) :: XT_ADV, XT_SOURCES, XT_DRAG
+ REAL(kind=MNHTIME), DIMENSION(2) :: XT_ADV, XT_SOURCES, XT_DRAG, XT_EOL
REAL(kind=MNHTIME), DIMENSION(2) :: XT_ADVUVW, XT_GRAV, XT_VISC
REAL(kind=MNHTIME), DIMENSION(2) :: XT_DIFF, XT_RELAX, XT_PARAM, XT_SPECTRA
REAL(kind=MNHTIME), DIMENSION(2) :: XT_HALO, XT_RAD_BOUND, XT_PRESS
@@ -68,7 +68,7 @@ integer, pointer :: nfile_backup_current => Null()
integer, pointer :: nfile_output_current => Null()
REAL(kind=MNHTIME), DIMENSION(:), POINTER :: XT_START=>NULL()
REAL(kind=MNHTIME), DIMENSION(:), POINTER :: XT_STORE=>NULL(), XT_BOUND=>NULL(), XT_GUESS=>NULL()
-REAL(kind=MNHTIME), DIMENSION(:), POINTER :: XT_ADV=>NULL(), XT_SOURCES=>NULL(), XT_DRAG=>NULL()
+REAL(kind=MNHTIME), DIMENSION(:), POINTER :: XT_ADV=>NULL(), XT_SOURCES=>NULL(), XT_DRAG=>NULL(), XT_EOL=>NULL()
REAL(kind=MNHTIME), DIMENSION(:), POINTER :: XT_ADVUVW=>NULL(), XT_GRAV=>NULL()
REAL(kind=MNHTIME), DIMENSION(:), POINTER :: XT_DIFF=>NULL(), XT_RELAX=>NULL(), XT_PARAM=>NULL(), XT_SPECTRA=>NULL()
REAL(kind=MNHTIME), DIMENSION(:), POINTER :: XT_HALO=>NULL(), XT_RAD_BOUND=>NULL(), XT_PRESS=>NULL()
@@ -127,6 +127,7 @@ XT_ADVUVW=>SUB_MODEL_MODEL(KTO)%XT_ADVUVW
XT_GRAV=>SUB_MODEL_MODEL(KTO)%XT_GRAV
XT_SOURCES=>SUB_MODEL_MODEL(KTO)%XT_SOURCES
XT_DRAG=>SUB_MODEL_MODEL(KTO)%XT_DRAG
+XT_EOL=>SUB_MODEL_MODEL(KTO)%XT_EOL
XT_DIFF=>SUB_MODEL_MODEL(KTO)%XT_DIFF
XT_RELAX=>SUB_MODEL_MODEL(KTO)%XT_RELAX
XT_PARAM=>SUB_MODEL_MODEL(KTO)%XT_PARAM
diff --git a/src/MNH/modeln.f90 b/src/MNH/modeln.f90
old mode 100644
new mode 100755
index 3ae4c24cab15fcd04948ff3a2aafb90575db8707..585b131a369aa4ed33c2a3a1f66b00946260d0e9
--- a/src/MNH/modeln.f90
+++ b/src/MNH/modeln.f90
@@ -303,6 +303,7 @@ USE MODD_DYN_n
USE MODD_DYNZD
USE MODD_DYNZD_n
USE MODD_ELEC_DESCR
+USE MODD_EOL_MAIN
USE MODD_FIELD_n
USE MODD_FRC
USE MODD_FRC_n
@@ -714,6 +715,7 @@ IF (KTCOUNT == 1) THEN
XT_TURB = 0.0_MNHTIME
XT_MAFL = 0.0_MNHTIME
XT_DRAG = 0.0_MNHTIME
+ XT_EOL = 0.0_MNHTIME
XT_TRACER = 0.0_MNHTIME
XT_SHADOWS = 0.0_MNHTIME
XT_ELEC = 0.0_MNHTIME
@@ -1335,10 +1337,10 @@ XT_RELAX = XT_RELAX + ZTIME2 - ZTIME1 &
!
ZTIME1 = ZTIME2
!
-CALL PHYS_PARAM_n( KTCOUNT, TZBAKFILE, &
- XT_RAD, XT_SHADOWS, XT_DCONV, XT_GROUND, &
- XT_MAFL, XT_DRAG, XT_TURB, XT_TRACER, &
- ZTIME, ZWETDEPAER, GMASKkids, GCLOUD_ONLY )
+CALL PHYS_PARAM_n( KTCOUNT, TZBAKFILE, &
+ XT_RAD, XT_SHADOWS, XT_DCONV, XT_GROUND, &
+ XT_MAFL, XT_DRAG, XT_EOL, XT_TURB, XT_TRACER, &
+ ZTIME, ZWETDEPAER, GMASKkids, GCLOUD_ONLY )
!
IF (CDCONV/='NONE') THEN
XPACCONV = XPACCONV + XPRCONV * XTSTEP
@@ -2152,6 +2154,7 @@ IF (OEXIT) THEN
CALL TIME_STAT_ll(XT_TURB,ZTOT, ' TURB = '//CTURB ,'-')
CALL TIME_STAT_ll(XT_MAFL,ZTOT, ' MAFL = '//CSCONV,'-')
CALL TIME_STAT_ll(XT_CHEM,ZTOT, ' CHIMIE' ,'-')
+ CALL TIME_STAT_ll(XT_EOL,ZTOT, ' WIND TURBINE' ,'-')
CALL TIMING_LEGEND()
CALL TIME_STAT_ll(XT_COUPL,ZTOT, ' SET_COUPLING','=')
CALL TIME_STAT_ll(XT_RAD_BOUND,ZTOT, ' RAD_BOUND','=')
diff --git a/src/MNH/modn_eol.f90 b/src/MNH/modn_eol.f90
new file mode 100755
index 0000000000000000000000000000000000000000..e9ded370931082c900a32f6f206af3130519129d
--- /dev/null
+++ b/src/MNH/modn_eol.f90
@@ -0,0 +1,39 @@
+!MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
+!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
+!MNH_LIC for details. version 1.
+!!
+!! #####################
+ MODULE MODN_EOL
+!! #####################
+!!
+!!*** *MODN_EOL*
+!!
+!! PURPOSE
+!! -------
+!! NAM_EOL activate the parameterization of wind turbines, and allows
+!! the selection of the aerodynamic method.
+!!
+!!** AUTHOR
+!! ------
+!! PA. Joulin *CNRM & IFPEN*
+!
+!! MODIFICATIONS
+!! -------------
+!! Original 24/01/17
+!! Modification 14/10/20 (PA. Joulin) Updated for a main version
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+USE MODD_EOL_MAIN
+!!
+!-----------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! -----------------
+IMPLICIT NONE
+SAVE
+NAMELIST /NAM_EOL/ &
+ LMAIN_EOL,CMETH_EOL,CSMEAR,NMODEL_EOL
+!
+END MODULE MODN_EOL
diff --git a/src/MNH/modn_eol_adnr.f90 b/src/MNH/modn_eol_adnr.f90
new file mode 100755
index 0000000000000000000000000000000000000000..b5767e37f73d97e98cde58d91dbe7cfe15d90b70
--- /dev/null
+++ b/src/MNH/modn_eol_adnr.f90
@@ -0,0 +1,43 @@
+!MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
+!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
+!MNH_LIC for details. version 1.
+!!
+!! #####################
+ MODULE MODN_EOL_ADNR
+!! #####################
+!!
+!!*** *MODN_EOL_ADNR*
+!!
+!! PURPOSE
+!! -------
+!! NAM_EOL activate the parameterization of wind turbines, and several
+!! models are available. One of the models is the Non-Rotating Actuator
+!! Disk Non Rotating (ADNR).
+!! The aim of NAM_EOL_ADNR is to specify ADNR parameters.
+!!
+!!** AUTHOR
+!! ------
+!! PA. Joulin *CNRM & IFPEN*
+!
+!! MODIFICATIONS
+!! -------------
+!! Original 04/16
+!! Modification 14/10/20 (PA. Joulin) Updated for a main version
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+USE MODD_EOL_ADNR
+USE MODD_EOL_SHARED_IO
+!!
+!-----------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! -----------------
+IMPLICIT NONE
+SAVE
+NAMELIST /NAM_EOL_ADNR/ &
+ CFARM_CSVDATA, CTURBINE_CSVDATA, &
+ CINTERP
+!
+END MODULE MODN_EOL_ADNR
diff --git a/src/MNH/modn_eol_alm.f90 b/src/MNH/modn_eol_alm.f90
new file mode 100755
index 0000000000000000000000000000000000000000..c7aa7e613f39b540d1d713d3998e21be74de39b4
--- /dev/null
+++ b/src/MNH/modn_eol_alm.f90
@@ -0,0 +1,45 @@
+!MNH_LIC Copyright 1994-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
+!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
+!MNH_LIC for details. version 1.
+!!
+!! #####################
+ MODULE MODN_EOL_ALM
+!! #####################
+!!
+!!*** *MODN_EOL_ALM*
+!!
+!! PURPOSE
+!! -------
+!! NAM_EOL activate the parameterization of wind turbines, and several
+!! models are available. One of the models is the Actuator Line
+!! Method (ALM).
+!! The aim of NAM_EOL_ALM is to specify ALM parameters.
+!!
+!!** AUTHOR
+!! ------
+!! PA. Joulin *CNRM & IFPEN*
+!
+!! MODIFICATIONS
+!! -------------
+!! Original 24/01/17
+!! Modification 14/10/20 (PA. Joulin) Updated for a main version
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+USE MODD_EOL_ALM
+USE MODD_EOL_SHARED_IO
+!!
+!-----------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! -----------------
+IMPLICIT NONE
+SAVE
+NAMELIST /NAM_EOL_ALM/ &
+ CFARM_CSVDATA, CTURBINE_CSVDATA, CBLADE_CSVDATA, CAIRFOIL_CSVDATA, &
+ NNB_BLAELT, &
+ CINTERP, LTIMESPLIT, LTIPLOSSG, &
+ LTECOUTPTS
+!
+END MODULE MODN_EOL_ALM
diff --git a/src/MNH/phys_paramn.f90 b/src/MNH/phys_paramn.f90
old mode 100644
new mode 100755
index 921afbd475c9005f1a2b21261aaf63b281a1418f..7bc1e8a89e726cae05875b8cce6aa74eebb2fee0
--- a/src/MNH/phys_paramn.f90
+++ b/src/MNH/phys_paramn.f90
@@ -11,8 +11,8 @@
INTERFACE
!
SUBROUTINE PHYS_PARAM_n( KTCOUNT, TPFILE, &
- PRAD, PSHADOWS, PKAFR, PGROUND, PMAFL, PDRAG, PTURB, PTRACER, &
- PTIME_BU, PWETDEPAER, OMASKkids, OCLOUD_ONLY )
+ PRAD, PSHADOWS, PKAFR, PGROUND, PMAFL, PDRAG,PEOL, PTURB, &
+ PTRACER, PTIME_BU, PWETDEPAER, OMASKkids, OCLOUD_ONLY )
!
USE MODD_IO, ONLY: TFILEDATA
use modd_precision, only: MNHTIME
@@ -20,7 +20,7 @@ use modd_precision, only: MNHTIME
INTEGER, INTENT(IN) :: KTCOUNT ! temporal iteration count
TYPE(TFILEDATA), INTENT(IN) :: TPFILE ! Synchronous output file
! advection schemes
-REAL(kind=MNHTIME), DIMENSION(2), INTENT(INOUT) :: PRAD,PSHADOWS,PKAFR,PGROUND,PTURB,PMAFL,PDRAG,PTRACER ! to store CPU
+REAL(kind=MNHTIME), DIMENSION(2), INTENT(INOUT) :: PRAD,PSHADOWS,PKAFR,PGROUND,PTURB,PMAFL,PDRAG,PTRACER,PEOL ! to store CPU
! time for computing time
REAL(kind=MNHTIME), DIMENSION(2), INTENT(INOUT) :: PTIME_BU ! time used in budget&LES budgets statistics
REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PWETDEPAER
@@ -36,8 +36,8 @@ END MODULE MODI_PHYS_PARAM_n
!
! ########################################################################################
SUBROUTINE PHYS_PARAM_n( KTCOUNT, TPFILE, &
- PRAD, PSHADOWS, PKAFR, PGROUND, PMAFL, PDRAG, PTURB, PTRACER, &
- PTIME_BU, PWETDEPAER, OMASKkids, OCLOUD_ONLY )
+ PRAD, PSHADOWS, PKAFR, PGROUND, PMAFL, PEOL, PDRAG, PTURB, &
+ PTRACER, PTIME_BU, PWETDEPAER, OMASKkids, OCLOUD_ONLY )
! ########################################################################################
!
!!**** *PHYS_PARAM_n * -monitor of the parameterizations used by model _n
@@ -263,6 +263,7 @@ USE MODD_DRAGTREE_n
USE MODD_DUST
USE MODD_DYN
USE MODD_DYN_n
+USE MODD_EOL_MAIN, ONLY: LMAIN_EOL, CMETH_EOL, NMODEL_EOL
USE MODD_FIELD_n
USE MODD_FRC
USE MODD_FRC_n
@@ -327,6 +328,7 @@ USE MODI_EDDY_FLUX_n ! Ajout PP
USE MODI_EDDY_FLUX_ONE_WAY_n ! Ajout PP
USE MODI_EDDYUV_FLUX_n ! Ajout PP
USE MODI_EDDYUV_FLUX_ONE_WAY_n ! Ajout PP
+USE MODI_EOL_MAIN
USE MODI_GROUND_PARAM_n
USE MODI_PASPOL
USE MODI_RADIATIONS
@@ -346,7 +348,7 @@ IMPLICIT NONE
INTEGER, INTENT(IN) :: KTCOUNT ! temporal iteration count
TYPE(TFILEDATA), INTENT(IN) :: TPFILE ! Synchronous output file
! advection schemes
-REAL(kind=MNHTIME), DIMENSION(2), INTENT(INOUT) :: PRAD,PSHADOWS,PKAFR,PGROUND,PTURB,PMAFL,PDRAG,PTRACER ! to store CPU
+REAL(kind=MNHTIME), DIMENSION(2), INTENT(INOUT) :: PRAD,PSHADOWS,PKAFR,PGROUND,PTURB,PMAFL,PDRAG,PTRACER,PEOL ! to store CPU
! time for computing time
REAL(kind=MNHTIME), DIMENSION(2), INTENT(INOUT) :: PTIME_BU ! time used in budget&LES budgets statistics
REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PWETDEPAER
@@ -1245,7 +1247,7 @@ CALL SECOND_MNH2(ZTIME2)
PTRACER = PTRACER + ZTIME2 - ZTIME1
!-----------------------------------------------------------------------------
!
-!* 5. Drag force
+!* 5a. Drag force
! ----------
!
ZTIME1 = ZTIME2
@@ -1265,6 +1267,29 @@ PDRAG = PDRAG + ZTIME2 - ZTIME1 &
!
PTIME_BU = PTIME_BU + XTIME_LES_BU_PROCESS + XTIME_BU_PROCESS
!
+!* 5b. Drag force from wind turbines
+! -----------------------
+!
+ZTIME1 = ZTIME2
+XTIME_BU_PROCESS = 0.
+XTIME_LES_BU_PROCESS = 0.
+!
+IF (LMAIN_EOL .AND. IMI == NMODEL_EOL) THEN
+ CALL EOL_MAIN(KTCOUNT,XTSTEP, &
+ XDXX,XDYY,XDZZ, &
+ XRHODJ, &
+ XUT,XVT,XWT, &
+ XRUS, XRVS, XRWS )
+END IF
+!
+CALL SECOND_MNH2(ZTIME2)
+!
+PEOL = PEOL + ZTIME2 - ZTIME1 &
+ - XTIME_LES_BU_PROCESS - XTIME_BU_PROCESS
+!
+PTIME_BU = PTIME_BU + XTIME_LES_BU_PROCESS + XTIME_BU_PROCESS
+!
+!*
!-----------------------------------------------------------------------------
!
!* 6. TURBULENCE SCHEME
diff --git a/src/MNH/read_exsegn.f90 b/src/MNH/read_exsegn.f90
old mode 100644
new mode 100755
index 170904e6019807d2d4e856d4ee487613e53e3a42..5d1bb14704c045a16db425d97a16d797de0e61ae
--- a/src/MNH/read_exsegn.f90
+++ b/src/MNH/read_exsegn.f90
@@ -356,6 +356,9 @@ USE MODN_DUST
USE MODN_DYN
USE MODN_DYN_n ! to avoid the duplication of this routine for each model.
USE MODN_ELEC
+USE MODN_EOL
+USE MODN_EOL_ADNR
+USE MODN_EOL_ALM
#ifdef MNH_FOREFIRE
USE MODN_FOREFIRE
#endif
@@ -529,6 +532,12 @@ CALL POSNAM(ILUSEG,'NAM_DRAGTREEN',GFOUND,ILUOUT)
IF (GFOUND) READ(UNIT=ILUSEG,NML=NAM_DRAGTREEn)
CALL POSNAM(ILUSEG,'NAM_DRAGBLDGN',GFOUND,ILUOUT)
IF (GFOUND) READ(UNIT=ILUSEG,NML=NAM_DRAGBLDGn)
+CALL POSNAM(ILUSEG,'NAM_EOL',GFOUND,ILUOUT)
+IF (GFOUND) READ(UNIT=ILUSEG,NML=NAM_EOL)
+CALL POSNAM(ILUSEG,'NAM_EOL_ADNR',GFOUND,ILUOUT)
+IF (GFOUND) READ(UNIT=ILUSEG,NML=NAM_EOL_ADNR)
+CALL POSNAM(ILUSEG,'NAM_EOL_ALM',GFOUND,ILUOUT)
+IF (GFOUND) READ(UNIT=ILUSEG,NML=NAM_EOL_ALM)
!
IF (KMI == 1) THEN
WRITE(UNIT=ILUOUT,FMT="(' namelists common to all the models ')")
diff --git a/src/MNH/write_lfin.f90 b/src/MNH/write_lfin.f90
old mode 100644
new mode 100755
index 11e3dc1298e5368a6836277ff1586b754886a24a..c6a48c561cfa99011aa6f6995cca0a6a73f2f6ce
--- a/src/MNH/write_lfin.f90
+++ b/src/MNH/write_lfin.f90
@@ -173,6 +173,7 @@ END MODULE MODI_WRITE_LFIFM_n
! S. Bielli 02/2019: Sea salt: significant sea wave height influences salt emission; 5 salt modes
! P. Wautelet 10/04/2019: replace ABORT and STOP calls by Print_msg
! P. Tulet 02/2020: correction for dust and sea salts
+! PA. Joulin 12/2020: add wind turbine outputs
! P. Wautelet 10/03/2021: use scalar variable names for dust and salt
! P. Wautelet 11/03/2021: bugfix: correct name for NSV_LIMA_IMM_NUCL
!-------------------------------------------------------------------------------
@@ -276,7 +277,12 @@ USE MODD_ADVFRC_n ! Modif PP ADV FRC
USE MODD_RELFRC_n
!
USE MODD_PARAM_C2R2
-!
+!
+USE MODD_EOL_MAIN
+USE MODD_EOL_SHARED_IO
+USE MODD_EOL_ADNR
+USE MODD_EOL_ALM
+!
IMPLICIT NONE
!
!* 0.1 Declarations of arguments
@@ -2079,6 +2085,200 @@ IF ( CPROGRAM=='REAL ' ) THEN
!
END IF
!
+!* 1.15 Wind turbine variables
+!
+! i) Main
+!
+IF (LMAIN_EOL .AND. IMI == NMODEL_EOL) THEN
+ TZFIELD%NGRID = 1
+ TZFIELD%NTYPE = TYPEREAL
+ TZFIELD%NDIMS = 3
+ TZFIELD%CDIR = 'XYZ'
+ TZFIELD%CUNITS = 'N'
+!
+ TZFIELD%CMNHNAME = 'FX_RG'
+ TZFIELD%CLONGNAME = 'FX_RG'
+ TZFIELD%CCOMMENT = 'X-component field of aerodynamic force (wind->rotor) in global frame (N)'
+ CALL IO_Field_write(TPFILE,TZFIELD,XFX_RG)
+!
+ TZFIELD%CMNHNAME = 'FY_RG'
+ TZFIELD%CLONGNAME = 'FY_RG'
+ TZFIELD%CCOMMENT = 'Y-component field of aerodynamic force (wind->rotor) in global frame (N)'
+ CALL IO_Field_write(TPFILE,TZFIELD,XFY_RG)
+!
+ TZFIELD%CMNHNAME = 'FZ_RG'
+ TZFIELD%CLONGNAME = 'FZ_RG'
+ TZFIELD%CCOMMENT = 'Z-component field of aerodynamic force (wind->rotor) in global frame (N)'
+ CALL IO_Field_write(TPFILE,TZFIELD,XFZ_RG)
+!
+ TZFIELD%CMNHNAME = 'FX_SMR_RG'
+ TZFIELD%CLONGNAME = 'FX_SMR_RG'
+ TZFIELD%CCOMMENT = 'X-component field of smeared aerodynamic force (wind->rotor) in global frame (N)'
+ TZFIELD%CCOMMENT = ''
+ CALL IO_Field_write(TPFILE,TZFIELD,XFX_SMR_RG)
+!
+ TZFIELD%CMNHNAME = 'FY_SMR_RG'
+ TZFIELD%CLONGNAME = 'FY_SMR_RG'
+ TZFIELD%CCOMMENT = 'Y-component field of smeared aerodynamic force (wind->rotor) in global frame (N)'
+ CALL IO_Field_write(TPFILE,TZFIELD,XFY_SMR_RG)
+!
+ TZFIELD%CMNHNAME = 'FZ_SMR_RG'
+ TZFIELD%CLONGNAME = 'FZ_SMR_RG'
+ TZFIELD%CCOMMENT = 'Z-component field of smeared aerodynamic force (wind->rotor) in global frame (N)'
+ CALL IO_Field_write(TPFILE,TZFIELD,XFZ_SMR_RG)
+!
+SELECT CASE(CMETH_EOL)
+!
+! ii) Actuator Disk without Rotation model
+!
+ CASE('ADNR') ! Actuator Disc Non-Rotating
+!
+ TZFIELD%NGRID = 1
+ TZFIELD%NTYPE = TYPEREAL
+ TZFIELD%NDIMS = 1
+ TZFIELD%CDIR = ''
+ TZFIELD%CUNITS = '-'
+!
+ TZFIELD%CMNHNAME = 'A_INDU'
+ TZFIELD%CLONGNAME = 'INDUCTION_FACTOR'
+ TZFIELD%CCOMMENT = 'Induction factor (-)'
+ CALL IO_Field_write(TPFILE,TZFIELD,XA_INDU)
+!
+ TZFIELD%CMNHNAME = 'CT_D'
+ TZFIELD%CLONGNAME = 'CTHRUST_D'
+ TZFIELD%CCOMMENT = 'Thrust coefficient at disk (-), &
+ used with wind speed at disk'
+ CALL IO_Field_write(TPFILE,TZFIELD,XCT_D)
+!
+ TZFIELD%CMNHNAME = 'THRUT'
+ TZFIELD%CLONGNAME = 'THRUSTT_EOL'
+ TZFIELD%CUNITS = 'N'
+ TZFIELD%CCOMMENT = 'RID instantaneous thrust of the wind turbines (N)'
+ CALL IO_Field_write(TPFILE,TZFIELD,XTHRUT)
+!
+ IF (MEAN_COUNT /= 0) THEN
+
+ TZFIELD%CMNHNAME = 'THRUMME'
+ TZFIELD%CLONGNAME = 'MEAN_THRUST_EOL'
+ TZFIELD%CUNITS = 'N'
+ TZFIELD%CCOMMENT = 'RID mean thrust of the wind turbines (N)'
+ CALL IO_Field_write(TPFILE,TZFIELD,XTHRU_SUM/MEAN_COUNT)
+!
+ END IF
+! iii) Actuator Line Model
+!
+ CASE('ALM') ! Actuator Line Method
+!
+ TZFIELD%NGRID = 1
+ TZFIELD%NTYPE = TYPEREAL
+ TZFIELD%CDIR = ''
+!
+ TZFIELD%NDIMS = 1
+!
+ TZFIELD%CMNHNAME = 'THRUT'
+ TZFIELD%CLONGNAME = 'THRUSTT_EOL'
+ TZFIELD%CUNITS = 'N'
+ TZFIELD%CCOMMENT = 'RID instantaneous thrust (N) of wind turbines'
+ CALL IO_Field_write(TPFILE,TZFIELD,XTHRUT)
+!
+ TZFIELD%CMNHNAME = 'TORQT'
+ TZFIELD%CLONGNAME = 'TORQUET_EOL'
+ TZFIELD%CUNITS = 'Nm'
+ TZFIELD%CCOMMENT = 'RID instantaneous torque (Nm) of wind turbines'
+ CALL IO_Field_write(TPFILE,TZFIELD,XTORQT)
+!
+ TZFIELD%CMNHNAME = 'POWT'
+ TZFIELD%CLONGNAME = 'POWERT_EOL'
+ TZFIELD%CUNITS = 'W'
+ TZFIELD%CCOMMENT = 'RID instantaneous power (W) of wind turbines'
+ CALL IO_Field_write(TPFILE,TZFIELD,XPOWT)
+!
+ TZFIELD%NDIMS = 3
+!
+ TZFIELD%CMNHNAME = 'ELT_RAD'
+ TZFIELD%CLONGNAME = 'ELT_RAD'
+ TZFIELD%CUNITS = 'm'
+ TZFIELD%CCOMMENT = 'RID_BID_EID radius (m) of wind turbine blade elements'
+ CALL IO_Field_write(TPFILE,TZFIELD,XELT_RAD)
+!
+ TZFIELD%CMNHNAME = 'AOA'
+ TZFIELD%CLONGNAME = 'ANGLE OF ATTACK'
+ TZFIELD%CUNITS = 'rad'
+ TZFIELD%CCOMMENT = 'RID_BID_EID instantaneous angle of attack (rad)'
+ CALL IO_Field_write(TPFILE,TZFIELD,XAOA_GLB)
+!
+ TZFIELD%CMNHNAME = 'FLIFT'
+ TZFIELD%CLONGNAME = 'LIFT FORCE'
+ TZFIELD%CUNITS = 'N'
+ TZFIELD%CCOMMENT = 'RID_BID_EID instantaneous lift (N) in relative frame'
+ CALL IO_Field_write(TPFILE,TZFIELD,XFLIFT_GLB)
+!
+ TZFIELD%CMNHNAME = 'FDRAG'
+ TZFIELD%CLONGNAME = 'DRAG FORCE'
+ TZFIELD%CUNITS = 'N'
+ TZFIELD%CCOMMENT = 'RID_BID_EID instantaneous drag (N) in relative frame'
+ CALL IO_Field_write(TPFILE,TZFIELD,XFDRAG_GLB)
+!
+ TZFIELD%NDIMS = 4
+!
+ TZFIELD%CMNHNAME = 'FAERO_RE'
+ TZFIELD%CLONGNAME = 'AERODYNAMIC FORCE RE'
+ TZFIELD%CUNITS = 'N'
+ TZFIELD%CCOMMENT = 'RID_BID_EID_XYZ instantaneous forces (N) in RE'
+ CALL IO_Field_write(TPFILE,TZFIELD,XFAERO_RE_GLB)
+!
+ TZFIELD%CMNHNAME = 'FAERO_RG'
+ TZFIELD%CLONGNAME = 'AERODYNAMIC FORCE RG'
+ TZFIELD%CUNITS = 'N'
+ TZFIELD%CCOMMENT = 'RID_BID_EID_XYZ instantaneous forces (N) in RG'
+ CALL IO_Field_write(TPFILE,TZFIELD,XFAERO_RG_GLB)
+!
+ IF (MEAN_COUNT /= 0) THEN
+!
+ TZFIELD%NGRID = 1
+ TZFIELD%NTYPE = TYPEREAL
+ TZFIELD%CDIR = ''
+!
+ TZFIELD%NDIMS = 1
+!
+ TZFIELD%CMNHNAME = 'THRUMME'
+ TZFIELD%CLONGNAME = 'MEAN_THRUST_EOL'
+ TZFIELD%CUNITS = 'N'
+ TZFIELD%CCOMMENT = 'RID mean thrust of the wind turbines (N)'
+ CALL IO_Field_write(TPFILE,TZFIELD,XTHRU_SUM/MEAN_COUNT)
+!
+ TZFIELD%CMNHNAME = 'TORQMME'
+ TZFIELD%CLONGNAME = 'MEAN_TORQUE_EOL'
+ TZFIELD%CUNITS = 'Nm'
+ TZFIELD%CCOMMENT = 'RID mean torque of the wind turbines (Nm)'
+ CALL IO_Field_write(TPFILE,TZFIELD,XTORQ_SUM/MEAN_COUNT)
+!
+ TZFIELD%CMNHNAME = 'POWMME'
+ TZFIELD%CLONGNAME = 'MEAN_POWER_EOL'
+ TZFIELD%CUNITS = 'W'
+ TZFIELD%CCOMMENT = 'RID mean power of the wind turbines (W)'
+ CALL IO_Field_write(TPFILE,TZFIELD,XPOW_SUM/MEAN_COUNT)
+!
+ TZFIELD%NDIMS = 3
+!
+ TZFIELD%CMNHNAME = 'AOAMME'
+ TZFIELD%CLONGNAME = 'MEAN_ANGLE_OF_ATTACK'
+ TZFIELD%CUNITS = 'rad'
+ TZFIELD%CCOMMENT = 'RID_BID_EID mean angle of attack (rad)'
+ CALL IO_Field_write(TPFILE,TZFIELD,XAOA_SUM/MEAN_COUNT)
+!
+ TZFIELD%NDIMS = 4
+!
+ TZFIELD%CMNHNAME = 'FAEROMME_RE'
+ TZFIELD%CLONGNAME = 'MEAN_AERODYNAMIC_FORCE_RE'
+ TZFIELD%CUNITS = 'N'
+ TZFIELD%CCOMMENT = 'RID_BID_EID_XYZ mean forces (N) in RE'
+ CALL IO_Field_write(TPFILE,TZFIELD,XFAERO_RE_SUM/MEAN_COUNT)
+!
+ END IF
+!
+ END SELECT
+END IF
!
DEALLOCATE(ZWORK2D,ZWORK3D)
!