source: LMDZ6/branches/Amaury_dev/libf/phylmd/coare30_flux_cnrm_mod.F90 @ 5135

Last change on this file since 5135 was 5119, checked in by abarral, 16 months ago

enforce PRIVATE by default in several modules, expose PUBLIC as needed
move eigen.f90 to obsolete/
(lint) aslong the way
File size: 19.9 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!----------Rajout Olive ---------
75USE dimphy
76USE indice_sol_mod
77USE coare_cp_mod, ONLY: PSIFCTT=>psit_30, PSIFCTU=>psiuo
78
79!--------------------------------
80
81USE MODD_CSTS,       ONLY: XKARMAN, XG, XSTEFAN, XRD, XRV, XPI, &
82                            XLVTT, XCL, XCPD, XCPV, XRHOLW, XTT, &
83                            XP00
84USE lmdz_abort_physic, ONLY: abort_physic
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
103
104IMPLICIT NONE
105
106!*      0.1    declarations of arguments
107
108
109
110!TYPE(SEAFLUX_t), INTENT(INOUT) :: S
111
112REAL, DIMENSION(klon), INTENT(IN)       :: PTA   ! air temperature at atm. level (K)
113REAL, DIMENSION(klon), INTENT(IN)       :: PQA   ! air humidity at atm. level (kg/kg)
114!REAL, DIMENSION(:), INTENT(IN)       :: PEXNA ! Exner function at atm. level
115!REAL, DIMENSION(:), INTENT(IN)       :: PRHOA ! air density at atm. level
116REAL, DIMENSION(klon), INTENT(IN)       :: PVMOD ! module of wind at atm. wind level (m/s)
117REAL, DIMENSION(klon), INTENT(IN)       :: PZREF ! atm. level for temp. and humidity (m)
118REAL, DIMENSION(klon), INTENT(IN)       :: PUREF ! atm. level for wind (m)
119REAL, DIMENSION(klon), INTENT(IN)       :: PSST  ! Sea Surface Temperature (K)
120!REAL, DIMENSION(:), INTENT(IN)       :: PEXNS ! Exner function at sea surface
121REAL, DIMENSION(klon), INTENT(IN)       :: PPS   ! air pressure at sea surface (Pa)
122REAL, DIMENSION(klon), INTENT(IN)       :: PRAIN !precipitation rate (kg/s/m2)
123REAL, DIMENSION(klon), INTENT(IN)       :: PPA        ! air pressure at atm level (Pa)
124REAL, DIMENSION(klon), INTENT(IN)       :: PQSATA        ! air pressure at atm level (Pa)
125
126REAL, DIMENSION(klon), INTENT(INOUT)    :: PZ0SEA! roughness length over the ocean
127
128!  surface fluxes : latent heat, sensible heat, friction fluxes
129REAL, DIMENSION(klon), INTENT(OUT)      :: PSFTH ! heat flux (W/m2)
130REAL, DIMENSION(klon), INTENT(OUT)      :: PSFTQ ! water flux (kg/m2/s)
131REAL, DIMENSION(klon), INTENT(OUT)      :: PUSTAR! friction velocity (m/s)
132
133! diagnostics
134REAL, DIMENSION(klon), INTENT(OUT)      :: PQSAT ! humidity at saturation
135REAL, DIMENSION(klon), INTENT(OUT)      :: PCD   ! heat drag coefficient
136REAL, DIMENSION(klon), INTENT(OUT)      :: PCDN  ! momentum drag coefficient
137REAL, DIMENSION(klon), INTENT(OUT)      :: PCH   ! neutral momentum drag coefficient
138REAL, DIMENSION(klon), INTENT(OUT)      :: PCE  !transfer coef. for latent heat flux
139REAL, DIMENSION(klon), INTENT(OUT)      :: PRI   ! Richardson number
140REAL, DIMENSION(klon), INTENT(OUT)      :: PRESA ! aerodynamical resistance
141REAL, DIMENSION(klon), INTENT(OUT)      :: PZ0HSEA ! heat roughness length
142
143LOGICAL,            INTENT(IN)    :: LPRECIP   !
144LOGICAL,            INTENT(IN)    :: LPWG    !
145REAL, DIMENSION(3), INTENT(INOUT)      :: coeffs
146
147
148!*      0.2    declarations of local variables
149
150REAL, DIMENSION(SIZE(PTA))      :: ZVMOD    ! wind intensity
151REAL, DIMENSION(SIZE(PTA))      :: ZPA      ! Pressure at atm. level
152REAL, DIMENSION(SIZE(PTA))      :: ZTA      ! Temperature at atm. level
153REAL, DIMENSION(SIZE(PTA))      :: ZQASAT   ! specific humidity at saturation  at atm. level (kg/kg)
154
155!rajout
156REAL, DIMENSION(SIZE(PTA))       :: PEXNA      ! Exner function at atm level
157REAL, DIMENSION(SIZE(PTA))       :: PEXNS      ! Exner function at atm level
158
159REAL, DIMENSION(SIZE(PTA))      :: ZO       ! rougness length ref
160REAL, DIMENSION(SIZE(PTA))      :: ZWG      ! gustiness factor (m/s)
161
162REAL, DIMENSION(SIZE(PTA))      :: ZDU,ZDT,ZDQ,ZDUWG !differences
163
164REAL, DIMENSION(SIZE(PTA))      :: ZUSR        !velocity scaling parameter "ustar" (m/s) = friction velocity
165REAL, DIMENSION(SIZE(PTA))      :: ZTSR        !temperature sacling parameter "tstar" (degC)
166REAL, DIMENSION(SIZE(PTA))      :: ZQSR        !humidity scaling parameter "qstar" (kg/kg)
167
168REAL, DIMENSION(SIZE(PTA))      :: ZU10,ZT10   !vertical profils (10-m height)
169REAL, DIMENSION(SIZE(PTA))      :: ZVISA       !kinematic viscosity of dry air
170REAL, DIMENSION(SIZE(PTA))      :: ZO10,ZOT10  !roughness length at 10m
171REAL, DIMENSION(SIZE(PTA))      :: ZCD,ZCT,ZCC
172REAL, DIMENSION(SIZE(PTA))      :: ZCD10,ZCT10 !transfer coef. at 10m
173REAL, DIMENSION(SIZE(PTA))      :: ZRIBU,ZRIBCU
174REAL, DIMENSION(SIZE(PTA))      :: ZETU,ZL10
175
176REAL, DIMENSION(SIZE(PTA))      :: ZCHARN                      !Charnock number depends on wind module
177REAL, DIMENSION(SIZE(PTA))      :: ZTWAVE,ZHWAVE,ZCWAVE,ZLWAVE !to compute gravity waves' impact
178
179REAL, DIMENSION(SIZE(PTA))      :: ZZL,ZZTL!,ZZQL    !Obukhovs stability
180                                                     !param. z/l for u,T,q
181REAL, DIMENSION(SIZE(PTA))      :: ZRR
182REAL, DIMENSION(SIZE(PTA))      :: ZOT,ZOQ           !rougness length ref
183REAL, DIMENSION(SIZE(PTA))      :: ZPUZ,ZPTZ,ZPQZ    !PHI funct. for u,T,q
184
185REAL, DIMENSION(SIZE(PTA))      :: ZBF               !constants to compute gustiness factor
186
187REAL, DIMENSION(SIZE(PTA))      :: ZTAU       !momentum flux (W/m2)
188REAL, DIMENSION(SIZE(PTA))      :: ZHF        !sensible heat flux (W/m2)
189REAL, DIMENSION(SIZE(PTA))      :: ZEF        !latent heat flux (W/m2)
190REAL, DIMENSION(SIZE(PTA))      :: ZWBAR      !diag for webb correction but not used here after
191REAL, DIMENSION(SIZE(PTA))      :: ZTAUR      !momentum flux due to rain (W/m2)
192REAL, DIMENSION(SIZE(PTA))      :: ZRF        !sensible heat flux due to rain (W/m2)
193REAL, DIMENSION(SIZE(PTA))      :: ZCHN,ZCEN  !neutral coef. for heat and vapor
194
195REAL, DIMENSION(SIZE(PTA))      :: ZLV      !latent heat constant
196
197REAL, DIMENSION(SIZE(PTA))      :: ZTAC,ZDQSDT,ZDTMP,ZDWAT,ZALFAC ! for precipitation impact
198REAL, DIMENSION(SIZE(PTA))      :: ZXLR                           ! vaporisation  heat  at a given temperature
199REAL, DIMENSION(SIZE(PTA))      :: ZCPLW                          ! specific heat for water at a given temperature
200
201REAL, DIMENSION(SIZE(PTA))      :: ZUSTAR2  ! square of friction velocity
202
203REAL, DIMENSION(SIZE(PTA))      :: ZDIRCOSZW! orography slope cosine (=1 on water!)
204REAL, DIMENSION(SIZE(PTA))      :: ZAC      ! Aerodynamical conductance
205
206
207INTEGER, DIMENSION(SIZE(PTA))   :: ITERMAX             ! maximum number of iterations
208
209REAL    :: ZRVSRDM1,ZRDSRV,ZR2 ! thermodynamic constants
210REAL    :: ZBETAGUST           !gustiness factor
211REAL    :: ZZBL                !atm. boundary layer depth (m)
212REAL    :: ZVISW               !m2/s kinematic viscosity of water
213REAL    :: ZS                  !height of rougness length ref
214REAL    :: ZCH10               !transfer coef. at 10m
215
216REAL    :: QSAT_SEAWATER
217REAL    :: QSATSEAW_1D
218
219INTEGER :: J, JLOOP    !loop indice
220!REAL(KIND=JPRB) :: ZHOOK_HANDLE
221
222!--------- Modif Olive -----------------
223REAL, DIMENSION(SIZE(PTA))        :: PRHOA
224REAL, PARAMETER                   :: XUNDEF = 1.E+20
225
226REAL       :: XVCHRNK = 0.021
227REAL       :: XVZ0CM = 1.0E-5
228!REAL       :: XRIMAX
229
230
231INTEGER :: PREF             ! reference pressure for exner function
232INTEGER :: NGRVWAVES        ! Pour le choix du z0
233
234INCLUDE "YOMCST.h"
235INCLUDE "clesphys.h"
236
237!--------------------------------------
238
239
240PRHOA(:) = PPS(:) / (287.1 * PTA(:) * (1.+.61*PQA(:)))
241
242PREF = 100000.                     ! = 1000 hPa
243NGRVWAVES = 1
244
245PEXNA = (PPA/PREF)**(RD/RCPD)
246PEXNS = (PPS/PREF)**(RD/RCPD)
247
248!-------------------------------------------------------------------------------
249
250!       1.     Initializations
251!              ---------------
252
253!       1.1   Constants and parameters
254
255!IF (LHOOK) CALL DR_HOOK('COARE30_FLUX',0,ZHOOK_HANDLE)
256
257ZRVSRDM1  = XRV/XRD-1. ! 0.607766
258ZRDSRV    = XRD/XRV    ! 0.62198
259ZR2       = 1.-ZRDSRV  ! pas utilisé dans cette routine
260ZBETAGUST = 1.2        ! value based on TOGA-COARE experiment
261ZZBL      = 600.       ! Set a default value for boundary layer depth
262ZS        = 10.        ! Standard heigth =10m
263ZCH10     = 0.00115
264
265ZVISW     = 1.E-6
266
267!       1.2   Array initialization by undefined values
268
269PSFTH (:)=XUNDEF
270PSFTQ (:)=XUNDEF
271PUSTAR(:)=XUNDEF
272
273PCD(:) = XUNDEF
274PCDN(:) = XUNDEF
275PCH(:) = XUNDEF
276PCE(:) =XUNDEF
277PRI(:) = XUNDEF
278
279PRESA(:)=XUNDEF
280
281!-------------------------------------------------------------------------------
282!       2. INITIAL GUESS FOR THE ITERATIVE METHOD
283!          -------------------------------------
284
285!       2.0     Temperature
286
287! Set a non-zero value for the temperature gradient
288
289WHERE((PTA(:)*PEXNS(:)/PEXNA(:)-PSST(:))==0.)
290      ZTA(:)=PTA(:)-1E-3
291ELSEWHERE
292      ZTA(:)=PTA(:)
293ENDWHERE
294
295!       2.1     Wind and humidity
296
297! Sea surface specific humidity
298
299!PQSAT(:)=QSAT_SEAWATER(PSST(:),PPS(:))
300PQSAT(:)=QSATSEAW_1D(PSST(:),PPS(:))
301
302! Set a minimum value to wind
303
304!ZVMOD(:) = WIND_THRESHOLD(PVMOD(:),PUREF(:))
305
306
307ZVMOD = MAX(PVMOD , 0.1 * MIN(10.,PUREF) )    !set a minimum value to wind
308!ZVMOD = PVMOD    !set a minimum value to wind
309
310! Specific humidity at saturation at the atm. level
311
312ZPA(:) = XP00* (PEXNA(:)**(XCPD/XRD))
313!ZQASAT(:) = QSAT_SEAWATER(ZTA(:),ZPA(:))
314ZQASAT = QSATSEAW_1D(ZTA(:),ZPA(:))
315
316
317ZO(:)  = 0.0001
318ZWG(:) = 0.
319IF (LPWG) ZWG(:) = 0.5
320
321ZCHARN(:) = 0.011
322
323DO J=1,SIZE(PTA)
324
325  !      2.2       initial guess
326
327  ZDU(J) = ZVMOD(J)   !wind speed difference with surface current(=0) (m/s)
328                      !initial guess for gustiness factor
329  ZDT(J) = -(ZTA(J)/PEXNA(J)) + (PSST(J)/PEXNS(J)) !potential temperature difference
330  ZDQ(J) = PQSAT(J)-PQA(J)                         !specific humidity difference
331
332  ZDUWG(J) = SQRT(ZDU(J)**2+ZWG(J)**2)     !wind speed difference including gustiness ZWG
333
334  !      2.3   initialization of neutral coefficients
335
336  ZU10(J)  = ZDUWG(J)*LOG(ZS/ZO(J))/LOG(PUREF(J)/ZO(J))
337  ZUSR(J)  = 0.035*ZU10(J)
338  ZVISA(J) = 1.326E-5*(1.+6.542E-3*(ZTA(J)-XTT)+&
339             8.301E-6*(ZTA(J)-XTT)**2-4.84E-9*(ZTA(J)-XTT)**3) !Andrea (1989) CRREL Rep. 89-11
340
341  ZO10(J) = ZCHARN(J)*ZUSR(J)*ZUSR(J)/XG+0.11*ZVISA(J)/ZUSR(J)
342  ZCD(J)  = (XKARMAN/LOG(PUREF(J)/ZO10(J)))**2  !drag coefficient
343  ZCD10(J)= (XKARMAN/LOG(ZS/ZO10(J)))**2
344  ZCT10(J)= ZCH10/SQRT(ZCD10(J))
345  ZOT10(J)= ZS/EXP(XKARMAN/ZCT10(J))
346
347  !-------------------------------------------------------------------------------
348  !             Grachev and Fairall (JAM, 1997)
349  ZCT(J) = XKARMAN/LOG(PZREF(J)/ZOT10(J))      !temperature transfer coefficient
350  ZCC(J) = XKARMAN*ZCT(J)/ZCD(J)               !z/L vs Rib linear coef.
351
352  ZRIBCU(J) = -PUREF(J)/(ZZBL*0.004*ZBETAGUST**3) !saturation or plateau Rib
353  !ZRIBU(J) =-XG*PUREF(J)*(ZDT(J)+ZRVSRDM1*(ZTA(J)-XTT)*ZDQ)/&
354  !     &((ZTA(J)-XTT)*ZDUWG(J)**2)
355  ZRIBU(J)  = -XG*PUREF(J)*(ZDT(J)+ZRVSRDM1*ZTA(J)*ZDQ(J))/&
356               (ZTA(J)*ZDUWG(J)**2)
357
358  IF (ZRIBU(J)<0.) THEN
359    ZETU(J) = ZCC(J)*ZRIBU(J)/(1.+ZRIBU(J)/ZRIBCU(J))    !Unstable G and F
360  ELSE
361    ZETU(J) = ZCC(J)*ZRIBU(J)/(1.+27./9.*ZRIBU(J)/ZCC(J))!Stable
362  ENDIF
363
364  ZL10(J) = PUREF(J)/ZETU(J) !MO length
365
366ENDDO
367
368!  First guess M-O stability dependent scaling params. (u*,T*,q*) to estimate ZO and z/L (ZZL)
369ZUSR(:) = ZDUWG(:)*XKARMAN/(LOG(PUREF(:)/ZO10(:))-PSIFCTU(PUREF(1)/ZL10(1)))
370ZTSR(:) = -ZDT(:)*XKARMAN/(LOG(PZREF(:)/ZOT10(:))-PSIFCTT(PZREF(1)/ZL10(1)))
371ZQSR(:) = -ZDQ(:)*XKARMAN/(LOG(PZREF(:)/ZOT10(:))-PSIFCTT(PZREF(1)/ZL10(1)))
372
373ZZL(:) = 0.0
374
375DO J=1,SIZE(PTA)
376
377  IF (ZETU(J)>50.) THEN
378    ITERMAX(J) = 1
379  ELSE
380    ITERMAX(J) = 3 !number of iterations
381  ENDIF
382
383  !then modify Charnork for high wind speeds Chris Fairall's data
384  IF (ZDUWG(J)>10.) ZCHARN(J) = 0.011 + (0.018-0.011)*(ZDUWG(J)-10.)/(18.-10.)
385  IF (ZDUWG(J)>18.) ZCHARN(J) = 0.018
386
387  !                3.  ITERATIVE LOOP TO COMPUTE USR, TSR, QSR
388  !                -------------------------------------------
389
390  ZHWAVE(J) = 0.018*ZVMOD(J)*ZVMOD(J)*(1.+0.015*ZVMOD(J))
391  ZTWAVE(J) = 0.729*ZVMOD(J)
392  ZCWAVE(J) = XG*ZTWAVE(J)/(2.*XPI)
393  ZLWAVE(J) = ZTWAVE(J)*ZCWAVE(J)
394
395ENDDO
396
397DO JLOOP=1,MAXVAL(ITERMAX) ! begin of iterative loop
398
399  DO J=1,SIZE(PTA)
400
401    IF (JLOOP>ITERMAX(J)) CYCLE
402
403    IF (NGRVWAVES==0) THEN
404      ZO(J) = ZCHARN(J)*ZUSR(J)*ZUSR(J)/XG + 0.11*ZVISA(J)/ZUSR(J) !Smith 1988
405    ELSE IF (NGRVWAVES==1) THEN
406      ZO(J) = (50./(2.*XPI))*ZLWAVE(J)*(ZUSR(J)/ZCWAVE(J))**4.5 &
407              + 0.11*ZVISA(J)/ZUSR(J)                       !Oost et al. 2002
408    ELSE IF (NGRVWAVES==2) THEN
409      ZO(J) = 1200.*ZHWAVE(J)*(ZHWAVE(J)/ZLWAVE(J))**4.5 &
410              + 0.11*ZVISA(J)/ZUSR(J)                       !Taulor and Yelland 2001
411    ENDIF
412
413    ZRR(J) = ZO(J)*ZUSR(J)/ZVISA(J)
414    ZOQ(J) = MIN(1.15E-4 , 5.5E-5/ZRR(J)**0.6)
415    ZOT(J) = ZOQ(J)
416
417    ZZL(J) = XKARMAN * XG * PUREF(J) * &
418              ( ZTSR(J)*(1.+ZRVSRDM1*PQA(J)) + ZRVSRDM1*ZTA(J)*ZQSR(J) ) / &
419              ( ZTA(J)*ZUSR(J)*ZUSR(J)*(1.+ZRVSRDM1*PQA(J)) )
420    ZZTL(J)= ZZL(J)*PZREF(J)/PUREF(J)  ! for T
421!    ZZQL(J)=ZZL(J)*PZREF(J)/PUREF(J)  ! for Q
422  ENDDO
423
424  ZPUZ(:) = PSIFCTU(ZZL(1))
425  ZPTZ(:) = PSIFCTT(ZZTL(1))
426
427  DO J=1,SIZE(PTA)
428
429    ! ZPQZ(J)=PSIFCTT(ZZQL(J))
430    ZPQZ(J) = ZPTZ(J)
431
432    !             3.1 scale parameters
433
434    ZUSR(J) = ZDUWG(J)*XKARMAN/(LOG(PUREF(J)/ZO(J)) -ZPUZ(J))
435    ZTSR(J) = -ZDT(J)  *XKARMAN/(LOG(PZREF(J)/ZOT(J))-ZPTZ(J))
436    ZQSR(J) = -ZDQ(J)  *XKARMAN/(LOG(PZREF(J)/ZOQ(J))-ZPQZ(J))
437
438    !             3.2 Gustiness factor (ZWG)
439
440    IF(LPWG) THEN
441      ZBF(J) = -XG/ZTA(J)*ZUSR(J)*(ZTSR(J)+ZRVSRDM1*ZTA(J)*ZQSR(J))
442      IF (ZBF(J)>0.) THEN
443        ZWG(J) = ZBETAGUST*(ZBF(J)*ZZBL)**(1./3.)
444      ELSE
445        ZWG(J) = 0.2
446      ENDIF
447    ENDIF
448    ZDUWG(J) = SQRT(ZVMOD(J)**2 + ZWG(J)**2)
449
450  ENDDO
451
452ENDDO
453!-------------------------------------------------------------------------------
454
455!            4.  COMPUTE transfer coefficients PCD, PCH, ZCE and SURFACE FLUXES
456!                --------------------------------------------------------------
457
458ZTAU(:) = XUNDEF
459ZHF(:)  = XUNDEF
460ZEF(:)  = XUNDEF
461
462ZWBAR(:) = 0.
463ZTAUR(:) = 0.
464ZRF(:)   = 0.
465
466DO J=1,SIZE(PTA)
467
468
469  !            4. transfert coefficients PCD, PCH and PCE
470  !                 and neutral PCDN, ZCHN, ZCEN
471
472  PCD(J) = (ZUSR(J)/ZDUWG(J))**2.
473  PCH(J) = ZUSR(J)*ZTSR(J)/(ZDUWG(J)*(ZTA(J)*PEXNS(J)/PEXNA(J)-PSST(J)))
474  PCE(J) = ZUSR(J)*ZQSR(J)/(ZDUWG(J)*(PQA(J)-PQSAT(J)))
475
476  PCDN(J) = (XKARMAN/LOG(ZS/ZO(J)))**2.
477  ZCHN(J) = (XKARMAN/LOG(ZS/ZO(J)))*(XKARMAN/LOG(ZS/ZOT(J)))
478  ZCEN(J) = (XKARMAN/LOG(ZS/ZO(J)))*(XKARMAN/LOG(ZS/ZOQ(J)))
479
480  ZLV(J) = XLVTT + (XCPV-XCL)*(PSST(J)-XTT)
481
482  !            4. 2 surface fluxes
483
484  IF (ABS(PCDN(J))>1.E-2) THEN   !!!! secure COARE3.0 CODE
485    WRITE(*,*) 'pb PCDN in COARE30: ',PCDN(J)
486    WRITE(*,*) 'point: ',J,"/",SIZE(PTA)
487    WRITE(*,*) 'roughness: ', ZO(J)
488    WRITE(*,*) 'ustar: ',ZUSR(J)
489    WRITE(*,*) 'wind: ',ZDUWG(J)
490    CALL abort_physic('COARE30',': PCDN too large -> no convergence',1)
491  ELSE
492    ZTSR(J) = -ZTSR(J)
493    ZQSR(J) = -ZQSR(J)
494    ZTAU(J) = -PRHOA(J)*ZUSR(J)*ZUSR(J)*ZVMOD(J)/ZDUWG(J)
495    ZHF(J)  =  PRHOA(J)*XCPD*ZUSR(J)*ZTSR(J)
496    ZEF(J)  =  PRHOA(J)*ZLV(J)*ZUSR(J)*ZQSR(J)
497
498    !           4.3 Contributions to surface  fluxes due to rainfall
499
500    ! SB: a priori, le facteur ZRDSRV=XRD/XRV est introduit pour
501    !     adapter la formule de Clausius-Clapeyron (pour l'air
502    !     sec) au cas humide.
503    IF (LPRECIP) THEN
504
505      ! heat surface  fluxes
506
507      ZTAC(J)  = ZTA(J)-XTT
508
509      ZXLR(J)  = XLVTT + (XCPV-XCL)* ZTAC(J)                            ! latent heat of rain vaporization
510      ZDQSDT(J)= ZQASAT(J) * ZXLR(J) / (XRD*ZTA(J)**2)                  ! Clausius-Clapeyron relation
511      ZDTMP(J) = (1.0 + 3.309e-3*ZTAC(J) -1.44e-6*ZTAC(J)*ZTAC(J)) * &  !heat diffusivity
512                  0.02411 / (PRHOA(J)*XCPD)
513
514      ZDWAT(J) = 2.11e-5 * (XP00/ZPA(J)) * (ZTA(J)/XTT)**1.94           ! water vapour diffusivity from eq (13.3)
515      !                                                                 ! of Pruppacher and Klett (1978)
516      ZALFAC(J)= 1.0 / (1.0 + &                                         ! Eq.11 in GoF95
517                   ZRDSRV*ZDQSDT(J)*ZXLR(J)*ZDWAT(J)/(ZDTMP(J)*XCPD))   ! ZALFAC=wet-bulb factor (sans dim)
518      ZCPLW(J) = 4224.8482 + ZTAC(J) * &
519                              ( -4.707 + ZTAC(J) * &
520                                (0.08499 + ZTAC(J) * &
521                                  (1.2826e-3 + ZTAC(J) * &
522                                    (4.7884e-5 - 2.0027e-6* ZTAC(J))))) ! specific heat
523
524      ZRF(J)   = PRAIN(J) * ZCPLW(J) * ZALFAC(J) * &                    !Eq.12 in GoF95 !SIGNE?
525                   (PSST(J) - ZTA(J) + (PQSAT(J)-PQA(J))*ZXLR(J)/XCPD )
526
527      ! Momentum flux due to rainfall
528
529      ZTAUR(J)=-0.85*(PRAIN(J) *ZVMOD(J)) !pp3752 in FBR96
530
531    ENDIF
532
533    !             4.4   Webb correction to latent heat flux
534
535    ZWBAR(J)=- (1./ZRDSRV)*ZUSR(J)*ZQSR(J) / (1.0+(1./ZRDSRV)*PQA(J)) &
536               - ZUSR(J)*ZTSR(J)/ZTA(J)                        ! Eq.21*rhoa in FBR96
537
538    !             4.5   friction velocity which contains correction du to rain
539
540    ZUSTAR2(J)= - (ZTAU(J) + ZTAUR(J)) / PRHOA(J)
541    PUSTAR(J) =  SQRT(ZUSTAR2(J))
542
543    !             4.6   Total surface fluxes
544
545    PSFTH (J) =  ZHF(J) + ZRF(J)
546    PSFTQ (J) =  ZEF(J) / ZLV(J)
547
548  ENDIF
549ENDDO
550
551
552coeffs = [PCD,&
553       PCE,&
554       PCH]
555
556!-------------------------------------------------------------------------------
557
558!       5.  FINAL STEP : TOTAL SURFACE FLUXES AND DERIVED DIAGNOSTICS
559!           -----------
560!       5.1    Richardson number
561
562
563!------------STOP LA --------------------
564!ZDIRCOSZW(:) = 1.
565! CALL SURFACE_RI(PSST,PQSAT,PEXNS,PEXNA,ZTA,ZQASAT,&
566!                PZREF,PUREF,ZDIRCOSZW,PVMOD,PRI   )
567!!
568!!       5.2     Aerodynamical conductance and resistance
569!!
570!ZAC(:) = PCH(:)*ZVMOD(:)
571!PRESA(:) = 1. / MAX(ZAC(:),XSURF_EPSILON)
572
573!!       5.3 Z0 and Z0H over sea
574!!
575!PZ0SEA(:) =  ZCHARN(:) * ZUSTAR2(:) / XG + XVZ0CM * PCD(:) / PCDN(:)
576!!
577!!PZ0HSEA(:) = PZ0SEA(:)
578!!
579!IF (LHOOK) CALL DR_HOOK('COARE30_FLUX',1,ZHOOK_HANDLE)
580
581!-------------------------------------------------------------------------------
582
583END SUBROUTINE COARE30_FLUX_CNRM
584
585END MODULE coare30_flux_cnrm_mod
Note: See TracBrowser for help on using the repository browser.