thermal_khione.f

Go to the documentation of this file.

!                     *********************
                      MODULE thermal_khione
!                     *********************
!
!***********************************************************************
! KHIONE   V7P3
!***********************************************************************
!
!brief    Module containing all subroutines to deal with the physics
!+        of thermal exchanges, at a node level.
!+        To be joined up with THERMAL_WAQTEL.
!
!history  F. SOUILLE (EDF)
!+        30/09/2019
!+        V8P0
!+        Added coefficients for full heat budget
!+        Added freezing temperature
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE declarations_special
      IMPLICIT NONE
!
      PRIVATE
      PUBLIC :: thermal_fluxes,waterice_heat_coef,icover_growth,leap,
     &          daynum
!
!=======================================================================
!
!       1) THERMAL FLUXES
!
      CONTAINS
!
!=======================================================================
!
!                   *************************
                    SUBROUTINE thermal_fluxes
!                   *************************
!
     &(tair,twat,tfrz,tdew,cc,visb,wind,pluie,sumph,phcl,phri,phps,
     & phib,phie,phih,phip,anfem,dt,at,mardat,martim,lambd0,cice)
!
!***********************************************************************
! RICE-2D    V7P2                                          11/11/2016
!***********************************************************************
!
!brief
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| ANFEM   |-->| CONCENTRATION OF SURFACE ICE PARTICLES
!| AT      |-->| CURRENT TIME
!| CC      |-->| CLOUD COVER
!| DN      |-->| CURRENT TIME
!| DT      |-->| TIME STEP
!| LAMBD0  |-->| LATITUDE OF ORIGIN POINT (KEYWORD, IN DEGREES)
!| MARDAT  |-->| DATE (YEAR, MONTH,DAY)
!| MARTIM  |-->| TIME (HOUR, MINUTE,SECOND)
!| PHCL    |<->| SOLAR RAD (FLUX) REACHING SURFACE, UNDER CLEAR SKY
!| PHIB    |<->| EFFECTIVE BACK RADIATION (OPEN WATER OR ICE)
!| PHIE    |<->| EVAPORATIVE HEAT TRANSFER
!| PHIH    |<->| CONVECTIVE HEAT TRANSFER
!| PHIP    |<->| HEAT TRANSFER DUE TO PRECIPITATION
!| PHPS    |<->| NET SOLAR RADIATION (FLUX) AFTER REFLEXION
!| PHRI    |<->| SOLAR RAD (FLUX) REACHING SURFACE, UNDER CLOUDY SKY
!| PLUIE   |-->| RAIN
!| SUMPH   |<->| NET SUM OF ALL THERMAL FLUXES
!| T1      |-->| STARTING HOUR FOR SOLAR RADIATION CALCULATION (HRS)
!| T2      |-->| ENDING HOUR FOR SOLAR RADIATION CALCULATION (HRS)
!| TAIR    |-->| AIR TEMPERATURE
!| TDEW    |-->| DEWPOINT TEMPERATURE
!| TFRZ    |-->| FREEZING TEMPERATURE
!| TWAT    |-->| WATER TEMPERATURE
!| VISB    |-->| VISIBILITY
!| WIND    |-->| WIND SPEED EFFECT ON ICE
!| CICE    |-->| EXCHANGES WITH OPEN WATER (0) OR ICE COVER (1)
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE declarations_khione,      ONLY: lh_ice,cp_ice,
     &  lin_watair,cst_watair,lin_iceair,cst_iceair,
     &  coef_phib,coef_phie,coef_phih,coef_phip,sgma,
     &  cp_eau, atmoexch
      USE meteo_telemac,            ONLY: windz
!
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN)             :: MARDAT(3),MARTIM(3),CICE
!
      DOUBLE PRECISION, INTENT(IN)    :: TWAT,ANFEM,TFRZ
      DOUBLE PRECISION, INTENT(IN)    :: TAIR,TDEW,CC,VISB,WIND,PLUIE
      DOUBLE PRECISION, INTENT(IN)    :: DT,AT
      DOUBLE PRECISION, INTENT(INOUT) :: SUMPH,PHCL,PHRI
      DOUBLE PRECISION, INTENT(INOUT) :: PHPS,PHIB,PHIE,PHIH,PHIP
      DOUBLE PRECISION, INTENT(IN)    :: LAMBD0
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
!     DAY NUMBER, ORBITAL CORRECTION
      INTEGER IYEAR,IMONTH,IDAY,IHOUR,IMIN,ISEC
      INTEGER NDLT,I
