!>\file funcphys.f90 !! This file includes API for basic thermodynamic physics. !>\defgroup func_phys GFS Physics Function Module !! This module provides API for computing basic thermodynamic physics !! functions. !> This module provides an Application Program Interface (API) for computing !! basic thermodynamic physics functions, in particular: !! -# saturation vapor pressure as a function of temperature; !! -# dewpoint temperature as a function of vapor pressure; !! -# equivalent potential temperature as a function of temperature and !! scaled pressure to the kappa power; !! -# temperature and specific humidity along a moist adiabat as functions !! of equivalent potential temperature and scaled pressure to the kappa power; !! -# scaled pressure to the kappa power as a function of pressure, and !! -# temperature at the lifting condensation level as a function of temperature !! and dewpoint depression. !! !! The entry points required to set up lookup tables start with a "g". !! All the other entry points are functions starting with an "f" or !! are subroutines starting with an "s". These other functions and subroutines !! are elemental; that is, they return a scalar if they are passed only scalars, !! but they return an array if they are passed an array. These other functions !! and subroutines can in inlined, too. module funcphys !$$$ Module Documentation Block ! ! Module: funcphys API for basic thermodynamic physics ! Author: Iredell Org: W/NX23 Date: 1999-03-01 ! ! Abstract: This module provides an Application Program Interface ! for computing basic thermodynamic physics functions, in particular ! (1) saturation vapor pressure as a function of temperature, ! (2) dewpoint temperature as a function of vapor pressure, ! (3) equivalent potential temperature as a function of temperature ! and scaled pressure to the kappa power, ! (4) temperature and specific humidity along a moist adiabat ! as functions of equivalent potential temperature and ! scaled pressure to the kappa power, ! (5) scaled pressure to the kappa power as a function of pressure, and ! (6) temperature at the lifting condensation level as a function ! of temperature and dewpoint depression. ! The entry points required to set up lookup tables start with a "g". ! All the other entry points are functions starting with an "f" or ! are subroutines starting with an "s". These other functions and ! subroutines are elemental; that is, they return a scalar if they ! are passed only scalars, but they return an array if they are passed ! an array. These other functions and subroutines can be inlined, too. ! ! Program History Log: ! 1999-03-01 Mark Iredell ! 1999-10-15 Mark Iredell SI unit for pressure (Pascals) ! 2001-02-26 Mark Iredell Ice phase changes of Hong and Moorthi ! ! Public Variables: ! krealfp Integer parameter kind or length of reals (=kind_phys) ! ! Public Subprograms: ! gpvsl Compute saturation vapor pressure over liquid table ! ! fpvsl Elementally compute saturation vapor pressure over liquid ! function result Real(krealfp) saturation vapor pressure in Pascals ! t Real(krealfp) temperature in Kelvin ! ! fpvslq Elementally compute saturation vapor pressure over liquid ! function result Real(krealfp) saturation vapor pressure in Pascals ! t Real(krealfp) temperature in Kelvin ! ! fpvslx Elementally compute saturation vapor pressure over liquid ! function result Real(krealfp) saturation vapor pressure in Pascals ! t Real(krealfp) temperature in Kelvin ! ! gpvsi Compute saturation vapor pressure over ice table ! ! fpvsi Elementally compute saturation vapor pressure over ice ! function result Real(krealfp) saturation vapor pressure in Pascals ! t Real(krealfp) temperature in Kelvin ! ! fpvsiq Elementally compute saturation vapor pressure over ice ! function result Real(krealfp) saturation vapor pressure in Pascals ! t Real(krealfp) temperature in Kelvin ! ! fpvsix Elementally compute saturation vapor pressure over ice ! function result Real(krealfp) saturation vapor pressure in Pascals ! t Real(krealfp) temperature in Kelvin ! ! gpvs Compute saturation vapor pressure table ! ! fpvs Elementally compute saturation vapor pressure ! function result Real(krealfp) saturation vapor pressure in Pascals ! t Real(krealfp) temperature in Kelvin ! ! fpvsq Elementally compute saturation vapor pressure ! function result Real(krealfp) saturation vapor pressure in Pascals ! t Real(krealfp) temperature in Kelvin ! ! fpvsx Elementally compute saturation vapor pressure ! function result Real(krealfp) saturation vapor pressure in Pascals ! t Real(krealfp) temperature in Kelvin ! ! gtdpl Compute dewpoint temperature over liquid table ! ! ftdpl Elementally compute dewpoint temperature over liquid ! function result Real(krealfp) dewpoint temperature in Kelvin ! pv Real(krealfp) vapor pressure in Pascals ! ! ftdplq Elementally compute dewpoint temperature over liquid ! function result Real(krealfp) dewpoint temperature in Kelvin ! pv Real(krealfp) vapor pressure in Pascals ! ! ftdplx Elementally compute dewpoint temperature over liquid ! function result Real(krealfp) dewpoint temperature in Kelvin ! pv Real(krealfp) vapor pressure in Pascals ! ! ftdplxg Elementally compute dewpoint temperature over liquid ! function result Real(krealfp) dewpoint temperature in Kelvin ! t Real(krealfp) guess dewpoint temperature in Kelvin ! pv Real(krealfp) vapor pressure in Pascals ! ! gtdpi Compute dewpoint temperature table over ice ! ! ftdpi Elementally compute dewpoint temperature over ice ! function result Real(krealfp) dewpoint temperature in Kelvin ! pv Real(krealfp) vapor pressure in Pascals ! ! ftdpiq Elementally compute dewpoint temperature over ice ! function result Real(krealfp) dewpoint temperature in Kelvin ! pv Real(krealfp) vapor pressure in Pascals ! ! ftdpix Elementally compute dewpoint temperature over ice ! function result Real(krealfp) dewpoint temperature in Kelvin ! pv Real(krealfp) vapor pressure in Pascals ! ! ftdpixg Elementally compute dewpoint temperature over ice ! function result Real(krealfp) dewpoint temperature in Kelvin ! t Real(krealfp) guess dewpoint temperature in Kelvin ! pv Real(krealfp) vapor pressure in Pascals ! ! gtdp Compute dewpoint temperature table ! ! ftdp Elementally compute dewpoint temperature ! function result Real(krealfp) dewpoint temperature in Kelvin ! pv Real(krealfp) vapor pressure in Pascals ! ! ftdpq Elementally compute dewpoint temperature ! function result Real(krealfp) dewpoint temperature in Kelvin ! pv Real(krealfp) vapor pressure in Pascals ! ! ftdpx Elementally compute dewpoint temperature ! function result Real(krealfp) dewpoint temperature in Kelvin ! pv Real(krealfp) vapor pressure in Pascals ! ! ftdpxg Elementally compute dewpoint temperature ! function result Real(krealfp) dewpoint temperature in Kelvin ! t Real(krealfp) guess dewpoint temperature in Kelvin ! pv Real(krealfp) vapor pressure in Pascals ! ! gthe Compute equivalent potential temperature table ! ! fthe Elementally compute equivalent potential temperature ! function result Real(krealfp) equivalent potential temperature in Kelvin ! t Real(krealfp) LCL temperature in Kelvin ! pk Real(krealfp) LCL pressure over 1e5 Pa to the kappa power ! ! ftheq Elementally compute equivalent potential temperature ! function result Real(krealfp) equivalent potential temperature in Kelvin ! t Real(krealfp) LCL temperature in Kelvin ! pk Real(krealfp) LCL pressure over 1e5 Pa to the kappa power ! ! fthex Elementally compute equivalent potential temperature ! function result Real(krealfp) equivalent potential temperature in Kelvin ! t Real(krealfp) LCL temperature in Kelvin ! pk Real(krealfp) LCL pressure over 1e5 Pa to the kappa power ! ! gtma Compute moist adiabat tables ! ! stma Elementally compute moist adiabat temperature and moisture ! the Real(krealfp) equivalent potential temperature in Kelvin ! pk Real(krealfp) pressure over 1e5 Pa to the kappa power ! tma Real(krealfp) parcel temperature in Kelvin ! qma Real(krealfp) parcel specific humidity in kg/kg ! ! stmaq Elementally compute moist adiabat temperature and moisture ! the Real(krealfp) equivalent potential temperature in Kelvin ! pk Real(krealfp) pressure over 1e5 Pa to the kappa power ! tma Real(krealfp) parcel temperature in Kelvin ! qma Real(krealfp) parcel specific humidity in kg/kg ! ! stmax Elementally compute moist adiabat temperature and moisture ! the Real(krealfp) equivalent potential temperature in Kelvin ! pk Real(krealfp) pressure over 1e5 Pa to the kappa power ! tma Real(krealfp) parcel temperature in Kelvin ! qma Real(krealfp) parcel specific humidity in kg/kg ! ! stmaxg Elementally compute moist adiabat temperature and moisture ! tg Real(krealfp) guess parcel temperature in Kelvin ! the Real(krealfp) equivalent potential temperature in Kelvin ! pk Real(krealfp) pressure over 1e5 Pa to the kappa power ! tma Real(krealfp) parcel temperature in Kelvin ! qma Real(krealfp) parcel specific humidity in kg/kg ! ! gpkap Compute pressure to the kappa table ! ! fpkap Elementally raise pressure to the kappa power. ! function result Real(krealfp) p over 1e5 Pa to the kappa power ! p Real(krealfp) pressure in Pascals ! ! fpkapq Elementally raise pressure to the kappa power. ! function result Real(krealfp) p over 1e5 Pa to the kappa power ! p Real(krealfp) pressure in Pascals ! ! fpkapo Elementally raise pressure to the kappa power. ! function result Real(krealfp) p over 1e5 Pa to the kappa power ! p Real(krealfp) surface pressure in Pascals ! ! fpkapx Elementally raise pressure to the kappa power. ! function result Real(krealfp) p over 1e5 Pa to the kappa power ! p Real(krealfp) pressure in Pascals ! ! grkap Compute pressure to the 1/kappa table ! ! frkap Elementally raise pressure to the 1/kappa power. ! function result Real(krealfp) pressure in Pascals ! pkap Real(krealfp) p over 1e5 Pa to the 1/kappa power ! ! frkapq Elementally raise pressure to the kappa power. ! function result Real(krealfp) pressure in Pascals ! pkap Real(krealfp) p over 1e5 Pa to the kappa power ! ! frkapx Elementally raise pressure to the kappa power. ! function result Real(krealfp) pressure in Pascals ! pkap Real(krealfp) p over 1e5 Pa to the kappa power ! ! gtlcl Compute LCL temperature table ! ! ftlcl Elementally compute LCL temperature. ! function result Real(krealfp) temperature at the LCL in Kelvin ! t Real(krealfp) temperature in Kelvin ! tdpd Real(krealfp) dewpoint depression in Kelvin ! ! ftlclq Elementally compute LCL temperature. ! function result Real(krealfp) temperature at the LCL in Kelvin ! t Real(krealfp) temperature in Kelvin ! tdpd Real(krealfp) dewpoint depression in Kelvin ! ! ftlclo Elementally compute LCL temperature. ! function result Real(krealfp) temperature at the LCL in Kelvin ! t Real(krealfp) temperature in Kelvin ! tdpd Real(krealfp) dewpoint depression in Kelvin ! ! ftlclx Elementally compute LCL temperature. ! function result Real(krealfp) temperature at the LCL in Kelvin ! t Real(krealfp) temperature in Kelvin ! tdpd Real(krealfp) dewpoint depression in Kelvin ! ! gfuncphys Compute all physics function tables ! ! Attributes: ! Language: Fortran 90 ! !$$$ use machine,only:kind_phys,r8=>kind_dbl_prec,r4=>kind_sngl_prec use physcons implicit none private ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! Public Variables ! integer,public,parameter:: krealfp=selected_real_kind(15,45) integer,public,parameter:: krealfp=kind_phys !< Integer parameter kind or length of reals ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! Private Variables real(krealfp),parameter:: psatb=con_psat*1.e-5 integer,parameter:: nxpvsl=7501 real(krealfp) c1xpvsl,c2xpvsl,tbpvsl(nxpvsl) integer,parameter:: nxpvsi=7501 real(krealfp) c1xpvsi,c2xpvsi,tbpvsi(nxpvsi) integer,parameter:: nxpvs=7501 real(krealfp) c1xpvs,c2xpvs,tbpvs(nxpvs) integer,parameter:: nxtdpl=5001 real(krealfp) c1xtdpl,c2xtdpl,tbtdpl(nxtdpl) integer,parameter:: nxtdpi=5001 real(krealfp) c1xtdpi,c2xtdpi,tbtdpi(nxtdpi) integer,parameter:: nxtdp=5001 real(krealfp) c1xtdp,c2xtdp,tbtdp(nxtdp) integer,parameter:: nxthe=241,nythe=151 real(krealfp) c1xthe,c2xthe,c1ythe,c2ythe,tbthe(nxthe,nythe) integer,parameter:: nxma=151,nyma=121 real(krealfp) c1xma,c2xma,c1yma,c2yma,tbtma(nxma,nyma),tbqma(nxma,nyma) ! integer,parameter:: nxpkap=5501 integer,parameter:: nxpkap=11001 real(krealfp) c1xpkap,c2xpkap,tbpkap(nxpkap) integer,parameter:: nxrkap=5501 real(krealfp) c1xrkap,c2xrkap,tbrkap(nxrkap) integer,parameter:: nxtlcl=151,nytlcl=61 real(krealfp) c1xtlcl,c2xtlcl,c1ytlcl,c2ytlcl,tbtlcl(nxtlcl,nytlcl) ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! Public Subprograms public gpvsl,fpvsl,fpvslq,fpvslx public gpvsi,fpvsi,fpvsiq,fpvsix public gpvs,fpvs,fpvsq,fpvsx public gtdpl,ftdpl,ftdplq,ftdplx,ftdplxg public gtdpi,ftdpi,ftdpiq,ftdpix,ftdpixg public gtdp,ftdp,ftdpq,ftdpx,ftdpxg public gthe,fthe,ftheq,fthex public gtma,stma,stmaq,stmax,stmaxg public gpkap,fpkap,fpkapq,fpkapo,fpkapx public grkap,frkap,frkapq,frkapx public gtlcl,ftlcl,ftlclq,ftlclo,ftlclx public gfuncphys interface fpvsl module procedure fpvsl_r4, fpvsl_r8 end interface fpvsl interface fpvsi module procedure fpvsi_r4, fpvsi_r8 end interface fpvsi contains !------------------------------------------------------------------------------- !> This subroutine computes saturation vapor pressure table as a function of !! temperature for the table lookup function fpval. Exact saturation vapor !! pressures are calculated in subprogram fpvslx(). The current implementation !! computes a table with a length of 7501 for temperature ranging from 180. to !! 330. Kelvin. subroutine gpvsl !$$$ Subprogram Documentation Block ! ! Subprogram: gpvsl Compute saturation vapor pressure table over liquid ! Author: N Phillips W/NMC2X2 Date: 30 dec 82 ! ! Abstract: Computes saturation vapor pressure table as a function of ! temperature for the table lookup function fpvsl. ! Exact saturation vapor pressures are calculated in subprogram fpvslx. ! The current implementation computes a table with a length ! of 7501 for temperatures ranging from 180. to 330. Kelvin. ! ! Program History Log: ! 91-05-07 Iredell ! 94-12-30 Iredell expand table ! 1999-03-01 Iredell f90 module ! ! Usage: call gpvsl ! ! Subprograms called: ! (fpvslx) inlinable function to compute saturation vapor pressure ! ! Attributes: ! Language: Fortran 90. ! !$$$ implicit none integer jx real(krealfp) xmin,xmax,xinc,x,t ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - xmin=180.0_krealfp xmax=330.0_krealfp xinc=(xmax-xmin)/(nxpvsl-1) ! c1xpvsl=1.-xmin/xinc c2xpvsl=1./xinc c1xpvsl=1.-xmin*c2xpvsl do jx=1,nxpvsl x=xmin+(jx-1)*xinc t=x tbpvsl(jx)=fpvslx(t) enddo ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end subroutine !------------------------------------------------------------------------------- !> This funtion computes saturation vapor pressure from the temperature. !! A linear interpolation is done between values in a lookup table computed !! in gpvsl(). See documentation for fpvslx() for details. Input values !! outside table range are reset to table extrema. !>\author N phillips elemental function fpvsl_r4(t) !$$$ Subprogram Documentation Block ! ! Subprogram: fpvsl Compute saturation vapor pressure over liquid ! Author: N Phillips w/NMC2X2 Date: 30 dec 82 ! ! Abstract: Compute saturation vapor pressure from the temperature. ! A linear interpolation is done between values in a lookup table ! computed in gpvsl. See documentation for fpvslx for details. ! Input values outside table range are reset to table extrema. ! The interpolation accuracy is almost 6 decimal places. ! On the Cray, fpvsl is about 4 times faster than exact calculation. ! This function should be expanded inline in the calling routine. ! ! Program History Log: ! 91-05-07 Iredell made into inlinable function ! 94-12-30 Iredell expand table ! 1999-03-01 Iredell f90 module ! ! Usage: pvsl=fpvsl(t) ! ! Input argument list: ! t Real(krealfp) temperature in Kelvin ! ! Output argument list: ! fpvsl Real(krealfp) saturation vapor pressure in Pascals ! ! Attributes: ! Language: Fortran 90. ! !$$$ implicit none real(r4) fpvsl_r4 real(r4),intent(in):: t integer jx real(r4) xj ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - xj=min(max(c1xpvsl+c2xpvsl*t,1._r4),real(nxpvsl,r4)) jx=min(xj,nxpvsl-1._r4) fpvsl_r4=tbpvsl(jx)+(xj-jx)*(tbpvsl(jx+1)-tbpvsl(jx)) ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end function fpvsl_r4 elemental function fpvsl_r8(t) !$$$ Subprogram Documentation Block ! ! Subprogram: fpvsl Compute saturation vapor pressure over liquid ! Author: N Phillips w/NMC2X2 Date: 30 dec 82 ! ! Abstract: Compute saturation vapor pressure from the temperature. ! A linear interpolation is done between values in a lookup table ! computed in gpvsl. See documentation for fpvslx for details. ! Input values outside table range are reset to table extrema. ! The interpolation accuracy is almost 6 decimal places. ! On the Cray, fpvsl is about 4 times faster than exact calculation. ! This function should be expanded inline in the calling routine. ! ! Program History Log: ! 91-05-07 Iredell made into inlinable function ! 94-12-30 Iredell expand table ! 1999-03-01 Iredell f90 module ! ! Usage: pvsl=fpvsl(t) ! ! Input argument list: ! t Real(krealfp) temperature in Kelvin ! ! Output argument list: ! fpvsl Real(krealfp) saturation vapor pressure in Pascals ! ! Attributes: ! Language: Fortran 90. ! !$$$ implicit none real(r8) fpvsl_r8 real(r8),intent(in):: t integer jx real(r8) xj ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - xj=min(max(c1xpvsl+c2xpvsl*t,1._r8),real(nxpvsl,r8)) jx=min(xj,nxpvsl-1._r8) fpvsl_r8=tbpvsl(jx)+(xj-jx)*(tbpvsl(jx+1)-tbpvsl(jx)) ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end function fpvsl_r8 !------------------------------------------------------------------------------- !> This function computes saturation vapor pressure from the temperature. !! A quadratic interpolation is done between values in a lookup table !! computed in gpvsl(). See documentaion for fpvslx() for details. !! Input values outside table range are reset to table extrema. elemental function fpvslq(t) !>\author N Phillips !$$$ Subprogram Documentation Block ! ! Subprogram: fpvslq Compute saturation vapor pressure over liquid ! Author: N Phillips w/NMC2X2 Date: 30 dec 82 ! ! Abstract: Compute saturation vapor pressure from the temperature. ! A quadratic interpolation is done between values in a lookup table ! computed in gpvsl. See documentation for fpvslx for details. ! Input values outside table range are reset to table extrema. ! The interpolation accuracy is almost 9 decimal places. ! On the Cray, fpvslq is about 3 times faster than exact calculation. ! This function should be expanded inline in the calling routine. ! ! Program History Log: ! 91-05-07 Iredell made into inlinable function ! 94-12-30 Iredell quadratic interpolation ! 1999-03-01 Iredell f90 module ! ! Usage: pvsl=fpvslq(t) ! ! Input argument list: ! t Real(krealfp) temperature in Kelvin ! ! Output argument list: ! fpvslq Real(krealfp) saturation vapor pressure in Pascals ! ! Attributes: ! Language: Fortran 90. ! !$$$ implicit none real(krealfp) fpvslq real(krealfp),intent(in):: t integer jx real(krealfp) xj,dxj,fj1,fj2,fj3 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - xj=min(max(c1xpvsl+c2xpvsl*t,1._krealfp),real(nxpvsl,krealfp)) jx=min(max(nint(xj),2),nxpvsl-1) dxj=xj-jx fj1=tbpvsl(jx-1) fj2=tbpvsl(jx) fj3=tbpvsl(jx+1) fpvslq=(((fj3+fj1)/2-fj2)*dxj+(fj3-fj1)/2)*dxj+fj2 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end function !------------------------------------------------------------------------------- !> This function exactly computes saturation vapor pressure from temperature. !! The water model assumes a perfect gas, constant specific heats !! for gas and liquid, and neglects the volume of the liquid. !! The model does account for the variation of the latent heat !! of condensation with temperature. The ice option is not included. !! The Clausius-Clapeyron equation is integrated from the triple point !! to get the formula: !!\n pvsl=con_psat*(tr**xa)*exp(xb*(1.-tr)) !!\n where tr is ttp/t and other values are physical constants. !! This function should be expanded inline in the calling routine. !>\author N Phillips !>\param[in] t real, temperature in Kelvin !\param[out] fpvslx real, saturation vapor pressure in Pascals elemental function fpvslx(t) !$$$ Subprogram Documentation Block ! ! Subprogram: fpvslx Compute saturation vapor pressure over liquid ! Author: N Phillips w/NMC2X2 Date: 30 dec 82 ! ! Abstract: Exactly compute saturation vapor pressure from temperature. ! The water model assumes a perfect gas, constant specific heats ! for gas and liquid, and neglects the volume of the liquid. ! The model does account for the variation of the latent heat ! of condensation with temperature. The ice option is not included. ! The Clausius-Clapeyron equation is integrated from the triple point ! to get the formula ! pvsl=con_psat*(tr**xa)*exp(xb*(1.-tr)) ! where tr is ttp/t and other values are physical constants. ! This function should be expanded inline in the calling routine. ! ! Program History Log: ! 91-05-07 Iredell made into inlinable function ! 94-12-30 Iredell exact computation ! 1999-03-01 Iredell f90 module ! ! Usage: pvsl=fpvslx(t) ! ! Input argument list: ! t Real(krealfp) temperature in Kelvin ! ! Output argument list: ! fpvslx Real(krealfp) saturation vapor pressure in Pascals ! ! Attributes: ! Language: Fortran 90. ! !$$$ implicit none real(krealfp) fpvslx real(krealfp),intent(in):: t real(krealfp),parameter:: dldt=con_cvap-con_cliq real(krealfp),parameter:: heat=con_hvap real(krealfp),parameter:: xpona=-dldt/con_rv real(krealfp),parameter:: xponb=-dldt/con_rv+heat/(con_rv*con_ttp) real(krealfp) tr ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - tr=con_ttp/t fpvslx=con_psat*(tr**xpona)*exp(xponb*(1.-tr)) ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end function !------------------------------------------------------------------------------- !> This subroutine computes saturation vapor pressure table as a function of !! temperature for the table lookup function fpvsi(). Exact saturation vapor !! pressures are calculated in subprogram fpvsix(). The current implementation !! computes a table with a length of 7501 for temperatures ranging from 180. !! to 330. Kelvin. !>\author N Phillips subroutine gpvsi !$$$ Subprogram Documentation Block ! ! Subprogram: gpvsi Compute saturation vapor pressure table over ice ! Author: N Phillips W/NMC2X2 Date: 30 dec 82 ! ! Abstract: Computes saturation vapor pressure table as a function of ! temperature for the table lookup function fpvsi. ! Exact saturation vapor pressures are calculated in subprogram fpvsix. ! The current implementation computes a table with a length ! of 7501 for temperatures ranging from 180. to 330. Kelvin. ! ! Program History Log: ! 91-05-07 Iredell ! 94-12-30 Iredell expand table ! 1999-03-01 Iredell f90 module ! 2001-02-26 Iredell ice phase ! ! Usage: call gpvsi ! ! Subprograms called: ! (fpvsix) inlinable function to compute saturation vapor pressure ! ! Attributes: ! Language: Fortran 90. ! !$$$ implicit none integer jx real(krealfp) xmin,xmax,xinc,x,t ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - xmin=180.0_krealfp xmax=330.0_krealfp xinc=(xmax-xmin)/(nxpvsi-1) ! c1xpvsi=1.-xmin/xinc c2xpvsi=1./xinc c1xpvsi=1.-xmin*c2xpvsi do jx=1,nxpvsi x=xmin+(jx-1)*xinc t=x tbpvsi(jx)=fpvsix(t) enddo ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end subroutine !------------------------------------------------------------------------------- !> This function computes saturation vapor pressure from the temperature. !! A linear interpolation is done between values in a lookup table !! computed in gpvsi(). See documentation for fpvsix() for details. !! Input values outside table range are reset to table extrema. !>\author N Phillips elemental function fpvsi_r4(t) !$$$ Subprogram Documentation Block ! ! Subprogram: fpvsi Compute saturation vapor pressure over ice ! Author: N Phillips w/NMC2X2 Date: 30 dec 82 ! ! Abstract: Compute saturation vapor pressure from the temperature. ! A linear interpolation is done between values in a lookup table ! computed in gpvsi. See documentation for fpvsix for details. ! Input values outside table range are reset to table extrema. ! The interpolation accuracy is almost 6 decimal places. ! On the Cray, fpvsi is about 4 times faster than exact calculation. ! This function should be expanded inline in the calling routine. ! ! Program History Log: ! 91-05-07 Iredell made into inlinable function ! 94-12-30 Iredell expand table ! 1999-03-01 Iredell f90 module ! 2001-02-26 Iredell ice phase ! ! Usage: pvsi=fpvsi(t) ! ! Input argument list: ! t Real(krealfp) temperature in Kelvin ! ! Output argument list: ! fpvsi Real(krealfp) saturation vapor pressure in Pascals ! ! Attributes: ! Language: Fortran 90. ! !$$$ implicit none real(r4) fpvsi_r4 real(r4),intent(in):: t integer jx real(r4) xj ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - xj=min(max(c1xpvsi+c2xpvsi*t,1._r4),real(nxpvsi,r4)) jx=min(xj,nxpvsi-1._r4) fpvsi_r4=tbpvsi(jx)+(xj-jx)*(tbpvsi(jx+1)-tbpvsi(jx)) ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end function fpvsi_r4 elemental function fpvsi_r8(t) !$$$ Subprogram Documentation Block ! ! Subprogram: fpvsi Compute saturation vapor pressure over ice ! Author: N Phillips w/NMC2X2 Date: 30 dec 82 ! ! Abstract: Compute saturation vapor pressure from the temperature. ! A linear interpolation is done between values in a lookup table ! computed in gpvsi. See documentation for fpvsix for details. ! Input values outside table range are reset to table extrema. ! The interpolation accuracy is almost 6 decimal places. ! On the Cray, fpvsi is about 4 times faster than exact calculation. ! This function should be expanded inline in the calling routine. ! ! Program History Log: ! 91-05-07 Iredell made into inlinable function ! 94-12-30 Iredell expand table ! 1999-03-01 Iredell f90 module ! 2001-02-26 Iredell ice phase ! ! Usage: pvsi=fpvsi(t) ! ! Input argument list: ! t Real(krealfp) temperature in Kelvin ! ! Output argument list: ! fpvsi Real(krealfp) saturation vapor pressure in Pascals ! ! Attributes: ! Language: Fortran 90. ! !$$$ implicit none real(r8) fpvsi_r8 real(r8),intent(in):: t integer jx real(r8) xj ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - xj=min(max(c1xpvsi+c2xpvsi*t,1._r8),real(nxpvsi,r8)) jx=min(xj,nxpvsi-1._r8) fpvsi_r8=tbpvsi(jx)+(xj-jx)*(tbpvsi(jx+1)-tbpvsi(jx)) ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end function fpvsi_r8 !------------------------------------------------------------------------------- !> This function computes saturation vapor pressure from the temperature. !! A quadratic interpolation is done between values in a lookup table !! computed in gpvsi(). See documentation for fpvsix() for details. !! Input values outside table range are reset to table extrema. !>\author N Phillips elemental function fpvsiq(t) !$$$ Subprogram Documentation Block ! ! Subprogram: fpvsiq Compute saturation vapor pressure over ice ! Author: N Phillips w/NMC2X2 Date: 30 dec 82 ! ! Abstract: Compute saturation vapor pressure from the temperature. ! A quadratic interpolation is done between values in a lookup table ! computed in gpvsi. See documentation for fpvsix for details. ! Input values outside table range are reset to table extrema. ! The interpolation accuracy is almost 9 decimal places. ! On the Cray, fpvsiq is about 3 times faster than exact calculation. ! This function should be expanded inline in the calling routine. ! ! Program History Log: ! 91-05-07 Iredell made into inlinable function ! 94-12-30 Iredell quadratic interpolation ! 1999-03-01 Iredell f90 module ! 2001-02-26 Iredell ice phase ! ! Usage: pvsi=fpvsiq(t) ! ! Input argument list: ! t Real(krealfp) temperature in Kelvin ! ! Output argument list: ! fpvsiq Real(krealfp) saturation vapor pressure in Pascals ! ! Attributes: ! Language: Fortran 90. ! !$$$ implicit none real(krealfp) fpvsiq real(krealfp),intent(in):: t integer jx real(krealfp) xj,dxj,fj1,fj2,fj3 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - xj=min(max(c1xpvsi+c2xpvsi*t,1._krealfp),real(nxpvsi,krealfp)) jx=min(max(nint(xj),2),nxpvsi-1) dxj=xj-jx fj1=tbpvsi(jx-1) fj2=tbpvsi(jx) fj3=tbpvsi(jx+1) fpvsiq=(((fj3+fj1)/2-fj2)*dxj+(fj3-fj1)/2)*dxj+fj2 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end function !------------------------------------------------------------------------------- !> This funtion exactly computes saturation vapor pressure from temperature. !! The water model assumes a perfect gas, constant specific heats !! for gas and ice, and neglects the volume of the ice. The model does !! account for the variation of the latent heat of condensation with temperature. !! The liquid option is not included. The Clausius- Clapeyron equation is !! integrated from the triple point to get the formula: !!\n pvsi=con_psat*(tr**xa)*exp(xb*(1.-tr)) !!\n where tr is ttp/t and other values are physical constants. !! This function should be expanded inline in the calling routine. !>\param[in] t real, temperature in Kelvin !\param[out] fpvsix real, saturation vapor pressure in Pascals elemental function fpvsix(t) !$$$ Subprogram Documentation Block ! ! Subprogram: fpvsix Compute saturation vapor pressure over ice ! Author: N Phillips w/NMC2X2 Date: 30 dec 82 ! ! Abstract: Exactly compute saturation vapor pressure from temperature. ! The water model assumes a perfect gas, constant specific heats ! for gas and ice, and neglects the volume of the ice. ! The model does account for the variation of the latent heat ! of condensation with temperature. The liquid option is not included. ! The Clausius-Clapeyron equation is integrated from the triple point ! to get the formula ! pvsi=con_psat*(tr**xa)*exp(xb*(1.-tr)) ! where tr is ttp/t and other values are physical constants. ! This function should be expanded inline in the calling routine. ! ! Program History Log: ! 91-05-07 Iredell made into inlinable function ! 94-12-30 Iredell exact computation ! 1999-03-01 Iredell f90 module ! 2001-02-26 Iredell ice phase ! ! Usage: pvsi=fpvsix(t) ! ! Input argument list: ! t Real(krealfp) temperature in Kelvin ! ! Output argument list: ! fpvsix Real(krealfp) saturation vapor pressure in Pascals ! ! Attributes: ! Language: Fortran 90. ! !$$$ implicit none real(krealfp) fpvsix real(krealfp),intent(in):: t real(krealfp),parameter:: dldt=con_cvap-con_csol real(krealfp),parameter:: heat=con_hvap+con_hfus real(krealfp),parameter:: xpona=-dldt/con_rv real(krealfp),parameter:: xponb=-dldt/con_rv+heat/(con_rv*con_ttp) real(krealfp) tr ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - tr=con_ttp/t fpvsix=con_psat*(tr**xpona)*exp(xponb*(1.-tr)) ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end function !------------------------------------------------------------------------------- !> This subroutine computes saturation vapor pressure table as a function of !! temperature for the table lookup function fpvs(). !! Exact saturation vapor pressures are calculated in subprogram fpvsx(). !! The current implementation computes a table with a length !! of 7501 for temperatures ranging from 180. to 330. Kelvin. subroutine gpvs !$$$ Subprogram Documentation Block ! ! Subprogram: gpvs Compute saturation vapor pressure table ! Author: N Phillips W/NMC2X2 Date: 30 dec 82 ! ! Abstract: Computes saturation vapor pressure table as a function of ! temperature for the table lookup function fpvs. ! Exact saturation vapor pressures are calculated in subprogram fpvsx. ! The current implementation computes a table with a length ! of 7501 for temperatures ranging from 180. to 330. Kelvin. ! ! Program History Log: ! 91-05-07 Iredell ! 94-12-30 Iredell expand table ! 1999-03-01 Iredell f90 module ! 2001-02-26 Iredell ice phase ! ! Usage: call gpvs ! ! Subprograms called: ! (fpvsx) inlinable function to compute saturation vapor pressure ! ! Attributes: ! Language: Fortran 90. ! !$$$ implicit none integer jx real(krealfp) xmin,xmax,xinc,x,t ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - xmin=180.0_krealfp xmax=330.0_krealfp xinc=(xmax-xmin)/(nxpvs-1) ! c1xpvs=1.-xmin/xinc c2xpvs=1./xinc c1xpvs=1.-xmin*c2xpvs do jx=1,nxpvs x=xmin+(jx-1)*xinc t=x tbpvs(jx)=fpvsx(t) enddo ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end subroutine !------------------------------------------------------------------------------- !> This function computes saturation vapor pressure from the temperature. !! A linear interpolation is done between values in a lookup table !! computed in gpvs(). See documentation for fpvsx() for details. !! Input values outside table range are reset to table extrema. !>\param[in] t real, temperature in Kelvin !\param[out] fpvs real, saturation vapor pressure in Pascals elemental function fpvs(t) !$$$ Subprogram Documentation Block ! ! Subprogram: fpvs Compute saturation vapor pressure ! Author: N Phillips w/NMC2X2 Date: 30 dec 82 ! ! Abstract: Compute saturation vapor pressure from the temperature. ! A linear interpolation is done between values in a lookup table ! computed in gpvs. See documentation for fpvsx for details. ! Input values outside table range are reset to table extrema. ! The interpolation accuracy is almost 6 decimal places. ! On the Cray, fpvs is about 4 times faster than exact calculation. ! This function should be expanded inline in the calling routine. ! ! Program History Log: ! 91-05-07 Iredell made into inlinable function ! 94-12-30 Iredell expand table ! 1999-03-01 Iredell f90 module ! 2001-02-26 Iredell ice phase ! ! Usage: pvs=fpvs(t) ! ! Input argument list: ! t Real(krealfp) temperature in Kelvin ! ! Output argument list: ! fpvs Real(krealfp) saturation vapor pressure in Pascals ! ! Attributes: ! Language: Fortran 90. ! !$$$ implicit none real(krealfp) fpvs real(krealfp),intent(in):: t integer jx real(krealfp) xj ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - xj=min(max(c1xpvs+c2xpvs*t,1._krealfp),real(nxpvs,krealfp)) jx=min(xj,nxpvs-1._krealfp) fpvs=tbpvs(jx)+(xj-jx)*(tbpvs(jx+1)-tbpvs(jx)) ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end function !------------------------------------------------------------------------------- !> This function computes saturation vapor pressure from the temperature. !! A quadratic interpolation is done between values in a lookup table !! computed in gpvs(). See documentation for fpvsx() for details. !! Input values outside table range are reset to table extrema. !>\param[in] t real, temperatue in Kelvin !\param[out] fpvsq real, saturation vapor pressure in Pascals elemental function fpvsq(t) !$$$ Subprogram Documentation Block ! ! Subprogram: fpvsq Compute saturation vapor pressure ! Author: N Phillips w/NMC2X2 Date: 30 dec 82 ! ! Abstract: Compute saturation vapor pressure from the temperature. ! A quadratic interpolation is done between values in a lookup table ! computed in gpvs. See documentation for fpvsx for details. ! Input values outside table range are reset to table extrema. ! The interpolation accuracy is almost 9 decimal places. ! On the Cray, fpvsq is about 3 times faster than exact calculation. ! This function should be expanded inline in the calling routine. ! ! Program History Log: ! 91-05-07 Iredell made into inlinable function ! 94-12-30 Iredell quadratic interpolation ! 1999-03-01 Iredell f90 module ! 2001-02-26 Iredell ice phase ! ! Usage: pvs=fpvsq(t) ! ! Input argument list: ! t Real(krealfp) temperature in Kelvin ! ! Output argument list: ! fpvsq Real(krealfp) saturation vapor pressure in Pascals ! ! Attributes: ! Language: Fortran 90. ! !$$$ implicit none real(krealfp) fpvsq real(krealfp),intent(in):: t integer jx real(krealfp) xj,dxj,fj1,fj2,fj3 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - xj=min(max(c1xpvs+c2xpvs*t,1._krealfp),real(nxpvs,krealfp)) jx=min(max(nint(xj),2),nxpvs-1) dxj=xj-jx fj1=tbpvs(jx-1) fj2=tbpvs(jx) fj3=tbpvs(jx+1) fpvsq=(((fj3+fj1)/2-fj2)*dxj+(fj3-fj1)/2)*dxj+fj2 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end function !------------------------------------------------------------------------------- !> This function exactly computes saturation vapor pressure from temperature. !! The saturation vapor pressure over either liquid and ice is computed !! over liquid for temperatures above the triple point, !! over ice for temperatures 20 degress below the triple point, !! and a linear combination of the two for temperatures in between. !! The water model assumes a perfect gas, constant specific heats !! for gas, liquid and ice, and neglects the volume of the condensate. !! The model does account for the variation of the latent heat !! of condensation and sublimation with temperature. !! The Clausius-Clapeyron equation is integrated from the triple point !! to get the formula: !!\n pvsl=con_psat*(tr**xa)*exp(xb*(1.-tr)) !!\n where tr is ttp/t and other values are physical constants. !! The reference for this computation is Emanuel(1994), pages 116-117. !! This function should be expanded inline in the calling routine. !!\param[in] t real, temperature in Kelvin !\param[out] fpvsx real, saturation vapor pressure in Pascals elemental function fpvsx(t) !$$$ Subprogram Documentation Block ! ! Subprogram: fpvsx Compute saturation vapor pressure ! Author: N Phillips w/NMC2X2 Date: 30 dec 82 ! ! Abstract: Exactly compute saturation vapor pressure from temperature. ! The saturation vapor pressure over either liquid and ice is computed ! over liquid for temperatures above the triple point, ! over ice for temperatures 20 degress below the triple point, ! and a linear combination of the two for temperatures in between. ! The water model assumes a perfect gas, constant specific heats ! for gas, liquid and ice, and neglects the volume of the condensate. ! The model does account for the variation of the latent heat ! of condensation and sublimation with temperature. ! The Clausius-Clapeyron equation is integrated from the triple point ! to get the formula ! pvsl=con_psat*(tr**xa)*exp(xb*(1.-tr)) ! where tr is ttp/t and other values are physical constants. ! The reference for this computation is Emanuel(1994), pages 116-117. ! This function should be expanded inline in the calling routine. ! ! Program History Log: ! 91-05-07 Iredell made into inlinable function ! 94-12-30 Iredell exact computation ! 1999-03-01 Iredell f90 module ! 2001-02-26 Iredell ice phase ! ! Usage: pvs=fpvsx(t) ! ! Input argument list: ! t Real(krealfp) temperature in Kelvin ! ! Output argument list: ! fpvsx Real(krealfp) saturation vapor pressure in Pascals ! ! Attributes: ! Language: Fortran 90. ! !$$$ implicit none real(krealfp) fpvsx real(krealfp),intent(in):: t real(krealfp),parameter:: tliq=con_ttp real(krealfp),parameter:: tice=con_ttp-20.0 real(krealfp),parameter:: dldtl=con_cvap-con_cliq real(krealfp),parameter:: heatl=con_hvap real(krealfp),parameter:: xponal=-dldtl/con_rv real(krealfp),parameter:: xponbl=-dldtl/con_rv+heatl/(con_rv*con_ttp) real(krealfp),parameter:: dldti=con_cvap-con_csol real(krealfp),parameter:: heati=con_hvap+con_hfus real(krealfp),parameter:: xponai=-dldti/con_rv real(krealfp),parameter:: xponbi=-dldti/con_rv+heati/(con_rv*con_ttp) real(krealfp) tr,w,pvl,pvi ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - tr=con_ttp/t if(t.ge.tliq) then fpvsx=con_psat*(tr**xponal)*exp(xponbl*(1.-tr)) elseif(t.lt.tice) then fpvsx=con_psat*(tr**xponai)*exp(xponbi*(1.-tr)) else w=(t-tice)/(tliq-tice) pvl=con_psat*(tr**xponal)*exp(xponbl*(1.-tr)) pvi=con_psat*(tr**xponai)*exp(xponbi*(1.-tr)) fpvsx=w*pvl+(1.-w)*pvi endif ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end function !------------------------------------------------------------------------------- !> This subroutine computes dewpoint temperature table as a function of !! vapor pressure for inlinable function ftdpl(). !! Exact dewpoint temperatures are calculated in subprogram ftdplxg(). !! The current implementation computes a table with a length !! of 5001 for vapor pressures ranging from 1 to 10001 Pascals !! giving a dewpoint temperature range of 208 to 319 Kelvin. subroutine gtdpl !$$$ Subprogram Documentation Block ! ! Subprogram: gtdpl Compute dewpoint temperature over liquid table ! Author: N Phillips w/NMC2X2 Date: 30 dec 82 ! ! Abstract: Compute dewpoint temperature table as a function of ! vapor pressure for inlinable function ftdpl. ! Exact dewpoint temperatures are calculated in subprogram ftdplxg. ! The current implementation computes a table with a length ! of 5001 for vapor pressures ranging from 1 to 10001 Pascals ! giving a dewpoint temperature range of 208 to 319 Kelvin. ! ! Program History Log: ! 91-05-07 Iredell ! 94-12-30 Iredell expand table ! 1999-03-01 Iredell f90 module ! ! Usage: call gtdpl ! ! Subprograms called: ! (ftdplxg) inlinable function to compute dewpoint temperature over liquid ! ! Attributes: ! Language: Fortran 90. ! !$$$ implicit none integer jx real(krealfp) xmin,xmax,xinc,t,x,pv ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - xmin=1 xmax=10001 xinc=(xmax-xmin)/(nxtdpl-1) c1xtdpl=1.-xmin/xinc c2xtdpl=1./xinc t=208.0 do jx=1,nxtdpl x=xmin+(jx-1)*xinc pv=x t=ftdplxg(t,pv) tbtdpl(jx)=t enddo ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end subroutine !------------------------------------------------------------------------------- !> This function Compute dewpoint temperature from vapor pressure. !! A linear interpolation is done between values in a lookup table !! computed in gtdpl(). See documentation for ftdplxg() for details. !! Input values outside table range are reset to table extrema. elemental function ftdpl(pv) !$$$ Subprogram Documentation Block ! ! Subprogram: ftdpl Compute dewpoint temperature over liquid ! Author: N Phillips w/NMC2X2 Date: 30 dec 82 ! ! Abstract: Compute dewpoint temperature from vapor pressure. ! A linear interpolation is done between values in a lookup table ! computed in gtdpl. See documentation for ftdplxg for details. ! Input values outside table range are reset to table extrema. ! The interpolation accuracy is better than 0.0005 Kelvin ! for dewpoint temperatures greater than 250 Kelvin, ! but decreases to 0.02 Kelvin for a dewpoint around 230 Kelvin. ! On the Cray, ftdpl is about 75 times faster than exact calculation. ! This function should be expanded inline in the calling routine. ! ! Program History Log: ! 91-05-07 Iredell made into inlinable function ! 94-12-30 Iredell expand table ! 1999-03-01 Iredell f90 module ! ! Usage: tdpl=ftdpl(pv) ! ! Input argument list: ! pv Real(krealfp) vapor pressure in Pascals ! ! Output argument list: ! ftdpl Real(krealfp) dewpoint temperature in Kelvin ! ! Attributes: ! Language: Fortran 90. ! !$$$ implicit none real(krealfp) ftdpl real(krealfp),intent(in):: pv integer jx real(krealfp) xj ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - xj=min(max(c1xtdpl+c2xtdpl*pv,1._krealfp),real(nxtdpl,krealfp)) jx=min(xj,nxtdpl-1._krealfp) ftdpl=tbtdpl(jx)+(xj-jx)*(tbtdpl(jx+1)-tbtdpl(jx)) ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end function !------------------------------------------------------------------------------- !> This function computes dewpoint temperature from vapor pressure. !! A quadratic interpolation is done between values in a lookup table !! computed in gtdpl(). See documentation for ftdplxg() for details. !! Input values outside table range are reset to table extrema. !!\param[in] pv real, vapor pressure in Pascals !\param[out] ftdplq real, dewpoint temperature in Kelvin elemental function ftdplq(pv) !$$$ Subprogram Documentation Block ! ! Subprogram: ftdplq Compute dewpoint temperature over liquid ! Author: N Phillips w/NMC2X2 Date: 30 dec 82 ! ! Abstract: Compute dewpoint temperature from vapor pressure. ! A quadratic interpolation is done between values in a lookup table ! computed in gtdpl. see documentation for ftdplxg for details. ! Input values outside table range are reset to table extrema. ! the interpolation accuracy is better than 0.00001 Kelvin ! for dewpoint temperatures greater than 250 Kelvin, ! but decreases to 0.002 Kelvin for a dewpoint around 230 Kelvin. ! On the Cray, ftdplq is about 60 times faster than exact calculation. ! This function should be expanded inline in the calling routine. ! ! Program History Log: ! 91-05-07 Iredell made into inlinable function ! 94-12-30 Iredell quadratic interpolation ! 1999-03-01 Iredell f90 module ! ! Usage: tdpl=ftdplq(pv) ! ! Input argument list: ! pv Real(krealfp) vapor pressure in Pascals ! ! Output argument list: ! ftdplq Real(krealfp) dewpoint temperature in Kelvin ! ! Attributes: ! Language: Fortran 90. ! !$$$ implicit none real(krealfp) ftdplq real(krealfp),intent(in):: pv integer jx real(krealfp) xj,dxj,fj1,fj2,fj3 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - xj=min(max(c1xtdpl+c2xtdpl*pv,1._krealfp),real(nxtdpl,krealfp)) jx=min(max(nint(xj),2),nxtdpl-1) dxj=xj-jx fj1=tbtdpl(jx-1) fj2=tbtdpl(jx) fj3=tbtdpl(jx+1) ftdplq=(((fj3+fj1)/2-fj2)*dxj+(fj3-fj1)/2)*dxj+fj2 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end function !------------------------------------------------------------------------------- !> This function exactly compute dewpoint temperature from vapor pressure. !! An approximate dewpoint temperature for function ftdplxg() !! is obtained using ftdpl() so gtdpl() must be already called. !! See documentation for ftdplxg() for details. !!\param[in] pv real, vapor pressure in Pascals !\param[out] ftdplx real, dewpoint temperature in Kelvin elemental function ftdplx(pv) !$$$ Subprogram Documentation Block ! ! Subprogram: ftdplx Compute dewpoint temperature over liquid ! Author: N Phillips w/NMC2X2 Date: 30 dec 82 ! ! Abstract: exactly compute dewpoint temperature from vapor pressure. ! An approximate dewpoint temperature for function ftdplxg ! is obtained using ftdpl so gtdpl must be already called. ! See documentation for ftdplxg for details. ! ! Program History Log: ! 91-05-07 Iredell made into inlinable function ! 94-12-30 Iredell exact computation ! 1999-03-01 Iredell f90 module ! ! Usage: tdpl=ftdplx(pv) ! ! Input argument list: ! pv Real(krealfp) vapor pressure in Pascals ! ! Output argument list: ! ftdplx Real(krealfp) dewpoint temperature in Kelvin ! ! Subprograms called: ! (ftdpl) inlinable function to compute dewpoint temperature over liquid ! (ftdplxg) inlinable function to compute dewpoint temperature over liquid ! ! Attributes: ! Language: Fortran 90. ! !$$$ implicit none real(krealfp) ftdplx real(krealfp),intent(in):: pv real(krealfp) tg ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - tg=ftdpl(pv) ftdplx=ftdplxg(tg,pv) ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end function !------------------------------------------------------------------------------- !> This function exactly computes dewpoint temperature from vapor pressure. !! A guess dewpoint temperature must be provided. !! The water model assumes a perfect gas, constant specific heats !! for gas and liquid, and neglects the volume of the liquid. !! The model does account for the variation of the latent heat !! of condensation with temperature. The ice option is not included. !! The Clausius-Clapeyron equation is integrated from the triple point !! to get the formula: !!\n pvs=con_psat*(tr**xa)*exp(xb*(1.-tr)) !!\n where tr is ttp/t and other values are physical constants. !!\n The formula is inverted by iterating Newtonian approximations !! for each pvs until t is found to within 1.e-6 Kelvin. !! This function can be expanded inline in the calling routine. !!\param[in] tg real, guess dewpoint temperature in Kelvin !!\param[in] pv real, vapor pressure in Pascals !\param[out] ftdplxg real, dewpoint temperature in Kelvin elemental function ftdplxg(tg,pv) !$$$ Subprogram Documentation Block ! ! Subprogram: ftdplxg Compute dewpoint temperature over liquid ! Author: N Phillips w/NMC2X2 Date: 30 dec 82 ! ! Abstract: Exactly compute dewpoint temperature from vapor pressure. ! A guess dewpoint temperature must be provided. ! The water model assumes a perfect gas, constant specific heats ! for gas and liquid, and neglects the volume of the liquid. ! The model does account for the variation of the latent heat ! of condensation with temperature. The ice option is not included. ! The Clausius-Clapeyron equation is integrated from the triple point ! to get the formula ! pvs=con_psat*(tr**xa)*exp(xb*(1.-tr)) ! where tr is ttp/t and other values are physical constants. ! The formula is inverted by iterating Newtonian approximations ! for each pvs until t is found to within 1.e-6 Kelvin. ! This function can be expanded inline in the calling routine. ! ! Program History Log: ! 91-05-07 Iredell made into inlinable function ! 94-12-30 Iredell exact computation ! 1999-03-01 Iredell f90 module ! ! Usage: tdpl=ftdplxg(tg,pv) ! ! Input argument list: ! tg Real(krealfp) guess dewpoint temperature in Kelvin ! pv Real(krealfp) vapor pressure in Pascals ! ! Output argument list: ! ftdplxg Real(krealfp) dewpoint temperature in Kelvin ! ! Attributes: ! Language: Fortran 90. ! !$$$ implicit none real(krealfp) ftdplxg real(krealfp),intent(in):: tg,pv real(krealfp),parameter:: terrm=1.e-6 real(krealfp),parameter:: dldt=con_cvap-con_cliq real(krealfp),parameter:: heat=con_hvap real(krealfp),parameter:: xpona=-dldt/con_rv real(krealfp),parameter:: xponb=-dldt/con_rv+heat/(con_rv*con_ttp) real(krealfp) t,tr,pvt,el,dpvt,terr integer i ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - t=tg do i=1,100 tr=con_ttp/t pvt=con_psat*(tr**xpona)*exp(xponb*(1.-tr)) el=heat+dldt*(t-con_ttp) dpvt=el*pvt/(con_rv*t**2) terr=(pvt-pv)/dpvt t=t-terr if(abs(terr).le.terrm) exit enddo ftdplxg=t ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end function !------------------------------------------------------------------------------- !> This subroutine computes dewpoint temperature table as a function of !! vapor pressure for inlinable function ftdpi(). !! Exact dewpoint temperatures are calculated in subprogram ftdpixg(). !! The current implementation computes a table with a length !! of 5001 for vapor pressures ranging from 0.1 to 1000.1 Pascals !! giving a dewpoint temperature range of 197 to 279 Kelvin. subroutine gtdpi !$$$ Subprogram Documentation Block ! ! Subprogram: gtdpi Compute dewpoint temperature over ice table ! Author: N Phillips w/NMC2X2 Date: 30 dec 82 ! ! Abstract: Compute dewpoint temperature table as a function of ! vapor pressure for inlinable function ftdpi. ! Exact dewpoint temperatures are calculated in subprogram ftdpixg. ! The current implementation computes a table with a length ! of 5001 for vapor pressures ranging from 0.1 to 1000.1 Pascals ! giving a dewpoint temperature range of 197 to 279 Kelvin. ! ! Program History Log: ! 91-05-07 Iredell ! 94-12-30 Iredell expand table ! 1999-03-01 Iredell f90 module ! 2001-02-26 Iredell ice phase ! ! Usage: call gtdpi ! ! Subprograms called: ! (ftdpixg) inlinable function to compute dewpoint temperature over ice ! ! Attributes: ! Language: Fortran 90. ! !$$$ implicit none integer jx real(krealfp) xmin,xmax,xinc,t,x,pv ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - xmin=0.1 xmax=1000.1 xinc=(xmax-xmin)/(nxtdpi-1) c1xtdpi=1.-xmin/xinc c2xtdpi=1./xinc t=197.0 do jx=1,nxtdpi x=xmin+(jx-1)*xinc pv=x t=ftdpixg(t,pv) tbtdpi(jx)=t enddo ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end subroutine !------------------------------------------------------------------------------- !> This subroutine computes dewpoint temperature from vapor pressure. !! A linear interpolation is done between values in a lookup table !! computed in gtdpi(). See documentation for ftdpixg for details. !! Input values outside table range are reset to table extrema. !>\param[in] pv real, vapor pressure in Pascals !\param[out] ftdpi real, dewpoint temperature in Kelvin elemental function ftdpi(pv) !$$$ Subprogram Documentation Block ! ! Subprogram: ftdpi Compute dewpoint temperature over ice ! Author: N Phillips w/NMC2X2 Date: 30 dec 82 ! ! Abstract: Compute dewpoint temperature from vapor pressure. ! A linear interpolation is done between values in a lookup table ! computed in gtdpi. See documentation for ftdpixg for details. ! Input values outside table range are reset to table extrema. ! The interpolation accuracy is better than 0.0005 Kelvin ! for dewpoint temperatures greater than 250 Kelvin, ! but decreases to 0.02 Kelvin for a dewpoint around 230 Kelvin. ! On the Cray, ftdpi is about 75 times faster than exact calculation. ! This function should be expanded inline in the calling routine. ! ! Program History Log: ! 91-05-07 Iredell made into inlinable function ! 94-12-30 Iredell expand table ! 1999-03-01 Iredell f90 module ! 2001-02-26 Iredell ice phase ! ! Usage: tdpi=ftdpi(pv) ! ! Input argument list: ! pv Real(krealfp) vapor pressure in Pascals ! ! Output argument list: ! ftdpi Real(krealfp) dewpoint temperature in Kelvin ! ! Attributes: ! Language: Fortran 90. ! !$$$ implicit none real(krealfp) ftdpi real(krealfp),intent(in):: pv integer jx real(krealfp) xj ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - xj=min(max(c1xtdpi+c2xtdpi*pv,1._krealfp),real(nxtdpi,krealfp)) jx=min(xj,nxtdpi-1._krealfp) ftdpi=tbtdpi(jx)+(xj-jx)*(tbtdpi(jx+1)-tbtdpi(jx)) ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end function !------------------------------------------------------------------------------- !> This function computes dewpoint temperature from vapor pressure. !! A quadratic interpolation is done between values in a lookup table !! computed in gtdpi(). see documentation for ftdpixg() for details. !! Input values outside table range are reset to table extrema. !>\param[in] pv real, vapor pressure in Pascals !\param[out] ftdpiq real, dewpoint temperature in Kelvin elemental function ftdpiq(pv) !$$$ Subprogram Documentation Block ! ! Subprogram: ftdpiq Compute dewpoint temperature over ice ! Author: N Phillips w/NMC2X2 Date: 30 dec 82 ! ! Abstract: Compute dewpoint temperature from vapor pressure. ! A quadratic interpolation is done between values in a lookup table ! computed in gtdpi. see documentation for ftdpixg for details. ! Input values outside table range are reset to table extrema. ! the interpolation accuracy is better than 0.00001 Kelvin ! for dewpoint temperatures greater than 250 Kelvin, ! but decreases to 0.002 Kelvin for a dewpoint around 230 Kelvin. ! On the Cray, ftdpiq is about 60 times faster than exact calculation. ! This function should be expanded inline in the calling routine. ! ! Program History Log: ! 91-05-07 Iredell made into inlinable function ! 94-12-30 Iredell quadratic interpolation ! 1999-03-01 Iredell f90 module ! 2001-02-26 Iredell ice phase ! ! Usage: tdpi=ftdpiq(pv) ! ! Input argument list: ! pv Real(krealfp) vapor pressure in Pascals ! ! Output argument list: ! ftdpiq Real(krealfp) dewpoint temperature in Kelvin ! ! Attributes: ! Language: Fortran 90. ! !$$$ implicit none real(krealfp) ftdpiq real(krealfp),intent(in):: pv integer jx real(krealfp) xj,dxj,fj1,fj2,fj3 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - xj=min(max(c1xtdpi+c2xtdpi*pv,1._krealfp),real(nxtdpi,krealfp)) jx=min(max(nint(xj),2),nxtdpi-1) dxj=xj-jx fj1=tbtdpi(jx-1) fj2=tbtdpi(jx) fj3=tbtdpi(jx+1) ftdpiq=(((fj3+fj1)/2-fj2)*dxj+(fj3-fj1)/2)*dxj+fj2 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end function !------------------------------------------------------------------------------- !> This function exactly computes dewpoint temperature from vapor pressure. !! An approximate dewpoint temperature for function ftdpixg() !! is obtained using ftdpi() so gtdpi() must be already called. !! See documentation for ftdpixg() for details. !>\param[in] pv real, vapor pressure in Pascals !\param[out] ftdpix real, dewpoint temperature in Kelvin elemental function ftdpix(pv) !$$$ Subprogram Documentation Block ! ! Subprogram: ftdpix Compute dewpoint temperature over ice ! Author: N Phillips w/NMC2X2 Date: 30 dec 82 ! ! Abstract: exactly compute dewpoint temperature from vapor pressure. ! An approximate dewpoint temperature for function ftdpixg ! is obtained using ftdpi so gtdpi must be already called. ! See documentation for ftdpixg for details. ! ! Program History Log: ! 91-05-07 Iredell made into inlinable function ! 94-12-30 Iredell exact computation ! 1999-03-01 Iredell f90 module ! 2001-02-26 Iredell ice phase ! ! Usage: tdpi=ftdpix(pv) ! ! Input argument list: ! pv Real(krealfp) vapor pressure in Pascals ! ! Output argument list: ! ftdpix Real(krealfp) dewpoint temperature in Kelvin ! ! Subprograms called: ! (ftdpi) inlinable function to compute dewpoint temperature over ice ! (ftdpixg) inlinable function to compute dewpoint temperature over ice ! ! Attributes: ! Language: Fortran 90. ! !$$$ implicit none real(krealfp) ftdpix real(krealfp),intent(in):: pv real(krealfp) tg ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - tg=ftdpi(pv) ftdpix=ftdpixg(tg,pv) ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end function !------------------------------------------------------------------------------- !> This function exactly computes dewpoint temperature from vapor pressure. !! A guess dewpoint temperature must be provided. !! The water model assumes a perfect gas, constant specific heats !! for gas and ice, and neglects the volume of the ice. !! The model does account for the variation of the latent heat !! of sublimation with temperature. The liquid option is not included. !! The Clausius-Clapeyron equation is integrated from the triple point !! to get the formula: !!\n pvs=con_psat*(tr**xa)*exp(xb*(1.-tr)) !!\n where tr is ttp/t and other values are physical constants. !! The formula is inverted by iterating Newtonian approximations !! for each pvs until t is found to within 1.e-6 Kelvin. !! This function can be expanded inline in the calling routine. !>\param[in] tg real, guess dewpoint temperature in Kelvin !>\param[in] pv real, vapor pressure in Pascals !\param[out] ftdpixg real, dewpoint temperature in Kelvin elemental function ftdpixg(tg,pv) !$$$ Subprogram Documentation Block ! ! Subprogram: ftdpixg Compute dewpoint temperature over ice ! Author: N Phillips w/NMC2X2 Date: 30 dec 82 ! ! Abstract: Exactly compute dewpoint temperature from vapor pressure. ! A guess dewpoint temperature must be provided. ! The water model assumes a perfect gas, constant specific heats ! for gas and ice, and neglects the volume of the ice. ! The model does account for the variation of the latent heat ! of sublimation with temperature. The liquid option is not included. ! The Clausius-Clapeyron equation is integrated from the triple point ! to get the formula ! pvs=con_psat*(tr**xa)*exp(xb*(1.-tr)) ! where tr is ttp/t and other values are physical constants. ! The formula is inverted by iterating Newtonian approximations ! for each pvs until t is found to within 1.e-6 Kelvin. ! This function can be expanded inline in the calling routine. ! ! Program History Log: ! 91-05-07 Iredell made into inlinable function ! 94-12-30 Iredell exact computation ! 1999-03-01 Iredell f90 module ! 2001-02-26 Iredell ice phase ! ! Usage: tdpi=ftdpixg(tg,pv) ! ! Input argument list: ! tg Real(krealfp) guess dewpoint temperature in Kelvin ! pv Real(krealfp) vapor pressure in Pascals ! ! Output argument list: ! ftdpixg Real(krealfp) dewpoint temperature in Kelvin ! ! Attributes: ! Language: Fortran 90. ! !$$$ implicit none real(krealfp) ftdpixg real(krealfp),intent(in):: tg,pv real(krealfp),parameter:: terrm=1.e-6 real(krealfp),parameter:: dldt=con_cvap-con_csol real(krealfp),parameter:: heat=con_hvap+con_hfus real(krealfp),parameter:: xpona=-dldt/con_rv real(krealfp),parameter:: xponb=-dldt/con_rv+heat/(con_rv*con_ttp) real(krealfp) t,tr,pvt,el,dpvt,terr integer i ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - t=tg do i=1,100 tr=con_ttp/t pvt=con_psat*(tr**xpona)*exp(xponb*(1.-tr)) el=heat+dldt*(t-con_ttp) dpvt=el*pvt/(con_rv*t**2) terr=(pvt-pv)/dpvt t=t-terr if(abs(terr).le.terrm) exit enddo ftdpixg=t ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end function !------------------------------------------------------------------------------- !> This subroutine computes dewpoint temperature table as a function of !! vapor pressure for inlinable function ftdp(). !! Exact dewpoint temperatures are calculated in subprogram ftdpxg(). !! The current implementation computes a table with a length !! of 5001 for vapor pressures ranging from 0.5 to 1000.5 Pascals !! giving a dewpoint temperature range of 208 to 319 Kelvin. subroutine gtdp !$$$ Subprogram Documentation Block ! ! Subprogram: gtdp Compute dewpoint temperature table ! Author: N Phillips w/NMC2X2 Date: 30 dec 82 ! ! Abstract: Compute dewpoint temperature table as a function of ! vapor pressure for inlinable function ftdp. ! Exact dewpoint temperatures are calculated in subprogram ftdpxg. ! The current implementation computes a table with a length ! of 5001 for vapor pressures ranging from 0.5 to 1000.5 Pascals ! giving a dewpoint temperature range of 208 to 319 Kelvin. ! ! Program History Log: ! 91-05-07 Iredell ! 94-12-30 Iredell expand table ! 1999-03-01 Iredell f90 module ! 2001-02-26 Iredell ice phase ! ! Usage: call gtdp ! ! Subprograms called: ! (ftdpxg) inlinable function to compute dewpoint temperature ! ! Attributes: ! Language: Fortran 90. ! !$$$ implicit none integer jx real(krealfp) xmin,xmax,xinc,t,x,pv ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - xmin=0.5 xmax=10000.5 xinc=(xmax-xmin)/(nxtdp-1) c1xtdp=1.-xmin/xinc c2xtdp=1./xinc t=208.0 do jx=1,nxtdp x=xmin+(jx-1)*xinc pv=x t=ftdpxg(t,pv) tbtdp(jx)=t enddo ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end subroutine !------------------------------------------------------------------------------- !> This function computes dewpoint temperature from vapor pressure. !! A linear interpolation is done between values in a lookup table !! computed in gtdp(). See documentation for ftdpxg() for details. !! Input values outside table range are reset to table extrema. !>\param[in] pv real, vapor pressure in Pascals !\param[out] ftdp real, dewpoint temperature in Kelvin elemental function ftdp(pv) !$$$ Subprogram Documentation Block ! ! Subprogram: ftdp Compute dewpoint temperature ! Author: N Phillips w/NMC2X2 Date: 30 dec 82 ! ! Abstract: Compute dewpoint temperature from vapor pressure. ! A linear interpolation is done between values in a lookup table ! computed in gtdp. See documentation for ftdpxg for details. ! Input values outside table range are reset to table extrema. ! The interpolation accuracy is better than 0.0005 Kelvin ! for dewpoint temperatures greater than 250 Kelvin, ! but decreases to 0.02 Kelvin for a dewpoint around 230 Kelvin. ! On the Cray, ftdp is about 75 times faster than exact calculation. ! This function should be expanded inline in the calling routine. ! ! Program History Log: ! 91-05-07 Iredell made into inlinable function ! 94-12-30 Iredell expand table ! 1999-03-01 Iredell f90 module ! 2001-02-26 Iredell ice phase ! ! Usage: tdp=ftdp(pv) ! ! Input argument list: ! pv Real(krealfp) vapor pressure in Pascals ! ! Output argument list: ! ftdp Real(krealfp) dewpoint temperature in Kelvin ! ! Attributes: ! Language: Fortran 90. ! !$$$ implicit none real(krealfp) ftdp real(krealfp),intent(in):: pv integer jx real(krealfp) xj ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - xj=min(max(c1xtdp+c2xtdp*pv,1._krealfp),real(nxtdp,krealfp)) jx=min(xj,nxtdp-1._krealfp) ftdp=tbtdp(jx)+(xj-jx)*(tbtdp(jx+1)-tbtdp(jx)) ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end function !------------------------------------------------------------------------------- !> This function computes dewpoint temperature from vapor pressure. !! A quadratic interpolation is done between values in a lookup table !! computed in gtdp(). See documentation for ftdpxg() for details. !! Input values outside table range are reset to table extrema. !>\param[in] pv real, vapor pressure in Pascals !\param[out] ftdpq real, dewpoint temperature in Kelvin elemental function ftdpq(pv) !$$$ Subprogram Documentation Block ! ! Subprogram: ftdpq Compute dewpoint temperature ! Author: N Phillips w/NMC2X2 Date: 30 dec 82 ! ! Abstract: Compute dewpoint temperature from vapor pressure. ! A quadratic interpolation is done between values in a lookup table ! computed in gtdp. see documentation for ftdpxg for details. ! Input values outside table range are reset to table extrema. ! the interpolation accuracy is better than 0.00001 Kelvin ! for dewpoint temperatures greater than 250 Kelvin, ! but decreases to 0.002 Kelvin for a dewpoint around 230 Kelvin. ! On the Cray, ftdpq is about 60 times faster than exact calculation. ! This function should be expanded inline in the calling routine. ! ! Program History Log: ! 91-05-07 Iredell made into inlinable function ! 94-12-30 Iredell quadratic interpolation ! 1999-03-01 Iredell f90 module ! 2001-02-26 Iredell ice phase ! ! Usage: tdp=ftdpq(pv) ! ! Input argument list: ! pv Real(krealfp) vapor pressure in Pascals ! ! Output argument list: ! ftdpq Real(krealfp) dewpoint temperature in Kelvin ! ! Attributes: ! Language: Fortran 90. ! !$$$ implicit none real(krealfp) ftdpq real(krealfp),intent(in):: pv integer jx real(krealfp) xj,dxj,fj1,fj2,fj3 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - xj=min(max(c1xtdp+c2xtdp*pv,1._krealfp),real(nxtdp,krealfp)) jx=min(max(nint(xj),2),nxtdp-1) dxj=xj-jx fj1=tbtdp(jx-1) fj2=tbtdp(jx) fj3=tbtdp(jx+1) ftdpq=(((fj3+fj1)/2-fj2)*dxj+(fj3-fj1)/2)*dxj+fj2 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end function !------------------------------------------------------------------------------- !> This function exactly computes dewpoint temperature from vapor pressure. !! An approximate dewpoint temperature for function ftdpxg() !! is obtained using ftdp() so gtdp() must be already called. !! See documentation for ftdpxg() for details. !>\param[in] pv real, vapor pressure in Pascals !\param[out] ftdpx real, dewpoint temperature in Kelvin elemental function ftdpx(pv) !$$$ Subprogram Documentation Block ! ! Subprogram: ftdpx Compute dewpoint temperature ! Author: N Phillips w/NMC2X2 Date: 30 dec 82 ! ! Abstract: exactly compute dewpoint temperature from vapor pressure. ! An approximate dewpoint temperature for function ftdpxg ! is obtained using ftdp so gtdp must be already called. ! See documentation for ftdpxg for details. ! ! Program History Log: ! 91-05-07 Iredell made into inlinable function ! 94-12-30 Iredell exact computation ! 1999-03-01 Iredell f90 module ! 2001-02-26 Iredell ice phase ! ! Usage: tdp=ftdpx(pv) ! ! Input argument list: ! pv Real(krealfp) vapor pressure in Pascals ! ! Output argument list: ! ftdpx Real(krealfp) dewpoint temperature in Kelvin ! ! Subprograms called: ! (ftdp) inlinable function to compute dewpoint temperature ! (ftdpxg) inlinable function to compute dewpoint temperature ! ! Attributes: ! Language: Fortran 90. ! !$$$ implicit none real(krealfp) ftdpx real(krealfp),intent(in):: pv real(krealfp) tg ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - tg=ftdp(pv) ftdpx=ftdpxg(tg,pv) ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end function !------------------------------------------------------------------------------- !> This function exactly computes dewpoint temperature from vapor pressure. !! A guess dewpoint temperature must be provided. !! The saturation vapor pressure over either liquid and ice is computed !! over liquid for temperatures above the triple point, !! over ice for temperatures 20 degress below the triple point, !! and a linear combination of the two for temperatures in between. !! The water model assumes a perfect gas, constant specific heats !! for gas, liquid and ice, and neglects the volume of the condensate. !! The model does account for the variation of the latent heat !! of condensation and sublimation with temperature. !! The Clausius-Clapeyron equation is integrated from the triple point !! to get the formula: !!\n pvsl=con_psat*(tr**xa)*exp(xb*(1.-tr)) !!\n where tr is ttp/t and other values are physical constants. !!\n The reference for this decision is Emanuel(1994), pages 116-117. !! The formula is inverted by iterating Newtonian approximations !! for each pvs until t is found to within 1.e-6 Kelvin. !! This function can be expanded inline in the calling routine. !>\param[in] tg real, guess dewpoint temperature in Kelvin !>\param[in] pv real, vapor pressure in Pascals !\param[out] ftdpxg real, dewpoint temperature in Kelvin elemental function ftdpxg(tg,pv) !$$$ Subprogram Documentation Block ! ! Subprogram: ftdpxg Compute dewpoint temperature ! Author: N Phillips w/NMC2X2 Date: 30 dec 82 ! ! Abstract: Exactly compute dewpoint temperature from vapor pressure. ! A guess dewpoint temperature must be provided. ! The saturation vapor pressure over either liquid and ice is computed ! over liquid for temperatures above the triple point, ! over ice for temperatures 20 degress below the triple point, ! and a linear combination of the two for temperatures in between. ! The water model assumes a perfect gas, constant specific heats ! for gas, liquid and ice, and neglects the volume of the condensate. ! The model does account for the variation of the latent heat ! of condensation and sublimation with temperature. ! The Clausius-Clapeyron equation is integrated from the triple point ! to get the formula ! pvsl=con_psat*(tr**xa)*exp(xb*(1.-tr)) ! where tr is ttp/t and other values are physical constants. ! The reference for this decision is Emanuel(1994), pages 116-117. ! The formula is inverted by iterating Newtonian approximations ! for each pvs until t is found to within 1.e-6 Kelvin. ! This function can be expanded inline in the calling routine. ! ! Program History Log: ! 91-05-07 Iredell made into inlinable function ! 94-12-30 Iredell exact computation ! 1999-03-01 Iredell f90 module ! 2001-02-26 Iredell ice phase ! ! Usage: tdp=ftdpxg(tg,pv) ! ! Input argument list: ! tg Real(krealfp) guess dewpoint temperature in Kelvin ! pv Real(krealfp) vapor pressure in Pascals ! ! Output argument list: ! ftdpxg Real(krealfp) dewpoint temperature in Kelvin ! ! Attributes: ! Language: Fortran 90. ! !$$$ implicit none real(krealfp) ftdpxg real(krealfp),intent(in):: tg,pv real(krealfp),parameter:: terrm=1.e-6 real(krealfp),parameter:: tliq=con_ttp real(krealfp),parameter:: tice=con_ttp-20.0 real(krealfp),parameter:: dldtl=con_cvap-con_cliq real(krealfp),parameter:: heatl=con_hvap real(krealfp),parameter:: xponal=-dldtl/con_rv real(krealfp),parameter:: xponbl=-dldtl/con_rv+heatl/(con_rv*con_ttp) real(krealfp),parameter:: dldti=con_cvap-con_csol real(krealfp),parameter:: heati=con_hvap+con_hfus real(krealfp),parameter:: xponai=-dldti/con_rv real(krealfp),parameter:: xponbi=-dldti/con_rv+heati/(con_rv*con_ttp) real(krealfp) t,tr,w,pvtl,pvti,pvt,ell,eli,el,dpvt,terr integer i ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - t=tg do i=1,100 tr=con_ttp/t if(t.ge.tliq) then pvt=con_psat*(tr**xponal)*exp(xponbl*(1.-tr)) el=heatl+dldtl*(t-con_ttp) dpvt=el*pvt/(con_rv*t**2) elseif(t.lt.tice) then pvt=con_psat*(tr**xponai)*exp(xponbi*(1.-tr)) el=heati+dldti*(t-con_ttp) dpvt=el*pvt/(con_rv*t**2) else w=(t-tice)/(tliq-tice) pvtl=con_psat*(tr**xponal)*exp(xponbl*(1.-tr)) pvti=con_psat*(tr**xponai)*exp(xponbi*(1.-tr)) pvt=w*pvtl+(1.-w)*pvti ell=heatl+dldtl*(t-con_ttp) eli=heati+dldti*(t-con_ttp) dpvt=(w*ell*pvtl+(1.-w)*eli*pvti)/(con_rv*t**2) endif terr=(pvt-pv)/dpvt t=t-terr if(abs(terr).le.terrm) exit enddo ftdpxg=t ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end function !------------------------------------------------------------------------------- !> This subroutine computes equivalent potential temperature table !! as a function of LCL temperature and pressure over 1e5 Pa !! to the kappa power for function fthe(). !! Equivalent potential temperatures are calculated in subprogram fthex() !! the current implementation computes a table with a first dimension !! of 241 for temperatures ranging from 183.16 to 303.16 Kelvin !! and a second dimension of 151 for pressure over 1e5 Pa !! to the kappa power ranging from \f$0.04^{rocp}\f$ to \f$1.10^{rocp}\f$. subroutine gthe !$$$ Subprogram Documentation Block ! ! Subprogram: gthe Compute equivalent potential temperature table ! Author: N Phillips w/NMC2X2 Date: 30 dec 82 ! ! Abstract: Compute equivalent potential temperature table ! as a function of LCL temperature and pressure over 1e5 Pa ! to the kappa power for function fthe. ! Equivalent potential temperatures are calculated in subprogram fthex ! the current implementation computes a table with a first dimension ! of 241 for temperatures ranging from 183.16 to 303.16 Kelvin ! and a second dimension of 151 for pressure over 1e5 Pa ! to the kappa power ranging from 0.04**rocp to 1.10**rocp. ! ! Program History Log: ! 91-05-07 Iredell ! 94-12-30 Iredell expand table ! 1999-03-01 Iredell f90 module ! ! Usage: call gthe ! ! Subprograms called: ! (fthex) inlinable function to compute equiv. pot. temperature ! ! Attributes: ! Language: Fortran 90. ! !$$$ implicit none integer jx,jy real(krealfp) xmin,xmax,ymin,ymax,xinc,yinc,x,y,pk,t ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - xmin=con_ttp-90._krealfp xmax=con_ttp+30._krealfp ymin=0.04_krealfp**con_rocp ymax=1.10_krealfp**con_rocp xinc=(xmax-xmin)/(nxthe-1) c1xthe=1.-xmin/xinc c2xthe=1./xinc yinc=(ymax-ymin)/(nythe-1) c1ythe=1.-ymin/yinc c2ythe=1./yinc do jy=1,nythe y=ymin+(jy-1)*yinc pk=y do jx=1,nxthe x=xmin+(jx-1)*xinc t=x tbthe(jx,jy)=fthex(t,pk) enddo enddo ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end subroutine !------------------------------------------------------------------------------- !> This function computes equivalent potential temperature at the LCL !! from temperature and pressure over 1e5 Pa to the kappa power. !! A bilinear interpolation is done between values in a lookup table !! computed in gthe(). see documentation for fthex() for details. !! Input values outside table range are reset to table extrema, !! except zero is returned for too cold or high LCLs. !>\param[in] t real, LCL temperature in Kelvin !>\param[in] pk real, LCL pressure over 1e5 Pa to the kappa power !\param[out] fthe real, equivalent potential temperature in Kelvin elemental function fthe(t,pk) !$$$ Subprogram Documentation Block ! ! Subprogram: fthe Compute equivalent potential temperature ! Author: N Phillips w/NMC2X2 Date: 30 dec 82 ! ! Abstract: Compute equivalent potential temperature at the LCL ! from temperature and pressure over 1e5 Pa to the kappa power. ! A bilinear interpolation is done between values in a lookup table ! computed in gthe. see documentation for fthex for details. ! Input values outside table range are reset to table extrema, ! except zero is returned for too cold or high LCLs. ! The interpolation accuracy is better than 0.01 Kelvin. ! On the Cray, fthe is almost 6 times faster than exact calculation. ! This function should be expanded inline in the calling routine. ! ! Program History Log: ! 91-05-07 Iredell made into inlinable function ! 94-12-30 Iredell expand table ! 1999-03-01 Iredell f90 module ! ! Usage: the=fthe(t,pk) ! ! Input argument list: ! t Real(krealfp) LCL temperature in Kelvin ! pk Real(krealfp) LCL pressure over 1e5 Pa to the kappa power ! ! Output argument list: ! fthe Real(krealfp) equivalent potential temperature in Kelvin ! ! Attributes: ! Language: Fortran 90. ! !$$$ implicit none real(krealfp) fthe real(krealfp),intent(in):: t,pk integer jx,jy real(krealfp) xj,yj,ftx1,ftx2 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - xj=min(c1xthe+c2xthe*t,real(nxthe,krealfp)) yj=min(c1ythe+c2ythe*pk,real(nythe,krealfp)) if(xj.ge.1..and.yj.ge.1.) then jx=min(xj,nxthe-1._krealfp) jy=min(yj,nythe-1._krealfp) ftx1=tbthe(jx,jy)+(xj-jx)*(tbthe(jx+1,jy)-tbthe(jx,jy)) ftx2=tbthe(jx,jy+1)+(xj-jx)*(tbthe(jx+1,jy+1)-tbthe(jx,jy+1)) fthe=ftx1+(yj-jy)*(ftx2-ftx1) else fthe=0. endif ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end function !------------------------------------------------------------------------------- !> This function computes equivalent potential temperature at the LCL !! from temperature and pressure over 1e5 Pa to the kappa power. !! A biquadratic interpolation is done between values in a lookup table !! computed in gthe(). see documentation for fthex() for details. !! Input values outside table range are reset to table extrema, !! except zero is returned for too cold or high LCLs. !>\param[in] t real, LCL temperature in Kelvin !>\param[in] pk real, LCL pressure over 1e5 Pa to the kappa power !\param[out] ftheq real, equivalent potential temperature in Kelvin elemental function ftheq(t,pk) !$$$ Subprogram Documentation Block ! ! Subprogram: ftheq Compute equivalent potential temperature ! Author: N Phillips w/NMC2X2 Date: 30 dec 82 ! ! Abstract: Compute equivalent potential temperature at the LCL ! from temperature and pressure over 1e5 Pa to the kappa power. ! A biquadratic interpolation is done between values in a lookup table ! computed in gthe. see documentation for fthex for details. ! Input values outside table range are reset to table extrema, ! except zero is returned for too cold or high LCLs. ! The interpolation accuracy is better than 0.0002 Kelvin. ! On the Cray, ftheq is almost 3 times faster than exact calculation. ! This function should be expanded inline in the calling routine. ! ! Program History Log: ! 91-05-07 Iredell made into inlinable function ! 94-12-30 Iredell quadratic interpolation ! 1999-03-01 Iredell f90 module ! ! Usage: the=ftheq(t,pk) ! ! Input argument list: ! t Real(krealfp) LCL temperature in Kelvin ! pk Real(krealfp) LCL pressure over 1e5 Pa to the kappa power ! ! Output argument list: ! ftheq Real(krealfp) equivalent potential temperature in Kelvin ! ! Attributes: ! Language: Fortran 90. ! !$$$ implicit none real(krealfp) ftheq real(krealfp),intent(in):: t,pk integer jx,jy real(krealfp) xj,yj,dxj,dyj real(krealfp) ft11,ft12,ft13,ft21,ft22,ft23,ft31,ft32,ft33 real(krealfp) ftx1,ftx2,ftx3 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - xj=min(c1xthe+c2xthe*t,real(nxthe,krealfp)) yj=min(c1ythe+c2ythe*pk,real(nythe,krealfp)) if(xj.ge.1..and.yj.ge.1.) then jx=min(max(nint(xj),2),nxthe-1) jy=min(max(nint(yj),2),nythe-1) dxj=xj-jx dyj=yj-jy ft11=tbthe(jx-1,jy-1) ft12=tbthe(jx-1,jy) ft13=tbthe(jx-1,jy+1) ft21=tbthe(jx,jy-1) ft22=tbthe(jx,jy) ft23=tbthe(jx,jy+1) ft31=tbthe(jx+1,jy-1) ft32=tbthe(jx+1,jy) ft33=tbthe(jx+1,jy+1) ftx1=(((ft31+ft11)/2-ft21)*dxj+(ft31-ft11)/2)*dxj+ft21 ftx2=(((ft32+ft12)/2-ft22)*dxj+(ft32-ft12)/2)*dxj+ft22 ftx3=(((ft33+ft13)/2-ft23)*dxj+(ft33-ft13)/2)*dxj+ft23 ftheq=(((ftx3+ftx1)/2-ftx2)*dyj+(ftx3-ftx1)/2)*dyj+ftx2 else ftheq=0. endif ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end function !------------------------------------------------------------------------------- ! elemental function fthex(t,pk) !> This function exactly computes equivalent potential temperature at the LCL !! from temperature and pressure over 1e5 Pa to the kappa power. !! Equivalent potential temperature is constant for a saturated parcel !! rising adiabatically up a moist adiabat when the heat and mass !! of the condensed water are neglected. Ice is also neglected. !! The formula for equivalent potential temperature (Holton) is !!\n the=t*(pd**(-rocp))*exp(el*eps*pv/(cp*t*pd)) !!\n where t is the temperature, pv is the saturated vapor pressure, !! pd is the dry pressure p-pv, el is the temperature dependent !! latent heat of condensation hvap+dldt*(t-ttp), and other values !! are physical constants defined in parameter statements in the code. !! Zero is returned if the input values make saturation impossible. !! This function should be expanded inline in the calling routine. !>\param[in] t real, LCL temperature in Kelvin !>\param[in] pk real, LCL pressure over 1e5 Pa to the kappa power !\param[out] fthex real, equivalent potential temperature in Kelvin function fthex(t,pk) !$$$ Subprogram Documentation Block ! ! Subprogram: fthex Compute equivalent potential temperature ! Author: N Phillips w/NMC2X2 Date: 30 dec 82 ! ! Abstract: Exactly compute equivalent potential temperature at the LCL ! from temperature and pressure over 1e5 Pa to the kappa power. ! Equivalent potential temperature is constant for a saturated parcel ! rising adiabatically up a moist adiabat when the heat and mass ! of the condensed water are neglected. Ice is also neglected. ! The formula for equivalent potential temperature (Holton) is ! the=t*(pd**(-rocp))*exp(el*eps*pv/(cp*t*pd)) ! where t is the temperature, pv is the saturated vapor pressure, ! pd is the dry pressure p-pv, el is the temperature dependent ! latent heat of condensation hvap+dldt*(t-ttp), and other values ! are physical constants defined in parameter statements in the code. ! Zero is returned if the input values make saturation impossible. ! This function should be expanded inline in the calling routine. ! ! Program History Log: ! 91-05-07 Iredell made into inlinable function ! 94-12-30 Iredell exact computation ! 1999-03-01 Iredell f90 module ! ! Usage: the=fthex(t,pk) ! ! Input argument list: ! t Real(krealfp) LCL temperature in Kelvin ! pk Real(krealfp) LCL pressure over 1e5 Pa to the kappa power ! ! Output argument list: ! fthex Real(krealfp) equivalent potential temperature in Kelvin ! ! Attributes: ! Language: Fortran 90. ! !$$$ implicit none real(krealfp) fthex real(krealfp),intent(in):: t,pk real(krealfp) p,tr,pv,pd,el,expo,expmax ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - p=pk**con_cpor tr=con_ttp/t pv=psatb*(tr**con_xpona)*exp(con_xponb*(1.-tr)) pd=p-pv if(pd.gt.pv) then el=con_hvap+con_dldt*(t-con_ttp) expo=el*con_eps*pv/(con_cp*t*pd) fthex=t*pd**(-con_rocp)*exp(expo) else fthex=0. endif ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end function !------------------------------------------------------------------------------- !> This subroutine computes temperature and specific humidity tables !! as a function of equivalent potential temperature and !! pressure over 1e5 Pa to the kappa power for subprogram stma(). !! Exact parcel temperatures are calculated in subprogram stmaxg(). !! The current implementation computes a table with a first dimension !! of 151 for equivalent potential temperatures ranging from 200 to 500 !! Kelvin and a second dimension of 121 for pressure over 1e5 Pa !! to the kappa power ranging from \f$0.01^{rocp}\f$ to \f$1.10^{rocp}\f$. subroutine gtma !$$$ Subprogram Documentation Block ! ! Subprogram: gtma Compute moist adiabat tables ! Author: N Phillips w/NMC2X2 Date: 30 dec 82 ! ! Abstract: Compute temperature and specific humidity tables ! as a function of equivalent potential temperature and ! pressure over 1e5 Pa to the kappa power for subprogram stma. ! Exact parcel temperatures are calculated in subprogram stmaxg. ! The current implementation computes a table with a first dimension ! of 151 for equivalent potential temperatures ranging from 200 to 500 ! Kelvin and a second dimension of 121 for pressure over 1e5 Pa ! to the kappa power ranging from 0.01**rocp to 1.10**rocp. ! ! Program History Log: ! 91-05-07 Iredell ! 94-12-30 Iredell expand table ! 1999-03-01 Iredell f90 module ! ! Usage: call gtma ! ! Subprograms called: ! (stmaxg) inlinable subprogram to compute parcel temperature ! ! Attributes: ! Language: Fortran 90. ! !$$$ implicit none integer jx,jy real(krealfp) xmin,xmax,ymin,ymax,xinc,yinc,x,y,pk,the,t,q,tg ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - xmin=200._krealfp xmax=500._krealfp ymin=0.01_krealfp**con_rocp ymax=1.10_krealfp**con_rocp xinc=(xmax-xmin)/(nxma-1) c1xma=1.-xmin/xinc c2xma=1./xinc yinc=(ymax-ymin)/(nyma-1) c1yma=1.-ymin/yinc c2yma=1./yinc do jy=1,nyma y=ymin+(jy-1)*yinc pk=y tg=xmin*y do jx=1,nxma x=xmin+(jx-1)*xinc the=x call stmaxg(tg,the,pk,t,q) tbtma(jx,jy)=t tbqma(jx,jy)=q tg=t enddo enddo ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end subroutine !------------------------------------------------------------------------------- !> This subroutine computes temperature and specific humidity of a parcel !! lifted up a moist adiabat from equivalent potential temperature !! at the LCL and pressure over 1e5 Pa to the kappa power. !! Bilinear interpolations are done between values in a lookup table !! computed in gtma(). See documentation for stmaxg() for details. !! Input values outside table range are reset to table extrema. !>\param[in] the real, equivalent potential temperature in Kelvin !>\param[in] pk real, pressure over 1e5 Pa to the kappa power !>\param[out] tma real, parcel temperature in Kelvin !>\param[out] qma real, parcel specific humidity in kg/kg elemental subroutine stma(the,pk,tma,qma) !$$$ Subprogram Documentation Block ! ! Subprogram: stma Compute moist adiabat temperature ! Author: N Phillips w/NMC2X2 Date: 30 dec 82 ! ! Abstract: Compute temperature and specific humidity of a parcel ! lifted up a moist adiabat from equivalent potential temperature ! at the LCL and pressure over 1e5 Pa to the kappa power. ! Bilinear interpolations are done between values in a lookup table ! computed in gtma. See documentation for stmaxg for details. ! Input values outside table range are reset to table extrema. ! The interpolation accuracy is better than 0.01 Kelvin ! and 5.e-6 kg/kg for temperature and humidity, respectively. ! On the Cray, stma is about 35 times faster than exact calculation. ! This subprogram should be expanded inline in the calling routine. ! ! Program History Log: ! 91-05-07 Iredell made into inlinable function ! 94-12-30 Iredell expand table ! 1999-03-01 Iredell f90 module ! ! Usage: call stma(the,pk,tma,qma) ! ! Input argument list: ! the Real(krealfp) equivalent potential temperature in Kelvin ! pk Real(krealfp) pressure over 1e5 Pa to the kappa power ! ! Output argument list: ! tma Real(krealfp) parcel temperature in Kelvin ! qma Real(krealfp) parcel specific humidity in kg/kg ! ! Attributes: ! Language: Fortran 90. ! !$$$ implicit none real(krealfp),intent(in):: the,pk real(krealfp),intent(out):: tma,qma integer jx,jy real(krealfp) xj,yj,ftx1,ftx2,qx1,qx2 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - xj=min(max(c1xma+c2xma*the,1._krealfp),real(nxma,krealfp)) yj=min(max(c1yma+c2yma*pk,1._krealfp),real(nyma,krealfp)) jx=min(xj,nxma-1._krealfp) jy=min(yj,nyma-1._krealfp) ftx1=tbtma(jx,jy)+(xj-jx)*(tbtma(jx+1,jy)-tbtma(jx,jy)) ftx2=tbtma(jx,jy+1)+(xj-jx)*(tbtma(jx+1,jy+1)-tbtma(jx,jy+1)) tma=ftx1+(yj-jy)*(ftx2-ftx1) qx1=tbqma(jx,jy)+(xj-jx)*(tbqma(jx+1,jy)-tbqma(jx,jy)) qx2=tbqma(jx,jy+1)+(xj-jx)*(tbqma(jx+1,jy+1)-tbqma(jx,jy+1)) qma=qx1+(yj-jy)*(qx2-qx1) ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end subroutine !------------------------------------------------------------------------------- !> This subroutine computes temperature and specific humidity of a parcel !! lifted up a moist adiabat from equivalent potential temperature !! at the LCL and pressure over 1e5 Pa to the kappa power. !! Biquadratic interpolations are done between values in a lookup table !! computed in gtma(). See documentation for stmaxg() for details. !! Input values outside table range are reset to table extrema. !>\param[in] the real, equivalent potential temperature in Kelvin !>\param[in] pk real, pressure over 1e5 Pa to the kappa power !>\param[out] tma real, parcel temperature in Kelvin !>\param[out] qma real, parcel specific humidity in kg/kg elemental subroutine stmaq(the,pk,tma,qma) !$$$ Subprogram Documentation Block ! ! Subprogram: stmaq Compute moist adiabat temperature ! Author: N Phillips w/NMC2X2 Date: 30 dec 82 ! ! Abstract: Compute temperature and specific humidity of a parcel ! lifted up a moist adiabat from equivalent potential temperature ! at the LCL and pressure over 1e5 Pa to the kappa power. ! Biquadratic interpolations are done between values in a lookup table ! computed in gtma. See documentation for stmaxg for details. ! Input values outside table range are reset to table extrema. ! the interpolation accuracy is better than 0.0005 Kelvin ! and 1.e-7 kg/kg for temperature and humidity, respectively. ! On the Cray, stmaq is about 25 times faster than exact calculation. ! This subprogram should be expanded inline in the calling routine. ! ! Program History Log: ! 91-05-07 Iredell made into inlinable function ! 94-12-30 Iredell quadratic interpolation ! 1999-03-01 Iredell f90 module ! ! Usage: call stmaq(the,pk,tma,qma) ! ! Input argument list: ! the Real(krealfp) equivalent potential temperature in Kelvin ! pk Real(krealfp) pressure over 1e5 Pa to the kappa power ! ! Output argument list: ! tmaq Real(krealfp) parcel temperature in Kelvin ! qma Real(krealfp) parcel specific humidity in kg/kg ! ! Attributes: ! Language: Fortran 90. ! !$$$ implicit none real(krealfp),intent(in):: the,pk real(krealfp),intent(out):: tma,qma integer jx,jy real(krealfp) xj,yj,dxj,dyj real(krealfp) ft11,ft12,ft13,ft21,ft22,ft23,ft31,ft32,ft33 real(krealfp) ftx1,ftx2,ftx3 real(krealfp) q11,q12,q13,q21,q22,q23,q31,q32,q33,qx1,qx2,qx3 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - xj=min(max(c1xma+c2xma*the,1._krealfp),real(nxma,krealfp)) yj=min(max(c1yma+c2yma*pk,1._krealfp),real(nyma,krealfp)) jx=min(max(nint(xj),2),nxma-1) jy=min(max(nint(yj),2),nyma-1) dxj=xj-jx dyj=yj-jy ft11=tbtma(jx-1,jy-1) ft12=tbtma(jx-1,jy) ft13=tbtma(jx-1,jy+1) ft21=tbtma(jx,jy-1) ft22=tbtma(jx,jy) ft23=tbtma(jx,jy+1) ft31=tbtma(jx+1,jy-1) ft32=tbtma(jx+1,jy) ft33=tbtma(jx+1,jy+1) ftx1=(((ft31+ft11)/2-ft21)*dxj+(ft31-ft11)/2)*dxj+ft21 ftx2=(((ft32+ft12)/2-ft22)*dxj+(ft32-ft12)/2)*dxj+ft22 ftx3=(((ft33+ft13)/2-ft23)*dxj+(ft33-ft13)/2)*dxj+ft23 tma=(((ftx3+ftx1)/2-ftx2)*dyj+(ftx3-ftx1)/2)*dyj+ftx2 q11=tbqma(jx-1,jy-1) q12=tbqma(jx-1,jy) q13=tbqma(jx-1,jy+1) q21=tbqma(jx,jy-1) q22=tbqma(jx,jy) q23=tbqma(jx,jy+1) q31=tbqma(jx+1,jy-1) q32=tbqma(jx+1,jy) q33=tbqma(jx+1,jy+1) qx1=(((q31+q11)/2-q21)*dxj+(q31-q11)/2)*dxj+q21 qx2=(((q32+q12)/2-q22)*dxj+(q32-q12)/2)*dxj+q22 qx3=(((q33+q13)/2-q23)*dxj+(q33-q13)/2)*dxj+q23 qma=(((qx3+qx1)/2-qx2)*dyj+(qx3-qx1)/2)*dyj+qx2 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end subroutine !------------------------------------------------------------------------------- !> This subroutine exactly computes temperature and humidity of a parcel !! lifted up a moist adiabat from equivalent potential temperature !! at the LCL and pressure over 1e5 Pa to the kappa power. !! An approximate parcel temperature for subprogram stmaxg() !! is obtained using stma() so gtma() must be already called. !! See documentation for stmaxg() for details. !>\param[in] the real, equivalent potential temperature in Kelvin !>\param[in] pk real, pressure over 1e5 Pa to the kappa power !>\param[out] tma real, parcel temperature in Kelvin !>\param[out] qma real, parcel specific humidity in kg/kg subroutine stmax(the,pk,tma,qma) !$$$ Subprogram Documentation Block ! ! Subprogram: stmax Compute moist adiabat temperature ! Author: N Phillips w/NMC2X2 Date: 30 dec 82 ! ! Abstract: Exactly compute temperature and humidity of a parcel ! lifted up a moist adiabat from equivalent potential temperature ! at the LCL and pressure over 1e5 Pa to the kappa power. ! An approximate parcel temperature for subprogram stmaxg ! is obtained using stma so gtma must be already called. ! See documentation for stmaxg for details. ! ! Program History Log: ! 91-05-07 Iredell made into inlinable function ! 94-12-30 Iredell exact computation ! 1999-03-01 Iredell f90 module ! ! Usage: call stmax(the,pk,tma,qma) ! ! Input argument list: ! the Real(krealfp) equivalent potential temperature in Kelvin ! pk Real(krealfp) pressure over 1e5 Pa to the kappa power ! ! Output argument list: ! tma Real(krealfp) parcel temperature in Kelvin ! qma Real(krealfp) parcel specific humidity in kg/kg ! ! Subprograms called: ! (stma) inlinable subprogram to compute parcel temperature ! (stmaxg) inlinable subprogram to compute parcel temperature ! ! Attributes: ! Language: Fortran 90. ! !$$$ implicit none real(krealfp),intent(in):: the,pk real(krealfp),intent(out):: tma,qma real(krealfp) tg,qg ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - call stma(the,pk,tg,qg) call stmaxg(tg,the,pk,tma,qma) ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end subroutine !------------------------------------------------------------------------------- !> This subroutine exactly computes temperature and humidity of a parcel !! lifted up a moist adiabat from equivalent potential temperature !! at the LCL and pressure over 1e5 Pa to the kappa power. !! A guess parcel temperature must be provided. !! Equivalent potential temperature is constant for a saturated parcel !! rising adiabatically up a moist adiabat when the heat and mass !! of the condensed water are neglected. Ice is also neglected. !! The formula for equivalent potential temperature (Holton) is !!\n the=t*(pd**(-rocp))*exp(el*eps*pv/(cp*t*pd)) !!\n where t is the temperature, pv is the saturated vapor pressure, !! pd is the dry pressure p-pv, el is the temperature dependent !! latent heat of condensation hvap+dldt*(t-ttp), and other values !! are physical constants defined in parameter statements in the code. !! The formula is inverted by iterating Newtonian approximations !! for each the and p until t is found to within 1.e-4 Kelvin. !! The specific humidity is then computed from pv and pd. !! This subprogram can be expanded inline in the calling routine. !>\param[in] tg real, guess parcel temperature in Kelvin !>\param[in] the real, equivalent potential temperature in Kelvin !>\param[in] pk real, pressure over 1e5 Pa to the kappa power !>\param[out] tma real, parcel temperature in Kelvin !>\param[out] qma real, parcel specific humidity in kg/kg subroutine stmaxg(tg,the,pk,tma,qma) !$$$ Subprogram Documentation Block ! ! Subprogram: stmaxg Compute moist adiabat temperature ! Author: N Phillips w/NMC2X2 Date: 30 dec 82 ! ! Abstract: exactly compute temperature and humidity of a parcel ! lifted up a moist adiabat from equivalent potential temperature ! at the LCL and pressure over 1e5 Pa to the kappa power. ! A guess parcel temperature must be provided. ! Equivalent potential temperature is constant for a saturated parcel ! rising adiabatically up a moist adiabat when the heat and mass ! of the condensed water are neglected. Ice is also neglected. ! The formula for equivalent potential temperature (Holton) is ! the=t*(pd**(-rocp))*exp(el*eps*pv/(cp*t*pd)) ! where t is the temperature, pv is the saturated vapor pressure, ! pd is the dry pressure p-pv, el is the temperature dependent ! latent heat of condensation hvap+dldt*(t-ttp), and other values ! are physical constants defined in parameter statements in the code. ! The formula is inverted by iterating Newtonian approximations ! for each the and p until t is found to within 1.e-4 Kelvin. ! The specific humidity is then computed from pv and pd. ! This subprogram can be expanded inline in the calling routine. ! ! Program History Log: ! 91-05-07 Iredell made into inlinable function ! 94-12-30 Iredell exact computation ! 1999-03-01 Iredell f90 module ! ! Usage: call stmaxg(tg,the,pk,tma,qma) ! ! Input argument list: ! tg Real(krealfp) guess parcel temperature in Kelvin ! the Real(krealfp) equivalent potential temperature in Kelvin ! pk Real(krealfp) pressure over 1e5 Pa to the kappa power ! ! Output argument list: ! tma Real(krealfp) parcel temperature in Kelvin ! qma Real(krealfp) parcel specific humidity in kg/kg ! ! Attributes: ! Language: Fortran 90. ! !$$$ implicit none real(krealfp),intent(in):: tg,the,pk real(krealfp),intent(out):: tma,qma real(krealfp),parameter:: terrm=1.e-4 real(krealfp) t,p,tr,pv,pd,el,expo,thet,dthet,terr integer i ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - t=tg p=pk**con_cpor do i=1,100 tr=con_ttp/t pv=psatb*(tr**con_xpona)*exp(con_xponb*(1.-tr)) pd=p-pv el=con_hvap+con_dldt*(t-con_ttp) expo=el*con_eps*pv/(con_cp*t*pd) thet=t*pd**(-con_rocp)*exp(expo) dthet=thet/t*(1.+expo*(con_dldt*t/el+el*p/(con_rv*t*pd))) terr=(thet-the)/dthet t=t-terr if(abs(terr).le.terrm) exit enddo tma=t tr=con_ttp/t pv=psatb*(tr**con_xpona)*exp(con_xponb*(1.-tr)) pd=p-pv qma=con_eps*pv/(pd+con_eps*pv) ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end subroutine !------------------------------------------------------------------------------- !> This subroutine computes pressure to the kappa table as a function of pressure !! for the table lookup function fpkap(). !! Exact pressure to the kappa values are calculated in subprogram fpkapx(). !! The current implementation computes a table with a length !! of 5501 for pressures ranging up to 110000 Pascals. subroutine gpkap !$$$ Subprogram documentation block ! ! Subprogram: gpkap Compute coefficients for p**kappa ! Author: Phillips org: w/NMC2X2 Date: 29 dec 82 ! ! Abstract: Computes pressure to the kappa table as a function of pressure ! for the table lookup function fpkap. ! Exact pressure to the kappa values are calculated in subprogram fpkapx. ! The current implementation computes a table with a length ! of 5501 for pressures ranging up to 110000 Pascals. ! ! Program History Log: ! 94-12-30 Iredell ! 1999-03-01 Iredell f90 module ! 1999-03-24 Iredell table lookup ! ! Usage: call gpkap ! ! Subprograms called: ! fpkapx function to compute exact pressure to the kappa ! ! Attributes: ! Language: Fortran 90. ! !$$$ implicit none integer jx real(krealfp) xmin,xmax,xinc,x,p ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - xmin=0._krealfp xmax=110000._krealfp xinc=(xmax-xmin)/(nxpkap-1) c1xpkap=1.-xmin/xinc c2xpkap=1./xinc do jx=1,nxpkap x=xmin+(jx-1)*xinc p=x tbpkap(jx)=fpkapx(p) enddo ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end subroutine !------------------------------------------------------------------------------- !> This subroutine raises pressure over 1e5 Pa to the kappa power. !! A linear interpolation is done between values in a lookup table !! computed in gpkap(). See documentation for fpkapx() for details. !! Input values outside table range are reset to table extrema. !>\param[in] p real, pressure in Pascals !\param[out] fpkap real, p over 1e5 Pa to the kappa power elemental function fpkap(p) !$$$ Subprogram Documentation Block ! ! Subprogram: fpkap raise pressure to the kappa power. ! Author: N Phillips w/NMC2X2 Date: 30 dec 82 ! ! Abstract: Raise pressure over 1e5 Pa to the kappa power. ! A linear interpolation is done between values in a lookup table ! computed in gpkap. See documentation for fpkapx for details. ! Input values outside table range are reset to table extrema. ! The interpolation accuracy ranges from 9 decimal places ! at 100000 Pascals to 5 decimal places at 1000 Pascals. ! On the Cray, fpkap is over 5 times faster than exact calculation. ! This function should be expanded inline in the calling routine. ! ! Program History Log: ! 91-05-07 Iredell made into inlinable function ! 94-12-30 Iredell standardized kappa, ! increased range and accuracy ! 1999-03-01 Iredell f90 module ! 1999-03-24 Iredell table lookup ! ! Usage: pkap=fpkap(p) ! ! Input argument list: ! p Real(krealfp) pressure in Pascals ! ! Output argument list: ! fpkap Real(krealfp) p over 1e5 Pa to the kappa power ! ! Attributes: ! Language: Fortran 90. ! !$$$ implicit none real(krealfp) fpkap real(krealfp),intent(in):: p integer jx real(krealfp) xj ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - xj=min(max(c1xpkap+c2xpkap*p,1._krealfp),real(nxpkap,krealfp)) jx=min(xj,nxpkap-1._krealfp) fpkap=tbpkap(jx)+(xj-jx)*(tbpkap(jx+1)-tbpkap(jx)) ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end function !------------------------------------------------------------------------------- !> This function raises pressure over 1e5 Pa to the kappa power. !! A quadratic interpolation is done between values in a lookup table !! computed in gpkap(). see documentation for fpkapx() for details. !! Input values outside table range are reset to table extrema. !>\param[in] p real, pressure in Pascals !\param[out] fpkapq real, p over 1e5 Pa to the kappa power elemental function fpkapq(p) !$$$ Subprogram Documentation Block ! ! Subprogram: fpkapq raise pressure to the kappa power. ! Author: N Phillips w/NMC2X2 Date: 30 dec 82 ! ! Abstract: Raise pressure over 1e5 Pa to the kappa power. ! A quadratic interpolation is done between values in a lookup table ! computed in gpkap. see documentation for fpkapx for details. ! Input values outside table range are reset to table extrema. ! The interpolation accuracy ranges from 12 decimal places ! at 100000 Pascals to 7 decimal places at 1000 Pascals. ! On the Cray, fpkap is over 4 times faster than exact calculation. ! This function should be expanded inline in the calling routine. ! ! Program History Log: ! 91-05-07 Iredell made into inlinable function ! 94-12-30 Iredell standardized kappa, ! increased range and accuracy ! 1999-03-01 Iredell f90 module ! 1999-03-24 Iredell table lookup ! ! Usage: pkap=fpkapq(p) ! ! Input argument list: ! p Real(krealfp) pressure in Pascals ! ! Output argument list: ! fpkapq Real(krealfp) p over 1e5 Pa to the kappa power ! ! Attributes: ! Language: Fortran 90. ! !$$$ implicit none real(krealfp) fpkapq real(krealfp),intent(in):: p integer jx real(krealfp) xj,dxj,fj1,fj2,fj3 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - xj=min(max(c1xpkap+c2xpkap*p,1._krealfp),real(nxpkap,krealfp)) jx=min(max(nint(xj),2),nxpkap-1) dxj=xj-jx fj1=tbpkap(jx-1) fj2=tbpkap(jx) fj3=tbpkap(jx+1) fpkapq=(((fj3+fj1)/2-fj2)*dxj+(fj3-fj1)/2)*dxj+fj2 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end function !------------------------------------------------------------------------------- !> This function raises surface pressure over 1e5 Pa to the kappa power !! using a rational weighted chebyshev approximation. !! The numerator is of order 2 and the denominator is of order 4. !! The pressure range is 40000-110000 Pa and kappa is defined in fpkapx(). !>\param[in] p real, surface pressure in Pascals p should be in the !! range 40000 to 110000 !\param[out] fpkapo real, p over 1e5 Pa to the kappa power function fpkapo(p) !$$$ Subprogram documentation block ! ! Subprogram: fpkapo raise surface pressure to the kappa power. ! Author: Phillips org: w/NMC2X2 Date: 29 dec 82 ! ! Abstract: Raise surface pressure over 1e5 Pa to the kappa power ! using a rational weighted chebyshev approximation. ! The numerator is of order 2 and the denominator is of order 4. ! The pressure range is 40000-110000 Pa and kappa is defined in fpkapx. ! The accuracy of this approximation is almost 8 decimal places. ! On the Cray, fpkap is over 10 times faster than exact calculation. ! This function should be expanded inline in the calling routine. ! ! Program History Log: ! 91-05-07 Iredell made into inlinable function ! 94-12-30 Iredell standardized kappa, ! increased range and accuracy ! 1999-03-01 Iredell f90 module ! ! Usage: pkap=fpkapo(p) ! ! Input argument list: ! p Real(krealfp) surface pressure in Pascals ! p should be in the range 40000 to 110000 ! ! Output argument list: ! fpkapo Real(krealfp) p over 1e5 Pa to the kappa power ! ! Attributes: ! Language: Fortran 90. ! !$$$ implicit none real(krealfp) fpkapo real(krealfp),intent(in):: p integer,parameter:: nnpk=2,ndpk=4 real(krealfp):: cnpk(0:nnpk)=(/3.13198449e-1,5.78544829e-2,& 8.35491871e-4/) real(krealfp):: cdpk(0:ndpk)=(/1.,8.15968401e-2,5.72839518e-4,& -4.86959812e-7,5.24459889e-10/) integer n real(krealfp) pkpa,fnpk,fdpk ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - pkpa=p*1.e-3_krealfp fnpk=cnpk(nnpk) do n=nnpk-1,0,-1 fnpk=pkpa*fnpk+cnpk(n) enddo fdpk=cdpk(ndpk) do n=ndpk-1,0,-1 fdpk=pkpa*fdpk+cdpk(n) enddo fpkapo=fnpk/fdpk ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end function !------------------------------------------------------------------------------- !> This function raises pressure over 1e5 Pa to the kappa power. !! Kappa is equal to rd/cp where rd and cp are physical constants. !>\param[in] p real, pressure in Pascals !\param[out] fpkapx real, p over 1e5 Pa to the kappa power elemental function fpkapx(p) !$$$ Subprogram documentation block ! ! Subprogram: fpkapx raise pressure to the kappa power. ! Author: Phillips org: w/NMC2X2 Date: 29 dec 82 ! ! Abstract: raise pressure over 1e5 Pa to the kappa power. ! Kappa is equal to rd/cp where rd and cp are physical constants. ! This function should be expanded inline in the calling routine. ! ! Program History Log: ! 94-12-30 Iredell made into inlinable function ! 1999-03-01 Iredell f90 module ! ! Usage: pkap=fpkapx(p) ! ! Input argument list: ! p Real(krealfp) pressure in Pascals ! ! Output argument list: ! fpkapx Real(krealfp) p over 1e5 Pa to the kappa power ! ! Attributes: ! Language: Fortran 90. ! !$$$ implicit none real(krealfp) fpkapx real(krealfp),intent(in):: p ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - fpkapx=(p/1.e5_krealfp)**con_rocp ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end function !------------------------------------------------------------------------------- !> This subroutine computes pressure to the 1/kappa table as a function of pressure !! for the table lookup function frkap(). !! Exact pressure to the 1/kappa values are calculated in subprogram frkapx(). !! The current implementation computes a table with a length !! of 5501 for pressures ranging up to 110000 Pascals. subroutine grkap !$$$ Subprogram documentation block ! ! Subprogram: grkap Compute coefficients for p**(1/kappa) ! Author: Phillips org: w/NMC2X2 Date: 29 dec 82 ! ! Abstract: Computes pressure to the 1/kappa table as a function of pressure ! for the table lookup function frkap. ! Exact pressure to the 1/kappa values are calculated in subprogram frkapx. ! The current implementation computes a table with a length ! of 5501 for pressures ranging up to 110000 Pascals. ! ! Program History Log: ! 94-12-30 Iredell ! 1999-03-01 Iredell f90 module ! 1999-03-24 Iredell table lookup ! ! Usage: call grkap ! ! Subprograms called: ! frkapx function to compute exact pressure to the 1/kappa ! ! Attributes: ! Language: Fortran 90. ! !$$$ implicit none integer jx real(krealfp) xmin,xmax,xinc,x,p ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - xmin=0._krealfp xmax=fpkapx(110000._krealfp) xinc=(xmax-xmin)/(nxrkap-1) c1xrkap=1.-xmin/xinc c2xrkap=1./xinc do jx=1,nxrkap x=xmin+(jx-1)*xinc p=x tbrkap(jx)=frkapx(p) enddo ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end subroutine !------------------------------------------------------------------------------- !> This subroutine raises pressure over 1e5 Pa to the 1/kappa power. !! A linear interpolation is done between values in a lookup table !! computed in grkap(). See documentation for frkapx() for details. !! Input values outside table range are reset to table extrema. !>\param[in] pkap real, p over 1e5 Pa to the kappa power !\param[out] frkap real, pressure in Pascals elemental function frkap(pkap) !$$$ Subprogram Documentation Block ! ! Subprogram: frkap raise pressure to the 1/kappa power. ! Author: N Phillips w/NMC2X2 Date: 30 dec 82 ! ! Abstract: Raise pressure over 1e5 Pa to the 1/kappa power. ! A linear interpolation is done between values in a lookup table ! computed in grkap. See documentation for frkapx for details. ! Input values outside table range are reset to table extrema. ! The interpolation accuracy is better than 7 decimal places. ! On the IBM, fpkap is about 4 times faster than exact calculation. ! This function should be expanded inline in the calling routine. ! ! Program History Log: ! 91-05-07 Iredell made into inlinable function ! 94-12-30 Iredell standardized kappa, ! increased range and accuracy ! 1999-03-01 Iredell f90 module ! 1999-03-24 Iredell table lookup ! ! Usage: p=frkap(pkap) ! ! Input argument list: ! pkap Real(krealfp) p over 1e5 Pa to the kappa power ! ! Output argument list: ! frkap Real(krealfp) pressure in Pascals ! ! Attributes: ! Language: Fortran 90. ! !$$$ implicit none real(krealfp) frkap real(krealfp),intent(in):: pkap integer jx real(krealfp) xj ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - xj=min(max(c1xrkap+c2xrkap*pkap,1._krealfp),real(nxrkap,krealfp)) jx=min(xj,nxrkap-1._krealfp) frkap=tbrkap(jx)+(xj-jx)*(tbrkap(jx+1)-tbrkap(jx)) ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end function !------------------------------------------------------------------------------- !> This function raises pressure over 1e5 Pa to the 1/kappa power. !! A quadratic interpolation is done between values in a lookup table !! computed in grkap(). see documentation for frkapx() for details. !! Input values outside table range are reset to table extrema. !>\param[in] pkap real, p over 1e5 Pa to the kappa power !\param[out] frkapq real, pressure in Pascals elemental function frkapq(pkap) !$$$ Subprogram Documentation Block ! ! Subprogram: frkapq raise pressure to the 1/kappa power. ! Author: N Phillips w/NMC2X2 Date: 30 dec 82 ! ! Abstract: Raise pressure over 1e5 Pa to the 1/kappa power. ! A quadratic interpolation is done between values in a lookup table ! computed in grkap. see documentation for frkapx for details. ! Input values outside table range are reset to table extrema. ! The interpolation accuracy is better than 11 decimal places. ! On the IBM, fpkap is almost 4 times faster than exact calculation. ! This function should be expanded inline in the calling routine. ! ! Program History Log: ! 91-05-07 Iredell made into inlinable function ! 94-12-30 Iredell standardized kappa, ! increased range and accuracy ! 1999-03-01 Iredell f90 module ! 1999-03-24 Iredell table lookup ! ! Usage: p=frkapq(pkap) ! ! Input argument list: ! pkap Real(krealfp) p over 1e5 Pa to the kappa power ! ! Output argument list: ! frkapq Real(krealfp) pressure in Pascals ! ! Attributes: ! Language: Fortran 90. ! !$$$ implicit none real(krealfp) frkapq real(krealfp),intent(in):: pkap integer jx real(krealfp) xj,dxj,fj1,fj2,fj3 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - xj=min(max(c1xrkap+c2xrkap*pkap,1._krealfp),real(nxrkap,krealfp)) jx=min(max(nint(xj),2),nxrkap-1) dxj=xj-jx fj1=tbrkap(jx-1) fj2=tbrkap(jx) fj3=tbrkap(jx+1) frkapq=(((fj3+fj1)/2-fj2)*dxj+(fj3-fj1)/2)*dxj+fj2 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end function !------------------------------------------------------------------------------- !> This function raise pressure over 1e5 Pa to the 1/kappa power. !! Kappa is equal to rd/cp where rd and cp are physical constants. !>\param[in] pkap real, p over 1e5 Pa to the kappa power !\param[out] frkapx real, pressure in Pascals elemental function frkapx(pkap) !$$$ Subprogram documentation block ! ! Subprogram: frkapx raise pressure to the 1/kappa power. ! Author: Phillips org: w/NMC2X2 Date: 29 dec 82 ! ! Abstract: raise pressure over 1e5 Pa to the 1/kappa power. ! Kappa is equal to rd/cp where rd and cp are physical constants. ! This function should be expanded inline in the calling routine. ! ! Program History Log: ! 94-12-30 Iredell made into inlinable function ! 1999-03-01 Iredell f90 module ! ! Usage: p=frkapx(pkap) ! ! Input argument list: ! pkap Real(krealfp) p over 1e5 Pa to the kappa power ! ! Output argument list: ! frkapx Real(krealfp) pressure in Pascals ! ! Attributes: ! Language: Fortran 90. ! !$$$ implicit none real(krealfp) frkapx real(krealfp),intent(in):: pkap ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - frkapx=pkap**(1/con_rocp)*1.e5_krealfp ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end function !------------------------------------------------------------------------------- !> This subroutine computes lifting condensation level temperature table !! as a function of temperature and dewpoint depression for function ftlcl(). !! Lifting condensation level temperature is calculated in subprogram ftlclx() !! The current implementation computes a table with a first dimension !! of 151 for temperatures ranging from 180.0 to 330.0 Kelvin !! and a second dimension of 61 for dewpoint depression ranging from !! 0 to 60 Kelvin. subroutine gtlcl !$$$ Subprogram Documentation Block ! ! Subprogram: gtlcl Compute equivalent potential temperature table ! Author: N Phillips w/NMC2X2 Date: 30 dec 82 ! ! Abstract: Compute lifting condensation level temperature table ! as a function of temperature and dewpoint depression for function ftlcl. ! Lifting condensation level temperature is calculated in subprogram ftlclx ! The current implementation computes a table with a first dimension ! of 151 for temperatures ranging from 180.0 to 330.0 Kelvin ! and a second dimension of 61 for dewpoint depression ranging from ! 0 to 60 Kelvin. ! ! Program History Log: ! 1999-03-01 Iredell f90 module ! ! Usage: call gtlcl ! ! Subprograms called: ! (ftlclx) inlinable function to compute LCL temperature ! ! Attributes: ! Language: Fortran 90. ! !$$$ implicit none integer jx,jy real(krealfp) xmin,xmax,ymin,ymax,xinc,yinc,x,y,tdpd,t ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - xmin=180._krealfp xmax=330._krealfp ymin=0._krealfp ymax=60._krealfp xinc=(xmax-xmin)/(nxtlcl-1) c1xtlcl=1.-xmin/xinc c2xtlcl=1./xinc yinc=(ymax-ymin)/(nytlcl-1) c1ytlcl=1.-ymin/yinc c2ytlcl=1./yinc do jy=1,nytlcl y=ymin+(jy-1)*yinc tdpd=y do jx=1,nxtlcl x=xmin+(jx-1)*xinc t=x tbtlcl(jx,jy)=ftlclx(t,tdpd) enddo enddo ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end subroutine !------------------------------------------------------------------------------- !> This function computes temperature at the lifting condensation level !! from temperature and dewpoint depression. !! A bilinear interpolation is done between values in a lookup table !! computed in gtlcl(). See documentation for ftlclx() for details. !! Input values outside table range are reset to table extrema. !>\param[in] t real, LCL temperature in Kelvin !>\param[in] tdpd real, dewpoint depression in Kelvin !\param[out] ftlcl real, temperature at the LCL in Kelvin elemental function ftlcl(t,tdpd) !$$$ Subprogram Documentation Block ! ! Subprogram: ftlcl Compute LCL temperature ! Author: N Phillips w/NMC2X2 Date: 30 dec 82 ! ! Abstract: Compute temperature at the lifting condensation level ! from temperature and dewpoint depression. ! A bilinear interpolation is done between values in a lookup table ! computed in gtlcl. See documentation for ftlclx for details. ! Input values outside table range are reset to table extrema. ! The interpolation accuracy is better than 0.0005 Kelvin. ! On the Cray, ftlcl is ? times faster than exact calculation. ! This function should be expanded inline in the calling routine. ! ! Program History Log: ! 1999-03-01 Iredell f90 module ! ! Usage: tlcl=ftlcl(t,tdpd) ! ! Input argument list: ! t Real(krealfp) LCL temperature in Kelvin ! tdpd Real(krealfp) dewpoint depression in Kelvin ! ! Output argument list: ! ftlcl Real(krealfp) temperature at the LCL in Kelvin ! ! Attributes: ! Language: Fortran 90. ! !$$$ implicit none real(krealfp) ftlcl real(krealfp),intent(in):: t,tdpd integer jx,jy real(krealfp) xj,yj,ftx1,ftx2 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - xj=min(max(c1xtlcl+c2xtlcl*t,1._krealfp),real(nxtlcl,krealfp)) yj=min(max(c1ytlcl+c2ytlcl*tdpd,1._krealfp),real(nytlcl,krealfp)) jx=min(xj,nxtlcl-1._krealfp) jy=min(yj,nytlcl-1._krealfp) ftx1=tbtlcl(jx,jy)+(xj-jx)*(tbtlcl(jx+1,jy)-tbtlcl(jx,jy)) ftx2=tbtlcl(jx,jy+1)+(xj-jx)*(tbtlcl(jx+1,jy+1)-tbtlcl(jx,jy+1)) ftlcl=ftx1+(yj-jy)*(ftx2-ftx1) ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end function !------------------------------------------------------------------------------- !> This function computes temperature at the lifting condensation level !! from temperature and dewpoint depression. !! A biquadratic interpolation is done between values in a lookup table !! computed in gtlcl(). see documentation for ftlclx() for details. !! Input values outside table range are reset to table extrema. !>\param[in] t real, LCL temperature in Kelvin !>\param[in] tdpd real, dowpoint depression in Kelvin !\param[out] ftlcl real, temperature at the LCL in Kelvin elemental function ftlclq(t,tdpd) !$$$ Subprogram Documentation Block ! ! Subprogram: ftlclq Compute LCL temperature ! Author: N Phillips w/NMC2X2 Date: 30 dec 82 ! ! Abstract: Compute temperature at the lifting condensation level ! from temperature and dewpoint depression. ! A biquadratic interpolation is done between values in a lookup table ! computed in gtlcl. see documentation for ftlclx for details. ! Input values outside table range are reset to table extrema. ! The interpolation accuracy is better than 0.000003 Kelvin. ! On the Cray, ftlclq is ? times faster than exact calculation. ! This function should be expanded inline in the calling routine. ! ! Program History Log: ! 1999-03-01 Iredell f90 module ! ! Usage: tlcl=ftlclq(t,tdpd) ! ! Input argument list: ! t Real(krealfp) LCL temperature in Kelvin ! tdpd Real(krealfp) dewpoint depression in Kelvin ! ! Output argument list: ! ftlcl Real(krealfp) temperature at the LCL in Kelvin ! ! Attributes: ! Language: Fortran 90. ! !$$$ implicit none real(krealfp) ftlclq real(krealfp),intent(in):: t,tdpd integer jx,jy real(krealfp) xj,yj,dxj,dyj real(krealfp) ft11,ft12,ft13,ft21,ft22,ft23,ft31,ft32,ft33 real(krealfp) ftx1,ftx2,ftx3 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - xj=min(max(c1xtlcl+c2xtlcl*t,1._krealfp),real(nxtlcl,krealfp)) yj=min(max(c1ytlcl+c2ytlcl*tdpd,1._krealfp),real(nytlcl,krealfp)) jx=min(max(nint(xj),2),nxtlcl-1) jy=min(max(nint(yj),2),nytlcl-1) dxj=xj-jx dyj=yj-jy ft11=tbtlcl(jx-1,jy-1) ft12=tbtlcl(jx-1,jy) ft13=tbtlcl(jx-1,jy+1) ft21=tbtlcl(jx,jy-1) ft22=tbtlcl(jx,jy) ft23=tbtlcl(jx,jy+1) ft31=tbtlcl(jx+1,jy-1) ft32=tbtlcl(jx+1,jy) ft33=tbtlcl(jx+1,jy+1) ftx1=(((ft31+ft11)/2-ft21)*dxj+(ft31-ft11)/2)*dxj+ft21 ftx2=(((ft32+ft12)/2-ft22)*dxj+(ft32-ft12)/2)*dxj+ft22 ftx3=(((ft33+ft13)/2-ft23)*dxj+(ft33-ft13)/2)*dxj+ft23 ftlclq=(((ftx3+ftx1)/2-ftx2)*dyj+(ftx3-ftx1)/2)*dyj+ftx2 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end function !------------------------------------------------------------------------------- !> This function computes temperature at the lifting condensation level !! from temperature and dewpoint depression. the formula used is !! a polynomial taken from Phillips mstadb routine which empirically !! approximates the original exact implicit relationship. !>\param[in] t real, temperature in Kelvin !>\param[in] tdpd real, dewpoint depression in Kelvin !\param[out] ftlclo real, temperature at the LCL in Kelvin function ftlclo(t,tdpd) !$$$ Subprogram documentation block ! ! Subprogram: ftlclo Compute LCL temperature. ! Author: Phillips org: w/NMC2X2 Date: 29 dec 82 ! ! Abstract: Compute temperature at the lifting condensation level ! from temperature and dewpoint depression. the formula used is ! a polynomial taken from Phillips mstadb routine which empirically ! approximates the original exact implicit relationship. ! (This kind of approximation is customary (inman, 1969), but ! the original source for this particular one is not yet known. -MI) ! Its accuracy is about 0.03 Kelvin for a dewpoint depression of 30. ! This function should be expanded inline in the calling routine. ! ! Program History Log: ! 91-05-07 Iredell made into inlinable function ! 1999-03-01 Iredell f90 module ! ! Usage: tlcl=ftlclo(t,tdpd) ! ! Input argument list: ! t Real(krealfp) temperature in Kelvin ! tdpd Real(krealfp) dewpoint depression in Kelvin ! ! Output argument list: ! ftlclo Real(krealfp) temperature at the LCL in Kelvin ! ! Attributes: ! Language: Fortran 90. ! !$$$ implicit none real(krealfp) ftlclo real(krealfp),intent(in):: t,tdpd real(krealfp),parameter:: clcl1= 0.954442e+0,clcl2= 0.967772e-3,& clcl3=-0.710321e-3,clcl4=-0.270742e-5 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ftlclo=t-tdpd*(clcl1+clcl2*t+tdpd*(clcl3+clcl4*t)) ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end function !------------------------------------------------------------------------------- !> This function computes temperature at the lifting condensation level !! from temperature and dewpoint depression. A parcel lifted !! adiabatically becomes saturated at the lifting condensation level. !! The water model assumes a perfect gas, constant specific heats !! for gas and liquid, and neglects the volume of the liquid. !! The model does account for the variation of the latent heat !! of condensation with temperature. The ice option is not included. !! The Clausius-Clapeyron equation is integrated from the triple point !! to get the formulas: !!\n pvlcl=con_psat*(trlcl**xa)*exp(xb*(1.-trlcl)) !!\n pvdew=con_psat*(trdew**xa)*exp(xb*(1.-trdew)) !!\n where pvlcl is the saturated parcel vapor pressure at the LCL, !! pvdew is the unsaturated parcel vapor pressure initially, !! trlcl is ttp/tlcl and trdew is ttp/tdew. The adiabatic lifting !! of the parcel is represented by the following formula !!\n pvdew=pvlcl*(t/tlcl)**(1/kappa) !!\n This formula is inverted by iterating Newtonian approximations !! until tlcl is found to within 1.e-6 Kelvin. Note that the minimum !! returned temperature is 180 Kelvin. !>\param[in] t real, temperature in Kelvin !>\param[in] tdpd real, dewpoint depression in Kelvin !\param[out] ftlclx real, temperature at the LCL in Kelvin elemental function ftlclx(t,tdpd) !$$$ Subprogram documentation block ! ! Subprogram: ftlclx Compute LCL temperature. ! Author: Iredell org: w/NMC2X2 Date: 25 March 1999 ! ! Abstract: Compute temperature at the lifting condensation level ! from temperature and dewpoint depression. A parcel lifted ! adiabatically becomes saturated at the lifting condensation level. ! The water model assumes a perfect gas, constant specific heats ! for gas and liquid, and neglects the volume of the liquid. ! The model does account for the variation of the latent heat ! of condensation with temperature. The ice option is not included. ! The Clausius-Clapeyron equation is integrated from the triple point ! to get the formulas ! pvlcl=con_psat*(trlcl**xa)*exp(xb*(1.-trlcl)) ! pvdew=con_psat*(trdew**xa)*exp(xb*(1.-trdew)) ! where pvlcl is the saturated parcel vapor pressure at the LCL, ! pvdew is the unsaturated parcel vapor pressure initially, ! trlcl is ttp/tlcl and trdew is ttp/tdew. The adiabatic lifting ! of the parcel is represented by the following formula ! pvdew=pvlcl*(t/tlcl)**(1/kappa) ! This formula is inverted by iterating Newtonian approximations ! until tlcl is found to within 1.e-6 Kelvin. Note that the minimum ! returned temperature is 180 Kelvin. ! ! Program History Log: ! 1999-03-25 Iredell ! ! Usage: tlcl=ftlclx(t,tdpd) ! ! Input argument list: ! t Real(krealfp) temperature in Kelvin ! tdpd Real(krealfp) dewpoint depression in Kelvin ! ! Output argument list: ! ftlclx Real(krealfp) temperature at the LCL in Kelvin ! ! Attributes: ! Language: Fortran 90. ! !$$$ implicit none real(krealfp) ftlclx real(krealfp),intent(in):: t,tdpd real(krealfp),parameter:: terrm=1.e-4,tlmin=180.,tlminx=tlmin-5. real(krealfp) tr,pvdew,tlcl,ta,pvlcl,el,dpvlcl,terr,terrp integer i ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - tr=con_ttp/(t-tdpd) pvdew=con_psat*(tr**con_xpona)*exp(con_xponb*(1.-tr)) tlcl=t-tdpd do i=1,100 tr=con_ttp/tlcl ta=t/tlcl pvlcl=con_psat*(tr**con_xpona)*exp(con_xponb*(1.-tr))*ta**(1/con_rocp) el=con_hvap+con_dldt*(tlcl-con_ttp) dpvlcl=(el/(con_rv*t**2)+1/(con_rocp*tlcl))*pvlcl terr=(pvlcl-pvdew)/dpvlcl tlcl=tlcl-terr if(abs(terr).le.terrm.or.tlcl.lt.tlminx) exit enddo ftlclx=max(tlcl,tlmin) ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end function !------------------------------------------------------------------------------- !> This subroutine computes all physics function tables. Lookup tables are !! set up for computing saturation vapor pressure, dewpoint temperature, !! equivalent potential temperature, moist adiabatic temperature and humidity, !! pressure to the kappa, and lifting condensation level temperature. subroutine gfuncphys !$$$ Subprogram Documentation Block ! ! Subprogram: gfuncphys Compute all physics function tables ! Author: N Phillips w/NMC2X2 Date: 30 dec 82 ! ! Abstract: Compute all physics function tables. Lookup tables are ! set up for computing saturation vapor pressure, dewpoint temperature, ! equivalent potential temperature, moist adiabatic temperature and humidity, ! pressure to the kappa, and lifting condensation level temperature. ! ! Program History Log: ! 1999-03-01 Iredell f90 module ! ! Usage: call gfuncphys ! ! Subprograms called: ! gpvsl compute saturation vapor pressure over liquid table ! gpvsi compute saturation vapor pressure over ice table ! gpvs compute saturation vapor pressure table ! gtdpl compute dewpoint temperature over liquid table ! gtdpi compute dewpoint temperature over ice table ! gtdp compute dewpoint temperature table ! gthe compute equivalent potential temperature table ! gtma compute moist adiabat tables ! gpkap compute pressure to the kappa table ! grkap compute pressure to the 1/kappa table ! gtlcl compute LCL temperature table ! ! Attributes: ! Language: Fortran 90. ! !$$$ implicit none ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - call gpvsl call gpvsi call gpvs call gtdpl call gtdpi call gtdp call gthe call gtma call gpkap call grkap call gtlcl ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end subroutine !------------------------------------------------------------------------------- end module