MODULE watersat_mod
 
!======================================================================================================================!
! Subject:
!---------
!   Module used to compute the saturation water vapor pressure and qsat
!----------------------------------------------------------------------------------------------------------------------!
! Reference:
!-----------
!  Laurent Li Parameterization documentation. Checked again Murphy and Koop 2005 (QJRMS)
!
!======================================================================================================================!


IMPLICIT NONE

CONTAINS
      
      SUBROUTINE watersat(gridsize,T,P,qsat)
      
      IMPLICIT NONE

!=======================================================================
!
!  Water mass mixing ratio at saturation (kg/kg) for a given pressure (Pa)
!  and Temperature (K) array of dimension gridsize
!
!  Moved in module by Deborah BARDET 02/07/18
!  Adapted in F90 by LL 09/02/2024
!=======================================================================

!   declarations:
!   -------------
!   arguments:
!   ----------

!   INPUT
      INTEGER, INTENT(in) :: gridsize              ! Dimension of the array
      REAL, INTENT(in) :: T(gridsize)    ! Temperature [K]
      REAL, INTENT(in) :: p(gridsize)              ! Pressure [Pa]
      
!   OUTPUT
      real, INTENT(out) :: qsat(gridsize)          ! qsat for water vapor [kg/kg]
      
!   local:
!   ------
      INTEGER :: i                                 ! Loop index
      REAL :: psat                                 ! Saturation water vapor over ice [Pa]
      REAL :: epsi                                 ! Rapport of water molar mass over CO2 molar mass [1]  
         
!    Code:
!    -------- 
      epsi = 18./44. !~0.41, as in other version. Perhaps should be recomputed properly is mu_dry changes ?
      do i=1,gridsize
          psat = 100.0 * 10**(2.07023 - 0.00320991*T(i) - 2484.896 / T(i) + 3.56654 * log10(T(i)))
          if(psat.gt.P(i)) then
              qsat(i) = 1.
          else
              qsat(i) = epsi*psat/(P(i)-(1.-epsi)*psat)
              qsat(i) = max(qsat(i), 1.e-30) ! in case of tiny value
          endif
      enddo

      END SUBROUTINE watersat
      
      END MODULE watersat_mod