Index: /NEMO/branches/2019/dev_ASINTER-01-05_merged/cfgs/SHARED/namelist_ref
===================================================================
--- /NEMO/branches/2019/dev_ASINTER-01-05_merged/cfgs/SHARED/namelist_ref	(revision 12080)
+++ /NEMO/branches/2019/dev_ASINTER-01-05_merged/cfgs/SHARED/namelist_ref	(revision 12081)
@@ -251,5 +251,5 @@
 /
 !-----------------------------------------------------------------------
-&namsbc_blk    !   namsbc_blk  generic Bulk formula                     (ln_blk =T)
+&namsbc_blk    !   namsbc_blk  generic Bulk formula          (ln_blk =T)
 !-----------------------------------------------------------------------
    !                    !  bulk algorithm :
@@ -259,18 +259,20 @@
    ln_ECMWF     = .false.   ! "ECMWF"     algorithm   (IFS cycle 31)
       !
-      rn_zqt      = 10.       !  Air temperature & humidity reference height (m)
-      rn_zu       = 10.       !  Wind vector reference height (m)
-      ln_Cd_L12   = .false.   !  air-ice drags = F(ice concentration) (Lupkes et al. 2012)
-      ln_Cd_L15   = .false.   !  air-ice drags = F(ice concentration) (Lupkes et al. 2015)
-      rn_pfac     = 1.        !  multiplicative factor for precipitation (total & snow)
-      rn_efac     = 1.        !  multiplicative factor for evaporation (0. or 1.)
-      rn_vfac     = 0.        !  multiplicative factor for ocean & ice velocity used to
-      !                       !  calculate the wind stress (0.=absolute or 1.=relative winds)
-      ln_skin_cs = .false.  !  use the cool-skin parameterization (only available in ECMWF and COARE algorithms) !LB
-      ln_skin_wl = .false.  !  use the warm-layer        "               "                    "
-      !
-      ln_humi_sph = .true.     !  humidity specified below in "sn_humi" is specific humidity     [kg/kg] if .true.
-      ln_humi_dpt = .false.    !  humidity specified below in "sn_humi" is dew-point temperature   [K]   if .true.
-      ln_humi_rlh = .false.    !  humidity specified below in "sn_humi" is relative humidity       [%]   if .true.
+      rn_zqt     = 10.      !  Air temperature & humidity reference height (m)
+      rn_zu      = 10.      !  Wind vector reference height (m)
+      ln_Cd_L12  = .false.  !  air-ice drags = F(ice conc.) (Lupkes et al. 2012)
+      ln_Cd_L15  = .false.  !  air-ice drags = F(ice conc.) (Lupkes et al. 2015)
+      !                     !  - module of the mean stress" data
+      rn_pfac    = 1.       !  multipl. factor for precipitation (total & snow)
+      rn_efac    = 1.       !  multipl. factor for evaporation (0. or 1.)
+      rn_vfac    = 0.       !  multipl. factor for ocean & ice velocity 
+      !                     !  used to calculate the wind stress
+      !                     ! (0. => absolute or 1. => relative winds)
+      ln_skin_cs = .false.  !  use the cool-skin parameterization
+      ln_skin_wl = .false.  !  use the warm-layer parameterization
+      !                     !   ==> only available in ECMWF and COARE algorithms
+      ln_humi_sph = .true.  !  humidity "sn_humi" is specific humidity  [kg/kg]
+      ln_humi_dpt = .false. !  humidity "sn_humi" is dew-point temperature [K]
+      ln_humi_rlh = .false. !  humidity "sn_humi" is relative humidity     [%]
    !
    cn_dir      = './'      !  root directory for the bulk data location
Index: /NEMO/branches/2019/dev_ASINTER-01-05_merged/doc/latex/NEMO/main/bibliography.bib
===================================================================
--- /NEMO/branches/2019/dev_ASINTER-01-05_merged/doc/latex/NEMO/main/bibliography.bib	(revision 12080)
+++ /NEMO/branches/2019/dev_ASINTER-01-05_merged/doc/latex/NEMO/main/bibliography.bib	(revision 12081)
@@ -188,4 +188,18 @@
 }
 
