exchange_with_atmosphere.f

Go to the documentation of this file.

!                     *******************************
                      MODULE exchange_with_atmosphere
!                     *******************************
!
!***********************************************************************
! WAQTEL   V8P2
!***********************************************************************
!
!brief    Module containing some subroutines to deal with heat exchange
!+        with atmosphere
!
!history  N. DURAND, A. GINEAU (EDF-LNHE)
!+        MAY 2011
!+        V6P0
!+   SOLRAD, SHORTRAD, EVAPO SUBROUTINES
!+   LEAP, DAYNUM FUNCTIONS FROM SOGREAH (NOW ARTELIA)
!
!history  C.-T. PHAM (EDF-LNHE)
!+        27/06/2012
!+        V6P2
!+   Creation of module EXCHANGE_WITH_ATMOSPHERE from previous
!+   subroutines
!
!history  N. DURAND, N. LORRAIN, C.-T. PHAM (EDF-LNHE)
!+        03/04/2014
!+        V7P0
!+   Update after Nicolas Lorrain's internship
!+
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE declarations_special
      IMPLICIT NONE
!
      PRIVATE
      PUBLIC :: solrad,shortrad,evapo,leap,daynum
!
!-----------------------------------------------------------------------
!
!  brought to declarations_waqtel
!  CP: SPECIFIC HEAT OF WATER AT CONSTANT PRESSURE
!      DOUBLE PRECISION, PARAMETER :: CP        = 4.18D3
!
!-----------------------------------------------------------------------
!
      CONTAINS
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
!                   *****************
                    SUBROUTINE solrad
!                   *****************
!
     &(ray_sol,nebu,mardat,martim,at,latitude,longitude)
!
!***********************************************************************
! TELEMAC-3D V7P0                             22/06/2012
!***********************************************************************
!
!brief    EVALUATES SOLAR RADIATION INCIDENT ON THE SEA SURFACE
!+          - CALCULATES SOLAR RADIATION AS FUNCTION OF DAY NUMBER
!+            OF THE YEAR AND GEOGRAPHICAL LOCATION
!+          - INCLUDES ATMOSPHERICAL ABSORPTION AND REFLECTION, CLOUD
!+            COVERAGE, SEA SURFACE ALBEDO
!+          - TIME EXPRESSED IN GMT
!+        SOURCES:
!+          - PERRIN DE BRICHAMBAUT (1975)
!+          - BERLIAND'S METHOD (1960)
!+          - COOPER'S FORMULA (1969)
!
!history  N. DURAND, A. GINEAU (EDF-LNHE)
!+        MAY 2011
!+        V6P0
!+
!
!history  C.-T. PHAM (EDF-LNHE)
!+        JUNE 2012
!+        V6P2
!+   Adding of MARDAT, MARTIM, LATITUDE, LONGITUDE, ALBEDO AND type of
!+   sky as new arguments + INTENT + cosmetics
!+
!
!history  C.-T. PHAM (EDF-LNHE)
!+        03/04/2014
!+        V7P0
!+   Update after Nicolas Lorrain's internship
!+   Albedo is not constant in time and is no more an argument
!+   The type of sky is no more an argument, but fixed in the subroutine
!+   and may be changed by the user
!+
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| AT             |-->| CURRENT TIME
!| NEBU           |-->| NEBULOSITY (IN OCTAS)
!| LATITUDE       |-->| LATITUDE
!| LONGITUDE      |-->| LONGITUDE
!| MARDAT         |-->| DATE (YEAR, MONTH,DAY)
!| MARTIM         |-->| TIME (HOUR, MINUTE,SECOND)
!| RAY_SOL        |<--| SOLAR RADIATION INCIDENT ON THE SEA SURFACE
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE declarations_waqtel,ONLY: iskytype
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN)           :: MARDAT(3),MARTIM(3)
      DOUBLE PRECISION, INTENT(IN)  :: AT,NEBU
      DOUBLE PRECISION, INTENT(IN)  :: LATITUDE,LONGITUDE
      DOUBLE PRECISION, INTENT(OUT) :: RAY_SOL
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER IYEAR,IMONTH,IDAY,IHOUR,IMIN,ISEC
!
      DOUBLE PRECISION DTR,PI,ALB
      DOUBLE PRECISION DAY,DAYREEL,NDAYS
