MODULE hdo_surfex_mod

      IMPLICIT NONE

      CONTAINS

      subroutine hdo_surfex(ngrid,nlay,nq,ptimestep,
     &                      zt,pplay,zq,pqsurf,
     &                      old_h2o_vap,qsat,pdqsdif,dwatercap_dif,
     &                      hdoflux)

      use tracer_mod, only: igcm_h2o_vap, igcm_h2o_ice,
     &                      igcm_hdo_vap, igcm_hdo_ice,
     &                      qperemin
      use surfdat_h, only: watercaptag
      use geometry_mod, only: longitude_deg,latitude_deg
      use comcstfi_h, only: pi
      use microphys_h, only: nav, kbz, mh2o, mco2, mhdo
      use microphys_h, only: molco2, molh2o, molhdo
      use write_output_mod, only: write_output

      implicit none
c------------------------------------------------------------------
c               Routine to compute the fluxes between air and surface
c               for HDO, based of the fluxes for H2O
c           L. Rossi.; M. Vals 2019
c------------------------------------------------------------------
      include "callkeys.h"
c------------------------------------------------------------------
c     Arguments:
c     ---------
c     Inputs:
      INTEGER, INTENT(IN) :: ngrid,nlay
      INTEGER, INTENT(IN) :: nq                 ! nombre de traceurs

      REAL, INTENT(IN) :: ptimestep             ! pas de temps physique (s)
      REAL, INTENT(IN) :: zt(ngrid,nlay)       ! local value of temperature
      REAL, INTENT(IN) :: pplay(ngrid,nlay)     ! pression au milieu des couches (Pa) 
      REAL, INTENT(IN) :: zq(ngrid,nlay,nq)    ! local value of tracers
      REAL, INTENT(IN) :: pqsurf(ngrid,nq)
      REAL, INTENT(IN) :: old_h2o_vap(ngrid)     ! traceur d'eau avant
                                           !traitement de l'eau (kg/kg)
      REAL, INTENT(IN) :: qsat(ngrid) ! saturation mixing ratio
      REAL, INTENT(IN) :: dwatercap_dif(ngrid)  ! trend related to permanent ice
      REAL, INTENT(INOUT) :: pdqsdif(ngrid,nq)    ! tendance towards surface 
                                 !   (kg/kg.s-1)

c     Output:
      REAL, INTENT(OUT) :: hdoflux(ngrid)       ! value of vapour flux of HDO

c------------------------------------------------------------------
c     Local variables:

      REAL alpha(ngrid)    ! equilibrium fractionation factor
      REAL alpha_c(ngrid)  ! real fractionation factor
      REAL extrasublim ! sublimation in excess of surface ice
      REAL tmpratio(ngrid)   ! D/H ratio in flux to surf
      REAL h2oflux(ngrid)       ! value of vapour flux of H2O
                                      ! same sign as pdqsdif
      REAL*8 satu(ngrid)          ! Water vapor saturation ratio over ice
      REAL zt1(ngrid),pplay1(ngrid)
      REAL Dv,Dv_hdo ! Water/HDO vapor diffusion coefficient
      INTEGER ig,l

      REAL DoH_vap(ngrid)

c-----------------------------------------------------------------------
c    Calculation of the fluxes for HDO
        !! Calculation of the saturation ratio in the layer above the surface
        satu(1:ngrid)=old_h2o_vap(1:ngrid) / qsat(1:ngrid)
        !! Initialisation of the fractionation coefficient
        alpha(1:ngrid)=1.        
        alpha_c(1:ngrid)=1.

        DO ig=1,ngrid
              
            h2oflux(ig) = pdqsdif(ig,igcm_h2o_ice) + 
     &          dwatercap_dif(ig)

            !! IF Sublimation
            if (h2oflux(ig).le.0.) then

               if (pqsurf(ig,igcm_h2o_ice).gt.qperemin) then
                pdqsdif(ig,igcm_hdo_ice) =
     &            pdqsdif(ig,igcm_h2o_ice)*
     &             (pqsurf(ig,igcm_hdo_ice)/
     &             pqsurf(ig,igcm_h2o_ice) )
               else
                pdqsdif(ig,igcm_hdo_ice) = 0.
               endif

                pdqsdif(ig,igcm_hdo_ice)=
     &             max(pdqsdif(ig,igcm_hdo_ice),
     &            -pqsurf(ig,igcm_hdo_ice)/ptimestep)

                hdoflux(ig) = pdqsdif(ig,igcm_hdo_ice)

             if(watercaptag(ig)) then

              !if we sublimate more than qsurf
              if ((-h2oflux(ig)*ptimestep)
     &           .gt.pqsurf(ig,igcm_h2o_ice)) then

