coare30_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.
!     #########
    SUBROUTINE coare30_flux (S,PZ0SEA,PTA,PEXNA,PRHOA,PSST,PEXNS,PQA,  &
            PVMOD,PZREF,PUREF,PPS,PQSAT,PSFTH,PSFTQ,PUSTAR,PCD,PCDN,PCH,PCE,PRI,&
            PRESA,PRAIN,PZ0HSEA)  
!     #######################################################################
!
!
!!****  *COARE25_FLUX*  
!!
!!    PURPOSE
!!    -------
!      Calculate the surface fluxes of heat, moisture, and momentum over
!      sea surface with bulk algorithm COARE3.0. 
!     
!!**  METHOD
!!    ------
!      transfer coefficients were obtained using a dataset which combined COARE
!      data with those from three other ETL field experiments, and reanalysis of
!      the HEXMAX data (DeCosmos et al. 1996). 
!      ITERMAX=3 
!      Take account of the surface gravity waves on the velocity roughness and 
!      hence the momentum transfer coefficient
!        NGRVWAVES=0 no gravity waves action (Charnock) !default value
!        NGRVWAVES=1 wave age parameterization of Oost et al. 2002
!        NGRVWAVES=2 model of Taylor and Yelland 2001
!
!!    EXTERNAL
!!    --------
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------ 
!!      
!!    REFERENCE
!!    ---------
!!      Fairall et al (2003), J. of Climate, vol. 16, 571-591
!!      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 C. Fairall's code)
!!
!!    MODIFICATIONS
!!    -------------
!!      Original     1/06/2006
!!      B. Decharme    06/2009 limitation of Ri
!!      B. Decharme    09/2012 Bug in Ri calculation and limitation of Ri in surface_ri.F90
!!      B. Decharme    06/2013 bug in z0 (output) computation 
!!      J.Escobar      06/2013  for REAL4/8 add EPSILON management
!!      C. Lebeaupin   03/2014 bug if PTA=PSST and PEXNA=PEXNS: set a minimum value
!!                             add abort if no convergence
!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!               ------------
!
!
USE modd_seaflux_n, ONLY : seaflux_t
!
USE modd_csts,       ONLY : xkarman, xg, xstefan, xrd, xrv, xpi, &
                            xlvtt, xcl, xcpd, xcpv, xrholw, xtt, &
                            xp00