+@article{ beljaars_QJRMS95,
+title = "The parametrization of surface fluxes in large-scale models under free convection",
+pages = "255--270",
+journal = "Quarterly Journal of the Royal Meteorological Society",
+volume = "121",
+number = "522",
+author = "Beljaars, Anton C. M.",
+year = "1995",
+month = "jan",
+ publisher     = "Wiley",
+issn = "00359009",
+doi = "10.1002/qj.49712152203"
+}
+
 @article{         bernie.guilyardi.ea_CD07,
   title         = "Impact of resolving the diurnal cycle in an
@@ -386,4 +400,18 @@
 }
 
+@article{         brodeau.barnier.ea_JPO17,
+  title         = "Climatologically Significant Effects of Some Approximations in the Bulk Parameterizations of Turbulent Air{\textendash}Sea Fluxes",
+  pages         = "5--28",
+  journal       = "Journal of Physical Oceanography",
+  volume        = "47",
+  number        = "1",
+  author        = "Brodeau, Laurent and Barnier, Bernard and Gulev, Sergey K. and Woods, Cian",
+  year          = "2017",
+  month         = "jan",
+  publisher     = "American Meteorological Society",
+  issn          = "0022-3670",
+  doi           = "10.1175/jpo-d-16-0169.1",
+}
+
 @article{         brown.campana_MWR78,
   title         = "An economical time-differencing system for numerical
@@ -745,4 +773,18 @@
 }
 
+@article{ edson.jampana.ea_JPO13,
+title = "On the Exchange of Momentum over the Open Ocean",
+pages = "1589--1610",
+journal = "Journal of Physical Oceanography",
+volume = "43",
+number = "8",
+author = "Edson, James B. and Jampana, Venkata and Weller, Robert A. and Bigorre, Sebastien P. and Plueddemann, Albert J. and Fairall, Christopher W. and Miller, Scott D. and Mahrt, Larry and Vickers, Dean and Hersbach, Hans",
+year = "2013",
+month = "aug",
+publisher = "American Meteorological Society",
+issn = "0022-3670",
+doi = "10.1175/JPO-D-12-0173.1"
+}
+
 @article{         egbert.ray_JGR01,
   title         = "Estimates of {M2} tidal energy dissipation from
@@ -819,4 +861,17 @@
 }
 
+@article{ fairall.bradley.ea_JC03,
+title = "Bulk parameterization of air-sea fluxes: Updates and verification for the COARE algorithm",
+pages = "571--591",
+journal = "Journal of Climate",
+volume = "16",
+number = "4",
+author = "Fairall, C. W. and Bradley, E. F. and Hare, J. E. and Grachev, A. A. and Edson, J. B.",
+year = "2003",
+publisher = "American Meteorological Society",
+issn = "08948755",
+doi = "10.1175/1520-0442(2003)016<0571:BPOASF>2.0.CO;2"
+}
+
 @phdthesis{       farge-coulombier_phd87,
   title         = "Dynamique non-lin\'{e}aire des ondes et des tourbillons
@@ -871,4 +926,18 @@
   publisher     = "UNESCO",
   url           = "https://unesdoc.unesco.org/ark:/48223/pf0000059832.locale=en"
+}
+
+@article{         foxkemper.ferrari_JPO08,
+  title         = "Parameterization of Mixed Layer Eddies. Part I: Theory and Diagnosis",
+  pages         = "1145--1165",
+  journal       = "Journal of Physical Oceanography",
+  volume        = "38",
+  number        = "6",
+  author        = "B. Fox-Kemper and R. Ferrari and B. Hallberg",
+  year          = "2008",
+  month         = "jun",
+  publisher     = "American Meteorological Society",
+  issn          = "1520-0485",
+  doi           = "10.1175/2007JPO3792.1"
 }
 
@@ -1545,4 +1614,24 @@
 }
 
