source: LMDZ6/trunk/libf/phylmd/coare30_flux_cnrm_mod.f90 @ 5279

Last change on this file since 5279 was 5274, checked in by abarral, 9 months ago

Replace yomcst.h by existing module
File size: 21.0 KB
Line 
1!SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
2!SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
3!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt 
4!SFX_LIC for details. version 1.
5!     #########
6
7module coare30_flux_cnrm_mod
8  implicit none
9  private
10  public COARE30_FLUX_CNRM
11
12contains
13
14
15SUBROUTINE COARE30_FLUX_CNRM(PZ0SEA,PTA,PSST,PQA,  &
16            PVMOD,PZREF,PUREF,PPS,PQSATA,PQSAT,PSFTH,PSFTQ,PUSTAR,PCD,PCDN,PCH,PCE,PRI,&
17            PRESA,PRAIN,PPA,PZ0HSEA,LPRECIP, LPWG,coeffs )
18!     #######################################################################
19!
20!
21!!****  *COARE25_FLUX*
22!!
23!!    PURPOSE
24!!    -------
25!      Calculate the surface fluxes of heat, moisture, and momentum over
26!      sea surface with bulk algorithm COARE3.0.
27!
28!!**  METHOD
29!!    ------
30!      transfer coefficients were obtained using a dataset which combined COARE
31!      data with those from three other ETL field experiments, and reanalysis of
32!      the HEXMAX data (DeCosmos et al. 1996).
33!      ITERMAX=3
34!      Take account of the surface gravity waves on the velocity roughness and
35!      hence the momentum transfer coefficient
36!        NGRVWAVES=0 no gravity waves action (Charnock) !default value
37!        NGRVWAVES=1 wave age parameterization of Oost et al. 2002
38!        NGRVWAVES=2 model of Taylor and Yelland 2001
39!
40!!    EXTERNAL
41!!    --------
42!!
43!!    IMPLICIT ARGUMENTS
44!!    ------------------
45!!
46!!    REFERENCE
47!!    ---------
48!!      Fairall et al (2003), J. of Climate, vol. 16, 571-591
49!!      Fairall et al (1996), JGR, 3747-3764
50!!      Gosnell et al (1995), JGR, 437-442
51!!      Fairall et al (1996), JGR, 1295-1308
52!!
53!!    AUTHOR
54!!    ------
55!!     C. Lebeaupin  *Météo-France* (adapted from C. Fairall's code)
56!!
57!!    MODIFICATIONS
58!!    -------------
59!!      Original     1/06/2006
60!!      B. Decharme    06/2009 limitation of Ri
61!!      B. Decharme    09/2012 Bug in Ri calculation and limitation of Ri in surface_ri.F90
62!!      B. Decharme    06/2013 bug in z0 (output) computation
63!!      J.Escobar      06/2013  for REAL4/8 add EPSILON management
64!!      C. Lebeaupin   03/2014 bug if PTA=PSST and PEXNA=PEXNS: set a minimum value
65!!                             add abort if no convergence
66!-------------------------------------------------------------------------------
67!
68!*       0.     DECLARATIONS
69!               ------------
70!
71!
72!USE MODD_SEAFLUX_n, ONLY : SEAFLUX_t
73!
74
75!----------Rajout Olive ---------
76USE dimphy
77USE indice_sol_mod
78USE coare_cp_mod, ONLY: PSIFCTT=>psit_30, PSIFCTU=>psiuo
79
80!--------------------------------
81
82USE MODD_CSTS,       ONLY : XKARMAN, XG, XSTEFAN, XRD, XRV, XPI, &
83                            XLVTT, XCL, XCPD, XCPV, XRHOLW, XTT, &
84                            XP00
85
86!USE MODD_SURF_ATM,   ONLY : XVZ0CM
87!
88!USE MODD_SURF_PAR,   ONLY : XUNDEF, XSURF_EPSILON
89!USE MODD_WATER_PAR
90!
91!USE MODI_SURFACE_RI
92!USE MODI_WIND_THRESHOLD
93!USE MODE_COARE30_PSI
94!
95!USE MODE_THERMOS
96!
97!
98!USE MODI_ABOR1_SFX
99!
100!
101!USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
102!USE PARKIND1  ,ONLY : JPRB
103USE yomcst_mod_h, ONLY: RPI, RCLUM, RHPLA, RKBOL, RNAVO                   &
104          , RDAY, REA, REPSM, RSIYEA, RSIDAY, ROMEGA                  &
105          , R_ecc, R_peri, R_incl                                      &
106          , RA, RG, R1SA                                         &
107          , RSIGMA                                                     &
108          , R, RMD, RMV, RD, RV, RCPD                    &
109          , RMO3, RMCO2, RMC, RMCH4, RMN2O, RMCFC11, RMCFC12        &
110          , RCPV, RCVD, RCVV, RKAPPA, RETV, eps_w                    &
111          , RCW, RCS                                                 &
112          , RLVTT, RLSTT, RLMLT, RTT, RATM                           &
113          , RESTT, RALPW, RBETW, RGAMW, RALPS, RBETS, RGAMS            &
114          , RALPD, RBETD, RGAMD
115IMPLICIT NONE
116!
117!*      0.1    declarations of arguments
118!
119!
120!
121!TYPE(SEAFLUX_t), INTENT(INOUT) :: S
122!
123REAL, DIMENSION(klon), INTENT(IN)       :: PTA   ! air temperature at atm. level (K)
124REAL, DIMENSION(klon), INTENT(IN)       :: PQA   ! air humidity at atm. level (kg/kg)
125!REAL, DIMENSION(:), INTENT(IN)       :: PEXNA ! Exner function at atm. level
126!REAL, DIMENSION(:), INTENT(IN)       :: PRHOA ! air density at atm. level
127REAL, DIMENSION(klon), INTENT(IN)       :: PVMOD ! module of wind at atm. wind level (m/s)
128REAL, DIMENSION(klon), INTENT(IN)       :: PZREF ! atm. level for temp. and humidity (m)
129REAL, DIMENSION(klon), INTENT(IN)       :: PUREF ! atm. level for wind (m)
130REAL, DIMENSION(klon), INTENT(IN)       :: PSST  ! Sea Surface Temperature (K)
131!REAL, DIMENSION(:), INTENT(IN)       :: PEXNS ! Exner function at sea surface
132REAL, DIMENSION(klon), INTENT(IN)       :: PPS   ! air pressure at sea surface (Pa)
133REAL, DIMENSION(klon), INTENT(IN)       :: PRAIN !precipitation rate (kg/s/m2)
134REAL, DIMENSION(klon), INTENT(IN)       :: PPA        ! air pressure at atm level (Pa)
135REAL, DIMENSION(klon), INTENT(IN)       :: PQSATA        ! air pressure at atm level (Pa)
136!
137REAL, DIMENSION(klon), INTENT(INOUT)    :: PZ0SEA! roughness length over the ocean
138!
139!  surface fluxes : latent heat, sensible heat, friction fluxes
140REAL, DIMENSION(klon), INTENT(OUT)      :: PSFTH ! heat flux (W/m2)
141REAL, DIMENSION(klon), INTENT(OUT)      :: PSFTQ ! water flux (kg/m2/s)
142REAL, DIMENSION(klon), INTENT(OUT)      :: PUSTAR! friction velocity (m/s)
143!
144! diagnostics
145REAL, DIMENSION(klon), INTENT(OUT)      :: PQSAT ! humidity at saturation
146REAL, DIMENSION(klon), INTENT(OUT)      :: PCD   ! heat drag coefficient
147REAL, DIMENSION(klon), INTENT(OUT)      :: PCDN  ! momentum drag coefficient
148REAL, DIMENSION(klon), INTENT(OUT)      :: PCH   ! neutral momentum drag coefficient
149REAL, DIMENSION(klon), INTENT(OUT)      :: PCE  !transfer coef. for latent heat flux
150REAL, DIMENSION(klon), INTENT(OUT)      :: PRI   ! Richardson number
151REAL, DIMENSION(klon), INTENT(OUT)      :: PRESA ! aerodynamical resistance
152REAL, DIMENSION(klon), INTENT(OUT)      :: PZ0HSEA ! heat roughness length
153
154LOGICAL,            INTENT(IN)    :: LPRECIP   !
155LOGICAL,            INTENT(IN)    :: LPWG    !
156real, dimension(3), intent(inout)      :: coeffs
157
158!*      0.2    declarations of local variables
159!
160REAL, DIMENSION(SIZE(PTA))      :: ZVMOD    ! wind intensity
161REAL, DIMENSION(SIZE(PTA))      :: ZPA      ! Pressure at atm. level
162REAL, DIMENSION(SIZE(PTA))      :: ZTA      ! Temperature at atm. level
163REAL, DIMENSION(SIZE(PTA))      :: ZQASAT   ! specific humidity at saturation  at atm. level (kg/kg)
164!
165!rajout
166REAL, DIMENSION(SIZE(PTA))       :: PEXNA      ! Exner function at atm level
167REAL, DIMENSION(SIZE(PTA))       :: PEXNS      ! Exner function at atm level
168!
169REAL, DIMENSION(SIZE(PTA))      :: ZO       ! rougness length ref
170REAL, DIMENSION(SIZE(PTA))      :: ZWG      ! gustiness factor (m/s)
171!
172REAL, DIMENSION(SIZE(PTA))      :: ZDU,ZDT,ZDQ,ZDUWG !differences
173!
174REAL, DIMENSION(SIZE(PTA))      :: ZUSR        !velocity scaling parameter "ustar" (m/s) = friction velocity
175REAL, DIMENSION(SIZE(PTA))      :: ZTSR        !temperature sacling parameter "tstar" (degC)
176REAL, DIMENSION(SIZE(PTA))      :: ZQSR        !humidity scaling parameter "qstar" (kg/kg)
177!
178REAL, DIMENSION(SIZE(PTA))      :: ZU10,ZT10   !vertical profils (10-m height)
179REAL, DIMENSION(SIZE(PTA))      :: ZVISA       !kinematic viscosity of dry air
180REAL, DIMENSION(SIZE(PTA))      :: ZO10,ZOT10  !roughness length at 10m
181REAL, DIMENSION(SIZE(PTA))      :: ZCD,ZCT,ZCC
182REAL, DIMENSION(SIZE(PTA))      :: ZCD10,ZCT10 !transfer coef. at 10m
183REAL, DIMENSION(SIZE(PTA))      :: ZRIBU,ZRIBCU
184REAL, DIMENSION(SIZE(PTA))      :: ZETU,ZL10
185!
186REAL, DIMENSION(SIZE(PTA))      :: ZCHARN                      !Charnock number depends on wind module
187REAL, DIMENSION(SIZE(PTA))      :: ZTWAVE,ZHWAVE,ZCWAVE,ZLWAVE !to compute gravity waves' impact
188!
189REAL, DIMENSION(SIZE(PTA))      :: ZZL,ZZTL!,ZZQL    !Obukhovs stability
190                                                     !param. z/l for u,T,q
191REAL, DIMENSION(SIZE(PTA))      :: ZRR
192REAL, DIMENSION(SIZE(PTA))      :: ZOT,ZOQ           !rougness length ref
193REAL, DIMENSION(SIZE(PTA))      :: ZPUZ,ZPTZ,ZPQZ    !PHI funct. for u,T,q
194!
195REAL, DIMENSION(SIZE(PTA))      :: ZBF               !constants to compute gustiness factor
196!
197REAL, DIMENSION(SIZE(PTA))      :: ZTAU       !momentum flux (W/m2)
198REAL, DIMENSION(SIZE(PTA))      :: ZHF        !sensible heat flux (W/m2)
199REAL, DIMENSION(SIZE(PTA))      :: ZEF        !latent heat flux (W/m2)
200REAL, DIMENSION(SIZE(PTA))      :: ZWBAR      !diag for webb correction but not used here after
201REAL, DIMENSION(SIZE(PTA))      :: ZTAUR      !momentum flux due to rain (W/m2)
202REAL, DIMENSION(SIZE(PTA))      :: ZRF        !sensible heat flux due to rain (W/m2)
203REAL, DIMENSION(SIZE(PTA))      :: ZCHN,ZCEN  !neutral coef. for heat and vapor
204!
205REAL, DIMENSION(SIZE(PTA))      :: ZLV      !latent heat constant
206!
207REAL, DIMENSION(SIZE(PTA))      :: ZTAC,ZDQSDT,ZDTMP,ZDWAT,ZALFAC ! for precipitation impact
208REAL, DIMENSION(SIZE(PTA))      :: ZXLR                           ! vaporisation  heat  at a given temperature
209REAL, DIMENSION(SIZE(PTA))      :: ZCPLW                          ! specific heat for water at a given temperature
210!
211REAL, DIMENSION(SIZE(PTA))      :: ZUSTAR2  ! square of friction velocity
212!
213REAL, DIMENSION(SIZE(PTA))      :: ZDIRCOSZW! orography slope cosine (=1 on water!)
214REAL, DIMENSION(SIZE(PTA))      :: ZAC      ! Aerodynamical conductance
215!
216!
217INTEGER, DIMENSION(SIZE(PTA))   :: ITERMAX             ! maximum number of iterations
218!
219REAL    :: ZRVSRDM1,ZRDSRV,ZR2 ! thermodynamic constants
220REAL    :: ZBETAGUST           !gustiness factor
221REAL    :: ZZBL                !atm. boundary layer depth (m)
222REAL    :: ZVISW               !m2/s kinematic viscosity of water
223REAL    :: ZS                  !height of rougness length ref
224REAL    :: ZCH10               !transfer coef. at 10m
225
226REAL    :: QSAT_SEAWATER
227REAL    :: QSATSEAW_1D
228!
229INTEGER :: J, JLOOP    !loop indice
230!REAL(KIND=JPRB) :: ZHOOK_HANDLE
231
232!--------- Modif Olive -----------------
233REAL, DIMENSION(SIZE(PTA))        :: PRHOA
234REAL, PARAMETER                   :: XUNDEF = 1.E+20
235
236REAL       :: XVCHRNK = 0.021
237REAL       :: XVZ0CM = 1.0E-5
238!REAL       :: XRIMAX
239
240
241INTEGER :: PREF             ! reference pressure for exner function
242INTEGER :: NGRVWAVES        ! Pour le choix du z0
243
244INCLUDE "clesphys.h"
245
246!--------------------------------------
247
248
249PRHOA(:) = PPS(:) / (287.1 * PTA(:) * (1.+.61*PQA(:)))
250
251PREF = 100000.                     ! = 1000 hPa
252NGRVWAVES = 1
253
254PEXNA = (PPA/PREF)**(RD/RCPD)
255PEXNS = (PPS/PREF)**(RD/RCPD)
256
257!
258!-------------------------------------------------------------------------------
259!
260!       1.     Initializations
261!              ---------------
262!
263!       1.1   Constants and parameters
264!
265!IF (LHOOK) CALL DR_HOOK('COARE30_FLUX',0,ZHOOK_HANDLE)
266!
267ZRVSRDM1  = XRV/XRD-1. ! 0.607766
268ZRDSRV    = XRD/XRV    ! 0.62198
269ZR2       = 1.-ZRDSRV  ! pas utilisé dans cette routine
270ZBETAGUST = 1.2        ! value based on TOGA-COARE experiment
271ZZBL      = 600.       ! Set a default value for boundary layer depth
272ZS        = 10.        ! Standard heigth =10m
273ZCH10     = 0.00115
274!
275ZVISW     = 1.E-6
276!
277!       1.2   Array initialization by undefined values
278!
279PSFTH (:)=XUNDEF
280PSFTQ (:)=XUNDEF
281PUSTAR(:)=XUNDEF
282!
283PCD(:) = XUNDEF
284PCDN(:) = XUNDEF
285PCH(:) = XUNDEF
286PCE(:) =XUNDEF
287PRI(:) = XUNDEF
288!
289PRESA(:)=XUNDEF
290!
291!-------------------------------------------------------------------------------
292!       2. INITIAL GUESS FOR THE ITERATIVE METHOD
293!          -------------------------------------
294!
295!       2.0     Temperature
296!
297! Set a non-zero value for the temperature gradient
298!
299WHERE((PTA(:)*PEXNS(:)/PEXNA(:)-PSST(:))==0.)
300      ZTA(:)=PTA(:)-1E-3
301ELSEWHERE
302      ZTA(:)=PTA(:)
303ENDWHERE
304
305!       2.1     Wind and humidity
306!
307! Sea surface specific humidity
308!
309!PQSAT(:)=QSAT_SEAWATER(PSST(:),PPS(:))
310PQSAT(:)=QSATSEAW_1D(PSST(:),PPS(:))
311
312
313
314!
315! Set a minimum value to wind
316!
317!ZVMOD(:) = WIND_THRESHOLD(PVMOD(:),PUREF(:))
318
319
320ZVMOD = MAX(PVMOD , 0.1 * MIN(10.,PUREF) )    !set a minimum value to wind
321!ZVMOD = PVMOD    !set a minimum value to wind
322
323!
324! Specific humidity at saturation at the atm. level
325!
326ZPA(:) = XP00* (PEXNA(:)**(XCPD/XRD))
327!ZQASAT(:) = QSAT_SEAWATER(ZTA(:),ZPA(:))
328ZQASAT = QSATSEAW_1D(ZTA(:),ZPA(:))
329
330
331!
332!
333ZO(:)  = 0.0001
334ZWG(:) = 0.
335IF (LPWG) ZWG(:) = 0.5
336!
337ZCHARN(:) = 0.011
338!
339DO J=1,SIZE(PTA)
340  !
341  !      2.2       initial guess
342  !
343  ZDU(J) = ZVMOD(J)   !wind speed difference with surface current(=0) (m/s)
344                      !initial guess for gustiness factor
345  ZDT(J) = -(ZTA(J)/PEXNA(J)) + (PSST(J)/PEXNS(J)) !potential temperature difference
346  ZDQ(J) = PQSAT(J)-PQA(J)                         !specific humidity difference
347  !
348  ZDUWG(J) = SQRT(ZDU(J)**2+ZWG(J)**2)     !wind speed difference including gustiness ZWG
349  !
350  !      2.3   initialization of neutral coefficients
351  !
352  ZU10(J)  = ZDUWG(J)*LOG(ZS/ZO(J))/LOG(PUREF(J)/ZO(J))
353  ZUSR(J)  = 0.035*ZU10(J)
354  ZVISA(J) = 1.326E-5*(1.+6.542E-3*(ZTA(J)-XTT)+&
355             8.301E-6*(ZTA(J)-XTT)**2-4.84E-9*(ZTA(J)-XTT)**3) !Andrea (1989) CRREL Rep. 89-11
356  !
357  ZO10(J) = ZCHARN(J)*ZUSR(J)*ZUSR(J)/XG+0.11*ZVISA(J)/ZUSR(J)
358  ZCD(J)  = (XKARMAN/LOG(PUREF(J)/ZO10(J)))**2  !drag coefficient
359  ZCD10(J)= (XKARMAN/LOG(ZS/ZO10(J)))**2
360  ZCT10(J)= ZCH10/SQRT(ZCD10(J))
361  ZOT10(J)= ZS/EXP(XKARMAN/ZCT10(J))
362  !
363  !-------------------------------------------------------------------------------
364  !             Grachev and Fairall (JAM, 1997)
365  ZCT(J) = XKARMAN/LOG(PZREF(J)/ZOT10(J))      !temperature transfer coefficient
366  ZCC(J) = XKARMAN*ZCT(J)/ZCD(J)               !z/L vs Rib linear coef.
367  !
368  ZRIBCU(J) = -PUREF(J)/(ZZBL*0.004*ZBETAGUST**3) !saturation or plateau Rib
369  !ZRIBU(J) =-XG*PUREF(J)*(ZDT(J)+ZRVSRDM1*(ZTA(J)-XTT)*ZDQ)/&
370  !     &((ZTA(J)-XTT)*ZDUWG(J)**2)
371  ZRIBU(J)  = -XG*PUREF(J)*(ZDT(J)+ZRVSRDM1*ZTA(J)*ZDQ(J))/&
372               (ZTA(J)*ZDUWG(J)**2)
373  !
374  IF (ZRIBU(J)<0.) THEN
375    ZETU(J) = ZCC(J)*ZRIBU(J)/(1.+ZRIBU(J)/ZRIBCU(J))    !Unstable G and F
376  ELSE
377    ZETU(J) = ZCC(J)*ZRIBU(J)/(1.+27./9.*ZRIBU(J)/ZCC(J))!Stable
378  ENDIF
379  !
380  ZL10(J) = PUREF(J)/ZETU(J) !MO length
381  !
382ENDDO
383!
384!  First guess M-O stability dependent scaling params. (u*,T*,q*) to estimate ZO and z/L (ZZL)
385ZUSR(:) = ZDUWG(:)*XKARMAN/(LOG(PUREF(:)/ZO10(:))-PSIFCTU(PUREF(1)/ZL10(1)))
386ZTSR(:) = -ZDT(:)*XKARMAN/(LOG(PZREF(:)/ZOT10(:))-PSIFCTT(PZREF(1)/ZL10(1)))
387ZQSR(:) = -ZDQ(:)*XKARMAN/(LOG(PZREF(:)/ZOT10(:))-PSIFCTT(PZREF(1)/ZL10(1)))
388!
389ZZL(:) = 0.0
390!
391DO J=1,SIZE(PTA)
392  !
393  IF (ZETU(J)>50.) THEN
394    ITERMAX(J) = 1
395  ELSE
396    ITERMAX(J) = 3 !number of iterations
397  ENDIF
398  !
399  !then modify Charnork for high wind speeds Chris Fairall's data
400  IF (ZDUWG(J)>10.) ZCHARN(J) = 0.011 + (0.018-0.011)*(ZDUWG(J)-10.)/(18.-10.)
401  IF (ZDUWG(J)>18.) ZCHARN(J) = 0.018
402  !
403  !                3.  ITERATIVE LOOP TO COMPUTE USR, TSR, QSR
404  !                -------------------------------------------
405  !
406  ZHWAVE(J) = 0.018*ZVMOD(J)*ZVMOD(J)*(1.+0.015*ZVMOD(J))
407  ZTWAVE(J) = 0.729*ZVMOD(J)
408  ZCWAVE(J) = XG*ZTWAVE(J)/(2.*XPI)
409  ZLWAVE(J) = ZTWAVE(J)*ZCWAVE(J)
410  !
411ENDDO
412!
413
414!
415DO JLOOP=1,MAXVAL(ITERMAX) ! begin of iterative loop
416  !
417  DO J=1,SIZE(PTA)
418    !
419    IF (JLOOP>ITERMAX(J)) CYCLE
420    !
421    IF (NGRVWAVES==0) THEN
422      ZO(J) = ZCHARN(J)*ZUSR(J)*ZUSR(J)/XG + 0.11*ZVISA(J)/ZUSR(J) !Smith 1988
423    ELSE IF (NGRVWAVES==1) THEN
424      ZO(J) = (50./(2.*XPI))*ZLWAVE(J)*(ZUSR(J)/ZCWAVE(J))**4.5 &
425              + 0.11*ZVISA(J)/ZUSR(J)                       !Oost et al. 2002
426    ELSE IF (NGRVWAVES==2) THEN
427      ZO(J) = 1200.*ZHWAVE(J)*(ZHWAVE(J)/ZLWAVE(J))**4.5 &
428              + 0.11*ZVISA(J)/ZUSR(J)                       !Taulor and Yelland 2001
429    ENDIF
430    !
431    ZRR(J) = ZO(J)*ZUSR(J)/ZVISA(J)
432    ZOQ(J) = MIN(1.15E-4 , 5.5E-5/ZRR(J)**0.6)
433    ZOT(J) = ZOQ(J)
434    !
435    ZZL(J) = XKARMAN * XG * PUREF(J) * &
436              ( ZTSR(J)*(1.+ZRVSRDM1*PQA(J)) + ZRVSRDM1*ZTA(J)*ZQSR(J) ) / &
437              ( ZTA(J)*ZUSR(J)*ZUSR(J)*(1.+ZRVSRDM1*PQA(J)) )
438    ZZTL(J)= ZZL(J)*PZREF(J)/PUREF(J)  ! for T
439!    ZZQL(J)=ZZL(J)*PZREF(J)/PUREF(J)  ! for Q
440  ENDDO
441  !
442  ZPUZ(:) = PSIFCTU(ZZL(1))
443  ZPTZ(:) = PSIFCTT(ZZTL(1))
444  !
445  DO J=1,SIZE(PTA)
446    !
447    ! ZPQZ(J)=PSIFCTT(ZZQL(J))
448    ZPQZ(J) = ZPTZ(J)
449    !
450    !             3.1 scale parameters
451    !
452    ZUSR(J) = ZDUWG(J)*XKARMAN/(LOG(PUREF(J)/ZO(J)) -ZPUZ(J))
453    ZTSR(J) = -ZDT(J)  *XKARMAN/(LOG(PZREF(J)/ZOT(J))-ZPTZ(J))
454    ZQSR(J) = -ZDQ(J)  *XKARMAN/(LOG(PZREF(J)/ZOQ(J))-ZPQZ(J))
455    !
456    !             3.2 Gustiness factor (ZWG)
457    !
458    IF(LPWG) THEN
459      ZBF(J) = -XG/ZTA(J)*ZUSR(J)*(ZTSR(J)+ZRVSRDM1*ZTA(J)*ZQSR(J))
460      IF (ZBF(J)>0.) THEN
461        ZWG(J) = ZBETAGUST*(ZBF(J)*ZZBL)**(1./3.)
462      ELSE
463        ZWG(J) = 0.2
464      ENDIF
465    ENDIF
466    ZDUWG(J) = SQRT(ZVMOD(J)**2 + ZWG(J)**2)
467    !
468  ENDDO
469  !
470ENDDO
471!-------------------------------------------------------------------------------
472!
473!            4.  COMPUTE transfer coefficients PCD, PCH, ZCE and SURFACE FLUXES
474!                --------------------------------------------------------------
475!
476ZTAU(:) = XUNDEF
477ZHF(:)  = XUNDEF
478ZEF(:)  = XUNDEF
479!
480ZWBAR(:) = 0.
481ZTAUR(:) = 0.
482ZRF(:)   = 0.
483!
484DO J=1,SIZE(PTA)
485  !
486  !
487  !            4. transfert coefficients PCD, PCH and PCE
488  !                 and neutral PCDN, ZCHN, ZCEN
489  !
490  PCD(J) = (ZUSR(J)/ZDUWG(J))**2.
491  PCH(J) = ZUSR(J)*ZTSR(J)/(ZDUWG(J)*(ZTA(J)*PEXNS(J)/PEXNA(J)-PSST(J)))
492  PCE(J) = ZUSR(J)*ZQSR(J)/(ZDUWG(J)*(PQA(J)-PQSAT(J)))
493  !
494  PCDN(J) = (XKARMAN/LOG(ZS/ZO(J)))**2.
495  ZCHN(J) = (XKARMAN/LOG(ZS/ZO(J)))*(XKARMAN/LOG(ZS/ZOT(J)))
496  ZCEN(J) = (XKARMAN/LOG(ZS/ZO(J)))*(XKARMAN/LOG(ZS/ZOQ(J)))
497  !
498  ZLV(J) = XLVTT + (XCPV-XCL)*(PSST(J)-XTT)
499  !
500  !            4. 2 surface fluxes
501  !
502  IF (ABS(PCDN(J))>1.E-2) THEN   !!!! secure COARE3.0 CODE
503    write(*,*) 'pb PCDN in COARE30: ',PCDN(J)
504    write(*,*) 'point: ',J,"/",SIZE(PTA)
505    write(*,*) 'roughness: ', ZO(J)
506    write(*,*) 'ustar: ',ZUSR(J)
507    write(*,*) 'wind: ',ZDUWG(J)
508    CALL abort_physic('COARE30',': PCDN too large -> no convergence',1)
509  ELSE
510    ZTSR(J) = -ZTSR(J)
511    ZQSR(J) = -ZQSR(J)
512    ZTAU(J) = -PRHOA(J)*ZUSR(J)*ZUSR(J)*ZVMOD(J)/ZDUWG(J)
513    ZHF(J)  =  PRHOA(J)*XCPD*ZUSR(J)*ZTSR(J)
514    ZEF(J)  =  PRHOA(J)*ZLV(J)*ZUSR(J)*ZQSR(J)
515    !
516    !           4.3 Contributions to surface  fluxes due to rainfall
517    !
518    ! SB: a priori, le facteur ZRDSRV=XRD/XRV est introduit pour
519    !     adapter la formule de Clausius-Clapeyron (pour l'air
520    !     sec) au cas humide.
521    IF (LPRECIP) THEN
522      !
523      ! heat surface  fluxes
524      !
525      ZTAC(J)  = ZTA(J)-XTT
526      !
527      ZXLR(J)  = XLVTT + (XCPV-XCL)* ZTAC(J)                            ! latent heat of rain vaporization
528      ZDQSDT(J)= ZQASAT(J) * ZXLR(J) / (XRD*ZTA(J)**2)                  ! Clausius-Clapeyron relation
529      ZDTMP(J) = (1.0 + 3.309e-3*ZTAC(J) -1.44e-6*ZTAC(J)*ZTAC(J)) * &  !heat diffusivity
530                  0.02411 / (PRHOA(J)*XCPD)
531      !
532      ZDWAT(J) = 2.11e-5 * (XP00/ZPA(J)) * (ZTA(J)/XTT)**1.94           ! water vapour diffusivity from eq (13.3)
533      !                                                                 ! of Pruppacher and Klett (1978)
534      ZALFAC(J)= 1.0 / (1.0 + &                                         ! Eq.11 in GoF95
535                   ZRDSRV*ZDQSDT(J)*ZXLR(J)*ZDWAT(J)/(ZDTMP(J)*XCPD))   ! ZALFAC=wet-bulb factor (sans dim)
536      ZCPLW(J) = 4224.8482 + ZTAC(J) * &
537                              ( -4.707 + ZTAC(J) * &
538                                (0.08499 + ZTAC(J) * &
539                                  (1.2826e-3 + ZTAC(J) * &
540                                    (4.7884e-5 - 2.0027e-6* ZTAC(J))))) ! specific heat
541      !
542      ZRF(J)   = PRAIN(J) * ZCPLW(J) * ZALFAC(J) * &                    !Eq.12 in GoF95 !SIGNE?
543                   (PSST(J) - ZTA(J) + (PQSAT(J)-PQA(J))*ZXLR(J)/XCPD )
544      !
545      ! Momentum flux due to rainfall
546      !
547      ZTAUR(J)=-0.85*(PRAIN(J) *ZVMOD(J)) !pp3752 in FBR96
548      !
549    ENDIF
550    !
551    !             4.4   Webb correction to latent heat flux
552    !
553    ZWBAR(J)=- (1./ZRDSRV)*ZUSR(J)*ZQSR(J) / (1.0+(1./ZRDSRV)*PQA(J)) &
554               - ZUSR(J)*ZTSR(J)/ZTA(J)                        ! Eq.21*rhoa in FBR96
555    !
556    !             4.5   friction velocity which contains correction du to rain
557    !
558    ZUSTAR2(J)= - (ZTAU(J) + ZTAUR(J)) / PRHOA(J)
559    PUSTAR(J) =  SQRT(ZUSTAR2(J))
560    !
561    !             4.6   Total surface fluxes
562    !
563    PSFTH (J) =  ZHF(J) + ZRF(J)
564    PSFTQ (J) =  ZEF(J) / ZLV(J)
565    !
566  ENDIF
567ENDDO
568
569
570coeffs = [PCD,&
571       PCE,&
572       PCH]
573
574!-------------------------------------------------------------------------------
575!
576!       5.  FINAL STEP : TOTAL SURFACE FLUXES AND DERIVED DIAGNOSTICS
577!           -----------
578!       5.1    Richardson number
579!
580!
581!------------STOP LA --------------------
582!ZDIRCOSZW(:) = 1.
583! CALL SURFACE_RI(PSST,PQSAT,PEXNS,PEXNA,ZTA,ZQASAT,&
584!                PZREF,PUREF,ZDIRCOSZW,PVMOD,PRI   )
585!!
586!!       5.2     Aerodynamical conductance and resistance
587!!
588!ZAC(:) = PCH(:)*ZVMOD(:)
589!PRESA(:) = 1. / MAX(ZAC(:),XSURF_EPSILON)
590!
591!!       5.3 Z0 and Z0H over sea
592!!
593!PZ0SEA(:) =  ZCHARN(:) * ZUSTAR2(:) / XG + XVZ0CM * PCD(:) / PCDN(:)
594!!
595!!PZ0HSEA(:) = PZ0SEA(:)
596!!
597!IF (LHOOK) CALL DR_HOOK('COARE30_FLUX',1,ZHOOK_HANDLE)
598!
599!-------------------------------------------------------------------------------
600!
601END SUBROUTINE COARE30_FLUX_CNRM
602
603end module coare30_flux_cnrm_mod
Note: See TracBrowser for help on using the repository browser.