!
      DOUBLE PRECISION DN, DAY,DAYREEL,NDAYS,DTHR
!
      DOUBLE PRECISION PHBA,PHBC,PHBR,PHBW
      DOUBLE PRECISION ES1,EA1,EPINA,AKN,VA,ASV
      DOUBLE PRECISION PCL,PRI,PPS,TAK,TSK,TDK
      DOUBLE PRECISION T1,T2,T10,T20
!
!-----------------------------------------------------------------------
!
      iyear  = mardat(1)
      imonth = mardat(2)
      iday   = mardat(3)
      ihour  = martim(1)
      imin   = martim(2)
      isec   = martim(3)
!
!-----------------------------------------------------------------------
!
!     PHCL: SOLAR RAD (FLUX) REACHING SURFACE, UNDER CLEAR SKY
!     PHRI: SOLAR RAD (FLUX) REACHING SURFACE, UNDER CLOUDY SKY
!     PHPS: NET SOLAR RAD (FLUX) AFTER REFLECTION
!     PHIB: BACK RADIATION HEAT TRANSFER
!     PHIE: EVAPORATIVE HEAT TRANSFER
!     PHIH: CONVECTIVE HEAT TRANSFER
!     PHIP: HEAT TRANSFER DUE TO PRECIPITATION
      sumph = 0.d0
!
!-----------------------------------------------------------------------
!
!     DAY NUMBER, ORBITAL CORRECTION
      day = daynum(iyear,imonth,iday,ihour,imin,isec)
     &    + floor(at/86400.d0)            ! DAY FROM JAN.1
      ndays = 365.d0 + REAL(leap(iyear))
      dayreel = modulo(day, ndays)        ! DAY NUMBER IN DATE
!
      dn = day       ! INT # OF DAYS FROM JAN 1 = DAYREEL
      dthr = dt / 3600.d0

      phcl = 0.d0
      phri = 0.d0
      phps = 0.d0

      phib = 0.d0
      phie = 0.d0
      phih = 0.d0
      phip = 0.d0

!     MSEC = NDAYS * 86400 + IHOUR * 3600 + IMIN * 60 + ISEC
!
!     T10 = FRACTION # OF HOURS FROM 0:00 @ CURRENT TIME (TSUM1)
!     T20 = FRACTION # OF HOURS FROM 0:00 AFTER DTHR
!
      t10 = ihour + modulo(at,86400.d0)/3600.d0
      t20 = t10 + dthr
      t1  = t10
!
      IF (dthr.GE.1.d0) THEN  ! IF TIME INTERVAL > 1.0 HR (RARE)
        ndlt = int(dthr + 0.001d0)
        DO i = 1,ndlt
          t2 = t1 + 1.d0
          IF (t2.GT.t20) THEN
            CALL solar(dn,t1,t20,cc,pcl,pri,pps,cice,lambd0)
            phcl = phcl + pcl  ! ADD FRACTION OF HOUR LEFT
            phri = phri + pri
            phps = phps + pps
            EXIT
          ELSE
            CALL solar(dn,t1,t2,cc,pcl,pri,pps,cice,lambd0)
            phcl = phcl + pcl  ! ADD HOUR INCREMENTS
            phri = phri + pri
            phps = phps + pps
          ENDIF
          t1 = t2
        ENDDO
      ELSE
        CALL solar(dn,t10,t20,cc,pcl,pri,pps,cice,lambd0)
        phcl = phcl + pcl  ! IF TIME INTERVAL A FRACTION OF AN HOUR
        phri = phri + pri
        phps = phps + pps
      ENDIF
!
      phps = phps / dthr  ! NET SOLAR RAD AT SURFACE, W/M^2
