 MODULE sflxrad
!!======================================================================
!!                   ***  MODULE  sflxrad  ***
!! Surface fluxes : estimate of atmospherical radiative fluxes
!!======================================================================
!! * Modules used
   USE parameters,ONLY  : ljmin,ljmax,ljminp1,rsh,rlg,valmanq,imin,imax
   USE comvars2d, ONLY  : h0,ierrorlog
   USE cometemp,  ONLY  : l_rapflux

   IMPLICIT NONE
   PRIVATE

!! * Accessibility
   PUBLIC sflx_ir,sflx_solar     ! function called by sflx_surf

!! * Interface
   INTERFACE sflx_solar

     MODULE PROCEDURE sflx_solar_luyten

   END INTERFACE
   INTERFACE sflx_ir

     MODULE PROCEDURE sflx_ir_swimbank

   END INTERFACE

!! * Shared module variables

!! * Private variables
   REAL(KIND=rsh), PARAMETER    :: kelvin=273.15_rsh, & ! shift to use kelvin unit
                                   stefan = 5.67e-8,  & ! Stefan constant (W/m2/KV4)
                                   em = 0.985_rsh       ! emissivity

 CONTAINS

!!======================================================================

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE sflx_solar  ***
!&E
!&E ** Purpose : short-wave heat flux estimation
!&E
!&E ** Description : the rate of absorption of solar radiation is usually
!&E              calculated as a product of (a) S_0, the net downward
!&E              flux of solar radiation just above the surface in
!&E              cloudless conditions (this is usually 0.7 - 0.75 of the
!&E              flux incident at the top of the atmosphere), (b) (1-alpha),
!&E              where alpha is the surface albedo, and (c) a correction
!&E              factor for cloud effect
!&E
!&E ** Options :
!&E      - Gill, 1982
!&E      - Luyten et al, 1992
!&E        see METEOCEAN report (2009)
!&E
!&E---------------------------------------------------------------------

!!======================================================================

  SUBROUTINE sflx_solar_gill(i,j,zfcc,zsinfiu,zcosfiu,zrad)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE sflx_solar_gill  ***
!&E
!&E ** Purpose : short-wave heat flux estimation
!&E
!&E ** Description : estimate of net incident solar flux
!&E              from POL bulf formula of NOMADS 2 program
!&E
!&E              Q_rad = S_0 sin(z) ( 1 - xhi + Tr)
!&E                      ( 1 - c_1 nu + c_2 beta) ( 1 - alpha )
!&E              - S_0      ! solar constant
!&E              - xhi      ! absorption coefficient by atmopherical layer
!&E              - Tr       ! diffusion coefficient in atmopherical layer
!&E              - c_1, c_2 ! constants
!&E              - nu       ! fractional cloud cover (0 < nu < 1)
!&E              - beta     ! sun altitude at midday
!&E              - alpha    ! albedo
!&E
!&E              sin(z) = sin(latitude) sin(declination) -
!&E                       cos(latitude) cos(declination) cos(2 pi t / T)
!&E              - T=24 ! number of hours in a day
!&E              - t    ! time in hours  tm=0 at midnight January, 1st
!&E
!&E              declination = 25°05min pi / 180 cos(2 pi t / T_1 - 2.950)
!&E              - T_1=365.24 days * 24 hours ! number of hours in a year
!&E
!&E
!&E ** Local variables
!&E
!&E              i,j    ! coordinate of the grid point
!&E              tm (=t)! time in hours  tm=0 at midnight January 1st
!&E              decln  !  maximum declination angle of the earth
!&E              w0     !  angular frequency of one year period  (/hours)
!&E              w1     !  angular frequency of one day period  (/hours)
!&E              sc (=S_0)     !  solar constant  1368.0  (W/m2)
!&E              snh (=sin(z))   !  sine of altitude angle of the sun
!&E              h1     !  altitude angle of the sun  (radians)
!&E              b1 (=1-xhi+Tr)    !  transmission coefficient  (b1 = 0.76)
!&E              alb (=alb) ! albedo (alb = 0.06)
!&E              c1 (=c_1) ! constant cc = 0.4
!&E              c2 (=c_2) ! constant    = 0.38
!&E              zfcc (=nu)  ! fractional cloud cover ( 0 < zfcc < 1 )
!&E              zfcc*zfcc (=beta)
!&E              jjulien number of the day in the year
!&E
!&E
!&E ** Called by : sflx_surf
!&E
!&E ** External calls :
!&E
!&E ** Reference : Gill (p34) + private modifications
!&E                METEOCEAN report (2009)
!&E
!&E ** History :
!&E       !  Original code
!&E       !  2006 (V. Garnier)  references
!&E       !  2011/01 (V. Garnier)  corrections :
!&E                  sin(z) = sin(latitude) sin(declination) +
!&E                           cos(latitude) cos(declination) cos(2 pi t / T)
!&E                  with 0<t<24h (TU + lag relative to the longitude)
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   USE parameters,ONLY : pi
   USE comvars2d, ONLY : jjulien,iheure,isec,iminu,lon2d

