source: LMDZ6/trunk/libf/phylmd/coare30_flux_cnrm.F90 @ 4912

Last change on this file since 4912 was 4722, checked in by Laurent Fairhead, 8 months ago

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