mode_psychro.F90

Go to the documentation of this file.

!SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
!SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
!SFX_LIC for details. version 1.
!##################
MODULE mode_psychro
!##################
!
!!****  *MODE_PSYCHRO* -
!!
!!    PURPOSE
!!    -------
!      
!
!!
!!**  IMPLICIT ARGUMENTS
!!    ------------------
!!       NONE          
!!
!!    REFERENCE
!!    ---------
!!
!!
!!    AUTHOR
!!    ------
!!      
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    12/04/11
!!      J.Escobar   11/13 :  remove space in ELSEWHERE statement
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
interface pe_from_pq
        module procedure pe_from_pq_0d
        module procedure pe_from_pq_1d
end interface
interface td_from_tq
        module procedure td_from_tq_0d
        module procedure td_from_tq_1d
end interface
interface rv_from_tptwb
        module procedure rv_from_tptwb_0d
        module procedure rv_from_tptwb_1d
end interface
interface twb_from_tpq
        module procedure twb_from_tpq_0d
        module procedure twb_from_tpq_1d
end interface
INTERFACE enth_fn_t_q
  MODULE PROCEDURE enth_fn_t_q
END INTERFACE
INTERFACE q_fn_t_enth
  MODULE PROCEDURE q_fn_t_enth
END INTERFACE

contains
!PE_FROM_PQ
!----------
function pe_from_pq_0d(PP, PQ) RESULT(PE) 
!arguments and result
REAL, INTENT(IN) :: pp !atmos. pressure (Pa)
REAL, INTENT(IN) :: pq !specific humidity (kg/kg)
REAL :: pe !water vapour pressure (Pa)
REAL(KIND=JPRB) :: zhook_handle
IF (lhook) CALL dr_hook('MODE_PSYCHRO:PE_FROM_PQ_0D',0,zhook_handle)
pe = pq * pp /(0.622 + 0.378 * pq)
IF (lhook) CALL dr_hook('MODE_PSYCHRO:PE_FROM_PQ_0D',1,zhook_handle)
end function pe_from_pq_0d

function pe_from_pq_1d(PP, PQ) RESULT(PE) 
!arguments and result
REAL, DIMENSION(:), INTENT(IN) :: pp !atmos. pressure (Pa)
REAL, DIMENSION(:), INTENT(IN) :: pq !specific humidity (kg/kg)
REAL, DIMENSION(SIZE(PQ))      :: pe !water vapour pressure (Pa)
REAL(KIND=JPRB) :: zhook_handle
IF (lhook) CALL dr_hook('MODE_PSYCHRO:PE_FROM_PQ_1D',0,zhook_handle)
pe(:) = pq(:) * pp(:) /(0.622 + 0.378 * pq(:))
IF (lhook) CALL dr_hook('MODE_PSYCHRO:PE_FROM_PQ_1D',1,zhook_handle)
end function pe_from_pq_1d
!-------------------------

!TD_FROM_TQ
function td_from_tq_0d(PT, PQ) RESULT(PTD)
USE modd_csts
USE modd_surf_par, ONLY: xundef
!arguments and result
REAL, INTENT(IN) :: pt !Air Temp. (K)
REAL, INTENT(IN) :: pq !Specific humidity (kg/kg)
REAL :: ptd !Dew Point Air Temp. (K)
!local variables
REAL :: alpha
REAL :: zpe !water vapour pressure
REAL(KIND=JPRB) :: zhook_handle

IF (lhook) CALL dr_hook('MODE_PSYCHRO:TD_FROM_TQ_0D',0,zhook_handle)
zpe = pe_from_pq(pt, pq)
alpha = log(zpe/1000.)
IF (pt .GE. xtt .AND. pt .GE. 93.+xtt) THEN
        ptd = xtt+6.54+14.526*alpha+0.7389*alpha*alpha+0.09486*alpha**3 &
              +0.4569*(zpe/1000.)**0.1984