!
!-----------------------------------------------------------------------
!
!     LINEAR HEAT TRANSFER FOR AIR-WATER INTERFACE ONLY
!
      IF( atmoexch.EQ.0 ) THEN
        IF(cice.EQ.1) THEN  ! ICE
          phih = cst_iceair + ( tair-tfrz )*lin_iceair
        ELSE              ! OPEN WATER
          phih = cst_watair + ( tair-twat )*lin_watair
        ENDIF
        sumph = phps + phih
!
!-----------------------------------------------------------------------
!
!     FULL BUDGET FOR AIR-WATER INTERFACE
!
      ELSEIF( atmoexch.EQ.1 ) THEN
!
        tak = tair + 273.16d0  ! TAK = AIR TEMPERATURE
        tdk = tdew + 273.16d0  ! TDK = DEW POINT
!
        IF (cice.EQ.0) THEN
          tsk = twat + 273.16d0  ! TSK = WATER TEMP
        ELSE
          tsk = 273.16d0
        ENDIF
!
!       ~~~~~~~~~~~~~~
!       BACK RADIATION
!       ~~~~~~~~~~~~~~
        IF (cice.EQ.0) THEN
!
!         FOR WATER SURFACE, ES1 FOUND USING WATER TEMP (TSK), (4.32)
          es1 = 7.95357242d10*exp((-18.1972839d0*373.16d0/tsk)+
     &          5.02808d0*log(373.16d0/tsk)-20242.1852d0*
     &          exp(-26.1205253d0*tsk/373.16d0)+
     &          58.0691913d0*exp(-8.039282d0*373.16d0/tsk))
          ea1 = 7.95357242d10*exp((-18.1972839d0*373.16d0/tdk)+
     &          5.02808d0*log(373.16d0/tdk)-20242.1852d0*
     &          exp(-26.1205253d0*tdk/373.16d0)+
     &          58.0691913d0*exp(-8.039282d0*373.16d0/tdk))
        ELSE
!
!         FOR ICE SURFACE, ES1 FOUND USING AIR TEMP (TAK), (4.33)
!         USUALLY CONSIDERED WHEN ANFEM(I) > 0.5,
!         MORE ICE THAN WATER AT SURFACE
          es1 = 5.75185606d10*exp((-20.947031d0*273.16d0/tak)-
     &          3.56654d0*log(273.16d0/tak)-2.01889049d0/273.16d0*tak)
          ea1 = 5.75185606d10*exp((-20.947031d0*273.16d0/tdk)-
     &        3.56654d0*log(273.16d0/tdk)-2.01889049d0/273.16d0*tdk)
        ENDIF
!
!       EMISSIVITY OF ATMOSPHERE = PHBA (INCL. EFFECTS OF CLOUDS)
        IF (tair.LT.0.d0) THEN
          ! (4.31) SATTERLUND (1979)
          epina = 1.08d0 * (1.d0 - exp(-ea1 ** (tak / 2016.d0)))
        ELSE
          ! TK 03-2010 IDSO-JACKSON(1969)
          epina = 1.d0-0.261d0*exp(-0.000777d0*(273.16d0-tak)**2)
        ENDIF
        ! (4.30)
        phbc = epina * sgma * tak ** 4
        ! = ATMOSPHERIC RADIATION UNDER CLOUDY SKY (TELEMAC 2D)
        phba = phbc * (1.d0 + 0.0017d0 * cc **2)
!
!       EMISSIVITY OF RIVER SURFACE = PHBW
        ! (4.29), TSK = WATER TEMP (ALWAYS)
        phbw = 0.97d0 * sgma * tsk ** 4

!       REFLECTED LONG WAVE RADIATION = PHBR
        ! (4.36)
        phbr = 0.03d0 * phba
!
!       EFFECTIVE BACK RADIATION
        phib = phba - phbr - phbw
        phib = coef_phib*phib
!
!       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!       EVAPORATIVE HEAT FLUX AND CONVECTIVE HEAT FLUX
!       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        IF (cice.EQ.0) THEN
          akn = 8.d0 + 0.35d0 * (twat - tair)
          va = (2.d0/windz) ** 0.15d0 * wind