!! * Declaration function

!! * Arguments
   INTEGER,        INTENT( IN )  :: i,j
   REAL(KIND=rsh), INTENT( IN )  :: zsinfiu, & !
                                    zcosfiu, & !
                                    zfcc       ! cloud cover fraction    [0<zfcc<1]
   REAL(KIND=rsh), INTENT( OUT ) :: zrad       ! solar flux

!! * Local declarations
   REAL(KIND=rsh), PARAMETER                    :: alb=0.06_rsh,    & ! albedo
                                                   sc=1368.0_rsh,   & ! solar constant
                                                   b1=0.76_rsh,     & ! transmission coefficient
                                                   c1=0.4_rsh,      & ! cst to estimate the solar extinction
                                                   c2=0.38_rsh
   REAL(KIND=rsh)                               :: decln,snh,w0,w1,xcosd,xsind,cosw1tm
   REAL(KIND=rlg)                               :: tm,th

!!----------------------------------------------------------------------
!! * Executable part

! time in hours
   tm=(jjulien-1)*24.0_rlg+iheure+iminu/60.0_rlg+isec/3600.0_rlg

   decln = 23.5_rsh*pi/180.0_rsh
   w0 = 2.0_rsh*pi/(365.24_rsh*24.0_rsh)
   w1 = 2.0_rsh*pi/24.0_rsh
   xcosd = cos(decln*cos(w0*tm - 2.950_rsh))
   xsind = sin(decln*cos(w0*tm - 2.950_rsh))
   th=REAL(iheure,rlg)+REAL(iminu,rlg)/60.0_rlg+REAL(isec,rlg)/3600.0_rlg+lon2d(i,j)/15.0_rlg
   cosw1tm = cos(w1*(th-12.0_rlg))


! Estimate the sine of the angle of the sun above the horizon, snh
! d is the declination angle - June 21st is the 171st day after tm=0
! ------------------------------------------------------------------

         snh = MAX(MIN(1.0_rsh,zcosfiu*xcosd*cosw1tm + zsinfiu*xsind),0.0_rsh)

! Take into account the effect of atmospheric absorption
! h1 = asin(snh)
! Defant 1961 suggests i(z)=i0*exp(-t*a/sin(h1))
! where  t             turbidity factor (unspecified)
!        a             0.128 to 0.054 times -ln(sin(h1))
!        b1=exp(-t*a)  transmission coefficient between 0.6 and 0.7
! Comparison with Bunker & Goldsmith averages gives b1=0.76 or a=0.274
! ----------------------------------------------------------------------

         zrad=sc*snh*b1*(1.0_rsh-alb)*(1.0_rsh-zfcc*c1-c2*zfcc**2)

  END SUBROUTINE sflx_solar_gill

!!======================================================================

  SUBROUTINE sflx_solar_luyten(i,j,zfcc,zsinfiu,zcosfiu,zrad)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE sflx_solar_luyten  ***
!&E
!&E ** Purpose : short-wave heat flux estimation
!&E
!&E ** Description : estimate of net incident solar flux
!&E              from bulf formula of MAST 3 program (Luyten et al, 1992)
!&E              (DEFAULT OPTION)
!&E
!&E              Q_rad = S_0 P sin(z) ( 1 - xhi + Tr) / 2
!&E                      ( 1 - c_1 nu + c_2 beta) ( 1 - alpha )
!&E              - S_0      ! solar constant
!&E              - P        ! orbital correction
!&E              - xhi      ! absorption coefficient by atmopherical layer
!&E              - Tr       ! diffusion coefficient in atmopherical layer
!&E              - c_1, c_2 ! constants
!&E              - nu       ! fractional cloud cover (0 < nu < 1)
!&E              - beta     ! sun altitude at midday
!&E              - alpha    ! albedo
!&E
!&E              sin(z) = sin(latitude) sin(declination) +
!&E                       cos(latitude) cos(declination) cos(2 pi t / T)
!&E              - T=24 ! number of hours in a day
!&E              - t    ! time in hours  tm=0 at midnight January, 1st
!&E
!&E              declination = d_0 + SUM_n=0_n=3 ( a_n cos(nJprime) + b_n sin(nJprime) )
!&E              - d_0=0.33281 ! mean inclination
!&E              - a_n (-22.984,-0.34990,-0.1398)
!&E              - b_n (3.7872,0.03205,0.07187)
!&E
!&E
!&E ** Local variables
!&E
!&E      aal (=xhi) : water vapor plus ozone absorption coefficient
!&E      alb (=alpha): albedo of sea surface
!&E      zfcc (=nu) : fractional cloud cover ( 0 < zfcc < 1 )
!&E      solc (=S_0) ! solar constant (1367 W/m2)
!&E
!&E      be (=beta) : sun s altitude at noon (deg)
!&E      dec : sun s declination (deg)
!&E      dgr : conversion factor (degrees > radians)
!&E      gam (=z) : sun s altitude (deg.)
!&E      ha : sun s hours angle (hours)
!&E      iheure : hour of the day (GMT)
!&E      iminu: minute (GMT)
!&E      jjulien (=J) : number of the day in the year
!&E      om (=m) optical air mass
!&E      orb (=P) : correction for the variation of the earth-sun distance
!&E      qd: direct component of solar radiation under cloudless
!&E               conditions (J/(m2*s))
!&E      qdif : diffuse sky radiation under cloudless conditions (J/(m2*s))
!&E      q0 (=S_0 P sin(z)): solar radiation uncorrected from atmospherical attenuation
!&E               (J/(m2*s))
!&E      rday : day number (radians)
!&E      rdec : sun s declination (radians)
!&E      r0 (=D_r) Raleigh s optical thickness
!&E      te  time equation
!&E      tf (=T_L) : Linke s total trouble factor (after Dogniaux)
!&E      var1 (=sin(z))
!&E
!&E ** Called by : sflx_surf
!&E
!&E ** External calls :
!&E
!&E ** Reference : (Luyten et al, 1992)
!&E      zrad=(qd+qdif)*(1.0-0.62*zfcc+0.0019*be)*(1.0-alb)
!&E                METEOCEAN report (2009)
!&E
!&E ** History :
!&E       !  02-05  ()  Original code
!&E       ! 2011/01 (V. Garnier) seconds introduced into the time
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   USE parameters,ONLY : pi
   USE comvars2d, ONLY : jjulien,iheure,iminu,isec,lat2d,lon2d

