C$$$ SUBPROGRAM DOCUMENTATION BLOCK C . . . . C SUBPROGRAM: HB2P INTERPOLATE FROM HYB-B TO ISOBARIC GRIDS C C PRGMMR: BENJAMIN, STAN ORG: ERL/PROFS DATE: 93-01-12 C C ABSTRACT: INTERPOLATE GRIDS FROM HYBRID-B LEVELS (NO INTERSECTING C SURFACES) TO ISOBARIC LEVELS C C PROGRAM HISTORY LOG: C 93-01-12 S. BENJAMIN ORIGINAL VERSION C C USAGE: CALL HB2P (g3,XT, C 1 H,T,RHP,UP,VP,XP,NX3D,PLEV,NZP) C C INPUT ARGUMENT LIST: C g3 - real C contains the following C PT - REAL HYBRID-B PRESSURE (Pa) C TH - REAL HYBRID-B VIRTUAL POTENTIAL TEMP (K) C M - REAL HYBRID-B height (m) C UT - REAL HYBRID-B U WIND COMPONENT (M/S) C VT - REAL HYBRID-B V WIND COMPONENT (M/S) C vertical velocity, cloud variables, TKE C XT - REAL other HYBRID-B variables C RHB - REAL HYBRID-B Relative humidity (0-1) C NX3D - INTEGER NUMBER OF other hybrid-b variables C PLEV - REAL ISOBARIC LEVELS IN MB C NZP - INTEGER NUMBER OF LEVELS TO INTERPOLATE TO C OUTPUT ARGUMENT LIST: C H - REAL ISOBARIC HEIGHT (M) C T - REAL ISOBARIC TEMPERATURE (K) (NON-VIRTUAL) C RHP - REAL ISOBARIC RELATIVE HUMIDITY (0-1) C UP - REAL ISOBARIC U WIND COMPONENT (M/S) C VP - REAL ISOBARIC V WIND COMPONENT (M/S) C XP - REAL other ISOBARIC variables C C REMARKS: SEE BENJAMIN ET AL., 1993, AMS CONFERENCE ON AVIATION C WEATHER FORECASTING, FOR DISCUSSION OF HYBRID-B COORDINATES. C THE MAXIMUM NUMBER OF ISOBARIC LEVELS IS SET AT 37. THIS CAN C BE INCREASED, IF NECESSARY. C C ATTRIBUTES: C LANGUAGE: FORTRAN-77 + EXTENSIONS C MACHINE: DEC - VAX, VMS C C$$$ SUBROUTINE HB2P (g3,XT,RHB, 1 H,T,RHP,UP,VP,XP,NX3D,PLEV,NZP) C *** This program converts a hybrid analysis/forecast in hybrid-b C coordinates (non-intersecting levels) into C typical isobaric variables on isobaric surfaces. C *** NZP is the number of pressure levels to which to interpolate. C The maximum value C for NZP is 39, which is the number of levels C spaced 25 mb apart between 1050 and 100 mb. If you need C to interpolate to more than 39 isobaric levels, change C NZP25_P in MAPSCON and recompile. C NZTN_P is the number of hybrid levels. An extra level C is needed for a higher isentropic level (500 K) with which C to extrapolate to higher isobaric levels. C NX3D is the number of extra 3d variables to be interpolated C to isobaric surfaces. Currently, these can be used for cloud C and TKE variables passed in thru g3. IMPLICIT NONE c INCLUDE 'MAPSCON' INTEGER 1 I,J,K,NZP,N,NX3D,IVAR integer nx_p,ny_p,nztn_p,nvar3_p parameter (nx_p=151,ny_p=113,nztn_p=40,nvar3_p=14) integer p_p,theta_p,u_p,v_p,h_p parameter (p_p=1,theta_p=3,u_p=5,v_p=6,h_p=2) integer nzp25_p,nzp_p parameter (nzp25_p=37,nzp_p=8) C *** HYBRID VARIABLES REAL G3 (NX_P,NY_P,NZTN_P,NVAR3_P) C pressure - Pa C virtual pot temp - K C height - m C u component of wind - m/s C v component of wind - m/s C other 3-d variables (e.g., pressure vertical velocity) 1 ,RHB(NX_P,NY_P,NZTN_P) C relative humidity (0-1) on hybrid levels 1 ,xt(NX_P,NY_P,NZTN_P) C extra variable - potential vorticity REAL 1 DUM (NX_P,NY_P,NZTN_P) 1 ,DUMMY (NX_P,NY_P) 1 ,TEMP,Q,PM1 1 ,EX1 C R*gamma/Cp 1 ,GAMMA 1 ,GAM,GAMD 1 ,PSFC,TSFC,ZSFC 1 ,T1,T5,Z1,Z5 C sfc values 1 ,C (4,NZTN_P+1) C work array for VERINT 1 ,XI(NZTN_P+1) C work array for VERINT C *** THETA VARIABLES REAL 1 EXNER (NX_P,NY_P,NZTN_P) C *** PRESSURE VARIABLES REAL H (NX_P,NY_P,NZP25_P) C height - m 1 ,T (NX_P,NY_P,NZP25_P) C temperature (non-virtual) C - K 1 ,RHP (NX_P,NY_P,NZP25_P) C relative humidity C - 0-1 scale 1 ,UP (NX_P,NY_P,NZP25_P) C u component of wind - m/s 1 ,VP (NX_P,NY_P,NZP25_P) C v component of wind - m/s 1 ,XP (NX_P,NY_P,NZP25_P,NX3D) C other 3-d variables, e.g., vertical velocity C *** PLEV are the pressure levels in mb REAL PLEV (NZP) 1 ,EXLEV (NZP25_P) C *** SMOOTHING VARIABLES REAL PMLEV (NZP_P) DATA PMLEV / 85000., 70000., 50000., 40000., 1 30000., 25000., 20000., 15000. / INTEGER KSM (NZP25_P) C - closest mandatory level INTEGER NSM (5, NZP_P) C - look-up table of recommended number of smoothing C passes, function of level and variable C H T RH U V DATA NSM / 3, 2, 1, 3, 3, 1 2, 1, 1, 3, 3, 1 2, 1, 1, 2, 2, 1 2, 1, 1, 1, 1, 1 2, 1, 1, 1, 1, 1 2, 1, 1, 1, 1, 1 2, 1, 1, 1, 1, 1 3, 1, 1, 2, 2 / C DATA EX1 /0.1903643/ C R*gamma/g DATA GAMMA /0.0065/ C 6.5 K/km DATA GAMD /0.0100/ C 10 K/km REAL CPD_P, ROVCP_P, G0_P, RD_P PARAMETER ( + CPD_P = 1004.686, C J/(kg*K) Specific heat at constant pressure C = 3.5* + ROVCP_P = .285714, C nd / + G0_P = 9.80665 C m/(s*s) Sea-level acceleration of gravity + ,RD_P = 287.04 1 ) INTEGER IG,K5 C IG - LEVEL USED FOR EXTRAPOLATING BELOW SURFACE C TYPE *,' ENTER LEVEL FOR EXTRAPOLATING BELOW SURFACE' C ACCEPT *,IG IG = 4 C TYPE *,' ENTER TOP LEVEL FOR COMPUTING LAPSE RATE' C ACCEPT *,K5 K5 = 5 DO 25 K = 1,NZTN_P DO 25 J = 1,NY_P DO 25 I = 1,NX_P EXNER(I,J,K) = CPD_P*(g3(I,J,K,p_p)/100000.)**ROVCP_P 25 CONTINUE do 27 k=1,nztn_p print *,' exner,pt,th,rhb',exner(15,15,k),g3(15,15,k,p_p), 1 g3(15,15,k,theta_p),rhb(15,15,k) 27 continue DO 10 K = 1,NZP EXLEV(K) = CPD_P*(PLEV(K)/100000.)**ROVCP_P PRINT *,' K,P,EX',K,PLEV(K),EXLEV(K) 10 CONTINUE CALL VERINT (0,1,EXNER,g3(1,1,1,u_p),DUM,NX_P,NY_P,NZTN_P,C,XI, 1 UP,NZP,EXLEV) CALL VERINT (0,1,EXNER,g3(1,1,1,v_p),DUM,NX_P,NY_P,NZTN_P,C,XI, 1 VP,NZP,EXLEV) DO 11 IVAR = 1,NX3D-1 CALL VERINT (0,1,EXNER,g3(1,1,1,6+IVAR),DUM,NX_P,NY_P, 1 NZTN_P,C,XI,XP(1,1,1,IVAR),NZP,EXLEV) 11 CONTINUE CALL VERINT (0,1,EXNER,RHB,DUM,NX_P,NY_P,NZTN_P,C,XI, 1 RHP,NZP,EXLEV) CALL VERINT (0,1,EXNER,xt,DUM,NX_P,NY_P,NZTN_P,C,XI, 1 xp(1,1,1,nx3d),NZP,EXLEV) CALL VERINT (0,1,EXNER,g3(1,1,1,theta_p),DUM,NX_P,NY_P,NZTN_P, 1 C,XI,T,NZP,EXLEV) CALL VERINT (0,1,EXNER,g3(1,1,1,h_p),DUM,NX_P,NY_P,NZTN_P, 1 C,XI,H,NZP,EXLEV) PRINT *,' H(15,15,5)',H(15,15,5) PRINT *,' T(15,15,3)',T(15,15,3) PRINT *,' RHP(15,15,5)',RHP(15,15,5) PRINT *,' UP(15,15,5)',UP(15,15,5) PRINT *,' VP(15,15,5)',VP(15,15,5) DO 12 IVAR = 1,NX3D PRINT *,' XP(15,15,5),IVAR',XP(15,15,5,IVAR) 12 CONTINUE C________________________________________________________ C________________________________________________________ C Extrapolate and interpolate near ground C________________________________________________________ C________________________________________________________ DO 50 J = 1,NY_P DO 50 I = 1,NX_P PSFC = 100000.*(EXNER(I,J,IG)/CPD_P)**3.5 TSFC = EXNER(I,J,IG)*g3(I,J,IG,theta_p)/CPD_P ZSFC = g3(I,J,IG,h_p) Cxxxxxx ERROR, should be theta at level 1 T1 = EXNER(I,J,1)*g3(I,J,IG,theta_p)/CPD_P Z1 = g3(I,J,1,h_p) T5 = EXNER(I,J,K5)*g3(I,J,K5,theta_p)/CPD_P Z5 = g3(I,J,K5,h_p) GAM = (T1-T5)/(Z5-Z1) GAM = MIN(GAMD,MAX(GAM,GAMMA)) EX1 = RD_P*GAM/G0_P DO 50 K = 1,NZP C--------------------------------- C The interpolation is considered in 3 layers: A,B,C C A = above IGth hybrid-b level (about 15-25 hPa above ground) C B = between IGth level and model surface C C = below model surface (beneath ground) C C ht temp RH u v x C ------ ------ ----- ----- ---- ----- C A interp interp inter inter inter inter C - - - - - - - - - - - - - - - - - - - - - - - - C B reduce reduce inter inter inter inter C _______________________________________________ C / / / / / / / / / surface / / / / / / / / / / / C C reduce reduce =k=1 =k=1 =k=1 =k=1 C C - layer A IF (EXLEV(K).LT.EXNER(I,J,IG)) THEN C --- get T from theta T(I,J,K) = T(I,J,K)*EXLEV(K)/CPD_P TEMP = T(I,J,K) C - layer B ELSE IF (EXLEV(K).GT.EXNER(I,J,IG) .AND. 1 EXLEV(K).LT.EXNER(I,J,1)) THEN C --- we need to do some extrapolating downward at this grid point. C Set RH, u, v = sfc values for isobaric levels beneath ground. C in pascals T(I,J,K) = TSFC*(PLEV(K)/PSFC)**EX1 H(I,J,K) = ZSFC - (T(I,J,K)-TSFC)/GAM ELSE C - layer C -- IF (EXLEV(K).GT.EXNER(I,J,1)) THEN T(I,J,K) = TSFC*(PLEV(K)/PSFC)**EX1 H(I,J,K) = ZSFC - (T(I,J,K)-TSFC)/GAM UP(I,J,K) = g3(I,J,1,u_p) VP(I,J,K) = g3(I,J,1,v_p) RHP(I,J,K)=RHB(I,J,1) DO 49 IVAR=1,NX3D-1 XP(I,J,K,IVAR) = g3(I,J,1,6+IVAR) 49 CONTINUE xp(i,j,k,nx3d) = xt(i,j,1) END IF c if (k.eq.1.and.(t(i,j,k).lt.200..or.t(i,j,k).gt.320.)) c 1 print *,' bad reduction',i,j,k,t(i,j,k),h(i,j,k), c 1 tsfc,zsfc,psfc,gam,gamd,gamma 50 CONTINUE C *** De-virtualize virtual temperatures DO 55 K = 1,NZP DO 55 J = 1,NY_P DO 55 I = 1,NX_P CALL TV2TQ (T(I,J,K),RHP(I,J,K)/100.,PLEV(K)/100., 1 TEMP,Q) c if (k.eq.1.and.(temp.lt.200..or.temp.gt.320.)) then c print *,' bad tv2tq',i,j,k,temp,t(i,j,k),rhp(i,j,k) c do k1=1,25 c print *,t(i,j,k1),rhp(i,j,k1),pt(i,j,k1),th(i,j,k1) c end do c end if T(I,J,K) = TEMP 55 CONTINUE C --- Find closest mandatory levels to each isobaric level DO 58 K = 1,NZP CALL FNDLVL (PLEV(K),PMLEV,NZP_P,KSM(K)) 58 CONTINUE C________________________________________________________ C________________________________________________________ C Extrapolate to upper isobaric levels, if necessary C Check from 175 to 100 mb. C________________________________________________________ C________________________________________________________ do 52 k=nzp-3,nzp do 52 j=1,ny_p do 52 i=1,nx_p if (t(i,j,k).lt.0.) then t(i,j,k) = g3(i,j,nztn_p,theta_p)*exner(i,j,nztn_p)/cpd_p h(i,j,k) = g3(i,j,nztn_p,h_p) + rd_p*t(i,j,k)/g0_p + *alog(g3(i,j,nztn_p,p_p)/plev(k)) rhp(i,j,k) = rhb(i,j,nztn_p) up (i,j,k) = g3 (i,j,nztn_p,u_p) vp (i,j,k) = g3 (i,j,nztn_p,v_p) do 515 ivar = 1,nx3d-1 xp(i,j,k,ivar) = g3(i,j,nztn_p,6+ivar) 515 continue xp(i,j,k,nx3d) = xt(i,j,nztn_p) end if 52 continue C --- Smooth by number of passes for recommended level C at nearest mandatory level (NSM) DO 60 K = 1,NZP DO 591 N=1,NSM(1,KSM(K)) CALL SMOOTH(H(1,1,K),DUMMY,NX_P,NY_P,0.5) 591 CONTINUE DO 592 N=1,NSM(2,KSM(K)) CALL SMOOTH(T(1,1,K),DUMMY,NX_P,NY_P,0.5) 592 CONTINUE DO 593 N=1,NSM(3,KSM(K)) CALL SMOOTH(RHP(1,1,K),DUMMY,NX_P,NY_P,0.5) DO 5925 IVAR=1,NX3D CALL SMOOTH(XP(1,1,K,IVAR),DUMMY,NX_P,NY_P,0.5) 5925 CONTINUE 593 CONTINUE DO 594 N=1,NSM(4,KSM(K)) CALL SMOOTH(UP(1,1,K),DUMMY,NX_P,NY_P,0.5) CALL SMOOTH(VP(1,1,K),DUMMY,NX_P,NY_P,0.5) 594 CONTINUE 60 CONTINUE DO 70 K = 1,NZP DO 70 J = 1,NY_P DO 70 I = 1,NX_P RHP(I,J,K) = MAX(0.,MIN(100.,RHP(I,J,K))) 70 CONTINUE RETURN END C$$$ SUBPROGRAM DOCUMENTATION BLOCK C . . . . C SUBPROGRAM: VERINT PERFORM VERTICAL INTERPOLATION C PRGMMR: STAN BENJAMIN ORG: FSL/PROFS DATE: 90-09-19 C C ABSTRACT: COMPUTE VERTICAL HERMITE OR LINEAR INTERPOLATION. C C PROGRAM HISTORY LOG: C C C USAGE: CALL VERINT (ITYPE, CALCDF, SUPORT, F, DF, C IDIN, JDIN, KDIN1, C, XI, C FNEW, KDIN2, NEWLEV) C INPUT ARGUMENT LIST: C--MODEV C ITYPE - INTEGER TYPE OF INTERPOLATION SCHEME. C 0 = LINEAR, 1 = HERMITE C CALCDF - INTEGER DECIDE WHETHER USE THE DERIVATIVE VALUE C FROM OUTSIDE OF ROUTINE OR NOT. C 0 = USE, 1 = CALCULATE IN THE ROUTINE. C--QI C SUPORT - REAL ARRAY SUPORT(IDIN,JDIN,KDIN1) C SUPPORT POINTS IN NEW VERTICAL C COORDINATES ON OLD LEVEL. E.G. THETA C--FI C F - REAL ARRAY F(IDIN,JDIN,KDIN1) C VALUES AT SUPPORT POINT TO INTERPOLATE. C--DFI C DF - REAL ARRAY DF(IDIN,JDIN,KDIN1) C DERIVATIVE OF F FOR HERMITE INTERPOLATION. C IDIN - INTEGER NUMBER OF GRID POINTS IN X-DIRECTION. C JDIN - INTEGER NUMBER OF GRID POINTS IN Y-DIRECTION. C KDIN1 - INTEGER NUMBER OF GRID POINTS IN VERTICAL LEVEL. C C - REAL ARRAY C(4,KDIN1) C COEFFICIENTS OF CUBIC POLYNOMIAL BETWEEN 2 POINTS. C XI - REAL ARRAY XI(KDIN1) C TEMPORARY STORAGE VALUE TO BE INTERPOLATED. C--KDOU C KDIN2 - INTEGER NUMBER OF GRID POINTS IN VERTICAL LEVEL C ON NEW COORDINATES(NEWLEV). C--QKO C NEWLEV - REAL ARRAY NEWLEV(KDIN2) C NEW COORDINATES LEVELS TO EVALUATE C INTERPOLATING POLYNOMIALS AT. C C OUTPUT ARGUMENT LIST: C--FO C FNEW - REAL ARRAY FNEW(IDIN,JDIN,KDIN2) C INTERPOLATED VALUES OF F ON NEW COORDINATES LEVEL. C C REMARKS: THIS PROGRAM WORKS FOR ONLY MONOTONICALLY DECREASING C SUPPORT POINT. ALSO, WHEN THE DERIVATIVE VALUE OF F (=DF) ARE C USED FROM PASSING ARGUEMENT ( WHEN CALCDF = 0), THE VALUES OF C DF CORRESPOND TO THE MONOTONICALLY DECREASING SUPPORT POINTS. C C ATTRIBUTES: C LANGUAGE: FORTRAN-77 C MACHINE: VAX-8800 C$$$ C SUBROUTINE VERINT (ITYPE, CALCDF, SUPORT, F, DF, - IDIN, JDIN, KDIN1, C, XI, - FNEW, KDIN2, NEWLEV) C IMPLICIT NONE INTEGER ITYPE, CALCDF, IDIN, JDIN, KDIN1, KDIN2 REAL SUPORT, F, FNEW, NEWLEV REAL DF, XI, C DIMENSION SUPORT(IDIN,JDIN,KDIN1), F(IDIN,JDIN,KDIN1), - DF(IDIN,JDIN,KDIN1) DIMENSION XI(KDIN1), C(4,KDIN1) DIMENSION FNEW(IDIN,JDIN,KDIN2), - NEWLEV(KDIN2) INTEGER LIN, HERM PARAMETER ( LIN = 0, HERM = 1 ) C C IF (ITYPE .EQ. LIN) THEN CALL LININT ( SUPORT, CALCDF, F, DF, IDIN, JDIN, KDIN1, - FNEW, KDIN2, NEWLEV, XI) ELSEIF (ITYPE .EQ. HERM) THEN CALL HERINT ( SUPORT, CALCDF, F, DF, IDIN, JDIN, KDIN1, - FNEW, KDIN2, NEWLEV, C, XI) ELSE PRINT*, ' INVALID INTERPOLATION TYPE. ITYPE =', ITYPE ENDIF RETURN END C$$$ SUBPROGRAM DOCUMENTATION BLOCK C . . . . C SUBPROGRAM: LININT PERFORM LINEAR INTERPOLATION. C PRGMMR: BENJAMIN, STAN ORG: ERL/PROFS DATE: 93-01-18 C C ABSTRACT: PERFORM LINEAR INTERPOLATION. C C PROGRAM HISTORY LOG: C 88/05/31 S. BENJAMIN ORIGINAL VERSION C C USAGE: CALL LININT ( SUPORT, CALCDF, F, DF, IDIN, JDIN, KDIN1, C FNEW, KDIN2, NEWLEV, XI) C C INPUT ARGUMENT LIST: C C OUTPUT ARGUMENT LIST: C C REMARKS: NONE C C ATTRIBUTES: C LANGUAGE: FORTRAN-77 + EXTENSIONS C MACHINE: DEC - VAX, VMS C C$$$ C------------------------------------------------------------------- C PERFORM LINEAR INTERPOLATION. C------------------------------------------------------------------- SUBROUTINE LININT ( SUPORT, CALCDF, F, DF, IDIN, JDIN, KDIN1, - FNEW, KDIN2, NEWLEV, XI) C IMPLICIT NONE INTEGER CALCDF, IDIN, JDIN, KDIN1, KDIN2, I, J, K, KK, PFIND REAL XBAR REAL SUPORT(IDIN,JDIN,KDIN1), F(IDIN,JDIN,KDIN1) REAL DF(IDIN,JDIN,KDIN1), - XI(KDIN1) C REAL FNEW(IDIN,JDIN,KDIN2), - NEWLEV(KDIN2) C IF (CALCDF.EQ.1) .THEN C C CALCULATE DERIVATIVE OF F. C DO 100 J=1, JDIN DO 100 I=1, IDIN DO 100 K=1, KDIN1-1 DF(I,J,K) = (F(I,J,K+1) - F(I,J,K))/ - (SUPORT(I,J,K+1) - SUPORT(I,J,K)) 100 CONTINUE C DO 110 J=1, JDIN DO 110 I=1, IDIN DF(I,J,KDIN1) = DF(I,J,KDIN1-1) 110 CONTINUE ENDIF C C GATHER SUPPORT POINT C DO 160 J=1, JDIN DO 160 I=1, IDIN DO 140 K=1, KDIN1 XI(K)= SUPORT(I,J,K) 140 CONTINUE C C EVALUATE INTERPOLATION. C DO 150 K=1, KDIN2 XBAR = NEWLEV(K) KK = PFIND(XBAR, XI, KDIN1-1) FNEW(I,J,K) = F(I,J,KK) + DF(I,J,KK)*(XBAR - XI(KK)) 150 CONTINUE 160 CONTINUE C RETURN END C$$$ SUBPROGRAM DOCUMENTATION BLOCK C . . . . C SUBPROGRAM: PFIND RETURN LINEAR INTERP INDEX C PRGMMR: BENJAMIN, STAN ORG: ERL/PROFS DATE: 93-01-18 C C ABSTRACT: RETURNS THE INDEX VALUE AT XBAR OF THE LINEAR INTERPOLATION C ON N INTERVALS WITH BREAKPOINT SEQUENCE XI. C C PROGRAM HISTORY LOG: C 88/05/31 S. BENJAMIN ORIGINAL VERSION C C USAGE: CALL PFIND(XBAR, XI, N) C C INPUT ARGUMENT LIST: C C OUTPUT ARGUMENT LIST: C C REMARKS: NONE C C ATTRIBUTES: C LANGUAGE: FORTRAN-77 + EXTENSIONS C MACHINE: DEC - VAX, VMS C C$$$ C------------------------------------------------------------------- C RETURNS THE INDEX VALUE AT XBAR OF THE LINEAR INTERPOLATION C ON N INTERVALS WITH BREAKPOINT SEQUENCE XI. C------------------------------------------------------------------- INTEGER FUNCTION PFIND (XBAR, XI, N) IMPLICIT NONE INTEGER N, I, J REAL XBAR, XI(N+1) DATA I /1/ IF (XBAR .LE. XI(I)) THEN DO 555 J=I, N IF (XBAR .GT. XI(J+1)) GO TO 777 555 CONTINUE J = N ELSE DO 666 J=I-1, 1, -1 IF (XBAR .GE. XI(J)) GO TO 777 666 CONTINUE J = 1 ENDIF 777 PFIND = J RETURN END C$$$ SUBPROGRAM DOCUMENTATION BLOCK C . . . . C SUBPROGRAM: HERINT PERFORM HERIMITE INTERPOLATION C PRGMMR: BENJAMIN, STAN ORG: ERL/PROFS DATE: 93-01-18 C C ABSTRACT: PERFORM HERMITE INTERPOLATION C C PROGRAM HISTORY LOG: C 88/05/31 S. BENJAMIN ORIGINAL VERSION C C USAGE: CALL HERINT ( SUPORT, CALCDF, F, DF, C 1 IDIN, JDIN, KDIN1, C 2 FNEW, KDIN2, NEWLEV, C, XI) C C INPUT ARGUMENT LIST: C C OUTPUT ARGUMENT LIST: C C REMARKS: NONE C C ATTRIBUTES: C LANGUAGE: FORTRAN-77 + EXTENSIONS C MACHINE: DEC - VAX, VMS C C$$$ C-------------------------------------------------------------------- C PERFORM HERMITE INTERPOLATION. C-------------------------------------------------------------------- SUBROUTINE HERINT ( SUPORT, CALCDF, F, DF, - IDIN, JDIN, KDIN1, - FNEW, KDIN2, NEWLEV, C, XI) C IMPLICIT NONE REAL XBAR, PCUBIC INTEGER CALCDF, IDIN, JDIN, KDIN1, KDIN2, I, J, K REAL SUPORT(IDIN,JDIN,KDIN1), F(IDIN,JDIN,KDIN1) REAL DF(IDIN,JDIN,KDIN1), - C(4,KDIN1), XI(KDIN1) C REAL FNEW(IDIN,JDIN,KDIN2), - NEWLEV(KDIN2) C IF (CALCDF.EQ.1) THEN C C CALCULATE THE DERIVATIVE OF F. C DO 200 K=2, KDIN1 - 1 DO 200 J=1, JDIN DO 200 I=1, IDIN DF(I,J,K) = (F(I,J,K+1) - F(I,J,K-1))/ - (SUPORT(I,J,K+1) - SUPORT(I,J,K-1)) 200 CONTINUE DO 210 J=1, JDIN DO 210 I=1, IDIN DF(I,J,1) = (F(I,J,2) - F(I,J,1))/ - (SUPORT(I,J,2) - SUPORT(I,J,1)) DF(I,J,KDIN1) = (F(I,J,KDIN1) - F(I,J,KDIN1-1))/ - (SUPORT(I,J,KDIN1) - SUPORT(I,J,KDIN1-1)) 210 CONTINUE ENDIF C C GATHER SUPPORT POINT AND CALCULATE THE CUBIC COEFFICIENT C OF EACH INTERVAL OF F. C DO 250 J=1, JDIN DO 250 I=1, IDIN DO 235 K=1,KDIN1 XI(K) = SUPORT(I,J,K) C(1,K) = F(I,J,K) C(2,K) = DF(I,J,K) 235 CONTINUE CALL CALCF( XI,C,KDIN1-1 ) DO 240 K=1,KDIN2 XBAR = NEWLEV(K) FNEW(I,J,K) = PCUBIC(XBAR,XI,C,KDIN1-1) 240 CONTINUE 250 CONTINUE RETURN END C$$$ SUBPROGRAM DOCUMENTATION BLOCK C . . . . C SUBPROGRAM: CALCF ??? C PRGMMR: BENJAMIN, STAN ORG: ERL/PROFS DATE: 93-01-18 C C ABSTRACT: ??? C C PROGRAM HISTORY LOG: C 88/05/31 S. BENJAMIN ORIGINAL VERSION C C USAGE: CALL CALCF (XI, C, N) C C INPUT ARGUMENT LIST: C C OUTPUT ARGUMENT LIST: C C REMARKS: C INPUT; C XI(1), ..., X(I+1) STRICTLY INCREASING SEQUENCE OF BREAKPOINTS. C C(1,I), C(2,I), VALUES AND FIRST DERIVATIVE AT XI(I), I=1,...,N+1 C C OUTPUT; C POLYNOMIAL COEFFICIENTS FOR C(1,I), C(2,I), C(3,I), C(4,I) OF C INTERVAL (XI(I), XI(I+1)). I= 1,..,N C------------------------------------------------------------------- C C ATTRIBUTES: C LANGUAGE: FORTRAN-77 + EXTENSIONS C MACHINE: DEC - VAX, VMS C C$$$ C------------------------------------------------------------------- C INPUT; C XI(1), ..., X(I+1) STRICTLY INCREASING SEQUENCE OF BREAKPOINTS. C C(1,I), C(2,I), VALUES AND FIRST DERIVATIVE AT XI(I), I=1,...,N+1 C OUTPUT; C POLYNOMIAL COEFFICIENTS FOR C(1,I), C(2,I), C(3,I), C(4,I) OF C INTERVAL (XI(I), XI(I+1)). I= 1,..,N C------------------------------------------------------------------- SUBROUTINE CALCF (XI, C, N) INTEGER N, I REAL C(4,N+1), XI(N+1), DIVDF1, DIVDF3, DX DO 12 I=1, N DX = XI(I+1) - XI(I) DIVDF1 = (C(1,I+1) - C(1,I))/DX DIVDF3 = C(2,I) + C(2,I+1) - 2.0*DIVDF1 C(3,I) = (DIVDF1 - C(2,I) - DIVDF3)/DX C(4,I) = DIVDF3/(DX*DX) 12 CONTINUE RETURN END C$$$ SUBPROGRAM DOCUMENTATION BLOCK C . . . . C SUBPROGRAM: PCUBIC PIECEWISE CUBIC FUNCTIONS C PRGMMR: BENJAMIN, STAN ORG: ERL/PROFS DATE: 93-01-18 C C ABSTRACT: PIECEWISE CUBIC FUNCTION C C PROGRAM HISTORY LOG: C 88/05/31 S. BENJAMIN ORIGINAL VERSION C C USAGE: PCUBIC (XBAR, XI, C, N) C C INPUT ARGUMENT LIST: C C OUTPUT ARGUMENT LIST: C C REMARKS: NONE C C ATTRIBUTES: C LANGUAGE: FORTRAN-77 + EXTENSIONS C MACHINE: DEC - VAX, VMS C C$$$ C------------------------------------------------------------------- C RETURNS THE VALUE AT XBAR OF THE PIECEWISE CUBIC FUNCTION C ON N INTERVALS WITH BREAKPOINT SEQUENCE XI AND COEFFICIENTS C. C------------------------------------------------------------------- REAL FUNCTION PCUBIC (XBAR, XI, C, N) IMPLICIT NONE INTEGER N, I, J REAL C(4, N), XBAR, XI(N+1), DX DATA I /1/ IF (XBAR .LE. XI(I)) THEN DO 222 J=I, N IF (XBAR .GT. XI(J+1)) GO TO 444 222 CONTINUE J = N ELSE DO 333 J=I-1, 1, -1 IF (XBAR .LE. XI(J)) GO TO 444 333 CONTINUE J = 1 ENDIF 444 I = J DX = XBAR - XI(I) PCUBIC = C(1,I) + DX*(C(2,I) + DX*(C(3,I) + DX*C(4,I))) RETURN END C$$$ SUBPROGRAM DOCUMENTATION BLOCK C . . . . C SUBPROGRAM: TV2TQ DEVIRTUALIZE TEMPERATURE C C PRGMMR: BENJAMIN, STAN ORG: ERL/PROFS DATE: 93-01-12 C C ABSTRACT: CALCULATES TEMPERATURE AND SPECIFIC HUMIDITY FROM C VIRTUAL TEMPERATURE, PRESSURE, AND RELATIVE HUMIDITY C C PROGRAM HISTORY LOG: C 92 S. BENJAMIN ORIGINAL VERSION C 94-05-16 S. BENJAMIN IMPROVE ACCURACY BY USING SPECIFIC C HUMIDITY INSTEAD OF MIXING RATIO. C ALSO HAS FASTER CONVERGENCE. C ADD ORDER OF ACCURACY TO 0.6078 C CONSTANT. C 94-05-23 S. BENJAMIN LIMIT ALGORITHM TO 3 ITERATIONS C FOR MORE SPEED. C 95-01-02 S. BENJAMIN NEW VERSION C 95-07-28 S. Benjamin Limit number of iterations C C USAGE: CALL TV2TQ(TV, RH, P, T, Q) C C INPUT ARGUMENT LIST: C TV - REAL VIRTUAL TEMPERATURE IN KELVIN C RH - REAL RELATIVE HUMIDITY C P - REAL PRESSURE IN MILIBARS C C OUTPUT ARGUMENT LIST: C T - REAL TEMPERATURE IN KELVIN C Q - REAL SPECIFIC HUMIDITY (G/G) C C REMARKS: NONE C C ATTRIBUTES: C LANGUAGE: FORTRAN-77 + EXTENSIONS C MACHINE: DEC - VAX, VMS C C$$$ SUBROUTINE TV2TQ(TV,RH,P,T,Q) C C THIS ROUTINE RETURNS TEMPERATURE (KELVIN) AND SPECIFIC HUMIDITY C (G/G) GIVEN VIRTUAL TEMPERATURE (KELVIN), PRESSURE (MB), AND C RELATIVE HUMIDITY (RANGE 0.0-1.0). C IT USES AN ITERATIVE NEWTON-RAPHSON TECHNIQUE. FOUR ITERATIONS ARE C GENERALLY ADEQUATE TO PROVIDE CONVERGENCE TO 5 DECIMAL PLACES. C THE WOBUS FUNCTION FOR SATURATION VAPOR PRESSURE OVER LIQUID WATER C IS USED C S.G. BENJAMIN 13-DEC-1990 Original version C S.G. BENJAMIN 02-JAN-1995 New version IMPLICIT NONE REAL TV,RH,P,T,Q REAL TG,QG,FAC,FAC1,TVG,DTVDT,TGNEW REAL ESW,E INTEGER IT IT = 0 C -- TG - guess at non-virtual temperature E = RH*ESW(TV) C 0.37803 = 0.6078 * 0.62197 TG = TV/(1.+0.37803*E/(P-0.37803*E)) 10 CONTINUE IT = IT + 1 C -- QG - guess at specific humidity E = RH*ESW(TG) QG = 0.62197*E/(P-0.37803*E) FAC = (1.+0.6078*QG) C -- TVG - guess at virtual temperature TVG = TG*FAC E = RH*ESW(TG+0.1) FAC1= (1.+0.37803*E/(P-0.37803*E)) C -- DTVDT - d(Tv)/dT estimated over 0.1 K interval C from 1st term Taylor series expansion DTVDT = FAC + TG*(FAC1-FAC)/0.1 TGNEW = TG + (TV-TVG)/DTVDT C PRINT *,',TV,TVG,TGNEW,QG,FAC,FAC1,DTVDT',TV,TVG,TGNEW, C 1 QG,FAC,FAC1,DTVDT IF (ABS(TV-TVG).LT.0.01) GO TO 20 IF (IT.GT.6) GO TO 20 TG = TGNEW GO TO 10 20 CONTINUE T = TGNEW E = RH*ESW(T) Q = 0.62197*E/(P-0.37803*E) RETURN END C$$$ SUBPROGRAM DOCUMENTATION BLOCK C . . . . C SUBPROGRAM: FNDLVL FIND THE CLOSEST ANALYSIS LEVEL TO OBS C PRGMMR: TRACY SMITH ORG: FSL/PROFS DATE: 90-07-19 C C ABSTRACT: FIND WHAT ANALYSIS LEVEL IS CLOSEST TO THE OBSERVATION C FOR SINGLE LEVEL DATA. C C PROGRAM HISTORY LOG: C C C USAGE: CALL FNDLVL(OB,VCOORD,NZ,LVLC) C INPUT ARGUMENT LIST: C OB - REAL OBSERVATION RELATING TO VERTICAL COORDINATES C VCOORD - REAL ARRAY VERTICAL COORDINATES (LEVELS) C NZ - INTEGER NUMBER OF LEVELS C C OUTPUT ARGUMENT LIST: C LVLC - INTEGER LEVEL THE OBSERVATION IS AT C C REMARKS: C C ATTRIBUTES: C LANGUAGE: FORTRAN-77 C MACHINE: NAS-9000 C$$$ SUBROUTINE FNDLVL(OB,VERT_COORD,NZ,LVLC) IMPLICIT NONE c INCLUDE 'MAPSCON' REAL VERT_COORD(*) C REAL OB,DIFF1,DIFF2 INTEGER LVLC,I,NZ C C*** Decide if the parameter increases or decreases IF(VERT_COORD(1).LT.VERT_COORD(2)) THEN C*** If the parameter decreases DO 70 I=1, NZ-1 IF(OB.LT.VERT_COORD(I)) GO TO 71 70 CONTINUE 71 IF(I.NE.1) THEN DIFF1=VERT_COORD(I)-OB DIFF2=OB-VERT_COORD(I-1) IF(DIFF1.GT.DIFF2) THEN LVLC=I-1 ELSE LVLC=I END IF ELSE LVLC=1 END IF C WRITE(6,915) OB, LVLC, VERT_COORD(LVLC) C 915 FORMAT(' CLOSEST LEVEL TO OBS=',F8.0,' IS ',I4, C , ' WHERE OBS=', F8.0) C*** If the variable increases ELSE DO 80 I=1, NZ-1 IF(OB.GT.VERT_COORD(I)) GO TO 81 80 CONTINUE 81 IF(I.NE.1) THEN DIFF1=OB-VERT_COORD(I) DIFF2=VERT_COORD(I-1)-OB IF(DIFF1.GT.DIFF2) THEN LVLC=I-1 ELSE LVLC=I END IF ELSE LVLC=1 END IF C WRITE(6,916) OB, LVLC, VERT_COORD(LVLC) C916 FORMAT(' CLOSEST LEVEL TO OBS=',F8.0,' IS ',I4, C , ' WHERE OBS=', F8.0) ENDIF C ! the variable increases or decreases RETURN END C$$$ SUBPROGRAM DOCUMENTATION BLOCK C . . . . C SUBPROGRAM: SMOOTH SMOOTH A METEOROLOGICAL FIELD C PRGMMR: STAN BENJAMIN ORG: FSL/PROFS DATE: 90-06-15 C C ABSTRACT: SHAPIRO SMOOTHER. C C PROGRAM HISTORY LOG: C 85-12-09 S. BENJAMIN ORIGINAL VERSION C C USAGE: CALL SMOOTH (FIELD,HOLD,IX,IY,SMTH) C INPUT ARGUMENT LIST: C FIELD - REAL ARRAY FIELD(IX,IY) C METEOROLOGICAL FIELD C HOLD - REAL ARRAY HOLD(IX,2) C HOLDING THE VALUE FOR FIELD C IX - INTEGER X COORDINATES OF FIELD C IY - INTEGER Y COORDINATES OF FIELD C SMTH - REAL C C OUTPUT ARGUMENT LIST: C FIELD - REAL ARRAY FIELD(IX,IY) C SMOOTHED METEOROLOGICAL FIELD C C REMARKS: REFERENCE: SHAPIRO, 1970: "SMOOTHING, FILTERING, AND C BOUNDARY EFFECTS", REV. GEOPHYS. SP. PHYS., 359-387. C THIS FILTER IS OF THE TYPE C Z(I) = (1-S)Z(I) + S(Z(I+1)+Z(I-1))/2 C FOR A FILTER WHICH IS SUPPOSED TO DAMP 2DX WAVES COMPLETELY C BUT LEAVE 4DX AND LONGER WITH LITTLE DAMPING, C IT SHOULD BE RUN WITH 2 PASSES USING SMTH (OR S) OF 0.5 C AND -0.5. C C ATTRIBUTES: C LANGUAGE: FORTRAN-77 + EXTENSIONS C MACHINE: NAS-9000, VAX, UNIX C$$$ C********************************************************************** C********************************************************************** SUBROUTINE SMOOTH (FIELD,HOLD,IX,IY,SMTH) C********************************************************************** C********************************************************************** REAL FIELD(IX,IY), 1 HOLD (IX,2) SMTH1 = 0.25 * SMTH * SMTH SMTH2 = 0.5 * SMTH * (1.-SMTH) SMTH3 = (1.-SMTH) * (1.-SMTH) SMTH4 = (1.-SMTH) SMTH5 = 0.5 * SMTH I1 = 2 I2 = 1 DO 30 J=2,IY-1 IT = I1 I1 = I2 I2 = IT DO 10 I = 2,IX-1 SUM1 = FIELD (I-1,J+1) + FIELD (I-1,J-1) 1 + FIELD (I+1,J+1) + FIELD (I+1,J-1) SUM2 = FIELD (I ,J+1) + FIELD (I+1,J ) 1 + FIELD (I ,J-1) + FIELD (I-1,J ) HOLD(I,I1) = SMTH1*SUM1 + SMTH2*SUM2 + SMTH3*FIELD(I,J) 10 CONTINUE IF (J.EQ.2) GO TO 200 DO 20 I=2,IX-1 FIELD(I,J-1) = HOLD(I,I2) 20 CONTINUE 200 CONTINUE 30 CONTINUE DO 40 I = 2,IX-1 FIELD (I,IY-1) = HOLD(I,I1) 40 CONTINUE DO 50 I = 2,IX-1 FIELD(I,1) = SMTH4* FIELD(I,1) 1 + SMTH5 * (FIELD(I-1,1) + FIELD(I+1,1)) FIELD(I,IY) = SMTH4* FIELD(I,IY) 1 + SMTH5 * (FIELD(I-1,IY) + FIELD(I+1,IY)) 50 CONTINUE DO 60 J = 2,IY-1 FIELD(1,J) = SMTH4* FIELD(1,J) 1 + SMTH5 * (FIELD(1,J-1) + FIELD(1,J+1)) FIELD(IX,J) = SMTH4* FIELD(IX,J) 1 + SMTH5 * (FIELD(IX,J-1) + FIELD(IX,J+1)) 60 CONTINUE RETURN END C$$$ SUBPROGRAM DOCUMENTATION BLOCK C . . . . C SUBPROGRAM: ESW COMPUTE THE SATURATION VAPOR PRESSURE C PRGMMR: STAN BENJAMIN ORG: FSL/PROFS DATE: 90-07-19 C C ABSTRACT: FUNCTION TO RETURNS THE SATURATION VAPOR PRESSURE C ESW (MILLIBARS) OVER LIQUID WATER GIVEN THE TEMPERATURE C T (CELSIUS OR KELVIN). C C PROGRAM HISTORY LOG: C 82-05-17 TOM SCHLATTER ORIGINAL VERSION C C USAGE: ESW(T) C INPUT ARGUMENT LIST: C T - REAL TEMPERATURE (CELSIUS OR KELVIN) C C OUTPUT ARGUMENT LIST: C C REMARKS: THE POLYNOMIAL APPROXIMATION BELOW IS DUE TO HERMAN WOBUS, C A MATHEMATICIAN WHO WORKED AT THE NAVY WEATHER RESEARCH FACILITY, C NORFOLK, VIRGINIA. THE COEFFICIENTS OF C THE POLYNOMIAL WERE CHOSEN TO FIT THE VALUES IN TABLE 94 ON C PP. 351-353 OF THE SMITHSONIAN METEOROLOGICAL TABLES BY ROLAND C LIST (6TH EDITION). THE APPROXIMATION IS VALID FOR C -50 < T < 100C. C C ATTRIBUTES: C LANGUAGE: FORTRAN-77 C MACHINE: NAS-9000 C$$$ FUNCTION ESW(T1) C C THIS FUNCTION RETURNS THE SATURATION VAPOR PRESSURE ESW (MILLIBARS) C OVER LIQUID WATER GIVEN THE TEMPERATURE T (CELSIUS). C C BAKER,SCHLATTER 17-MAY-1982 Original version C C THE POLYNOMIAL APPROXIMATION BELOW IS DUE TO HERMAN WOBUS, A C MATHEMATICIAN WHO WORKED AT THE NAVY WEATHER RESEARCH FACILITY, C NORFOLK, VIRGINIA. THE COEFFICIENTS OF THE C POLYNOMIAL WERE CHOSEN TO FIT THE VALUES IN TABLE 94 ON PP. 351-353 C C OF THE SMITHSONIAN METEOROLOGICAL TABLES BY ROLAND LIST (6TH C EDITION). THE APPROXIMATION IS VALID FOR -50 < T < 100C. C C ES0 = SATURATION VAPOR RESSURE OVER LIQUID WATER AT 0C DATA ES0/6.1078/ T = T1 IF(T.GT.150.) T=T-273.15 POL = 0.99999683 + T*(-0.90826951E-02 + 1 T*(0.78736169E-04 + T*(-0.61117958E-06 + 2 T*(0.43884187E-08 + T*(-0.29883885E-10 + 3 T*(0.21874425E-12 + T*(-0.17892321E-14 + 4 T*(0.11112018E-16 + T*(-0.30994571E-19))))))))) ESW = ES0/POL**8 RETURN END