From d5cdbc3e353d91bea857a7d055af0dfdcecad025 Mon Sep 17 00:00:00 2001
From: Gabriele Farina <gabriele.farina@unibs.it>
Date: Fri, 21 Feb 2025 10:45:05 +0000
Subject: [PATCH] New vegetation friction law in TELEMAC-2D based on Cui et al.
 (2023)

---
 NEWS.txt                                      |   2 +
 documentation/data/biblio.bib                 |  11 ++
 .../telemac2d/user/latex/appendix5.tex        |  14 +-
 examples/telemac2d/vegetation/f2d_veg10.slf   |   3 +
 examples/telemac2d/vegetation/friction10.tbl  | 163 ++++++++++++++++++
 examples/telemac2d/vegetation/friction9.tbl   |   6 +-
 .../vegetation/user_fortran/friction_zones.f  |  11 +-
 .../telemac2d/vegetation/vnv_vegetation.py    |  43 ++++-
 sources/api/api.cmdf                          |   1 +
 sources/telemac2d/friction_calc.f             |   2 +-
 sources/telemac2d/friction_cui.f              |  69 ++++++++
 sources/telemac2d/friction_def.f              |   4 +-
 sources/telemac2d/friction_read.f             |  16 +-
 sources/telemac2d/friction_zones.f            |  12 +-
 sources/telemac2d/interface_telemac2d.f       |  12 ++
 sources/telemac2d/telemac2d.cmdf              |   1 +
 16 files changed, 348 insertions(+), 22 deletions(-)
 create mode 100644 examples/telemac2d/vegetation/f2d_veg10.slf
 create mode 100644 examples/telemac2d/vegetation/friction10.tbl
 create mode 100644 sources/telemac2d/friction_cui.f

diff --git a/NEWS.txt b/NEWS.txt
index 8e936ce4b9..915d473204 100644
--- a/NEWS.txt
+++ b/NEWS.txt
@@ -1,6 +1,8 @@
 Latest changes
 ==============
 
+TELEMAC-2D: Add Cui vegetation law
+
 TELEMAC-2D: Improvements of the friction laws of telemac2d
 - Fix the dico to allow Colebrook-White law to be used
 - Improved the fixed point algorithm of the Colebrook-White for better convergence