!! * Declaration function

!! * Arguments
   INTEGER,        INTENT( IN )  :: i,j
   REAL(KIND=rsh), INTENT( IN )  :: zsinfiu, & !
                                    zcosfiu, & !
                                    zfcc       ! cloud cover fraction    [0<zfcc<1]
   REAL(KIND=rsh), INTENT( OUT ) :: zrad       ! solar flux

!! * Local declarations
   INTEGER                     :: k
   REAL(KIND=rsh)              :: var1,gam
   REAL(KIND=rsh)              :: aal,alb,solc
   REAL(KIND=rsh), DIMENSION(3) :: a,b,c,d
   REAL(KIND=rlg)    :: orb,rday,dgr,dec,rdec,te,ha,q0,om,ro,w1,w2,tf,qd,qdif,be

!!----------------------------------------------------------------------
!! * Executable part

   aal=0.09_rsh
   alb=0.07_rsh
   solc=1367.0_rsh

   a(1)=-22.984_rsh
   a(2)=-0.34990_rsh
   a(3)=-0.1398_rsh
   b(1)=3.7872_rsh
   b(2)=0.03205_rsh
   b(3)=0.07187_rsh
   c(1)=0.0072_rsh
   c(2)=-0.0528_rsh
   c(3)=-0.0012_rsh
   d(1)=-0.1299_rsh
   d(2)=-0.1565_rsh
   d(3)=-0.0041_rsh
   dgr=pi/180.0_rlg

! orbital correction (P factor)
! -----------------------------

! Jprime = rday = K_2 J, K_2=0.9856
   rday=0.9856_rlg*REAL(jjulien,rlg)*dgr
! P = orb = 1 + K_1 cos(Jprime - 2.8), K_1=0.03344
   orb=1.0_rlg+0.03344_rlg*cos(rday-2.8_rlg*dgr)   ! orbital correction

! declination = d_0 + SUM_n=1_n=3 ( a_n cos(nJprime) + b_n sin(nJprime) )
! ---------------------------------------------------------------

   dec=0.32281_rlg
   DO k=1,3
     dec=dec+a(k)*cos(k*rday)+b(k)*sin(k*rday) ! sun s declination
   END DO
   rdec=dec*dgr

! t = ha = t_h -12 + TE + theta
! t_h, real hour in the day
! TE, time equation
! theta, longitude in hours divided by 15
! ---------------------------------------

! TE = SUM_n=1_n=3 ( c_n cos(nJ') + d_n sin(nJ') )
   te=0.0_rlg
   DO k=1,3
     te=te+c(k)*cos(k*rday)+d(k)*sin(k*rday)   ! time equation
   END DO

! sin(z)  sun s hour angle and altitude
!--------------------------------

! time in hours
         ha=(REAL(iheure,rlg)+REAL(iminu,rlg)/60.0_rlg+REAL(isec,rlg)/3600.0_rlg-12.0_rlg+te+lon2d(i,j)/15.0_rlg)*pi/12.0_rlg
! sin(z)
         var1=zsinfiu*sin(rdec)+zcosfiu*cos(rdec)*cos(ha)
         gam=asin(var1)/dgr

! S_0 P sin(z) radiation at the top of the atmosphere
!---------------------------------------

         q0=solc*orb*var1

