source: LMDZ6/trunk/libf/phylmd/ecumev6_flux.F90 @ 5202

Last change on this file since 5202 was 4722, checked in by Laurent Fairhead, 13 months ago

Modification by O. Torres to the cdrag routines to include different bulk formulae
to calculate cdrag coefficients over ocean as well as an iteration of that
calculation.
The iteration is controlled by flag ok_cdrag_iter which if set to FALSE by default
to converge with previous results.
The choice of bulk formulae is set with the choix_bulk parameter
The number of iterations to run is set with nit_bulk
OT, PB, CD, LF
File size: 28.3 KB
RevLine 
[4661]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.