ELSE IF (pt .LT. xtt) THEN
        ptd = xtt+6.09+12.608*alpha+0.4959*alpha*alpha
ELSE
        ptd = xundef
ENDIF
ptd = min(ptd, pt)
IF (lhook) CALL dr_hook('MODE_PSYCHRO:TD_FROM_TQ_0D',1,zhook_handle)
end function td_from_tq_0d

function td_from_tq_1d(PT, PQ) RESULT(PTD)
USE modd_csts
USE modd_surf_par, ONLY: xundef
!arguments and result
REAL, DIMENSION(:), INTENT(IN) :: pt !Air Temp. (K)
REAL, DIMENSION(:), INTENT(IN) :: pq !Specific humidity (kg/kg)
REAL, DIMENSION(SIZE(PQ))      :: ptd !Dew Point Air Temp. (K)
!local variables
REAL, DIMENSION(SIZE(PQ)) :: alpha
REAL, DIMENSION(SIZE(PQ)) :: zpe !water vapour pressure
REAL(KIND=JPRB) :: zhook_handle

IF (lhook) CALL dr_hook('MODE_PSYCHRO:TD_FROM_TQ_1D',0,zhook_handle)
zpe = pe_from_pq(pt, pq)
alpha(:) = log(zpe(:)/1000.)
WHERE (pt .GE. xtt .AND. pt .GE. 93.+xtt)
        ptd = xtt+6.54+14.526*alpha+0.7389*alpha*alpha+0.09486*alpha**3 &
              +0.4569*(zpe/1000.)**0.1984
      ELSEWHERE (pt .LT. xtt)
        ptd = xtt+6.09+12.608*alpha+0.4959*alpha*alpha
ELSEWHERE
        ptd = xundef
END WHERE
ptd(:) = min(ptd(:), pt(:))
IF (lhook) CALL dr_hook('MODE_PSYCHRO:TD_FROM_TQ_1D',1,zhook_handle)
end function td_from_tq_1d
!-------------------------

!RV_FROM_TPTWB
function rv_from_tptwb_0d(PT, PP, PTWB) RESULT(PRV)
USE mode_thermos
USE modd_csts
!arguments and result
REAL, INTENT(IN) :: pt !Air temperature (K)
REAL, INTENT(IN) :: pp !Atmos. Pressure (Pa)
REAL, INTENT(IN) :: ptwb !Wet Bulb Temp. (K)
REAL :: prv !water vapor mixing ratio (kg/kg)
REAL :: zrvsat !saturation water vapor mixing ratio
REAL(KIND=JPRB) :: zhook_handle

IF (lhook) CALL dr_hook('MODE_PSYCHRO:RV_FROM_TPTWB_0D',0,zhook_handle)
zrvsat = qsat(pt, pp) / (1 - qsat(pt, pp))
prv = ((2501. - 2.326*(ptwb-xtt))*zrvsat - 1.006*(pt - ptwb)) &
       / (2501. + 1.86*(pt - xtt) -4.186*(ptwb - xtt))
IF (lhook) CALL dr_hook('MODE_PSYCHRO:RV_FROM_TPTWB_0D',1,zhook_handle)
end function rv_from_tptwb_0d

function rv_from_tptwb_1d(PT, PP, PTWB) RESULT(PRV)
USE mode_thermos
USE modd_csts
!arguments and result
REAL, DIMENSION(:), INTENT(IN) :: pt !Air temperature (K)
REAL, DIMENSION(:),INTENT(IN) :: pp !Atmos. Pressure (Pa)
REAL, DIMENSION(:),INTENT(IN) :: ptwb !Wet Bulb Temp. (K)
REAL, DIMENSION(SIZE(PT)) :: prv !water vapor mixing ratio (kg/kg)
REAL, DIMENSION(SIZE(PT)) :: zrvsat !saturation water vapor mixing ratio
REAL(KIND=JPRB) :: zhook_handle