USE modd_surf_atm,   ONLY : xvz0cm
!
USE modd_surf_par,   ONLY : xundef, xsurf_epsilon
USE modd_water_par
!
USE modi_surface_ri
USE modi_wind_threshold
USE mode_coare30_psi
!
USE mode_thermos
!
!
USE modi_abor1_sfx
!
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
!
!
TYPE(seaflux_t), INTENT(INOUT) :: S
!
REAL, DIMENSION(:), INTENT(IN)       :: PTA   ! air temperature at atm. level (K)
REAL, DIMENSION(:), INTENT(IN)       :: PQA   ! air humidity at atm. level (kg/kg)
REAL, DIMENSION(:), INTENT(IN)       :: PEXNA ! Exner function at atm. level
REAL, DIMENSION(:), INTENT(IN)       :: PRHOA ! air density at atm. level
REAL, DIMENSION(:), INTENT(IN)       :: PVMOD ! module of wind at atm. wind level (m/s)
REAL, DIMENSION(:), INTENT(IN)       :: PZREF ! atm. level for temp. and humidity (m)
REAL, DIMENSION(:), INTENT(IN)       :: PUREF ! atm. level for wind (m)
REAL, DIMENSION(:), INTENT(IN)       :: PSST  ! Sea Surface Temperature (K)
REAL, DIMENSION(:), INTENT(IN)       :: PEXNS ! Exner function at sea surface
REAL, DIMENSION(:), INTENT(IN)       :: PPS   ! air pressure at sea surface (Pa)
REAL, DIMENSION(:), INTENT(IN)       :: PRAIN !precipitation rate (kg/s/m2)
!
REAL, DIMENSION(:), INTENT(INOUT)    :: PZ0SEA! roughness length over the ocean
!                                                                                 
!  surface fluxes : latent heat, sensible heat, friction fluxes
REAL, DIMENSION(:), INTENT(OUT)      :: PSFTH ! heat flux (W/m2)
REAL, DIMENSION(:), INTENT(OUT)      :: PSFTQ ! water flux (kg/m2/s)
REAL, DIMENSION(:), INTENT(OUT)      :: PUSTAR! friction velocity (m/s)
!
! diagnostics
REAL, DIMENSION(:), INTENT(OUT)      :: PQSAT ! humidity at saturation
REAL, DIMENSION(:), INTENT(OUT)      :: PCD   ! heat drag coefficient
REAL, DIMENSION(:), INTENT(OUT)      :: PCDN  ! momentum drag coefficient
REAL, DIMENSION(:), INTENT(OUT)      :: PCH   ! neutral momentum drag coefficient
REAL, DIMENSION(:), INTENT(OUT)      :: PCE  !transfer coef. for latent heat flux
REAL, DIMENSION(:), INTENT(OUT)      :: PRI   ! Richardson number
REAL, DIMENSION(:), INTENT(OUT)      :: PRESA ! aerodynamical resistance
REAL, DIMENSION(:), INTENT(OUT)      :: PZ0HSEA ! heat roughness length
!
!
!*      0.2    declarations of local variables
!
REAL, DIMENSION(SIZE(PTA))      :: ZVMOD    ! wind intensity
REAL, DIMENSION(SIZE(PTA))      :: ZPA      ! Pressure at atm. level
REAL, DIMENSION(SIZE(PTA))      :: ZTA      ! Temperature at atm. level
REAL, DIMENSION(SIZE(PTA))      :: ZQASAT   ! specific humidity at saturation  at atm. level (kg/kg)
!
REAL, DIMENSION(SIZE(PTA))      :: ZO       ! rougness length ref 
REAL, DIMENSION(SIZE(PTA))      :: ZWG      ! gustiness factor (m/s)
!
REAL, DIMENSION(SIZE(PTA))      :: ZDU,ZDT,ZDQ,ZDUWG !differences
!
REAL, DIMENSION(SIZE(PTA))      :: ZUSR        !velocity scaling parameter "ustar" (m/s) = friction velocity
REAL, DIMENSION(SIZE(PTA))      :: ZTSR        !temperature sacling parameter "tstar" (degC)
REAL, DIMENSION(SIZE(PTA))      :: ZQSR        !humidity scaling parameter "qstar" (kg/kg)
!
REAL, DIMENSION(SIZE(PTA))      :: ZU10,ZT10   !vertical profils (10-m height) 
REAL, DIMENSION(SIZE(PTA))      :: ZVISA       !kinematic viscosity of dry air
REAL, DIMENSION(SIZE(PTA))      :: ZO10,ZOT10  !roughness length at 10m
REAL, DIMENSION(SIZE(PTA))      :: ZCD,ZCT,ZCC
REAL, DIMENSION(SIZE(PTA))      :: ZCD10,ZCT10 !transfer coef. at 10m
REAL, DIMENSION(SIZE(PTA))      :: ZRIBU,ZRIBCU
REAL, DIMENSION(SIZE(PTA))      :: ZETU,ZL10
!
REAL, DIMENSION(SIZE(PTA))      :: ZCHARN                      !Charnock number depends on wind module
REAL, DIMENSION(SIZE(PTA))      :: ZTWAVE,ZHWAVE,ZCWAVE,ZLWAVE !to compute gravity waves' impact
!
REAL, DIMENSION(SIZE(PTA))      :: ZZL,ZZTL!,ZZQL    !Obukhovs stability 
                                                     !param. z/l for u,T,q
