source: LMDZ6/trunk/libf/phylmd/ecumev6_flux.f90 @ 5319

Last change on this file since 5319 was 5285, checked in by abarral, 10 days ago

As discussed internally, remove generic ONLY: ... for new _mod_h modules
File size: 28.3 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    SUBROUTINE ECUMEV6_FLUX(PZ0SEA,PTA,PSST,PQA,PQSAT,PVMOD, &
7                            PZREF,PSSS,PUREF,PPS,PPA,OPRECIP,OPWEBB,        &
8                            PSFTH,PSFTQ,PUSTAR,PCD,PCDN,PCH,PCE,        &
9                            PRI,PRESA,PRAIN,PZ0HSEA,OPERTFLUX,coeffs   )
10!###############################################################################
11!!
12!!****  *ECUMEV6_FLUX*
13!!
14!!    PURPOSE
15!!    -------
16!       Calculate the surface turbulent fluxes of heat, moisture, and momentum
17!       over sea surface + corrections due to rainfall & Webb effect.
18!!
19!!**  METHOD
20!!    ------
21!       The estimation of the transfer coefficients relies on the iterative
22!       computation of the scaling parameters U*/Teta*/q*. The convergence is
23!       supposed to be reached in NITERFL iterations maximum.
24!       Neutral transfer coefficients for momentum/temperature/humidity
25!       are computed as a function of the 10m-height neutral wind speed using
26!       the ECUME_V6 formulation based on the multi-campaign (POMME,FETCH,CATCH,
27!       SEMAPHORE,EQUALANT) ALBATROS dataset.
28!!
29!!    EXTERNAL
30!!    --------
31!!
32!!    IMPLICIT ARGUMENTS
33!!    ------------------
34!!
35!!    REFERENCE
36!!    ---------
37!!      Fairall et al (1996), JGR, 3747-3764
38!!      Gosnell et al (1995), JGR, 437-442
39!!      Fairall et al (1996), JGR, 1295-1308
40!!
41!!    AUTHOR
42!!    ------
43!!      C. Lebeaupin  *Météo-France* (adapted from S. Belamari's code)
44!!
45!!    MODIFICATIONS
46!!    -------------
47!!      Original     15/03/2005
48!!      Modified        01/2006  C. Lebeaupin (adapted from  A. Pirani's code)
49!!      Modified        08/2009  B. Decharme: limitation of Ri
50!!      Modified        09/2012  B. Decharme: CD correction
51!!      Modified        09/2012  B. Decharme: limitation of Ri in surface_ri.F90
52!!      Modified        10/2012  P. Le Moigne: extra inputs for FLake use
53!!      Modified        06/2013  B. Decharme: bug in z0 (output) computation
54!!      Modified        12/2013  S. Belamari: ZRF computation updated:
55!!                                1. ZP00/PPA in ZDWAT, ZLVA in ZDQSDT/ZBULB/ZRF
56!!                                2. ZDWAT/ZDTMP in ZBULB/ZRF (Gosnell et al 95)
57!!                                3. cool skin correction included
58!!      Modified        01/2014  S. Belamari: salinity impact on latent heat of
59!!                                vaporization of seawater included
60!!      Modified        01/2014  S. Belamari: new formulation for pure water
61!!                                specific heat (ZCPWA)
62!!      Modified        01/2014  S. Belamari: 4 choices for PZ0SEA computation
63!!      Modified        12/2015  S. Belamari: ECUME now provides parameterisations
64!!                                for:  U10n*sqrt(CDN)          instead of CDN
65!!                                      U10n*CHN/sqrt(CDN)         "       CHN
66!!                                      U10n*CEN/sqrt(CDN)         "       CEN
67!!      Modified        01/2016  S. Belamari: New ECUME formulation
68!!
69!!      To be done:
70!!      include gustiness computation following Mondon & Redelsperger (1998)
71!!!
72!-------------------------------------------------------------------------------
73!!
74!!    MODIFICATIONS RELATED TO SST CORRECTION COMPUTATION
75!!    ---------------------------------------------------
76!!      Modified        09/2013  S. Belamari: use 0.98 for the ocean emissivity
77!!                                following up to date satellite measurements in
78!!                                the 8-14 μm range (obtained values range from
79!!                                0.98 to 0.99).
80!!!
81!-------------------------------------------------------------------------------
82!
83!       0.   DECLARATIONS
84!            ------------
85!
86USE dimphy
87USE indice_sol_mod
88USE MODD_CSTS,             ONLY : XPI, XDAY, XKARMAN, XG, XP00, XSTEFAN, XRD, XRV,   &
89                                  XCPD, XCPV, XCL, XTT, XLVTT
90
91
92!USE MODD_SURF_PAR,         ONLY : XUNDEF
93!USE MODD_SURF_ATM,         ONLY : XVCHRNK, XVZ0CM
94!USE MODD_REPROD_OPER,      ONLY : CCHARNOCK
95!
96!USE MODE_THERMOS
97!USE MODI_WIND_THRESHOLD
98!USE MODI_SURFACE_RI
99!
100!USE YOMHOOK,   ONLY : LHOOK,   DR_HOOK
101!USE PARKIND1,  ONLY : JPRB
102!
103!USE MODI_ABOR1_SFX
104USE yomcst_mod_h
105USE clesphys_mod_h
106IMPLICIT NONE
107!
108!       0.1. Declarations of arguments
109!
110REAL, DIMENSION(klon), INTENT(IN)    :: PVMOD      ! module of wind at atm level (m/s)
111REAL, DIMENSION(klon), INTENT(IN)    :: PTA        ! air temperature at atm level (K)
112REAL, DIMENSION(klon), INTENT(IN)    :: PQA        ! air spec. hum. at atm level (kg/kg)
113REAL, DIMENSION(klon), INTENT(IN)    :: PQSAT      ! sea surface spec. hum. (kg/kg)
114REAL, DIMENSION(klon), INTENT(IN)    :: PPA        ! air pressure at atm level (Pa)
115!REAL, DIMENSION(:), INTENT(IN)    :: PRHOA      ! air density at atm level (kg/m3)
116!REAL, DIMENSION(:), INTENT(IN)    :: PEXNA      ! Exner function at atm level
117REAL, DIMENSION(klon), INTENT(IN)    :: PUREF      ! atm level for wind (m)
118REAL, DIMENSION(klon), INTENT(IN)    :: PZREF      ! atm level for temp./hum. (m)
119REAL, DIMENSION(klon), INTENT(IN)    :: PSSS       ! Sea Surface Salinity (g/kg)
120REAL, DIMENSION(klon), INTENT(IN)    :: PPS        ! air pressure at sea surface (Pa)
121!REAL, DIMENSION(:), INTENT(IN)    :: PEXNS      ! Exner function at sea surface
122!REAL, DIMENSION(:), INTENT(IN)    :: PPERTFLUX  ! stochastic flux perturbation pattern
123! for correction
124!REAL,               INTENT(IN)    :: PICHCE    !
125LOGICAL,            INTENT(IN)    :: OPRECIP   !
126LOGICAL,            INTENT(IN)    :: OPWEBB    !
127LOGICAL,            INTENT(IN)    :: OPERTFLUX
128REAL, DIMENSION(klon), INTENT(IN)    :: PRAIN     ! precipitation rate (kg/s/m2)
129!
130!INTEGER,            INTENT(IN)    :: KZ0
131!
132REAL, DIMENSION(klon), INTENT(INOUT) :: PSST       ! Sea Surface Temperature (K)
133REAL, DIMENSION(klon), INTENT(INOUT) :: PZ0SEA     ! roughness length over sea
134REAL, DIMENSION(klon), INTENT(OUT)   :: PZ0HSEA    ! heat roughness length over sea
135
136! surface fluxes : latent heat, sensible heat, friction fluxes
137REAL, DIMENSION(klon), INTENT(OUT)   :: PUSTAR     ! friction velocity (m/s)
138REAL, DIMENSION(klon), INTENT(OUT)   :: PSFTH      ! heat flux (W/m2)
139REAL, DIMENSION(klon), INTENT(OUT)   :: PSFTQ      ! water flux (kg/m2/s)
140
141! diagnostics
142REAL, DIMENSION(klon), INTENT(OUT)   :: PCD        ! transfer coef. for momentum
143REAL, DIMENSION(klon), INTENT(OUT)   :: PCH        ! transfer coef. for temperature
144REAL, DIMENSION(klon), INTENT(OUT)   :: PCE        ! transfer coef. for humidity
145REAL, DIMENSION(klon), INTENT(OUT)   :: PCDN       ! neutral coef. for momentum
146REAL, DIMENSION(klon), INTENT(OUT)   :: PRI        ! Richardson number
147REAL, DIMENSION(klon), INTENT(OUT)   :: PRESA      ! aerodynamical resistance
148real, dimension(3), intent(out)   :: coeffs
149
150!       0.2. Declarations of local variables
151!
152! specif SB
153INTEGER, DIMENSION(SIZE(PTA))     :: JCV        ! convergence index
154INTEGER, DIMENSION(SIZE(PTA))     :: JITER      ! nb of iterations to converge
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))        :: ZTAU       ! momentum flux (N/m2)
160REAL, DIMENSION(SIZE(PTA))        :: ZHF        ! sensible heat flux (W/m2)
161REAL, DIMENSION(SIZE(PTA))        :: ZEF        ! latent heat flux (W/m2)
162REAL, DIMENSION(SIZE(PTA))        :: ZTAUR      ! momentum flx due to rain (N/m2)
163REAL, DIMENSION(SIZE(PTA))        :: ZRF        ! sensible flx due to rain (W/m2)
164REAL, DIMENSION(SIZE(PTA))        :: ZEFWEBB    ! Webb corr. on latent flx (W/m2)
165
166REAL, DIMENSION(SIZE(PTA))        :: ZVMOD      ! wind intensity at atm level (m/s)
167REAL, DIMENSION(SIZE(PTA))        :: ZQSATA     ! sat.spec.hum. at atm level (kg/kg)
168REAL, DIMENSION(SIZE(PTA))        :: ZLVA       ! vap.heat of pure water at atm level (J/kg)
169REAL, DIMENSION(SIZE(PTA))        :: ZLVS       ! vap.heat of seawater at sea surface (J/kg)
170REAL, DIMENSION(SIZE(PTA))        :: ZCPA       ! specif.heat moist air (J/kg/K)
171REAL, DIMENSION(SIZE(PTA))        :: ZVISA      ! kinemat.visc. of dry air (m2/s)
172REAL, DIMENSION(SIZE(PTA))        :: ZDU        ! U   vert.grad. (real atm)
173REAL, DIMENSION(SIZE(PTA))        :: ZDT,ZDQ    ! T,Q vert.grad. (real atm)
174REAL, DIMENSION(SIZE(PTA))        :: ZDDU       ! U   vert.grad. (real atm + gust)
175REAL, DIMENSION(SIZE(PTA))        :: ZDDT,ZDDQ  ! T,Q vert.grad. (real atm + WL/CS)
176REAL, DIMENSION(SIZE(PTA))        :: ZUSR       ! velocity scaling param. (m/s)
177                                                ! =friction velocity
178REAL, DIMENSION(SIZE(PTA))        :: ZTSR       ! temperature scaling param. (K)
179REAL, DIMENSION(SIZE(PTA))        :: ZQSR       ! humidity scaling param. (kg/kg)
180REAL, DIMENSION(SIZE(PTA))        :: ZDELTAU10N,ZDELTAT10N,ZDELTAQ10N
181                                                ! U,T,Q vert.grad. (10m, neutral atm)
182REAL, DIMENSION(SIZE(PTA))        :: ZUSR0,ZTSR0,ZQSR0    ! ITERATIVE PROCESS
183REAL, DIMENSION(SIZE(PTA))        :: ZDUSTO,ZDTSTO,ZDQSTO ! ITERATIVE PROCESS
184REAL, DIMENSION(SIZE(PTA))        :: ZPSIU,ZPSIT! PSI funct for U, T/Q (Z0 comp)
185REAL, DIMENSION(SIZE(PTA))        :: ZCHARN     ! Charnock parameter   (Z0 comp)
186
187REAL, DIMENSION(SIZE(PTA))        :: ZUSTAR2    ! square of friction velocity
188REAL, DIMENSION(SIZE(PTA))        :: ZAC        ! aerodynamical conductance
189REAL, DIMENSION(SIZE(PTA))        :: ZDIRCOSZW  ! orography slope cosine
190                                                ! (=1 on water!)
191REAL, DIMENSION(SIZE(PTA))        :: ZPARUN,ZPARTN,ZPARQN ! neutral parameter for U,T,Q
192
193!-- rajout pour la pression saturante
194REAL, DIMENSION(SIZE(PPA))        :: ZFOES                                 ! [OPWEBB]
195REAL, DIMENSION(SIZE(PPA))        :: ZWORK1
196REAL, DIMENSION(SIZE(PPA))        :: ZWORK2 
197REAL, DIMENSION(SIZE(PPS))        :: ZWORK1A
198REAL, DIMENSION(SIZE(PPS))        :: ZWORK2A
199!#####################
200
201REAL, DIMENSION(0:5)              :: ZCOEFU,ZCOEFT,ZCOEFQ
202
203!--------- Modif Olive -----------------
204REAL, DIMENSION(SIZE(PTA))        :: PRHOA
205REAL, PARAMETER                   :: XUNDEF = 1.E+20
206
207
208REAL       :: XVCHRNK = 0.021
209REAL       :: XVZ0CM = 1.0E-5
210!REAL       :: XRIMAX
211
212CHARACTER  :: CCHARNOCK = 'NEW'
213
214
215!--------------------------------------
216
217
218! local constants
219LOGICAL :: OPCVFLX              ! to force convergence
220INTEGER :: NITERMAX             ! nb of iterations to get free convergence
221INTEGER :: NITERSUP             ! nb of additional iterations if OPCVFLX=.TRUE.
222INTEGER :: NITERFL              ! maximum number of iterations
223REAL    :: ZETV,ZRDSRV          ! thermodynamic constants
224REAL    :: ZSQR3
225REAL    :: ZLMOMIN,ZLMOMAX      ! min/max value of Obukhovs stability param. z/l
226REAL    :: ZBTA,ZGMA            ! parameters of the stability functions
227REAL    :: ZDUSR0,ZDTSR0,ZDQSR0 ! maximum gap for USR/TSR/QSR between 2 steps
228REAL    :: ZP00                 ! [OPRECIP] - water vap. diffusiv.ref.press.(Pa)
229REAL    :: ZUTU,ZUTT,ZUTQ       ! U10n threshold in ECUME parameterisation
230REAL    :: ZCDIRU,ZCDIRT,ZCDIRQ ! coef directeur pour fonction affine U,T,Q
231REAL    :: ZORDOU,ZORDOT,ZORDOQ ! ordonnee a l'origine pour fonction affine U,T,Q
232
233INTEGER :: JJ                                   ! for ITERATIVE PROCESS
234INTEGER :: JLON,JK
235REAL    :: ZLMOU,ZLMOT                          ! Obukhovs param. z/l for U, T/Q
236REAL    :: ZPSI_U,ZPSI_T                        ! PSI funct. for U, T/Q
237REAL    :: Z0TSEA,Z0QSEA                        ! roughness length for T, Q
238REAL    :: ZCHIC,ZCHIK,ZPSIC,ZPSIK,ZLOGUS10,ZLOGTS10
239REAL    :: ZTAC,ZCPWA,ZDQSDT,ZDWAT,ZDTMP,ZBULB  ! [OPRECIP]
240REAL    :: ZWW                                  ! [OPWEBB]
241
242
243INTEGER :: PREF             ! reference pressure for exner function
244REAL, DIMENSION(klon)    :: PQSATA      ! sea surface spec. hum. (kg/kg)
245
246REAL    :: qsat_seawater2,qsat_seawater
247
248!REAL(KIND=JPRB) :: ZHOOK_HANDLE
249!
250!-------------------------------------------------------------------------------
251!----------------------- Modif Olive calcul de PRHOA ---------------------------
252
253!write(*,*) "PZ0SEA ",PZ0SEA
254!write(*,*) "PTA ",PTA
255!write(*,*) "PSST ",PSST
256!write(*,*) "PQA ",PQA
257!write(*,*) "PVMOD ",PVMOD
258!write(*,*) "PZREF ",PZREF
259!write(*,*) "PUREF ",PUREF
260!write(*,*) "PPS ",PPS
261!write(*,*) "PPA ",PPA
262!write(*,*) "OPRECIP ",OPRECIP
263!write(*,*) "PZ0HSEA ",PZ0HSEA
264!write(*,*) "PRAIN ",PRAIN
265
266
267PRHOA(:) = PPS(:) / (287.1 * PTA(:) * (1.+.61*PQA(:)))
268!write(*,*) "klon klon ",klon,PTA
269!write(*,*) "PRHOA ",SIZE(PRHOA),PRHOA
270
271PREF = 100900.                    ! = 1000 hPa
272
273!PEXNA = (PPA/PPS)**RKAPPA
274!PEXNS = (PPS/PPS)**RKAPPA
275
276PEXNA = (PPA/PREF)**(RD/RCPD)
277PEXNS = (PPS/PREF)**(RD/RCPD)
278
279!IF (LHOOK) CALL DR_HOOK('ECUMEV6_FLUX',0,ZHOOK_HANDLE)
280!
281ZDUSR0   = 1.E-06
282ZDTSR0   = 1.E-06
283ZDQSR0   = 1.E-09
284!
285NITERMAX = 5
286NITERSUP = 5
287OPCVFLX  = .TRUE.
288!
289NITERFL = NITERMAX
290IF (OPCVFLX) NITERFL = NITERMAX+NITERSUP
291!
292ZCOEFU = (/ 1.00E-03, 3.66E-02, -1.92E-03, 2.32E-04, -7.02E-06,  6.40E-08 /)
293ZCOEFT = (/ 5.36E-03, 2.90E-02, -1.24E-03, 4.50E-04, -2.06E-05,       0.0 /)
294ZCOEFQ = (/ 1.00E-03, 3.59E-02, -2.87E-04,      0.0,       0.0,       0.0 /)
295!
296ZUTU = 40.0
297ZUTT = 14.4
298ZUTQ = 10.0
299!
300ZCDIRU = ZCOEFU(1) + 2.0*ZCOEFU(2)*ZUTU + 3.0*ZCOEFU(3)*ZUTU**2   &
301                   + 4.0*ZCOEFU(4)*ZUTU**3 + 5.0*ZCOEFU(5)*ZUTU**4
302ZCDIRT = ZCOEFT(1) + 2.0*ZCOEFT(2)*ZUTT + 3.0*ZCOEFT(3)*ZUTT**2   &
303                   + 4.0*ZCOEFT(4)*ZUTT**3
304ZCDIRQ = ZCOEFQ(1) + 2.0*ZCOEFQ(2)*ZUTQ
305!
306ZORDOU = ZCOEFU(0) + ZCOEFU(1)*ZUTU + ZCOEFU(2)*ZUTU**2 + ZCOEFU(3)*ZUTU**3   &
307                   + ZCOEFU(4)*ZUTU**4 + ZCOEFU(5)*ZUTU**5
308ZORDOT = ZCOEFT(0) + ZCOEFT(1)*ZUTT + ZCOEFT(2)*ZUTT**2 + ZCOEFT(3)*ZUTT**3   &
309                   + ZCOEFT(4)*ZUTT**4
310ZORDOQ = ZCOEFQ(0) + ZCOEFQ(1)*ZUTQ + ZCOEFQ(2)*ZUTQ**2
311!
312!-------------------------------------------------------------------------------
313!
314!       1.   AUXILIARY CONSTANTS & ARRAY INITIALISATION BY UNDEFINED VALUES.
315!       --------------------------------------------------------------------
316!
317ZDIRCOSZW(:) = 1.0
318!
319ZETV    = XRV/XRD-1.0   !~0.61 (cf Liu et al. 1979)
320ZRDSRV  = XRD/XRV       !~0.622
321ZSQR3   = SQRT(3.0)
322ZLMOMIN = -200.0
323ZLMOMAX = 0.25
324ZBTA    = 16.0
325ZGMA    = 7.0           !initially =4.7, modified to 7.0 following G. Caniaux
326!
327ZP00    = 1013.25E+02
328!
329PCD  = XUNDEF
330PCH  = XUNDEF
331PCE  = XUNDEF
332PCDN = XUNDEF
333ZUSR = XUNDEF
334ZTSR = XUNDEF
335ZQSR = XUNDEF
336ZTAU = XUNDEF
337ZHF  = XUNDEF
338ZEF  = XUNDEF
339!
340PSFTH  = XUNDEF
341PSFTQ  = XUNDEF
342PUSTAR = XUNDEF
343PRESA  = XUNDEF
344PRI    = XUNDEF
345!
346ZTAUR   = 0.0
347ZRF     = 0.0
348ZEFWEBB = 0.0
349!
350!-------------------------------------------------------------------------------
351!
352!       2.   INITIALISATIONS BEFORE ITERATIVE LOOP.
353!       -------------------------------------------
354!
355!ZVMOD(:) = WIND_THRESHOLD(PVMOD(:),PUREF(:))    !set a minimum value to wind
356ZVMOD = MAX(PVMOD , 0.1 * MIN(10.,PUREF) )    !set a minimum value to wind
357
358write(*,*) "ZVMOD ",SIZE(ZVMOD)
359
360!
361!       2.0. Radiative fluxes - For warm layer & cool skin
362!
363!       2.0b. Warm Layer correction
364!
365!       2.1. Specific humidity at saturation
366!
367WHERE(PSSS(:)>0.0.AND.PSSS(:)/=XUNDEF)
368  PQSATA = QSAT_SEAWATER2 (PSST(:),PPS(:),PSSS(:))    !at sea surface
369ELSEWHERE
370  PQSATA (:) = QSAT_SEAWATER (PSST(:),PPS(:))            !at sea surface
371ENDWHERE
372
373!ZQSATA(:) = QSAT(PTA(:),PPA(:))                         !at atm level
374
375!### OLIVIER POUR PRESSION SATURANTE #####
376!-------------------------------------------------------------------------------
377!
378ZFOES  = 1 !PSAT(PT(:))
379ZFOES  = 0.98*ZFOES
380!
381ZWORK1    = ZFOES/PPS
382ZWORK2    = XRD/XRV
383
384ZWORK1A    = ZFOES/PPA
385ZWORK2A    = XRD/XRV
386
387!write(*,*) "ZFOES ",ZFOES
388!write(*,*) "PPS ",PPS
389!write(*,*) "ZWORK1 ",ZWORK1
390!write(*,*) "XRD ",XRD
391!write(*,*) "XRV ",XRV
392!write(*,*) "PPA ",PPA
393!write(*,*) "ZWORK1A ",ZWORK1A
394
395write(*,*) "PQSAT : ",PQSAT
396write(*,*) "PQSATA : ",PQSATA
397
398
399!
400!*       2.    COMPUTE SATURATION HUMIDITY
401!              ---------------------------
402!
403!PQSAT  = ZWORK2*ZWORK1 / (1.+(ZWORK2-1.)*ZWORK1)
404!ZQSATA = ZWORK2A*ZWORK1A / (1.+(ZWORK2A-1.)*ZWORK1A)
405ZQSATA = QSAT_SEAWATER (PTA(:),PPA(:))            !at sea surface
406
407
408!
409!       2.2. Gradients at the air-sea interface
410!
411ZDU(:) = ZVMOD(:)               !one assumes u is measured / sea surface current
412ZDT(:) = PTA(:)/PEXNA(:)-PSST(:)/PEXNS(:)
413ZDQ(:) = PQA(:)-PQSATA(:)
414
415write(*,*) "PQA ",PQA(:)
416write(*,*) "PQSAT",PQSAT(:)
417write(*,*) "ZDQ",ZDQ(:)
418!
419!       2.3. Latent heat of vaporisation
420!
421ZLVA(:) = XLVTT+(XCPV-XCL)*(PTA (:)-XTT)                !of pure water at atm level
422ZLVS(:) = XLVTT+(XCPV-XCL)*(PSST(:)-XTT)                !of pure water at sea surface
423
424
425
426write(*,*) "ZLVA ",ZLVA
427write(*,*) "ZLVS ",ZLVS
428
429
430WHERE(PSSS(:)>0.0.AND.PSSS(:)/=XUNDEF)
431  ZLVS(:) = ZLVS(:)*(1.0-1.00472E-3*PSSS(:))            !of seawater at sea surface
432ENDWHERE
433!
434!       2.4. Specific heat of moist air (Businger 1982)
435!
436!ZCPA(:) = XCPD*(1.0+(XCPV/XCPD-1.0)*PQA(:))
437ZCPA(:) = XCPD
438!
439!       2.4b Kinematic viscosity of dry air (Andreas 1989, CRREL Rep. 89-11)
440!
441ZVISA(:) = 1.326E-05*(1.0+6.542E-03*(PTA(:)-XTT)+8.301E-06*(PTA(:)-XTT)**2   &
442           -4.84E-09*(PTA(:)-XTT)**3)
443!
444!       2.4c Coefficients for warm layer and/or cool skin correction
445!
446!       2.5. Initial guess
447!
448ZDDU(:) = ZDU(:)
449ZDDT(:) = ZDT(:)
450ZDDQ(:) = ZDQ(:)
451ZDDU(:) = SIGN(MAX(ABS(ZDDU(:)),10.0*ZDUSR0),ZDDU(:))
452ZDDT(:) = SIGN(MAX(ABS(ZDDT(:)),10.0*ZDTSR0),ZDDT(:))
453ZDDQ(:) = SIGN(MAX(ABS(ZDDQ(:)),10.0*ZDQSR0),ZDDQ(:))
454
455write(*,*) "ZDDU ",ZDDU
456write(*,*) "ZDDQ ",ZDDQ
457write(*,*) "ZDDT ",ZDDT
458
459!
460JCV (:) = -1
461ZUSR(:) = 0.04*ZDDU(:)
462ZTSR(:) = 0.04*ZDDT(:)
463ZQSR(:) = 0.04*ZDDQ(:)
464ZDELTAU10N(:) = ZDDU(:)
465ZDELTAT10N(:) = ZDDT(:)
466ZDELTAQ10N(:) = ZDDQ(:)
467JITER(:) = 99
468!
469! In the following, we suppose that Richardson number PRI < XRIMAX
470! If not true, Monin-Obukhov theory can't (and therefore shouldn't) be applied !
471!-------------------------------------------------------------------------------
472!
473!       3.   ITERATIVE LOOP TO COMPUTE U*, T*, Q*.
474!       ------------------------------------------
475!
476DO JJ=1,NITERFL
477  DO JLON=1,SIZE(PTA)
478!
479  IF (JCV(JLON) == -1) THEN
480    ZUSR0(JLON)=ZUSR(JLON)
481    ZTSR0(JLON)=ZTSR(JLON)
482    ZQSR0(JLON)=ZQSR(JLON)
483    IF (JJ == NITERMAX+1 .OR. JJ == NITERMAX+NITERSUP) THEN
484      ZDELTAU10N(JLON) = 0.5*(ZDUSTO(JLON)+ZDELTAU10N(JLON))    !forced convergence
485      ZDELTAT10N(JLON) = 0.5*(ZDTSTO(JLON)+ZDELTAT10N(JLON))
486      ZDELTAQ10N(JLON) = 0.5*(ZDQSTO(JLON)+ZDELTAQ10N(JLON))
487      IF (JJ == NITERMAX+NITERSUP) JCV(JLON)=3
488    ENDIF
489    ZDUSTO(JLON) = ZDELTAU10N(JLON)
490    ZDTSTO(JLON) = ZDELTAT10N(JLON)
491    ZDQSTO(JLON) = ZDELTAQ10N(JLON)
492!
493!       3.1. Neutral parameter for wind speed (ECUME_V6 formulation)
494!
495    IF (ZDELTAU10N(JLON) <= ZUTU) THEN
496      ZPARUN(JLON) = ZCOEFU(0) + ZCOEFU(1)*ZDELTAU10N(JLON)      &
497                               + ZCOEFU(2)*ZDELTAU10N(JLON)**2   &
498                               + ZCOEFU(3)*ZDELTAU10N(JLON)**3   &
499                               + ZCOEFU(4)*ZDELTAU10N(JLON)**4   &
500                               + ZCOEFU(5)*ZDELTAU10N(JLON)**5
501    ELSE
502      ZPARUN(JLON) = ZCDIRU*(ZDELTAU10N(JLON)-ZUTU) + ZORDOU
503    ENDIF
504    PCDN(JLON) = (ZPARUN(JLON)/ZDELTAU10N(JLON))**2
505!
506!       3.2. Neutral parameter for temperature (ECUME_V6 formulation)
507!
508    IF (ZDELTAU10N(JLON) <= ZUTT) THEN
509      ZPARTN(JLON) = ZCOEFT(0) + ZCOEFT(1)*ZDELTAU10N(JLON)      &
510                               + ZCOEFT(2)*ZDELTAU10N(JLON)**2   &
511                               + ZCOEFT(3)*ZDELTAU10N(JLON)**3   &
512                               + ZCOEFT(4)*ZDELTAU10N(JLON)**4
513    ELSE
514      ZPARTN(JLON) = ZCDIRT*(ZDELTAU10N(JLON)-ZUTT) + ZORDOT
515    ENDIF
516!
517!       3.3. Neutral parameter for humidity (ECUME_V6 formulation)
518!
519    IF (ZDELTAU10N(JLON) <= ZUTQ) THEN
520      ZPARQN(JLON) = ZCOEFQ(0) + ZCOEFQ(1)*ZDELTAU10N(JLON)      &
521                               + ZCOEFQ(2)*ZDELTAU10N(JLON)**2
522    ELSE
523      ZPARQN(JLON) = ZCDIRQ*(ZDELTAU10N(JLON)-ZUTQ) + ZORDOQ
524    ENDIF
525!
526!       3.4. Scaling parameters U*, T*, Q*
527!
528    ZUSR(JLON) = ZPARUN(JLON)
529    ZTSR(JLON) = ZPARTN(JLON)*ZDELTAT10N(JLON)/ZDELTAU10N(JLON)
530    ZQSR(JLON) = ZPARQN(JLON)*ZDELTAQ10N(JLON)/ZDELTAU10N(JLON)
531!
532!       3.4b Gustiness factor (Deardorff 1970)
533!
534!       3.4c Cool skin correction
535!
536!       3.5. Obukhovs stability param. z/l following Liu et al. (JAS, 1979)
537!
538! For U
539    ZLMOU = PUREF(JLON)*XG*XKARMAN*(ZTSR(JLON)/PTA(JLON)   &
540            +ZETV*ZQSR(JLON)/(1.0+ZETV*PQA(JLON)))/ZUSR(JLON)**2
541! For T/Q
542    ZLMOT = ZLMOU*(PZREF(JLON)/PUREF(JLON))
543    ZLMOU = MAX(MIN(ZLMOU,ZLMOMAX),ZLMOMIN)
544    ZLMOT = MAX(MIN(ZLMOT,ZLMOMAX),ZLMOMIN)
545!
546!       3.6. Stability function psi (see Liu et al, 1979 ; Dyer and Hicks, 1970)
547!            Modified to include convective form following Fairall (unpublished)
548!
549! For U
550    IF (ZLMOU == 0.0) THEN
551      ZPSI_U = 0.0
552    ELSEIF (ZLMOU > 0.0) THEN
553      ZPSI_U = -ZGMA*ZLMOU
554    ELSE
555      ZCHIK  = (1.0-ZBTA*ZLMOU)**0.25
556      ZPSIK  = 2.0*LOG((1.0+ZCHIK)/2.0)  &
557                +LOG((1.0+ZCHIK**2)/2.0) &
558                -2.0*ATAN(ZCHIK)+0.5*XPI
559      ZCHIC  = (1.0-12.87*ZLMOU)**(1.0/3.0)       !for very unstable conditions
560      ZPSIC  = 1.5*LOG((ZCHIC**2+ZCHIC+1.0)/3.0)   &
561                -ZSQR3*ATAN((2.0*ZCHIC+1.0)/ZSQR3) &
562                +XPI/ZSQR3
563      ZPSI_U = ZPSIC+(ZPSIK-ZPSIC)/(1.0+ZLMOU**2) !match Kansas & free-conv. forms
564    ENDIF
565    ZPSIU(JLON) = ZPSI_U
566! For T/Q
567    IF (ZLMOT == 0.0) THEN
568      ZPSI_T = 0.0
569    ELSEIF (ZLMOT > 0.0) THEN
570      ZPSI_T = -ZGMA*ZLMOT
571    ELSE
572      ZCHIK  = (1.0-ZBTA*ZLMOT)**0.25
573      ZPSIK  = 2.0*LOG((1.0+ZCHIK**2)/2.0)
574      ZCHIC  = (1.0-12.87*ZLMOT)**(1.0/3.0)       !for very unstable conditions
575      ZPSIC  = 1.5*LOG((ZCHIC**2+ZCHIC+1.0)/3.0)   &
576                -ZSQR3*ATAN((2.0*ZCHIC+1.0)/ZSQR3) &
577                +XPI/ZSQR3
578      ZPSI_T = ZPSIC+(ZPSIK-ZPSIC)/(1.0+ZLMOT**2) !match Kansas & free-conv. forms
579    ENDIF
580    ZPSIT(JLON) = ZPSI_T
581!
582!       3.7. Update air-sea gradients
583!
584    ZDDU(JLON) = ZDU(JLON)
585    ZDDT(JLON) = ZDT(JLON)
586    ZDDQ(JLON) = ZDQ(JLON)
587    ZDDU(JLON) = SIGN(MAX(ABS(ZDDU(JLON)),10.0*ZDUSR0),ZDDU(JLON))
588    ZDDT(JLON) = SIGN(MAX(ABS(ZDDT(JLON)),10.0*ZDTSR0),ZDDT(JLON))
589    ZDDQ(JLON) = SIGN(MAX(ABS(ZDDQ(JLON)),10.0*ZDQSR0),ZDDQ(JLON))
590    ZLOGUS10   = LOG(PUREF(JLON)/10.0)
591    ZLOGTS10   = LOG(PZREF(JLON)/10.0)
592    ZDELTAU10N(JLON) = ZDDU(JLON)-ZUSR(JLON)*(ZLOGUS10-ZPSI_U)/XKARMAN
593    ZDELTAT10N(JLON) = ZDDT(JLON)-ZTSR(JLON)*(ZLOGTS10-ZPSI_T)/XKARMAN
594    ZDELTAQ10N(JLON) = ZDDQ(JLON)-ZQSR(JLON)*(ZLOGTS10-ZPSI_T)/XKARMAN
595    ZDELTAU10N(JLON) = SIGN(MAX(ABS(ZDELTAU10N(JLON)),10.0*ZDUSR0),   &
596                            ZDELTAU10N(JLON))
597    ZDELTAT10N(JLON) = SIGN(MAX(ABS(ZDELTAT10N(JLON)),10.0*ZDTSR0),   &
598                            ZDELTAT10N(JLON))
599    ZDELTAQ10N(JLON) = SIGN(MAX(ABS(ZDELTAQ10N(JLON)),10.0*ZDQSR0),   &
600                            ZDELTAQ10N(JLON))
601!
602!       3.8. Test convergence for U*, T*, Q*
603!
604    IF (ABS(ZUSR(JLON)-ZUSR0(JLON)) < ZDUSR0 .AND.   &
605        ABS(ZTSR(JLON)-ZTSR0(JLON)) < ZDTSR0 .AND.   &
606        ABS(ZQSR(JLON)-ZQSR0(JLON)) < ZDQSR0) THEN
607      JCV(JLON) = 1                                     !free convergence
608      IF (JJ >= NITERMAX+1) JCV(JLON) = 2               !leaded convergence
609    ENDIF
610    JITER(JLON) = JJ
611  ENDIF
612!
613  ENDDO
614ENDDO
615!
616!-------------------------------------------------------------------------------
617!
618!       4.   COMPUTATION OF TURBULENT FLUXES AND EXCHANGE COEFFICIENTS.
619!       ---------------------------------------------------------------
620!
621DO JLON=1,SIZE(PTA)
622!
623!       4.1. Surface turbulent fluxes
624!            (ATM CONV.: ZTAU<<0 ; ZHF,ZEF<0 if atm looses heat)
625!
626  ZTAU(JLON) = -PRHOA(JLON)*ZUSR(JLON)**2
627  ZHF(JLON)  = -PRHOA(JLON)*ZCPA(JLON)*ZUSR(JLON)*ZTSR(JLON)
628  ZEF(JLON)  = -PRHOA(JLON)*ZLVS(JLON)*ZUSR(JLON)*ZQSR(JLON)
629
630  write(*,*) "ZTAU = ",ZTAU(JLON)
631  write(*,*) "SENS = ",ZHF(JLON)
632  write(*,*) "LAT = ",ZEF(JLON)
633
634!
635!       4.2. Exchange coefficients PCD, PCH, PCE
636!
637  PCD(JLON) = (ZUSR(JLON)/ZDDU(JLON))**2
638  PCH(JLON) = (ZUSR(JLON)*ZTSR(JLON))/(ZDDU(JLON)*ZDDT(JLON))
639  PCE(JLON) = (ZUSR(JLON)*ZQSR(JLON))/(ZDDU(JLON)*ZDDQ(JLON))
640
641  write(*,*) "ZUSR = ",ZUSR(JLON)
642  write(*,*) "ZTSR = ",ZTSR(JLON)
643  write(*,*) "ZQSR = ",ZQSR(JLON)
644
645!
646!       4.3. Stochastic perturbation of turbulent fluxes
647!
648!  IF( OPERTFLUX )THEN
649!    ZTAU(JLON) = ZTAU(JLON)* ( 1. + PPERTFLUX(JLON) / 2. )
650!    ZHF (JLON) = ZHF(JLON)*  ( 1. + PPERTFLUX(JLON) / 2. )
651!    ZEF (JLON) = ZEF(JLON)*  ( 1. + PPERTFLUX(JLON) / 2. )
652!  ENDIF
653!
654ENDDO
655!
656!-------------------------------------------------------------------------------
657!
658!       5.   COMPUTATION OF FLUX CORRECTIONS DUE TO RAINFALL.
659!            (ATM conv: ZRF<0 if atm. looses heat, ZTAUR<<0)
660!       -----------------------------------------------------
661!
662IF (OPRECIP) THEN
663  DO JLON=1,SIZE(PTA)
664!
665!       5.1. Momentum flux due to rainfall (ZTAUR, N/m2)
666!
667! See pp3752 in FBR96.
668    ZTAUR(JLON) = -0.85*PRAIN(JLON)*PVMOD(JLON)
669!
670!       5.2. Sensible heat flux due to rainfall (ZRF, W/m2)
671!
672! See Eq.12 in GoF95 with ZCPWA as specific heat of water at atm level (J/kg/K),
673! ZDQSDT from Clausius-Clapeyron relation, ZDWAT as water vapor diffusivity
674! (Eq.13-3 of Pruppacher and Klett, 1978), ZDTMP as heat diffusivity, and ZBULB
675! as wet-bulb factor (Eq.11 in GoF95).
676!
677    ZTAC   = PTA(JLON)-XTT
678    ZCPWA  = 4217.51 -3.65566*ZTAC +0.1381*ZTAC**2       &
679              -2.8309E-03*ZTAC**3 +3.42061E-05*ZTAC**4   &
680              -2.18107E-07*ZTAC**5 +5.74535E-10*ZTAC**6
681    ZDQSDT = (ZLVA(JLON)*ZQSATA(JLON))/(XRV*PTA(JLON)**2)
682    ZDWAT  = 2.11E-05*(ZP00/PPA(JLON))*(PTA(JLON)/XTT)**1.94
683    ZDTMP  = (1.0+3.309E-03*ZTAC-1.44E-06*ZTAC**2)   &
684              *0.02411/(PRHOA(JLON)*ZCPA(JLON))
685    ZBULB  = 1.0/(1.0+ZDQSDT*(ZLVA(JLON)*ZDWAT)/(ZCPA(JLON)*ZDTMP))
686    ZRF(JLON) = PRAIN(JLON)*ZCPWA*ZBULB*((PSST(JLON)-PTA(JLON))   &
687                +(PQSATA(JLON)-PQA(JLON))*(ZLVA(JLON)*ZDWAT)/(ZCPA(JLON)*ZDTMP))
688!
689  ENDDO
690ENDIF
691!
692!-------------------------------------------------------------------------------
693!
694!       6.   COMPUTATION OF WEBB CORRECTION TO LATENT HEAT FLUX (ZEFWEBB, W/m2).
695!       ------------------------------------------------------------------------
696!
697! See Eq.21 and Eq.22 in FBR96.
698IF (OPWEBB) THEN
699  DO JLON=1,SIZE(PTA)
700    ZWW = (1.0+ZETV)*(ZUSR(JLON)*ZQSR(JLON))   &
701           +(1.0+(1.0+ZETV)*PQA(JLON))*(ZUSR(JLON)*ZTSR(JLON))/PTA(JLON)
702    ZEFWEBB(JLON) = -PRHOA(JLON)*ZLVS(JLON)*ZWW*PQA(JLON)
703  ENDDO
704ENDIF
705!
706!-------------------------------------------------------------------------------
707!
708!       7.   FINAL STEP : TOTAL SURFACE FLUXES AND DERIVED DIAGNOSTICS.
709!       ---------------------------------------------------------------
710!
711!       7.1. Richardson number
712!
713! CALL SURFACE_RI(PSST,PQSAT,PEXNS,PEXNA,PTA,PQA,   &
714!                PZREF,PUREF,ZDIRCOSZW,PVMOD,PRI)
715!
716!       7.2. Friction velocity which contains correction due to rain
717!
718ZUSTAR2(:) = -(ZTAU(:)+ZTAUR(:))/PRHOA(:)       !>>0 as ZTAU<<0 & ZTAUR<=0
719!
720IF (OPRECIP) THEN
721  PCD(:) = ZUSTAR2(:)/ZDDU(:)**2
722ENDIF
723!
724PUSTAR(:) = SQRT(ZUSTAR2(:))                    !>>0
725!
726!       7.3. Aerodynamical conductance and resistance
727!
728ZAC  (:) = PCH(:)*ZDDU(:)
729PRESA(:) = 1.0/ZAC(:)
730!
731!       7.4. Total surface fluxes
732!
733PSFTH(:) =  ZHF(:)+ZRF(:)
734PSFTQ(:) = (ZEF(:)+ZEFWEBB(:))/ZLVS(:)
735!
736!       7.5. Charnock number
737!
738IF (CCHARNOCK == 'OLD') THEN
739  ZCHARN(:) = XVCHRNK
740ELSE            !modified for moderate wind speed as in COARE3.0
741  ZCHARN(:) = MIN(0.018,MAX(0.011,0.011+(0.007/8.0)*(ZDDU(:)-10.0)))
742ENDIF
743!
744!       7.6. Roughness lengths Z0 and Z0H over sea
745!
746!IF (KZ0 == 0) THEN      ! ARPEGE formulation
747!  PZ0SEA (:) = (ZCHARN(:)/XG)*ZUSTAR2(:) + XVZ0CM*PCD(:)/PCDN(:)
748!  PZ0HSEA(:) = PZ0SEA (:)
749!ELSEIF (KZ0 == 1) THEN  ! Smith (1988) formulation
750!  PZ0SEA (:) = (ZCHARN(:)/XG)*ZUSTAR2(:) + 0.11*ZVISA(:)/PUSTAR(:)
751!  PZ0HSEA(:) = PZ0SEA (:)
752!ELSEIF (KZ0 == 2) THEN  ! Direct computation using the stability functions
753!  DO JLON=1,SIZE(PTA)
754!    PZ0SEA (JLON) = PUREF(JLON)/EXP(XKARMAN*ZDDU(JLON)/PUSTAR(JLON)+ZPSIU(JLON))
755!    Z0TSEA        = PZREF(JLON)/EXP(XKARMAN*ZDDT(JLON)/ZTSR  (JLON)+ZPSIT(JLON))
756!    Z0QSEA        = PZREF(JLON)/EXP(XKARMAN*ZDDQ(JLON)/ZQSR  (JLON)+ZPSIT(JLON))
757!    PZ0HSEA(JLON) = 0.5*(Z0TSEA+Z0QSEA)
758!  ENDDO
759!ENDIF
760
761write(*,*) "JLON ",JLON
762write(*,*) "PTA ",klon,PTA
763write(*,*) "PCD ",SIZE(PCD),PCD
764write(*,*) "PCQ ",SIZE(PCE),PCE
765write(*,*) "PCH ",SIZE(PCH),PCH
766
767coeffs = [PCD,&
768       PCE,&
769       PCH]
770!
771!
772!IF (LHOOK) CALL DR_HOOK('ECUMEV6_FLUX',1,ZHOOK_HANDLE)
773!
774!-------------------------------------------------------------------------------
775   END SUBROUTINE ECUMEV6_FLUX
Note: See TracBrowser for help on using the repository browser.