!         WINDZ = HEIGHT WIND VELOCITY IS MEASURED (*.PAR)
!
!         EVAPORATION
          phie = -(1.56d0*akn + 6.08d0 * va) * (es1-ea1)*4.1855d0/8.64d0
          phie = coef_phie*phie
        ELSE
          akn = 8.d0 + 0.35d0 * (0.d0 - tair)
          va = (2.d0/windz) ** 0.15d0 * wind
        ENDIF
!
!       CONDUCTIVE HEAT TRANSFER
        phih = (akn + 3.9d0*va) * (tak-tsk)*4.1855d0/8.64d0
        phih = coef_phih*phih
!
!       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!       HEAT TRANSFER DUE TO PRECIPITATION
!       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        IF (visb.GT.0.d0.AND.visb.LT.1.d0) THEN
!         MIN VALUE = 1.0 KM, PHIP = 380.28 W/M^2
          asv = 78.5d0 / 86400.d0
        ELSE
!         AT VISB = 10 KM, PHIP = 1.604 W/M^2
          asv = 78.5d0 / 86400.d0 * visb ** (-2.375d0)
        ENDIF
!
!       NO SNOW
        IF(visb.LT.1.d-5) THEN ! NO SNOW
          asv = 0.d0
        ENDIF
!
!       SNOW FALL
        IF(cice.EQ.0) THEN
          phip = asv * (cp_ice * (tair - twat) - lh_ice)
        ELSE
          phip = 0.d0
        ENDIF
!
!       RAIN FALL
        IF(pluie.GT.0.d0) THEN
          asv = pluie/3600.d0
          IF(cice.EQ.0) THEN   ! WATER
            IF(tair.LT.0.d0) THEN
              phip = - asv*cp_eau*(0.d0 - twat) ! HEAT LOSS
            ELSE
              phip = - asv*cp_eau*(tair - twat) ! HEAT GAIN
            ENDIF
          ELSE              ! ICE
            IF(tair.GT.0) THEN
              phip = - asv*cp_eau*(tair - 0.d0) ! HEAT GAIN
            ELSE
              phip = 0.d0
            ENDIF
          ENDIF
        ENDIF
        phip = coef_phip*phip
!
!       ~~~~~~~~~~~~~~~~~~~~
!       SUMMATION AND OUTPUT
!       ~~~~~~~~~~~~~~~~~~~~
        sumph = phps + phib + phie + phih + phip
!
      ENDIF
!
!-----------------------------------------------------------------------
!
      RETURN
      END SUBROUTINE thermal_fluxes
!
!=======================================================================
!
!                   ************************
                    SUBROUTINE icover_growth
!                   ************************
!
     &(twat,tmelt,sumph,anfem,thifems,hwi,dt)
!
!***********************************************************************
! RICE-2D    V7P3                                          11/11/2016
!***********************************************************************
!
!brief
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| DT      |-->| TIME STEP
!| HWI     |-->| WATERICE HEAT COEFFICIENT
!| SUMPH   |-->| NET HEAT EXCHANGE FLUX WITH THE ATMOSPHERE
!| ANFEM   |<->| CONCENTRATION OF SURFACE ICE PARTICLES
!| THIFEMS |<->| SOLID ICE THICKNESS
!| TMELT   |-->| FREEZING POINT OF WATER
!| TWAT    |-->| WATER TEMPERATURE
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE declarations_khione, ONLY: rho_ice,lh_ice,tc_bi,lin_iceair
!
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      DOUBLE PRECISION, INTENT(IN)    :: HWI
      DOUBLE PRECISION, INTENT(IN)    :: TWAT,TMELT,SUMPH
      DOUBLE PRECISION, INTENT(IN)    :: DT
      DOUBLE PRECISION, INTENT(INOUT) :: ANFEM,THIFEMS
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
!
      DOUBLE PRECISION AUX,DH
!
!-----------------------------------------------------------------------
!
      aux = dt/(rho_ice*lh_ice)
      dh = - aux*sumph/(thifems*lin_iceair/tc_bi + 1.d0)
     &     - aux*hwi*(twat-tmelt)
!
      IF( dh.GT.(-thifems) ) THEN
        thifems = thifems + dh
      ELSE
        thifems = 0.d0
      ENDIF