! atmospherical attenuation
!---------------------------

         IF (var1 < 0.0_rsh) THEN

           q0=0.0_rlg
           om=0.0_rlg
           tf=0.0_rlg
           ro=0.0_rlg
           w1=0.0_rlg
           w2=0.0_rlg

         ELSE

           om=1.0_rlg/REAL((var1+0.15_rsh*(gam+3.885_rsh)**(-1.253_rsh)),rlg)
           ro=1.0_rlg/(0.9_rlg*om+9.4_rlg)
           w1=4.0_rlg
           w2=0.1_rlg
           tf=(gam+85.0_rlg)/(39.5_rlg*exp(-w1)+47.4_rlg)+0.1_rlg+(16.0_rlg+0.22_rlg*w1)*w2

         END IF

! solar irradiance
!-----------------

! T_r = exp(-om*ro*tf)
         qd=q0*exp(-om*ro*tf)
! S_0 P sin(z) ( 1 - xhi - Tr) / 2
         qdif=((1.0_rlg-aal)*q0-qd)/2.0_rlg

! beta
         be=90.0_rlg-lat2d(i,j)+dec
! qd+qdif = S_0 P sin(z) ( 1 - xhi + Tr) / 2
! rad = S_0 P sin(z) ( 1 - xhi + Tr) / 2 ( 1 - c_1 nu + c_2 beta) ( 1 - alpha )
         zrad=(qd+qdif)*(1.0_rlg-0.62_rlg*zfcc+0.0019_rlg*be)*(1.0_rlg-alb)

  END SUBROUTINE sflx_solar_luyten

!!======================================================================

  SUBROUTINE sflx_solar_empty(i,j,zfcc,zsinfiu,zcosfiu,zrad)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE sflx_solar_empty  ***
!&E
!&E ** Purpose : warning
!&E
!&E---------------------------------------------------------------------
!! * Modules used

!! * Declaration function

!! * Arguments
   INTEGER,        INTENT( IN )  :: i,j
   REAL(KIND=rsh), INTENT( IN )  :: zsinfiu, & !
                                    zcosfiu, & !
                                    zfcc       ! cloud cover fraction    [0<zfcc<1]
   REAL(KIND=rsh), INTENT( OUT ) :: zrad       ! solar flux

!! * Local declarations

!!----------------------------------------------------------------------
!! * Executable part

   zrad=0.0_rsh
   WRITE(ierrorlog,*) 'You do not read the solar atmospherical flux'
   WRITE(ierrorlog,*) 'and you have not chosen any formulation either'
   WRITE(ierrorlog,*) 'The simulation is stopped'
   STOP

  END SUBROUTINE sflx_solar_empty


!!======================================================================

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE sflx_ir  ***
!&E
!&E ** Purpose : long-wave heat flux estimation
!&E
!&E ** Description : The net upward flux of long-wave radiation from the
!&E              ocean is usually calculated as a product of (a) sigma T**4,
!&E              which is the flux emitted by a blackbody at the temperature T,
!&E              sigma being Stefan s constant, (b) a correction factor, 0.985
!&E              for departure of the ocean surface from blackbody behavior,
!&E              (c) a correction factor for back radiation in the absence of
!&E              cloud, and (d) a correction factor for cloud effect.
!&E
!&E
!&E ** Options :
!&E      - Clark et al, 1995
!&E      - Luyten et al, 1992
!&E      - Swimbank, dans thèse A. Agoumi
!&E      - Bigmani et al, 1995
!&E      - Berliand and Berliand, 1952
!&E      - Hastenrath and Lamb, 1978
!&E      - Ayina, 2008
!&E
!&E        see METEOCEAN report (2009)
!&E
!&E---------------------------------------------------------------------

!!======================================================================
 
  SUBROUTINE sflx_ir_swimbank(ztemp,zsat,zrh,zfcc,zirlos,zdirlos)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE sflx_ir_swimbank  ***
!&E
!&E ** Purpose : long-wave heat flux estimation
!&E
!&E ** Description :
!&E      The flux emitted by the ocean (corrected from blackbody
!&E      theory) is corrected from the atmosphere radiation and
!&E      cloud effect.
!&E
!&E      ztemp          -  sea surface temperature (degre C)
!&E      zsat           -  air temperature (degre C)
!&E      st             -  sea surface temperature (Kelvin)
!&E      at             -  air temperature (Kelvin)
!&E      em             -  emissivity  (0.985)
!&E      zirlos         -  infra-red heat flux [W/m2]
!&E
!&E ** Called by :
!&E
!&E ** External calls :
!&E
!&E ** Reference :     Swimbank (Agoumi thesis)
!&E                    METEOCEAN report (2009)
!&E
!&E ** History :
!&E       !          Original code
!&E       !  06-05  (V. Garnier) 0 < rh fcc < 1
!&E       !  2009-06  (V. Garnier) correction according to METEOCEAN report
!&E
!&E---------------------------------------------------------------------
!! * Modules used

!! * Declaration subroutine

!! * Arguments
   REAL(KIND=rsh), INTENT( IN )  :: ztemp, & ! sea surface temperature [C]
                                    zsat,  & ! air temperature         [C]
                                    zrh,   & ! relative humidity       [0<zrh<1]
                                    zfcc     ! cloud cover fraction    [0<zfcc<1]
   REAL(KIND=rsh), INTENT( OUT ) :: zirlos, & ! infra-red flux lost by the ocean
                                    zdirlos