diff --git a/documentation/data/biblio.bib b/documentation/data/biblio.bib
index 186d1dea90..97a018ca02 100644
--- a/documentation/data/biblio.bib
+++ b/documentation/data/biblio.bib
@@ -1,3 +1,14 @@
+@article{cui2023predicting,
+  title={Predicting flow resistance in open-channel flows with submerged vegetation},
+  author={Cui, Hanwen and Felder, Stefan and Kramer, Matthias},
+  journal={Environmental Fluid Mechanics},
+  volume={23},
+  number={4},
+  pages={757--778},
+  year={2023},
+  publisher={Springer}
+}
+
 @InProceedings{Quetin1977,
   author = 	 {GAUTHIER M. and QUETIN B.},
   title = 	 {Modèles mathématiques de calcul des écoulements induits par le vent},
diff --git a/documentation/telemac2d/user/latex/appendix5.tex b/documentation/telemac2d/user/latex/appendix5.tex
index 63c855b7a9..cfd20d2d8e 100644
--- a/documentation/telemac2d/user/latex/appendix5.tex
+++ b/documentation/telemac2d/user/latex/appendix5.tex
@@ -112,7 +112,8 @@ Huthoff et al & 5 & HUTH & Cd, mD, Hp, sp \\ \hline
 Van Velzen et al & 6 & VANV & Cd, mD, Hp \\ \hline
 Luhar \& Nepf & 7 & LUHE & Cd, Cv, a, Hp \\ \hline
 Vastila \& Jaervelae & 8 & VAST & Cdf, LAI, Ureff, Vogelf, CDs,  SAI, Vogels, Urefs, Hp \\ \hline
-Folke et al & 1 & HYBR & Cdx, LAI, Uref, Vogel, Hp \\ \hline
+Folke et al & 9 & HYBR & Cdx, LAI, Uref, Vogel, Hp \\ \hline
+Cui et al \cite{cui2023predicting} & 10 & VCUI & Cd, A, HdVeg, Zeta, Y, Alfa, Beta, CPI  \\ \hline
 \end{tabular}
 
 Explanation of the used parameters: 
@@ -133,18 +134,21 @@ Vogelf: & Foliage Vogel exponent \\
 Vogels: & Stem Vogel exponent \\
 Hp: & Plant height \\
 D: & Vegetation diameter \\
-sp: & Vegetation spacing  \\
-NOTE: Huthoff et al defines sp as spacing between two trees (bark to bark)\\
+sp: & Vegetation spacing (Huthoff et al defines sp as spacing between two trees (bark to bark)) \\
 Ap0: & Initial projected area \\
 EI: & Flexural rigidity \\
 mD: & m: vegetation density (1/sp**2) D: vegetation diameter 0.5*LAI/Hp (Finnigan 2000) \\
 Cv: & Friction coefficient on top of the vegetation layer, for emerged case Cv=0\\
 a: & Frontal area of vegetation per volume (=mD) \\
+HdVeg: & Height of deflected vegetation \\
+Zeta: & Adimensional parameter ui/uUD  \\
+Y: & Adimensional parameter y0/Hdeg  \\
+Alfa: & Adimensional parameter Le/HdVeg  \\
+Beta: & Adimensional parameter yi/HdVeg  \\
+CPI: & Wake function coefficient \\
 \end{tabular}
 
 
-
-
 \textbf{III -- Steering file}
 
 In order to use a friction computation by domains, the next keywords have to be
diff --git a/examples/telemac2d/vegetation/f2d_veg10.slf b/examples/telemac2d/vegetation/f2d_veg10.slf
new file mode 100644
index 0000000000..bb6b4d44b9
--- /dev/null
+++ b/examples/telemac2d/vegetation/f2d_veg10.slf
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:b3fd8cb81e004aca5f10fc09825496ab6981478f61783c63609ba820487bcac5
+size 719116
diff --git a/examples/telemac2d/vegetation/friction10.tbl b/examples/telemac2d/vegetation/friction10.tbl
new file mode 100644
index 0000000000..e26edc3748
--- /dev/null
+++ b/examples/telemac2d/vegetation/friction10.tbl
@@ -0,0 +1,163 @@
+* -----------------------------------------------------------------------------
+*  Description of friction of the zones
+*
+*  Bed roughness laws
+*    NOFR : no friction         (no value)
+*    HAAL : Haaland   law       (1 Value  : rB)
+*    CHEZ : Chezy     law       (1 Value  : rB)
+*    STRI : Strickler law       (1 Value  : rB)
+*    MANN : Manning   law       (1 Value  : rB)
+*    NIKU : Nikuradse law       (1 Value  : rB)
+*    LOGW : Log Wall  law       (1 Value  : rB)
+*    COWH : Colebrook-White law (1 Value  : rB. nDef)
+*
+*  Vegetation roughness laws
+*   NULL : no vegetation        (no value)
+*   LIND : Lindner              (1982)  approach  (2 Values: dp, sp)
+*   JAER : Jarvela              (2004)  approach  (5 Values: Cd, LAI, Uref, VOGEL, hp)
+*   WHIT : Whittaker et al      (2015)  approach  (6 Values: CDo, Apo, EI, VOGEL, Sp, hp)  
+*   BAPT : Baptist et al        (2007)  approach  (3 Values: Cd, mD, hp)
+*   HUTH : Huthoff et al        (2007)  approach  (4 values: Cd, mD, hp, Sp)
+*   VANV : Van Velzen et al     (2003)  approach  (3 Values: Cd, mD,  hp)
+*   LUNE : Luhar & Nepf         (2013)  approach  (4 values: Cd, Cv, a, hp)
+*   VAST : Vastila and Jarvela  (2014)  approach  (8 values: Cdf, CDs, LAI, SAI, VOGELf, VOGELs, Ureff, Urefs)
+*   HYBR : Folke et al.         (2021)  approach  (5 Values: Cd, LAI, Uref, VOGEL, hp)
+*   VCUI : Cui et al            (2023)  approach  (8 values: Cd, A, HdVeg, Zeta, Y, Alfa, Beta, CPI)
+*
+*  no             : Number of zone
+*
+*  Bed setup
+*  -------------
+*  typeB          : Roughness law for sole boundary
+*  rB             : Roughness coefficient for bottom boundary
+*  nDefB          : Manning coefficient for flat overflow areas (option)
+*
+*  Walls setup (for k-epsilon model)
+*  -----------------------------------------
+*  typeS          : roughness law for lateral boundary (option)
+*  rS             : roughness factor for lateral boundary (option)
+*  nDefS          : Manning coefficient for flat overflow areas (option)
+*
+*****************************************************************************************************
+*  Vegetation parameters(if needed)
+*  ------------------------
+*  typeV          : Roughness law for vegetation
+*
+* Case Lindner:
+*  dp             : average diameter (option)
+*  sp             : mean distance from ideal roughness elements (option)
+*                   from center to center
+*
+*
+* Case Jarvela:
+*  Cdx            : vegetation drag coefficient (Specie specific)
+*  LAI            : Leaf area index
+*  Uref           : (Specie-specific) lowest velocity used to determine the Vogel exponent
+*  Vogel          : (Specie-specific) Vogel exponent
+*  hp             : Vegetation height
+*
+*
+* Case Whittaker:
+*  CDo            : Specie-specific rigid drag coefficient
+*  Apo            : Initial-projected area (Rigid)
+*  EI             : Flexural Rigidity
+*  VOGEL          : Vogel exponent
+*  Sp             : Spacing between roughness elements (same as Lindner)
+*  hp             : Undeflected plant height
+*
+* Case Baptist:
+*  Cd             : vegetation bulk drag coefficient (Normally 1.0)
+*  mD             : m: vegetation density=1/spacing**2
+*                   D: Vegetation diameter when modeled as cylinder      
+*                   Note: can be estimated by mD=0.5*LAI/hp (Finnigan 2000)
+*  hp             : vegetation height
+*
+* Case Huthoff (2007):
+*  Cd             : vegetation bulk drag coefficient
+*  mD             : m: vegetation density=1/spacing**2
+*                   D: Vegetation diameter when modeled as cylinder      
+*                   Note: can be estimated by mD=0.5*LAI/hp (Finnigan 2000)
+*  hp             : Vegetation height
+*  Sp             : mean spacing between ideal roughness elements
+*                   from bark to bark
+*
+* Case Van Velzen (2003):
+*  Cd             : vegetation drag coefficient
+*  mD             : m: vegetation density=1/spacing**2
+*                   D: Vegetation diameter when modeled as cylinder      
+*                   Note: can be estimated by mD=0.5*LAI/hp (Finnigan 2000)
+*  hp             : Vegetation height
+*
+*
+* Case Luhar and Nepf (2013):
+*  Cd             : vegetation drag coefficient
+*  Cv             : friction coefficient on top of the vegetation layer
+*                   Note: Case of emerged, input 0
+*  a              : frontal area of vegetation per volume (=mD)
+*                   Note: can be estimated by 0.5 LAI/H (Finnigan 2000)
+*  hp             : Vegetation height             
+*
+*
+* Case Vastila and Jarvela (2014):
+*  CDf            : Foliage drag coefficient
+*  CDs            : Stem drag coefficient
+*  LAI            : Leaf area index
+*  SAI            : Stem Area index
+*  VOGELf         : Foliage Vogel exponent
+*  VOGELs         : Stem Vogel exponent
+*  Ureff          : Foliage reference velocity (lowest velocity used to determine the vogel exponent)
+*  Urefs          : Stem reference velocity
+*
+*
+* Case Hybrid (Folke et al, 2021: Combination of Jarvela and Baptist):
+*  Cdx            : vegetation drag coefficient (Specie specific)
+*  LAI            : Leaf area index
+*  Uref           : (Specie-specific) lowest velocity used to determine the Vogel exponent
+*  Vogel          : (Specie-specific) Vogel exponent
+*  hp             : Vegetation height
+*
+*
+* Case Cui et al (2023):
+* Cd		  : Bulk drag coefficient for vegetation
+* a               : Frontal area per unit volume
+* HdVeg           : Height of deflected vegetation
+* Zeta            : Adimensional parameter ui/uUD (u_i: time-average velocity at the inflection point, u_UD: quasi-constant velocity distribution)
+* Y               : Adimensional parameter y0/Hdeg (y_0: roughness height)
+* Alfa            : Adimensional parameter Le/HdVeg (Le: characteristic length scale describing the penetration depth of vortices)
+* Beta            : Adimensional parameter yi/HdVeg (y_i: bed normal elevation of the inflection point)
+* CPI             : Wake function coefficient
+*
+* --------------------------------------------------------------------------------------------------------------------
+* Lindner:
+*     no    typeB    rB    NDefB   typeS   rS    NDefS   typeV   dp    sp
+* Jarvela:
+*     no    typeB    rB    NDefB   typeS   rS    NDefS   typeV   CDx    LAI   Ux   VOGEL   hp
+* Whittaker
+*     no    typeB    rB    NDefB   typeS   rS    NDefS   typeV   CDo   Apo   EI    VOGEL   Sp   hp
+* Baptist:
+*     no    typeB    rB    NDefB   typeS   rS    NDefS   typeV   CD    mD   hp 
+* Huthoff:
+*     no    typeB    rB    NDefB   typeS   rS    NDefS   typeV   CD    mD   hp  Sp
+* Van Velzen:
+*     no    typeB    rB    NDefB   typeS   rS    NDefS   typeV   CD    mD   hp
+* Luhar and Nepf:
+*     no    typeB    rB    NDefB   typeS   rS    NDefS   typeV   CD    Cv   a   hp
+* Vaestilae:
+*     no    typeB    rB    NDefB   typeV   CDf   LAI  Ureff   VOGELf   CDs  SAI  Urefs  VOGELs hp
+* Hybrid:
+*     no    typeB    rB    NDefB   typeS   rS    NDefS   typeV   CDx    LAI   Ux   VOGEL   hp
+* Cui et al:
+*     no    typeB    rB    NDefB   typeV   CD   a   HdVeg   Zeta   Y   Alfa   Beta   CPI 
+* --------------------------------------------------------------------------------------------------------------------
+*
+* main channel
+100         NIKU   0.0005   NULL    
+*
+* bank
+200         NIKU   0.0005   NULL
+*
+* floodplain
+300         NIKU   0.0030   VCUI         0.5  2.8  0.03    1.21  0.65  0.27  0.76  0.04
+*
+*
+END
diff --git a/examples/telemac2d/vegetation/friction9.tbl b/examples/telemac2d/vegetation/friction9.tbl
index e1b8c92f4c..1a7830880b 100644
--- a/examples/telemac2d/vegetation/friction9.tbl
+++ b/examples/telemac2d/vegetation/friction9.tbl
@@ -135,13 +135,13 @@
 *     no    typeB    rB    NDefB   typeS   rS    NDefS   typeV   CDx    LAI   Ux   VOGEL   hp
 * --------------------------------------------------------------------------------------------------------------------
 *
-* Flussschlauch
+* main channel
 100         NIKU   0.0005   NULL    
 *
-* Ufer
+* bank
 200         NIKU   0.0005   NULL
 *
-* Vorland
+* floodplain
 300         NIKU   0.0030   HYBR    0.5  0.1 0.1  -0.9    1.0
 *
 *
diff --git a/examples/telemac2d/vegetation/user_fortran/friction_zones.f b/examples/telemac2d/vegetation/user_fortran/friction_zones.f
index c0d6219b3a..a9c14504db 100644
--- a/examples/telemac2d/vegetation/user_fortran/friction_zones.f
+++ b/examples/telemac2d/vegetation/user_fortran/friction_zones.f
@@ -214,16 +214,23 @@
      &       VCOEFF%ADR(4)%P%R(I), VCOEFF%ADR(5)%P%R(I),
      &       KARMAN,CP)
           CF%R(I) =CF%R(I)+CP
