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

Last change on this file since 5144 was 5144, checked in by abarral, 3 months ago

Put YOMCST.h into modules

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