!  HA  : SUN'S HOUR ANGLE         [rad]
!  HR  : TIME OF THE DAY IN HOURS (GMT)
!  RDEC: SUN'S DECLINATION        [rad]
!  RLAT: LATITUDE                 [rad]
!  SING: SIN(GAMMA)
!  TE  : TIME EQUATION            [hours]
      DOUBLE PRECISION HA,HR,RDEC,RLAT,SING,TE
!  AA,BB: COEFFICIENTS DEALING WITH LUMINOSITY AND SKY COLOUR
      DOUBLE PRECISION AA,BB
!
!     INTEGER  LEAP,DAYNUM
!     EXTERNAL LEAP,DAYNUM
!
!-----------------------------------------------------------------------
!
!
      IF(iskytype.EQ.1) THEN
!  VERY PURE SKY
        aa = 1130.d0
        bb = 1.15d0
      ELSEIF(iskytype.EQ.2) THEN
!  MEAN PURE SKY
        aa = 1080.d0
        bb = 1.22d0
      ELSEIF(iskytype.EQ.3) THEN
!  INDUSTRIAL AREA
        aa = 995.d0
        bb = 1.25d0
      ENDIF
!
      iyear  = mardat(1)
      imonth = mardat(2)
      iday   = mardat(3)
      ihour  = martim(1)
      imin   = martim(2)
      isec   = martim(3)
!
!-----------------------------------------------------------------------
!
      pi  = 4.d0*atan(1.d0)
      dtr = pi/180.d0
!
!  DAY NUMBER, ORBITAL CORRECTION
      day = daynum(iyear,imonth,iday,ihour,imin,isec)
     &    + floor(at/86400.d0)
      ndays = 365.d0 + REAL(leap(iyear))
      dayreel = modulo(day, ndays)

!  ALBEDO WITH RESPECT OF THE MONTH
      IF(dayreel.GE.0.d0.AND.dayreel.LE.31.d0) THEN
        alb = 0.11d0
      ELSEIF(dayreel.GT.31.d0.AND.dayreel.LE.59.d0) THEN
        alb = 0.10d0
      ELSEIF(dayreel.GT.59.d0.AND.dayreel.LE.90.d0) THEN
        alb = 0.08d0
      ELSEIF(dayreel.GT.90.d0.AND.dayreel.LE.120.d0) THEN
        alb = 0.07d0
      ELSEIF(dayreel.GT.120.d0.AND.dayreel.LE.151.d0) THEN
        alb = 0.06d0
      ELSEIF(dayreel.GT.151.d0.AND.dayreel.LE.181.d0) THEN
        alb = 0.06d0
      ELSEIF(dayreel.GT.181.d0.AND.dayreel.LE.212.d0) THEN
        alb = 0.06d0
      ELSEIF(dayreel.GT.212.d0.AND.dayreel.LE.243.d0) THEN
        alb = 0.07d0
      ELSEIF(dayreel.GT.243.d0.AND.dayreel.LE.273.d0) THEN
        alb = 0.07d0
      ELSEIF(dayreel.GT.273.d0.AND.dayreel.LE.304.d0) THEN
        alb = 0.08d0
      ELSEIF(dayreel.GT.304.d0.AND.dayreel.LE.334.d0) THEN
        alb = 0.11d0
      ELSEIF(dayreel.GT.334.d0.AND.dayreel.LE.365.d0) THEN
        alb = 0.12d0
      ENDIF

