source: LMDZ6/branches/Amaury_dev/libf/phylmd/ecumev6_flux.F90 @ 5136

Last change on this file since 5136 was 5117, checked in by abarral, 2 months ago

rename modules properly lmdz_*
move some unused files to obsolete/
(lint) uppercase fortran keywords
File size: 28.2 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!       2.0. Radiative fluxes - For warm layer & cool skin
363
364!       2.0b. Warm Layer correction
365
366!       2.1. Specific humidity at saturation
367
368WHERE(PSSS(:)>0.0.AND.PSSS(:)/=XUNDEF)
369  PQSATA = QSAT_SEAWATER2 (PSST(:),PPS(:),PSSS(:))    !at sea surface
370ELSEWHERE
371  PQSATA (:) = QSAT_SEAWATER (PSST(:),PPS(:))            !at sea surface
372ENDWHERE
373
374!ZQSATA(:) = QSAT(PTA(:),PPA(:))                         !at atm level
375
376!### OLIVIER POUR PRESSION SATURANTE #####
377!-------------------------------------------------------------------------------
378
379ZFOES  = 1 !PSAT(PT(:))
380ZFOES  = 0.98*ZFOES
381
382ZWORK1    = ZFOES/PPS
383ZWORK2    = XRD/XRV
384
385ZWORK1A    = ZFOES/PPA
386ZWORK2A    = XRD/XRV
387
388!WRITE(*,*) "ZFOES ",ZFOES
389!WRITE(*,*) "PPS ",PPS
390!WRITE(*,*) "ZWORK1 ",ZWORK1
391!WRITE(*,*) "XRD ",XRD
392!WRITE(*,*) "XRV ",XRV
393!WRITE(*,*) "PPA ",PPA
394!WRITE(*,*) "ZWORK1A ",ZWORK1A
395
396WRITE(*,*) "PQSAT : ",PQSAT
397WRITE(*,*) "PQSATA : ",PQSATA
398
399!*       2.    COMPUTE SATURATION HUMIDITY
400!              ---------------------------
401
402!PQSAT  = ZWORK2*ZWORK1 / (1.+(ZWORK2-1.)*ZWORK1)
403!ZQSATA = ZWORK2A*ZWORK1A / (1.+(ZWORK2A-1.)*ZWORK1A)
404ZQSATA = QSAT_SEAWATER (PTA(:),PPA(:))            !at sea surface
405
406!       2.2. Gradients at the air-sea interface
407
408ZDU(:) = ZVMOD(:)               !one assumes u is measured / sea surface current
409ZDT(:) = PTA(:)/PEXNA(:)-PSST(:)/PEXNS(:)
410ZDQ(:) = PQA(:)-PQSATA(:)
411
412WRITE(*,*) "PQA ",PQA(:)
413WRITE(*,*) "PQSAT",PQSAT(:)
414WRITE(*,*) "ZDQ",ZDQ(:)
415
416!       2.3. Latent heat of vaporisation
417
418ZLVA(:) = XLVTT+(XCPV-XCL)*(PTA (:)-XTT)                !of pure water at atm level
419ZLVS(:) = XLVTT+(XCPV-XCL)*(PSST(:)-XTT)                !of pure water at sea surface
420
421
422
423WRITE(*,*) "ZLVA ",ZLVA
424WRITE(*,*) "ZLVS ",ZLVS
425
426
427WHERE(PSSS(:)>0.0.AND.PSSS(:)/=XUNDEF)
428  ZLVS(:) = ZLVS(:)*(1.0-1.00472E-3*PSSS(:))            !of seawater at sea surface
429ENDWHERE
430
431!       2.4. Specific heat of moist air (Businger 1982)
432
433!ZCPA(:) = XCPD*(1.0+(XCPV/XCPD-1.0)*PQA(:))
434ZCPA(:) = XCPD
435
436!       2.4b Kinematic viscosity of dry air (Andreas 1989, CRREL Rep. 89-11)
437
438ZVISA(:) = 1.326E-05*(1.0+6.542E-03*(PTA(:)-XTT)+8.301E-06*(PTA(:)-XTT)**2   &
439           -4.84E-09*(PTA(:)-XTT)**3)
440
441!       2.4c Coefficients for warm layer and/or cool skin correction
442
443!       2.5. Initial guess
444
445ZDDU(:) = ZDU(:)
446ZDDT(:) = ZDT(:)
447ZDDQ(:) = ZDQ(:)
448ZDDU(:) = SIGN(MAX(ABS(ZDDU(:)),10.0*ZDUSR0),ZDDU(:))
449ZDDT(:) = SIGN(MAX(ABS(ZDDT(:)),10.0*ZDTSR0),ZDDT(:))
450ZDDQ(:) = SIGN(MAX(ABS(ZDDQ(:)),10.0*ZDQSR0),ZDDQ(:))
451
452WRITE(*,*) "ZDDU ",ZDDU
453WRITE(*,*) "ZDDQ ",ZDDQ
454WRITE(*,*) "ZDDT ",ZDDT
455
456JCV (:) = -1
457ZUSR(:) = 0.04*ZDDU(:)
458ZTSR(:) = 0.04*ZDDT(:)
459ZQSR(:) = 0.04*ZDDQ(:)
460ZDELTAU10N(:) = ZDDU(:)
461ZDELTAT10N(:) = ZDDT(:)
462ZDELTAQ10N(:) = ZDDQ(:)
463JITER(:) = 99
464
465! In the following, we suppose that Richardson number PRI < XRIMAX
466! If not true, Monin-Obukhov theory can't (and therefore shouldn't) be applied !
467!-------------------------------------------------------------------------------
468
469!       3.   ITERATIVE LOOP TO COMPUTE U*, T*, Q*.
470!       ------------------------------------------
471
472DO JJ=1,NITERFL
473  DO JLON=1,SIZE(PTA)
474
475  IF (JCV(JLON) == -1) THEN
476    ZUSR0(JLON)=ZUSR(JLON)
477    ZTSR0(JLON)=ZTSR(JLON)
478    ZQSR0(JLON)=ZQSR(JLON)
479    IF (JJ == NITERMAX+1 .OR. JJ == NITERMAX+NITERSUP) THEN
480      ZDELTAU10N(JLON) = 0.5*(ZDUSTO(JLON)+ZDELTAU10N(JLON))    !forced convergence
481      ZDELTAT10N(JLON) = 0.5*(ZDTSTO(JLON)+ZDELTAT10N(JLON))
482      ZDELTAQ10N(JLON) = 0.5*(ZDQSTO(JLON)+ZDELTAQ10N(JLON))
483      IF (JJ == NITERMAX+NITERSUP) JCV(JLON)=3
484    ENDIF
485    ZDUSTO(JLON) = ZDELTAU10N(JLON)
486    ZDTSTO(JLON) = ZDELTAT10N(JLON)
487    ZDQSTO(JLON) = ZDELTAQ10N(JLON)
488
489!       3.1. Neutral parameter for wind speed (ECUME_V6 formulation)
490
491    IF (ZDELTAU10N(JLON) <= ZUTU) THEN
492      ZPARUN(JLON) = ZCOEFU(0) + ZCOEFU(1)*ZDELTAU10N(JLON)      &
493                               + ZCOEFU(2)*ZDELTAU10N(JLON)**2   &
494                               + ZCOEFU(3)*ZDELTAU10N(JLON)**3   &
495                               + ZCOEFU(4)*ZDELTAU10N(JLON)**4   &
496                               + ZCOEFU(5)*ZDELTAU10N(JLON)**5
497    ELSE
498      ZPARUN(JLON) = ZCDIRU*(ZDELTAU10N(JLON)-ZUTU) + ZORDOU
499    ENDIF
500    PCDN(JLON) = (ZPARUN(JLON)/ZDELTAU10N(JLON))**2
501
502!       3.2. Neutral parameter for temperature (ECUME_V6 formulation)
503
504    IF (ZDELTAU10N(JLON) <= ZUTT) THEN
505      ZPARTN(JLON) = ZCOEFT(0) + ZCOEFT(1)*ZDELTAU10N(JLON)      &
506                               + ZCOEFT(2)*ZDELTAU10N(JLON)**2   &
507                               + ZCOEFT(3)*ZDELTAU10N(JLON)**3   &
508                               + ZCOEFT(4)*ZDELTAU10N(JLON)**4
509    ELSE
510      ZPARTN(JLON) = ZCDIRT*(ZDELTAU10N(JLON)-ZUTT) + ZORDOT
511    ENDIF
512
513!       3.3. Neutral parameter for humidity (ECUME_V6 formulation)
514
515    IF (ZDELTAU10N(JLON) <= ZUTQ) THEN
516      ZPARQN(JLON) = ZCOEFQ(0) + ZCOEFQ(1)*ZDELTAU10N(JLON)      &
517                               + ZCOEFQ(2)*ZDELTAU10N(JLON)**2
518    ELSE
519      ZPARQN(JLON) = ZCDIRQ*(ZDELTAU10N(JLON)-ZUTQ) + ZORDOQ
520    ENDIF
521
522!       3.4. Scaling parameters U*, T*, Q*
523
524    ZUSR(JLON) = ZPARUN(JLON)
525    ZTSR(JLON) = ZPARTN(JLON)*ZDELTAT10N(JLON)/ZDELTAU10N(JLON)
526    ZQSR(JLON) = ZPARQN(JLON)*ZDELTAQ10N(JLON)/ZDELTAU10N(JLON)
527
528!       3.4b Gustiness factor (Deardorff 1970)
529
530!       3.4c Cool skin correction
531
532!       3.5. Obukhovs stability param. z/l following Liu et al. (JAS, 1979)
533
534! For U
535    ZLMOU = PUREF(JLON)*XG*XKARMAN*(ZTSR(JLON)/PTA(JLON)   &
536            +ZETV*ZQSR(JLON)/(1.0+ZETV*PQA(JLON)))/ZUSR(JLON)**2
537! For T/Q
538    ZLMOT = ZLMOU*(PZREF(JLON)/PUREF(JLON))
539    ZLMOU = MAX(MIN(ZLMOU,ZLMOMAX),ZLMOMIN)
540    ZLMOT = MAX(MIN(ZLMOT,ZLMOMAX),ZLMOMIN)
541
542!       3.6. Stability function psi (see Liu et al, 1979 ; Dyer and Hicks, 1970)
543!            Modified to include convective form following Fairall (unpublished)
544
545! For U
546    IF (ZLMOU == 0.0) THEN
547      ZPSI_U = 0.0
548    ELSEIF (ZLMOU > 0.0) THEN
549      ZPSI_U = -ZGMA*ZLMOU
550    ELSE
551      ZCHIK  = (1.0-ZBTA*ZLMOU)**0.25
552      ZPSIK  = 2.0*LOG((1.0+ZCHIK)/2.0)  &
553                +LOG((1.0+ZCHIK**2)/2.0) &
554                -2.0*ATAN(ZCHIK)+0.5*XPI
555      ZCHIC  = (1.0-12.87*ZLMOU)**(1.0/3.0)       !for very unstable conditions
556      ZPSIC  = 1.5*LOG((ZCHIC**2+ZCHIC+1.0)/3.0)   &
557                -ZSQR3*ATAN((2.0*ZCHIC+1.0)/ZSQR3) &
558                +XPI/ZSQR3
559      ZPSI_U = ZPSIC+(ZPSIK-ZPSIC)/(1.0+ZLMOU**2) !match Kansas & free-conv. forms
560    ENDIF
561    ZPSIU(JLON) = ZPSI_U
562! For T/Q
563    IF (ZLMOT == 0.0) THEN
564      ZPSI_T = 0.0
565    ELSEIF (ZLMOT > 0.0) THEN
566      ZPSI_T = -ZGMA*ZLMOT
567    ELSE
568      ZCHIK  = (1.0-ZBTA*ZLMOT)**0.25
569      ZPSIK  = 2.0*LOG((1.0+ZCHIK**2)/2.0)
570      ZCHIC  = (1.0-12.87*ZLMOT)**(1.0/3.0)       !for very unstable conditions
571      ZPSIC  = 1.5*LOG((ZCHIC**2+ZCHIC+1.0)/3.0)   &
572                -ZSQR3*ATAN((2.0*ZCHIC+1.0)/ZSQR3) &
573                +XPI/ZSQR3
574      ZPSI_T = ZPSIC+(ZPSIK-ZPSIC)/(1.0+ZLMOT**2) !match Kansas & free-conv. forms
575    ENDIF
576    ZPSIT(JLON) = ZPSI_T
577
578!       3.7. Update air-sea gradients
579
580    ZDDU(JLON) = ZDU(JLON)
581    ZDDT(JLON) = ZDT(JLON)
582    ZDDQ(JLON) = ZDQ(JLON)
583    ZDDU(JLON) = SIGN(MAX(ABS(ZDDU(JLON)),10.0*ZDUSR0),ZDDU(JLON))
584    ZDDT(JLON) = SIGN(MAX(ABS(ZDDT(JLON)),10.0*ZDTSR0),ZDDT(JLON))
585    ZDDQ(JLON) = SIGN(MAX(ABS(ZDDQ(JLON)),10.0*ZDQSR0),ZDDQ(JLON))
586    ZLOGUS10   = LOG(PUREF(JLON)/10.0)
587    ZLOGTS10   = LOG(PZREF(JLON)/10.0)
588    ZDELTAU10N(JLON) = ZDDU(JLON)-ZUSR(JLON)*(ZLOGUS10-ZPSI_U)/XKARMAN
589    ZDELTAT10N(JLON) = ZDDT(JLON)-ZTSR(JLON)*(ZLOGTS10-ZPSI_T)/XKARMAN
590    ZDELTAQ10N(JLON) = ZDDQ(JLON)-ZQSR(JLON)*(ZLOGTS10-ZPSI_T)/XKARMAN
591    ZDELTAU10N(JLON) = SIGN(MAX(ABS(ZDELTAU10N(JLON)),10.0*ZDUSR0),   &
592                            ZDELTAU10N(JLON))
593    ZDELTAT10N(JLON) = SIGN(MAX(ABS(ZDELTAT10N(JLON)),10.0*ZDTSR0),   &
594                            ZDELTAT10N(JLON))
595    ZDELTAQ10N(JLON) = SIGN(MAX(ABS(ZDELTAQ10N(JLON)),10.0*ZDQSR0),   &
596                            ZDELTAQ10N(JLON))
597
598!       3.8. Test convergence for U*, T*, Q*
599
600    IF (ABS(ZUSR(JLON)-ZUSR0(JLON)) < ZDUSR0 .AND.   &
601        ABS(ZTSR(JLON)-ZTSR0(JLON)) < ZDTSR0 .AND.   &
602        ABS(ZQSR(JLON)-ZQSR0(JLON)) < ZDQSR0) THEN
603      JCV(JLON) = 1                                     !free convergence
604      IF (JJ >= NITERMAX+1) JCV(JLON) = 2               !leaded convergence
605    ENDIF
606    JITER(JLON) = JJ
607  ENDIF
608
609  ENDDO
610ENDDO
611
612!-------------------------------------------------------------------------------
613
614!       4.   COMPUTATION OF TURBULENT FLUXES AND EXCHANGE COEFFICIENTS.
615!       ---------------------------------------------------------------
616
617DO JLON=1,SIZE(PTA)
618
619!       4.1. Surface turbulent fluxes
620!            (ATM CONV.: ZTAU<<0 ; ZHF,ZEF<0 if atm looses heat)
621
622  ZTAU(JLON) = -PRHOA(JLON)*ZUSR(JLON)**2
623  ZHF(JLON)  = -PRHOA(JLON)*ZCPA(JLON)*ZUSR(JLON)*ZTSR(JLON)
624  ZEF(JLON)  = -PRHOA(JLON)*ZLVS(JLON)*ZUSR(JLON)*ZQSR(JLON)
625
626  WRITE(*,*) "ZTAU = ",ZTAU(JLON)
627  WRITE(*,*) "SENS = ",ZHF(JLON)
628  WRITE(*,*) "LAT = ",ZEF(JLON)
629
630!       4.2. Exchange coefficients PCD, PCH, PCE
631
632  PCD(JLON) = (ZUSR(JLON)/ZDDU(JLON))**2
633  PCH(JLON) = (ZUSR(JLON)*ZTSR(JLON))/(ZDDU(JLON)*ZDDT(JLON))
634  PCE(JLON) = (ZUSR(JLON)*ZQSR(JLON))/(ZDDU(JLON)*ZDDQ(JLON))
635
636  WRITE(*,*) "ZUSR = ",ZUSR(JLON)
637  WRITE(*,*) "ZTSR = ",ZTSR(JLON)
638  WRITE(*,*) "ZQSR = ",ZQSR(JLON)
639
640!       4.3. Stochastic perturbation of turbulent fluxes
641
642!  IF( OPERTFLUX )THEN
643!    ZTAU(JLON) = ZTAU(JLON)* ( 1. + PPERTFLUX(JLON) / 2. )
644!    ZHF (JLON) = ZHF(JLON)*  ( 1. + PPERTFLUX(JLON) / 2. )
645!    ZEF (JLON) = ZEF(JLON)*  ( 1. + PPERTFLUX(JLON) / 2. )
646!  ENDIF
647
648ENDDO
649
650!-------------------------------------------------------------------------------
651
652!       5.   COMPUTATION OF FLUX CORRECTIONS DUE TO RAINFALL.
653!            (ATM conv: ZRF<0 if atm. looses heat, ZTAUR<<0)
654!       -----------------------------------------------------
655
656IF (OPRECIP) THEN
657  DO JLON=1,SIZE(PTA)
658
659!       5.1. Momentum flux due to rainfall (ZTAUR, N/m2)
660
661! See pp3752 in FBR96.
662    ZTAUR(JLON) = -0.85*PRAIN(JLON)*PVMOD(JLON)
663
664!       5.2. Sensible heat flux due to rainfall (ZRF, W/m2)
665
666! See Eq.12 in GoF95 with ZCPWA as specific heat of water at atm level (J/kg/K),
667! ZDQSDT from Clausius-Clapeyron relation, ZDWAT as water vapor diffusivity
668! (Eq.13-3 of Pruppacher and Klett, 1978), ZDTMP as heat diffusivity, and ZBULB
669! as wet-bulb factor (Eq.11 in GoF95).
670
671    ZTAC   = PTA(JLON)-XTT
672    ZCPWA  = 4217.51 -3.65566*ZTAC +0.1381*ZTAC**2       &
673              -2.8309E-03*ZTAC**3 +3.42061E-05*ZTAC**4   &
674              -2.18107E-07*ZTAC**5 +5.74535E-10*ZTAC**6
675    ZDQSDT = (ZLVA(JLON)*ZQSATA(JLON))/(XRV*PTA(JLON)**2)
676    ZDWAT  = 2.11E-05*(ZP00/PPA(JLON))*(PTA(JLON)/XTT)**1.94
677    ZDTMP  = (1.0+3.309E-03*ZTAC-1.44E-06*ZTAC**2)   &
678              *0.02411/(PRHOA(JLON)*ZCPA(JLON))
679    ZBULB  = 1.0/(1.0+ZDQSDT*(ZLVA(JLON)*ZDWAT)/(ZCPA(JLON)*ZDTMP))
680    ZRF(JLON) = PRAIN(JLON)*ZCPWA*ZBULB*((PSST(JLON)-PTA(JLON))   &
681                +(PQSATA(JLON)-PQA(JLON))*(ZLVA(JLON)*ZDWAT)/(ZCPA(JLON)*ZDTMP))
682
683  ENDDO
684ENDIF
685
686!-------------------------------------------------------------------------------
687
688!       6.   COMPUTATION OF WEBB CORRECTION TO LATENT HEAT FLUX (ZEFWEBB, W/m2).
689!       ------------------------------------------------------------------------
690
691! See Eq.21 and Eq.22 in FBR96.
692IF (OPWEBB) THEN
693  DO JLON=1,SIZE(PTA)
694    ZWW = (1.0+ZETV)*(ZUSR(JLON)*ZQSR(JLON))   &
695           +(1.0+(1.0+ZETV)*PQA(JLON))*(ZUSR(JLON)*ZTSR(JLON))/PTA(JLON)
696    ZEFWEBB(JLON) = -PRHOA(JLON)*ZLVS(JLON)*ZWW*PQA(JLON)
697  ENDDO
698ENDIF
699
700!-------------------------------------------------------------------------------
701
702!       7.   FINAL STEP : TOTAL SURFACE FLUXES AND DERIVED DIAGNOSTICS.
703!       ---------------------------------------------------------------
704
705!       7.1. Richardson number
706
707! CALL SURFACE_RI(PSST,PQSAT,PEXNS,PEXNA,PTA,PQA,   &
708!                PZREF,PUREF,ZDIRCOSZW,PVMOD,PRI)
709
710!       7.2. Friction velocity which contains correction due to rain
711
712ZUSTAR2(:) = -(ZTAU(:)+ZTAUR(:))/PRHOA(:)       !>>0 as ZTAU<<0 & ZTAUR<=0
713
714IF (OPRECIP) THEN
715  PCD(:) = ZUSTAR2(:)/ZDDU(:)**2
716ENDIF
717
718PUSTAR(:) = SQRT(ZUSTAR2(:))                    !>>0
719
720!       7.3. Aerodynamical conductance and resistance
721
722ZAC  (:) = PCH(:)*ZDDU(:)
723PRESA(:) = 1.0/ZAC(:)
724
725!       7.4. Total surface fluxes
726
727PSFTH(:) =  ZHF(:)+ZRF(:)
728PSFTQ(:) = (ZEF(:)+ZEFWEBB(:))/ZLVS(:)
729
730!       7.5. Charnock number
731
732IF (CCHARNOCK == 'OLD') THEN
733  ZCHARN(:) = XVCHRNK
734ELSE            !modified for moderate wind speed as in COARE3.0
735  ZCHARN(:) = MIN(0.018,MAX(0.011,0.011+(0.007/8.0)*(ZDDU(:)-10.0)))
736ENDIF
737
738!       7.6. Roughness lengths Z0 and Z0H over sea
739
740!IF (KZ0 == 0) THEN      ! ARPEGE formulation
741!  PZ0SEA (:) = (ZCHARN(:)/XG)*ZUSTAR2(:) + XVZ0CM*PCD(:)/PCDN(:)
742!  PZ0HSEA(:) = PZ0SEA (:)
743!ELSEIF (KZ0 == 1) THEN  ! Smith (1988) formulation
744!  PZ0SEA (:) = (ZCHARN(:)/XG)*ZUSTAR2(:) + 0.11*ZVISA(:)/PUSTAR(:)
745!  PZ0HSEA(:) = PZ0SEA (:)
746!ELSEIF (KZ0 == 2) THEN  ! Direct computation using the stability functions
747!  DO JLON=1,SIZE(PTA)
748!    PZ0SEA (JLON) = PUREF(JLON)/EXP(XKARMAN*ZDDU(JLON)/PUSTAR(JLON)+ZPSIU(JLON))
749!    Z0TSEA        = PZREF(JLON)/EXP(XKARMAN*ZDDT(JLON)/ZTSR  (JLON)+ZPSIT(JLON))
750!    Z0QSEA        = PZREF(JLON)/EXP(XKARMAN*ZDDQ(JLON)/ZQSR  (JLON)+ZPSIT(JLON))
751!    PZ0HSEA(JLON) = 0.5*(Z0TSEA+Z0QSEA)
752!  ENDDO
753!ENDIF
754
755WRITE(*,*) "JLON ",JLON
756WRITE(*,*) "PTA ",klon,PTA
757WRITE(*,*) "PCD ",SIZE(PCD),PCD
758WRITE(*,*) "PCQ ",SIZE(PCE),PCE
759WRITE(*,*) "PCH ",SIZE(PCH),PCH
760
761coeffs = [PCD,&
762       PCE,&
763       PCH]
764
765
766!IF (LHOOK) CALL DR_HOOK('ECUMEV6_FLUX',1,ZHOOK_HANDLE)
767
768!-------------------------------------------------------------------------------
769   END SUBROUTINE ECUMEV6_FLUX
Note: See TracBrowser for help on using the repository browser.