+@article{         kraus.turner_tellus67,
+  author = {Kraus, E.B. and Turner, J.},
+  journal = {Tellus},
+  pages = {98--106},
+  title = {A one dimensional model of the seasonal thermocline {II}. {T}he general theory and its consequences},
+  volume = {19},
+  year = {1967}
+}
+
+@article{        large.ea_RG97,
+  author        = "Large, W. G. and McWilliams, J. C. and Doney, S. C.",
+  doi           = "10.1029/94RG01872",
+  journal       = "Reviews of Geophysics",
+  number        = {4},
+  pages         = {363--403},
+  publisher     = {AGU},
+  title         = "Oceanic vertical mixing: {A} review and a model with a nonlocal boundary layer parameterization",
+  year          = "1994"
+}
+
 @techreport{      large.yeager_rpt04,
   title         = "Diurnal to decadal global forcing for ocean and sea-ice
@@ -1839,4 +1928,28 @@
   author        = "F. Lott and G. Madec and J. Verron",
   year          = "1990"
+}
+
+@article{	  lupkes.gryanik.ea_JGR12,
+  author 	= "L{\"{u}}pkes, Christof and Gryanik, Vladimir M. and Hartmann, J{\"{o}}rg and Andreas, Edgar L.",
+  doi 		= "10.1029/2012JD017630",
+  issn 		= "01480227",
+  journal 	= "Journal of Geophysical Research Atmospheres",
+  number 	= "13",
+  pages 	= "1--18",
+  title 	= "A parametrization, based on sea ice morphology, of the neutral atmospheric drag coefficients for weather prediction and climate models",
+  volume 	= "117",
+  year 		= "2012"
+}
+
+@article{	  lupkes.gryanik_JGR15,
+  author 	= "L{\"{u}}pkes, Christof and Gryanik, Vladimir M.",
+  doi 		= "10.1002/2014JD022418",
+  issn 		= "21562202",
+  journal 	= "Journal of Geophysical Research",
+  number 	= "2",
+  pages 	= "552--581",
+  title 	= "A stability-dependent parametrization of transfer coefficients formomentum and heat over polar sea ice to be used in climate models",
+  volume 	= "120",
+  year 		= "2015"
 }
 
@@ -2091,4 +2204,13 @@
 }
 
+@article{mcwilliams.ea_JFM97,
+	author = {McWilliams, James C. and Sullivan, Peter P. and Moeng, Chin-Hoh},
+	doi = {10.1017/S0022112096004375},
+	journal = {Journal of Fluid Mechanics},
+	pages = {1--30},
+	title = {Langmuir turbulence in the ocean},
+	volume = {334},
+	year = {1997},
+}
 @article{         mellor.blumberg_JPO04,
   title         = "Wave Breaking and Ocean Surface Layer Thermal Response",
@@ -2304,4 +2426,20 @@
   issn          = "0002-9947",
   doi           = "10.1090/s0002-9947-1959-0104919-1"
+}
+
+@article{         qiao.yuan.ea_OD10,
+  title         = "A three-dimensional surface wave–ocean circulation coupled 
+                  model and its initial testing",
+  pages         = "1339--1335",
+  journal       = "Ocean Dynamics",
+  volume        = "60",
+  number        = "5",
+  author        = "F. Qiao and Y. Yuan and T. Ezer and C. Xia and 
+                   Y. Yang and X. Lu and Z. Song ",
+  year          = "2010",
+  month         = "oct",
+  publisher     = "Springer-Verlag",
+  issn          = "1616-7341",
+  doi           = "10.1007/s10236-010-0326-y"
 }
 
@@ -2544,4 +2682,18 @@
 }
 
+@article{         smagorinsky_MW63,
+  title         = "General circulation experiments with the primitive equations: I. The basic experiment ",
+  pages         = "99--164",
+  journal       = "Monthly Weather Review",
+  volume        = "91",
+  number        = "3",
+  author        = "J. Smagorinsky",
+  year          = "1963",
+  month         = "mar",
+  publisher     = "American Meteorological Society",
+  issn          = "1520-0493",
+  doi           = "10.1175/1520-0493(1963)091<0099:GCEWTP>2.3.CO;2"
+}
+
 @article{         song.haidvogel_JCP94,
   title         = "A semi-implicit ocean circulation model using a
@@ -2896,4 +3048,15 @@
 }
 
+@article{         white.hoskins.ea_QJRMS05,
+  title         = "Consistent approximate models of the global atmosphere: shallow, deep,
+                  hydrostatic, quasi-hydrostatic and non-hydrostatic",
+  pages         = "2081--2107",
+  journal       = "Quarterly Journal of the Royal Meteorological Society",
+  volume        = "131",
+  author        = "A. A. White and B. J. Hoskins and I. Roulstone and A. Staniforth",
+  year          = "2005",
+  doi           = "10.1256/qj.04.49"
+}
+
 @article{         white.adcroft.ea_JCP09,
   title         = "High-order regridding-remapping schemes for continuous
@@ -2971,2 +3134,73 @@
   doi           = "10.1029/92jc00911"
 }
