C----------------------------------------------------------------------- C----------------------------------------------------------------------- PROGRAM RRVENTS C$$$ MAIN PROGRAM DOCUMENTATION BLOCK C . . . . C MAIN PROGRAM: PREPFITS QM MARKS TO PREPBUFR FILE C PRGMMR: JANJIC ORG: NP22 DATE: 99-01-20 C C ABSTRACT: PROGRAM ADDS QM MARKS TO PREPBUFR FILE (LIKE GLOBAL) C C PROGRAM HISTORY LOG: C 12-12-07 J. WOOLLEN added tolerences C 13-01 W. EBISUZAKI conversion to linux C 17-02 W. EBISUZAKI program failed because the tolerences C for the location test were too tight. (New obs_proc) C text file(unit=20) lat/lon to xxx.xx degrees C old prepbufr: lat/lon to xxx.xx degrees C test: abs(lat1-lat2) .le. 0.001 C new prepbufr: lat/lon to xxx.xxxxxx C test can fail C so new test abs(lat1-lat2) .le. 0.006 same for lon C C changed tol(x1-x2) -> tolx(x1,x2) C changed tol(y1-y2) -> tolx(y2,y2) C changed PARAMETER(MPP=50000,MXX=200000) to PARAMETER(MPP=90000,MXX=300000) C MPP/MXX were not a problem but increase number to be safe C C USAGE: MAIN PROGRAM C C INPUT FILES: C UNIT20 TEXT FILE WITH QM MARKS C UNIT21 INPUT BUFR file C OUTPUT FILES: C UNIT51 Output BUFR file with obs & forecast profiles C C SUBPROGRAMS CALLED: C UFBCPY C UFBINT C CLOSBF C OPENMB C IREADMG C IREADSB C WRITSB C C EXIT STATES: C COND = 0 - NORMAL EXIT C COND = 1 - ERROR C C ATTRIBUTES: C LANGUAGE: FORTRAN 90 C MACHINE : IBM SP C C$$$ PARAMETER(MPP=90000,MXX=300000) CHARACTER*80 HDSTR,OBSTR,PVSTR,QVSTR,TVSTR,UVSTR,VVSTR CHARACTER*8 SUBSET,VAR,SID CHARACTER*8 PID(MPP),QID(MXX),TID(MXX),UID(MXX),VID(MXX) DIMENSION EVN(11) DIMENSION PCD(MPP),POO(MPP),PFC(MPP),PAN(MPP),PXX(MPP),PYY(MPP) DIMENSION PEL(MPP),PDH(MPP),PPP(MPP),PQM(MPP),PER(MPP),PRT(MPP) DIMENSION QCD(MXX),QOO(MXX),QFC(MXX),QAN(MXX),QXX(MXX),QYY(MXX) DIMENSION QEL(MXX),QDH(MXX),QPR(MXX),QQM(MXX),QER(MXX),QRT(MXX) DIMENSION TCD(MXX),TOO(MXX),TFC(MXX),TAN(MXX),TXX(MXX),TYY(MXX) DIMENSION TEL(MXX),TDH(MXX),TPR(MXX),TQM(MXX),TER(MXX),TRT(MXX) DIMENSION UCD(MXX),UOO(MXX),UFC(MXX),UAN(MXX),UXX(MXX),UYY(MXX) DIMENSION UEL(MXX),UDH(MXX),UPR(MXX),UQM(MXX),UER(MXX),URT(MXX) DIMENSION VCD(MXX),VOO(MXX),VFC(MXX),VAN(MXX),VXX(MXX),VYY(MXX) DIMENSION VEL(MXX),VDH(MXX),VPR(MXX),VQM(MXX),VER(MXX),VRT(MXX) EQUIVALENCE (HDR(1),SID) logical tol, tolx REAL*8 HDR(10) REAL*8 OBS(10,255),PEV(10,255),QEV(10,255),QMS(10,255) REAL*8 TEV(10,255),UEV(10,255),VEV(10,255) CHARACTER*50 string DATA HDSTR /'SID XOB YOB DHR ELV TYP '/ DATA OBSTR /'POB QOB TOB UOB VOB CAT ZOB'/ DATA PVSTR /'POB PQM PRC PPC PFC PAN PER'/ DATA QVSTR /'QOB QQM QRC QPC QFC QAN QER'/ DATA TVSTR /'TOB TQM TRC TPC TFC TAN TER'/ DATA UVSTR /'UOB WQM WRC WPC UFC UAN WER'/ DATA VVSTR /'VOB NUL NUL NUL VFC VAN NUL'/ DATA ANLCODE /10/ DATA LUNEV /20/ DATA LUNIN /21/ DATA LUOUT /51/ DATA BMISS /10E10/ C----------------------------------------------------------------------- tol(a,b) = abs(a-b).le..001 tolx(a,b) = abs(a-b).le..006 C----------------------------------------------------------------------- C CLEAR THE SUMMARY COUNTERS NPP = 0 ; NPOB = 0 ; NP3DX = 0 ; NPR2X = 0 NQQ = 0 ; NQOB = 0 ; NQ3DX = 0 ; NQR2X = 0 NTT = 0 ; NTOB = 0 ; NT3DX = 0 ; NTR2X = 0 NUU = 0 ; NUOB = 0 ; NU3DX = 0 ; NUR2X = 0 NVV = 0 ; NVOB = 0 ; NV3DX = 0 ; NVR2X = 0 C READ AND STORE THE ASCII TEXT EVENTS C ------------------------------------ write(*,*) '>> read unit=',lunev 10 READ(LUNEV,11,END=20,err=101) VAR,ICD,EVN,SID 11 FORMAT(A3,I3.3,11E13.5,2X,A8) C WNE 11/2012 III = LJUST(SID) write(*,*) 'prior: SID=',SID CALL JSTCHR(SID,III) write(*,*) 'post: SID=',SID write(*,*) 'var=',var IF(VAR.EQ.'PP.') THEN IF(NPP+1.GT.MPP) CALL BORT('TOO MANY P EVENTS') NPP = NPP+1 PCD(NPP) = ICD POO(NPP) = EVN( 1) PFC(NPP) = EVN( 2) PAN(NPP) = EVN( 3) PXX(NPP) = EVN( 4) PYY(NPP) = EVN( 5) PEL(NPP) = EVN( 6) PDH(NPP) = EVN( 7) PPP(NPP) = EVN( 8) PQM(NPP) = EVN( 9) PER(NPP) = EVN(10) PRT(NPP) = EVN(11) PID(NPP) = SID ELSEIF(VAR.EQ.'TT.') THEN IF(NTT+1.GT.MXX) CALL BORT('TOO MANY T EVENTS') NTT = NTT+1 TCD(NTT) = ICD TOO(NTT) = EVN( 1)-273.15 TFC(NTT) = EVN( 2)-273.15 TAN(NTT) = EVN( 3)-273.15 TXX(NTT) = EVN( 4) TYY(NTT) = EVN( 5) TPR(NTT) = EVN( 6) TEL(NTT) = EVN( 7) TDH(NTT) = EVN( 8) TQM(NTT) = EVN( 9) TER(NTT) = EVN(10) TRT(NTT) = EVN(11) TID(NTT) = SID ELSEIF(VAR.EQ.'QQ.') THEN IF(NQQ+1.GT.MXX) CALL BORT('TOO MANY Q EVENTS') NQQ = NQQ+1 QCD(NQQ) = ICD QOO(NQQ) = EVN( 1)*1.E6; qoo(nqq) = nint(qoo(nqq)) QFC(NQQ) = EVN( 2)*1.E6; qfc(nqq) = nint(qfc(nqq)) QAN(NQQ) = EVN( 3)*1.E6; qan(nqq) = nint(qan(nqq)) QXX(NQQ) = EVN( 4) QYY(NQQ) = EVN( 5) QPR(NQQ) = EVN( 6) QEL(NQQ) = EVN( 7) QDH(NQQ) = EVN( 8) QQM(NQQ) = EVN( 9) QER(NQQ) = EVN(10) QRT(NQQ) = EVN(11) QID(NQQ) = SID if(qrt(nqq).eq.119.) nqq = nqq-1 ELSEIF(VAR.EQ.'UU.') THEN IF(NUU+1.GT.MXX) CALL BORT('TOO MANY U EVENTS') NUU = NUU+1 UCD(NUU) = ICD UOO(NUU) = EVN( 1) UFC(NUU) = EVN( 2) UAN(NUU) = EVN( 3) UXX(NUU) = EVN( 4) UYY(NUU) = EVN( 5) UPR(NUU) = EVN( 6) UEL(NUU) = EVN( 7) UDH(NUU) = EVN( 8) UQM(NUU) = EVN( 9) UER(NUU) = EVN(10) URT(NUU) = EVN(11) UID(NUU) = SID ELSEIF(VAR.EQ.'VV.') THEN IF(NVV+1.GT.MXX) CALL BORT('TOO MANY V EVENTS') NVV = NVV+1 VCD(NVV) = ICD VOO(NVV) = EVN( 1) VFC(NVV) = EVN( 2) VAN(NVV) = EVN( 3) VXX(NVV) = EVN( 4) VYY(NVV) = EVN( 5) VPR(NVV) = EVN( 6) VEL(NVV) = EVN( 7) VDH(NVV) = EVN( 8) VQM(NVV) = EVN( 9) VER(NVV) = EVN(10) VRT(NVV) = EVN(11) VID(NVV) = SID c ELSE c CALL BORT('UNKNOWN EVENT TYPE '//VAR) ENDIF GOTO 10 101 continue READ(LUNEV,'(a50)',END=20,err=20) string write(*,*) 'error found in 10-read: ', string goto 10 C OPEN INPUT AND OUTPUT PREPBUFR FILES C ------------------------------------ 20 CONTINUE write(*,*) 'open prepbufr file' CALL OPENBF(LUNIN,'IN ',LUNIN) CALL OPENBF(LUOUT,'OUT',LUNIN) write(*,*) 'finished open prepbufr' C NOW GO THROUGH A PREPBUFR FILE AND ENCODE THESE EVENTS C ------------------------------------------------------ DO WHILE(IREADMG(LUNIN,SUBSET,IDATE).EQ.0) DO WHILE(IREADSB(LUNIN).EQ.0) CALL UFBINT(LUNIN,HDR,10,1,IRET,HDSTR) C 11-2012 intel compiler reads binary 0s C DO I=1,8 C IF(ICHAR(SID(I:I)).EQ.00) SID(I:I) = ' ' C ENDDO C WNE 11/2012 III = LJUST(SID) CALL JSTCHR(SID,III) XOB = HDR(2) YOB = HDR(3) DHR = HDR(4) ELV = HDR(5) RTP = HDR(6) C INTERROGATE VARIABLES FOR EACH OB AT EACH LEVEL - STORE EVENTS C -------------------------------------------------------------- CALL UFBINT(LUNIN,OBS,10,255,NLEV,OBSTR) CALL UFBINT(LUNIN,QMS,10,255,NLEV,'PQM QQM TQM WQM') PEV = BMISS QEV = BMISS TEV = BMISS UEV = BMISS VEV = BMISS DO L=1,NLEV POB = OBS(1,L) QOB = OBS(2,L) TOB = OBS(3,L) UOB = OBS(4,L) VOB = OBS(5,L) CAT = OBS(6,L) ELV = OBS(7,L) QMP = QMS(1,L) QMQ = QMS(2,L) QMT = QMS(3,L) QMW = QMS(4,L) C CHECK FOR A SURFACE PRESSURE OB EVENT C ------------------------------------- IF(CAT.EQ.0) THEN NPOB = NPOB+1 DO N=1,NPP IF(SID.EQ.PID(N).AND.RTP.EQ.PRT(N).AND.PCD(N).GE.0) THEN IF(TOLX(XOB,PXX(N)) .AND. TOLX(YOB,PYY(N)) .AND. . TOL(DHR,PDH(N))) THEN IF(PCD(N).NE.0) THEN ! STORE 3DVAR QC REJECTION CODES PEV(1,L) = POB ! OBSERVED VALUE PEV(2,L) = 10 ! 3DVAR QC REJECT CODE PEV(3,L) = PCD(N) ! 3DVAR QC REASON CODE PEV(4,L) = ANLCODE ! 3DVAR PROGRAM CODE ENDIF DP = POO(N)/POB PEV(5,L) = PFC(N)/DP ! FIRST GUESS PEV(6,L) = PAN(N)/DP ! ANALYSIS PEV(7,L) = PER(N) ! OB ERROR PCD(N) = -1 GOTO 21 ENDIF ENDIF ENDDO ! NO EVENT - EITHER REJECTED BY R2 QC OR 3DVAR APRIORI FILTER IF(QMP.LT.4) THEN ! STORE FILTER EVENT CODES PEV(1,L) = POB ! OBSERVED VALUE PEV(2,L) = 9 ! 3DVAR FILTER CODE PEV(3,L) = 9 ! 3DVAR REASON CODE PEV(4,L) = ANLCODE ! 3DVAR PROGRAM CODE NP3DX = NP3DX+1 ELSE NPR2X = NPR2X+1 ENDIF ENDIF C CHECK FOR A SPECIFIC HUMIDITY OB EVENT C -------------------------------------- 21 IF(QOB.LT.BMISS) THEN NQOB = NQOB+1 DO N=1,NQQ IF(SID.EQ.QID(N).AND.RTP.EQ.QRT(N).AND.QCD(N).GE.0) THEN IF(TOLX(XOB,QXX(N)) .AND. TOLX(YOB,QYY(N)) .AND. . TOL(DHR,QDH(N)) .AND. TOL(QOB,QOO(N))) THEN IF(QCD(N).NE.0) THEN ! STORE 3DVAR QC REJECTION CODES QEV(1,L) = QOB ! OBSERVED VALUE QEV(2,L) = 10 ! 3DVAR QC REJECT CODE QEV(3,L) = QCD(N) ! 3DVAR QC REASON CODE QEV(4,L) = ANLCODE ! 3DVAR PROGRAM CODE ENDIF QEV(5,L) = QFC(N) ! FIRST GUESS QEV(6,L) = QAN(N) ! ANALYSIS QEV(7,L) = QER(N) ! OB ERROR QCD(N) = -1 GOTO 22 ENDIF ENDIF ENDDO ! NO EVENT - EITHER REJECTED BY R2 QC OR 3DVAR APRIORI FILTER IF(QMQ.LT.4) THEN ! STORE FILTER EVENT CODES QEV(1,L) = QOB ! OBSERVED VALUE QEV(2,L) = 9 ! 3DVAR FILTER MARK QEV(3,L) = 9 ! 3DVAR REASON CODE QEV(4,L) = ANLCODE ! 3DVAR PROGRAM CODE NQ3DX = NQ3DX+1 ELSE NQR2X = NQR2X+1 ENDIF ENDIF C CHECK FOR A TEMPERATURE OB EVENT C -------------------------------- 22 IF(TOB.LT.BMISS) THEN NTOB = NTOB+1 DO N=1,NTT IF(SID.EQ.TID(N).AND.RTP.EQ.TRT(N).AND.TCD(N).GE.0) THEN IF(TOLX(XOB,TXX(N)) .AND. TOLX(YOB,TYY(N)) .AND. . TOL(DHR,TDH(N)) .AND. TOL(POB,TPR(N))) THEN IF(TCD(N).NE.0) THEN ! STORE 3DVAR QC REJECTION CODES TEV(1,L) = TOB ! OBSERVED VALUE TEV(2,L) = 10 ! 3DVAR QC REJECT CODE TEV(3,L) = TCD(N) ! 3DVAR QC REASON CODE TEV(4,L) = ANLCODE ! 3DVAR PROGRAM CODE ENDIF TEV(5,L) = TFC(N) ! FIRST GUESS TEV(6,L) = TAN(N) ! ANALYSIS TEV(7,L) = TER(N) ! OB ERROR TCD(N) = -1 GOTO 23 ENDIF ENDIF ENDDO ! NO EVENT - EITHER REJECTED BY R2 QC OR 3DVAR APRIORI FILTER IF(QMT.LT.4) THEN ! STORE FILTER EVENT CODES TEV(1,L) = TOB ! OBSERVED VALUE TEV(2,L) = 9 ! 3DVAR FILTER MARK TEV(3,L) = 9 ! 3DVAR REASON CODE TEV(4,L) = ANLCODE ! 3DVAR PROGRAM CODE NT3DX = NT3DX+1 ELSE NTR2X = NTR2X+1 ENDIF ENDIF C CHECK FOR A U COMPONENT OB C -------------------------- 23 IF(UOB.LT.BMISS) THEN NUOB = NUOB+1 DO N=1,NUU IF(SID.EQ.UID(N).AND.RTP.EQ.URT(N).AND.UCD(N).GE.0) THEN IF(TOLX(XOB,UXX(N)) .AND. TOLX(YOB,UYY(N)) .AND. . TOL(DHR,UDH(N)) .AND. TOL(POB,UPR(N)) .AND. . TOL(UOB,UOO(N))) THEN IF(UCD(N).NE.0) THEN ! STORE 3DVAR QC REJECTION CODES UEV(1,L) = UOB ! OBSERVED VALUE UEV(2,L) = 10 ! 3DVAR QC REJECT CODE UEV(3,L) = UCD(N) ! 3DVAR QC REASON CODE UEV(4,L) = ANLCODE ! 3DVAR PROGRAM CODE ENDIF UEV(5,L) = UFC(N) ! FIRST GUESS UEV(6,L) = UAN(N) ! ANALYSIS UEV(7,L) = UER(N) ! OB ERROR UCD(N) = -1 GOTO 24 ENDIF ENDIF ENDDO ! NO EVENT - EITHER REJECTED BY R2 QC OR 3DVAR APRIORI FILTER IF(QMW.LT.4) THEN ! STORE FILTER EVENT CODES UEV(1,L) = UOB ! OBSERVED VALUE UEV(2,L) = 9 ! 3DVAR FILTER MARK UEV(3,L) = 9 ! 3DVAR REASON CODE UEV(4,L) = ANLCODE ! 3DVAR PROGRAM CODE NU3DX = NU3DX+1 ELSE NUR2X = NUR2X+1 ENDIF ENDIF C CHECK FOR A V COMPONENT OB C -------------------------- 24 IF(VOB.LT.BMISS) THEN NVOB = NVOB+1 DO N=1,NVV IF(SID.EQ.VID(N).AND.RTP.EQ.VRT(N).AND.VCD(N).GE.0) THEN IF(TOLX(XOB,VXX(N)) .AND. TOLX(YOB,VYY(N)) .AND. . TOL(DHR,VDH(N)) .AND. TOL(POB,VPR(N)) .AND. . TOL(VOB,VOO(N))) THEN IF(VCD(N).NE.0) THEN ! STORE 3DVAR QC REJECTION CODES VEV(1,L) = VOB ! OBSERVED SPECIFIC HUMIDITY VEV(2,L) = 10 ! 3DVAR QC REJECT CODE VEV(3,L) = VCD(N) ! 3DVAR QC REASON CODE VEV(4,L) = ANLCODE ! 3DVAR PROGRAM CODE ENDIF VEV(5,L) = VFC(N) ! FIRST GUESS VEV(6,L) = VAN(N) ! ANALYSIS VEV(7,L) = VER(N) ! OB ERROR VCD(N) = -1 GOTO 25 ENDIF ENDIF ENDDO ! NO EVENT - EITHER REJECTED BY R2 QC OR 3DVAR APRIORI FILTER IF(QMW.LT.4) THEN ! STORE FILTER EVENT CODES VEV(1,L) = VOB ! OBSERVED VALUE VEV(2,L) = 9 ! 3DVAR FILTER MARK VEV(3,L) = 9 ! 3DVAR REASON CODE VEV(4,L) = ANLCODE ! 3DVAR PROGRAM CODE NV3DX = NV3DX+1 ELSE NVR2X = NVR2X+1 ENDIF ENDIF 25 CONTINUE C MAKE SURE UEV AND VEV ARE COORDINATED PROPERLY C ---------------------------------------------- IF(UEV(1,L).LT.BMISS .AND. VEV(1,L).GE.BMISS) THEN VEV(1,L) = VOB ELSEIF(UEV(1,L).GE.BMISS .AND. VEV(1,L).LT.BMISS) THEN UEV(1,L) = UOB UEV(2,L) = VEV(2,L) UEV(3,L) = VEV(3,L) UEV(4,L) = VEV(4,L) ENDIF ENDDO ! END OF REPORT LEVEL LOOP C WRITE THIS REPORT WITH EVENTS INTO THE OUTPUT FILE C -------------------------------------------------- CALL OPENMB(LUOUT,SUBSET,IDATE) CALL UFBCPY(LUNIN,LUOUT) CALL UFBINT(LUOUT,PEV,10,NLEV,IRET,PVSTR) CALL UFBINT(LUOUT,QEV,10,NLEV,IRET,QVSTR) CALL UFBINT(LUOUT,TEV,10,NLEV,IRET,TVSTR) CALL UFBINT(LUOUT,UEV,10,NLEV,IRET,UVSTR) CALL UFBINT(LUOUT,VEV,10,NLEV,IRET,VVSTR) CALL WRITSB(LUOUT) C END OF READ/WRITE LOOP C ---------------------- ENDDO ! END OF READSB LOOP ENDDO ! END OF READMG LOOP CALL CLOSBF(LUOUT) C WRITE A SUMMARY OF EVENTS FOR EACH VARIABLE C ------------------------------------------- IPSUM = NPOB-(NPP+NP3DX+NPR2X) IQSUM = NQOB-(NQQ+NQ3DX+NQR2X) ITSUM = NTOB-(NTT+NT3DX+NTR2X) IUSUM = NUOB-(NUU+NU3DX+NUR2X) IVSUM = NVOB-(NVV+NV3DX+NVR2X) print'(1x, 5a10)','totobs','totevn','3dfltr','r2fltr','cksum' PRINT'("p",5I10)',NPOB,NPP,NP3DX,NPR2X,IPSUM PRINT'("q",5I10)',NQOB,NQQ,NQ3DX,NQR2X,IQSUM PRINT'("t",5I10)',NTOB,NTT,NT3DX,NTR2X,ITSUM PRINT'("u",5I10)',NUOB,NUU,NU3DX,NUR2X,IUSUM PRINT'("v",5I10)',NVOB,NVV,NV3DX,NVR2X,IVSUM !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! this block for diagnosing checksum errors (example given for t, similar for p,q,u,v) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! do n=1,nuu C if(tcd(n).ge.0) then if(ucd(n).ge.0) then C print*,tid(n),txx(n),tyy(n),tdh(n),trt(n),tpr(n) print*,'n=',n,uid(n),uxx(n),uyy(n),udh(n),urt(n),upr(n) endif enddo !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! if(ipsum.ne.0.or.iqsum.ne.0.or.itsum.ne.0.or. . iusum.ne.0.or.ivsum.ne.0) then print*,'>>>>>>>>>>>rrvent checksum error<<<<<<<<<<' stop 1 endif if (nuu.ne.nvv) then print*,'>>>>>>>>>>>rrvent uv event error<<<<<<<<<<' stop 1 endif STOP END