source: LMDZ6/branches/LMDZ_ECRad/libf/phylmd/coare30_flux_cnrm.F90 @ 5452

Last change on this file since 5452 was 4727, checked in by idelkadi, 15 months ago

Merged trunk changes -r4488:4726 LMDZ_ECRad branch
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.