SURFEX v7.3
General documentation of Surfex
 All Classes Files Functions Variables Typedefs
/home/dasprezs/EXPORT_v7_3/src/SURFEX/coare30_flux.F90
Go to the documentation of this file.
00001 !     #########
00002     SUBROUTINE COARE30_FLUX(PZ0SEA,PTA,PEXNA,PRHOA,PSST,PEXNS,PQA,  &
00003             PVMOD,PZREF,PUREF,PPS,PQSAT,PSFTH,PSFTQ,PUSTAR,PCD,PCDN,PCH,PCE,PRI,&
00004             PRESA,PRAIN,PZ0HSEA)  
00005 !     #######################################################################
00006 !
00007 !
00008 !!****  *COARE25_FLUX*  
00009 !!
00010 !!    PURPOSE
00011 !!    -------
00012 !      Calculate the surface fluxes of heat, moisture, and momentum over
00013 !      sea surface with bulk algorithm COARE3.0. 
00014 !     
00015 !!**  METHOD
00016 !!    ------
00017 !      transfer coefficients were obtained using a dataset which combined COARE
00018 !      data with those from three other ETL field experiments, and reanalysis of
00019 !      the HEXMAX data (DeCosmos et al. 1996). 
00020 !      ITERMAX=3 
00021 !      Take account of the surface gravity waves on the velocity roughness and 
00022 !      hence the momentum transfer coefficient
00023 !        NGRVWAVES=0 no gravity waves action (Charnock) !default value
00024 !        NGRVWAVES=1 wave age parameterization of Oost et al. 2002
00025 !        NGRVWAVES=2 model of Taylor and Yelland 2001
00026 !
00027 !!    EXTERNAL
00028 !!    --------
00029 !!
00030 !!    IMPLICIT ARGUMENTS
00031 !!    ------------------ 
00032 !!      
00033 !!    REFERENCE
00034 !!    ---------
00035 !!      Fairall et al (2003), J. of Climate, vol. 16, 571-591
00036 !!      Fairall et al (1996), JGR, 3747-3764
00037 !!      Gosnell et al (1995), JGR, 437-442
00038 !!      Fairall et al (1996), JGR, 1295-1308
00039 !!      
00040 !!    AUTHOR
00041 !!    ------
00042 !!     C. Lebeaupin  *Météo-France* (adapted from C. Fairall's code)
00043 !!
00044 !!    MODIFICATIONS
00045 !!    -------------
00046 !!      Original     1/06/2006
00047 !!      B. Decharme    06/2009 limitation of Ri
00048 !!      B. Decharme    09/2012 Bug in Ri calculation and limitation of Ri in surface_ri.F90
00049 !-------------------------------------------------------------------------------
00050 !
00051 !*       0.     DECLARATIONS
00052 !               ------------
00053 !
00054 USE MODD_CSTS,       ONLY : XKARMAN, XG, XSTEFAN, XRD, XRV, XPI, &
00055                               XLVTT, XCL, XCPD, XCPV, XRHOLW, XTT,XP00  
00056 USE MODD_SEAFLUX_n
00057 USE MODD_SURF_PAR,   ONLY : XUNDEF
00058 USE MODD_WATER_PAR
00059 !
00060 USE MODI_SURFACE_RI
00061 USE MODI_WIND_THRESHOLD
00062 USE MODE_COARE30_PSI
00063 !
00064 USE MODE_THERMOS
00065 !
00066 !
00067 !
00068 USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
00069 USE PARKIND1  ,ONLY : JPRB
00070 !
00071 IMPLICIT NONE
00072 !
00073 !*      0.1    declarations of arguments
00074 !
00075 !
00076 REAL, DIMENSION(:), INTENT(IN)       :: PTA   ! air temperature at atm. level (K)
00077 REAL, DIMENSION(:), INTENT(IN)       :: PQA   ! air humidity at atm. level (kg/kg)
00078 REAL, DIMENSION(:), INTENT(IN)       :: PEXNA ! Exner function at atm. level
00079 REAL, DIMENSION(:), INTENT(IN)       :: PRHOA ! air density at atm. level
00080 REAL, DIMENSION(:), INTENT(IN)       :: PVMOD ! module of wind at atm. wind level (m/s)
00081 REAL, DIMENSION(:), INTENT(IN)       :: PZREF ! atm. level for temp. and humidity (m)
00082 REAL, DIMENSION(:), INTENT(IN)       :: PUREF ! atm. level for wind (m)
00083 REAL, DIMENSION(:), INTENT(IN)       :: PSST  ! Sea Surface Temperature (K)
00084 REAL, DIMENSION(:), INTENT(IN)       :: PEXNS ! Exner function at sea surface
00085 REAL, DIMENSION(:), INTENT(IN)       :: PPS   ! air pressure at sea surface (Pa)
00086 REAL, DIMENSION(:), INTENT(IN)       :: PRAIN !precipitation rate (kg/s/m2)
00087 !
00088 REAL, DIMENSION(:), INTENT(INOUT)    :: PZ0SEA! roughness length over the ocean
00089 !                                                                                 
00090 !  surface fluxes : latent heat, sensible heat, friction fluxes
00091 REAL, DIMENSION(:), INTENT(OUT)      :: PSFTH ! heat flux (W/m2)
00092 REAL, DIMENSION(:), INTENT(OUT)      :: PSFTQ ! water flux (kg/m2/s)
00093 REAL, DIMENSION(:), INTENT(OUT)      :: PUSTAR! friction velocity (m/s)
00094 !
00095 ! diagnostics
00096 REAL, DIMENSION(:), INTENT(OUT)      :: PQSAT ! humidity at saturation
00097 REAL, DIMENSION(:), INTENT(OUT)      :: PCD   ! heat drag coefficient
00098 REAL, DIMENSION(:), INTENT(OUT)      :: PCDN  ! momentum drag coefficient
00099 REAL, DIMENSION(:), INTENT(OUT)      :: PCH   ! neutral momentum drag coefficient
00100 REAL, DIMENSION(:), INTENT(OUT)      :: PCE  !transfer coef. for latent heat flux
00101 REAL, DIMENSION(:), INTENT(OUT)      :: PRI   ! Richardson number
00102 REAL, DIMENSION(:), INTENT(OUT)      :: PRESA ! aerodynamical resistance
00103 REAL, DIMENSION(:), INTENT(OUT)      :: PZ0HSEA ! heat roughness length
00104 !
00105 !
00106 !*      0.2    declarations of local variables
00107 !
00108 REAL, DIMENSION(SIZE(PTA))      :: ZVMOD    ! wind intensity
00109 REAL, DIMENSION(SIZE(PTA))      :: ZPA      ! Pressure at atm. level
00110 REAL, DIMENSION(SIZE(PTA))      :: ZQASAT   ! specific humidity at saturation  at atm. level (kg/kg)
00111 !
00112 REAL, DIMENSION(SIZE(PTA))      :: ZO       ! rougness length ref 
00113 REAL, DIMENSION(SIZE(PTA))      :: ZWG      ! gustiness factor (m/s)
00114 !
00115 REAL, DIMENSION(SIZE(PTA))      :: ZDU,ZDT,ZDQ,ZDUWG !differences
00116 !
00117 REAL, DIMENSION(SIZE(PTA))      :: ZUSR        !velocity scaling parameter "ustar" (m/s) = friction velocity
00118 REAL, DIMENSION(SIZE(PTA))      :: ZTSR        !temperature sacling parameter "tstar" (degC)
00119 REAL, DIMENSION(SIZE(PTA))      :: ZQSR        !humidity scaling parameter "qstar" (kg/kg)
00120 !
00121 REAL, DIMENSION(SIZE(PTA))      :: ZU10,ZT10   !vertical profils (10-m height) 
00122 REAL, DIMENSION(SIZE(PTA))      :: ZVISA       !kinematic viscosity of dry air
00123 REAL, DIMENSION(SIZE(PTA))      :: ZO10,ZOT10  !roughness length at 10m
00124 REAL, DIMENSION(SIZE(PTA))      :: ZCD,ZCT,ZCC
00125 REAL, DIMENSION(SIZE(PTA))      :: ZCD10,ZCT10 !transfer coef. at 10m
00126 REAL, DIMENSION(SIZE(PTA))      :: ZRIBU,ZRIBCU
00127 REAL, DIMENSION(SIZE(PTA))      :: ZETU,ZL10
00128 !
00129 REAL, DIMENSION(SIZE(PTA))      :: ZCHARN                      !Charnock number depends on wind module
00130 REAL, DIMENSION(SIZE(PTA))      :: ZTWAVE,ZHWAVE,ZCWAVE,ZLWAVE !to compute gravity waves' impact
00131 !
00132 REAL, DIMENSION(SIZE(PTA))      :: ZZL,ZZTL!,ZZQL    !Obukhovs stability 
00133                                                      !param. z/l for u,T,q
00134 REAL, DIMENSION(SIZE(PTA))      :: ZRR
00135 REAL, DIMENSION(SIZE(PTA))      :: ZOT,ZOQ           !rougness length ref
00136 REAL, DIMENSION(SIZE(PTA))      :: ZPUZ,ZPTZ,ZPQZ    !PHI funct. for u,T,q 
00137 !
00138 REAL, DIMENSION(SIZE(PTA))      :: ZBF               !constants to compute gustiness factor
00139 !
00140 REAL, DIMENSION(SIZE(PTA))      :: ZTAU       !momentum flux (W/m2)
00141 REAL, DIMENSION(SIZE(PTA))      :: ZHF        !sensible heat flux (W/m2)
00142 REAL, DIMENSION(SIZE(PTA))      :: ZEF        !latent heat flux (W/m2)
00143 REAL, DIMENSION(SIZE(PTA))      :: ZWBAR      !diag for webb correction but not used here after
00144 REAL, DIMENSION(SIZE(PTA))      :: ZTAUR      !momentum flux due to rain (W/m2)
00145 REAL, DIMENSION(SIZE(PTA))      :: ZRF        !sensible heat flux due to rain (W/m2)
00146 REAL, DIMENSION(SIZE(PTA))      :: ZCHN,ZCEN  !neutral coef. for heat and vapor
00147 !
00148 REAL, DIMENSION(SIZE(PTA))      :: ZLV      !latent heat constant
00149 !
00150 REAL, DIMENSION(SIZE(PTA))      :: ZTAC,ZDQSDT,ZDTMP,ZDWAT,ZALFAC ! for precipitation impact
00151 REAL, DIMENSION(SIZE(PTA))      :: ZXLR                           ! vaporisation  heat  at a given temperature
00152 REAL, DIMENSION(SIZE(PTA))      :: ZCPLW                          ! specific heat for water at a given temperature 
00153 !
00154 REAL, DIMENSION(SIZE(PTA))      :: ZUSTAR2  ! square of friction velocity
00155 !
00156 REAL, DIMENSION(SIZE(PTA))      :: ZDIRCOSZW! orography slope cosine (=1 on water!)
00157 REAL, DIMENSION(SIZE(PTA))      :: ZAC      ! Aerodynamical conductance
00158 !
00159 !
00160 INTEGER, DIMENSION(SIZE(PTA))   :: ITERMAX             ! maximum number of iterations
00161 !
00162 REAL    :: ZRVSRDM1,ZRDSRV,ZR2 ! thermodynamic constants
00163 REAL    :: ZBETAGUST           !gustiness factor
00164 REAL    :: ZZBL                !atm. boundary layer depth (m)
00165 REAL    :: ZVISW               !m2/s kinematic viscosity of water
00166 REAL    :: ZS                  !height of rougness length ref
00167 REAL    :: ZCH10               !transfer coef. at 10m
00168 !
00169 INTEGER :: J, JLOOP    !loop indice
00170 REAL(KIND=JPRB) :: ZHOOK_HANDLE
00171 !
00172 !-------------------------------------------------------------------------------
00173 !
00174 !       1.     Initializations
00175 !              ---------------
00176 !
00177 !       1.1   Constants and parameters
00178 !
00179 IF (LHOOK) CALL DR_HOOK('COARE30_FLUX',0,ZHOOK_HANDLE)
00180 !
00181 ZRVSRDM1  = XRV/XRD-1. ! 0.607766
00182 ZRDSRV    = XRD/XRV    ! 0.62198
00183 ZR2       = 1.-ZRDSRV  ! pas utilisé dans cette routine
00184 ZBETAGUST = 1.2        ! value based on TOGA-COARE experiment
00185 ZZBL      = 600.       ! Set a default value for boundary layer depth
00186 ZS        = 10.        ! Standard heigth =10m
00187 ZCH10     = 0.00115
00188 !
00189 ZVISW     = 1.E-6
00190 !
00191 !       1.2   Array initialization by undefined values
00192 !
00193 PSFTH (:)=XUNDEF
00194 PSFTQ (:)=XUNDEF
00195 PUSTAR(:)=XUNDEF
00196 !
00197 PCD(:) = XUNDEF
00198 PCDN(:) = XUNDEF
00199 PCH(:) = XUNDEF
00200 PCE(:) =XUNDEF
00201 PRI(:) = XUNDEF
00202 !
00203 PRESA(:)=XUNDEF
00204 !
00205 !-------------------------------------------------------------------------------
00206 !       2. INITIAL GUESS FOR THE ITERATIVE METHOD 
00207 !          -------------------------------------
00208 !
00209 !       2.1     Wind and humidity 
00210 !
00211 ! Sea surface specific humidity 
00212 !
00213 PQSAT(:)=QSAT_SEAWATER(PSST(:),PPS(:))         
00214 !              
00215 ! Set a minimum value to wind 
00216 !
00217 ZVMOD(:) = WIND_THRESHOLD(PVMOD(:),PUREF(:))
00218 !
00219 ! Specific humidity at saturation at the atm. level 
00220 !
00221 ZPA(:) = XP00* (PEXNA(:)**(XCPD/XRD))
00222 ZQASAT(:) = QSAT(PTA(:),ZPA(:)) 
00223 !
00224 !
00225 ZO(:)  = 0.0001
00226 ZWG(:) = 0.
00227 IF (LPWG) ZWG(:) = 0.5
00228 !
00229 DO J=1,SIZE(PTA)
00230   !
00231   !      2.2       initial guess
00232   !    
00233   ZDU(J) = ZVMOD(J)   !wind speed difference with surface current(=0) (m/s)
00234                       !initial guess for gustiness factor
00235   ZDT(J) = -(PTA(J)/PEXNA(J)) + (PSST(J)/PEXNS(J)) !potential temperature difference
00236   ZDQ(J) = PQSAT(J)-PQA(J)                         !specific humidity difference
00237   !
00238   ZDUWG(J) = (ZDU(J)*ZDU(J)+ZWG(J)*ZWG(J))**0.5    !wind speed difference including gustiness ZWG
00239   !
00240   !      2.3   initialization of neutral coefficients
00241   !
00242   ZU10(J)  = ZDUWG(J)*LOG(ZS/ZO(J))/LOG(PUREF(J)/ZO(J))
00243   ZUSR(J)  = 0.035*ZU10(J)
00244   ZVISA(J) = 1.326E-5*(1.+6.542E-3*(PTA(J)-XTT)+&
00245              8.301E-6*(PTA(J)-XTT)**2.-4.84E-9*(PTA(J)-XTT)**3.) !Andrea (1989) CRREL Rep. 89-11
00246   ! 
00247   ZO10(J) = 0.011*ZUSR(J)*ZUSR(J)/XG+0.11*ZVISA(J)/ZUSR(J)
00248   ZCD(J)  = (XKARMAN/LOG(PUREF(J)/ZO10(J)))**2  !drag coefficient
00249   ZCD10(J)= (XKARMAN/LOG(ZS/ZO10(J)))**2
00250   ZCT10(J)= ZCH10/SQRT(ZCD10(J))
00251   ZOT10(J)= ZS/EXP(XKARMAN/ZCT10(J))
00252   !
00253   !-------------------------------------------------------------------------------
00254   !             Grachev and Fairall (JAM, 1997)
00255   ZCT(J) = XKARMAN/LOG(PZREF(J)/ZOT10(J))      !temperature transfer coefficient
00256   ZCC(J) = XKARMAN*ZCT(J)/ZCD(J)               !z/L vs Rib linear coef.
00257   !
00258   ZRIBCU(J) = -PUREF(J)/(ZZBL*0.004*ZBETAGUST**3) !saturation or plateau Rib
00259   !ZRIBU(J) =-XG*PUREF(J)*(ZDT(J)+ZRVSRDM1*(PTA(J)-XTT)*ZDQ)/&
00260   !     &((PTA(J)-XTT)*ZDUWG(J)**2)
00261   ZRIBU(J)  = -XG*PUREF(J)*(ZDT(J)+ZRVSRDM1*PTA(J)*ZDQ(J))/&
00262                (PTA(J)*ZDUWG(J)**2)  
00263   !
00264   IF (ZRIBU(J)<0.) THEN
00265     ZETU(J) = ZCC(J)*ZRIBU(J)/(1.+ZRIBU(J)/ZRIBCU(J))    !Unstable G and F
00266   ELSE
00267     ZETU(J) = ZCC(J)*ZRIBU(J)/(1.+27./9.*ZRIBU(J)/ZCC(J))!Stable 
00268   ENDIF
00269   !
00270   ZL10(J) = PUREF(J)/ZETU(J) !MO length
00271   !
00272 ENDDO
00273 !
00274 !  First guess M-O stability dependent scaling params. (u*,T*,q*) to estimate ZO and z/L (ZZL)
00275 ZUSR(:) = ZDUWG(:)*XKARMAN/(LOG(PUREF(:)/ZO10(:))-PSIFCTU(PUREF(:)/ZL10(:)))
00276 ZTSR(:) = -ZDT(:)*XKARMAN/(LOG(PZREF(:)/ZOT10(:))-PSIFCTT(PZREF(:)/ZL10(:)))
00277 ZQSR(:) = -ZDQ(:)*XKARMAN/(LOG(PZREF(:)/ZOT10(:))-PSIFCTT(PZREF(:)/ZL10(:)))
00278 !
00279 ZCHARN(:) = 0.011  
00280 !
00281 ZZL(:) = 0.0
00282 !
00283 DO J=1,SIZE(PTA)
00284   !
00285   IF (ZETU(J)>50.) THEN
00286     ITERMAX(J) = 1
00287   ELSE
00288     ITERMAX(J) = 3 !number of iterations
00289   ENDIF
00290   !
00291   !then modify Charnork for high wind speeds Chris Fairall's data
00292   IF (ZDUWG(J)>10.) ZCHARN(J) = 0.011 + (0.018-0.011)*(ZDUWG(J)-10.)/(18-10)
00293   IF (ZDUWG(J)>18.) ZCHARN(J) = 0.018
00294   !
00295   !                3.  ITERATIVE LOOP TO COMPUTE USR, TSR, QSR 
00296   !                -------------------------------------------
00297   !
00298   ZHWAVE(J) = 0.018*PVMOD(J)*PVMOD(J)*(1.+0.015*PVMOD(J))
00299   ZTWAVE(J) = 0.729*PVMOD(J)
00300   ZCWAVE(J) = XG*ZTWAVE(J)/(2.*XPI)
00301   ZLWAVE(J) = ZTWAVE(J)*ZCWAVE(J)
00302   !
00303 ENDDO
00304 !
00305    
00306 !
00307 DO JLOOP=1,MAXVAL(ITERMAX) ! begin of iterative loop
00308   !
00309   DO J=1,SIZE(PTA)
00310     !
00311     IF (JLOOP.GT.ITERMAX(J)) CYCLE
00312     !
00313     IF (NGRVWAVES==0) THEN
00314       ZO(J) = ZCHARN(J)*ZUSR(J)*ZUSR(J)/XG + 0.11*ZVISA(J)/ZUSR(J) !Smith 1988
00315     ELSE IF (NGRVWAVES==1) THEN
00316       ZO(J) = (50./(2.*XPI))*ZLWAVE(J)*(ZUSR(J)/ZCWAVE(J))**4.5 &
00317               + 0.11*ZVISA(J)/ZUSR(J)                       !Oost et al. 2002  
00318     ELSE IF (NGRVWAVES==2) THEN
00319       ZO(J) = 1200.*ZHWAVE(J)*(ZHWAVE(J)/ZLWAVE(J))**4.5 &
00320               + 0.11*ZVISA(J)/ZUSR(J)                       !Taulor and Yelland 2001  
00321     ENDIF
00322     !
00323     ZRR(J) = ZO(J)*ZUSR(J)/ZVISA(J)
00324     ZOQ(J) = MIN(1.15E-4 , 5.5E-5/ZRR(J)**0.6)
00325     ZOT(J) = ZOQ(J)
00326     !
00327     ZZL(J) = XKARMAN * XG * PUREF(J) * &
00328               ( ZTSR(J)*(1.+ZRVSRDM1*PQA(J)) + ZRVSRDM1*PTA(J)*ZQSR(J) ) / &
00329               ( PTA(J)*ZUSR(J)*ZUSR(J)*(1.+ZRVSRDM1*PQA(J)) )  
00330     ZZTL(J)= ZZL(J)*PZREF(J)/PUREF(J)  ! for T 
00331 !    ZZQL(J)=ZZL(J)*PZREF(J)/PUREF(J)  ! for Q
00332   ENDDO
00333   !
00334   ZPUZ(:) = PSIFCTU(ZZL(:))     
00335   ZPTZ(:) = PSIFCTT(ZZTL(:))
00336   !
00337   DO J=1,SIZE(PTA)
00338     !
00339     ! ZPQZ(J)=PSIFCTT(ZZQL(J))    
00340     ZPQZ(J) = ZPTZ(J)
00341     !
00342     !             3.1 scale parameters
00343     !
00344     ZUSR(J) = ZDUWG(J)*XKARMAN/(LOG(PUREF(J)/ZO(J)) -ZPUZ(J))
00345     ZTSR(J) = -ZDT(J)  *XKARMAN/(LOG(PZREF(J)/ZOT(J))-ZPTZ(J))
00346     ZQSR(J) = -ZDQ(J)  *XKARMAN/(LOG(PZREF(J)/ZOQ(J))-ZPQZ(J))
00347     !
00348     !             3.2 Gustiness factor (ZWG)
00349     !
00350     IF(LPWG) THEN
00351       ZBF(J) = -XG/PTA(J)*ZUSR(J)*(ZTSR(J)+ZRVSRDM1*PTA(J)*ZQSR(J))
00352       IF (ZBF(J)>0.) THEN
00353         ZWG(J) = ZBETAGUST*(ZBF(J)*ZZBL)**(1./3.)
00354       ELSE
00355         ZWG(J) = 0.2
00356       ENDIF
00357     ENDIF  
00358     ZDUWG(J) = SQRT(ZVMOD(J)**2. + ZWG(J)**2.)
00359     !
00360   ENDDO
00361   !
00362 ENDDO
00363 !-------------------------------------------------------------------------------
00364 !
00365 !            4.  COMPUTE transfer coefficients PCD, PCH, ZCE and SURFACE FLUXES
00366 !                --------------------------------------------------------------
00367 !
00368 ZTAU(:) = XUNDEF
00369 ZHF(:)  = XUNDEF
00370 ZEF(:)  = XUNDEF
00371 !
00372 ZWBAR(:) = 0.
00373 ZTAUR(:) = 0.
00374 ZRF(:)   = 0.
00375 !
00376 DO J=1,SIZE(PTA)
00377   !
00378   !            4. 0  roughness length over the ocean PZOSEA
00379   !
00380   IF(ZUSR(J)/=0.0) PZ0SEA(J) = 10./EXP(XKARMAN*ZU10(J)/ZUSR(J))
00381   !
00382   !            4. 1 transfert coefficients PCD, PCH and PCE 
00383   !                 and neutral PCDN, ZCHN, ZCEN 
00384   !
00385   PCD(J) = (ZUSR(J)/ZDUWG(J))**2.
00386   PCH(J) = ZUSR(J)*ZTSR(J)/(ZDUWG(J)*(PTA(J)*PEXNS(J)/PEXNA(J)-PSST(J)))
00387   PCE(J) = ZUSR(J)*ZQSR(J)/(ZDUWG(J)*(PQA(J)-PQSAT(J)))
00388   !
00389   PCDN(J) = (XKARMAN/LOG(ZS/ZO(J)))**2.
00390   ZCHN(J) = (XKARMAN/LOG(ZS/ZO(J)))*(XKARMAN/LOG(ZS/ZOT(J)))
00391   ZCEN(J) = (XKARMAN/LOG(ZS/ZO(J)))*(XKARMAN/LOG(ZS/ZOQ(J)))
00392   !
00393   ZLV(J) = XLVTT + (XCPV-XCL)*(PSST(J)-XTT)
00394   !
00395   !            4. 2 surface fluxes 
00396   !
00397   IF (ABS(PCDN(J))<=1.E-2) THEN   !!!! secure COARE3.0 CODE 
00398     ZTSR(J) = -ZTSR(J)
00399     ZQSR(J) = -ZQSR(J)
00400     ZTAU(J) = -PRHOA(J)*ZUSR(J)*ZUSR(J)*ZVMOD(J)/ZDUWG(J)
00401     ZHF(J)  =  PRHOA(J)*XCPD*ZUSR(J)*ZTSR(J)
00402     ZEF(J)  =  PRHOA(J)*ZLV(J)*ZUSR(J)*ZQSR(J)
00403     !    
00404     !           4.3 Contributions to surface  fluxes due to rainfall
00405     !
00406     ! SB: a priori, le facteur ZRDSRV=XRD/XRV est introduit pour
00407     !     adapter la formule de Clausius-Clapeyron (pour l'air
00408     !     sec) au cas humide.
00409     IF (LPRECIP) THEN
00410       ! 
00411       ! heat surface  fluxes
00412       !
00413       ZTAC(J)  = PTA(J)-XTT
00414       !
00415       ZXLR(J)  = XLVTT + (XCPV-XCL)* ZTAC(J)                            ! latent heat of rain vaporization
00416       ZDQSDT(J)= ZQASAT(J) * ZXLR(J) / (XRD*PTA(J)**2)                  ! Clausius-Clapeyron relation
00417       ZDTMP(J) = (1.0 + 3.309e-3*ZTAC(J) -1.44e-6*ZTAC(J)*ZTAC(J)) * &  !heat diffusivity
00418                   0.02411 / (PRHOA(J)*XCPD)
00419       !
00420       ZDWAT(J) = 2.11e-5 * (XP00/ZPA(J)) * (PTA(J)/XTT)**1.94           ! water vapour diffusivity from eq (13.3)
00421       !                                                                 ! of Pruppacher and Klett (1978)      
00422       ZALFAC(J)= 1.0 / (1.0 + &                                         ! Eq.11 in GoF95
00423                    ZRDSRV*ZDQSDT(J)*ZXLR(J)*ZDWAT(J)/(ZDTMP(J)*XCPD))   ! ZALFAC=wet-bulb factor (sans dim)     
00424       ZCPLW(J) = 4224.8482 + ZTAC(J) * &
00425                               ( -4.707 + ZTAC(J) * &
00426                                 (0.08499 + ZTAC(J) * &
00427                                   (1.2826e-3 + ZTAC(J) * &
00428                                     (4.7884e-5 - 2.0027e-6* ZTAC(J))))) ! specific heat  
00429       !       
00430       ZRF(J)   = PRAIN(J) * ZCPLW(J) * ZALFAC(J) * &                    !Eq.12 in GoF95 !SIGNE?
00431                    (PSST(J) - PTA(J) + (PQSAT(J)-PQA(J))*ZXLR(J)/XCPD )
00432       !
00433       ! Momentum flux due to rainfall  
00434       !
00435       ZTAUR(J)=-0.85*(PRAIN(J) *ZVMOD(J)) !pp3752 in FBR96
00436       !
00437     ENDIF
00438     !
00439     !             4.4   Webb correction to latent heat flux
00440     ! 
00441     ZWBAR(J)=- (1./ZRDSRV)*ZUSR(J)*ZQSR(J) / (1.0+(1./ZRDSRV)*PQA(J)) &
00442                - ZUSR(J)*ZTSR(J)/PTA(J)                        ! Eq.21*rhoa in FBR96    
00443     !
00444     !             4.5   friction velocity which contains correction du to rain            
00445     !
00446     ZUSTAR2(J)= - (ZTAU(J) + ZTAUR(J)) / PRHOA(J)
00447     PUSTAR(J) =  SQRT(ZUSTAR2(J))
00448     !
00449     !             4.6   Total surface fluxes
00450     !           
00451     PSFTH (J) =  ZHF(J) + ZRF(J)
00452     PSFTQ (J) =  ZEF(J) / ZLV(J)
00453     ! 
00454   ENDIF
00455 ENDDO                      
00456 !-------------------------------------------------------------------------------
00457 !
00458 !       5.  FINAL STEP : TOTAL SURFACE FLUXES AND DERIVED DIAGNOSTICS 
00459 !           -----------
00460 !       5.1    Richardson number
00461 !             
00462 !
00463 ZDIRCOSZW(:) = 1.
00464  CALL SURFACE_RI(PSST,PQSAT,PEXNS,PEXNA,PTA,ZQASAT,  &
00465                   PZREF,PUREF,ZDIRCOSZW,PVMOD,PRI)  
00466 !
00467 !       5.2     Aerodynamical conductance and resistance
00468 !             
00469 ZAC(:) = PCH(:)*ZVMOD(:)
00470 PRESA(:) = 1. / ZAC(:)
00471 !
00472 PZ0HSEA(:) = PZ0SEA(:)
00473 IF (LHOOK) CALL DR_HOOK('COARE30_FLUX',1,ZHOOK_HANDLE)
00474 !
00475 !-------------------------------------------------------------------------------
00476 !
00477 END SUBROUTINE COARE30_FLUX