+        CASE (10)
+          CALL FRICTION_CUI
+     &       (T1%R(I), VCOEFF%ADR(1)%P%R(I),
+     &       VCOEFF%ADR(2)%P%R(I), VCOEFF%ADR(3)%P%R(I),
+     &       VCOEFF%ADR(4)%P%R(I), VCOEFF%ADR(5)%P%R(I),
+     &       VCOEFF%ADR(6)%P%R(I), VCOEFF%ADR(7)%P%R(I),
+     &       VCOEFF%ADR(8)%P%R(I),KARMAN,CP)
+          CF%R(I) =CP
         CASE DEFAULT
 !
           WRITE(LU,2) VEGLAW%I(I)
- 2             FORMAT(I5,' : UNKNOWN FRICTION LAW')
+ 2             FORMAT(I5,' : UNKNOWN VEGETATION FRICTION LAW ')
         CALL PLANTE(1)
         STOP
         END SELECT
 
 
-
         ENDIF !VEGETATION
 !
       ENDDO
diff --git a/examples/telemac2d/vegetation/vnv_vegetation.py b/examples/telemac2d/vegetation/vnv_vegetation.py
index 79afb8a51f..fd107bfbf8 100644
--- a/examples/telemac2d/vegetation/vnv_vegetation.py
+++ b/examples/telemac2d/vegetation/vnv_vegetation.py
@@ -44,11 +44,20 @@ class VnvStudy(AbstractVnvStudy):
                        cas=cas)
         del cas
 