C               dwatercap_dif is how much we sublimate in excess of
C               pqsurf for H2O                        
C               hdoflux(ig) is the flux of HDO from atm. to surf.
c               The D/H of the old ice is supposed to be 5 SMOW
c               We need D/H of the flux to be 5, so we need
c               dwatercap_dif* 5 * 2 * 155.76e-6 (=1 SMOW)
                    hdoflux(ig)= hdoflux(ig)
     &                   +(dwatercap_dif(ig)*(2.*155.76e-6)*5.)
                endif
             endif ! watercap

            else ! condensation

               if (hdofrac) then !do we use fractionation?
                !! Calculation of the H2O vapor diffusion coefficient
                Dv = 1./3. * sqrt( 8*kbz*zt(ig,1)/(pi*mh2o/nav) )
     &             * kbz * zt(ig,1) / 
     &            ( pi * pplay(ig,1) * (molco2+molh2o)*(molco2+molh2o) 
     &             * sqrt(1.+mh2o/mco2) )
                !! Calculation of the HDO vapor diffusion coefficient
                Dv_hdo = 1./3. * sqrt( 8*kbz*zt(ig,1)/(pi*mhdo/nav) )
     &             * kbz * zt(ig,1) / 
     &            ( pi * pplay(ig,1) * (molco2+molhdo)*(molco2+molhdo) 
     &             * sqrt(1.+mhdo/mco2) )
                !! Calculation of the "equilibrium" fractionation coefficient
c               alpha(ig) = exp(16288./zt(ig,1)**2.-9.34e-2)
                alpha(ig) = exp(13525./zt(ig,1)**2.-5.59e-2) !Lamb
                !! Calculation of the 'real' fractionnation coefficient (effect of kinetics, see Jouzel and Merlivat, 1984)
                alpha_c(ig) = (alpha(ig)*satu(ig))/
     &             ( (alpha(ig)*(Dv/Dv_hdo)*(satu(ig)-1.)) + 1.)
               else
                alpha_c(ig) = 1.
               endif

               if (old_h2o_vap(ig).gt.qperemin) then
                         pdqsdif(ig,igcm_hdo_ice)=
     &                      alpha_c(ig)*pdqsdif(ig,igcm_h2o_ice)*
     &                      (zq(ig,1,igcm_hdo_vap)/
     &                          old_h2o_vap(ig))
                else
                 pdqsdif(ig,igcm_hdo_ice)= 0.
                endif

               if (hdofrac) then !do we use fractionation?
                pdqsdif(ig,igcm_hdo_ice)= 
     &      min(  pdqsdif(ig,igcm_hdo_ice), 
     &           (zq(ig,1,igcm_hdo_vap)/ptimestep) )
               endif

                hdoflux(ig)=pdqsdif(ig,igcm_hdo_ice)

            endif !sublimation

        ENDDO ! of DO ig=1,ngrid

c           CALL write_output('extrasublim',
c    &                       'extrasublimation',
c    &                       ' ',tmpratio)
c           CALL write_output('alpha_c_s',
c    &                       'alpha_c_s',
c    &                       ' ',alpha_c)

       return
      end subroutine hdo_surfex
c------------------------------------------------------------------

      end module hdo_surfex_mod