SURFEX  V8_0
Surfex V8_0 release
 All Classes Files Functions Variables
ecumev6_flux.F90
Go to the documentation of this file.
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,PEXNA,PRHOA,PSST,PSSS,PEXNS,PQA,PVMOD, &
7  pzref,puref,pps,ppa,pichce,oprecip,opwebb, &
8  pqsat,psfth,psftq,pustar,pcd,pcdn,pch,pce, &
9  pri,presa,prain,kz0,pz0hsea,opertflux,ppertflux )
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 !
86 USE modd_csts, ONLY : xpi, xday, xkarman, xg, xp00, xstefan, xrd, xrv, &
87  xcpd, xcpv, xcl, xtt, xlvtt
88 USE modd_surf_par, ONLY : xundef
89 USE modd_surf_atm, ONLY : xvchrnk, xvz0cm
90 USE modd_reprod_oper, ONLY : ccharnock
91 !
92 USE mode_thermos
93 USE modi_wind_threshold
94 USE modi_surface_ri
95 !
96 USE yomhook, ONLY : lhook, dr_hook
97 USE parkind1, ONLY : jprb
98 !
99 USE modi_abor1_sfx
100 !
101 IMPLICIT NONE
102 !
103 ! 0.1. Declarations of arguments
104 !
105 REAL, DIMENSION(:), INTENT(IN) :: pvmod ! module of wind at atm level (m/s)
106 REAL, DIMENSION(:), INTENT(IN) :: pta ! air temperature at atm level (K)
107 REAL, DIMENSION(:), INTENT(IN) :: pqa ! air spec. hum. at atm level (kg/kg)
108 REAL, DIMENSION(:), INTENT(IN) :: ppa ! air pressure at atm level (Pa)
109 REAL, DIMENSION(:), INTENT(IN) :: prhoa ! air density at atm level (kg/m3)
110 REAL, DIMENSION(:), INTENT(IN) :: pexna ! Exner function at atm level
111 REAL, DIMENSION(:), INTENT(IN) :: puref ! atm level for wind (m)
112 REAL, DIMENSION(:), INTENT(IN) :: pzref ! atm level for temp./hum. (m)
113 REAL, DIMENSION(:), INTENT(IN) :: psss ! Sea Surface Salinity (g/kg)
114 REAL, DIMENSION(:), INTENT(IN) :: pps ! air pressure at sea surface (Pa)
115 REAL, DIMENSION(:), INTENT(IN) :: pexns ! Exner function at sea surface
116 REAL, DIMENSION(:), INTENT(IN) :: ppertflux ! stochastic flux perturbation pattern
117 ! for correction
118 REAL, INTENT(IN) :: pichce !
119 LOGICAL, INTENT(IN) :: oprecip !
120 LOGICAL, INTENT(IN) :: opwebb !
121 LOGICAL, INTENT(IN) :: opertflux
122 REAL, DIMENSION(:), INTENT(IN) :: prain ! precipitation rate (kg/s/m2)
123 !
124 INTEGER, INTENT(IN) :: kz0
125 !
126 REAL, DIMENSION(:), INTENT(INOUT) :: psst ! Sea Surface Temperature (K)
127 REAL, DIMENSION(:), INTENT(INOUT) :: pz0sea ! roughness length over sea
128 REAL, DIMENSION(:), INTENT(OUT) :: pz0hsea ! heat roughness length over sea
129 
130 ! surface fluxes : latent heat, sensible heat, friction fluxes
131 REAL, DIMENSION(:), INTENT(OUT) :: pustar ! friction velocity (m/s)
132 REAL, DIMENSION(:), INTENT(OUT) :: psfth ! heat flux (W/m2)
133 REAL, DIMENSION(:), INTENT(OUT) :: psftq ! water flux (kg/m2/s)
134 
135 ! diagnostics
136 REAL, DIMENSION(:), INTENT(OUT) :: pqsat ! sea surface spec. hum. (kg/kg)
137 REAL, DIMENSION(:), INTENT(OUT) :: pcd ! transfer coef. for momentum
138 REAL, DIMENSION(:), INTENT(OUT) :: pch ! transfer coef. for temperature
139 REAL, DIMENSION(:), INTENT(OUT) :: pce ! transfer coef. for humidity
140 REAL, DIMENSION(:), INTENT(OUT) :: pcdn ! neutral coef. for momentum
141 REAL, DIMENSION(:), INTENT(OUT) :: pri ! Richardson number
142 REAL, DIMENSION(:), INTENT(OUT) :: presa ! aerodynamical resistance
143 !
144 ! 0.2. Declarations of local variables
145 !
146 ! specif SB
147 INTEGER, DIMENSION(SIZE(PTA)) :: jcv ! convergence index
148 INTEGER, DIMENSION(SIZE(PTA)) :: jiter ! nb of iterations to converge
149 
150 REAL, DIMENSION(SIZE(PTA)) :: ztau ! momentum flux (N/m2)
151 REAL, DIMENSION(SIZE(PTA)) :: zhf ! sensible heat flux (W/m2)
152 REAL, DIMENSION(SIZE(PTA)) :: zef ! latent heat flux (W/m2)
153 REAL, DIMENSION(SIZE(PTA)) :: ztaur ! momentum flx due to rain (N/m2)
154 REAL, DIMENSION(SIZE(PTA)) :: zrf ! sensible flx due to rain (W/m2)
155 REAL, DIMENSION(SIZE(PTA)) :: zefwebb ! Webb corr. on latent flx (W/m2)
156 
157 REAL, DIMENSION(SIZE(PTA)) :: zvmod ! wind intensity at atm level (m/s)
158 REAL, DIMENSION(SIZE(PTA)) :: zqsata ! sat.spec.hum. at atm level (kg/kg)
159 REAL, DIMENSION(SIZE(PTA)) :: zlva ! vap.heat of pure water at atm level (J/kg)
160 REAL, DIMENSION(SIZE(PTA)) :: zlvs ! vap.heat of seawater at sea surface (J/kg)
161 REAL, DIMENSION(SIZE(PTA)) :: zcpa ! specif.heat moist air (J/kg/K)
162 REAL, DIMENSION(SIZE(PTA)) :: zvisa ! kinemat.visc. of dry air (m2/s)
163 REAL, DIMENSION(SIZE(PTA)) :: zdu ! U vert.grad. (real atm)
164 REAL, DIMENSION(SIZE(PTA)) :: zdt,zdq ! T,Q vert.grad. (real atm)
165 REAL, DIMENSION(SIZE(PTA)) :: zddu ! U vert.grad. (real atm + gust)
166 REAL, DIMENSION(SIZE(PTA)) :: zddt,zddq ! T,Q vert.grad. (real atm + WL/CS)
167 REAL, DIMENSION(SIZE(PTA)) :: zusr ! velocity scaling param. (m/s)
168  ! =friction velocity
169 REAL, DIMENSION(SIZE(PTA)) :: ztsr ! temperature scaling param. (K)
170 REAL, DIMENSION(SIZE(PTA)) :: zqsr ! humidity scaling param. (kg/kg)
171 REAL, DIMENSION(SIZE(PTA)) :: zdeltau10n,zdeltat10n,zdeltaq10n
172  ! U,T,Q vert.grad. (10m, neutral atm)
173 REAL, DIMENSION(SIZE(PTA)) :: zusr0,ztsr0,zqsr0 ! ITERATIVE PROCESS
174 REAL, DIMENSION(SIZE(PTA)) :: zdusto,zdtsto,zdqsto ! ITERATIVE PROCESS
175 REAL, DIMENSION(SIZE(PTA)) :: zpsiu,zpsit! PSI funct for U, T/Q (Z0 comp)
176 REAL, DIMENSION(SIZE(PTA)) :: zcharn ! Charnock parameter (Z0 comp)
177 
178 REAL, DIMENSION(SIZE(PTA)) :: zustar2 ! square of friction velocity
179 REAL, DIMENSION(SIZE(PTA)) :: zac ! aerodynamical conductance
180 REAL, DIMENSION(SIZE(PTA)) :: zdircoszw ! orography slope cosine
181  ! (=1 on water!)
182 REAL, DIMENSION(SIZE(PTA)) :: zparun,zpartn,zparqn ! neutral parameter for U,T,Q
183 REAL, DIMENSION(0:5) :: zcoefu,zcoeft,zcoefq
184 
185 ! local constants
186 LOGICAL :: opcvflx ! to force convergence
187 INTEGER :: nitermax ! nb of iterations to get free convergence
188 INTEGER :: nitersup ! nb of additional iterations if OPCVFLX=.TRUE.
189 INTEGER :: niterfl ! maximum number of iterations
190 REAL :: zetv,zrdsrv ! thermodynamic constants
191 REAL :: zsqr3
192 REAL :: zlmomin,zlmomax ! min/max value of Obukhovs stability param. z/l
193 REAL :: zbta,zgma ! parameters of the stability functions
194 REAL :: zdusr0,zdtsr0,zdqsr0 ! maximum gap for USR/TSR/QSR between 2 steps
195 REAL :: zp00 ! [OPRECIP] - water vap. diffusiv.ref.press.(Pa)
196 REAL :: zutu,zutt,zutq ! U10n threshold in ECUME parameterisation
197 REAL :: zcdiru,zcdirt,zcdirq ! coef directeur pour fonction affine U,T,Q
198 REAL :: zordou,zordot,zordoq ! ordonnee a l'origine pour fonction affine U,T,Q
199 
200 INTEGER :: jj ! for ITERATIVE PROCESS
201 INTEGER :: jlon,jk
202 REAL :: zlmou,zlmot ! Obukhovs param. z/l for U, T/Q
203 REAL :: zpsi_u,zpsi_t ! PSI funct. for U, T/Q
204 REAL :: z0tsea,z0qsea ! roughness length for T, Q
205 REAL :: zchic,zchik,zpsic,zpsik,zlogus10,zlogts10
206 REAL :: ztac,zcpwa,zdqsdt,zdwat,zdtmp,zbulb ! [OPRECIP]
207 REAL :: zww ! [OPWEBB]
208 
209 REAL(KIND=JPRB) :: zhook_handle
210 !
211 !-------------------------------------------------------------------------------
212 IF (lhook) CALL dr_hook('ECUMEV6_FLUX',0,zhook_handle)
213 !
214 zdusr0 = 1.e-06
215 zdtsr0 = 1.e-06
216 zdqsr0 = 1.e-09
217 !
218 nitermax = 5
219 nitersup = 5
220 opcvflx = .true.
221 !
222 niterfl = nitermax
223 IF (opcvflx) niterfl = nitermax+nitersup
224 !
225 zcoefu = (/ 1.00e-03, 3.66e-02, -1.92e-03, 2.32e-04, -7.02e-06, 6.40e-08 /)
226 zcoeft = (/ 5.36e-03, 2.90e-02, -1.24e-03, 4.50e-04, -2.06e-05, 0.0 /)
227 zcoefq = (/ 1.00e-03, 3.59e-02, -2.87e-04, 0.0, 0.0, 0.0 /)
228 !
229 zutu = 40.0
230 zutt = 14.4
231 zutq = 10.0
232 !
233 zcdiru = zcoefu(1) + 2.0*zcoefu(2)*zutu + 3.0*zcoefu(3)*zutu**2 &
234  + 4.0*zcoefu(4)*zutu**3 + 5.0*zcoefu(5)*zutu**4
235 zcdirt = zcoeft(1) + 2.0*zcoeft(2)*zutt + 3.0*zcoeft(3)*zutt**2 &
236  + 4.0*zcoeft(4)*zutt**3
237 zcdirq = zcoefq(1) + 2.0*zcoefq(2)*zutq
238 !
239 zordou = zcoefu(0) + zcoefu(1)*zutu + zcoefu(2)*zutu**2 + zcoefu(3)*zutu**3 &
240  + zcoefu(4)*zutu**4 + zcoefu(5)*zutu**5
241 zordot = zcoeft(0) + zcoeft(1)*zutt + zcoeft(2)*zutt**2 + zcoeft(3)*zutt**3 &
242  + zcoeft(4)*zutt**4
243 zordoq = zcoefq(0) + zcoefq(1)*zutq + zcoefq(2)*zutq**2
244 !
245 !-------------------------------------------------------------------------------
246 !
247 ! 1. AUXILIARY CONSTANTS & ARRAY INITIALISATION BY UNDEFINED VALUES.
248 ! --------------------------------------------------------------------
249 !
250 zdircoszw(:) = 1.0
251 !
252 zetv = xrv/xrd-1.0 !~0.61 (cf Liu et al. 1979)
253 zrdsrv = xrd/xrv !~0.622
254 zsqr3 = sqrt(3.0)
255 zlmomin = -200.0
256 zlmomax = 0.25
257 zbta = 16.0
258 zgma = 7.0 !initially =4.7, modified to 7.0 following G. Caniaux
259 !
260 zp00 = 1013.25e+02
261 !
262 pcd(:) = xundef
263 pch(:) = xundef
264 pce(:) = xundef
265 pcdn(:) = xundef
266 zusr(:) = xundef
267 ztsr(:) = xundef
268 zqsr(:) = xundef
269 ztau(:) = xundef
270 zhf(:) = xundef
271 zef(:) = xundef
272 !
273 psfth(:) = xundef
274 psftq(:) = xundef
275 pustar(:) = xundef
276 presa(:) = xundef
277 pri(:) = xundef
278 !
279 ztaur(:) = 0.0
280 zrf(:) = 0.0
281 zefwebb(:) = 0.0
282 !
283 !-------------------------------------------------------------------------------
284 !
285 ! 2. INITIALISATIONS BEFORE ITERATIVE LOOP.
286 ! -------------------------------------------
287 !
288 zvmod(:) = wind_threshold(pvmod(:),puref(:)) !set a minimum value to wind
289 !
290 ! 2.0. Radiative fluxes - For warm layer & cool skin
291 !
292 ! 2.0b. Warm Layer correction
293 !
294 ! 2.1. Specific humidity at saturation
295 !
296 WHERE(psss(:)>0.0.AND.psss(:)/=xundef)
297  pqsat(:) = qsat_seawater2(psst(:),pps(:),psss(:)) !at sea surface
298 ELSEWHERE
299  pqsat(:) = qsat_seawater(psst(:),pps(:)) !at sea surface
300 ENDWHERE
301 zqsata(:) = qsat(pta(:),ppa(:)) !at atm level
302 !
303 ! 2.2. Gradients at the air-sea interface
304 !
305 zdu(:) = zvmod(:) !one assumes u is measured / sea surface current
306 zdt(:) = pta(:)/pexna(:)-psst(:)/pexns(:)
307 zdq(:) = pqa(:)-pqsat(:)
308 !
309 ! 2.3. Latent heat of vaporisation
310 !
311 zlva(:) = xlvtt+(xcpv-xcl)*(pta(:)-xtt) !of pure water at atm level
312 zlvs(:) = xlvtt+(xcpv-xcl)*(psst(:)-xtt) !of pure water at sea surface
313 WHERE(psss(:)>0.0.AND.psss(:)/=xundef)
314  zlvs(:) = zlvs(:)*(1.0-1.00472e-3*psss(:)) !of seawater at sea surface
315 ENDWHERE
316 !
317 ! 2.4. Specific heat of moist air (Businger 1982)
318 !
319 !ZCPA(:) = XCPD*(1.0+(XCPV/XCPD-1.0)*PQA(:))
320 zcpa(:) = xcpd
321 !
322 ! 2.4b Kinematic viscosity of dry air (Andreas 1989, CRREL Rep. 89-11)
323 !
324 zvisa(:) = 1.326e-05*(1.0+6.542e-03*(pta(:)-xtt)+8.301e-06*(pta(:)-xtt)**2 &
325  -4.84e-09*(pta(:)-xtt)**3)
326 !
327 ! 2.4c Coefficients for warm layer and/or cool skin correction
328 !
329 ! 2.5. Initial guess
330 !
331 zddu(:) = zdu(:)
332 zddt(:) = zdt(:)
333 zddq(:) = zdq(:)
334 zddu(:) = sign(max(abs(zddu(:)),10.0*zdusr0),zddu(:))
335 zddt(:) = sign(max(abs(zddt(:)),10.0*zdtsr0),zddt(:))
336 zddq(:) = sign(max(abs(zddq(:)),10.0*zdqsr0),zddq(:))
337 !
338 jcv(:) = -1
339 zusr(:) = 0.04*zddu(:)
340 ztsr(:) = 0.04*zddt(:)
341 zqsr(:) = 0.04*zddq(:)
342 zdeltau10n(:) = zddu(:)
343 zdeltat10n(:) = zddt(:)
344 zdeltaq10n(:) = zddq(:)
345 jiter(:) = 99
346 !
347 ! In the following, we suppose that Richardson number PRI < XRIMAX
348 ! If not true, Monin-Obukhov theory can't (and therefore shouldn't) be applied !
349 !-------------------------------------------------------------------------------
350 !
351 ! 3. ITERATIVE LOOP TO COMPUTE U*, T*, Q*.
352 ! ------------------------------------------
353 !
354 DO jj=1,niterfl
355  DO jlon=1,SIZE(pta)
356 !
357  IF (jcv(jlon) == -1) THEN
358  zusr0(jlon)=zusr(jlon)
359  ztsr0(jlon)=ztsr(jlon)
360  zqsr0(jlon)=zqsr(jlon)
361  IF (jj == nitermax+1 .OR. jj == nitermax+nitersup) THEN
362  zdeltau10n(jlon) = 0.5*(zdusto(jlon)+zdeltau10n(jlon)) !forced convergence
363  zdeltat10n(jlon) = 0.5*(zdtsto(jlon)+zdeltat10n(jlon))
364  zdeltaq10n(jlon) = 0.5*(zdqsto(jlon)+zdeltaq10n(jlon))
365  IF (jj == nitermax+nitersup) jcv(jlon)=3
366  ENDIF
367  zdusto(jlon) = zdeltau10n(jlon)
368  zdtsto(jlon) = zdeltat10n(jlon)
369  zdqsto(jlon) = zdeltaq10n(jlon)
370 !
371 ! 3.1. Neutral parameter for wind speed (ECUME_V6 formulation)
372 !
373  IF (zdeltau10n(jlon) <= zutu) THEN
374  zparun(jlon) = zcoefu(0) + zcoefu(1)*zdeltau10n(jlon) &
375  + zcoefu(2)*zdeltau10n(jlon)**2 &
376  + zcoefu(3)*zdeltau10n(jlon)**3 &
377  + zcoefu(4)*zdeltau10n(jlon)**4 &
378  + zcoefu(5)*zdeltau10n(jlon)**5
379  ELSE
380  zparun(jlon) = zcdiru*(zdeltau10n(jlon)-zutu) + zordou
381  ENDIF
382  pcdn(jlon) = (zparun(jlon)/zdeltau10n(jlon))**2
383 !
384 ! 3.2. Neutral parameter for temperature (ECUME_V6 formulation)
385 !
386  IF (zdeltau10n(jlon) <= zutt) THEN
387  zpartn(jlon) = zcoeft(0) + zcoeft(1)*zdeltau10n(jlon) &
388  + zcoeft(2)*zdeltau10n(jlon)**2 &
389  + zcoeft(3)*zdeltau10n(jlon)**3 &
390  + zcoeft(4)*zdeltau10n(jlon)**4
391  ELSE
392  zpartn(jlon) = zcdirt*(zdeltau10n(jlon)-zutt) + zordot
393  ENDIF
394 !
395 ! 3.3. Neutral parameter for humidity (ECUME_V6 formulation)
396 !
397  IF (zdeltau10n(jlon) <= zutq) THEN
398  zparqn(jlon) = zcoefq(0) + zcoefq(1)*zdeltau10n(jlon) &
399  + zcoefq(2)*zdeltau10n(jlon)**2
400  ELSE
401  zparqn(jlon) = zcdirq*(zdeltau10n(jlon)-zutq) + zordoq
402  ENDIF
403 !
404 ! 3.4. Scaling parameters U*, T*, Q*
405 !
406  zusr(jlon) = zparun(jlon)
407  ztsr(jlon) = zpartn(jlon)*zdeltat10n(jlon)/zdeltau10n(jlon)
408  zqsr(jlon) = zparqn(jlon)*zdeltaq10n(jlon)/zdeltau10n(jlon)
409 !
410 ! 3.4b Gustiness factor (Deardorff 1970)
411 !
412 ! 3.4c Cool skin correction
413 !
414 ! 3.5. Obukhovs stability param. z/l following Liu et al. (JAS, 1979)
415 !
416 ! For U
417  zlmou = puref(jlon)*xg*xkarman*(ztsr(jlon)/pta(jlon) &
418  +zetv*zqsr(jlon)/(1.0+zetv*pqa(jlon)))/zusr(jlon)**2
419 ! For T/Q
420  zlmot = zlmou*(pzref(jlon)/puref(jlon))
421  zlmou = max(min(zlmou,zlmomax),zlmomin)
422  zlmot = max(min(zlmot,zlmomax),zlmomin)
423 !
424 ! 3.6. Stability function psi (see Liu et al, 1979 ; Dyer and Hicks, 1970)
425 ! Modified to include convective form following Fairall (unpublished)
426 !
427 ! For U
428  IF (zlmou == 0.0) THEN
429  zpsi_u = 0.0
430  ELSEIF (zlmou > 0.0) THEN
431  zpsi_u = -zgma*zlmou
432  ELSE
433  zchik = (1.0-zbta*zlmou)**0.25
434  zpsik = 2.0*log((1.0+zchik)/2.0) &
435  +log((1.0+zchik**2)/2.0) &
436  -2.0*atan(zchik)+0.5*xpi
437  zchic = (1.0-12.87*zlmou)**(1.0/3.0) !for very unstable conditions
438  zpsic = 1.5*log((zchic**2+zchic+1.0)/3.0) &
439  -zsqr3*atan((2.0*zchic+1.0)/zsqr3) &
440  +xpi/zsqr3
441  zpsi_u = zpsic+(zpsik-zpsic)/(1.0+zlmou**2) !match Kansas & free-conv. forms
442  ENDIF
443  zpsiu(jlon) = zpsi_u
444 ! For T/Q
445  IF (zlmot == 0.0) THEN
446  zpsi_t = 0.0
447  ELSEIF (zlmot > 0.0) THEN
448  zpsi_t = -zgma*zlmot
449  ELSE
450  zchik = (1.0-zbta*zlmot)**0.25
451  zpsik = 2.0*log((1.0+zchik**2)/2.0)
452  zchic = (1.0-12.87*zlmot)**(1.0/3.0) !for very unstable conditions
453  zpsic = 1.5*log((zchic**2+zchic+1.0)/3.0) &
454  -zsqr3*atan((2.0*zchic+1.0)/zsqr3) &
455  +xpi/zsqr3
456  zpsi_t = zpsic+(zpsik-zpsic)/(1.0+zlmot**2) !match Kansas & free-conv. forms
457  ENDIF
458  zpsit(jlon) = zpsi_t
459 !
460 ! 3.7. Update air-sea gradients
461 !
462  zddu(jlon) = zdu(jlon)
463  zddt(jlon) = zdt(jlon)
464  zddq(jlon) = zdq(jlon)
465  zddu(jlon) = sign(max(abs(zddu(jlon)),10.0*zdusr0),zddu(jlon))
466  zddt(jlon) = sign(max(abs(zddt(jlon)),10.0*zdtsr0),zddt(jlon))
467  zddq(jlon) = sign(max(abs(zddq(jlon)),10.0*zdqsr0),zddq(jlon))
468  zlogus10 = log(puref(jlon)/10.0)
469  zlogts10 = log(pzref(jlon)/10.0)
470  zdeltau10n(jlon) = zddu(jlon)-zusr(jlon)*(zlogus10-zpsi_u)/xkarman
471  zdeltat10n(jlon) = zddt(jlon)-ztsr(jlon)*(zlogts10-zpsi_t)/xkarman
472  zdeltaq10n(jlon) = zddq(jlon)-zqsr(jlon)*(zlogts10-zpsi_t)/xkarman
473  zdeltau10n(jlon) = sign(max(abs(zdeltau10n(jlon)),10.0*zdusr0), &
474  zdeltau10n(jlon))
475  zdeltat10n(jlon) = sign(max(abs(zdeltat10n(jlon)),10.0*zdtsr0), &
476  zdeltat10n(jlon))
477  zdeltaq10n(jlon) = sign(max(abs(zdeltaq10n(jlon)),10.0*zdqsr0), &
478  zdeltaq10n(jlon))
479 !
480 ! 3.8. Test convergence for U*, T*, Q*
481 !
482  IF (abs(zusr(jlon)-zusr0(jlon)) < zdusr0 .AND. &
483  abs(ztsr(jlon)-ztsr0(jlon)) < zdtsr0 .AND. &
484  abs(zqsr(jlon)-zqsr0(jlon)) < zdqsr0) THEN
485  jcv(jlon) = 1 !free convergence
486  IF (jj >= nitermax+1) jcv(jlon) = 2 !leaded convergence
487  ENDIF
488  jiter(jlon) = jj
489  ENDIF
490 !
491  ENDDO
492 ENDDO
493 !
494 !-------------------------------------------------------------------------------
495 !
496 ! 4. COMPUTATION OF TURBULENT FLUXES AND EXCHANGE COEFFICIENTS.
497 ! ---------------------------------------------------------------
498 !
499 DO jlon=1,SIZE(pta)
500 !
501 ! 4.1. Surface turbulent fluxes
502 ! (ATM CONV.: ZTAU<<0 ; ZHF,ZEF<0 if atm looses heat)
503 !
504  ztau(jlon) = -prhoa(jlon)*zusr(jlon)**2
505  zhf(jlon) = -prhoa(jlon)*zcpa(jlon)*zusr(jlon)*ztsr(jlon)
506  zef(jlon) = -prhoa(jlon)*zlvs(jlon)*zusr(jlon)*zqsr(jlon)
507 !
508 ! 4.2. Exchange coefficients PCD, PCH, PCE
509 !
510  pcd(jlon) = (zusr(jlon)/zddu(jlon))**2
511  pch(jlon) = (zusr(jlon)*ztsr(jlon))/(zddu(jlon)*zddt(jlon))
512  pce(jlon) = (zusr(jlon)*zqsr(jlon))/(zddu(jlon)*zddq(jlon))
513 !
514 ! 4.3. Stochastic perturbation of turbulent fluxes
515 !
516  IF( opertflux )THEN
517  ztau(jlon) = ztau(jlon)* ( 1. + ppertflux(jlon) / 2. )
518  zhf(jlon) = zhf(jlon)* ( 1. + ppertflux(jlon) / 2. )
519  zef(jlon) = zef(jlon)* ( 1. + ppertflux(jlon) / 2. )
520  ENDIF
521 !
522 ENDDO
523 !
524 !-------------------------------------------------------------------------------
525 !
526 ! 5. COMPUTATION OF FLUX CORRECTIONS DUE TO RAINFALL.
527 ! (ATM conv: ZRF<0 if atm. looses heat, ZTAUR<<0)
528 ! -----------------------------------------------------
529 !
530 IF (oprecip) THEN
531  DO jlon=1,SIZE(pta)
532 !
533 ! 5.1. Momentum flux due to rainfall (ZTAUR, N/m2)
534 !
535 ! See pp3752 in FBR96.
536  ztaur(jlon) = -0.85*prain(jlon)*pvmod(jlon)
537 !
538 ! 5.2. Sensible heat flux due to rainfall (ZRF, W/m2)
539 !
540 ! See Eq.12 in GoF95 with ZCPWA as specific heat of water at atm level (J/kg/K),
541 ! ZDQSDT from Clausius-Clapeyron relation, ZDWAT as water vapor diffusivity
542 ! (Eq.13-3 of Pruppacher and Klett, 1978), ZDTMP as heat diffusivity, and ZBULB
543 ! as wet-bulb factor (Eq.11 in GoF95).
544 !
545  ztac = pta(jlon)-xtt
546  zcpwa = 4217.51 -3.65566*ztac +0.1381*ztac**2 &
547  -2.8309e-03*ztac**3 +3.42061e-05*ztac**4 &
548  -2.18107e-07*ztac**5 +5.74535e-10*ztac**6
549  zdqsdt = (zlva(jlon)*zqsata(jlon))/(xrv*pta(jlon)**2)
550  zdwat = 2.11e-05*(zp00/ppa(jlon))*(pta(jlon)/xtt)**1.94
551  zdtmp = (1.0+3.309e-03*ztac-1.44e-06*ztac**2) &
552  *0.02411/(prhoa(jlon)*zcpa(jlon))
553  zbulb = 1.0/(1.0+zdqsdt*(zlva(jlon)*zdwat)/(zcpa(jlon)*zdtmp))
554  zrf(jlon) = prain(jlon)*zcpwa*zbulb*((psst(jlon)-pta(jlon)) &
555  +(pqsat(jlon)-pqa(jlon))*(zlva(jlon)*zdwat)/(zcpa(jlon)*zdtmp))
556 !
557  ENDDO
558 ENDIF
559 !
560 !-------------------------------------------------------------------------------
561 !
562 ! 6. COMPUTATION OF WEBB CORRECTION TO LATENT HEAT FLUX (ZEFWEBB, W/m2).
563 ! ------------------------------------------------------------------------
564 !
565 ! See Eq.21 and Eq.22 in FBR96.
566 IF (opwebb) THEN
567  DO jlon=1,SIZE(pta)
568  zww = (1.0+zetv)*(zusr(jlon)*zqsr(jlon)) &
569  +(1.0+(1.0+zetv)*pqa(jlon))*(zusr(jlon)*ztsr(jlon))/pta(jlon)
570  zefwebb(jlon) = -prhoa(jlon)*zlvs(jlon)*zww*pqa(jlon)
571  ENDDO
572 ENDIF
573 !
574 !-------------------------------------------------------------------------------
575 !
576 ! 7. FINAL STEP : TOTAL SURFACE FLUXES AND DERIVED DIAGNOSTICS.
577 ! ---------------------------------------------------------------
578 !
579 ! 7.1. Richardson number
580 !
581  CALL surface_ri(psst,pqsat,pexns,pexna,pta,pqa, &
582  pzref,puref,zdircoszw,pvmod,pri)
583 !
584 ! 7.2. Friction velocity which contains correction due to rain
585 !
586 zustar2(:) = -(ztau(:)+ztaur(:))/prhoa(:)
587 !
588 IF (oprecip) THEN
589  pcd(:) = zustar2(:)/zddu(:)**2
590 ENDIF
591 !
592 pustar(:) = sqrt(zustar2(:))
593 !
594 ! 7.3. Aerodynamical conductance and resistance
595 !
596 zac(:) = pch(:)*zddu(:)
597 presa(:) = 1.0/zac(:)
598 !
599 ! 7.4. Total surface fluxes
600 !
601 psfth(:) = zhf(:)+zrf(:)
602 psftq(:) = (zef(:)+zefwebb(:))/zlvs(:)
603 !
604 ! 7.5. Charnock number
605 !
606 IF (ccharnock == 'OLD') THEN
607  zcharn(:) = xvchrnk
608 ELSE !modified for moderate wind speed as in COARE3.0
609  zcharn(:) = min(0.018,max(0.011,0.011+(0.007/8.0)*(zddu(:)-10.0)))
610 ENDIF
611 !
612 ! 7.6. Roughness lengths Z0 and Z0H over sea
613 !
614 IF (kz0 == 0) THEN ! ARPEGE formulation
615  pz0sea(:) = (zcharn(:)/xg)*zustar2(:) + xvz0cm*pcd(:)/pcdn(:)
616  pz0hsea(:) = pz0sea(:)
617 ELSEIF (kz0 == 1) THEN ! Smith (1988) formulation
618  pz0sea(:) = (zcharn(:)/xg)*zustar2(:) + 0.11*zvisa(:)/pustar(:)
619  pz0hsea(:) = pz0sea(:)
620 ELSEIF (kz0 == 2) THEN ! Direct computation using the stability functions
621  DO jlon=1,SIZE(pta)
622  pz0sea(jlon) = puref(jlon)/exp(xkarman*zddu(jlon)/pustar(jlon)+zpsiu(jlon))
623  z0tsea = pzref(jlon)/exp(xkarman*zddt(jlon)/ztsr(jlon)+zpsit(jlon))
624  z0qsea = pzref(jlon)/exp(xkarman*zddq(jlon)/zqsr(jlon)+zpsit(jlon))
625  pz0hsea(jlon) = 0.5*(z0tsea+z0qsea)
626  ENDDO
627 ENDIF
628 !
629 !
630 IF (lhook) CALL dr_hook('ECUMEV6_FLUX',1,zhook_handle)
631 !
632 !-------------------------------------------------------------------------------
633  END SUBROUTINE ecumev6_flux
mode_thermos::qsat_seawater
Definition: mode_thermos.F90:62
wind_threshold
real function, dimension(size(pwind)) wind_threshold(PWIND, PUREF)
Definition: wind_threshold.F90:6
surface_ri
subroutine surface_ri(PTG, PQS, PEXNS, PEXNA, PTA, PQA, PZREF, PUREF, PDIRCOSZW, PVMOD, PRI)
Definition: surface_ri.F90:6
mode_thermos
Definition: mode_thermos.F90:6
ecumev6_flux
subroutine ecumev6_flux(PZ0SEA, PTA, PEXNA, PRHOA, PSST, PSSS, PEXNS, PQA, PVMOD, PZREF, PUREF, PPS, PPA, PICHCE, OPRECIP, OPWEBB, PQSAT, PSFTH, PSFTQ, PUSTAR, PCD, PCDN, PCH, PCE, PRI, PRESA, PRAIN, KZ0, PZ0HSEA, OPERTFLUX, PPERTFLUX)
Definition: ecumev6_flux.F90:6
modd_surf_atm
Definition: modd_surf_atm.F90:6
mode_thermos::qsat_seawater2
Definition: mode_thermos.F90:65
mode_thermos::qsat
Definition: mode_thermos.F90:57
modd_surf_par
Definition: modd_surf_par.F90:6
modd_reprod_oper
Definition: modd_reprod_oper.F90:6
modd_csts
Definition: modd_csts.F90:6
Generated on Fri May 13 2016 09:31:04 for SURFEX by   doxygen 1.8.5