!  DECLINATION OF SUN (COOPER'S FORMULA)
!      RDEC = (23.45D0*SIN(2.D0*PI*(DAYREEL+284.D0)/NDAYS))*DTR
      rdec = (23.45d0*cos(2.d0*pi*(172.d0-dayreel)/ndays))*dtr

!  TIME EQUATION
      te = ( 450.68d0*sin(2.d0*pi*dayreel/ndays-0.026903d0)
     &      +595.40d0*sin(4.d0*pi*dayreel/ndays+0.352835d0))/3600.d0
!  SOLAR ALTITUDE
      hr = ihour+modulo(at,86400.d0)/3600.d0
!
      rlat = latitude*dtr
      ha   = (hr-te-12.d0 + longitude/15.d0)*pi/12.d0
      sing = sin(rlat)*sin(rdec) + cos(rlat)*cos(rdec)*cos(ha)
!  SOLAR RADIATION
      IF(sing.LE.0.d0) THEN
        ray_sol = 0.d0
      ELSE
!  THE NEBULOSITY IS GIVEN IN OCTAS
        ray_sol = aa*(sing**bb)*(1.d0-0.65d0*(nebu/8.d0)**2)
     &              *(1.d0-alb)
      ENDIF
!
!-----------------------------------------------------------------------
!
      RETURN
      END SUBROUTINE solrad

!                       *********************
                        INTEGER FUNCTION leap
!                       *********************
!
     &(iyear)
!
!***********************************************************************
! TELEMAC-3D V6P2                             27/06/2012
!***********************************************************************
!
!brief    DETERMINES WHETHER IYEAR IS A LEAP YEAR
!+        DESCRIPTION - RETURNS 1 IF IYEAR IS A LEAP YEAR, 0 OTHERWISE
!+
!
!history  C. GUILBAUD (SOGREAH)
!+        JUNE 2001
!+        V6P0?
!+
!
!history  C.-T. PHAM (EDF-LNHE)
!+        27/06/2012
!+        V6P2
!+   Introduction into EXCHANGE_WITH_ATMOSPHERE module + INTENT
!+
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| IYEAR          |-->| INDEX OF YEAR
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN) :: IYEAR
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      IF( mod(iyear,4).EQ.0.AND.
     &   (mod(iyear,100).NE.0.OR.mod(iyear,400).EQ.0)) THEN
        leap = 1
      ELSE
        leap = 0
      ENDIF
!
!-----------------------------------------------------------------------
!
      RETURN
      END FUNCTION leap

!                   ********************************
                    DOUBLE PRECISION FUNCTION daynum
!                   ********************************
!
     &(iyear,imonth,iday,ihour,imin,isec)
!
!***********************************************************************
! TELEMAC-3D V6P2                             27/06/2012
!***********************************************************************
!
!brief    RETURNS DAY NUMBER OF THE YEAR (FRACTIONAL)
!+
!
!history  C. GUILBAUD (SOGREAH)
!+        JUNE 2001
!+        V6P0?
!+
!
!history  C.-T. PHAM (EDF-LNHE)
!+        27/06/2012
!+        V6P2
!+   Introduction into EXCHANGE_WITH_ATMOSPHERE module
!+   Change for type of result (double precision, not integer)
!+   + REAL conversion + addition of seconds ISEC
!+
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| IDAY           |-->| INDEX OF DAY
!| IHOUR          |-->| INDEX OF HOUR
!| IMIN           |-->| INDEX OF MINUTE
!| IMONTH         |-->| INDEX OF MONTH
!| ISEC           |-->| INDEX OF SECOND
!| IYEAR          |-->| INDEX OF YEAR
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN) :: IYEAR,IMONTH,IDAY,IHOUR,IMIN,ISEC
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
!
!     INTEGER  LEAP
!     EXTERNAL LEAP
!
      INTEGER MONTH(12)
      parameter( month=(/0,31,59,90,120,151,181,212,243,273,304,334/) )
!
!-----------------------------------------------------------------------
!
      daynum = REAL(month(imonth)+iday)
     &       + REAL(ihour)/24.d0+REAL(imin)/1440.d0+REAL(isec)/86400.d0
      IF(imonth.GT.2) daynum = daynum + REAL(leap(iyear))
!
!-----------------------------------------------------------------------
!
      RETURN
      END FUNCTION daynum

!                   *******************
                    SUBROUTINE shortrad
!                   *******************
!
     &(treel,tair,nebu,hrel,ray_atm,ray_eau)
