SUBROUTINE CALRH_GFS(P1,T1,Q1,RH)
!$$$ SUBPROGRAM DOCUMENTATION BLOCK
! . . .
! SUBPROGRAM: CALRH COMPUTES RELATIVE HUMIDITY
! PRGRMMR: TREADON ORG: W/NP2 DATE: 92-12-22
!
! ABSTRACT:
! THIS ROUTINE COMPUTES RELATIVE HUMIDITY GIVEN PRESSURE,
! TEMPERATURE, SPECIFIC HUMIDITY. AN UPPER AND LOWER BOUND
! OF 100 AND 1 PERCENT RELATIVE HUMIDITY IS ENFORCED. WHEN
! THESE BOUNDS ARE APPLIED THE PASSED SPECIFIC HUMIDITY
! ARRAY IS ADJUSTED AS NECESSARY TO PRODUCE THE SET RELATIVE
! HUMIDITY.
! .
!
! PROGRAM HISTORY LOG:
! ??-??-?? DENNIS DEAVEN
! 92-12-22 RUSS TREADON - MODIFIED AS DESCRIBED ABOVE.
! 98-06-08 T BLACK - CONVERSION FROM 1-D TO 2-D
! 98-08-18 MIKE BALDWIN - MODIFY TO COMPUTE RH OVER ICE AS IN MODEL
! 98-12-16 GEOFF MANIKIN - UNDO RH COMPUTATION OVER ICE
! 00-01-04 JIM TUCCILLO - MPI VERSION
! 02-06-11 MIKE BALDWIN - WRF VERSION
!
! USAGE: CALL CALRH(P1,T1,Q1,RH)
! INPUT ARGUMENT LIST:
! P1 - PRESSURE (PA)
! T1 - TEMPERATURE (K)
! Q1 - SPECIFIC HUMIDITY (KG/KG)
!
! OUTPUT ARGUMENT LIST:
! RH - RELATIVE HUMIDITY (DECIMAL FORM)
! Q1 - ADJUSTED SPECIFIC HUMIDITY (KG/KG)
!
! OUTPUT FILES:
! NONE
!
! SUBPROGRAMS CALLED:
! UTILITIES:
! LIBRARY:
! NONE
!
! ATTRIBUTES:
! LANGUAGE: FORTRAN
! MACHINE : CRAY C-90
!$$$
!
use params_mod, only: rhmin
use ctlblk_mod, only: jsta, jend, spval, im, jm
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
implicit none
!
real,parameter:: con_rd =2.8705e+2 ! gas constant air (J/kg/K)
real,parameter:: con_rv =4.6150e+2 ! gas constant H2O
real,parameter:: con_eps =con_rd/con_rv
real,parameter:: con_epsm1 =con_rd/con_rv-1
! real,external::FPVSNEW
INTERFACE
ELEMENTAL FUNCTION FPVSNEW (t)
REAL FPVSNEW
REAL, INTENT(IN) :: t
END FUNCTION FPVSNEW
END INTERFACE
!
REAL,dimension(IM,JM),intent(in):: P1,T1
REAL,dimension(IM,JM),intent(inout):: Q1,RH
REAL ES,QC
integer :: I,J
!***************************************************************
!
! START CALRH.
!
DO J=JSTA,JEND
DO I=1,IM
IF (T1(I,J).LT.SPVAL .AND. P1(I,J).LT.SPVAL.AND.Q1(I,J)/=SPVAL) THEN
IF (ABS(P1(I,J)).GT.1) THEN
ES=FPVSNEW(T1(I,J))
ES=MIN(ES,P1(I,J))
QC=CON_EPS*ES/(P1(I,J)+CON_EPSM1*ES)
! QC=PQ0/P1(I,J)
! 1 *EXP(A2*(T1(I,J)-A3)/(T1(I,J)-A4))
RH(I,J)=Q1(I,J)/QC
! BOUNDS CHECK
!
IF (RH(I,J).GT.1.0) THEN
RH(I,J)=1.0
Q1(I,J)=RH(I,J)*QC
ENDIF
! IF (RH(I,J).LT.0.01) THEN
IF (RH(I,J).LT.RHmin) THEN !use smaller RH limit for stratosphere
RH(I,J)=RHmin
Q1(I,J)=RH(I,J)*QC
ENDIF
ENDIF
ELSE
RH(I,J)=SPVAL
ENDIF
ENDDO
ENDDO
RETURN
END
elemental function fpvsnew(t)
!$$$ Subprogram Documentation Block
!
! Subprogram: fpvsnew 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=fpvsnew(t)
!
! Input argument list:
! t Real(krealfp) temperature in Kelvin
!
! Output argument list:
! fpvsnew Real(krealfp) saturation vapor pressure in Pascals
!
! Attributes:
! Language: Fortran 90.
!
!$$$
implicit none
integer,parameter:: nxpvs=7501
real,parameter:: con_ttp =2.7316e+2 ! temp at H2O 3pt
real,parameter:: con_psat =6.1078e+2 ! pres at H2O 3pt
real,parameter:: con_cvap =1.8460e+3 ! spec heat H2O gas (J/kg/K)
real,parameter:: con_cliq =4.1855e+3 ! spec heat H2O liq
real,parameter:: con_hvap =2.5000e+6 ! lat heat H2O cond
real,parameter:: con_rv =4.6150e+2 ! gas constant H2O
real,parameter:: con_csol =2.1060e+3 ! spec heat H2O ice
real,parameter:: con_hfus =3.3358e+5 ! lat heat H2O fusion
real,parameter:: tliq=con_ttp
real,parameter:: tice=con_ttp-20.0
real,parameter:: dldtl=con_cvap-con_cliq
real,parameter:: heatl=con_hvap
real,parameter:: xponal=-dldtl/con_rv
real,parameter:: xponbl=-dldtl/con_rv+heatl/(con_rv*con_ttp)
real,parameter:: dldti=con_cvap-con_csol
real,parameter:: heati=con_hvap+con_hfus
real,parameter:: xponai=-dldti/con_rv
real,parameter:: xponbi=-dldti/con_rv+heati/(con_rv*con_ttp)
real tr,w,pvl,pvi
real fpvsnew
real,intent(in):: t
integer jx
real xj,x,tbpvs(nxpvs),xp1
real xmin,xmax,xinc,c2xpvs,c1xpvs
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
xmin=180.0
xmax=330.0
xinc=(xmax-xmin)/(nxpvs-1)
! c1xpvs=1.-xmin/xinc
c2xpvs=1./xinc
c1xpvs=1.-xmin*c2xpvs
! xj=min(max(c1xpvs+c2xpvs*t,1.0),real(nxpvs,krealfp))
xj=min(max(c1xpvs+c2xpvs*t,1.0),float(nxpvs))
jx=min(xj,float(nxpvs)-1.0)
x=xmin+(jx-1)*xinc
tr=con_ttp/x
if(x.ge.tliq) then
tbpvs(jx)=con_psat*(tr**xponal)*exp(xponbl*(1.-tr))
elseif(x.lt.tice) then
tbpvs(jx)=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))
tbpvs(jx)=w*pvl+(1.-w)*pvi
endif
xp1=xmin+(jx-1+1)*xinc
tr=con_ttp/xp1
if(xp1.ge.tliq) then
tbpvs(jx+1)=con_psat*(tr**xponal)*exp(xponbl*(1.-tr))
elseif(xp1.lt.tice) then
tbpvs(jx+1)=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))
tbpvs(jx+1)=w*pvl+(1.-w)*pvi
endif
fpvsnew=tbpvs(jx)+(xj-jx)*(tbpvs(jx+1)-tbpvs(jx))
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
end function fpvsnew