+
+@article{large.yeager_CD09,
+author="Large, W. G. and Yeager, S. G.",
+title="The Global Climatology of an Interannually Varying Air-Sea Flux Data Set",
+pages = "341--364",
+journal="Climate Dynamics",
+volume = "33",
+number = "2-3",
+year="2009",
+month = "aug",
+publisher = "Springer Science and Business Media LLC",
+doi="10.1007/s00382-008-0441-3"
+}
+
+@book{sverdrup.johnson.ea_1942,
+author = {H. U. Sverdrup and Martin W. Johnson and Richard H. Fleming},
+title = {The Oceans, Their Physics, Chemistry, and General Biology},
+publisher = {Prentice-Hall},
+address = {New York},
+year = {1942},
+pages = {1087},
+}
+
+@article{kraus.businger_QJRMS96,
+author = "E. B. Kraus and J. A. Businger",
+title = "Atmosphere-ocean interaction.",
+journal="Quarterly Journal of the Royal Meteorological Society",,
+year = "1996",
+volume = "122",
+number = "529",
+pages = "324-325",
+publisher = "John Wiley & Sons, Ltd",
+issn = "1477-870X",
+doi = "10.1002/qj.49712252914"
+}
+
+@article{josey.gulev.ea_2013,
+title = "Exchanges through the ocean surface",
+journal = "Ocean Circulation and Climate - A 21st Century Perspective, Int. Geophys. Ser.",
+year = "2013",
+author = "S. A. Josey and S. Gulev and L. Yu",
+pages = "115-140, edited by G. Siedler et al., Academic Press, Oxford",
+volume = "103",
+doi = "10.1016/B978-0-12-391851-2.00005-2"
+}
+
+@article{fairall.bradley.ea_JGR96,
+  year = "1996",
+ journal = "Journal of Geophysical Research: Oceans",
+  month = "jan",
+  publisher = "American Geophysical Union",
+  volume = "101",
+  number = "C1",
+  pages = "1295-1308",
+  author = "C. W. Fairall and E. F. Bradley and J. S. Godfrey and G. A. Wick and J. B. Edson and G. S. Young",
+  title = "Cool-skin and warm-layer effects on sea surface temperature",
+  doi = "10.1029/95jc03190"
+}
+
+@article{zeng.beljaars_GRL05,
+  year = "2005",
+  month = "jul",
+  publisher = "American Geophysical Union",
+  volume = "32",
+  number = "14",
+  author = "Xubin Zeng and Anton Beljaars",
+  title = "A prognostic scheme of sea surface skin temperature for modeling and data assimilation",
+  journal = "Geophysical Research Letters",
+  doi = "10.1029/2005gl023030"
+}
+
Index: /NEMO/branches/2019/dev_ASINTER-01-05_merged/doc/namelists/namsbc_blk
===================================================================
--- /NEMO/branches/2019/dev_ASINTER-01-05_merged/doc/namelists/namsbc_blk	(revision 12080)
+++ /NEMO/branches/2019/dev_ASINTER-01-05_merged/doc/namelists/namsbc_blk	(revision 12081)
@@ -3,5 +3,5 @@
 !-----------------------------------------------------------------------
    !                    !  bulk algorithm :
-   ln_NCAR      = .false.   ! "NCAR"      algorithm   (Large and Yeager 2008)
+   ln_NCAR      = .true.    ! "NCAR"      algorithm   (Large and Yeager 2008)
    ln_COARE_3p0 = .false.   ! "COARE 3.0" algorithm   (Fairall et al. 2003)
    ln_COARE_3p6 = .false.   ! "COARE 3.6" algorithm   (Edson et al. 2013)
@@ -12,5 +12,4 @@
       ln_Cd_L12  = .false.  !  air-ice drags = F(ice conc.) (Lupkes et al. 2012)
       ln_Cd_L15  = .false.  !  air-ice drags = F(ice conc.) (Lupkes et al. 2015)