IF (lhook) CALL dr_hook('MODE_PSYCHRO:RV_FROM_TPTWB_1D',0,zhook_handle)
zrvsat = qsat(pt, pp) / (1 - qsat(pt, pp))
prv(:) = ((2501. - 2.326*(ptwb(:)-xtt))*zrvsat(:) - 1.006*(pt(:) - ptwb(:))) &
       / (2501. + 1.86*(pt(:) - xtt) -4.186*(ptwb(:) - xtt))
IF (lhook) CALL dr_hook('MODE_PSYCHRO:RV_FROM_TPTWB_1D',1,zhook_handle)       
end function rv_from_tptwb_1d
!----------------------------

!TWB_FROM_TPQ
!------------
function twb_from_tpq_0d(PT, PP, PQ) RESULT(PTWB)
!arguments and results
REAL, INTENT(IN) :: pt !air temperature (K)
REAL, INTENT(IN) :: pq !mixing ratio (kg/kg)
REAL, INTENT(IN) :: pp !atmos. pressure (Pa)
REAL :: ptwb !Wet Bulb Temp. (K)
!local variable
REAL :: ztd !Dew Point Temp. (K)
REAL :: ztwbinf, ztwbsup, zrv
INTEGER :: jiter
REAL(KIND=JPRB) :: zhook_handle
IF (lhook) CALL dr_hook('MODE_PSYCHRO:TWB_FROM_TPQ_0D',0,zhook_handle)
jiter = 1
ztd = td_from_tq(pt, pq) 
!initial guess
ztwbsup = pt
ztwbinf = ztd
ptwb = 0.5 * (ztwbsup + ztwbinf)
DO WHILE (ztwbsup - ztwbinf > 0.001 .OR. jiter .LE. 50)
   zrv = rv_from_tptwb(pt, pp, ptwb)
   IF (zrv .GT. pq/(1 - pq)) THEN
           ztwbsup = ptwb
   ELSE
           ztwbinf = ptwb
   ENDIF
   ptwb = 0.5 * (ztwbinf + ztwbsup)
   jiter = jiter + 1
END DO
IF (lhook) CALL dr_hook('MODE_PSYCHRO:TWB_FROM_TPQ_0D',1,zhook_handle)
end function twb_from_tpq_0d

function twb_from_tpq_1d(PT, PP, PQ) RESULT(PTWB)
!arguments and results
REAL, DIMENSION(:), INTENT(IN) :: pt !air temperature (K)
REAL, DIMENSION(:), INTENT(IN) :: pq !humidity content (kg/kg)
REAL, DIMENSION(:), INTENT(IN) :: pp !atmos. pressure (Pa)
REAL, DIMENSION(SIZE(PT)) :: ptwb !Wet Bulb Temp. (K)
!local variable
REAL, DIMENSION(SIZE(PT)) :: ztd !Dew Point Temp. (K)
REAL, DIMENSION(SIZE(PT)) :: ztwbinf, ztwbsup, zrv
INTEGER :: jiter, ji
REAL(KIND=JPRB) :: zhook_handle
IF (lhook) CALL dr_hook('MODE_PSYCHRO:TWB_FROM_TPQ_1D',0,zhook_handle)
ztd = td_from_tq(pt, pq) 
!initial guess
ztwbsup = pt
ztwbinf = ztd
ptwb = 0.5 * (ztwbsup + ztwbinf)
DO ji=1,SIZE(pt)
   jiter = 1
   DO WHILE (ztwbsup(ji) - ztwbinf(ji) > 0.001 .OR. jiter .LE. 50)
      zrv(ji) = rv_from_tptwb(pt(ji), pp(ji), ptwb(ji))
      IF (zrv(ji) .GT. pq(ji)/(1 - pq(ji))) THEN
              ztwbsup(ji) = ptwb(ji)
      ELSE
              ztwbinf(ji) = ptwb(ji)
      ENDIF
      ptwb(ji) = 0.5 * (ztwbinf(ji) + ztwbsup(ji))
   END DO
