The TELEMAC-MASCARET system  trunk
exchange_with_atmosphere.f
Go to the documentation of this file.
1 ! *******************************
3 ! *******************************
4 !
5 !***********************************************************************
6 ! WAQTEL V8P2
7 !***********************************************************************
8 !
9 !brief Module containing some subroutines to deal with heat exchange
10 !+ with atmosphere
11 !
12 !history N. DURAND, A. GINEAU (EDF-LNHE)
13 !+ MAY 2011
14 !+ V6P0
15 !+ SOLRAD, SHORTRAD, EVAPO SUBROUTINES
16 !+ LEAP, DAYNUM FUNCTIONS FROM SOGREAH (NOW ARTELIA)
17 !
18 !history C.-T. PHAM (EDF-LNHE)
19 !+ 27/06/2012
20 !+ V6P2
21 !+ Creation of module EXCHANGE_WITH_ATMOSPHERE from previous
22 !+ subroutines
23 !
24 !history N. DURAND, N. LORRAIN, C.-T. PHAM (EDF-LNHE)
25 !+ 03/04/2014
26 !+ V7P0
27 !+ Update after Nicolas Lorrain's internship
28 !+
29 !
30 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
31 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
32 !
34  IMPLICIT NONE
35 !
36  PRIVATE
37  PUBLIC :: solrad,shortrad,evapo,leap,daynum
38 !
39 !-----------------------------------------------------------------------
40 !
41 ! brought to declarations_waqtel
42 ! CP: SPECIFIC HEAT OF WATER AT CONSTANT PRESSURE
43 ! DOUBLE PRECISION, PARAMETER :: CP = 4.18D3
44 !
45 !-----------------------------------------------------------------------
46 !
47  CONTAINS
48 !
49 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
50 !
51 ! *****************
52  SUBROUTINE solrad
53 ! *****************
54 !
55  &(ray_sol,nebu,mardat,martim,at,latitude,longitude)
56 !
57 !***********************************************************************
58 ! TELEMAC-3D V7P0 22/06/2012
59 !***********************************************************************
60 !
61 !brief EVALUATES SOLAR RADIATION INCIDENT ON THE SEA SURFACE
62 !+ - CALCULATES SOLAR RADIATION AS FUNCTION OF DAY NUMBER
63 !+ OF THE YEAR AND GEOGRAPHICAL LOCATION
64 !+ - INCLUDES ATMOSPHERICAL ABSORPTION AND REFLECTION, CLOUD
65 !+ COVERAGE, SEA SURFACE ALBEDO
66 !+ - TIME EXPRESSED IN GMT
67 !+ SOURCES:
68 !+ - PERRIN DE BRICHAMBAUT (1975)
69 !+ - BERLIAND'S METHOD (1960)
70 !+ - COOPER'S FORMULA (1969)
71 !
72 !history N. DURAND, A. GINEAU (EDF-LNHE)
73 !+ MAY 2011
74 !+ V6P0
75 !+
76 !
77 !history C.-T. PHAM (EDF-LNHE)
78 !+ JUNE 2012
79 !+ V6P2
80 !+ Adding of MARDAT, MARTIM, LATITUDE, LONGITUDE, ALBEDO AND type of
81 !+ sky as new arguments + INTENT + cosmetics
82 !+
83 !
84 !history C.-T. PHAM (EDF-LNHE)
85 !+ 03/04/2014
86 !+ V7P0
87 !+ Update after Nicolas Lorrain's internship
88 !+ Albedo is not constant in time and is no more an argument
89 !+ The type of sky is no more an argument, but fixed in the subroutine
90 !+ and may be changed by the user
91 !+
92 !
93 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
94 !| AT |-->| CURRENT TIME
95 !| NEBU |-->| NEBULOSITY (IN OCTAS)
96 !| LATITUDE |-->| LATITUDE
97 !| LONGITUDE |-->| LONGITUDE
98 !| MARDAT |-->| DATE (YEAR, MONTH,DAY)
99 !| MARTIM |-->| TIME (HOUR, MINUTE,SECOND)
100 !| RAY_SOL |<--| SOLAR RADIATION INCIDENT ON THE SEA SURFACE
101 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
102 !
103  USE declarations_waqtel,ONLY: iskytype
104  IMPLICIT NONE
105 !
106 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
107 !
108  INTEGER, INTENT(IN) :: MARDAT(3),MARTIM(3)
109  DOUBLE PRECISION, INTENT(IN) :: AT,NEBU
110  DOUBLE PRECISION, INTENT(IN) :: LATITUDE,LONGITUDE
111  DOUBLE PRECISION, INTENT(OUT) :: RAY_SOL
112 !
113 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
114 !
115  INTEGER IYEAR,IMONTH,IDAY,IHOUR,IMIN,ISEC
116 !
117  DOUBLE PRECISION DTR,PI,ALB
118  DOUBLE PRECISION DAY,DAYREEL,NDAYS
119 ! HA : SUN'S HOUR ANGLE [rad]
120 ! HR : TIME OF THE DAY IN HOURS (GMT)
121 ! RDEC: SUN'S DECLINATION [rad]
122 ! RLAT: LATITUDE [rad]
123 ! SING: SIN(GAMMA)
124 ! TE : TIME EQUATION [hours]
125  DOUBLE PRECISION HA,HR,RDEC,RLAT,SING,TE
126 ! AA,BB: COEFFICIENTS DEALING WITH LUMINOSITY AND SKY COLOUR
127  DOUBLE PRECISION AA,BB
128 !
129 ! INTEGER LEAP,DAYNUM
130 ! EXTERNAL LEAP,DAYNUM
131 !
132 !-----------------------------------------------------------------------
133 !
134 !
135  IF(iskytype.EQ.1) THEN
136 ! VERY PURE SKY
137  aa = 1130.d0
138  bb = 1.15d0
139  ELSEIF(iskytype.EQ.2) THEN
140 ! MEAN PURE SKY
141  aa = 1080.d0
142  bb = 1.22d0
143  ELSEIF(iskytype.EQ.3) THEN
144 ! INDUSTRIAL AREA
145  aa = 995.d0
146  bb = 1.25d0
147  ENDIF
148 !
149  iyear = mardat(1)
150  imonth = mardat(2)
151  iday = mardat(3)
152  ihour = martim(1)
153  imin = martim(2)
154  isec = martim(3)
155 !
156 !-----------------------------------------------------------------------
157 !
158  pi = 4.d0*atan(1.d0)
159  dtr = pi/180.d0
160 !
161 ! DAY NUMBER, ORBITAL CORRECTION
162  day = daynum(iyear,imonth,iday,ihour,imin,isec)
163  & + floor(at/86400.d0)
164  ndays = 365.d0 + REAL(leap(iyear))
165  dayreel = modulo(day, ndays)
166 
167 ! ALBEDO WITH RESPECT OF THE MONTH
168  IF(dayreel.GE.0.d0.AND.dayreel.LE.31.d0) THEN
169  alb = 0.11d0
170  ELSEIF(dayreel.GT.31.d0.AND.dayreel.LE.59.d0) THEN
171  alb = 0.10d0
172  ELSEIF(dayreel.GT.59.d0.AND.dayreel.LE.90.d0) THEN
173  alb = 0.08d0
174  ELSEIF(dayreel.GT.90.d0.AND.dayreel.LE.120.d0) THEN
175  alb = 0.07d0
176  ELSEIF(dayreel.GT.120.d0.AND.dayreel.LE.151.d0) THEN
177  alb = 0.06d0
178  ELSEIF(dayreel.GT.151.d0.AND.dayreel.LE.181.d0) THEN
179  alb = 0.06d0
180  ELSEIF(dayreel.GT.181.d0.AND.dayreel.LE.212.d0) THEN
181  alb = 0.06d0
182  ELSEIF(dayreel.GT.212.d0.AND.dayreel.LE.243.d0) THEN
183  alb = 0.07d0
184  ELSEIF(dayreel.GT.243.d0.AND.dayreel.LE.273.d0) THEN
185  alb = 0.07d0
186  ELSEIF(dayreel.GT.273.d0.AND.dayreel.LE.304.d0) THEN
187  alb = 0.08d0
188  ELSEIF(dayreel.GT.304.d0.AND.dayreel.LE.334.d0) THEN
189  alb = 0.11d0
190  ELSEIF(dayreel.GT.334.d0.AND.dayreel.LE.365.d0) THEN
191  alb = 0.12d0
192  ENDIF
193 
194 ! DECLINATION OF SUN (COOPER'S FORMULA)
195 ! RDEC = (23.45D0*SIN(2.D0*PI*(DAYREEL+284.D0)/NDAYS))*DTR
196  rdec = (23.45d0*cos(2.d0*pi*(172.d0-dayreel)/ndays))*dtr
197 
198 ! TIME EQUATION
199  te = ( 450.68d0*sin(2.d0*pi*dayreel/ndays-0.026903d0)
200  & +595.40d0*sin(4.d0*pi*dayreel/ndays+0.352835d0))/3600.d0
201 ! SOLAR ALTITUDE
202  hr = ihour+modulo(at,86400.d0)/3600.d0
203 !
204  rlat = latitude*dtr
205  ha = (hr-te-12.d0 + longitude/15.d0)*pi/12.d0
206  sing = sin(rlat)*sin(rdec) + cos(rlat)*cos(rdec)*cos(ha)
207 ! SOLAR RADIATION
208  IF(sing.LE.0.d0) THEN
209  ray_sol = 0.d0
210  ELSE
211 ! THE NEBULOSITY IS GIVEN IN OCTAS
212  ray_sol = aa*(sing**bb)*(1.d0-0.65d0*(nebu/8.d0)**2)
213  & *(1.d0-alb)
214  ENDIF
215 !
216 !-----------------------------------------------------------------------
217 !
218  RETURN
219  END SUBROUTINE solrad
220 
221 ! *********************
222  INTEGER FUNCTION leap
223 ! *********************
224 !
225  &(iyear)
226 !
227 !***********************************************************************
228 ! TELEMAC-3D V6P2 27/06/2012
229 !***********************************************************************
230 !
231 !brief DETERMINES WHETHER IYEAR IS A LEAP YEAR
232 !+ DESCRIPTION - RETURNS 1 IF IYEAR IS A LEAP YEAR, 0 OTHERWISE
233 !+
234 !
235 !history C. GUILBAUD (SOGREAH)
236 !+ JUNE 2001
237 !+ V6P0?
238 !+
239 !
240 !history C.-T. PHAM (EDF-LNHE)
241 !+ 27/06/2012
242 !+ V6P2
243 !+ Introduction into EXCHANGE_WITH_ATMOSPHERE module + INTENT
244 !+
245 !
246 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
247 !| IYEAR |-->| INDEX OF YEAR
248 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
249 !
250  IMPLICIT NONE
251 !
252 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
253 !
254  INTEGER, INTENT(IN) :: IYEAR
255 !
256 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
257 !
258  IF( mod(iyear,4).EQ.0.AND.
259  & (mod(iyear,100).NE.0.OR.mod(iyear,400).EQ.0)) THEN
260  leap = 1
261  ELSE
262  leap = 0
263  ENDIF
264 !
265 !-----------------------------------------------------------------------
266 !
267  RETURN
268  END FUNCTION leap
270 ! ********************************
271  DOUBLE PRECISION FUNCTION daynum
272 ! ********************************
273 !
274  &(iyear,imonth,iday,ihour,imin,isec)
275 !
276 !***********************************************************************
277 ! TELEMAC-3D V6P2 27/06/2012
278 !***********************************************************************
279 !
280 !brief RETURNS DAY NUMBER OF THE YEAR (FRACTIONAL)
281 !+
282 !
283 !history C. GUILBAUD (SOGREAH)
284 !+ JUNE 2001
285 !+ V6P0?
286 !+
287 !
288 !history C.-T. PHAM (EDF-LNHE)
289 !+ 27/06/2012
290 !+ V6P2
291 !+ Introduction into EXCHANGE_WITH_ATMOSPHERE module
292 !+ Change for type of result (double precision, not integer)
293 !+ + REAL conversion + addition of seconds ISEC
294 !+
295 !
296 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
297 !| IDAY |-->| INDEX OF DAY
298 !| IHOUR |-->| INDEX OF HOUR
299 !| IMIN |-->| INDEX OF MINUTE
300 !| IMONTH |-->| INDEX OF MONTH
301 !| ISEC |-->| INDEX OF SECOND
302 !| IYEAR |-->| INDEX OF YEAR
303 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
304 !
305  IMPLICIT NONE
306 !
307 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
308 !
309  INTEGER, INTENT(IN) :: IYEAR,IMONTH,IDAY,IHOUR,IMIN,ISEC
310 !
311 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
312 !
313 !
314 ! INTEGER LEAP
315 ! EXTERNAL LEAP
316 !
317  INTEGER MONTH(12)
318  parameter( month=(/0,31,59,90,120,151,181,212,243,273,304,334/) )
319 !
320 !-----------------------------------------------------------------------
321 !
322  daynum = REAL(month(imonth)+iday)
323  & + REAL(ihour)/24.d0+REAL(imin)/1440.d0+REAL(isec)/86400.d0
324  IF(imonth.GT.2) daynum = daynum + REAL(leap(iyear))
325 !
326 !-----------------------------------------------------------------------
327 !
328  RETURN
329  END FUNCTION daynum
330 
331 ! *******************
332  SUBROUTINE shortrad
333 ! *******************
334 !
335  &(treel,tair,nebu,hrel,ray_atm,ray_eau)
336 !
337 !***********************************************************************
338 ! TELEMAC-3D V7P0 22/06/2012
339 !***********************************************************************
340 !
341 !brief CALCULATES ATMOSPHERIC AND WATER RADIATIONS
342 !+ SOURCES:
343 !+ - SWINBANK'S METHOD
344 !+ - T.V.A. 1972
345 !+
346 !
347 !history N. DURAND, A. GINEAU (EDF-LNHE)
348 !+ MAY 2011
349 !+ V6P0
350 !+
351 !
352 !history C.-T. PHAM (EDF-LNHE)
353 !+ JUNE 2012
354 !+ V6P2
355 !+ Type of cloud taken into account + INTENT + cosmetics
356 !+
357 !
358 !history C.-T. PHAM (EDF-LNHE)
359 !+ 03/04/2014
360 !+ V7P0
361 !+ Update after Nicolas Lorrain's internship
362 !+ The type of cloud is no more an argument, but fixed in the subroutine
363 !+ and may be changed by the user
364 !+
365 !
366 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
367 !| NEBU |-->| NEBULOSITY (IN OCTAS)
368 !| RAY_ATM |<--| ATMOSPHERIC RADIATION
369 !| RAY_EAU |<--| WATER RADIATION
370 !| TAIR |-->| AIR TEMPERATURE
371 !| TREEL |-->| REAL WATER TEMPERATURE
372 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
373 !
376  IMPLICIT NONE
377 !
378 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
379 !
380  DOUBLE PRECISION, INTENT(IN) :: TREEL,TAIR,NEBU,HREL
381  DOUBLE PRECISION, INTENT(OUT) :: RAY_ATM,RAY_EAU
382 !
383 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
384 !
385 ! NUA: COEFFICIENT DEPENDING ON THE TYPE OF CLOUDS
386 ! REPLACED BY COEF_K (KEYWORD COEFFICIENT OF CLOUDING RATE)
387 ! ICLOUDTYPE FOR THE TYPE OF CLOUD IS THEN USELESS
388 !
389 ! TMP TO USE THE SAME FORMULAE AS IN GLM
390  DOUBLE PRECISION EPS_STAR
391  DOUBLE PRECISION QSAT_AIR
392 !
393 !-----------------------------------------------------------------------
394 !
395 ! ATMOSPHERIC RADIATION
396  IF(iray_atm.EQ.1) THEN
397 ! CASE 1: IDSO AND JACKSON (1969)
398  eps_star = (1.d0+0.275d0*(nebu/8.d0))
399  & *(1.d0-0.261d0*exp(-7.77d-8*tair**2))
400  ELSEIF(iray_atm.EQ.2) THEN
401 ! CASE 2: SWINBANK (1963)
402 ! DEFAULT, LIKE VERSIONS V7P1 AND BEFORE
403  eps_star = (1.d0+coef_k*(nebu/8.d0)**2)*0.937d-5
404  & *(tair+273.15d0)**2
405  ELSEIF(iray_atm.EQ.3) THEN
406 ! CASE 3: BRUTSAERT (1975)
407  qsat_air = 10.d0**(9.28603523d0-(2322.37885d0/(tair+273.15d0)))
408  eps_star = (1.d0+0.275d0*(nebu/8.d0))*0.642d0
409  & *((hrel/100.d0)*qsat_air/(tair+273.15d0))**(1.d0/7.d0)
410  ELSEIF(iray_atm.EQ.4) THEN
411 ! CASE 4: YAJIMA TONO DAM (2014)
412  qsat_air = 10.d0**(9.28603523d0-(2322.37885d0/(tair+273.15d0)))
413  eps_star = (1.d0-(nebu/8.d0)**2.796d0)*1.24d0
414  & *((hrel/100.d0)*qsat_air/(tair+273.15d0))**(1.d0/7.d0)
415  & + 0.955d0*(nebu/8.d0)**2.796d0
416  ELSE
417  WRITE(lu,*) 'FORMULA OF ATMOSPHERIC RADIATION NOT POSSIBLE'
418  CALL plante(1)
419  ENDIF
420 !
421  ray_atm = ema*eps_star*boltz*(tair+273.15d0)**4
422 !
423 ! WATER RADIATION
424  ray_eau = emi_eau*boltz*(treel+273.15d0)**4
425 !
426 !-----------------------------------------------------------------------
427 !
428  RETURN
429  END SUBROUTINE shortrad
430 
431 ! ****************
432  SUBROUTINE evapo
433 ! ****************
434 !
435  &(treel,tair,w2,patm,hrel,ro,flux_evap,flux_sens,debevap,c_atmos,
436  & c1_atmos,c2_atmos)
437 !
438 !***********************************************************************
439 ! WAQTEL V8P2
440 !***********************************************************************
441 !
442 !brief CALCULATES FLUX OF LATENT HEAT (W/M^2)
443 !+ CALCULATES SENSIBLE FLUX (W/M^2)
444 !+ CALCULAGES EVAPORATED WATER FLOWRATE (M/S)
445 !+ SOURCES:
446 !+ - BOLTON 1980 FOR SATURATION VAPOUR PRESSURE
447 !+
448 !
449 !history N. DURAND, A. GINEAU (EDF-LNHE)
450 !+ MAY 2011
451 !+ V6P0
452 !+
453 !
454 !history C.-T. PHAM (EDF-LNHE)
455 !+ JUNE 2012
456 !+ V6P2
457 !+ Parameter B to calibrate in arguments + cosmetics
458 !+
459 !
460 !history C.-T. PHAM (EDF-LNHE)
461 !+ 03/04/2014
462 !+ V7P0
463 !+ Update after Nicolas Lorrain's internship
464 !+ The salinity is not an argument anymore.
465 !+ Every salinity correction has been removed
466 !+
467 !
468 !history M. GANT, L. ABBAS, C.-T. PHAM (EDF-LNHE)
469 !+ 07/12/2015
470 !+ V7P1
471 !+ DEBEVAP and FLUX_EVAP have to be greater than or equal to zero.
472 !+ Otherwise, rain is read through files
473 !+
474 !
475 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
476 !| C_ATMOS |-->| PARAMETER TO CALIBRATE (FWW = C_ATMOS*(1.D0+W2))
477 !| C1_ATMOS |-->| PARAMETER TO CALIBRATE (FWW = C1_ATMOS+C2_ATMOS*W)
478 !| C2_ATMOS |-->| PARAMETER TO CALIBRATE (FWW = C1_ATMOS+C2_ATMOS*W)
479 !| DEB_EVAP |<--| EVAPORATION FLOWRATE AT THE SURFACE
480 !| FLUX_EVAP |<--| ENERGY FLUX DUE TO EVAPORATED WATER
481 !| FLUX_SENS |<--| HEAT FLUX BY CONVECTION
482 !| HREL |-->| RELATIVE HUMIDITY
483 !| PATM |-->| ATMOSPHERIC PRESSURE
484 !| RO |-->| DENSITY
485 !| TAIR |-->| AIR TEMPERATURE
486 !| TREEL |-->| REAL WATER TEMPERATURE
487 !| W2 |-->| RELATIVE MAGNITUDE OF WIND AT 2 M
488 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
489 !
491  IMPLICIT NONE
492 !
493 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
494 !
495  DOUBLE PRECISION, INTENT(IN) :: TREEL,TAIR,W2,PATM,HREL,RO
496  DOUBLE PRECISION, INTENT(IN) :: C_ATMOS,C1_ATMOS,C2_ATMOS
497  DOUBLE PRECISION, INTENT(OUT) :: FLUX_EVAP,FLUX_SENS,DEBEVAP
498 !
499 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
500 !
501  DOUBLE PRECISION Q_SAT_EAU,Q_SAT_AIR,HUMI_EAU,HUMI_AIR,FWW,ROAIR
502 !
503 !-----------------------------------------------------------------------
504 !
505 ! SATURATION VAPOUR PRESSURE (MAGNUS TETENS)
506  q_sat_eau = exp(2.3026d0*(7.5d0*treel/(treel+237.3d0)+0.7858d0))
507  q_sat_air = exp(2.3026d0*(7.5d0*tair/(tair+237.3d0)+0.7858d0))
508 !
509 ! AIR DENSITY : IDEAL GAZ LAW
510  roair = patm*100.d0/(287.d0*(tair+273.15d0))
511  & - 1.32d-5*hrel*q_sat_air/(tair+273.15d0)
512 !
513 ! HUMIDITY
514 ! 0.378D0 = 1.D0-0.622D0
515  humi_eau = 0.622d0*q_sat_eau/(patm-0.378d0*q_sat_eau)
516  humi_air = 0.622d0*(hrel/100.d0)*q_sat_air
517  & / (patm-(0.378d0*(hrel/100.d0)*q_sat_air))
518 ! HEAT FLUX BY EVAPORATION (SALENCON)
519 !
520  IF(n_c_atmos.EQ.2) THEN
521  fww = c1_atmos+c2_atmos*w2
522  ELSE
523 ! N_C_ATMOS.EQ.1 IS THE DEFAULT OPTION, NO OTHER OPTION THAN 1 OR 2
524  fww = c_atmos*(1.d0+w2)
525  ENDIF
526 !
527  flux_evap = roair*(2500.9d3-treel*2.365d3)*fww
528  & *(humi_eau-humi_air)
529  flux_evap = max(flux_evap,0.d0)
530 ! HEAT FLUX BY CONVECTION
531  flux_sens = cp_air*roair*fww*(treel-tair)
532 ! EVAPORATION FLOWRATE AT THE SURFACE
533  debevap = roair*fww/ro*(humi_eau-humi_air)
534  debevap = max(debevap,0.d0)
535 !
536 !-----------------------------------------------------------------------
537 !
538  RETURN
539  END SUBROUTINE evapo
540 !
541 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
542 !
543  END MODULE exchange_with_atmosphere
double precision cp_air
double precision emi_eau
double precision coef_k
double precision, parameter boltz