!! * Local declarations
   REAL(KIND=rsh), PARAMETER  :: cst = 0.937e-5,    &
                                 landa = 0.3_rsh      ! cloud influence
   REAL(KIND=rsh)             :: st,at,ema

!!----------------------------------------------------------------------
!! * Executable part

! sea temperature (Kelvin)
         st=ztemp + kelvin
! air temperature (Kelvin)
         at=zsat + kelvin
! ??????????????????
         ema=cst*at**2

! thermal flux (infra-red)
         zirlos=em*stefan*at**4*ema*(1.0_rsh+landa*zfcc**2)
         zirlos=zirlos-em*stefan*st**4
         IF (l_rapflux) zdirlos=-4.0_rsh*em*stefan*st**3


  END SUBROUTINE sflx_ir_swimbank

!!======================================================================
 
  SUBROUTINE sflx_ir_luyten(ztemp,zsat,zrh,zfcc,zirlos,zdirlos)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE sflx_ir_luyten  ***
!&E
!&E ** Purpose : long-wave heat flux estimation
!&E
!&E ** Description :
!&E      The flux emitted by the ocean (corrected from blackbody
!&E      theory) is corrected from the back atmosphere radiation
!&E      in cloudless condition and from cloud effect.
!&E
!&E      ztemp          -  sea surface temperature (degre C)
!&E      zsat           -  air temperature (degre C)
!&E      st             -  sea surface temperature (Kelvin)
!&E      at             -  air temperature (Kelvin)
!&E      em             -  emissivity  (0.985)
!&E      ew_a           -  saturation vapour pressure
!&E      zirlos         -  infra-red heat flux [W/m2]
!&E
!&E ** Called by :
!&E
!&E ** External calls :
!&E
!&E ** Reference :  Luyten and De Mulder, 1992
!&E                 METEOCEAN report (2009)
!&E
!&E ** History :
!&E       !  2009-06  (V. Garnier) Original code according to METEOCEAN report
!&E
!&E---------------------------------------------------------------------
!! * Modules used

!! * Declaration subroutine

!! * Arguments
   REAL(KIND=rsh), INTENT( IN )  :: ztemp, & ! sea surface temperature [C]
                                    zsat,  & ! air temperature         [C]
                                    zrh,   & ! relative humidity       [0<zrh<1]
                                    zfcc     ! cloud cover fraction    [0<zfcc<1]
   REAL(KIND=rsh), INTENT( OUT ) :: zirlos, & ! infra-red flux lost by the ocean
                                    zdirlos

!! * Local declarations
   INTEGER                    :: i,j
   REAL(KIND=rsh), PARAMETER  :: cst1 = 0.39_rsh,    &
                                 cst2 = 0.05_rsh,    &
                                 landa = 0.6_rsh      ! cloud influence
   REAL(KIND=rsh)             :: st,at,ew_a

!!----------------------------------------------------------------------
!! * Executable part

! sea temperature (Kelvin)
         st=ztemp + kelvin
! air temperature (Celsius for ew_gill)
         at=zsat
! saturation vapor pressure of air
         ew_a=zrh*ew_gill(at)

! thermal flux (infra-red)
         zirlos=-em*stefan*st**4
         zirlos=zirlos * ( cst1-cst2*SQRT(ew_a) ) * (1.0_rsh-landa*zfcc**2)
         IF (l_rapflux) zdirlos=-4.0_rsh*em*stefan*st**3 *      &
                        ( cst1-cst2*SQRT(ew_a) ) * (1.0_rsh-landa*zfcc**2)

  END SUBROUTINE sflx_ir_luyten

!!======================================================================
 
  SUBROUTINE sflx_ir_clark(ztemp,zsat,zrh,zfcc,zirlos,zdirlos)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE sflx_ir_clark  ***
!&E
!&E ** Purpose : long-wave heat flux estimation
!&E
!&E ** Description :
!&E      The flux emitted by the ocean (corrected from blackbody
!&E      theory) is corrected from the back atmosphere radiation
!&E      in cloudless condition. The presence of cloud is not
!&E      taken into account.
!&E
!&E      ztemp          -  sea surface temperature (degre C)
!&E      zsat           -  air temperature (degre C)
!&E      st             -  sea surface temperature (Kelvin)
!&E      at             -  air temperature (Kelvin)
!&E      em             -  emissivity  (0.985)
!&E      ew_a           -  saturation vapour pressure
!&E      zirlos         -  infra-red heat flux [W/m2]
!&E
!&E ** Called by :
!&E
!&E ** External calls :
!&E
!&E ** Reference :  Clark et al, 1995
!&E                 METEOCEAN report (2009)
!&E
!&E ** History :
!&E       !  2009-06  (V. Garnier) Original code according to METEOCEAN report
!&E
!&E---------------------------------------------------------------------
!! * Modules used

