! ! !################################################################## !################################################################## !###### ###### !###### FUNCTION F_ES ###### !###### ###### !###### Developed by ###### !###### Center for Analysis and Prediction of Storms ###### !###### University of Oklahoma ###### !###### ###### !################################################################## !################################################################## ! FUNCTION f_es( p, t ) 7,2 ! !----------------------------------------------------------------------- ! ! PURPOSE: ! ! Calculate the saturation specific humidity using enhanced Teten's ! formula. ! !----------------------------------------------------------------------- ! ! AUTHOR: Yuhe Liu ! 01/08/1998 ! ! MODIFICATION HISTORY: ! !----------------------------------------------------------------------- ! ! INPUT : ! ! p Pressure (Pascal) ! t Temperature (K) ! ! OUTPUT: ! ! f_es Saturation water vapor pressure (Pa) ! !----------------------------------------------------------------------- ! ! !----------------------------------------------------------------------- ! ! Variable Declarations. ! !----------------------------------------------------------------------- ! IMPLICIT NONE REAL :: p ! Pressure (Pascal) REAL :: t ! Temperature (K) REAL :: f_es ! Saturation water vapor pressure (Pa) ! !----------------------------------------------------------------------- ! ! Function f_es and inline directive for Cray PVP ! !----------------------------------------------------------------------- ! REAL :: f_esl, f_esi !fpp$ expand (f_esl) !fpp$ expand (f_esi) !dir$ inline always f_esl, f_esi !*$* inline routine (f_esl, f_esi) ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! IF ( t >= 273.15 ) THEN ! for water f_es = f_esl( p,t ) ELSE ! for ice f_es = f_esi( p,t ) END IF RETURN END FUNCTION f_es ! !----------------------------------------------------------------------- ! ! Calculate the saturation water vapor over liquid water using ! enhanced Teten's formula. ! !----------------------------------------------------------------------- ! FUNCTION f_esl( p, t ) 3 IMPLICIT NONE REAL :: p ! Pressure (Pascal) REAL :: t ! Temperature (K) REAL :: f_esl ! Saturation water vapor pressure over liquid water REAL :: f INCLUDE 'globcst.inc' ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! f = satfwa + satfwb * p f_esl = f * satewa * EXP( satewb*(t-273.15)/(t-satewc) ) RETURN END FUNCTION f_esl ! !----------------------------------------------------------------------- ! ! Calculate the saturation water vapor over ice using enhanced ! Teten's formula. ! !----------------------------------------------------------------------- ! FUNCTION f_esi( p, t ) 3 IMPLICIT NONE REAL :: p ! Pressure (Pascal) REAL :: t ! Temperature (K) REAL :: f_esi ! Saturation water vapor pressure over ice (Pa) REAL :: f INCLUDE 'globcst.inc' ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! f = satfia + satfib * p f_esi = f * sateia * EXP( sateib*(t-273.15)/(t-sateic) ) RETURN END FUNCTION f_esi ! ! !################################################################## !################################################################## !###### ###### !###### FUNCTION F_DESDT ###### !###### ###### !###### Developed by ###### !###### Center for Analysis and Prediction of Storms ###### !###### University of Oklahoma ###### !###### ###### !################################################################## !################################################################## ! FUNCTION f_desdt( t ) 2,2 ! !----------------------------------------------------------------------- ! ! PURPOSE: ! ! Calculate d(es)/dt/es using enhanced Teten's formula. See function ! f_es. ! !----------------------------------------------------------------------- ! ! AUTHOR: Yuhe Liu ! 01/12/1998 ! ! MODIFICATION HISTORY: ! !----------------------------------------------------------------------- ! ! INPUT : ! ! t Temperature (K) ! ! OUTPUT: ! ! f_desdt d(es)/dt/es (1/K) ! !----------------------------------------------------------------------- ! ! !----------------------------------------------------------------------- ! ! Variable Declarations. ! !----------------------------------------------------------------------- ! IMPLICIT NONE REAL :: t ! Temperature (K) REAL :: f_desdt ! d(es)/dt/es (1/K) ! !----------------------------------------------------------------------- ! ! Function f_desdtl and f_desdti and inline directive for Cray PVP ! !----------------------------------------------------------------------- ! REAL :: f_desdtl, f_desdti !fpp$ expand (f_desdtl) !fpp$ expand (f_desdti) !dir$ inline always f_desdtl, f_desdti !*$* inline routine (f_desdtl, f_desdti) ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! IF ( t >= 273.15 ) THEN ! for water f_desdt = f_desdtl(t) ELSE ! for ice f_desdt = f_desdti(t) END IF RETURN END FUNCTION f_desdt ! !----------------------------------------------------------------------- ! ! Calculate d(esl)/dt/esl using enhanced Teten's formula. ! See function f_esl. ! !----------------------------------------------------------------------- ! FUNCTION f_desdtl( t ) 7 IMPLICIT NONE REAL :: t ! Temperature (K) REAL :: f_desdtl ! d(esl)/dt/esl (1/K) INCLUDE 'globcst.inc' ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! f_desdtl = satewb*(273.15-satewc)/(t-satewc)**2 RETURN END FUNCTION f_desdtl ! !----------------------------------------------------------------------- ! ! Calculate d(esi)/dt/esi using enhanced Teten's formula. ! See function f_esi. ! !----------------------------------------------------------------------- ! FUNCTION f_desdti( t ) 4 IMPLICIT NONE REAL :: t ! Temperature (K) REAL :: f_desdti ! d(esi)/dt/esi (1/K) INCLUDE 'globcst.inc' ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! f_desdti = sateib*(273.15-sateic)/(t-sateic)**2 RETURN END FUNCTION f_desdti ! !----------------------------------------------------------------------- ! ! Calculate d(qvs)/dt/qvs using enhanced Teten's formula. ! See function f_es. ! ! Rd ! rddrv = ---- ! Rv ! ! qvs = rddrv * es / (p - ( 1.0 - rddrv ) * es ) ! ! 1 d(qvs) 1 - rddrv 1 d(es) ! ----- -------- = ( 1 + ----------- * qvs ) * ---- ------- ! qvs dT rddrv es dT ! !----------------------------------------------------------------------- ! FUNCTION f_dqvsdt( p, t ),2 IMPLICIT NONE REAL :: p ! Pressure (Pa) REAL :: t ! Temperature (K) REAL :: f_dqvsdt ! d(esi)/dt/esi (1/K) REAL :: qvsat, desdt ! !----------------------------------------------------------------------- ! ! Function f_desdtl and f_desdti and inline directive for Cray PVP ! !----------------------------------------------------------------------- ! REAL :: f_desdt, f_qvsat !fpp$ expand (f_desdt) !fpp$ expand (f_qvsat) !dir$ inline always f_desdt, f_qvsat !*$* inline routine (f_desdt, f_qvsat) INCLUDE 'phycst.inc' ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! qvsat = f_qvsat( p,t ) desdt = f_desdt( t ) f_dqvsdt = ( 1. + qvsat*(1.-rddrv)/rddrv ) * desdt RETURN END FUNCTION f_dqvsdt ! ! !################################################################## !################################################################## !###### ###### !###### FUNCTION F_QVSAT ###### !###### ###### !###### Developed by ###### !###### Center for Analysis and Prediction of Storms ###### !###### University of Oklahoma ###### !###### ###### !################################################################## !################################################################## ! FUNCTION f_qvsat( p, t ) 38,2 ! !----------------------------------------------------------------------- ! ! PURPOSE: ! ! Calculate the saturation specific humidity using enhanced Teten's ! formula. ! !----------------------------------------------------------------------- ! ! AUTHOR: Yuhe Liu ! 01/08/1998 ! ! MODIFICATION HISTORY: ! !----------------------------------------------------------------------- ! ! INPUT : ! ! p Pressure (Pascal) ! t Temperature (K) ! ! OUTPUT: ! ! f_qvsat Saturation water vapor specific humidity (kg/kg). ! !----------------------------------------------------------------------- ! ! !----------------------------------------------------------------------- ! ! Variable Declarations. ! !----------------------------------------------------------------------- ! IMPLICIT NONE REAL :: p ! Pressure (Pascal) REAL :: t ! Temperature (K) REAL :: f_qvsat ! Saturation water vapor specific humidity (kg/kg) ! !----------------------------------------------------------------------- ! ! Include files: ! !----------------------------------------------------------------------- ! INCLUDE 'phycst.inc' ! !----------------------------------------------------------------------- ! ! Function f_es and inline directive for Cray PVP ! !----------------------------------------------------------------------- ! REAL :: f_es !fpp$ expand (f_es) !dir$ inline always f_es !*$* inline routine (f_es) ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! f_qvsat = rddrv * f_es(p,t) / (p-(1.0-rddrv)*f_es(p,t)) RETURN END FUNCTION f_qvsat ! !----------------------------------------------------------------------- ! ! Calculate the saturation specific humidity over liquid water using ! enhanced Teten's formula. ! !----------------------------------------------------------------------- ! FUNCTION f_qvsatl( p, t ) 13,1 IMPLICIT NONE REAL :: p ! Pressure (Pascal) REAL :: t ! Temperature (K) REAL :: f_qvsatl ! Saturation specific humidity over liquid ! water (kg/kg) REAL :: fesl INCLUDE 'phycst.inc' REAL :: f_esl !fpp$ expand (f_esl) !dir$ inline always f_esl !*$* inline routine (f_esl) ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! fesl = f_esl( p, t ) f_qvsatl = rddrv * fesl / (p-(1.0-rddrv)*fesl) RETURN END FUNCTION f_qvsatl ! !----------------------------------------------------------------------- ! ! Calculate the saturation specific humidity over ice using ! enhanced Teten's formula. ! !----------------------------------------------------------------------- ! FUNCTION f_qvsati( p, t ) 5,1 IMPLICIT NONE REAL :: p ! Pressure (Pascal) REAL :: t ! Temperature (K) REAL :: f_qvsati ! Saturation specific humidity over ice (kg/kg) REAL :: fesi INCLUDE 'phycst.inc' REAL :: f_esi !fpp$ expand (f_esi) !dir$ inline always f_esi !*$* inline routine (f_esi) ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! fesi = f_esi( p, t ) f_qvsati = rddrv * fesi / (p-(1.0-rddrv)*fesi) RETURN END FUNCTION f_qvsati ! ! !################################################################## !################################################################## !###### ###### !###### FUNCTION F_MRSAT ###### !###### ###### !###### Developed by ###### !###### Center for Analysis and Prediction of Storms ###### !###### University of Oklahoma ###### !###### ###### !################################################################## !################################################################## ! FUNCTION f_mrsat( p, t ) 9,1 ! !----------------------------------------------------------------------- ! ! PURPOSE: ! ! Calculate the saturation water vapor mixing ratio using enhanced ! Teten's formula. ! !----------------------------------------------------------------------- ! ! AUTHOR: Yuhe Liu ! 01/08/1998 ! ! MODIFICATION HISTORY: ! !----------------------------------------------------------------------- ! ! INPUT : ! ! p Pressure (Pascal) ! t Temperature (K) ! ! OUTPUT: ! ! f_mrsat Saturation water vapor mixing ratio (kg/kg). ! !----------------------------------------------------------------------- ! ! !----------------------------------------------------------------------- ! ! Variable Declarations. ! !----------------------------------------------------------------------- ! IMPLICIT NONE REAL :: p ! Pressure (Pascal) REAL :: t ! Temperature (K) REAL :: f_mrsat ! Saturation water vapor mixing ratio (kg/kg) ! !----------------------------------------------------------------------- ! ! Misc. local variables: ! !----------------------------------------------------------------------- ! REAL :: fes ! !----------------------------------------------------------------------- ! ! Include files: ! !----------------------------------------------------------------------- ! INCLUDE 'phycst.inc' ! !----------------------------------------------------------------------- ! ! Function f_es and inline directive for Cray PVP ! !----------------------------------------------------------------------- ! REAL :: f_es !fpp$ expand (f_es) !dir$ inline always f_es !*$* inline routine (f_es) ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! fes = f_es( p,t ) f_mrsat = rddrv * fes / (p-fes) RETURN END FUNCTION f_mrsat ! !----------------------------------------------------------------------- ! ! Calculate the saturated water vapor mixing ratio over liquid water ! using enhanced Teten's formula. ! !----------------------------------------------------------------------- ! FUNCTION f_mrsatl( p, t ),1 IMPLICIT NONE REAL :: p ! Pressure (Pascal) REAL :: t ! Temperature (K) REAL :: f_mrsatl ! Saturation water vapor mixing ratio over liquid ! water (kg/kg) REAL :: fesl INCLUDE 'phycst.inc' REAL :: f_esl !fpp$ expand (f_esl) !dir$ inline always f_esl !*$* inline routine (f_esl) ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! fesl = f_esl( p,t ) f_mrsatl = rddrv * fesl / (p-fesl) RETURN END FUNCTION f_mrsatl ! !----------------------------------------------------------------------- ! ! Calculate the saturated water vapor mixing ratio over ice ! using enhanced Teten's formula. ! !----------------------------------------------------------------------- ! FUNCTION f_mrsati( p, t ),1 IMPLICIT NONE REAL :: p ! Pressure (Pascal) REAL :: t ! Temperature (K) REAL :: f_mrsati ! Saturation water vapor mixing ratio over ice ! (kg/kg) REAL :: fesi INCLUDE 'phycst.inc' REAL :: f_esi !fpp$ expand (f_esi) !dir$ inline always f_esi !*$* inline routine (f_esi) ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! fesi = f_esi( p,t ) f_mrsati = rddrv * fesi / (p-fesi) RETURN END FUNCTION f_mrsati ! !################################################################## !################################################################## !###### ###### !###### FUNCTION F_TDEW ###### !###### ###### !###### Developed by ###### !###### Center for Analysis and Prediction of Storms ###### !###### University of Oklahoma ###### !###### ###### !################################################################## !################################################################## ! FUNCTION f_tdew( p, t, qv ) 2,3 ! !----------------------------------------------------------------------- ! ! PURPOSE: ! ! Calculate the dew point temperature using reversed enhanced ! Tetan's formula. ! !----------------------------------------------------------------------- ! ! AUTHOR: Yuhe Liu ! 01/09/1998 ! ! MODIFICATION HISTORY: ! !----------------------------------------------------------------------- ! ! INPUT : ! ! p Pressure (Pascal) ! t Temperature (K) ! qv Specific humidity (kg/kg) ! ! OUTPUT: ! ! f_tdew Dew point temperature (K) ! !----------------------------------------------------------------------- ! ! !----------------------------------------------------------------------- ! ! Variable Declarations. ! !----------------------------------------------------------------------- ! IMPLICIT NONE REAL :: p ! Pressure (Pascal) REAL :: t ! Temperature (K) REAL :: qv ! Specific humidity (kg/kg) REAL :: f_tdew ! Dew point temperature (K) ! !----------------------------------------------------------------------- ! ! Include file ! !----------------------------------------------------------------------- ! INCLUDE 'globcst.inc' ! !----------------------------------------------------------------------- ! ! Function f_tdewl and f_tdewi and inline directive for Cray PVP ! !----------------------------------------------------------------------- ! REAL :: f_tdewl, f_tdewi !fpp$ expand (f_tdewl) !fpp$ expand (f_tdewi) !dir$ inline always f_tdewl, f_tdewi !*$* inline routine (f_tdewl, f_tdewi) ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! IF ( t >= 273.15 ) THEN ! for water f_tdew = f_tdewl( p,qv ) IF ( f_tdew < 273.15 ) THEN ! using ice formula f_tdew = f_tdewi( p,qv ) END IF ELSE ! for ice f_tdew = f_tdewi( p,qv ) END IF RETURN END FUNCTION f_tdew ! !----------------------------------------------------------------------- ! ! Calculate the dew point temperature over liquid water using the ! reversed enhanced Teten's formula. ! !----------------------------------------------------------------------- ! FUNCTION f_tdewl( p, qv ) 1 IMPLICIT NONE REAL :: p ! Pressure (Pascal) REAL :: qv ! Specific humidity (kg/kg) REAL :: f_tdewl ! Dew point temperature over liquid water (K) REAL :: qvs, esl, f, ln INCLUDE 'globcst.inc' INCLUDE 'phycst.inc' ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! qvs = MAX( qv, 1.e-8 ) f = satfwa + satfwb * p esl = p*qvs/(rddrv+(1-rddrv)*qvs) ln = ALOG( esl/(f*satewa) ) f_tdewl = ( ln*satewc - 273.15*satewb ) / ( ln - satewb ) RETURN END FUNCTION f_tdewl ! !----------------------------------------------------------------------- ! ! Calculate the dew point temperature over ice using the ! reversed enhanced Teten's formula. ! !----------------------------------------------------------------------- ! FUNCTION f_tdewi( p, qv ) 2 IMPLICIT NONE REAL :: p ! Pressure (Pascal) REAL :: qv ! Specific humidity (kg/kg) REAL :: f_tdewi ! Dew point temperature over ice (K) REAL :: qvs, esi, f, ln INCLUDE 'globcst.inc' INCLUDE 'phycst.inc' ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! qvs = MAX( qv, 1.e-8 ) f = satfia + satfib * p esi = p*qvs/(rddrv+(1-rddrv)*qvs) ln = ALOG( esi/(f*sateia) ) f_tdewi = ( ln*sateic - 273.15*sateib ) / ( ln - sateib ) RETURN END FUNCTION f_tdewi ! !################################################################## !################################################################## !###### ###### !###### SUBROUTINE F_PT2PTE ###### !###### ###### !###### Developed by ###### !###### Center for Analysis and Prediction of Storms ###### !###### University of Oklahoma ###### !###### ###### !################################################################## !################################################################## ! FUNCTION f_pt2pte( p, pt, qv ) 1,1 ! !----------------------------------------------------------------------- ! ! PURPOSE: ! ! Calculate the equivalent potential temperature of an air ! parcel that is either saturated or unsaturated. ! !----------------------------------------------------------------------- ! ! AUTHOR: Yuhe Liu ! 01/14/1998 Re-created from subroutine EQUIPT ! ! MODIFICATION HISTORY: ! !----------------------------------------------------------------------- ! ! INPUT: ! ! p Base state pressure (Pascals) ! pt Potential temperature (degrees Kelvin) ! qv Water vapor specific humidity (kg/kg) ! ! OUTPUT: ! ! f_pt2pte Equivalent potential temperature (K) ! !----------------------------------------------------------------------- ! ! Variable Declarations. ! !----------------------------------------------------------------------- ! IMPLICIT NONE REAL :: p ! Pressure (Pascals) REAL :: pt ! Potential temperature (degrees Kelvin) REAL :: qv ! Water vapor specific humidity (kg/kg) REAL :: f_pt2pte ! Equivalent potential temperature (K) ! !----------------------------------------------------------------------- ! ! Misc. local variables ! !----------------------------------------------------------------------- ! REAL :: t ! Temperature REAL :: tdew ! Dew point temperature ! !----------------------------------------------------------------------- ! ! Include files: ! !----------------------------------------------------------------------- ! INCLUDE 'globcst.inc' INCLUDE 'phycst.inc' ! !----------------------------------------------------------------------- ! ! Function f_tdew and inline directive for Cray PVP ! !----------------------------------------------------------------------- ! REAL :: f_tdew !fpp$ expand (f_tdew) !dir$ inline always f_tdew !*$* inline routine (f_tdew) ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! t = pt * (p/p0)**rddcp tdew = f_tdew( p,t,qv ) ! Dew point temperture. f_pt2pte = pt * EXP( lathv*qv/(cp*tdew) ) RETURN END FUNCTION f_pt2pte ! !################################################################## !################################################################## !###### ###### !###### SUBROUTINE GETQVS ###### !###### ###### !###### Developed by ###### !###### Center for Analysis and Prediction of Storms ###### !###### University of Oklahoma ###### !###### ###### !################################################################## !################################################################## ! SUBROUTINE getqvs( nx,ny,nz, ibgn,iend,jbgn,jend,kbgn,kend, & 14,1 p, t, qvsat ) ! !----------------------------------------------------------------------- ! ! PURPOSE: ! ! Calculate the saturation specific humidity using Teten's formula. ! !----------------------------------------------------------------------- ! ! AUTHOR: Ming Xue ! 3/17/1991. ! ! MODIFICATION HISTORY: ! ! 5/02/92 (M. Xue) ! Added full documentation. ! ! 5/03/92 (M. Xue) ! Further documentation. ! ! 2/10/93 (K. Droegemeier) ! Cleaned up documentation. ! ! 9/10/94 (Weygandt & Y. Lu) ! Cleaned up documentation. ! ! 12/11/1997 (Yuhe Liu) ! Rewrote the subroutine calling function F_QVSAT to calculate the ! saturated specific humidity for array. ! !----------------------------------------------------------------------- ! ! INPUT : ! ! nx Array dimension in x-direction ! ny Array dimension in y-direction ! nz Array dimension in z-direction ! ibgn Starting index in x-direction ! jbgn Starting index in y-direction ! kbgn Starting index in z-direction ! iend Last index in x-direction ! jend Last index in y-direction ! kend Last index in z-direction ! ! p Pressure (Pascal) ! t Temperature (K) ! ! OUTPUT: ! ! qvsat Saturation water vapor specific humidity (kg/kg). ! !----------------------------------------------------------------------- ! ! !----------------------------------------------------------------------- ! ! Variable Declarations. ! !----------------------------------------------------------------------- ! IMPLICIT NONE INTEGER :: nx,ny,nz ! Input array dimension INTEGER :: ibgn ! Starting index in x-direction INTEGER :: jbgn ! Starting index in y-direction INTEGER :: kbgn ! Starting index in z-direction INTEGER :: iend ! Last index in x-direction INTEGER :: jend ! Last index in y-direction INTEGER :: kend ! Last index in z-direction REAL :: p (nx,ny,nz) ! Pressure (Pascal) REAL :: t (nx,ny,nz) ! Temperature (K) REAL :: qvsat(nx,ny,nz) ! Saturation water vapor specific ! humidity (kg/kg) ! !----------------------------------------------------------------------- ! ! Include file ! !----------------------------------------------------------------------- ! INCLUDE 'globcst.inc' ! !----------------------------------------------------------------------- ! ! Misc. local variables: ! !----------------------------------------------------------------------- ! INTEGER :: i,j,k ! !----------------------------------------------------------------------- ! ! Function f_qvsat and inline directives for Cray PVP ! !----------------------------------------------------------------------- ! REAL :: f_qvsat !fpp$ expand (f_qvsat) !dir$ inline always f_qvsat !*$* inline routine (f_qvsat) ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! DO k=kbgn,kend DO j=jbgn,jend DO i=ibgn,iend qvsat(i,j,k) = f_qvsat( p(i,j,k), t(i,j,k) ) END DO END DO END DO RETURN END SUBROUTINE getqvs ! ! !################################################################## !################################################################## !###### ###### !###### SUBROUTINE GETMRS ###### !###### ###### !###### Developed by ###### !###### Center for Analysis and Prediction of Storms ###### !###### University of Oklahoma ###### !###### ###### !################################################################## !################################################################## ! SUBROUTINE getmrs( nx,ny,nz, ibgn,iend,jbgn,jend,kbgn,kend, &,1 p, t, mrsat ) ! !----------------------------------------------------------------------- ! ! PURPOSE: ! ! Calculate the saturation mixing ratio using Teten's formula. ! !----------------------------------------------------------------------- ! ! AUTHOR: Yuhe Liu ! 01/09/1998 ! ! MODIFICATION HISTORY: ! !----------------------------------------------------------------------- ! ! INPUT : ! ! nx Array dimension in x-direction ! ny Array dimension in y-direction ! nz Array dimension in z-direction ! ibgn Starting index in x-direction ! jbgn Starting index in y-direction ! kbgn Starting index in z-direction ! iend Last index in x-direction ! jend Last index in y-direction ! kend Last index in z-direction ! ! p Pressure (Pascal) ! t Temperature (K) ! ! OUTPUT: ! ! mrsat Saturation water vapor Mixing ratio (kg/kg). ! !----------------------------------------------------------------------- ! ! !----------------------------------------------------------------------- ! ! Variable Declarations. ! !----------------------------------------------------------------------- ! IMPLICIT NONE INTEGER :: nx,ny,nz ! Input array dimension INTEGER :: ibgn ! Starting index in x-direction INTEGER :: jbgn ! Starting index in y-direction INTEGER :: kbgn ! Starting index in z-direction INTEGER :: iend ! Last index in x-direction INTEGER :: jend ! Last index in y-direction INTEGER :: kend ! Last index in z-direction REAL :: p (nx,ny,nz) ! Pressure (Pascal) REAL :: t (nx,ny,nz) ! Temperature (K) REAL :: mrsat(nx,ny,nz) ! Saturated water vapor mixing ratio (kg/kg) ! !----------------------------------------------------------------------- ! ! Include file ! !----------------------------------------------------------------------- ! INCLUDE 'globcst.inc' ! !----------------------------------------------------------------------- ! ! Misc. local variables: ! !----------------------------------------------------------------------- ! INTEGER :: i,j,k ! !----------------------------------------------------------------------- ! ! Function f_mrsat and inline directive for Cray PVP ! !----------------------------------------------------------------------- ! REAL :: f_mrsat !fpp$ expand (f_mrsat) !dir$ inline always f_mrsat !*$* inline routine (f_mrsat) ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! DO k=kbgn,kend DO j=jbgn,jend DO i=ibgn,iend mrsat(i,j,k) = f_mrsat( p(i,j,k), t(i,j,k) ) END DO END DO END DO RETURN END SUBROUTINE getmrs ! ! !################################################################## !################################################################## !###### ###### !###### SUBROUTINE GETDEW ###### !###### ###### !###### Developed by ###### !###### Center for Analysis and Prediction of Storms ###### !###### University of Oklahoma ###### !###### ###### !################################################################## !################################################################## ! SUBROUTINE getdew( nx,ny,nz, ibgn,iend,jbgn,jend,kbgn,kend, & 2,1 p, t, qv, tdew ) ! !----------------------------------------------------------------------- ! ! PURPOSE: ! ! Calculate the dew point temperature (K) using enhanced Teten's ! formula. ! !----------------------------------------------------------------------- ! ! AUTHOR: Yuhe Liu ! 01/09/1998 ! ! MODIFICATION HISTORY: ! !----------------------------------------------------------------------- ! ! INPUT : ! ! nx Array dimension in x-direction ! ny Array dimension in y-direction ! nz Array dimension in z-direction ! ibgn Starting index in x-direction ! jbgn Starting index in y-direction ! kbgn Starting index in z-direction ! iend Last index in x-direction ! jend Last index in y-direction ! kend Last index in z-direction ! ! p Pressure (Pascal) ! t Temperature (K) ! qv Specific humidity (kg/kg) ! ! OUTPUT: ! ! tdew Dew point temperature (K) ! !----------------------------------------------------------------------- ! ! !----------------------------------------------------------------------- ! ! Variable Declarations. ! !----------------------------------------------------------------------- ! IMPLICIT NONE INTEGER :: nx,ny,nz ! Input array dimension INTEGER :: ibgn ! Starting index in x-direction INTEGER :: jbgn ! Starting index in y-direction INTEGER :: kbgn ! Starting index in z-direction INTEGER :: iend ! Last index in x-direction INTEGER :: jend ! Last index in y-direction INTEGER :: kend ! Last index in z-direction REAL :: p (nx,ny,nz) ! Pressure (Pascal) REAL :: t (nx,ny,nz) ! Temperature (K) REAL :: qv (nx,ny,nz) ! Specific humidity (kg/kg) REAL :: tdew(nx,ny,nz) ! Dew point temperature (K) ! !----------------------------------------------------------------------- ! ! Misc. local variables: ! !----------------------------------------------------------------------- ! INTEGER :: i,j,k ! !----------------------------------------------------------------------- ! ! Function f_mrsat and inline directive for Cray PVP ! !----------------------------------------------------------------------- ! REAL :: f_tdew !fpp$ expand (f_tdew) !dir$ inline always f_tdew !*$* inline routine (f_tdew) ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! DO k=kbgn,kend DO j=jbgn,jend DO i=ibgn,iend tdew(i,j,k) = f_tdew( p(i,j,k), t(i,j,k), qv(i,j,k) ) END DO END DO END DO RETURN END SUBROUTINE getdew ! ! !################################################################## !################################################################## !###### ###### !###### SUBROUTINE PT2PTE ###### !###### ###### !###### Developed by ###### !###### Center for Analysis and Prediction of Storms ###### !###### University of Oklahoma ###### !###### ###### !################################################################## !################################################################## ! SUBROUTINE pt2pte( nx,ny,nz, ibgn,iend,jbgn,jend,kbgn,kend, &,1 p, pt, qv, pte ) ! !----------------------------------------------------------------------- ! ! PURPOSE: ! ! Calculate the dew point temperature (K) using enhanced Teten's ! formula. ! !----------------------------------------------------------------------- ! ! AUTHOR: Yuhe Liu ! 01/09/1998 ! ! MODIFICATION HISTORY: ! !----------------------------------------------------------------------- ! ! INPUT : ! ! nx Array dimension in x-direction ! ny Array dimension in y-direction ! nz Array dimension in z-direction ! ibgn Starting index in x-direction ! jbgn Starting index in y-direction ! kbgn Starting index in z-direction ! iend Last index in x-direction ! jend Last index in y-direction ! kend Last index in z-direction ! ! p Pressure (Pascal) ! pt Temperature (K) ! qv Specific humidity (kg/kg) ! ! OUTPUT: ! ! pte Equivalent potential temperature (K) ! !----------------------------------------------------------------------- ! ! !----------------------------------------------------------------------- ! ! Variable Declarations. ! !----------------------------------------------------------------------- ! IMPLICIT NONE INTEGER :: nx,ny,nz ! Input array dimension INTEGER :: ibgn ! Starting index in x-direction INTEGER :: jbgn ! Starting index in y-direction INTEGER :: kbgn ! Starting index in z-direction INTEGER :: iend ! Last index in x-direction INTEGER :: jend ! Last index in y-direction INTEGER :: kend ! Last index in z-direction REAL :: p (nx,ny,nz) ! Pressure (Pascal) REAL :: pt (nx,ny,nz) ! Temperature (K) REAL :: qv (nx,ny,nz) ! Specific humidity (kg/kg) REAL :: pte(nx,ny,nz) ! Equivalent potential temperature (K) ! !----------------------------------------------------------------------- ! ! Misc. local variables: ! !----------------------------------------------------------------------- ! INTEGER :: i,j,k ! !----------------------------------------------------------------------- ! ! Function f_pt2pte and inline directive for Cray PVP ! !----------------------------------------------------------------------- ! REAL :: f_pt2pte !fpp$ expand (f_pt2pte) !dir$ inline always f_pt2pte !*$* inline routine (f_pt2pte) ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! DO k=kbgn,kend DO j=jbgn,jend DO i=ibgn,iend pte(i,j,k) = f_pt2pte( p(i,j,k), pt(i,j,k), qv(i,j,k) ) END DO END DO END DO RETURN END SUBROUTINE pt2pte ! !################################################################## !################################################################## !###### ###### !###### FUNCTION F_PCCL ###### !###### ###### !###### Developed by ###### !###### Center for Analysis and Prediction of Storms ###### !###### University of Oklahoma ###### !###### ###### !################################################################## !################################################################## ! FUNCTION f_pccl(pm,p,t,td,mrbar,n),7 ! !----------------------------------------------------------------------- ! ! PURPOSE: ! ! this function returns the pressure at the convective condensation ! level given the appropriate sounding data. ! ! the algorithm is decribed on p.17 of stipanuk, g.s.,1973: ! "algorithms for generating a skew-t log p diagram and computing ! selected meteorological quantities," atmospheric sciences labora- ! tory, u.s. army electronics command, white sands missile range, new ! mexico 88002. ! !----------------------------------------------------------------------- ! ! AUTHOR: ! g.s. stipanuk 1973 original version. ! reference stipanuk paper entitled: ! "algorithms for generating a skew-t, log p ! diagram and computing selected meteorological ! quantities." ! atmospheric sciences laboratory ! u.s. army electronics command ! white sands missile range, new mexico 88002 ! 33 pages ! baker, schlatter 17-may-1982 ! ! ! MODIFICATIONS: ! ! 01/30/1998 (Yuhe Liu) ! Modified to use Kelvin and Pascal as units of temperature and ! pressure ! !----------------------------------------------------------------------- ! ! INPUT: ! ! p pressure (Pa). note that p(i).gt.p(i+1). ! t temperature (K) ! td dew point (K) ! n number of levels in the sounding and the dimension of ! p, t and td ! pm pressure (Pa) at upper boundary of the layer for ! computing the mean mixing ratio. p(1) is the lower ! boundary. ! ! OUTPUT: ! ! mrbar mean mixing ratio (kg/kg) in the layer bounded by ! pressures p(1) at the bottom and pm at the top ! f_pccl pressure (Pa) at the convective condensation level ! !----------------------------------------------------------------------- ! ! Variable Declarations. ! !----------------------------------------------------------------------- ! IMPLICIT NONE INTEGER :: n REAL :: t(n),p(n),td(n) REAL :: pm REAL :: mrbar REAL :: f_pccl ! !----------------------------------------------------------------------- ! ! Misc. local variables ! !----------------------------------------------------------------------- ! INTEGER :: i,j,k,l REAL :: pc,tq,del,x,a REAL :: mrsat1,mrsat2 ! !----------------------------------------------------------------------- ! ! Function f_mrsat and inline directive for Cray PVP ! !----------------------------------------------------------------------- ! REAL :: f_mrsat REAL :: f_tmr !fpp$ expand (f_mrsat) !fpp$ expand (f_tmr) !dir$ inline always f_mrsat, f_tmr !*$* inline routine (f_mrsat, f_tmr) ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! IF (pm /= p(1)) GO TO 5 mrbar = f_mrsat( p(1),td(1) ) pc= pm IF (ABS( GO TO 25 5 CONTINUE mrbar = 0. k = 0 10 CONTINUE k = k+1 IF (p(k) > pm) GO TO 10 k = k-1 j = k-1 IF (j < 1) GO TO 20 ! compute the average mixing ratio....alog = natural log DO i= 1,j l = i+1 mrsat1 = f_mrsat( p(i),td(i) ) mrsat2 = f_mrsat( p(l),td(l) ) mrbar = (mrsat1+mrsat2)*ALOG( END DO 20 CONTINUE l= k+1 ! estimate the dew point at pressure pm. mrsat1 = f_mrsat( p(k),td(k) ) tq= td(k)+(td(l)-td(k))*(ALOG(pm/p(k)))/(ALOG( mrsat2 = f_mrsat( pm,tq ) mrbar= mrbar+(mrsat1+mrsat2 )*ALOG( mrbar= mrbar/(2.*ALOG( ! find level at which the mixing ratio line mrbar crosses the ! environmental temperature profile. 25 CONTINUE DO j= 1,n i = n-j+1 IF (p(i) < 300.) CYCLE ! f_tmr = temperature (celsius) at pressure p (mb) along a mixing ! ratio line given by mrbar (g/kg) x = f_tmr(mrbar,p(i))-t(i) IF (x <= 0.) GO TO 35 END DO f_pccl = 0.0 RETURN ! set up bisection routine 35 l = i i = i+1 del = p(l)-p(i) pc = p(i)+.5*del a = (t(i)-t(l))/ALOG( DO j = 1,10 del = del/2. x = f_tmr(mrbar,pc)-t(l)-a*(ALOG( ! the sign function replaces the sign of the first argument ! with that of the second. pc = pc+SIGN(del,x) END DO 45 f_pccl = pc RETURN END FUNCTION f_pccl ! ! !################################################################## !################################################################## !###### ###### !###### FUNCTION F_CT ###### !###### ###### !###### Developed by ###### !###### Center for Analysis and Prediction of Storms ###### !###### University of Oklahoma ###### !###### ###### !################################################################## !################################################################## ! FUNCTION f_ct(mrbar,pc,ps),1 ! !----------------------------------------------------------------------- ! ! PURPOSE: ! ! This function returns the convective temperature ct (K) ! given the mean mixing ratio mrbar (kg/kg) in the surface layer, ! the pressure pc (Pa) at the convective condensation level (ccl) ! and the surface pressure ps (Pa). ! !----------------------------------------------------------------------- ! ! AUTHOR: G.S. Stipanuk 1973 original version. ! ! reference stipanuk paper entitled: ! "algorithms for generating a skew-t, log p ! diagram and computing selected meteorological ! quantities." ! atmospheric sciences laboratory ! u.s. army electronics command ! white sands missile range, new mexico 88002 ! 33 pages ! baker, schlatter 17-may-1982 ! ! MODIFICATION HISTORY: ! ! 01/30/1998 (Yuhe Liu) ! Modified for ARPS to use the units of Kelvin, Pascal and kg/kg for ! temperature, pressure and mixing ratio respectively ! !----------------------------------------------------------------------- ! ! INPUT : ! ! mrbar Mean mixing ratio (kg/kg) ! pc Pressure at CCL (Pa) ! ps Pressure at surface (Pa) ! ! OUTPUT: ! ! f_ct Convective temperature (K) ! !----------------------------------------------------------------------- ! ! Variable Declarations. ! !----------------------------------------------------------------------- ! IMPLICIT NONE REAL :: mrbar REAL :: pc REAL :: ps REAL :: f_ct REAL :: tc ! !----------------------------------------------------------------------- ! ! Include file ! !----------------------------------------------------------------------- ! INCLUDE 'phycst.inc' ! !----------------------------------------------------------------------- ! ! Function f_tmr and inline directive for Cray PVP ! !----------------------------------------------------------------------- ! REAL :: f_tmr !fpp$ expand (f_tmr) !dir$ inline always f_tmr !*$* inline routine (f_tmr) ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! compute the temperature (K) at the ccl. tc= f_tmr(mrbar,pc) ! !----------------------------------------------------------------------- ! ! compute the potential temperature (K) at ccl, i.e., the dry ! adiabat through the ccl. ! ! ptc = tc*((100000./pc)**rddcp) ! ! and then compute the surface temperature on the same dry adiabat. ! ! f_ct = ptc*((ps/100000.)**rddcp) ! !----------------------------------------------------------------------- ! f_ct = tc*((ps/pc)**rddcp) RETURN END FUNCTION f_ct ! ! !################################################################## !################################################################## !###### ###### !###### FUNCTION F_TMR ###### !###### ###### !###### Developed by ###### !###### Center for Analysis and Prediction of Storms ###### !###### University of Oklahoma ###### !###### ###### !################################################################## !################################################################## ! FUNCTION f_tmr(p,mr) 3 ! !----------------------------------------------------------------------- ! ! PURPOSE: ! ! This function returns the temperature (Kelvin) on a mixing ! ratio line mr (kg/kg) at pressure p (Pa). The formula is given in ! table 1 on page 7 of Stipanuk (1973). ! !----------------------------------------------------------------------- ! ! AUTHOR: G.S. Stipanuk, 1973 original version. ! ! reference stipanuk paper entitled: ! "algorithms for generating a skew-t, log p ! diagram and computing selected meteorological ! quantities." ! atmospheric sciences laboratory ! u.s. army electronics command ! white sands missile range, new mexico 88002 ! 33 pages ! baker, schlatter 17-may-1982 ! ! MODIFICATIONS: ! ! 01/30/1998 (Yuhe Liu) ! Modified to use the units of Kelvin, Pascal and kg/kg for ! temperature, pressure and mixing ratio respectively. ! !----------------------------------------------------------------------- ! ! INPUT: ! ! p Pressure (Pa) ! mr Mixing ratio (kg/kg) ! ! OUTPUT: ! ! f_tmr Temperature (K) ! !----------------------------------------------------------------------- ! ! Variable Declarations. ! !----------------------------------------------------------------------- ! IMPLICIT NONE REAL :: p ! Pressure (Pa) REAL :: mr ! Mixing ratio (kg/kg) REAL :: f_tmr ! Temperature (K) ! !----------------------------------------------------------------------- ! ! Local variables ! !----------------------------------------------------------------------- ! REAL :: x, pmb, mrg REAL :: c1,c2,c3,c4,c5,c6 DATA c1/.0498646455/ DATA c2/2.4082965/ DATA c3/7.07475/ DATA c4/38.9114/ DATA c5/.0915/ DATA c6/1.2035/ ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! pmb = p/100. mrg = mr*1000. x= ALOG10(mrg*pmb/(622.+mrg)) f_tmr= 10.**(c1*x+c2)-c3+c4*((10.**(c5*x)-c6)**2.) RETURN END FUNCTION f_tmr