SUBROUTINE ANNATV(IUNIT,IM,JM,IGRID,KPDS,KGDS,OPCP,MASK) ! ! M BALDWIN 2/6/97 ! MODIFIED VERSION OF PROGRAM PCPANA ! 1. READ IN RAINFALL DATA - DO QC ! 2. DO A BOX AVERAGE ANALYSIS TO THE GRID DEFINED BY KGDS ! ! PARAMETER (MXSIZE=150000,NPTS=20000) DIMENSION OPCP(IM*JM),MASK(IM*JM) DIMENSION ICOUNT(MXSIZE),XPTS(NPTS),YPTS(NPTS) DIMENSION SROT(NPTS),CROT(NPTS) DIMENSION OLAT(NPTS),OLON(NPTS),PCP(NPTS) INTEGER KPDS(25),KGDS(20),MNTH(12) CHARACTER*9 STATID, STID(NPTS) CHARACTER fn2*80 CHARACTER*3 month(12) CHARACTER END*5,DAYOW*5,STI*5 ! data month/'jan','feb','mar','apr','may','jun','jul'& ,'aug','sep','oct','nov','dec'/ DATA MNTH/31,28,31,30,31,30,31,31,30,31,30,31/ ! IMJM=IM*JM RADDEG=3.14159265/180.0 ! DO 100 J=1,IMJM ICOUNT(J) = 0 MASK(J) = 0 OPCP(J) = 0.0 100 CONTINUE ! ! THESE ARE THE GEODETIC COORDS TO DEFINE THE GRID, FROM KGDS ! IF (KGDS(1).EQ.0) THEN ALAT1=KGDS(4)*1.E-3 ALON1=KGDS(5)*1.E-3 ALAT2=KGDS(7)*1.E-3 ALON2=KGDS(8)*1.E-3 DX=(ALON2-ALON1)/(IM-1) DY=(ALAT2-ALAT1)/(JM-1) ELSEIF (KGDS(1).EQ.1) THEN ALAT1=KGDS(4)*1.E-3 ALON1=KGDS(5)*1.E-3 ALAT2=KGDS(7)*1.E-3 ALON2=KGDS(8)*1.E-3 ALATIN=KGDS(9)*1.E-3 DX=KGDS(12) ELSEIF (KGDS(1).EQ.3) THEN ALAT1=KGDS(4)*1.E-3 ALON1=KGDS(5)*1.E-3 ELONV=KGDS(7)*1.E-3 DX=KGDS(8) ALATAN=KGDS(12)*1.E-3 ELSEIF (KGDS(1).EQ.4) THEN ALAT1=KGDS(4)*1.E-3 ALON1=KGDS(5)*1.E-3 ALAT2=KGDS(7)*1.E-3 ALON2=KGDS(8)*1.E-3 NLATC=KGDS(10) ELSEIF (KGDS(1).EQ.5) THEN ALAT1=KGDS(4)*1.E-3 ALON1=KGDS(5)*1.E-3 ELONV=KGDS(7)*1.E-3 DX=KGDS(8) ELSEIF (KGDS(1).EQ.201) THEN IF (IGRID.EQ.90) THEN CLAT0=52. CLON0=111. DLATD=14./26. DLOND=15./26. ELSEIF (IGRID.EQ.92) THEN CLAT0=50. CLON0=107. DLATD=8./39. DLOND=2./9. ELSEIF (IGRID.EQ.94) THEN CLAT0=41. CLON0=97. DLATD=5./27. DLOND=7./36. ELSEIF (IGRID.EQ.96) THEN CLAT0=50. CLON0=111. DLATD=5./16. DLOND=41./124. ENDIF WBD=-(IM-1)*DLOND SBD=-((JM-1)/2)*DLATD ALAT1=KGDS(4)*1.E-3 ALON1=KGDS(5)*1.E-3 ENDIF ! ! READ PRECIP OBS ! K=0 READ(IUNIT,7773) IHR,IYR,IMN,IDA 330 READ(IUNIT,7778,END=1005,ERR=980) ALAT,ALON,PCPK,STI K=K+1 STID(K)=STI OLAT(K)=ALAT OLON(K)=-ALON PCP(K)=PCPK*25.4 ! CHANGE UNITS TO MM 980 CONTINUE GOTO 330 7773 FORMAT(29X,I2,5X,I4,2I2) 7778 FORMAT(F6.2,F9.2,F7.2,1X,A8) 1005 PRINT 101,K 101 FORMAT(1X,'FINISHED READING DATA',I8,' OBSERVATIONS') NOBSK=K ! ! CALC AVERAGE AND STANDARD DEVIATION FOR RFC DATA ! MAXIMUM ALLOWABLE OBSERVATION WILL BE AVG + 3*STDV ! NEED TO CHANGE THIS MAYBE? ! SUMK=0. SSQSK=0. DO K=1,NOBSK FIPCP = PCP(K) SUMK = SUMK + FIPCP SSQSK =SSQSK + FIPCP ** 2 ENDDO XN = FLOAT(NOBSK) IF (XN.GT.0.) THEN AVG = SUMK / XN ASUMSQ = SSQSK / XN YX = ASUMSQ - AVG ** 2 QCMAX = AVG + 6.* SQRT( YX ) ENDIF WRITE(6,*) 'OBS DATA avg variance qcmax',avg,yx,qcmax ! ! NOW READ THROUGH DATA AND TOSS OBS > AVG + 6 STD DEV ! THEN ANALYZE ! ! ! GET X,Y COORDS OF OBS LAT/LON ! CALL GDSWIZ(KGDS,-1,NOBSK,-9999.,XPTS,YPTS,OLON,OLAT,NRET,& 0,CROT,SROT) ! ! SUBROUTINE GDSWIZ(KGDS,IOPT,NPTS,FILL,XPTS,YPTS,RLON,RLAT,NRET,& ! LROT,CROT,SROT) ! INPUT ARGUMENT LIST: ! KGDS - INTEGER (200) GDS PARAMETERS AS DECODED BY W3FI63 ! IOPT - INTEGER OPTION FLAG ! ( 0 TO COMPUTE EARTH COORDS OF ALL THE GRID POINTS) ! (+1 TO COMPUTE EARTH COORDS OF SELECTED GRID COORDS) ! (-1 TO COMPUTE GRID COORDS OF SELECTED EARTH COORDS) ! NPTS - INTEGER MAXIMUM NUMBER OF COORDINATES ! FILL - REAL FILL VALUE TO SET INVALID OUTPUT DATA ! (MUST BE IMPOSSIBLE VALUE; SUGGESTED VALUE: -9999.) ! XPTS - REAL (NPTS) GRID X POINT COORDINATES IF IOPT>0 ! YPTS - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT>0 ! RLON - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT<0 ! (ACCEPTABLE RANGE: -360. TO 360.) ! RLAT - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT<0 ! (ACCEPTABLE RANGE: -90. TO 90.) ! LROT - INTEGER FLAG TO RETURN VECTOR ROTATIONS IF 1 ! ! OUTPUT ARGUMENT LIST: ! XPTS - REAL (NPTS) GRID X POINT COORDINATES IF IOPT<=0 ! YPTS - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT<=0 ! RLON - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT>=0 ! RLAT - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT>=0 ! NRET - INTEGER NUMBER OF VALID POINTS COMPUTED ! (-1 IF PROJECTION UNRECOGNIZED) ! CROT - REAL (NPTS) CLOCKWISE VECTOR ROTATION COSINES IF LROT=1 ! SROT - REAL (NPTS) CLOCKWISE VECTOR ROTATION SINES IF LROT=1 ! (UGRID=CROT*UEARTH-SROT*VEARTH; ! VGRID=SROT*UEARTH+CROT*VEARTH) ! DO K=1,NOBSK FIPCP = PCP(K) IF (FIPCP.GT.QCMAX) THEN WRITE(6,3301) STID(K),PCP(K),QCMAX ELSE ! ! FIND KINDEX OF OBSERVATION ! II = NINT(XPTS(K)) JJ = NINT(YPTS(K)) IF (II.GE.1.AND.II.LE.IM.AND.JJ.GE.1.AND.JJ.LE.JM) THEN KINDEX=(JJ-1)*IM+II ENDIF ! IF (KGDS(1).EQ.201) THEN ! ALN = -OLON(K) ! APH = OLAT(K) ! CALL ETA2IJ (APH,ALN,IM,CLON0,CLAT0,DLOND,DLATD, ! & WBD,SBD,KINDEX) ! ENDIF ! WRITE (6,*) 'II = ',ii,' JJ = ',jj,' KINDEX = ',kindex IF (KINDEX.GT.0.AND.KINDEX.LE.IMJM) THEN OPCP(KINDEX)=OPCP(KINDEX)+FIPCP ICOUNT(KINDEX)=ICOUNT(KINDEX)+1 ENDIF ENDIF ENDDO 3301 FORMAT(' QC TOSS STATION ',A5,' AMOUNT ',F6.2,' QCMAX ',F8.2) ! ! DIVIDE BY SUM OF WEIGHTS ! AMAX=0.0 DO K=1,IMJM IF (ICOUNT(K).GT.0) THEN RCOUNT=1./FLOAT(ICOUNT(K)) OPCP(K)=OPCP(K)*RCOUNT AMAX=AMAX1(OPCP(K),AMAX) MASK(K)=1 ELSE MASK(K)=0 ENDIF ENDDO WRITE(6,*) 'MAX ANALYSIS VALUE = ',AMAX ! ! SET UP DATE INFO ! IF(MOD(IYR,4).EQ.0) THEN MNTH(2) = 29 END IF IYR1=IYR IHR1=IHR IDA1=IDA-1 IMN1=IMN IF (IHR1.LT.0) THEN IHR1=IHR1+24 IDA1=IDA1-1 ENDIF IF (IDA1.LE.0) THEN IF(IMN1.EQ.1) THEN IMN1 = 12 IYR1 = IYR1 -1 ELSE IMN1 = IMN1 -1 ENDIF IDA1 = MNTH(IMN1) + IDA1 ENDIF ICENT=(IYR1-1)/100 + 1 ! ! SET UP PDS ! KPDS(1)=7 ! ID OF CENTER KPDS(2)=154 ! GENERATING PROCESS I MADE UP KPDS(3)=IGRID ! GRID NUMBER KPDS(4)=128+64 ! TABLE 1 FLAG 192 MEANS BMS IS HERE KPDS(5)=61 ! PARAMETER KPDS(6)=1 ! TYPE OF LEVEL KPDS(7)=0 ! VALUE OF LEVEL KPDS(8)=IYR1-(ICENT-1)*100 ! YEAR (START TIME) KPDS(9)=IMN1 ! MONTH KPDS(10)=IDA1 ! DAY KPDS(11)=IHR1 ! HOUR KPDS(12)=0 ! MINUTE KPDS(13)=1 ! TIME UNIT KPDS(14)=0 ! P1 KPDS(15)=24 ! P2 KPDS(16)=4 ! TIME RANGE INDICATOR KPDS(17)=0 ! NO INCLUDED IN AVERAGE KPDS(18)=1 ! GRIB VERSION NO KPDS(19)=2 ! PARAMETER TABLE VERSION NO KPDS(20)=0 ! NO MISSING KPDS(21)=ICENT ! CENTURY KPDS(22)=2 ! DECIMAL SCALE FACTOR KPDS(23)=4 ! SUB CENTER NUM KPDS(24)=0 ! ENSEMBLE STUFF KPDS(25)=0 ! RESERVED ! RETURN END