!
!***********************************************************************
! TELEMAC-3D V7P0                             22/06/2012
!***********************************************************************
!
!brief    CALCULATES ATMOSPHERIC AND WATER RADIATIONS
!+        SOURCES:
!+          - SWINBANK'S METHOD
!+          - T.V.A. 1972
!+
!
!history  N. DURAND, A. GINEAU (EDF-LNHE)
!+        MAY 2011
!+        V6P0
!+
!
!history  C.-T. PHAM (EDF-LNHE)
!+        JUNE 2012
!+        V6P2
!+   Type of cloud taken into account + INTENT + cosmetics
!+
!
!history  C.-T. PHAM (EDF-LNHE)
!+        03/04/2014
!+        V7P0
!+   Update after Nicolas Lorrain's internship
!+   The type of cloud is no more an argument, but fixed in the subroutine
!+   and may be changed by the user
!+
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| NEBU           |-->| NEBULOSITY (IN OCTAS)
!| RAY_ATM        |<--| ATMOSPHERIC RADIATION
!| RAY_EAU        |<--| WATER RADIATION
!| TAIR           |-->| AIR TEMPERATURE
!| TREEL          |-->| REAL WATER TEMPERATURE
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE declarations_waqtel, ONLY: emi_eau,boltz,ema,iray_atm,coef_k
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      DOUBLE PRECISION, INTENT(IN)  :: TREEL,TAIR,NEBU,HREL
      DOUBLE PRECISION, INTENT(OUT) :: RAY_ATM,RAY_EAU
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
!  NUA: COEFFICIENT DEPENDING ON THE TYPE OF CLOUDS
!  REPLACED BY COEF_K (KEYWORD COEFFICIENT OF CLOUDING RATE)
!  ICLOUDTYPE FOR THE TYPE OF CLOUD IS THEN USELESS
!
!  TMP TO USE THE SAME FORMULAE AS IN GLM
      DOUBLE PRECISION EPS_STAR
      DOUBLE PRECISION QSAT_AIR
!
!-----------------------------------------------------------------------
!
!  ATMOSPHERIC RADIATION
      IF(iray_atm.EQ.1) THEN
!       CASE 1: IDSO AND JACKSON (1969)
        eps_star = (1.d0+0.275d0*(nebu/8.d0))
     &            *(1.d0-0.261d0*exp(-7.77d-8*tair**2))
      ELSEIF(iray_atm.EQ.2) THEN
!       CASE 2: SWINBANK (1963)
!       DEFAULT, LIKE VERSIONS V7P1 AND BEFORE
        eps_star = (1.d0+coef_k*(nebu/8.d0)**2)*0.937d-5
     &            *(tair+273.15d0)**2
      ELSEIF(iray_atm.EQ.3) THEN
!       CASE 3: BRUTSAERT (1975)
        qsat_air = 10.d0**(9.28603523d0-(2322.37885d0/(tair+273.15d0)))
        eps_star = (1.d0+0.275d0*(nebu/8.d0))*0.642d0
     &            *((hrel/100.d0)*qsat_air/(tair+273.15d0))**(1.d0/7.d0)
      ELSEIF(iray_atm.EQ.4) THEN
!       CASE 4: YAJIMA TONO DAM (2014)
        qsat_air = 10.d0**(9.28603523d0-(2322.37885d0/(tair+273.15d0)))
        eps_star = (1.d0-(nebu/8.d0)**2.796d0)*1.24d0
     &            *((hrel/100.d0)*qsat_air/(tair+273.15d0))**(1.d0/7.d0)
     &           + 0.955d0*(nebu/8.d0)**2.796d0
      ELSE
        WRITE(lu,*) 'FORMULA OF ATMOSPHERIC RADIATION NOT POSSIBLE'
        CALL plante(1)
      ENDIF
!
      ray_atm = ema*eps_star*boltz*(tair+273.15d0)**4
!
!  WATER RADIATION
      ray_eau = emi_eau*boltz*(treel+273.15d0)**4
!
!-----------------------------------------------------------------------
!
      RETURN
      END SUBROUTINE shortrad

!                   ****************
                    SUBROUTINE evapo
!                   ****************
!
     &(treel,tair,w2,patm,hrel,ro,flux_evap,flux_sens,debevap,c_atmos,
     & c1_atmos,c2_atmos)