REAL, DIMENSION(SIZE(PTA))      :: ZRR
REAL, DIMENSION(SIZE(PTA))      :: ZOT,ZOQ           !rougness length ref
REAL, DIMENSION(SIZE(PTA))      :: ZPUZ,ZPTZ,ZPQZ    !PHI funct. for u,T,q 
!
REAL, DIMENSION(SIZE(PTA))      :: ZBF               !constants to compute gustiness factor
!
REAL, DIMENSION(SIZE(PTA))      :: ZTAU       !momentum flux (W/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))      :: ZWBAR      !diag for webb correction but not used here after
REAL, DIMENSION(SIZE(PTA))      :: ZTAUR      !momentum flux due to rain (W/m2)
REAL, DIMENSION(SIZE(PTA))      :: ZRF        !sensible heat flux due to rain (W/m2)
REAL, DIMENSION(SIZE(PTA))      :: ZCHN,ZCEN  !neutral coef. for heat and vapor
!
REAL, DIMENSION(SIZE(PTA))      :: ZLV      !latent heat constant
!
REAL, DIMENSION(SIZE(PTA))      :: ZTAC,ZDQSDT,ZDTMP,ZDWAT,ZALFAC ! for precipitation impact
REAL, DIMENSION(SIZE(PTA))      :: ZXLR                           ! vaporisation  heat  at a given temperature
REAL, DIMENSION(SIZE(PTA))      :: ZCPLW                          ! specific heat for water at a given temperature 
!
REAL, DIMENSION(SIZE(PTA))      :: ZUSTAR2  ! square of friction velocity
!
REAL, DIMENSION(SIZE(PTA))      :: ZDIRCOSZW! orography slope cosine (=1 on water!)
REAL, DIMENSION(SIZE(PTA))      :: ZAC      ! Aerodynamical conductance
!
!
INTEGER, DIMENSION(SIZE(PTA))   :: ITERMAX             ! maximum number of iterations
!
REAL    :: ZRVSRDM1,ZRDSRV,ZR2 ! thermodynamic constants
REAL    :: ZBETAGUST           !gustiness factor
REAL    :: ZZBL                !atm. boundary layer depth (m)
REAL    :: ZVISW               !m2/s kinematic viscosity of water
REAL    :: ZS                  !height of rougness length ref
REAL    :: ZCH10               !transfer coef. at 10m
!
INTEGER :: J, JLOOP    !loop indice
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
!-------------------------------------------------------------------------------
!
!       1.     Initializations
!              ---------------
!
!       1.1   Constants and parameters
!
IF (lhook) CALL dr_hook('COARE30_FLUX',0,zhook_handle)
!
zrvsrdm1  = xrv/xrd-1. ! 0.607766
zrdsrv    = xrd/xrv    ! 0.62198
zr2       = 1.-zrdsrv  ! pas utilisé dans cette routine
zbetagust = 1.2        ! value based on TOGA-COARE experiment
zzbl      = 600.       ! Set a default value for boundary layer depth
zs        = 10.        ! Standard heigth =10m
zch10     = 0.00115
!
zvisw     = 1.e-6
!
!       1.2   Array initialization by undefined values
!
psfth(:)=xundef
psftq(:)=xundef
pustar(:)=xundef
!
pcd(:) = xundef
pcdn(:) = xundef
pch(:) = xundef
pce(:) =xundef
pri(:) = xundef
!
presa(:)=xundef
!
!-------------------------------------------------------------------------------
!       2. INITIAL GUESS FOR THE ITERATIVE METHOD 
!          -------------------------------------
!
!       2.0     Temperature 
!
! Set a non-zero value for the temperature gradient
!
WHERE((pta(:)*pexns(:)/pexna(:)-psst(:))==0.) 
      zta(:)=pta(:)-1e-3
ELSEWHERE
      zta(:)=pta(:)      
ENDWHERE