!! * Declaration subroutine

!! * Arguments
   REAL(KIND=rsh), INTENT( IN )  :: ztemp, & ! sea surface temperature [C]
                                    zsat,  & ! air temperature         [C]
                                    zrh,   & ! relative humidity       [0<zrh<1]
                                    zfcc     ! cloud cover fraction    [0<zfcc<1]
   REAL(KIND=rsh), INTENT( OUT ) :: zirlos, & ! infra-red flux lost by the ocean
                                    zdirlos

!! * Local declarations
   INTEGER                    :: i,j
   REAL(KIND=rsh), PARAMETER  :: cst1 = 0.39_rsh,    &
                                 cst2 = 0.05_rsh
   REAL(KIND=rsh)             :: st,at,ew_a

!!----------------------------------------------------------------------
!! * Executable part

! sea temperature (Kelvin)
         st=ztemp + kelvin
! air temperature (Celsius for ew_gill)
         at=zsat
! saturation vapor pressure of air
         ew_a=zrh*ew_gill(at)

! thermal flux (infra-red)
         zirlos=-em*stefan*st**4
         zirlos=zirlos * ( cst1-cst2*SQRT(ew_a) )
         IF (l_rapflux) zdirlos=-4.0_rsh*em*stefan*st**3 *      &
                        ( cst1-cst2*SQRT(ew_a) )

  END SUBROUTINE sflx_ir_clark

!!======================================================================
 
  SUBROUTINE sflx_ir_berliand(ztemp,zsat,zrh,zfcc,zirlos,zdirlos)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE sflx_ir_berliand  ***
!&E
!&E ** Purpose : long-wave heat flux estimation
!&E
!&E ** Description :
!&E      The flux emitted by the ocean (corrected from blackbody
!&E      theory) is corrected from the back atmosphere radiation
!&E      in cloudless condition, from the cloud effect, and from
!&E      the air-sea temperature gradient in atmospheric layer
!&E      boundary.
!&E
!&E      ztemp          -  sea surface temperature (degre C)
!&E      zsat           -  air temperature (degre C)
!&E      st             -  sea surface temperature (Kelvin)
!&E      at             -  air temperature (Kelvin)
!&E      em             -  emissivity  (0.985)
!&E      ew_a           -  saturation vapour pressure
!&E      zirlos         -  infra-red heat flux [W/m2]
!&E
!&E ** Called by :
!&E
!&E ** External calls :
!&E
!&E ** Reference :  Berliand and Berliand (1952)
!&E                 METEOCEAN report (2009)
!&E
!&E ** History :
!&E       !  2009-06  (V. Garnier) Original code according to METEOCEAN report
!&E
!&E---------------------------------------------------------------------
!! * Modules used

!! * Declaration subroutine

!! * Arguments
   REAL(KIND=rsh), INTENT( IN )  :: ztemp, & ! sea surface temperature [C]
                                    zsat,  & ! air temperature         [C]
                                    zrh,   & ! relative humidity       [0<zrh<1]
                                    zfcc     ! cloud cover fraction    [0<zfcc<1]
   REAL(KIND=rsh), INTENT( OUT ) :: zirlos, & ! infra-red flux lost by the ocean
                                    zdirlos

!! * Local declarations
   INTEGER                    :: i,j
   REAL(KIND=rsh), PARAMETER  :: cst1 = 0.39_rsh,    &
                                 cst2 = 0.05_rsh,    &
                                 cst3 = 4.00_rsh,    &
                                 landa = 0.6_rsh      ! cloud influence
   REAL(KIND=rsh)             :: st,at,ew_a

!!----------------------------------------------------------------------
!! * Executable part

! sea temperature (Kelvin)
         st=ztemp + kelvin
! air temperature
         at=zsat + kelvin
! saturation vapor pressure of air
         ew_a=zrh*ew_gill(zsat)   ! Celsius for ew_gill

! thermal flux (infra-red)
         zirlos=-em*stefan*st**4
         zirlos=zirlos * ( cst1-cst2*SQRT(ew_a) ) * (1.0_rsh-landa*zfcc**2)
         zirlos=zirlos - cst3*em*stefan*(st-at)*st**3
         IF (l_rapflux) zdirlos=-4.0_rsh*em*stefan*st**3 *                  &
                        ( cst1-cst2*SQRT(ew_a) ) * (1.0_rsh-landa*zfcc**2)  &
                        - 16.0_rsh*em*stefan*st**3 - 12.0_rsh*em*stefan*at*st**2 

  END SUBROUTINE sflx_ir_berliand
 
!!======================================================================
 
  SUBROUTINE sflx_ir_hastenrath(ztemp,zsat,zrh,zfcc,zirlos,zdirlos)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE sflx_ir_hastenrath  ***