!
!-----------------------------------------------------------------------
!
!     NOT ALLOWING NEGATIVE ICE COVER THICKNESS
!
      IF( thifems.LE.0.d0 ) thifems = 0.d0
      IF( thifems.LE.0.d0 ) anfem = 0.d0
!
!-----------------------------------------------------------------------
!
      RETURN
      END SUBROUTINE icover_growth
!
!=======================================================================
!
!                     ****************
                      SUBROUTINE solar
!                     ****************
     &(dn,t1,t2,cc,phcl,phri,phps,cice,lambd0)
!
!***********************************************************************
! KHIONE   V7P3
!***********************************************************************
!
!brief    Same as SOLRAD within EXCHANGE_WITH_ATMOSPHERE (WAQTEL)
!+        TODO: Combined with SOLAR
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| CC      |-->| CLOUD COVER, IN THENTHS, 0-10 (/!\ NOT IN OCTAS)
!| CICE    |-->| ICE CONDITION, OPEN WATER : 0 ; ICE : 1
!| DN      |-->| CURRENT TIME (THE NUMBER OF DAYS FROM JAN. 1)
!| T1      |-->| STARTING HOUR FOR SOLAR RADIATION CALCULATION (HRS)
!| T2      |-->| ENDING HOUR FOR SOLAR RADIATION CALCULATION (HRS)
!| PHCL    |<->| SOLAR RAD (FLUX) REACHING SURFACE, UNDER CLEAR SKY
!| PHPS    |<->| NET SOLAR RADIATION (FLUX) AFTER REFLEXION
!| PHRI    |<->| SOLAR RAD (FLUX) REACHING SURFACE, UNDER CLOUDY SKY
!| LAMBD0  |-->| LATITUDE OF ORIGIN POINT (NORTH +, SOUTH -, DEGREES)
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE declarations_khione, ONLY: alsm,allm,etadir,sio,albe
      USE meteo_telemac,       ONLY: modelz,alphsd,alphrd
!
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER         , INTENT(IN)    :: CICE
      DOUBLE PRECISION, INTENT(IN)    :: CC
      DOUBLE PRECISION, INTENT(IN)    :: DN,T1,T2
      DOUBLE PRECISION, INTENT(INOUT) :: PHCL,PHRI,PHPS
      DOUBLE PRECISION, INTENT(IN)    :: LAMBD0
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      DOUBLE PRECISION :: PI
      DOUBLE PRECISION :: AA, ALPHA,ALPHR,ALPHS,ALHA
      DOUBLE PRECISION R,BB,AMO,AM,CR,DELTA,EO,ET,HRIS,HSS,
     &  papo,phiso,phrr,phs,rhs1,rhs2,sia,alr,hs1,hs2,hset
!
!-----------------------------------------------------------------------
!
! DELTA = SOLAR DECLINATION OF THE SUN;
! HSS = LOCAL HOUR ANGLE OF THE SUNSET, RADIANS;
! HSR = LOCAL HOUR ANGLE OF THE SUNRISE, RADIANS;
! ALMST = LOCAL MEAN SOLAR TIME IN HRS;
! ALST = STANDARD TIME OF THE TIME ZONE, IN HRS, COUNTED FROM MIDNIGHT, 0-24.0;
! ALSM = LONGITUDE OF THE STANDARD MERIDEAN, IN DEGREE; (*.PAR FILE)
! ALLM = LONGITUDE OF THE LOCAL MERIDEAN, IN DEGREES; (*.PAR FILE)
! ETA = -1, FOR WEST; + FOR EAST; (*.PAR FILE)
! MODELZ = ELEVATION ABOVE SEA LEVEL, M. (*.PAR FILE)
! ALPHSD = SUN EXIT ANGLE, 180 FOR HORIZONTAL; (*.PAR FILE)
! ALPHRD = SUN EMISION ANGLE, 0 DEGREE FOR HORIZONTAL; (*.PAR FILE)
! NM = NUMBER OF MONTH
! ND = DAY NUMBER IN THE DATE
!
!-----------------------------------------------------------------------
!
      phcl = 0.d0
      phri = 0.d0
      phps = 0.d0