+        # veg_law 
+        # 1: Lindner, 
+        # 2: Jaervelae, 
+        # 3: Whittaker, 
+        # 4: Baptist, 
+        # 5: Huthoff
+        # 6: van Velzen, 
+        # 7: Luhar Nepf, 
+        # 8: Vaestilae, 
+        # 9: hybrid
+        # 10: Cui
 
-# veg_law 1: Lindner, 2: Jaervelae, 3: Whittaker, 4: Baptist, 5: Huthoff
-# 6: van Velzen, 7: Luhar Nepf, 8: Vaestilae, 9: hybrid
         # vegetation T2D scalar mode
-        for veg_law in [2, 3, 4, 5, 6, 7, 8, 9]:
+        for veg_law in [2, 3, 4, 5, 6, 7, 8, 9, 10]:
             t2d_steering_file = 't2d_veg{}.cas'.format(veg_law)
             veg_file = 'friction{}.tbl'.format(veg_law)
 
@@ -214,6 +223,16 @@ class VnvStudy(AbstractVnvStudy):
                             'vnv_veg9_par:T2DRES',
                             eps=[2e-3])
 
+        # Comparison with the last time frame of the reference file.
+        self.check_epsilons('vnv_veg10_seq:T2DRES',
+                            'f2d_veg10.slf',
+                            eps=[3e-3])
+
+       # Comparison between sequential and parallel run.
+        self.check_epsilons('vnv_veg10_seq:T2DRES',
+                            'vnv_veg10_par:T2DRES',
+                            eps=[1e-1])
+
     def _post(self):
         """
         Post-treatment processes
@@ -229,9 +248,20 @@ class VnvStudy(AbstractVnvStudy):
         res, _ = self.get_study_res('vnv_veg1_seq:T2DRES')
 
         # Load all results as a list:
-        res_list, res_labels = self.get_study_res(module='T2D')
-#        print ("res_list", res_list)
-        print ("res_labels", res_labels)
+        res_list, res_labels = self.get_study_res(module='T2D', whitelist=['seq'])
+
+        res_labels = [
+          "Lindner", 
+          "Jaervelae", 
+          "Whittaker", 
+          "Baptist", 
+          "Huthoff",
+          "van Velzen", 
+          "Luhar Nepf", 
+          "Vaestilae", 
+          "hybrid",
+          "Cui"]
+
         # Plot bathy longitudinal section:
         vnv_plot1d_polylines(\
             'BOTTOM',
@@ -246,6 +276,7 @@ class VnvStudy(AbstractVnvStudy):
         vnv_plot1d_polylines(\
             'BOTTOM',
             res,
+            x_label='y (m)',
             poly=[[30, 0], [30, 5]],
             fig_size=(8, 2),
             record=0,
diff --git a/sources/api/api.cmdf b/sources/api/api.cmdf
index 9e0dd89df4..aaf7fbef2b 100644
--- a/sources/api/api.cmdf
+++ b/sources/api/api.cmdf
@@ -1075,6 +1075,7 @@ files: friction_def.f
   friction_calc.f
   friction_unif.f
   friction_baptist.f
+  friction_cui.f
   friction_huthoff.f
   friction_hybrid.f
   friction_jaervelae.f
diff --git a/sources/telemac2d/friction_calc.f b/sources/telemac2d/friction_calc.f
index 45b5d50227..642038ca51 100644
--- a/sources/telemac2d/friction_calc.f
+++ b/sources/telemac2d/friction_calc.f
@@ -220,7 +220,7 @@
       CASE DEFAULT
 !
         WRITE(LU,2) KFROT
-2       FORMAT(I5,' : UNKNOWN FRICTION LAW')
+2       FORMAT(I5,' : UNKNOWN BOTTOM FRICTION LAW')
         CALL PLANTE(1)
         STOP
 !
diff --git a/sources/telemac2d/friction_cui.f b/sources/telemac2d/friction_cui.f
new file mode 100644
index 0000000000..aad8ecab47
--- /dev/null
+++ b/sources/telemac2d/friction_cui.f
@@ -0,0 +1,69 @@
+!                   ***********************
+                    SUBROUTINE FRICTION_CUI
+!                   ***********************
+!
+     &(HA,CD,A,HDVEG,ZETA,Y,ALFA,BETA,CPI,KARMAN,CP)
+!
+!***********************************************************************
+! TELEMAC2D
+!***********************************************************************
+!
+!brief    COMPUTES FRICTION COEFFICIENT FOR SUBMERGED VEGETATION
+!+        FROM PARAMETERS WITH CUI ET AL.(2023) APPROACH
+!
+!history  GABRIELE FARINA (UNIBS) and VITO BACCHI (LNHE)
+!+        04/03/2024
+!+        V9.1
+!+
+!+   The algorithm was developed by
+!
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!| CD             |-->| BULK DRAG COEFFICIENT FOR VEGETATION
+!| CP             |<--| VEGETATION FRICTION COEFFICIENT
+!| HA             |-->| WATER DEPTH
+!| HDVEG          |-->| HEIGHT OF DEFLECTED VEGETATION
+!| KARMAN         |-->| VON KARMAN CONSTANT
+!| A              |-->| FRONTAL AREA PER UNIT VOLUME
+!| ZETA           |-->| ADIMENSIONAL PARAMETER ui/uUD
+!| Y              |-->| ADIMENSIONAL PARAMETER y0/HDVEG
+!| ALFA           |-->| ADIMENSIONAL PARAMETER Le/HDVEG
+!| BETA           |-->| ADIMENSIONAL PARAMETER yi/HDVEG
+!| CPI            |-->| WAKE FUNCTION COEFFICIENT
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!
+      USE INTERFACE_TELEMAC2D, EX_FRICTION_CUI => FRICTION_CUI
+      IMPLICIT NONE
+!
+!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
+!
+      DOUBLE PRECISION, INTENT(IN)  :: HA, CD, A, HDVEG, ZETA, Y
+      DOUBLE PRECISION, INTENT(IN)  :: ALFA,BETA,CPI, KARMAN
+      DOUBLE PRECISION, INTENT(OUT) :: CP
+!
+!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
+!
+      DOUBLE PRECISION TERM1,TERM2,TERM3,TERM4,TERM5,TERM6,TERM7,TERM8,C
+!
+!-----------------------------------------------------------------------
+!
+      IF(MIN(HA,HDVEG).LE.0.D0) THEN
+        CP = 0.D0
+      ELSE IF (MIN(HA,HDVEG).GT.0.D0.AND.HA.LE.HDVEG) THEN
+        CP = CD*A*HA
+      ELSE
+        TERM1 = SQRT(2.D0/(CD*A))
+        TERM2 = (ALFA*HDVEG)/HA
+        TERM3 = LOG(COSH(BETA/ALFA-HA/(ALFA*HDVEG))/(COSH(BETA/ALFA)))
+        TERM4 = SQRT(HDVEG*(HA/HDVEG-1.D0))
+        TERM5 = (HA-BETA*HDVEG)/HA
+        TERM6 = (HA-BETA*HDVEG)/HDVEG
+        C = -1/KARMAN*LOG(Y)
+        TERM7 = TERM1*(1.D0+(ZETA-1.D0)*(1.D0+TERM2*TERM3))
+        TERM8=TERM4*(TERM5/KARMAN*(LOG(TERM6)+KARMAN*C-1.D0)+CPI/KARMAN)
+        CP = 2.D0*HA/((TERM7+TERM8)*(TERM7+TERM8))
+      ENDIF
+!
+!-----------------------------------------------------------------------
+!
+      RETURN
+      END
diff --git a/sources/telemac2d/friction_def.f b/sources/telemac2d/friction_def.f
index f97cdeeb1b..fe7dc2ed78 100644
--- a/sources/telemac2d/friction_def.f
+++ b/sources/telemac2d/friction_def.f
@@ -36,12 +36,12 @@
 !         SEQUENCE
         INTEGER          :: GNUMB(2) ! GLOBAL NUMBER OF THE ZONE
         INTEGER          :: RTYPE ! TYPE OF LAW USED
-        INTEGER          :: VTYPE
+        INTEGER          :: VTYPE ! TYPE OF LAW FOR VEGETATION
         ! USE REAL BECAUSE CHESTR IS SAVED AS SIMPLE PRECISION IN SELAFIN DATA
         ! --------------------------------------------------------------------
         DOUBLE PRECISION :: RCOEF ! FRICTION PARAMETER
         DOUBLE PRECISION :: NDEF  ! DEFAULT MANNING (FOR C-W LAW)
-        DOUBLE PRECISION :: VCOEF(15)    ! 15 COEFFICIENTS 
+        DOUBLE PRECISION :: VCOEF(15) ! VEGETATION COEF: 15 COEFFICIENTS 
         
         TYPE(POINTER_TO_FRICTION), POINTER, DIMENSION(:) :: ADR
       END TYPE FRICTION_OBJ
diff --git a/sources/telemac2d/friction_read.f b/sources/telemac2d/friction_read.f
index 9cfebe2529..b734973091 100644
--- a/sources/telemac2d/friction_read.f
+++ b/sources/telemac2d/friction_read.f
@@ -6,7 +6,7 @@
      & SB) 
 !
 !***********************************************************************
-! TELEMAC2D   V8P2
+! TELEMAC2D
 !***********************************************************************
 !
 !brief    FRICTION FILE READ.
@@ -375,6 +375,20 @@
      &                 FRTAB%ADR(IZONE)%P%NDEF ,
      &                 LAWVEG,
      &                 (FRTAB%ADR(IZONE)%P%VCOEF(J),J=1,5)
+        CASE ('VCUI')
+          READ(NCOF,*,END=999,ERR=888)
+     &        (CHAINE(J),J=1,N),LAWVEG,(R3(J),J=1,8)
+          FRTAB%ADR(IZONE)%P%VTYPE = 10
+          DO J=1,8
+            FRTAB%ADR(IZONE)%P%VCOEF(J) = R3(J)
+          ENDDO
+          WRITE(LU,11) FRTAB%ADR(IZONE)%P%GNUMB(1),
+     &                 FRTAB%ADR(IZONE)%P%GNUMB(2),
+     &                 LAW,
+     &                 FRTAB%ADR(IZONE)%P%RCOEF,
+     &                 FRTAB%ADR(IZONE)%P%NDEF ,
+     &                 LAWVEG,
+     &                 (FRTAB%ADR(IZONE)%P%VCOEF(J),J=1,8)
         CASE DEFAULT
           WRITE(LU,10) LINE, IZ1, IZ2, LAWVEG
           CALL PLANTE(1)
diff --git a/sources/telemac2d/friction_zones.f b/sources/telemac2d/friction_zones.f
index 8b936884f8..42a99d6ee6 100644
--- a/sources/telemac2d/friction_zones.f
+++ b/sources/telemac2d/friction_zones.f
@@ -7,7 +7,7 @@
      & KARMAN, GRAV, T1, T2, CF, CFBOR)
 !
 !***********************************************************************
-! TELEMAC2D   V8P4
+! TELEMAC2D
 !***********************************************************************
 !
 !brief    COMPUTES FRICTION FOR EACH NODE AND ZONE.
@@ -214,10 +214,18 @@
      &       VCOEFF%ADR(4)%P%R(I), VCOEFF%ADR(5)%P%R(I),
      &       KARMAN,CP)
           CF%R(I) =CF%R(I)+CP
+        CASE (10)
+          CALL FRICTION_CUI
+     &       (T1%R(I), VCOEFF%ADR(1)%P%R(I),
+     &       VCOEFF%ADR(2)%P%R(I), VCOEFF%ADR(3)%P%R(I),
+     &       VCOEFF%ADR(4)%P%R(I), VCOEFF%ADR(5)%P%R(I),
+     &       VCOEFF%ADR(6)%P%R(I), VCOEFF%ADR(7)%P%R(I),
+     &       VCOEFF%ADR(8)%P%R(I),KARMAN,CP)
+          CF%R(I) =CP
         CASE DEFAULT
 !
           WRITE(LU,2) VEGLAW%I(I)
- 2             FORMAT(I5,' : UNKNOWN FRICTION LAW')
+ 2             FORMAT(I5,' : UNKNOWN VEGETATION FRICTION LAW ')
         CALL PLANTE(1)
         STOP
         END SELECT
diff --git a/sources/telemac2d/interface_telemac2d.f b/sources/telemac2d/interface_telemac2d.f
index da96f897c7..4b71e9f060 100644
--- a/sources/telemac2d/interface_telemac2d.f
+++ b/sources/telemac2d/interface_telemac2d.f
@@ -1163,6 +1163,18 @@
       END INTERFACE
 !
 !-----------------------------------------------------------------------
+!
+      INTERFACE
+        SUBROUTINE FRICTION_CUI(HA,CD,A,HDVEG,ZETA,Y,ALFA,BETA,
+     &       CPI,KARMAN,CP)
+      IMPLICIT NONE
+      DOUBLE PRECISION, INTENT(IN)  :: HA,CD,A,HDVEG,ZETA,Y,ALFA
+      DOUBLE PRECISION, INTENT(IN)  :: BETA,CPI,KARMAN
+      DOUBLE PRECISION, INTENT(OUT) :: CP
+        END SUBROUTINE
+      END INTERFACE
+!
+!-----------------------------------------------------------------------
 !
       INTERFACE
         SUBROUTINE FRICTION_READ
diff --git a/sources/telemac2d/telemac2d.cmdf b/sources/telemac2d/telemac2d.cmdf
index 58d5200407..ceeeafe72e 100644
--- a/sources/telemac2d/telemac2d.cmdf
+++ b/sources/telemac2d/telemac2d.cmdf
@@ -971,6 +971,7 @@ files: friction_def.f
   friction_calc.f
   friction_unif.f
   friction_baptist.f
+  friction_cui.f
   friction_huthoff.f
   friction_hybrid.f
   friction_jaervelae.f
-- 
GitLab