!       2.1     Wind and humidity 
!
! Sea surface specific humidity 
!
pqsat(:)=qsat_seawater(psst(:),pps(:))         
!              
! Set a minimum value to wind 
!
zvmod(:) = wind_threshold(pvmod(:),puref(:))
!
! Specific humidity at saturation at the atm. level 
!
zpa(:) = xp00* (pexna(:)**(xcpd/xrd))
zqasat(:) = qsat(zta(:),zpa(:)) 
!
!
zo(:)  = 0.0001
zwg(:) = 0.
IF (s%LPWG) zwg(:) = 0.5
!
zcharn(:) = 0.011  
!
DO j=1,SIZE(pta)
  !
  !      2.2       initial guess
  !    
  zdu(j) = zvmod(j)   !wind speed difference with surface current(=0) (m/s)
                      !initial guess for gustiness factor
  zdt(j) = -(zta(j)/pexna(j)) + (psst(j)/pexns(j)) !potential temperature difference
  zdq(j) = pqsat(j)-pqa(j)                         !specific humidity difference
  !
  zduwg(j) = sqrt(zdu(j)**2+zwg(j)**2)     !wind speed difference including gustiness ZWG
  !
  !      2.3   initialization of neutral coefficients
  !
  zu10(j)  = zduwg(j)*log(zs/zo(j))/log(puref(j)/zo(j))
  zusr(j)  = 0.035*zu10(j)
  zvisa(j) = 1.326e-5*(1.+6.542e-3*(zta(j)-xtt)+&
             8.301e-6*(zta(j)-xtt)**2-4.84e-9*(zta(j)-xtt)**3) !Andrea (1989) CRREL Rep. 89-11
  ! 
  zo10(j) = zcharn(j)*zusr(j)*zusr(j)/xg+0.11*zvisa(j)/zusr(j)
  zcd(j)  = (xkarman/log(puref(j)/zo10(j)))**2  !drag coefficient
  zcd10(j)= (xkarman/log(zs/zo10(j)))**2
  zct10(j)= zch10/sqrt(zcd10(j))
  zot10(j)= zs/exp(xkarman/zct10(j))
  !
  !-------------------------------------------------------------------------------
  !             Grachev and Fairall (JAM, 1997)
  zct(j) = xkarman/log(pzref(j)/zot10(j))      !temperature transfer coefficient
  zcc(j) = xkarman*zct(j)/zcd(j)               !z/L vs Rib linear coef.
  !
  zribcu(j) = -puref(j)/(zzbl*0.004*zbetagust**3) !saturation or plateau Rib
  !ZRIBU(J) =-XG*PUREF(J)*(ZDT(J)+ZRVSRDM1*(ZTA(J)-XTT)*ZDQ)/&
  !     &((ZTA(J)-XTT)*ZDUWG(J)**2)
  zribu(j)  = -xg*puref(j)*(zdt(j)+zrvsrdm1*zta(j)*zdq(j))/&
               (zta(j)*zduwg(j)**2)  
  !
  IF (zribu(j)<0.) THEN
    zetu(j) = zcc(j)*zribu(j)/(1.+zribu(j)/zribcu(j))    !Unstable G and F
  ELSE
    zetu(j) = zcc(j)*zribu(j)/(1.+27./9.*zribu(j)/zcc(j))!Stable 
  ENDIF
  !
  zl10(j) = puref(j)/zetu(j) !MO length
  !
ENDDO
!
!  First guess M-O stability dependent scaling params. (u*,T*,q*) to estimate ZO and z/L (ZZL)
zusr(:) = zduwg(:)*xkarman/(log(puref(:)/zo10(:))-psifctu(puref(:)/zl10(:)))
ztsr(:) = -zdt(:)*xkarman/(log(pzref(:)/zot10(:))-psifctt(pzref(:)/zl10(:)))
zqsr(:) = -zdq(:)*xkarman/(log(pzref(:)/zot10(:))-psifctt(pzref(:)/zl10(:)))
!
zzl(:) = 0.0
!
DO j=1,SIZE(pta)
  !
  IF (zetu(j)>50.) THEN
    itermax(j) = 1
  ELSE
    itermax(j) = 3 !number of iterations
  ENDIF
  !
  !then modify Charnork for high wind speeds Chris Fairall's data
  IF (zduwg(j)>10.) zcharn(j) = 0.011 + (0.018-0.011)*(zduwg(j)-10.)/(18.-10.)
  IF (zduwg(j)>18.) zcharn(j) = 0.018
  !
  !                3.  ITERATIVE LOOP TO COMPUTE USR, TSR, QSR 
  !                -------------------------------------------
  !
  zhwave(j) = 0.018*zvmod(j)*zvmod(j)*(1.+0.015*zvmod(j))
  ztwave(j) = 0.729*zvmod(j)
  zcwave(j) = xg*ztwave(j)/(2.*xpi)
  zlwave(j) = ztwave(j)*zcwave(j)
  !
ENDDO
!
   
