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 ! 13-08-13 S. Moorthi - Threading ! ! 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 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 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,jsta:jend),intent(in) :: P1,T1 REAL,dimension(IM,jsta:jend),intent(inout):: Q1,RH REAL ES,QC integer :: I,J !*************************************************************** ! ! START CALRH. ! !$omp parallel do private(i,j,es,qc) DO J=JSTA,JEND DO I=1,IM IF (T1(I,J) < SPVAL .AND. P1(I,J) < SPVAL.AND.Q1(I,J)/=SPVAL) THEN ! IF (ABS(P1(I,J)) > 1.0) THEN IF (P1(I,J) > 1.0) THEN ES = MIN(FPVSNEW(T1(I,J)),P1(I,J)) QC = CON_EPS*ES/(P1(I,J)+CON_EPSM1*ES) ! QC=PQ0/P1(I,J)*EXP(A2*(T1(I,J)-A3)/(T1(I,J)-A4)) RH(I,J) = min(1.0,max(Q1(I,J)/QC,rhmin)) q1(i,j) = rh(i,j)*qc ! BOUNDS CHECK ! ! IF (RH(I,J) > 1.0) THEN ! RH(I,J) = 1.0 ! Q1(I,J) = RH(I,J)*QC ! ELSEIF (RH(I,J) < 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