MODULE qsat_mod




CONTAINS

  SUBROUTINE qsat(f_temp,f_p,f_qs)
  USE icosa
  IMPLICIT NONE
    TYPE(t_field), POINTER :: f_temp(:)  ! IN    : temperature
    TYPE(t_field), POINTER :: f_p(:)     ! IN    : pressure at mid-levels
    TYPE(t_field), POINTER :: f_qs(:)    ! OUT   : vapeur d'eau saturante en kg/kg
    
    REAL(rstd),POINTER :: temp(:,:),  p(:,:), qs(:,:)
    INTEGER :: ind

    DO ind=1,ndomain
       IF (.NOT. assigned_domain(ind)) CYCLE
       CALL swap_dimensions(ind)
       CALL swap_geometry(ind)
       temp=f_temp(ind)
       p=f_p(ind)
       qs=f_qs(ind)
       CALL compute_qsat(temp,p,qs)
    END DO
  
  END SUBROUTINE qsat

  
  SUBROUTINE compute_qsat(temp,p,qsat)
  USE icosa
  USE omp_para
  IMPLICIT NONE
  
!======================================================================
! Autheur(s): Z.X. Li (LMD/CNRS)
!  reecriture vectorisee par F. Hourdin.
! Objet: calculer la vapeur d'eau saturante (formule Centre Euro.)
!======================================================================
! Arguments:
! kelvin---input-R: temperature en Kelvin
! millibar--input-R: pression en mb
!
! q_sat----output-R: vapeur d'eau saturante en kg/kg
!======================================================================
!
      REAL,INTENT(IN)  :: temp(iim*jjm,llm)
      REAL,INTENT(IN)  :: p   (iim*jjm,llm+1)
      REAL,INTENT(OUT) :: qsat(iim*jjm,llm)

      REAL, PARAMETER  ::  r2es=611.14 *18.0153/28.9644

      REAL :: r3es
      REAL, PARAMETER ::  r3les=17.269
      REAL, PARAMETER ::  r3ies=21.875

      REAL :: r4es
      REAL, PARAMETER :: r4les=35.86
      REAL, PARAMETER :: r4ies=7.66

      REAL, PARAMETER :: rtt=273.16
      REAL, PARAMETER :: retv=28.9644/18.0153 - 1.0

      REAL :: zqsat, pmid
      INTEGER :: l,i,j,ij
!
!     ------------------------------------------------------------------
!
!

    DO l=ll_begin,ll_end
      DO j=jj_begin,jj_end
        DO i=ii_begin,ii_end
          ij=(j-1)*iim+i

          IF (temp(ij,l) .LE. rtt) THEN
            r3es = r3ies
            r4es = r4ies
          ELSE
            r3es = r3les
            r4es = r4les
          ENDIF
          pmid=0.5*(p(ij,l)+p(ij,l+1))
          zqsat=r2es/pmid*EXP(r3es*(temp(ij,l)-rtt)/(temp(ij,l)-r4es))
          zqsat=MIN(0.5,zqsat)
          zqsat=zqsat/(1.-retv *zqsat)

          qsat(ij,l)= zqsat

        ENDDO
      ENDDO
    ENDDO
  
  
  END SUBROUTINE compute_qsat

END MODULE qsat_mod