!&E
!&E ** Purpose : long-wave heat flux estimation
!&E
!&E ** Description :
!&E      The flux absorbed by the ocean but emitted by the atmosphere
!&E      in cloudless condition, then corrected from the cloud effect
!&E      and from the air-sea temperature gradient in atmospheric layer
!&E      boundary.
!&E
!&E      ztemp          -  sea surface temperature (degre C)
!&E      zsat           -  air temperature (degre C)
!&E      st             -  sea surface temperature (Kelvin)
!&E      at             -  air temperature (Kelvin)
!&E      em             -  emissivity  (0.985)
!&E      ew_a           -  saturation vapour pressure
!&E      zirlos         -  infra-red heat flux [W/m2]
!&E
!&E ** Called by :
!&E
!&E ** External calls :
!&E
!&E ** Reference :  Hastenrath and Lamb (1978)
!&E                 METEOCEAN report (2009)
!&E
!&E ** History :
!&E       !  2009-06  (V. Garnier) Original code according to METEOCEAN report
!&E
!&E---------------------------------------------------------------------
!! * Modules used

!! * Declaration subroutine

!! * Arguments
   REAL(KIND=rsh), INTENT( IN )  :: ztemp, & ! sea surface temperature [C]
                                    zsat,  & ! air temperature         [C]
                                    zrh,   & ! relative humidity       [0<zrh<1]
                                    zfcc     ! cloud cover fraction    [0<zfcc<1]
   REAL(KIND=rsh), INTENT( OUT ) :: zirlos, & ! infra-red flux lost by the ocean
                                    zdirlos

!! * Local declarations
   INTEGER                    :: i,j
   REAL(KIND=rsh), PARAMETER  :: cst1 = 0.39_rsh,    &
                                 cst2 = 0.05_rsh,    &
                                 cst3 = 4.00_rsh,    &
                                 landa = 0.6_rsh      ! cloud influence
   REAL(KIND=rsh)             :: st,at,ew_a

!!----------------------------------------------------------------------
!! * Executable part

! sea temperature (Kelvin)
         st=ztemp + kelvin
! air temperature
         at=zsat + kelvin
! saturation vapor pressure of air
         ew_a=zrh*ew_gill(zsat)   ! Celsius for ew_gill

! thermal flux (infra-red)
         zirlos=-em*stefan*at**4
         zirlos=zirlos * ( cst1-cst2*SQRT(ew_a) ) * (1.0_rsh-landa*zfcc**2)
         zirlos=zirlos - cst3*em*stefan*(st-at)*at**3
         IF (l_rapflux) zdirlos=-cst3*em*stefan*at**3

  END SUBROUTINE sflx_ir_hastenrath


!!======================================================================
 
  SUBROUTINE sflx_ir_bignami(ztemp,zsat,zrh,zfcc,zirlos,zdirlos)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE sflx_ir_bignami  ***
!&E
!&E ** Purpose : long-wave heat flux estimation
!&E
!&E ** Description :
!&E      The flux emitted by the ocean (corrected from blackbody
!&E      theory) is corrected from the back atmosphere flux (radiation
!&E      in cloudless condition and cloud effect).
!&E
!&E      ztemp          -  sea surface temperature (degre C)
!&E      zsat           -  air temperature (degre C)
!&E      st             -  sea surface temperature (Kelvin)
!&E      at             -  air temperature (Kelvin)
!&E      em             -  emissivity  (0.985)
!&E      ew_a           -  saturation vapour pressure
!&E      zirlos         -  infra-red heat flux [W/m2]
!&E
!&E ** Called by :
!&E
!&E ** External calls :
!&E
!&E ** Reference :  Bignami et al (1995)
!&E                 METEOCEAN report (2009)
!&E
!&E ** History :
!&E       !  2009-06  (V. Garnier) Original code according to METEOCEAN report
!&E
!&E---------------------------------------------------------------------
!! * Modules used

!! * Declaration subroutine

!! * Arguments
   REAL(KIND=rsh), INTENT( IN )  :: ztemp, & ! sea surface temperature [C]
                                    zsat,  & ! air temperature         [C]
                                    zrh,   & ! relative humidity       [0<zrh<1]
                                    zfcc     ! cloud cover fraction    [0<zfcc<1]
   REAL(KIND=rsh), INTENT( OUT ) :: zirlos, & ! infra-red flux lost by the ocean
                                    zdirlos

!! * Local declarations
   INTEGER                    :: i,j
   REAL(KIND=rsh), PARAMETER  :: cst1 = 0.653_rsh,    &
                                 cst2 = 0.00535_rsh,  &
                                 landa = 0.1762_rsh      ! cloud influence
   REAL(KIND=rsh)             :: st,at,ew_a

!!----------------------------------------------------------------------
!! * Executable part

! sea temperature (Kelvin)
         st=ztemp + kelvin
! air temperature
         at=zsat + kelvin
! saturation vapor pressure of air
         ew_a=zrh*ew_gill(zsat)   ! Celsius for ew_gill

! thermal flux (infra-red)
         zirlos=stefan*at**4 * ( cst1+cst2*ew_a ) * (1.0_rsh+landa*zfcc**2)
         zirlos=zirlos - em*stefan*st**4
         IF (l_rapflux) zdirlos=-4.0_rsh*em*stefan*st**3

  END SUBROUTINE sflx_ir_bignami