!
!***********************************************************************
! WAQTEL   V8P2
!***********************************************************************
!
!brief    CALCULATES FLUX OF LATENT HEAT (W/M^2)
!+        CALCULATES SENSIBLE FLUX (W/M^2)
!+        CALCULAGES EVAPORATED WATER FLOWRATE (M/S)
!+        SOURCES:
!+          - BOLTON 1980 FOR SATURATION VAPOUR PRESSURE
!+
!
!history  N. DURAND, A. GINEAU (EDF-LNHE)
!+        MAY 2011
!+        V6P0
!+
!
!history  C.-T. PHAM (EDF-LNHE)
!+        JUNE 2012
!+        V6P2
!+   Parameter B to calibrate in arguments + cosmetics
!+
!
!history  C.-T. PHAM (EDF-LNHE)
!+        03/04/2014
!+        V7P0
!+   Update after Nicolas Lorrain's internship
!+   The salinity is not an argument anymore.
!+   Every salinity correction has been removed
!+
!
!history  M. GANT, L. ABBAS, C.-T. PHAM (EDF-LNHE)
!+        07/12/2015
!+        V7P1
!+   DEBEVAP and FLUX_EVAP have to be greater than or equal to zero.
!+   Otherwise, rain is read through files
!+
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| C_ATMOS        |-->| PARAMETER TO CALIBRATE (FWW = C_ATMOS*(1.D0+W2))
!| C1_ATMOS       |-->| PARAMETER TO CALIBRATE (FWW = C1_ATMOS+C2_ATMOS*W)
!| C2_ATMOS       |-->| PARAMETER TO CALIBRATE (FWW = C1_ATMOS+C2_ATMOS*W)
!| DEB_EVAP       |<--| EVAPORATION FLOWRATE AT THE SURFACE
!| FLUX_EVAP      |<--| ENERGY FLUX DUE TO EVAPORATED WATER
!| FLUX_SENS      |<--| HEAT FLUX BY CONVECTION
!| HREL           |-->| RELATIVE HUMIDITY
!| PATM           |-->| ATMOSPHERIC PRESSURE
!| RO             |-->| DENSITY
!| TAIR           |-->| AIR TEMPERATURE
!| TREEL          |-->| REAL WATER TEMPERATURE
!| W2             |-->| RELATIVE MAGNITUDE OF WIND AT 2 M
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE declarations_waqtel, ONLY: cp_air,n_c_atmos
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      DOUBLE PRECISION, INTENT(IN)  :: TREEL,TAIR,W2,PATM,HREL,RO
      DOUBLE PRECISION, INTENT(IN)  :: C_ATMOS,C1_ATMOS,C2_ATMOS
      DOUBLE PRECISION, INTENT(OUT) :: FLUX_EVAP,FLUX_SENS,DEBEVAP
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      DOUBLE PRECISION Q_SAT_EAU,Q_SAT_AIR,HUMI_EAU,HUMI_AIR,FWW,ROAIR
!
!-----------------------------------------------------------------------
!
!  SATURATION VAPOUR PRESSURE (MAGNUS TETENS)
      q_sat_eau  = exp(2.3026d0*(7.5d0*treel/(treel+237.3d0)+0.7858d0))
      q_sat_air  = exp(2.3026d0*(7.5d0*tair/(tair+237.3d0)+0.7858d0))
!
!  AIR DENSITY : IDEAL GAZ LAW
      roair =    patm*100.d0/(287.d0*(tair+273.15d0))
     &      - 1.32d-5*hrel*q_sat_air/(tair+273.15d0)
!
!  HUMIDITY
!  0.378D0 = 1.D0-0.622D0
      humi_eau  = 0.622d0*q_sat_eau/(patm-0.378d0*q_sat_eau)
      humi_air  =        0.622d0*(hrel/100.d0)*q_sat_air
     &          / (patm-(0.378d0*(hrel/100.d0)*q_sat_air))
!  HEAT FLUX BY EVAPORATION (SALENCON)
!
      IF(n_c_atmos.EQ.2) THEN
        fww     = c1_atmos+c2_atmos*w2
      ELSE
!     N_C_ATMOS.EQ.1 IS THE DEFAULT OPTION, NO OTHER OPTION THAN 1 OR 2
        fww     = c_atmos*(1.d0+w2)
      ENDIF
!
      flux_evap = roair*(2500.9d3-treel*2.365d3)*fww
     &                 *(humi_eau-humi_air)
      flux_evap = max(flux_evap,0.d0)
!  HEAT FLUX BY CONVECTION
      flux_sens = cp_air*roair*fww*(treel-tair)
!  EVAPORATION FLOWRATE AT THE SURFACE
      debevap   = roair*fww/ro*(humi_eau-humi_air)
      debevap   = max(debevap,0.d0)
!
!-----------------------------------------------------------------------
!
      RETURN
      END SUBROUTINE evapo
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      END MODULE exchange_with_atmosphere