END DO
IF (lhook) CALL dr_hook('MODE_PSYCHRO:TWB_FROM_TPQ_1D',1,zhook_handle)
end function twb_from_tpq_1d
!-------------------------------------------------------------------------------
!
!     ######################################
      FUNCTION enth_fn_t_q(PT,PQ) RESULT(PENTH)
!     ######################################
!
!!
!!    PURPOSE
!!    -------
!       The purpose of this function is to compute the enthalpy function
!       of temperature and humidity content
!      
!
!!**  METHOD
!!    ------
!!  
!!
!!    EXTERNAL
!!    --------
!!      NONE
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------  
!!      
!!    REFERENCE
!!    ---------
!!      
!!
!!
!!    AUTHOR
!!    ------
!!      
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    12/04/11 
!!      A. Alias    01/2013   compi. on Bull : must be 1.0E-5 instead of 1.0D-5
!-------------------------------------------------------------------------------
!
!*       0.    DECLARATIONS
!              ------------
!
IMPLICIT NONE
!
!*       0.1   Declarations of arguments and results
!
!
REAL, INTENT(IN)                :: pt     ! Temperature (K)
REAL, INTENT(IN)                :: pq     ! Humidity content (kg/kg)
REAL                           :: penth  ! Enthalpy (J/kg)
!
!*       0.2   Declarations of local variables
!
REAL        :: zt                          ! Temperature (C)
REAL        :: zrv                         ! Mixing ratio (kg/kg_da)
REAL(KIND=JPRB) :: zhook_handle
!
      IF (lhook) CALL dr_hook('MODE_PSYCHRO:ENTH_FN_T_Q',0,zhook_handle)
! calculate enthalpy
      zt = pt - 273.15
      zrv=max(pq/(1-pq),1.0e-5)
      penth=1.00484d3*zt+zrv*(2.50094d6+1.85895d3*zt)
!
      IF (lhook) CALL dr_hook('MODE_PSYCHRO:ENTH_FN_T_Q',1,zhook_handle)
!-------------------------------------------------------------------------------
!
END FUNCTION enth_fn_t_q
!
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
!
!     ######################################
      FUNCTION q_fn_t_enth(PT,PENTH) RESULT(PQ)
!     ######################################
!
!!
!!    PURPOSE
!!    -------
!       The purpose of this function is to compute the humidity content
!       as a function of temperature and enthalpy
!      
!
!!**  METHOD
!!    ------
!!  
!!
!!    EXTERNAL
!!    --------
!!      NONE
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------  
!!      
!!    REFERENCE
!!    ---------
!!      
!!
!!
!!    AUTHOR
!!    ------
!!      
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    12/04/11 
!-------------------------------------------------------------------------------
!
!*       0.    DECLARATIONS
!              ------------
!
IMPLICIT NONE
!
!*       0.1   Declarations of arguments and results
!
!
REAL, INTENT(IN)                :: pt     ! Temperature (K)
REAL, INTENT(IN)                :: penth  ! Enthalpy (J/kg)
REAL                           :: pq     ! Humidity content (kg/kg)
!
!*       0.2   Declarations of local variables
!
REAL        :: zt                          ! Temperature (C)
REAL        :: zrv                         ! Mixing ratio (kg/kg_da)
REAL(KIND=JPRB) :: zhook_handle
!
      IF (lhook) CALL dr_hook('MODE_PSYCHRO:Q_FN_T_ENTH',0,zhook_handle)
!
      zt = pt - 273.15
!
!    calculate mixing ratio
      zrv=(penth-1.00484d3*zt)/(2.50094d6+1.85895d3*zt)
!
!    validity test
      IF (zrv < 0.0d0) THEN
        zrv=1.d-5
      ENDIF
!
!    calculate humidity content
      pq = zrv/(1+zrv)
!
     IF (lhook) CALL dr_hook('MODE_PSYCHRO:Q_FN_T_ENTH',1,zhook_handle)
!
!-------------------------------------------------------------------------------
!
END FUNCTION q_fn_t_enth

END MODULE mode_psychro