!
      pi  = 4.d0*atan(1.d0)

!     LOCAL GEOGRAPHIC LATITUDE, CONVERT TO RADIANS
      phs = lambd0 * pi / 180.d0  ! LATITUDE
      alphs = alphsd * pi / 180.d0  ! EXIT ANGLE, RADIANS
      alphr = alphrd * pi / 180.d0  ! EMISSION ANGLE, RADIANS

!     COOPER 1969 SOLAR DECLINATION, IN RADIANS  (4.6)
      delta = 23.45d0*pi/180.d0*sin(360.d0*(284.d0+dn)/365.d0*pi/180.d0)

!     DIFFERENCE BETWEEN TRUE SOLAR TIME AND MEAN SOLAR TIME
      ! (4.3)
      r = 2.d0 * pi * (dn - 1.d0) / 365.d0
!     DUFFIE AND BECKMAN 1959, ECCENTRICITY CORRECTION FACTION OF THE EARTH ORBIT
      ! (4.5)
      eo = 1.d0 + 0.033d0 * cos(2.d0 * pi / 365.d0 * dn)
!     EQUATION OF TIME, IN HRS  (4.2)
      et = 3.8197d0 * (0.000075d0 + 0.001868d0 * cos(r)
     &  -0.032077d0*sin(r)-0.014615d0*cos(2.d0*r)-0.04089d0*sin(2.d0*r))

!     HOUR ANGLE AT SUNRISE, RADIANS  (4.9)
      hss = acos( -tan(phs) * tan(delta) )
      IF (hss.LT.0.d0) THEN
        hss=-hss
      ENDIF
!     SUN SET, HRS  (4.11)
      hset = 12.d0 + hss * 12.d0 / pi - (pi - alphs) / pi * 12.d0
!     SUN RISE, HRS  (4.10)
      hris = 12.d0 - hss * 12.d0 / pi + alphr / pi * 12.d0
!     CALCULATE HOUR ANGLE, TIME CORRECTION  (4.4)
      hs1 = t1 - etadir / 15.d0 * (alsm - allm) + et  ! HOURS
      rhs1 = (12.d0 - hs1) * pi / 12.d0   ! RADIANS
      hs2 = t2 - etadir / 15.d0 * (alsm - allm) + et  ! HOURS
      rhs2 = (12.d0 - hs2) * pi / 12.d0   ! RADIANS
!     RHS1/RHS2 IS # OF RADIANS OFF OF NOON, + BEFORE, - AFTER

!     NET SOLAR RADIATION, PHPS
      phcl = 0.d0  ! CLEAR SKY SOLAR RADIATION
      phri = 0.d0  ! INCL. CLOUD EFFECTS, IF ANY
      phps = 0.d0  ! NET SOLAR, AFTER REFLECTION AT EARTH SURFACE

!     NO SUN, BEFORE SUNRISE OR AFTER SUNSET
      IF (hs1.GT.hset.AND.hs2.GT.hset) THEN
        RETURN
      ENDIF
      IF (hs1.LT.hris.AND.hs2.LT.hris) THEN
        RETURN
      ENDIF

      IF (hs1.LT.hris .AND. hs2.GT.hris) THEN
        rhs1=(12.d0-hris)*pi/12.d0
        hs1=hris
      ENDIF
      IF (hs2.GT.hset.AND.hs1.LT.hset) THEN
        rhs2=(12.0-hset)*pi/12.d0
        hs2=hset
      ENDIF

!     AVERAGE SOLAR TIME ANGLE
      alha = (rhs1 + rhs2) / 2.d0  ! RADIANS
!     TOTAL SOLAR RADIATION T1 TO T2, PER UNIT AREA, (4.13)-RADIANS, (4.14)-HOURS
      phiso = 12.d0/pi*sio*eo*((rhs1-rhs2)*sin(delta)*sin(phs)+
     &        (sin(rhs1)-sin(rhs2))*cos(delta)*cos(phs))

      sia = sin(delta) * sin(phs) + cos(delta) * cos(phs) * cos(alha)