!
DO jloop=1,maxval(itermax) ! begin of iterative loop
  !
  DO j=1,SIZE(pta)
    !
    IF (jloop.GT.itermax(j)) cycle
    !
    IF (s%NGRVWAVES==0) THEN
      zo(j) = zcharn(j)*zusr(j)*zusr(j)/xg + 0.11*zvisa(j)/zusr(j) !Smith 1988
    ELSE IF (s%NGRVWAVES==1) THEN
      zo(j) = (50./(2.*xpi))*zlwave(j)*(zusr(j)/zcwave(j))**4.5 &
              + 0.11*zvisa(j)/zusr(j)                       !Oost et al. 2002  
    ELSE IF (s%NGRVWAVES==2) THEN
      zo(j) = 1200.*zhwave(j)*(zhwave(j)/zlwave(j))**4.5 &
              + 0.11*zvisa(j)/zusr(j)                       !Taulor and Yelland 2001  
    ENDIF
    !
    zrr(j) = zo(j)*zusr(j)/zvisa(j)
    zoq(j) = min(1.15e-4 , 5.5e-5/zrr(j)**0.6)
    zot(j) = zoq(j)
    !
    zzl(j) = xkarman * xg * puref(j) * &
              ( ztsr(j)*(1.+zrvsrdm1*pqa(j)) + zrvsrdm1*zta(j)*zqsr(j) ) / &
              ( zta(j)*zusr(j)*zusr(j)*(1.+zrvsrdm1*pqa(j)) )  
    zztl(j)= zzl(j)*pzref(j)/puref(j)  ! for T 
!    ZZQL(J)=ZZL(J)*PZREF(J)/PUREF(J)  ! for Q
  ENDDO
  !
  zpuz(:) = psifctu(zzl(:))     
  zptz(:) = psifctt(zztl(:))
  !
  DO j=1,SIZE(pta)
    !
    ! ZPQZ(J)=PSIFCTT(ZZQL(J))    
    zpqz(j) = zptz(j)
    !
    !             3.1 scale parameters
    !
    zusr(j) = zduwg(j)*xkarman/(log(puref(j)/zo(j)) -zpuz(j))
    ztsr(j) = -zdt(j)  *xkarman/(log(pzref(j)/zot(j))-zptz(j))
    zqsr(j) = -zdq(j)  *xkarman/(log(pzref(j)/zoq(j))-zpqz(j))
    !
    !             3.2 Gustiness factor (ZWG)
    !
    IF(s%LPWG) THEN
      zbf(j) = -xg/zta(j)*zusr(j)*(ztsr(j)+zrvsrdm1*zta(j)*zqsr(j))
      IF (zbf(j)>0.) THEN
        zwg(j) = zbetagust*(zbf(j)*zzbl)**(1./3.)
      ELSE
        zwg(j) = 0.2
      ENDIF
    ENDIF  
    zduwg(j) = sqrt(zvmod(j)**2 + zwg(j)**2)
    !
  ENDDO
  !