-      ln_taudif  = .false.  !  HF tau contribution: use "mean of stress module 
       !                     !  - module of the mean stress" data
       rn_pfac    = 1.       !  multipl. factor for precipitation (total & snow)
@@ -36,7 +35,8 @@
    sn_tair     = 't_10.15JUNE2009_fill'       ,    6.        , 'T_10_MOD',   .false.   , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc' , ''       , ''
    sn_humi     = 'q_10.15JUNE2009_fill'       ,    6.        , 'Q_10_MOD',   .false.   , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc' , ''       , ''
+   sn_hpgi     = 'NONE'                       ,   24.        , 'uhpg'    ,   .false.   , .false., 'monthly' , 'weights_ERAI3D_F128_2_ORCA2_bicubic', 'UG'     , ''
+   sn_hpgj     = 'NONE'                       ,   24.        , 'vhpg'    ,   .false.   , .false., 'monthly' , 'weights_ERAI3D_F128_2_ORCA2_bicubic', 'VG'     , ''
    sn_prec     = 'ncar_precip.15JUNE2009_fill',   -1.        , 'PRC_MOD1',   .false.   , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc' , ''       , ''
    sn_snow     = 'ncar_precip.15JUNE2009_fill',   -1.        , 'SNOW'    ,   .false.   , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc' , ''       , ''
    sn_slp      = 'slp.15JUNE2009_fill'        ,    6.        , 'SLP'     ,   .false.   , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc' , ''       , ''
-   sn_tdif     = 'taudif_core'                ,   24         , 'taudif'  ,   .false.   , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc' , ''       , ''
 /
Index: /NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/SBC/sbcblk_phy.F90
===================================================================
--- /NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/SBC/sbcblk_phy.F90	(revision 12080)
+++ /NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/SBC/sbcblk_phy.F90	(revision 12081)
@@ -533,5 +533,5 @@
             CALL BULK_FORMULA( pzu, pTs(ji,jj), pqs(ji,jj), pTa(ji,jj), pqa(ji,jj), &
                &              zCd, zCh, zCe,                                        &
-               &              pwnd(ji,jj), pUb(ji,jj), pslp(ji,jj), &
+               &              pwnd(ji,jj), pUb(ji,jj), pslp(ji,jj),                 &
                &              pTau(ji,jj), zQsen, zQlat )
 
@@ -568,35 +568,35 @@
       REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pQlat !  [W/m^2]
       !!
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(out), OPTIONAL :: pEvap  ! Evaporation [kg/m^2/s]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(out), OPTIONAL :: prhoa  ! Air density at z=pzu [kg/m^3]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(out), OPTIONAL :: pEvap ! Evaporation [kg/m^2/s]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(out), OPTIONAL :: prhoa ! Air density at z=pzu [kg/m^3]
       !!
       REAL(wp) :: ztaa, zgamma, zrho, zUrho, zevap
       INTEGER  :: ji, jj, jq     ! dummy loop indices
       !!----------------------------------------------------------------------------------
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-   
-               !! Need ztaa, absolute temperature at pzu (formula to estimate rho_air needs absolute temperature, not the potential temperature "pTa")
-               ztaa = pTa(ji,jj) ! first guess...
-               DO jq = 1, 4
-                  zgamma = gamma_moist( 0.5*(ztaa+pTs(ji,jj)) , pqa(ji,jj) )
-                  ztaa = pTa(ji,jj) - zgamma*pzu   ! Absolute temp. is slightly colder...
-               END DO
-               zrho = rho_air(ztaa, pqa(ji,jj), pslp(ji,jj))
-               zrho = rho_air(ztaa, pqa(ji,jj), pslp(ji,jj)-zrho*grav*pzu) ! taking into account that we are pzu m above the sea level where SLP is given!
-   
-               zUrho = pUb(ji,jj)*MAX(zrho, 1._wp)     ! rho*U10
-   
-               pTau(ji,jj) = zUrho * pCd(ji,jj) * pwnd(ji,jj) ! Wind stress module
-   
-               zevap        = zUrho * pCe(ji,jj) * (pqa(ji,jj) - pqs(ji,jj))
-               pQsen(ji,jj) = zUrho * pCh(ji,jj) * (pTa(ji,jj) - pTs(ji,jj)) * cp_air(pqa(ji,jj))
-               pQlat(ji,jj) = L_vap(pTs(ji,jj)) * zevap
-   
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+
+            !! Need ztaa, absolute temperature at pzu (formula to estimate rho_air needs absolute temperature, not the potential temperature "pTa")
+            ztaa = pTa(ji,jj) ! first guess...
+            DO jq = 1, 4
+               zgamma = gamma_moist( 0.5*(ztaa+pTs(ji,jj)) , pqa(ji,jj) )
+               ztaa = pTa(ji,jj) - zgamma*pzu   ! Absolute temp. is slightly colder...
+            END DO
+            zrho = rho_air(ztaa, pqa(ji,jj), pslp(ji,jj))
+            zrho = rho_air(ztaa, pqa(ji,jj), pslp(ji,jj)-zrho*grav*pzu) ! taking into account that we are pzu m above the sea level where SLP is given!
+
+            zUrho = pUb(ji,jj)*MAX(zrho, 1._wp)     ! rho*U10
+
+            pTau(ji,jj) = zUrho * pCd(ji,jj) * pwnd(ji,jj) ! Wind stress module
+
+            zevap        = zUrho * pCe(ji,jj) * (pqa(ji,jj) - pqs(ji,jj))
+            pQsen(ji,jj) = zUrho * pCh(ji,jj) * (pTa(ji,jj) - pTs(ji,jj)) * cp_air(pqa(ji,jj))
+            pQlat(ji,jj) = L_vap(pTs(ji,jj)) * zevap
+
             IF( PRESENT(pEvap) ) pEvap(ji,jj) = - zevap
             IF( PRESENT(prhoa) ) prhoa(ji,jj) = zrho
