source: LMDZ6/branches/LMDZ_cdrag_LSCE/libf/phylmd/ecumev6_flux.F90

Last change on this file was 4661, checked in by Laurent Fairhead, 14 months ago

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