ecumev6_flux.F90

Go to the documentation of this file.

!SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
!SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
!SFX_LIC for details. version 1.
!SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
!SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
!SFX_LIC for details. version 1.
!     #########
    SUBROUTINE ecumev6_flux(PZ0SEA,PTA,PEXNA,PRHOA,PSST,PSSS,PEXNS,PQA,PVMOD, &
                            PZREF,PUREF,PPS,PPA,PICHCE,OPRECIP,OPWEBB,        &
                            PQSAT,PSFTH,PSFTQ,PUSTAR,PCD,PCDN,PCH,PCE,        &
                            PRI,PRESA,PRAIN,KZ0,PZ0HSEA,OPERTFLUX,PPERTFLUX   )
!###############################################################################
!!
!!****  *ECUMEV6_FLUX*
!!
!!    PURPOSE
!!    -------
!       Calculate the surface turbulent fluxes of heat, moisture, and momentum 
!       over sea surface + corrections due to rainfall & Webb effect.
!!
!!**  METHOD
!!    ------
!       The estimation of the transfer coefficients relies on the iterative 
!       computation of the scaling parameters U*/Teta*/q*. The convergence is
!       supposed to be reached in NITERFL iterations maximum.
!       Neutral transfer coefficients for momentum/temperature/humidity
!       are computed as a function of the 10m-height neutral wind speed using
!       the ECUME_V6 formulation based on the multi-campaign (POMME,FETCH,CATCH,
!       SEMAPHORE,EQUALANT) ALBATROS dataset.
!!
!!    EXTERNAL
!!    --------
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!    REFERENCE
!!    ---------
!!      Fairall et al (1996), JGR, 3747-3764
!!      Gosnell et al (1995), JGR, 437-442
!!      Fairall et al (1996), JGR, 1295-1308
!!
!!    AUTHOR
!!    ------
!!      C. Lebeaupin  *Météo-France* (adapted from S. Belamari's code)
!!
!!    MODIFICATIONS
!!    -------------
!!      Original     15/03/2005
!!      Modified        01/2006  C. Lebeaupin (adapted from  A. Pirani's code)
!!      Modified        08/2009  B. Decharme: limitation of Ri
!!      Modified        09/2012  B. Decharme: CD correction
!!      Modified        09/2012  B. Decharme: limitation of Ri in surface_ri.F90
!!      Modified        10/2012  P. Le Moigne: extra inputs for FLake use
!!      Modified        06/2013  B. Decharme: bug in z0 (output) computation 
!!      Modified        12/2013  S. Belamari: ZRF computation updated:
!!                                1. ZP00/PPA in ZDWAT, ZLVA in ZDQSDT/ZBULB/ZRF
!!                                2. ZDWAT/ZDTMP in ZBULB/ZRF (Gosnell et al 95)
!!                                3. cool skin correction included
!!      Modified        01/2014  S. Belamari: salinity impact on latent heat of
!!                                vaporization of seawater included
!!      Modified        01/2014  S. Belamari: new formulation for pure water
!!                                specific heat (ZCPWA)
!!      Modified        01/2014  S. Belamari: 4 choices for PZ0SEA computation
!!      Modified        12/2015  S. Belamari: ECUME now provides parameterisations
!!                                for:  U10n*sqrt(CDN)          instead of CDN
!!                                      U10n*CHN/sqrt(CDN)         "       CHN
!!                                      U10n*CEN/sqrt(CDN)         "       CEN
!!      Modified        01/2016  S. Belamari: New ECUME formulation
!!
!!      To be done:
!!      include gustiness computation following Mondon & Redelsperger (1998)
!!!
!-------------------------------------------------------------------------------
!!
!!    MODIFICATIONS RELATED TO SST CORRECTION COMPUTATION
!!    ---------------------------------------------------
!!      Modified        09/2013  S. Belamari: use 0.98 for the ocean emissivity
!!                                following up to date satellite measurements in
!!                                the 8-14 μm range (obtained values range from
!!                                0.98 to 0.99).
!!!
!-------------------------------------------------------------------------------
!
!       0.   DECLARATIONS
!            ------------
!
USE modd_csts,             ONLY : xpi, xday, xkarman, xg, xp00, xstefan, xrd, xrv,   &
                                  xcpd, xcpv, xcl, xtt, xlvtt
USE modd_surf_par,         ONLY : xundef
USE modd_surf_atm,         ONLY : xvchrnk, xvz0cm
USE modd_reprod_oper,      ONLY : ccharnock
!
USE mode_thermos
USE modi_wind_threshold
USE modi_surface_ri
!
USE yomhook,   ONLY : lhook,   dr_hook
USE parkind1,  ONLY : jprb
!
USE modi_abor1_sfx
!
IMPLICIT NONE
!
!       0.1. Declarations of arguments
!
REAL, DIMENSION(:), INTENT(IN)    :: PVMOD      ! module of wind at atm level (m/s)
REAL, DIMENSION(:), INTENT(IN)    :: PTA        ! air temperature at atm level (K)
REAL, DIMENSION(:), INTENT(IN)    :: PQA        ! air spec. hum. at atm level (kg/kg)
REAL, DIMENSION(:), INTENT(IN)    :: PPA        ! air pressure at atm level (Pa)
REAL, DIMENSION(:), INTENT(IN)    :: PRHOA      ! air density at atm level (kg/m3)
REAL, DIMENSION(:), INTENT(IN)    :: PEXNA      ! Exner function at atm level
REAL, DIMENSION(:), INTENT(IN)    :: PUREF      ! atm level for wind (m)
REAL, DIMENSION(:), INTENT(IN)    :: PZREF      ! atm level for temp./hum. (m)
REAL, DIMENSION(:), INTENT(IN)    :: PSSS       ! Sea Surface Salinity (g/kg)
REAL, DIMENSION(:), INTENT(IN)    :: PPS        ! air pressure at sea surface (Pa)
REAL, DIMENSION(:), INTENT(IN)    :: PEXNS      ! Exner function at sea surface
REAL, DIMENSION(:), INTENT(IN)    :: PPERTFLUX  ! stochastic flux perturbation pattern
! for correction
REAL,               INTENT(IN)    :: PICHCE    !
LOGICAL,            INTENT(IN)    :: OPRECIP   !
LOGICAL,            INTENT(IN)    :: OPWEBB    !
LOGICAL,            INTENT(IN)    :: OPERTFLUX
REAL, DIMENSION(:), INTENT(IN)    :: PRAIN     ! precipitation rate (kg/s/m2)
!
INTEGER,            INTENT(IN)    :: KZ0
!
REAL, DIMENSION(:), INTENT(INOUT) :: PSST       ! Sea Surface Temperature (K)
REAL, DIMENSION(:), INTENT(INOUT) :: PZ0SEA     ! roughness length over sea
REAL, DIMENSION(:), INTENT(OUT)   :: PZ0HSEA    ! heat roughness length over sea

! surface fluxes : latent heat, sensible heat, friction fluxes
REAL, DIMENSION(:), INTENT(OUT)   :: PUSTAR     ! friction velocity (m/s)
REAL, DIMENSION(:), INTENT(OUT)   :: PSFTH      ! heat flux (W/m2)
REAL, DIMENSION(:), INTENT(OUT)   :: PSFTQ      ! water flux (kg/m2/s)

! diagnostics
REAL, DIMENSION(:), INTENT(OUT)   :: PQSAT      ! sea surface spec. hum. (kg/kg)
REAL, DIMENSION(:), INTENT(OUT)   :: PCD        ! transfer coef. for momentum
REAL, DIMENSION(:), INTENT(OUT)   :: PCH        ! transfer coef. for temperature
REAL, DIMENSION(:), INTENT(OUT)   :: PCE        ! transfer coef. for humidity
REAL, DIMENSION(:), INTENT(OUT)   :: PCDN       ! neutral coef. for momentum
REAL, DIMENSION(:), INTENT(OUT)   :: PRI        ! Richardson number
REAL, DIMENSION(:), INTENT(OUT)   :: PRESA      ! aerodynamical resistance
!
!       0.2. Declarations of local variables
!
! specif SB
INTEGER, DIMENSION(SIZE(PTA))     :: JCV        ! convergence index
INTEGER, DIMENSION(SIZE(PTA))     :: JITER      ! nb of iterations to converge

REAL, DIMENSION(SIZE(PTA))        :: ZTAU       ! momentum flux (N/m2)
REAL, DIMENSION(SIZE(PTA))        :: ZHF        ! sensible heat flux (W/m2)
REAL, DIMENSION(SIZE(PTA))        :: ZEF        ! latent heat flux (W/m2)
REAL, DIMENSION(SIZE(PTA))        :: ZTAUR      ! momentum flx due to rain (N/m2)
REAL, DIMENSION(SIZE(PTA))        :: ZRF        ! sensible flx due to rain (W/m2)
REAL, DIMENSION(SIZE(PTA))        :: ZEFWEBB    ! Webb corr. on latent flx (W/m2)

REAL, DIMENSION(SIZE(PTA))        :: ZVMOD      ! wind intensity at atm level (m/s)
REAL, DIMENSION(SIZE(PTA))        :: ZQSATA     ! sat.spec.hum. at atm level (kg/kg)
REAL, DIMENSION(SIZE(PTA))        :: ZLVA       ! vap.heat of pure water at atm level (J/kg)
REAL, DIMENSION(SIZE(PTA))        :: ZLVS       ! vap.heat of seawater at sea surface (J/kg)
REAL, DIMENSION(SIZE(PTA))        :: ZCPA       ! specif.heat moist air (J/kg/K)
REAL, DIMENSION(SIZE(PTA))        :: ZVISA      ! kinemat.visc. of dry air (m2/s)
REAL, DIMENSION(SIZE(PTA))        :: ZDU        ! U   vert.grad. (real atm)
REAL, DIMENSION(SIZE(PTA))        :: ZDT,ZDQ    ! T,Q vert.grad. (real atm)
REAL, DIMENSION(SIZE(PTA))        :: ZDDU       ! U   vert.grad. (real atm + gust)
REAL, DIMENSION(SIZE(PTA))        :: ZDDT,ZDDQ  ! T,Q vert.grad. (real atm + WL/CS)
REAL, DIMENSION(SIZE(PTA))        :: ZUSR       ! velocity scaling param. (m/s)
                                                ! =friction velocity
REAL, DIMENSION(SIZE(PTA))        :: ZTSR       ! temperature scaling param. (K)
REAL, DIMENSION(SIZE(PTA))        :: ZQSR       ! humidity scaling param. (kg/kg)
REAL, DIMENSION(SIZE(PTA))        :: ZDELTAU10N,ZDELTAT10N,ZDELTAQ10N
                                                ! U,T,Q vert.grad. (10m, neutral atm)
REAL, DIMENSION(SIZE(PTA))        :: ZUSR0,ZTSR0,ZQSR0    ! ITERATIVE PROCESS
REAL, DIMENSION(SIZE(PTA))        :: ZDUSTO,ZDTSTO,ZDQSTO ! ITERATIVE PROCESS
REAL, DIMENSION(SIZE(PTA))        :: ZPSIU,ZPSIT! PSI funct for U, T/Q (Z0 comp)
REAL, DIMENSION(SIZE(PTA))        :: ZCHARN     ! Charnock parameter   (Z0 comp)

REAL, DIMENSION(SIZE(PTA))        :: ZUSTAR2    ! square of friction velocity
REAL, DIMENSION(SIZE(PTA))        :: ZAC        ! aerodynamical conductance
REAL, DIMENSION(SIZE(PTA))        :: ZDIRCOSZW  ! orography slope cosine
                                                ! (=1 on water!)
REAL, DIMENSION(SIZE(PTA))        :: ZPARUN,ZPARTN,ZPARQN ! neutral parameter for U,T,Q
REAL, DIMENSION(0:5)              :: ZCOEFU,ZCOEFT,ZCOEFQ

! local constants
LOGICAL :: OPCVFLX              ! to force convergence
INTEGER :: NITERMAX             ! nb of iterations to get free convergence
INTEGER :: NITERSUP             ! nb of additional iterations if OPCVFLX=.TRUE.
INTEGER :: NITERFL              ! maximum number of iterations
REAL    :: ZETV,ZRDSRV          ! thermodynamic constants
REAL    :: ZSQR3
REAL    :: ZLMOMIN,ZLMOMAX      ! min/max value of Obukhovs stability param. z/l
REAL    :: ZBTA,ZGMA            ! parameters of the stability functions
REAL    :: ZDUSR0,ZDTSR0,ZDQSR0 ! maximum gap for USR/TSR/QSR between 2 steps
REAL    :: ZP00                 ! [OPRECIP] - water vap. diffusiv.ref.press.(Pa)
REAL    :: ZUTU,ZUTT,ZUTQ       ! U10n threshold in ECUME parameterisation
REAL    :: ZCDIRU,ZCDIRT,ZCDIRQ ! coef directeur pour fonction affine U,T,Q
REAL    :: ZORDOU,ZORDOT,ZORDOQ ! ordonnee a l'origine pour fonction affine U,T,Q

INTEGER :: JJ                                   ! for ITERATIVE PROCESS
INTEGER :: JLON,JK
REAL    :: ZLMOU,ZLMOT                          ! Obukhovs param. z/l for U, T/Q
REAL    :: ZPSI_U,ZPSI_T                        ! PSI funct. for U, T/Q
REAL    :: Z0TSEA,Z0QSEA                        ! roughness length for T, Q
REAL    :: ZCHIC,ZCHIK,ZPSIC,ZPSIK,ZLOGUS10,ZLOGTS10
REAL    :: ZTAC,ZCPWA,ZDQSDT,ZDWAT,ZDTMP,ZBULB  ! [OPRECIP]
REAL    :: ZWW                                  ! [OPWEBB]

REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
!-------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('ECUMEV6_FLUX',0,zhook_handle)
!
zdusr0   = 1.e-06
zdtsr0   = 1.e-06
zdqsr0   = 1.e-09
!
nitermax = 5
nitersup = 5
opcvflx  = .true.
!
niterfl = nitermax
IF (opcvflx) niterfl = nitermax+nitersup
!
zcoefu = (/ 1.00e-03, 3.66e-02, -1.92e-03, 2.32e-04, -7.02e-06,  6.40e-08 /)
zcoeft = (/ 5.36e-03, 2.90e-02, -1.24e-03, 4.50e-04, -2.06e-05,       0.0 /)
zcoefq = (/ 1.00e-03, 3.59e-02, -2.87e-04,      0.0,       0.0,       0.0 /)
!
zutu = 40.0
zutt = 14.4
zutq = 10.0
!
zcdiru = zcoefu(1) + 2.0*zcoefu(2)*zutu + 3.0*zcoefu(3)*zutu**2   &
                   + 4.0*zcoefu(4)*zutu**3 + 5.0*zcoefu(5)*zutu**4
zcdirt = zcoeft(1) + 2.0*zcoeft(2)*zutt + 3.0*zcoeft(3)*zutt**2   &
                   + 4.0*zcoeft(4)*zutt**3
zcdirq = zcoefq(1) + 2.0*zcoefq(2)*zutq
!
zordou = zcoefu(0) + zcoefu(1)*zutu + zcoefu(2)*zutu**2 + zcoefu(3)*zutu**3   &
                   + zcoefu(4)*zutu**4 + zcoefu(5)*zutu**5
zordot = zcoeft(0) + zcoeft(1)*zutt + zcoeft(2)*zutt**2 + zcoeft(3)*zutt**3   &
                   + zcoeft(4)*zutt**4
zordoq = zcoefq(0) + zcoefq(1)*zutq + zcoefq(2)*zutq**2
!
!-------------------------------------------------------------------------------
!
!       1.   AUXILIARY CONSTANTS & ARRAY INITIALISATION BY UNDEFINED VALUES.
!       --------------------------------------------------------------------
!
zdircoszw(:) = 1.0
!
zetv    = xrv/xrd-1.0   !~0.61 (cf Liu et al. 1979)
zrdsrv  = xrd/xrv       !~0.622
zsqr3   = sqrt(3.0)
zlmomin = -200.0
zlmomax = 0.25
zbta    = 16.0
zgma    = 7.0           !initially =4.7, modified to 7.0 following G. Caniaux
!
zp00    = 1013.25e+02
!
pcd(:) = xundef
pch(:) = xundef
pce(:) = xundef
pcdn(:) = xundef
zusr(:) = xundef
ztsr(:) = xundef
zqsr(:) = xundef
ztau(:) = xundef
zhf(:) = xundef
zef(:) = xundef
!
psfth(:) = xundef
psftq(:) = xundef
pustar(:) = xundef
presa(:) = xundef
pri(:) = xundef
!
ztaur(:) = 0.0
zrf(:) = 0.0
zefwebb(:) = 0.0
!
!-------------------------------------------------------------------------------
!
!       2.   INITIALISATIONS BEFORE ITERATIVE LOOP.
!       -------------------------------------------
!
zvmod(:) = wind_threshold(pvmod(:),puref(:))    !set a minimum value to wind
!
!       2.0. Radiative fluxes - For warm layer & cool skin
!
!       2.0b. Warm Layer correction
!
!       2.1. Specific humidity at saturation
!
WHERE(psss(:)>0.0.AND.psss(:)/=xundef)
  pqsat(:) = qsat_seawater2(psst(:),pps(:),psss(:))    !at sea surface
ELSEWHERE
  pqsat(:) = qsat_seawater(psst(:),pps(:))            !at sea surface
ENDWHERE
zqsata(:) = qsat(pta(:),ppa(:))                         !at atm level
!
!       2.2. Gradients at the air-sea interface
!
zdu(:) = zvmod(:)               !one assumes u is measured / sea surface current
zdt(:) = pta(:)/pexna(:)-psst(:)/pexns(:)
zdq(:) = pqa(:)-pqsat(:)
!
!       2.3. Latent heat of vaporisation
!
zlva(:) = xlvtt+(xcpv-xcl)*(pta(:)-xtt)                !of pure water at atm level
zlvs(:) = xlvtt+(xcpv-xcl)*(psst(:)-xtt)                !of pure water at sea surface
WHERE(psss(:)>0.0.AND.psss(:)/=xundef)
  zlvs(:) = zlvs(:)*(1.0-1.00472e-3*psss(:))            !of seawater at sea surface
ENDWHERE
!
!       2.4. Specific heat of moist air (Businger 1982)
!
!ZCPA(:) = XCPD*(1.0+(XCPV/XCPD-1.0)*PQA(:))
zcpa(:) = xcpd
!
!       2.4b Kinematic viscosity of dry air (Andreas 1989, CRREL Rep. 89-11)
!
zvisa(:) = 1.326e-05*(1.0+6.542e-03*(pta(:)-xtt)+8.301e-06*(pta(:)-xtt)**2   &
           -4.84e-09*(pta(:)-xtt)**3)
!
!       2.4c Coefficients for warm layer and/or cool skin correction
!
!       2.5. Initial guess
!
zddu(:) = zdu(:)
zddt(:) = zdt(:)
zddq(:) = zdq(:)
zddu(:) = sign(max(abs(zddu(:)),10.0*zdusr0),zddu(:))
zddt(:) = sign(max(abs(zddt(:)),10.0*zdtsr0),zddt(:))
zddq(:) = sign(max(abs(zddq(:)),10.0*zdqsr0),zddq(:))
!
jcv(:) = -1
zusr(:) = 0.04*zddu(:)
ztsr(:) = 0.04*zddt(:)
zqsr(:) = 0.04*zddq(:)
zdeltau10n(:) = zddu(:)
zdeltat10n(:) = zddt(:)
zdeltaq10n(:) = zddq(:)
jiter(:) = 99
!
! In the following, we suppose that Richardson number PRI < XRIMAX
! If not true, Monin-Obukhov theory can't (and therefore shouldn't) be applied !
!-------------------------------------------------------------------------------
!
!       3.   ITERATIVE LOOP TO COMPUTE U*, T*, Q*.
!       ------------------------------------------
!
DO jj=1,niterfl
  DO jlon=1,SIZE(pta)
!
  IF (jcv(jlon) == -1) THEN
    zusr0(jlon)=zusr(jlon)
    ztsr0(jlon)=ztsr(jlon)
    zqsr0(jlon)=zqsr(jlon)
    IF (jj == nitermax+1 .OR. jj == nitermax+nitersup) THEN
      zdeltau10n(jlon) = 0.5*(zdusto(jlon)+zdeltau10n(jlon))    !forced convergence
      zdeltat10n(jlon) = 0.5*(zdtsto(jlon)+zdeltat10n(jlon))
      zdeltaq10n(jlon) = 0.5*(zdqsto(jlon)+zdeltaq10n(jlon))
      IF (jj == nitermax+nitersup) jcv(jlon)=3
    ENDIF
    zdusto(jlon) = zdeltau10n(jlon)
    zdtsto(jlon) = zdeltat10n(jlon)
    zdqsto(jlon) = zdeltaq10n(jlon)
!
!       3.1. Neutral parameter for wind speed (ECUME_V6 formulation)
!
    IF (zdeltau10n(jlon) <= zutu) THEN
      zparun(jlon) = zcoefu(0) + zcoefu(1)*zdeltau10n(jlon)      &
                               + zcoefu(2)*zdeltau10n(jlon)**2   &
                               + zcoefu(3)*zdeltau10n(jlon)**3   &
                               + zcoefu(4)*zdeltau10n(jlon)**4   &
                               + zcoefu(5)*zdeltau10n(jlon)**5
    ELSE
      zparun(jlon) = zcdiru*(zdeltau10n(jlon)-zutu) + zordou
    ENDIF
    pcdn(jlon) = (zparun(jlon)/zdeltau10n(jlon))**2
!
!       3.2. Neutral parameter for temperature (ECUME_V6 formulation)
!
    IF (zdeltau10n(jlon) <= zutt) THEN
      zpartn(jlon) = zcoeft(0) + zcoeft(1)*zdeltau10n(jlon)      &
                               + zcoeft(2)*zdeltau10n(jlon)**2   &
                               + zcoeft(3)*zdeltau10n(jlon)**3   &
                               + zcoeft(4)*zdeltau10n(jlon)**4
    ELSE
      zpartn(jlon) = zcdirt*(zdeltau10n(jlon)-zutt) + zordot
    ENDIF
!
!       3.3. Neutral parameter for humidity (ECUME_V6 formulation)
!
    IF (zdeltau10n(jlon) <= zutq) THEN
      zparqn(jlon) = zcoefq(0) + zcoefq(1)*zdeltau10n(jlon)      &
                               + zcoefq(2)*zdeltau10n(jlon)**2
    ELSE
      zparqn(jlon) = zcdirq*(zdeltau10n(jlon)-zutq) + zordoq
    ENDIF
!
!       3.4. Scaling parameters U*, T*, Q*
!
    zusr(jlon) = zparun(jlon)
    ztsr(jlon) = zpartn(jlon)*zdeltat10n(jlon)/zdeltau10n(jlon)
    zqsr(jlon) = zparqn(jlon)*zdeltaq10n(jlon)/zdeltau10n(jlon)
!
!       3.4b Gustiness factor (Deardorff 1970)
!
!       3.4c Cool skin correction
!
!       3.5. Obukhovs stability param. z/l following Liu et al. (JAS, 1979)
!
! For U
    zlmou = puref(jlon)*xg*xkarman*(ztsr(jlon)/pta(jlon)   &
            +zetv*zqsr(jlon)/(1.0+zetv*pqa(jlon)))/zusr(jlon)**2
! For T/Q
    zlmot = zlmou*(pzref(jlon)/puref(jlon))
    zlmou = max(min(zlmou,zlmomax),zlmomin)
    zlmot = max(min(zlmot,zlmomax),zlmomin)
!
!       3.6. Stability function psi (see Liu et al, 1979 ; Dyer and Hicks, 1970)
!            Modified to include convective form following Fairall (unpublished)
!
! For U
    IF (zlmou == 0.0) THEN
      zpsi_u = 0.0
    ELSEIF (zlmou > 0.0) THEN
      zpsi_u = -zgma*zlmou
    ELSE
      zchik  = (1.0-zbta*zlmou)**0.25
      zpsik  = 2.0*log((1.0+zchik)/2.0)  &
                +log((1.0+zchik**2)/2.0) &
                -2.0*atan(zchik)+0.5*xpi
      zchic  = (1.0-12.87*zlmou)**(1.0/3.0)       !for very unstable conditions
      zpsic  = 1.5*log((zchic**2+zchic+1.0)/3.0)   &
                -zsqr3*atan((2.0*zchic+1.0)/zsqr3) &
                +xpi/zsqr3
      zpsi_u = zpsic+(zpsik-zpsic)/(1.0+zlmou**2) !match Kansas & free-conv. forms
    ENDIF
    zpsiu(jlon) = zpsi_u
! For T/Q
    IF (zlmot == 0.0) THEN
      zpsi_t = 0.0
    ELSEIF (zlmot > 0.0) THEN
      zpsi_t = -zgma*zlmot
    ELSE
      zchik  = (1.0-zbta*zlmot)**0.25
      zpsik  = 2.0*log((1.0+zchik**2)/2.0)
      zchic  = (1.0-12.87*zlmot)**(1.0/3.0)       !for very unstable conditions
      zpsic  = 1.5*log((zchic**2+zchic+1.0)/3.0)   &
                -zsqr3*atan((2.0*zchic+1.0)/zsqr3) &
                +xpi/zsqr3
      zpsi_t = zpsic+(zpsik-zpsic)/(1.0+zlmot**2) !match Kansas & free-conv. forms
    ENDIF
    zpsit(jlon) = zpsi_t
!
!       3.7. Update air-sea gradients
!
    zddu(jlon) = zdu(jlon)
    zddt(jlon) = zdt(jlon)
    zddq(jlon) = zdq(jlon)
    zddu(jlon) = sign(max(abs(zddu(jlon)),10.0*zdusr0),zddu(jlon))
    zddt(jlon) = sign(max(abs(zddt(jlon)),10.0*zdtsr0),zddt(jlon))
    zddq(jlon) = sign(max(abs(zddq(jlon)),10.0*zdqsr0),zddq(jlon))
    zlogus10   = log(puref(jlon)/10.0)
    zlogts10   = log(pzref(jlon)/10.0)
    zdeltau10n(jlon) = zddu(jlon)-zusr(jlon)*(zlogus10-zpsi_u)/xkarman
    zdeltat10n(jlon) = zddt(jlon)-ztsr(jlon)*(zlogts10-zpsi_t)/xkarman
    zdeltaq10n(jlon) = zddq(jlon)-zqsr(jlon)*(zlogts10-zpsi_t)/xkarman
    zdeltau10n(jlon) = sign(max(abs(zdeltau10n(jlon)),10.0*zdusr0),   &
                            zdeltau10n(jlon))
    zdeltat10n(jlon) = sign(max(abs(zdeltat10n(jlon)),10.0*zdtsr0),   &
                            zdeltat10n(jlon))
    zdeltaq10n(jlon) = sign(max(abs(zdeltaq10n(jlon)),10.0*zdqsr0),   &
                            zdeltaq10n(jlon))
!
!       3.8. Test convergence for U*, T*, Q*
!
    IF (abs(zusr(jlon)-zusr0(jlon)) < zdusr0 .AND.   &
        abs(ztsr(jlon)-ztsr0(jlon)) < zdtsr0 .AND.   &
        abs(zqsr(jlon)-zqsr0(jlon)) < zdqsr0) THEN
      jcv(jlon) = 1                                     !free convergence
      IF (jj >= nitermax+1) jcv(jlon) = 2               !leaded convergence
    ENDIF
    jiter(jlon) = jj
  ENDIF
!
  ENDDO
ENDDO
!
!-------------------------------------------------------------------------------
!
!       4.   COMPUTATION OF TURBULENT FLUXES AND EXCHANGE COEFFICIENTS.
!       ---------------------------------------------------------------
!
DO jlon=1,SIZE(pta)
!
!       4.1. Surface turbulent fluxes
!            (ATM CONV.: ZTAU<<0 ; ZHF,ZEF<0 if atm looses heat)
!
  ztau(jlon) = -prhoa(jlon)*zusr(jlon)**2
  zhf(jlon) = -prhoa(jlon)*zcpa(jlon)*zusr(jlon)*ztsr(jlon)
  zef(jlon) = -prhoa(jlon)*zlvs(jlon)*zusr(jlon)*zqsr(jlon)
!
!       4.2. Exchange coefficients PCD, PCH, PCE
!
  pcd(jlon) = (zusr(jlon)/zddu(jlon))**2
  pch(jlon) = (zusr(jlon)*ztsr(jlon))/(zddu(jlon)*zddt(jlon))
  pce(jlon) = (zusr(jlon)*zqsr(jlon))/(zddu(jlon)*zddq(jlon))
!
!       4.3. Stochastic perturbation of turbulent fluxes
!
  IF( opertflux )THEN
    ztau(jlon) = ztau(jlon)* ( 1. + ppertflux(jlon) / 2. )
    zhf(jlon) = zhf(jlon)*  ( 1. + ppertflux(jlon) / 2. )
    zef(jlon) = zef(jlon)*  ( 1. + ppertflux(jlon) / 2. )
  ENDIF
!
ENDDO
!
!-------------------------------------------------------------------------------
!
!       5.   COMPUTATION OF FLUX CORRECTIONS DUE TO RAINFALL.
!            (ATM conv: ZRF<0 if atm. looses heat, ZTAUR<<0)
!       -----------------------------------------------------
!
IF (oprecip) THEN
  DO jlon=1,SIZE(pta)
!
!       5.1. Momentum flux due to rainfall (ZTAUR, N/m2)
!
! See pp3752 in FBR96.
    ztaur(jlon) = -0.85*prain(jlon)*pvmod(jlon)
!
!       5.2. Sensible heat flux due to rainfall (ZRF, W/m2)
!
! See Eq.12 in GoF95 with ZCPWA as specific heat of water at atm level (J/kg/K),
! ZDQSDT from Clausius-Clapeyron relation, ZDWAT as water vapor diffusivity 
! (Eq.13-3 of Pruppacher and Klett, 1978), ZDTMP as heat diffusivity, and ZBULB
! as wet-bulb factor (Eq.11 in GoF95).
!
    ztac   = pta(jlon)-xtt
    zcpwa  = 4217.51 -3.65566*ztac +0.1381*ztac**2       &
              -2.8309e-03*ztac**3 +3.42061e-05*ztac**4   &
              -2.18107e-07*ztac**5 +5.74535e-10*ztac**6
    zdqsdt = (zlva(jlon)*zqsata(jlon))/(xrv*pta(jlon)**2)
    zdwat  = 2.11e-05*(zp00/ppa(jlon))*(pta(jlon)/xtt)**1.94
    zdtmp  = (1.0+3.309e-03*ztac-1.44e-06*ztac**2)   &
              *0.02411/(prhoa(jlon)*zcpa(jlon))
    zbulb  = 1.0/(1.0+zdqsdt*(zlva(jlon)*zdwat)/(zcpa(jlon)*zdtmp))
    zrf(jlon) = prain(jlon)*zcpwa*zbulb*((psst(jlon)-pta(jlon))   &
                +(pqsat(jlon)-pqa(jlon))*(zlva(jlon)*zdwat)/(zcpa(jlon)*zdtmp))
!
  ENDDO
ENDIF
!
!-------------------------------------------------------------------------------
!
!       6.   COMPUTATION OF WEBB CORRECTION TO LATENT HEAT FLUX (ZEFWEBB, W/m2).
!       ------------------------------------------------------------------------
!
! See Eq.21 and Eq.22 in FBR96.
IF (opwebb) THEN
  DO jlon=1,SIZE(pta)
    zww = (1.0+zetv)*(zusr(jlon)*zqsr(jlon))   &
           +(1.0+(1.0+zetv)*pqa(jlon))*(zusr(jlon)*ztsr(jlon))/pta(jlon)
    zefwebb(jlon) = -prhoa(jlon)*zlvs(jlon)*zww*pqa(jlon)
  ENDDO
ENDIF
!
!-------------------------------------------------------------------------------
!
!       7.   FINAL STEP : TOTAL SURFACE FLUXES AND DERIVED DIAGNOSTICS. 
!       ---------------------------------------------------------------
!
!       7.1. Richardson number
!
 CALL surface_ri(psst,pqsat,pexns,pexna,pta,pqa,   &
                pzref,puref,zdircoszw,pvmod,pri)
!
!       7.2. Friction velocity which contains correction due to rain
!
zustar2(:) = -(ztau(:)+ztaur(:))/prhoa(:)
!
IF (oprecip) THEN
  pcd(:) = zustar2(:)/zddu(:)**2
ENDIF
!
pustar(:) = sqrt(zustar2(:))
!
!       7.3. Aerodynamical conductance and resistance
!
zac(:) = pch(:)*zddu(:)
presa(:) = 1.0/zac(:)
!
!       7.4. Total surface fluxes
!
psfth(:) =  zhf(:)+zrf(:)
psftq(:) = (zef(:)+zefwebb(:))/zlvs(:)
!
!       7.5. Charnock number
!
IF (ccharnock == 'OLD') THEN
  zcharn(:) = xvchrnk
ELSE            !modified for moderate wind speed as in COARE3.0
  zcharn(:) = min(0.018,max(0.011,0.011+(0.007/8.0)*(zddu(:)-10.0)))
ENDIF
!
!       7.6. Roughness lengths Z0 and Z0H over sea
!
IF (kz0 == 0) THEN      ! ARPEGE formulation
  pz0sea(:) = (zcharn(:)/xg)*zustar2(:) + xvz0cm*pcd(:)/pcdn(:)
  pz0hsea(:) = pz0sea(:)
ELSEIF (kz0 == 1) THEN  ! Smith (1988) formulation
  pz0sea(:) = (zcharn(:)/xg)*zustar2(:) + 0.11*zvisa(:)/pustar(:)
  pz0hsea(:) = pz0sea(:)
ELSEIF (kz0 == 2) THEN  ! Direct computation using the stability functions
  DO jlon=1,SIZE(pta)
    pz0sea(jlon) = puref(jlon)/exp(xkarman*zddu(jlon)/pustar(jlon)+zpsiu(jlon))
    z0tsea        = pzref(jlon)/exp(xkarman*zddt(jlon)/ztsr(jlon)+zpsit(jlon))
    z0qsea        = pzref(jlon)/exp(xkarman*zddq(jlon)/zqsr(jlon)+zpsit(jlon))
    pz0hsea(jlon) = 0.5*(z0tsea+z0qsea)
  ENDDO
ENDIF
!
!
IF (lhook) CALL dr_hook('ECUMEV6_FLUX',1,zhook_handle)
!
!-------------------------------------------------------------------------------
   END SUBROUTINE ecumev6_flux