-   
-            END DO
-         END DO
+
+         END DO
+      END DO
    END SUBROUTINE BULK_FORMULA_VCTR
 
@@ -630,23 +630,23 @@
       !!----------------------------------------------------------------------------------
 
-         !! Need ztaa, absolute temperature at pzu (formula to estimate rho_air needs absolute temperature, not the potential temperature "pTa")
-         ztaa = pTa ! first guess...
-         DO jq = 1, 4
-            zgamma = gamma_moist( 0.5*(ztaa+pTs) , pqa )
-            ztaa = pTa - zgamma*pzu   ! Absolute temp. is slightly colder...
-         END DO
-         zrho = rho_air(ztaa, pqa, pslp)
-         zrho = rho_air(ztaa, pqa, pslp-zrho*grav*pzu) ! taking into account that we are pzu m above the sea level where SLP is given!
-   
-         zUrho = pUb*MAX(zrho, 1._wp)     ! rho*U10
-   
-         pTau = zUrho * pCd * pwnd ! Wind stress module
-   
-         zevap = zUrho * pCe * (pqa - pqs)
-         pQsen = zUrho * pCh * (pTa - pTs) * cp_air(pqa)
-         pQlat = L_vap(pTs) * zevap
-   
-         IF( PRESENT(pEvap) ) pEvap = - zevap
-         IF( PRESENT(prhoa) ) prhoa = zrho
+      !! Need ztaa, absolute temperature at pzu (formula to estimate rho_air needs absolute temperature, not the potential temperature "pTa")
+      ztaa = pTa ! first guess...
+      DO jq = 1, 4
+         zgamma = gamma_moist( 0.5*(ztaa+pTs) , pqa )
+         ztaa = pTa - zgamma*pzu   ! Absolute temp. is slightly colder...
+      END DO
+      zrho = rho_air(ztaa, pqa, pslp)
+      zrho = rho_air(ztaa, pqa, pslp-zrho*grav*pzu) ! taking into account that we are pzu m above the sea level where SLP is given!
+
+      zUrho = pUb*MAX(zrho, 1._wp)     ! rho*U10
+
+      pTau = zUrho * pCd * pwnd ! Wind stress module
+
+      zevap = zUrho * pCe * (pqa - pqs)
+      pQsen = zUrho * pCh * (pTa - pTs) * cp_air(pqa)
+      pQlat = L_vap(pTs) * zevap
+
+      IF( PRESENT(pEvap) ) pEvap = - zevap
+      IF( PRESENT(prhoa) ) prhoa = zrho
 
    END SUBROUTINE BULK_FORMULA_SCLR
