Index: NEMO/branches/2020/r12581_ticket2418/src/OCE/SBC/sbcblk_phy.F90
===================================================================
--- NEMO/branches/2020/r12581_ticket2418/src/OCE/SBC/sbcblk_phy.F90	(revision 12582)
+++ NEMO/branches/2020/r12581_ticket2418/src/OCE/SBC/sbcblk_phy.F90	(revision 12623)
@@ -520,9 +520,9 @@
          zCe = zz0*pqst(ji,jj)/zdq
 
-         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),                 &
-            &              pTau(ji,jj), zQsen, zQlat )
-
+         CALL BULK_FORMULA_SCLR( 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),                &
+            &                    pTau(ji,jj), zQsen, zQlat )
+         
          zTs2  = pTs(ji,jj)*pTs(ji,jj)
          zQlw  = emiss_w*(prlw(ji,jj) - stefan*zTs2*zTs2) ! Net longwave flux
@@ -535,8 +535,63 @@
 
 
-   SUBROUTINE BULK_FORMULA_VCTR( pzu, pTs, pqs, pTa, pqa,  &
-      &                          pCd, pCh, pCe,            &
-      &                          pwnd, pUb, pslp,          &
-      &                          pTau, pQsen, pQlat,  pEvap, prhoa )
+   SUBROUTINE BULK_FORMULA_SCLR( pzu, pTs, pqs, pTa, pqa, &
+      &                          pCd, pCh, pCe,           &
+      &                          pwnd, pUb, pslp,         &
+      &                          pTau, pQsen, pQlat,      &
+      &                          pEvap, prhoa, pfact_evap )
+      !!----------------------------------------------------------------------------------
+      REAL(wp),                     INTENT(in)  :: pzu  ! height above the sea-level where all this takes place (normally 10m)
+      REAL(wp), INTENT(in)  :: pTs  ! water temperature at the air-sea interface [K]
+      REAL(wp), INTENT(in)  :: pqs  ! satur. spec. hum. at T=pTs   [kg/kg]
+      REAL(wp), INTENT(in)  :: pTa  ! potential air temperature at z=pzu [K]
+      REAL(wp), INTENT(in)  :: pqa  ! specific humidity at z=pzu [kg/kg]
+      REAL(wp), INTENT(in)  :: pCd
+      REAL(wp), INTENT(in)  :: pCh
+      REAL(wp), INTENT(in)  :: pCe
+      REAL(wp), INTENT(in)  :: pwnd ! wind speed module at z=pzu [m/s]
+      REAL(wp), INTENT(in)  :: pUb  ! bulk wind speed at z=pzu (inc. pot. effect of gustiness etc) [m/s]
+      REAL(wp), INTENT(in)  :: pslp ! sea-level atmospheric pressure [Pa]
+      !!
+      REAL(wp), INTENT(out) :: pTau  ! module of the wind stress [N/m^2]
+      REAL(wp), INTENT(out) :: pQsen !  [W/m^2]
+      REAL(wp), INTENT(out) :: pQlat !  [W/m^2]
+      !!
+      REAL(wp), INTENT(out), OPTIONAL :: pEvap ! Evaporation [kg/m^2/s]
+      REAL(wp), INTENT(out), OPTIONAL :: prhoa ! Air density at z=pzu [kg/m^3]
+      REAL(wp), INTENT(in) , OPTIONAL :: pfact_evap  ! ABOMINATION: corrective factor for evaporation (doing this against my will! /laurent)
+      !!
+      REAL(wp) :: ztaa, zgamma, zrho, zUrho, zevap, zfact_evap
+      INTEGER  :: jq
+      !!----------------------------------------------------------------------------------
+      zfact_evap = 1._wp
+      IF( PRESENT(pfact_evap) ) zfact_evap = pfact_evap
+      
+      !! 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 )  !LOLO: why not "0.5*(pqs+pqa)" rather then "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 = - zfact_evap * zevap
+      IF( PRESENT(prhoa) ) prhoa = zrho
+
+   END SUBROUTINE BULK_FORMULA_SCLR
+
+   SUBROUTINE BULK_FORMULA_VCTR( pzu, pTs, pqs, pTa, pqa, &
+      &                          pCd, pCh, pCe,           &
+      &                          pwnd, pUb, pslp,         &
+      &                          pTau, pQsen, pQlat,      & 
+      &                          pEvap, prhoa, pfact_evap )      
       !!----------------------------------------------------------------------------------
       REAL(wp),                     INTENT(in)  :: pzu  ! height above the sea-level where all this takes place (normally 10m)