!     LINE ABOVE -- (4.7)
      alpha = asin( sia ) * 180.d0 / pi  ! CONVERT TO DEGREES
      ! (4.17)
      amo = 1.d0 / (sia + 0.15d0 * (alpha + 3.885d0) ** (-1.253d0))
      ! (4.18)
      papo = exp(-0.0001184d0 * modelz)
      ! (4.16)
      am = amo * papo
      ! (4.15)
      am = 0.99d0 - 0.17d0 * am
      ! ENERGY FLUX, REACHING GROUND UNDER CLEAR SKY
      phcl = am * phiso

      IF (phcl.LT.0.d0) phcl=0.d0
      ! (4.19)
      phri = phcl * ( 1.d0 - 0.0065d0 * cc **2 )
!     SOLAR RADIATION REACHING THE EARTH UNDER CLOUDY SKIES (ABOVE)

      IF (cice.EQ.0) THEN  ! IF ALBE = 0, "OPEN WATER" CONDITIONS
        IF (phcl.EQ.0.d0) THEN
          phrr=0.d0
        ELSE  ! ESTIMATE REFLECTIVITY OF OPEN WATER SURFACE
          ! (4.25) - CLOUDINESS RATIO
          cr = max((1.d0 - phri / phcl),0.d0)
          aa =  2.2d0 +cr**0.7d0/4.d0 -  (cr**0.7d0 - 0.4d0)**2 / 0.16d0
          bb = -1.02d0+cr**0.7d0/16.d0 + (cr**0.7d0 - 0.4d0)**2 / 0.64d0
          ! ESTIMATE OPEN WATER ALBEDO (4.22-4.24)
          alr = aa * alpha ** bb
          ! AMOUNT OF SOLAR RADIATION REFLECTED
          phrr = alr * phri
        ENDIF

!     USE THE ALBEDO SPECIFIED TO INCLUDE REFLECTIVE EFFECTS OF ICE PRESENT
      ELSE
        ! (4.20)
        phrr = albe * phri
        IF (phcl.EQ.0.) phrr=0.d0
      ENDIF
      ! NET SOLAR RADIATION BETWEEN T1-T2, (4.20)
      phps = phri - phrr

      RETURN
      END SUBROUTINE solar
!
!=======================================================================
!
!             ********************************************
              DOUBLE PRECISION FUNCTION waterice_heat_coef
!             ********************************************
     &(h,u,twat,tfrz)
!
!***********************************************************************
! KHIONE   V7P3
!***********************************************************************
!
!brief    Computes the heat exchange coefficient with the ice cover
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| H       |-->| PROPAGATION DEPTH
!| U       |-->| VELOCITY MAGNITUDE
!| TWAT    |-->| WATER TEMPERATURE
!| TFRZ    |-->| FREEZING POINT FOR WATER
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE declarations_khione, ONLY: tc_wt,xnu,cwi1,ciw1,ata
!
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      DOUBLE PRECISION, INTENT(IN) :: H,U,TWAT,TFRZ
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      DOUBLE PRECISION RE, DH
!
!-----------------------------------------------------------------------
!
!     REYNOLDS NUMBER
      dh = 4.d0*h
      re = u * dh / xnu
!
!     DRY BANKS
      IF(dh.LE.0.001d0.AND.dh.GT.0.d0) THEN
        waterice_heat_coef = 0.d0 !1394.D0
        RETURN
      ENDIF
!
!     LAMINAR FLOW
      IF(re.LT.2200.d0) THEN
        waterice_heat_coef = ata*tc_wt/dh
!
!     TURBULENT FLOW
      ELSEIF( twat.GT.tfrz ) THEN ! WATER TEMP.>TFRZ
        waterice_heat_coef = cwi1*u**0.8d0/dh**0.2d0 ! CWI1 = 1448
      ELSE ! WATER TEMP.<TFRZ
        waterice_heat_coef = ciw1*u**0.8d0/dh**0.2d0 ! CIW1 = 1118
      ENDIF
!
      RETURN
!
!-----------------------------------------------------------------------
!
      END FUNCTION waterice_heat_coef
!
!                       *********************
                        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

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      END MODULE thermal_khione