!!======================================================================
 
  SUBROUTINE sflx_ir_empty(ztemp,zsat,zrh,zfcc,zirlos,zdirlos)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE sflx_ir_empty  ***
!&E
!&E ** Purpose : warning
!&E
!&E---------------------------------------------------------------------
!! * Modules used

!! * Declaration function

!! * Arguments
   REAL(KIND=rsh), INTENT( IN )  :: ztemp, & ! sea surface temperature [C]
                                    zsat,  & ! air temperature         [C]
                                    zrh,   & ! relative humidity       [0<zrh<1]
                                    zfcc     ! cloud cover fraction    [0<zfcc<1]
   REAL(KIND=rsh), INTENT( OUT ) :: zirlos, & ! infra-red flux lost by the ocean
                                    zdirlos

!! * Local declarations

!!----------------------------------------------------------------------
!! * Executable part

   zirlos=0.0_rsh
   zdirlos=0.0_rsh

   WRITE(ierrorlog,*) 'You do not read the thermal (infra-red) atmospherical flux'
   WRITE(ierrorlog,*) 'and you have not chosen any formulation either'
   WRITE(ierrorlog,*) 'The simulation is stopped'
   STOP

  END SUBROUTINE sflx_ir_empty

!!======================================================================

  FUNCTION ew_gill(tcel)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE ew_gill  ***
!&E
!&E ** Purpose : estimate of the saturation vapor pressure
!&E
!&E ** Description : For temperatures t between +- 40°C, the value in millibars
!&E                  of the saturation vapor pressure is given correct to
!&E                  1 part in 500 by
!&E                  log_10 ew(t) = (0.7859 + 0.03477t) / (1 + 0.00412t)
!&E
!&E                  tcel - temperature in degrees Celsius
!&E                  ew   - saturation vapor pressure
!&E
!&E ** Called by : sflx_ir_*
!&E
!&E ** External calls :
!&E
!&E ** Reference :  Gill, 1982
!&E
!&E ** History :
!&E       !  2009-06  (V. Garnier)  Original code
!&E
!&E---------------------------------------------------------------------
!! * Modules used

!! * Declaration function

!! * Arguments
   REAL(KIND=rsh), INTENT( IN ) :: tcel           ! temperature

!! * Local declarations
   REAL(KIND=rsh), PARAMETER    :: cst1 = 0.7859_rsh,     &
                                   cst2 = 0.03477_rsh,    &
                                   cst3 = 0.00412_rsh
   REAL(KIND=rsh)               :: ew_gill        ! saturation vapor pressure

!!----------------------------------------------------------------------
!! * Executable part

   ew_gill = 10.0_rsh**( (cst1 + cst2*tcel) / (1.0_rsh + cst3*tcel) )

  END FUNCTION ew_gill

!!======================================================================

  FUNCTION sh_gill(tcel,rh,pr)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE sh_gill  ***
!&E
!&E ** Purpose : estimate of the saturation vapor pressure
!&E
!&E ** Description : The specific humidity q (or sh) is the mass of vapor per unit
!&E                  mass of moist air and is function of the vapor
!&E                  pressure e and the pressure p, according to the
!&E                  relation :
!&E                  e / p = q / (mmwar + (1-mmwar)q )
!&E                  with mmwar = 0.62197
!&E                            = molecular mass of water /
!&E                              molecular mass of air
!&E
!&E                  tcel - temperature in degrees Celsius
!&E                  pr   - atmospherical pressure in millibars
!&E                  ew   - saturation vapor pressure
!&E                  rh   - relative humidity
!&E                  sh   - specific humidity
!&E
!&E ** Called by : sflx_ir_*
!&E
!&E ** External calls :
!&E
!&E ** Reference :  Gill, 1982
!&E
!&E ** History :
!&E       !  2009-06  (V. Garnier)  Original code
!&E
!&E---------------------------------------------------------------------
!! * Modules used

!! * Declaration function

!! * Arguments
   REAL(KIND=rsh), INTENT( IN ) :: tcel,     &  ! temperature (degrees Celcius)
                                   rh,       &  ! relative humidity (0 < rh < 1)
                                   pr           ! pressure (mbar)

!! * Local declarations
   REAL(KIND=rsh), PARAMETER    :: cst1 = 0.7859_rsh,     &
                                   cst2 = 0.03477_rsh,    &
                                   cst3 = 0.00412_rsh,    &
                                   mmwar = 0.62197_rsh
   REAL(KIND=rsh)               :: ew,       &  ! saturation vapor pressure
                                   sh_gill      ! specific humidity

!!----------------------------------------------------------------------
!! * Executable part

   ew = rh*10.0_rsh**( (cst1 + cst2*tcel) / (1.0_rsh + cst3*tcel) )
   sh_gill = mmwar*ew / ( pr - (1.0_rsh-mmwar)*ew )

  END FUNCTION sh_gill

!!======================================================================

 END MODULE sflxrad