@@ -558,85 +613,25 @@
       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_2D_11_11
-
-         !! 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
+      REAL(wp),                     INTENT(in) , OPTIONAL :: pfact_evap  ! ABOMINATION: corrective factor for evaporation (doing this against my will! /laurent)
+      !!
+      REAL(wp) :: ztaa, zgamma, zrho, zUrho, zevap, zfact_evap
+      INTEGER  :: ji, jj
+      !!----------------------------------------------------------------------------------
+      zfact_evap = 1._wp
+      IF( PRESENT(pfact_evap) ) zfact_evap = pfact_evap
+
+      DO_2D_11_11
+
+         CALL BULK_FORMULA_SCLR( pzu, pTs(ji,jj), pqs(ji,jj), pTa(ji,jj), pqa(ji,jj), &
+            &                    pCd(ji,jj), pCh(ji,jj), pCe(ji,jj),                  &
+            &                    pwnd(ji,jj), pUb(ji,jj), pslp(ji,jj),                &
+            &                    pTau(ji,jj), pQsen(ji,jj), pQlat(ji,jj),             &
+            &                    pEvap=zevap, prhoa=zrho, pfact_evap=zfact_evap       )
+
+         IF( PRESENT(pEvap) ) pEvap(ji,jj) = zevap
          IF( PRESENT(prhoa) ) prhoa(ji,jj) = zrho
-
+   
       END_2D
    END SUBROUTINE BULK_FORMULA_VCTR
-
-
-   SUBROUTINE BULK_FORMULA_SCLR( pzu, pTs, pqs, pTa, pqa, &
-      &                          pCd, pCh, pCe,           &
-      &                          pwnd, pUb, pslp,         &
-      &                          pTau, pQsen, pQlat,  pEvap, prhoa )
-      !!----------------------------------------------------------------------------------
-      REAL(wp),                     INTENT(in)  :: pzu  ! height above the sea-level where all this takes place (normally 10m)
-      REAL(wp), INTENT(in)  :: pTs  ! water temperature at the air-sea interface [K]
-      REAL(wp), INTENT(in)  :: pqs  ! satur. spec. hum. at T=pTs   [kg/kg]
-      REAL(wp), INTENT(in)  :: pTa  ! potential air temperature at z=pzu [K]
-      REAL(wp), INTENT(in)  :: pqa  ! specific humidity at z=pzu [kg/kg]
-      REAL(wp), INTENT(in)  :: pCd
-      REAL(wp), INTENT(in)  :: pCh
-      REAL(wp), INTENT(in)  :: pCe
-      REAL(wp), INTENT(in)  :: pwnd ! wind speed module at z=pzu [m/s]
-      REAL(wp), INTENT(in)  :: pUb  ! bulk wind speed at z=pzu (inc. pot. effect of gustiness etc) [m/s]
-      REAL(wp), INTENT(in)  :: pslp ! sea-level atmospheric pressure [Pa]
-      !!
-      REAL(wp), INTENT(out) :: pTau  ! module of the wind stress [N/m^2]
-      REAL(wp), INTENT(out) :: pQsen !  [W/m^2]
-      REAL(wp), INTENT(out) :: pQlat !  [W/m^2]
-      !!
-      REAL(wp), INTENT(out), OPTIONAL :: pEvap ! Evaporation [kg/m^2/s]
-      REAL(wp), INTENT(out), OPTIONAL :: prhoa ! Air density at z=pzu [kg/m^3]
-      !!
-      REAL(wp) :: ztaa, zgamma, zrho, zUrho, zevap
-      INTEGER  :: jq
-      !!----------------------------------------------------------------------------------
-
-      !! 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
-
-
 
 