ENDDO
!-------------------------------------------------------------------------------
!
!            4.  COMPUTE transfer coefficients PCD, PCH, ZCE and SURFACE FLUXES
!                --------------------------------------------------------------
!
ztau(:) = xundef
zhf(:)  = xundef
zef(:)  = xundef
!
zwbar(:) = 0.
ztaur(:) = 0.
zrf(:)   = 0.
!
DO j=1,SIZE(pta)
  !
  !
  !            4. transfert coefficients PCD, PCH and PCE 
  !                 and neutral PCDN, ZCHN, ZCEN 
  !
  pcd(j) = (zusr(j)/zduwg(j))**2.
  pch(j) = zusr(j)*ztsr(j)/(zduwg(j)*(zta(j)*pexns(j)/pexna(j)-psst(j)))
  pce(j) = zusr(j)*zqsr(j)/(zduwg(j)*(pqa(j)-pqsat(j)))
  !
  pcdn(j) = (xkarman/log(zs/zo(j)))**2.
  zchn(j) = (xkarman/log(zs/zo(j)))*(xkarman/log(zs/zot(j)))
  zcen(j) = (xkarman/log(zs/zo(j)))*(xkarman/log(zs/zoq(j)))
  !
  zlv(j) = xlvtt + (xcpv-xcl)*(psst(j)-xtt)
  !
  !            4. 2 surface fluxes 
  !
  IF (abs(pcdn(j))>1.e-2) THEN   !!!! secure COARE3.0 CODE 
    write(*,*) 'pb PCDN in COARE30: ',pcdn(j)
    write(*,*) 'point: ',j,"/",SIZE(pta)
    write(*,*) 'roughness: ', zo(j)
    write(*,*) 'ustar: ',zusr(j)
    write(*,*) 'wind: ',zduwg(j)
    CALL abor1_sfx('COARE30: PCDN too large -> no convergence')
  ELSE
    ztsr(j) = -ztsr(j)
    zqsr(j) = -zqsr(j)
    ztau(j) = -prhoa(j)*zusr(j)*zusr(j)*zvmod(j)/zduwg(j)
    zhf(j)  =  prhoa(j)*xcpd*zusr(j)*ztsr(j)
    zef(j)  =  prhoa(j)*zlv(j)*zusr(j)*zqsr(j)
    !    
    !           4.3 Contributions to surface  fluxes due to rainfall
    !
    ! SB: a priori, le facteur ZRDSRV=XRD/XRV est introduit pour
    !     adapter la formule de Clausius-Clapeyron (pour l'air
    !     sec) au cas humide.
    IF (s%LPRECIP) THEN
      ! 
      ! heat surface  fluxes
      !
      ztac(j)  = zta(j)-xtt
      !
      zxlr(j)  = xlvtt + (xcpv-xcl)* ztac(j)                            ! latent heat of rain vaporization
      zdqsdt(j)= zqasat(j) * zxlr(j) / (xrd*zta(j)**2)                  ! Clausius-Clapeyron relation
      zdtmp(j) = (1.0 + 3.309e-3*ztac(j) -1.44e-6*ztac(j)*ztac(j)) * &  !heat diffusivity
                  0.02411 / (prhoa(j)*xcpd)
      !
      zdwat(j) = 2.11e-5 * (xp00/zpa(j)) * (zta(j)/xtt)**1.94           ! water vapour diffusivity from eq (13.3)
      !                                                                 ! of Pruppacher and Klett (1978)      
      zalfac(j)= 1.0 / (1.0 + &                                         ! Eq.11 in GoF95
                   zrdsrv*zdqsdt(j)*zxlr(j)*zdwat(j)/(zdtmp(j)*xcpd))   ! ZALFAC=wet-bulb factor (sans dim)     
      zcplw(j) = 4224.8482 + ztac(j) * &
                              ( -4.707 + ztac(j) * &
                                (0.08499 + ztac(j) * &
                                  (1.2826e-3 + ztac(j) * &
                                    (4.7884e-5 - 2.0027e-6* ztac(j))))) ! specific heat  
      !       
      zrf(j)   = prain(j) * zcplw(j) * zalfac(j) * &                    !Eq.12 in GoF95 !SIGNE?
                   (psst(j) - zta(j) + (pqsat(j)-pqa(j))*zxlr(j)/xcpd )
      !
      ! Momentum flux due to rainfall  
      !
      ztaur(j)=-0.85*(prain(j) *zvmod(j)) !pp3752 in FBR96
      !
    ENDIF
    !
    !             4.4   Webb correction to latent heat flux
    ! 
    zwbar(j)=- (1./zrdsrv)*zusr(j)*zqsr(j) / (1.0+(1./zrdsrv)*pqa(j)) &
               - zusr(j)*ztsr(j)/zta(j)                        ! Eq.21*rhoa in FBR96    
    !
    !             4.5   friction velocity which contains correction du to rain            
    !
    zustar2(j)= - (ztau(j) + ztaur(j)) / prhoa(j)
    pustar(j) =  sqrt(zustar2(j))
    !
    !             4.6   Total surface fluxes
    !           
    psfth(j) =  zhf(j) + zrf(j)
    psftq(j) =  zef(j) / zlv(j)
    ! 
  ENDIF
ENDDO                      
!-------------------------------------------------------------------------------
!
!       5.  FINAL STEP : TOTAL SURFACE FLUXES AND DERIVED DIAGNOSTICS 
!           -----------
!       5.1    Richardson number
!             
!
zdircoszw(:) = 1.
 CALL surface_ri(psst,pqsat,pexns,pexna,zta,zqasat,&
                pzref,puref,zdircoszw,pvmod,pri   )  
!
!       5.2     Aerodynamical conductance and resistance
!             
zac(:) = pch(:)*zvmod(:)
presa(:) = 1. / max(zac(:),xsurf_epsilon)
!
!       5.3 Z0 and Z0H over sea
!
pz0sea(:) =  zcharn(:) * zustar2(:) / xg + xvz0cm * pcd(:) / pcdn(:)
!
pz0hsea(:) = pz0sea(:)
!
IF (lhook) CALL dr_hook('COARE30_FLUX',1,zhook_handle)
!
!-------------------------------------------------------------------------------
!
END SUBROUTINE coare30_flux