c updated from operations 4/2001 C$$$ MAIN PROGRAM DOCUMENTATION BLOCK C . . . . C MAIN PROGRAM: OIQCTEMP OPTIMAL INTERPOLATION QUALITY CONTROL C PRGMMR: WOOLLEN ORG: W/NMC2 DATE: 90-11-06 C C ABSTRACT: QUALITY CONTROLS OBSERVATIONS TO BE USED IN THE GLOBAL C ANALYSIS-FORECAST SYSTEM. OBSERVATIONAL DATASET GENERATED BY C THE GESTEMP PREPROCESSOR IS INPUT TO OIQC WHICH THEN PERFORMS C A "COMPLEX" QUALITY CONTROL (FROM GANDIN,1985) EMPLOYING C MULTIVARIATE OI TO CHECK ALL OBSERVATIONS AGAINST NEARBY C NEIGHBORS. SEVERAL ITERATIONS OF THE BASIC SCHEME ARE REQUIRED C TO COMPLETE THE PROCESS. DURING EACH ITERATION ALL OBSERVATIONS C ARE SUBJETED TO FIVE INTERPOLATION CHECKS. CHECK 1 INTERPOLATES C COMPARISON VALUES FROM NEARBY TEMPERATURE DATA. CHECK 2 AND CHECK C 3 INTERPOLATE COMPARISON VALUES FROM NEARBY ZONAL AND LONGITUDINAL C WIND COMPONENTS RESPECTIVELY. CHECK 4 IS AN INTERPOLATION C FROM THE NEAREST OBS IN THE SAME PROFILE AS THE ONE BEING CHECKED. C A WEIGHTED COMBINATION OF THE OUTCOME OF THESE FOUR CHECKS IS C PRODUCED ON THE BASIS OF FACTORS WHICH DETERMINE WHICH C OF THE INTERPOLATIONS (IF ANY) MIGHT BE EXPECTED TO BE MORE C RELEVANT THAN THE OTHERS; FOR EXAMPLE THE NUMBER AND LOCATION C IN SPACE AND TIME OF OBSERVATIONS USED IN A GIVEN CHECK C AS WELL AS RELATIVE QUALITY OF THE OBSERVATIONS INVOLVED IS USED TO C INDICATE WHICH OF THE OUTCOMES SHOULD BE RELIED MORE HEAVILY UPON C TO MAKE A DECISION IN A PARTICULAR CASE. CHECK 5 INTERPOLATES C A VALUE FROM THE MULTIVARIATE SET OF OBS USED IN CHECKS 1-3, C WHICH IS USED TO MEASURE HOW THE ANALYSIS WOULD DRAW (OR NOT DRAW) C TO THE OB BEING CHECKED. THE TOLERANCE ALLOWED FOR DEVIATION FROM C THE EXPECTED OUTCOME OF THE INTERPOLATION CHECKS IS PROPORTIONAL C TO A MEASURE OF THE DRAW COMPUTED USING THE RESULT OF CHECK 5. C DURING EACH ITERATION EACH OBSEVATION EITHER PASSES OR FAILS THE C COMPLEX OF CHECKS JUST DESCRIBED. A "PASS" MEANS THAT OB MAY BE C USED FOR CHECKING OBS DURING THE SUBSEQUENT ITERATION. A "FAIL" C INDICATES THE OPPOSITE. THE SYSTEM IS ITERATED UNTIL IDENTICAL C RESULTS IN TERMS OF OBS THAT PASS AND FAIL ARE OBTAINED IN TWO C CONSEQUETIVE ITERATIONS UP TO FOUR COMPLETE ITERATIONS AT WHICH C POINT A FINAL ARBITRATION PROCEDURE IS INVOKED TO RESOLVE A (SMALL) C NUMBER OF AMBIGUOUS CASES WHICH ARE PREVENTING CONVERGENCE OF THE C SYSTEM. C C PROGRAM HISTORY LOG: C 90-11-06 J. WOOLLEN ORIGINAL VERSION FOR IMPLEMENTATION C 90-04-23 J. WOOLLEN VERSION FOR IMPLEMENTATION W/ SPECTRAL OI C C USAGE: C INPUT FILES: C FORT.11 - PREPDA - OBSERVATION FILE C FORT.12 - PREPQC - OBSERVATION INCREMENT FILE FROM GEST C FORT.14 - OBERRS - OBSERVATION ERROR FILE C FORT.16 - NMCDATE - NMC DATE FILE C FORT.18 - OBPRT.IN - LIST OF IPOINTS FOR DETAILED PRINT C FORT.20 - TOLLS.IN - TOLERANCE LEVELS (FOR EXPERIMENTATION) C C OUTPUT FILES: C FORT.60 - OBCNT.OUT - SUMMARY OF OIQC RUN C FORT.61 - TOSS.SFZ - DETAILED TOSSLIST FOR SURFACE HEIGHTS C FORT.62 - TOSS.OUT - DETAILED TOSSLIST FOR TEMPS AND WINDS C FORT.63 - TOSS.SAT - DETAILED TOSSLIST FOR SAT TEMPS C FORT.64 - TOSS.SMI - DETAILED TOSSLIST FOR SSMI WIND SPEEDS C FORT.65 - TOSSLIST - OFFICIAL OPERATIONAL TOSSLIST C FORT.66 - PSFILE.OUT - LIST OF OBS TOSSED FOR BEING UNDERGROUND C FORT.70 - PREPQM - OBSERVATIONAL DATA FOR OI INGEST C FORT.71 - 80 - RESERVED FOR DIAGNOSTIC TESTING C FORT.81 - OBOGRAM.OUT - SUMMARY OF QC RESULTS (ASCII) C FORT.82 - OBOGRAM.BIN - SUMMARY OF QC RESULTS (BINARY) C C SUBPROGRAMS CALLED: (*AUTOTASKED) C UNIQUE: - SETCON,RDOER,STORE,STOERR,PSFILE,MAKEAMAP, C SEARCH,CHDIST,QCLOOP,QCPTS,SATPTS, C *SELDAT,*QCOI,*COORS,*DRCTSL,*VSOLVE,*OUTPUT,*SATPUT, C OBOGRAM,TOSSLIST,PREPQM C C LIBRARY: C COMMON - W3LIB C C EXIT STATES: C COND = 0 - SUCCESSFUL RUN C C IF ANY DISASTORIUS PROBLEM SITUATIONS ARE DETECTED THE PROGRAM C WILL ABORT WITH A MESSAGE DETAILING THE CALAMITY. C C C REMARKS: C C ATTRIBUTES: C LANGUAGE: CRAY FORTRAN C MACHINE: CRAY C C$$$ C----------------------------------------------------------------------- CTREE GETFCER CTREE PILNLNP CTREE SATPTS CTREE SATPUT C----------------------------------------------------------------------- PROGRAM OIQCTEMP COMMON /OBS/ NREP,NLEV,SID( 200000 ),INOB( 200000 ),IMAP(360,181) COMMON /CKLIST / NDD,INDD( 300000 ) COMMON /TOSSUNIT/ ITUNIT CHARACTER*3 FOC EXTERNAL QCPTS,SATPTS DATA LUDAT /11/ DATA LUBFI /14/ DATA LUQCE /17/ DATA LUBFQ /70/ C----------------------------------------------------------------------- C UNIFIED QC ANALYSIS SETUP PROCEDURE IS AS FOLLOWS: C ------------------------------------------------- C 0) W3LOG - WRITE AN OPERATIONS START MESSAGE C 1) SETCON - SET CONSTANTS AND PARAMETERS C 2) STORE - READ IN FERR (OB/INCREMENT) FILE AND STORE IT C 3) RDOBER - READ IN THE OBSERVATION ERRORS C 4) STOERR - STORE THE OB AND FC ERRORS WITH EACH REPORT LEVEL C 5) MAKEAMAP - MAKE AN OB POINTER LIST SORTED BY LAT-LON C----------------------------------------------------------------------- C CALL W3LOG('$S','$M') CALL SETCON CALL STORE(LUDAT,LUBFI) CALL RDOBER (LUQCE) CALL STOERR CALL MAKEAMAP C OIQC CHECKS OF SURFACE HEIGHT/PRESSURE DATA C ------------------------------------------- NDD = 0 DO IREP=1,NREP USEFLG = ABS(FQL(IREP)) REPOFF = WQL(IREP,-USEFLG) IF(FOC(IREP).EQ.'SFC') THEN SET1 = WNL(IREP,1) ILEV = FLV(IREP)+1 IF(FOE(ILEV).GT.0.) THEN NDD = NDD+1 INDD(NDD) = ILEV REPON = WQL(IREP,USEFLG) SETZ = WLV(IREP,ILEV) ENDIF ENDIF ENDDO IF(NDD.GT.0) THEN WRITE(6 ,*)' CHECKING SURFACE HEIGHT DATA ',NDD,' OBS' WRITE(60,*)' CHECKING SURFACE HEIGHT DATA ',NDD,' OBS' ITUNIT = 61 CALL SEARCH CALL QCLOOP(QCPTS,4) DO N=1,NDD IREP = FHD(INDD(N)) SETT = WLV(IREP,INDD(N)-1) ENDDO ENDIF C OIQC CHECKS OF UPA AND SFC TEMP AND WIND DATA C ------------------------------------------- NDD = 0 DO IREP=1,NREP USEFLG = ABS(FQL(IREP)) REPOFF = WQL(IREP,-USEFLG) IF(FOC(IREP).EQ.'UPA' .OR. FOC(IREP).EQ.'SFC') THEN REPON = WQL(IREP,USEFLG) ILEV = FLV(IREP) KLEV = ILEV+FNL(IREP)-1 DO JLEV=ILEV,KLEV IF(FOE(JLEV).GT.0.) THEN NDD = NDD+1 INDD(NDD) = JLEV ENDIF ENDDO ENDIF ENDDO IF(NDD.GT.0) THEN WRITE(6 ,*)' CHECKING TEMP AND WIND DATA ',NDD,' OBS' WRITE(60,*)' CHECKING TEMP AND WIND DATA ',NDD,' OBS' ITUNIT = 62 CALL SEARCH CALL QCLOOP(QCPTS,4) ENDIF C OIQC CHECKS OF SATEM DATA C ------------------------- NDD = 0 DO IREP=1,NREP USEFLG = ABS(FQL(IREP)) IF(FOC(IREP).EQ.'SAT') THEN REPON = WQL(IREP,1.) ILEV = FLV(IREP) KLEV = ILEV+FNL(IREP)-1 DO JLEV=ILEV,KLEV IF(FOE(JLEV).GT.0.) THEN NDD = NDD+1 INDD(NDD) = JLEV ENDIF ENDDO ELSE REPOFF = WQL(IREP,-USEFLG) ENDIF ENDDO IF(NDD.GT.0) THEN WRITE(6 ,*)' CHECKING SAT TEMP DATA ',NDD,' OBS' WRITE(60,*)' CHECKING SAT TEMP DATA ',NDD,' OBS' ITUNIT = 63 CALL SEARCH CALL QCLOOP(SATPTS,1) ENDIF C DIRECTION ASSIGNMENT AND OIQC CHECKS FOR SSMI DATA C -------------------------------------------------- NDD = 0 DO IREP=1,NREP USEFLG = ABS(FQL(IREP)) IF(FOC(IREP).EQ.'SMI') THEN ILEV = FLV(IREP) IF(FOE(ILEV).GT.0.) THEN NDD = NDD+1 INDD(NDD) = ILEV REPON = WQL(IREP,USEFLG) ENDIF ELSE IF(FOC(IREP).EQ.'SAT') THEN REPOFF = WQL(IREP,-USEFLG) ELSE REPOFF = WQL(IREP,0.) ENDIF ENDDO IF(NDD.GT.0) THEN WRITE(6 ,*)' CHECKING SSMI DATA ',NDD,' OBS' WRITE(60,*)' CHECKING SSMI DATA ',NDD,' OBS' ITUNIT = 64 CALL SSMIPTS CALL SEARCH CALL QCLOOP(QCPTS,4) ENDIF C COMPUTE STATISTICS, WRITE LISTS, UPDATE QUALITY MARKS C ----------------------------------------------------- CALL OBOGRAM (11,81,82) CALL TOSSLIST (65) CALL PREPQM (LUBFI,LUBFQ) C WRITE AN OPERATIONS END MESSAGE C ------------------------------- C CALL W3LOG('$E') STOP END C$$$ SUBPROGRAM DOCUMENTATION BLOCK C C SUBPROGRAM: CHDIST C PRGMMR: WOOLLEN ORG: NMC22 DATE: 90-11-06 C C ABSTRACT: C COMPUTES CHORD LENGTH DISTANCE FROM A FIXED C POINT TO AN ARRAY OF POINTS USING THE FORMULA: C S**2/2 = 1 - COS(Y1-Y2) + COS(Y1)*COS(Y2)*(1-COS(X1-X2)). C ALSO THE DIRECTION OF EACH POINT IS COMPUTED WITH RESPECT TO C THE FIXED POINT AND RETURNED AS WELL. C C PROGRAM HISTORY LOG: C 90-11-06 J. WOOLLEN C C USAGE: C INPUT ARGUMENTS: C X1 - X-COORDINATE (LONGITUDE) OF FIXED POINT C Y1 - Y-COORDINATE (LATITUDE ) OF FIXED POINT C X2 - X-COORDINATES (LONGITUDE) FOR SET OF POINTS C Y2 - Y-COORDINATES (LATITUDE ) FOR SET OF POINTS C NP - NUMBER OF OUTPUTS REQUESTED C C OUTPUT ARGUMENTS: C DIST - CHORD LENGTH DISTANCES FROM FIXED POINT (KM) C DIRN - DIRECTION FROM FIXED POINT (DEG) C C SUBPROGRAMS CALLED: NONE C C EXIT STATES: NONE C C REMARKS: C C C ATTRIBUTES: C LANGUAGE: CRAY FORTRAN C MACHINE: CRAY C C$$$ SUBROUTINE CHDIST(X1,Y1,X2,Y2,DIST,DIRN,NP) DIMENSION X2(NP),Y2(NP),DIST(NP),DIRN(NP) DATA PI180/.0174532 /,RADE/6371./ C---------------------------------------------------------------------- C---------------------------------------------------------------------- IF(NP.EQ.0) RETURN C COMPUTE THE DISTANCE C -------------------- DO 10 I=1,NP COSY1 = COS(Y1*PI180) COSY2 = COS(Y2(I)*PI180) COSDX = COS((X1-X2(I))*PI180) COSDY = COS((Y1-Y2(I))*PI180) S = 1.0-COSDY+COSY1*COSY2*(1.0-COSDX) S = SQRT(2.*S) IF(S.LE..002) S = 0. DIST(I) = S*RADE 10 CONTINUE C COMPUTE DIRECTIONS C ------------------ DO 20 I=1,NP DX = (X2(I)-X1)*COS(Y1) DY = (Y2(I)-Y1) IF(DX.GT.0.) THEN DIRN(I) = 0. + ATAN(DY/DX)/PI180 ELSE IF(DX.LT.0.) THEN DIRN(I) = 180. + ATAN(DY/DX)/PI180 ELSE IF(DX.EQ.0.) THEN DIRN(I) = SIGN(90.,DY) ENDIF IF(DIRN(I).LT.0.) DIRN(I) = DIRN(I) + 360. 20 CONTINUE RETURN END C$$$ SUBPROGRAM DOCUMENTATION BLOCK C C SUBPROGRAM: COORS C PRGMMR: WOOLLEN ORG: NMC22 DATE: 90-11-06 C C ABSTRACT: COMPUTES CHORD LENGTH DISTANCE FOR A MATRIX OF C LOCATIONS <(X,Y)I,(X,Y)J> USING THE NORMAL ANGLE FORMULA: C S**2/2 = 1 - COS(Y1-Y2) + COS(Y1)*COS(Y2)*(1-COS(X1-X2)) C DERIVATIVES OF DISTANCE WITH RESPECT TO LAT AND LON ARE ALSO C COMPUTED AND COMBINED WITH APPROPRIATE CORRELATION FUNCTIONS C AND DERIVATIVES WITH RESPECT TO DISTANCE TO FORM THE MULTI- C VARIATE CORRELATIONS FOR THE HEIGHT WIND ANALYSIS. C C PROGRAM HISTORY LOG: C 90-11-06 J. WOOLLEN C C USAGE: C INPUT ARGUMENTS: C NP - NUMBER OF OUTPUTS REQUESTED C CH - VECTOR OF CONSTANTS FOR THE GAUSSIAN LENGTH SCALE C X1 - X-COORDINATES (LONGITUDE) FOR POINT 1 C Y1 - Y-COORDINATES (LATITUDE ) FOR POINT 1 C K1 - VARIABLE TYPE AT POINT 1 (I.E. 1-Z , 2-U , 3-V) C X2 - X-COORDINATES (LONGITUDE) FOR POINT 2 C Y2 - Y-COORDINATES (LATITUDE ) FOR POINT 2 C K2 - VARIABLE TYPE AT POINT 2 (I.E. 1-Z , 2-U , 3-V) C FL - DECOUPLING FACTOR C TZ - HEIGHT-TEMPERATURE FACTOR C C OUTPUT ARGUMENTS: C D - DISTANCES BETWEEN ARRAYS OF POINTS IN RADIANS C C1 - CORRELATION COEFFICIENTS AS SPECIFIED C C2 - CORRELATION COEFFICIENTS AS SPECIFIED C C SUBPROGRAMS CALLED: NONE C C EXIT STATES: NONE C C REMARKS: THIS PROGRAM COMPUTES A GAUSSIAN CORRELATION WITH C LENGTH SCALE GIVEN BY THE CHTWV INPUT ARGUMENT. C C C ATTRIBUTES: C LANGUAGE: CRAY FORTRAN C MACHINE: CRAY C C$$$ SUBROUTINE COOR1(NP,CH,FL,TZ,X1,Y1,K1,X2,Y2,K2,D,C1,C2) DIMENSION CH(NP),FL(NP),TZ(NP) DIMENSION X1(NP),Y1(NP),K1(NP) DIMENSION X2(NP),Y2(NP),K2(NP) DIMENSION D (NP),C1(NP),C2(NP) LOGICAL ZZ,ZU,ZV,UZ,UU,UV,VZ,VU,VV DATA PI180 /.0174532 /, SMIN/.002/ C---------------------------------------------------------------------- C---------------------------------------------------------------------- IF(NP.EQ.0) RETURN C LOOP OVER SET OF INPUT POINTS C ----------------------------- DO 20 I=1,NP C1(I) = 0. C2(I) = 0. C DETERMINE CORRELATION TYPE C -------------------------- ZZ = K1(I).EQ.1 .AND. K2(I).EQ.1 ZU = K1(I).EQ.1 .AND. K2(I).EQ.2 ZV = K1(I).EQ.1 .AND. K2(I).EQ.3 UZ = K1(I).EQ.2 .AND. K2(I).EQ.1 UU = K1(I).EQ.2 .AND. K2(I).EQ.2 UV = K1(I).EQ.2 .AND. K2(I).EQ.3 VZ = K1(I).EQ.3 .AND. K2(I).EQ.1 VU = K1(I).EQ.3 .AND. K2(I).EQ.2 VV = K1(I).EQ.3 .AND. K2(I).EQ.3 C COMPUTE THE MATRIX OF SINES AND COSINES C --------------------------------------- COSY1 = COS(Y1(I)*PI180) COSY2 = COS(Y2(I)*PI180) SINY1 = SIN(Y1(I)*PI180) SINY2 = SIN(Y2(I)*PI180) COSDX = COS((X1(I)-X2(I))*PI180) COSDY = COS((Y1(I)-Y2(I))*PI180) SINDX = SIN((X1(I)-X2(I))*PI180) SINDY = SIN((Y1(I)-Y2(I))*PI180) C COMPUTE THE NORMAL ANGLE C ------------------------ S = 1.0-COSDY+COSY1*COSY2*(1.0-COSDX) S = SQRT(2.*S) C CALCULATE THE FIRST PARTIALS WITH RESPECT TO X1 AND Y1 C ------------------------------------------------------ SI = MAX(S,SMIN) SI = 1./SI DY1 = SI*(SINDY-SINY1*COSY2*(1.0-COSDX)) DX1 = SI*(COSY1*COSY2*SINDX) DY2 = -SI*(SINDY+COSY1*SINY2*(1.0-COSDX)) DX2 = -DX1 C CALCULATE THE MIXED PARTIAL DERIVATIVES C --------------------------------------- DY1Y2= SI*(SINY1*SINY2*(1.0-COSDX)-COSDY-DY1*DY2) DY1X2= SI*(SINY1*COSY2*SINDX-DY1*DX2) DX1X2= -SI*(COSY1*COSY2*COSDX+DX1*DX2) DX1Y2= -SI*(COSY1*SINY2*SINDX+DX1*DY2) C COMPUTE THE VARIOUS CORRELATIONS C -------------------------------- RHO = EXP(-CH(I)*S*S) DRHO = -2.*CH(I)*S*RHO D2RHO = -2.*CH(I)*(RHO+S*DRHO) WW0 = 2.*CH(I) SRWW0 = SQRT(WW0) D(I) = S IF(S.LT.SMIN) THEN IF(ZZ) C1(I) = 1. IF(UU) C1(I) = 1. IF(VV) C1(I) = 1. D(I) = 0. ELSE IF(ZZ) THEN C1(I) = RHO ELSE IF(ZU) THEN C1(I) = FL(I)*TZ(I) * DRHO*DY2*SIGN(1.,-SINY2)/SRWW0 ELSE IF(ZV) THEN C1(I) = FL(I)*TZ(I) * DRHO*DX2*SIGN(1., SINY2)/(SRWW0*COSY2) ELSE IF(UZ) THEN C1(I) = FL(I)*TZ(I) * DRHO*DY1*SIGN(1.,-SINY1)/SRWW0 ELSE IF(UU) THEN C1(I) = DY1*DY2*D2RHO+DY1Y2*DRHO C1(I) = C1(I)*SIGN(1., SINY1*SINY2)/WW0 C1(I) = C1(I)*FL(I) + RHO*(1.-FL(I)) ELSE IF(UV) THEN C1(I) = DY1*DX2*D2RHO+DY1X2*DRHO C1(I) = FL(I) * C1(I)*SIGN(1.,-SINY1*SINY2)/(WW0*COSY2) ELSE IF(VZ) THEN C1(I) = FL(I)*TZ(I) * DRHO*DX1*SIGN(1., SINY1)/(SRWW0*COSY1) ELSE IF(VU) THEN C1(I) = DX1*DY2*D2RHO+DX1Y2*DRHO C1(I) = FL(I) * C1(I)*SIGN(1.,-SINY1*SINY2)/(WW0*COSY1) ELSE IF(VV) THEN C1(I) = DX1*DX2*D2RHO+DX1X2*DRHO C1(I) = C1(I)*SIGN(1., SINY1*SINY2)/(WW0*COSY1*COSY2) C1(I) = C1(I)*FL(I) + RHO*(1.-FL(I)) ENDIF 20 CONTINUE RETURN END SUBROUTINE COOR2(NP,CH,FL,TZ,X1,Y1,K1,X2,Y2,K2,D,C1,C2) DIMENSION CH(NP),FL(NP),TZ(NP) DIMENSION X1(NP),Y1(NP),K1(NP) DIMENSION X2(NP),Y2(NP),K2(NP) DIMENSION D (NP),C1(NP),C2(NP) LOGICAL ZZ,ZU,ZV,UZ,UU,UV,VZ,VU,VV DATA PI180 /.0174532 /, SMIN/.002/ C---------------------------------------------------------------------- C---------------------------------------------------------------------- IF(NP.EQ.0) RETURN C LOOP OVER SET OF INPUT POINTS C ----------------------------- DO 20 I=1,NP C1(I) = 0. C2(I) = 0. C DETERMINE CORRELATION TYPE C -------------------------- ZZ = K1(I).EQ.1 .AND. K2(I).EQ.1 ZU = K1(I).EQ.1 .AND. K2(I).EQ.2 ZV = K1(I).EQ.1 .AND. K2(I).EQ.3 UZ = K1(I).EQ.2 .AND. K2(I).EQ.1 UU = K1(I).EQ.2 .AND. K2(I).EQ.2 UV = K1(I).EQ.2 .AND. K2(I).EQ.3 VZ = K1(I).EQ.3 .AND. K2(I).EQ.1 VU = K1(I).EQ.3 .AND. K2(I).EQ.2 VV = K1(I).EQ.3 .AND. K2(I).EQ.3 C COMPUTE THE MATRIX OF SINES AND COSINES C --------------------------------------- COSY1 = COS(Y1(I)*PI180) COSY2 = COS(Y2(I)*PI180) SINY1 = SIN(Y1(I)*PI180) SINY2 = SIN(Y2(I)*PI180) COSDX = COS((X1(I)-X2(I))*PI180) COSDY = COS((Y1(I)-Y2(I))*PI180) SINDX = SIN((X1(I)-X2(I))*PI180) SINDY = SIN((Y1(I)-Y2(I))*PI180) C COMPUTE THE NORMAL ANGLE C ------------------------ S = 1.0-COSDY+COSY1*COSY2*(1.0-COSDX) S = SQRT(2.*S) C CALCULATE THE FIRST PARTIALS WITH RESPECT TO X1 AND Y1 C ------------------------------------------------------ SI = MAX(S,SMIN) SI = 1./SI DY1 = SI*(SINDY-SINY1*COSY2*(1.0-COSDX)) DX1 = SI*(COSY1*COSY2*SINDX) DY2 = -SI*(SINDY+COSY1*SINY2*(1.0-COSDX)) DX2 = -DX1 C CALCULATE THE MIXED PARTIAL DERIVATIVES C --------------------------------------- DY1Y2= SI*(SINY1*SINY2*(1.0-COSDX)-COSDY-DY1*DY2) DY1X2= SI*(SINY1*COSY2*SINDX-DY1*DX2) DX1X2= -SI*(COSY1*COSY2*COSDX+DX1*DX2) DX1Y2= -SI*(COSY1*SINY2*SINDX+DX1*DY2) C COMPUTE THE VARIOUS CORRELATIONS C -------------------------------- RHO = EXP(-CH(I)*S*S) DRHO = -2.*CH(I)*S*RHO D2RHO = -2.*CH(I)*(RHO+S*DRHO) WW0 = 2.*CH(I) SRWW0 = SQRT(WW0) D(I) = S IF(S.LT.SMIN) THEN IF(ZZ) C1(I) = 1. IF(ZZ) C2(I) = 1. IF(UU) C1(I) = 1. IF(VV) C1(I) = 1. IF(VU) C2(I) = 1. D(I) = 0. ELSE IF(ZZ) THEN C1(I) = RHO C2(I) = C1(I) ELSE IF(ZU) THEN C1(I) = FL(I)*TZ(I) * DRHO*DY2*SIGN(1.,-SINY2)/SRWW0 C2(I) = FL(I)*TZ(I) * DRHO*DX2*SIGN(1., SINY2)/(SRWW0*COSY2) ELSE IF(ZV) THEN C1(I) = FL(I)*TZ(I) * DRHO*DX2*SIGN(1., SINY2)/(SRWW0*COSY2) C2(I) = C1(I) ELSE IF(UZ) THEN C1(I) = FL(I)*TZ(I) * DRHO*DY1*SIGN(1.,-SINY1)/SRWW0 C2(I) = C1(I) ELSE IF(UU) THEN C1(I) = DY1*DY2*D2RHO+DY1Y2*DRHO C1(I) = C1(I)*SIGN(1., SINY1*SINY2)/WW0 C1(I) = C1(I)*FL(I) + RHO*(1.-FL(I)) C2(I) = DY1*DX2*D2RHO+DY1X2*DRHO C2(I) = FL(I) * C2(I)*SIGN(1.,-SINY1*SINY2)/(WW0*COSY2) ELSE IF(UV) THEN C1(I) = DY1*DX2*D2RHO+DY1X2*DRHO C1(I) = FL(I) * C1(I)*SIGN(1.,-SINY1*SINY2)/(WW0*COSY2) C2(I) = C1(I) ELSE IF(VZ) THEN C1(I) = FL(I)*TZ(I) * DRHO*DX1*SIGN(1., SINY1)/(SRWW0*COSY1) C2(I) = C1(I) ELSE IF(VU) THEN C1(I) = DX1*DY2*D2RHO+DX1Y2*DRHO C1(I) = FL(I) * C1(I)*SIGN(1.,-SINY1*SINY2)/(WW0*COSY1) C2(I) = DX1*DX2*D2RHO+DX1X2*DRHO C2(I) = C2(I)*SIGN(1., SINY1*SINY2)/(WW0*COSY1*COSY2) C2(I) = C2(I)*FL(I) + RHO*(1.-FL(I)) ELSE IF(VV) THEN C1(I) = DX1*DX2*D2RHO+DX1X2*DRHO C1(I) = C1(I)*SIGN(1., SINY1*SINY2)/(WW0*COSY1*COSY2) C1(I) = C1(I)*FL(I) + RHO*(1.-FL(I)) C2(I) = C1(I) ENDIF 20 CONTINUE RETURN END C$$$ SUBPROGRAM DOCUMENTATION BLOCK C C SUBPROGRAM: DRCTSL C PRGMMR: WOOLLEN ORG: NMC22 DATE: 90-11-06 C C ABSTRACT: DRIVER FOR CHOLESKY TYPE LINEAR EQUATION SOLVER. C C PROGRAM HISTORY LOG: C 90-11-06 J. WOOLLEN C C USAGE: C INPUT ARGUMENTS: C FAALL - ARRAY OF SYMMETRIC MATRIXES C RAALL - ARRAY OF MATRIX RIGHT HAND SIDES C NDIM - ARRAY OF MATRIX RANKS C MAXDIM - MAXIMUM RANK MATRIX IN MATRIX ARRAY C NXXYY - NUMBER OF MATRIXES IN MATRIX ARRAY C NFT - NUMBER OF RIGHT HAND SIDE VECTORS PER MATRIX C OUTPUT ARGUMENTS: C RAALL - ARRAY OF MATRIX SOLUTIONS C DOTPRD - ARRAY OF DOT PRODUCTS OF RHS VECTORS WITH MATRIX C SOLUTIONS C C SUBPROGRAMS CALLED: C UNIQUE: - VSOLVE C C REMARKS: ILL CONDITIONED OR NON POSITIVE DEFINITE MATRIXES ARE C IDENTIFIED BY DOT PRODUCTS GT 1 OR LT 0 OR BY A MESSAGE C FROM VSOLVE. FIVE ATTEMPTS ARE MADE TO SOLVE BAD ONES, C BY RIDGE REGRESSION, AFTER WHICH A NULL SOLUTION IS RETURNED. C C ATTRIBUTES: C LANGUAGE: CRAY FORTRAN C MACHINE: CRAY C C$$$ C----------------------------------------------------------------------- SUBROUTINE DRCTSL(FA,RA,DP,NDIM,NXXYY) PARAMETER( NMX = 64 , . NN = 10 *( 10 +1)/2 , . N = 10 ) DIMENSION FA(NMX,NN), RA(NMX,N,2), DP(NMX,2), NDIM(NXXYY) DIMENSION A (NMX,NN), B (NMX,N,2) DIMENSION BAD(NMX), SMOOTH(6) LOGICAL BAD DATA SMOOTH /1.00,1.01,1.02,1.05,1.10,2.00/ C---------------------------------------------------------------------- C---------------------------------------------------------------------- C LOOP FOR POSSIBILITY OF BAD MATRIXS C ----------------------------------- DO 50 ITRY=1,6 NBAD = 0 C FIND THE LARGEST MATRIX IN THE STACK C ------------------------------------ MAXDIM = 0 DO 5 I=1,NXXYY MAXDIM = MAX(MAXDIM,NDIM(I)) 5 CONTINUE IF(MAXDIM.EQ.0) RETURN C INITIALIZE THE WORKING ARRAYS C ----------------------------- DO 10 J=1,MAXDIM*(MAXDIM+1)/2 DO 10 I=1,NXXYY A(I,J) = FA(I,J) 10 CONTINUE DO 11 J=1,MAXDIM DO 11 I=1,NXXYY B(I,J,1) = RA(I,J,1) B(I,J,2) = RA(I,J,2) 11 CONTINUE DO 12 I=1,NXXYY DP(I,1) = 0. DP(I,2) = 0. 12 CONTINUE DO 13 J=1,MAXDIM JJ = J*(J+1)/2 DO 13 I=1,NXXYY A(I,JJ) = FA(I,JJ)*SMOOTH(ITRY) 13 CONTINUE C CALL THE MATRIX SOLVER C ---------------------- CALL VSOLVE (A,B,NDIM,BAD,NXXYY,MAXDIM) C MAKE THE DOT PRODUCTS OF SOLUTIONS WITH RIGHT HAND SIDES C -------------------------------------------------------- DO 20 J=1,MAXDIM DO 20 I=1,NXXYY DP(I,1) = DP(I,1) + RA(I,J,1)*B(I,J,1) DP(I,2) = DP(I,2) + RA(I,J,2)*B(I,J,2) 20 CONTINUE C NORMALIZE THE WEIGHTS IF THEY ARE TOO BIG C ----------------------------------------- DO 25 I=1,NXXYY IF(DP(I,1).GT.1. .OR. DP(I,2).GT.1.) THEN DO 24 J=1,MAXDIM B(I,J,1) = B(I,J,1)/DP(I,1) B(I,J,2) = B(I,J,2)/DP(I,2) 24 CONTINUE DP(I,1) = 1. DP(I,2) = 1. ENDIF 25 CONTINUE C CHECK FOR BAD ONES - DO IT AGAIN IF NECESSARY C --------------------------------------------- DO 30 I=1,NXXYY BAD(I) = BAD(I) .OR. DP(I,1).LT.0. .OR. DP(I,2).LT.0. BAD(I) = BAD(I) .OR. DP(I,1).GT.1. .OR. DP(I,2).GT.1. 30 CONTINUE DO 40 I=1,NXXYY IF(BAD(I)) THEN DP(I,1) = 0. DP(I,2) = 0. DO 35 J=1,MAXDIM B(I,J,1) = 0. B(I,J,2) = 0. 35 CONTINUE NBAD = NBAD + 1 ENDIF 40 CONTINUE IF(NBAD.NE.0) PRINT*, 'NBAD=',NBAD,' ON TRY ',ITRY IF(NBAD.EQ.0) GOTO 55 50 CONTINUE C COPY SOLUTIONS INTO OUTPUT ARRAY - ZERO BAD ONES OUT C ---------------------------------------------------- 55 DO 60 J=1,MAXDIM DO 60 I=1,NXXYY RA(I,J,1) = B(I,J,1) RA(I,J,2) = B(I,J,2) 60 CONTINUE RETURN END C$$$ SUBPROGRAM DOCUMENTATION BLOCK C . . . . C SUBPROGRAM: GETFCER GET FORECAST ERRORS C PRGMMR: WOOLLEN ORG: NMC22 DATE: 92-07-29 C C ABSTRACT: FUNCTION WHICH RETURNS A FORECAST ERROR IN THE FORCAST C HEIGHT OR WIND FIELD FOR A GIVEN LAT, LON, PRESSURE. ENTRY POINTS C GETZER AND GETWER RETURN ERROR VALUES; GETFCER RETURNS 0. C C PROGRAM HISTORY LOG: C 92-07-29 J. WOOLLEN C C USAGE: CALL GETZER(SLAT,SLON,PRES) - (RETURNS HEIGHT ERROR) C CALL GETWER(SLAT,SLON,PRES) - (RETURNS WIND ERROR) C C INPUT ARGUMENTS: C SLAT - LATITUDE C SLON - LONGITUDE C PRES - PRESSURE C C FUNCTION RETURN VALUE: C GETZER - FORECAST HEIGHT ERROR C GETWER - FORECAST WIND ERROR C C SUBPROGRAMS CALLED: C C EXIT STATES: C C REMARKS: C C C ATTRIBUTES: C LANGUAGE: CRAY FORTRAN C MACHINE: CRAY C C----------------------------------------------------------------------- FUNCTION GETFCER(SLAT,SLON,PRES) COMMON /OBERRS/ QIKE(21,300),PMAND(21),P50(21),CHL(21) COMMON /CORP / CHLP(1500),CVC,FFLAT(91),IFA( 10 , 10 ) COMMON /FCRATS/ RATION(21),RATIOS(21),TROPUV(21) LOGICAL WIND DATA PI180 / .0174532 / DATA OMEGA / .72921E-4 / DATA RADE / 6.371E6 / DATA GRAV / 9.8 / C--------------------------------------------------------------------- C--------------------------------------------------------------------- GETFCER = 0. RETURN ENTRY GETZER(SLAT,SLON,PRES) WIND = .FALSE. GOTO 10 ENTRY GETWER(SLAT,SLON,PRES) WIND = .TRUE. C COMPUTE A HEIGHT ERROR PROPORTIONAL TO RAWINSONDE ERROR C ------------------------------------------------------- 10 RAWZER = PILNLNP(PRES,PMAND,QIKE(1,120),21) IF(SLAT.LT.0.) THEN ZER = RAWZER/SQRT(PILNLNP(PRES,PMAND,RATIOS,21)) * 1.1 ELSE ZER = RAWZER/SQRT(PILNLNP(PRES,PMAND,RATION,21)) ENDIF IF(.NOT.WIND) THEN GETZER = ZER RETURN ENDIF C COMPUTE THE GEOSTROPHIC WIND ERROR FOR GETWER C --------------------------------------------- TRPWER = PILNLNP(PRES,PMAND,TROPUV,21) IPRS = MIN(1,MAX(NINT(PRES),1000)) SRWW0 = SQRT(2.0*CHLP(IPRS)) IF(SLAT.EQ.0.) ALPHA = 0. IF(SLAT.NE.0.) ALPHA = GRAV/(2.*OMEGA*SIN(ABS(SLAT)*PI180)*RADE) IY = ABS(SLAT)+1.5 G = FFLAT(IY) GETWER = G*ZER*ALPHA*SRWW0 + (1-G)*TRPWER RETURN END C----------------------------------------------------------------------- C C----------------------------------------------------------------------- SUBROUTINE MAKEAMAP COMMON /OBS/ NREP,NLEV,SID( 200000 ),INOB( 200000 ),IMAP(360,181) COMMON /CKLIST/ NDD,INDD( 300000 ) C----------------------------------------------------------------------- C----------------------------------------------------------------------- C MAKE SURE NONE OF THE LONGITUDES EQUAL 360. C ------------------------------------------- DO 1 N=1,NREP IF(FLN(N).EQ.360.) PACK = WLN(N,0.) 1 CONTINUE C LONGITUDE SORT C -------------- NN = 0 DO 10 NX=1,360 DO 10 N=1,NREP IX = FLN(N)+1. IF(IX.EQ.NX) THEN NN = NN+1 INDD(NN) = N ENDIF 10 CONTINUE C LATITUDE SORT C ------------- NOB = 0 DO 20 NY=1,181 DO 20 N=1,NN IY = FLT(INDD(N))+91. IF(IY.EQ.NY) THEN NOB = NOB+1 INOB(NOB) = INDD(N) ENDIF 20 CONTINUE C CHECK TO MAKE SURE ALL REPS ARE ACCOUNTED FOR C --------------------------------------------- IF(NREP.NE.NN .OR. NREP.NE.NOB) THEN PRINT* PRINT*,'*****************************************************' PRINT*,'** WARNING*WARNING*WARNING*WARNING*WARNING*WARNING **' PRINT*,'*****************************************************' PRINT*,'** IN MAKEAMAP ALL REPORTS HAVE NOT BEEN PROCESSED **' PRINT*,'*****************************************************' PRINT*,'** WARNING*WARNING*WARNING*WARNING*WARNING*WARNING **' PRINT*,'*****************************************************' PRINT* ENDIF C INITIALIZE IMAP ARRAY C --------------------- C-CRA IMAP = 0 DO IJ=1,360*181 IMAP(IJ,1) = 0 ENDDO C DO 30 N=1,NOB NX = FLN(INOB(N)) + 1. NY = FLT(INOB(N)) + 91. IF(IMAP(NX,NY).EQ.0) IMAP(NX,NY) = N 30 CONTINUE C FILL GAPS IN IMAP FOR CONTINUITY (BACKWARDS) C -------------------------------------------- LASTN = NOB DO 40 NY=181,1,-1 DO 40 NX=360,1,-1 IF(IMAP(NX,NY).EQ.0) THEN IMAP(NX,NY) = LASTN ELSE LASTN = IMAP(NX,NY) ENDIF 40 CONTINUE RETURN END C$$$ SUBPROGRAM DOCUMENTATION BLOCK C C SUBPROGRAM: OBOGRAM C PRGMMR: WOOLLEN ORG: NMC22 DATE: 90-11-06 C C ABSTRACT: EXAMINES THE RESULTS OF FINAL QC DECISIONS AND MAKES C TABLES OF STATISTICS FOR ACCEPTED AND REJECTED REPORTS BY C REPORT TYPE AND LEVEL. C C PROGRAM HISTORY LOG: C 90-11-06 J. WOOLLEN C C USAGE: CALL OBOGRAM(LASC,LBIN) C C INPUT ARGUMENTS: C LASC - LOGICAL UNIT FOR ASCII OUTPUT C LBIN - LOGICAL UNIT FOR BINARY OUTPUT C C OUTPUT ARGUMENTS: NONE C C OUTPUT FILES: C FORT.LASC - OBOGRAM.OUT - SUMMARY OF QC RESULTS (ASCII) C FORT.LBIN - OBOGRAM.BIN - SUMMARY OF QC RESULTS (BINARY) C C SUBPROGRAMS CALLED: C LIBRARY: C COMMON - ABORT C C REMARKS: C C C ATTRIBUTES: C LANGUAGE: CRAY FORTRAN C MACHINE: CRAY C C$$$ C----------------------------------------------------------------------- SUBROUTINE OBOGRAM(LDAT,LASC,LBIN) COMMON /OBS/ NREP,NLEV,SID( 200000 ),INOB( 200000 ),IMAP(360,181) COMMON /OBERRS/ QIKE(21,300),PMAND(21),P50(21),CHL(21) CHARACTER*3 AR(2),TAG,FOC DIMENSION STAT(300,15,6,2,2,3),IOBS(300),RZONE(6) LOGICAL QOBS(300,6) DATA AR/'ACC','REJ'/ DATA XMSG/32767./ C---------------------------------------------------------------------- C----------------------------------------------------------------------- REWIND (LDAT) READ(LDAT,1) NMCDATE 1 FORMAT(8X,I8) WRITE(LASC,*) WRITE(LASC,*)' STATS FOR OBS CHECKED ',NMCDATE WRITE(LASC,*) C START UP C -------- C-CRA QOBS = .FALSE. C LOGICAL QOBS(300,6) DO IJ=1,6*300 QOBS(IJ,1) = .FALSE. ENDDO C DIMENSION STAT(300,15,6,2,2,3),IOBS(300),RZONE(6) C-CRA STAT = 0 DO IJ=1,300*15*6*2*2*3 STAT(IJ,1,1,1,1,1) = 0 ENDDO C GO THROUGH THE BIG BUFFER TO COUNT OBS FOR OBOGRAM C -------------------------------------------------- DO 10 IREP=1,NREP X = FLN(IREP) Y = FLT(IREP) KX = FKX(IREP) KK = KX/100 LB = MIN(IFIX((Y+90.)/30.+1.),6) TAG = FOC(IREP) IF(TAG.EQ.'SFC') SFC2 = WNL(IREP,2) ILEV = FLV(IREP) KLEV = ILEV+FNL(IREP)-1 C GO THROUGH LEVELS FOR ACCEPT REJECT STATISTICS C ---------------------------------------------- DO 10 JLEV=ILEV,KLEV IF(FOI(JLEV).NE.0) THEN IF(FOI(JLEV).EQ.1) THEN IAR = 1 ELSE IF(FOI(JLEV).EQ.-1) THEN IAR = 2 ELSE GOTO 10 ENDIF BUFF11 = FOB(JLEV,1) BUFF12 = FOB(JLEV,2) IF(TAG.EQ.'SFC' .AND. JLEV.EQ.KLEV) THEN KX = KX - 100 ELSE IF(KK.EQ.1) BUFF12 = BUFF11 ENDIF QOBS(KX,LB) = .TRUE. PP = FPR(JLEV) DO 5 L=1,15 IF(L.EQ.1) THEN IF(PP.GT.(PMAND(L)+PMAND(L+1))/2.) GOTO 6 ELSE IF(L.EQ.15) THEN IF(PP.LE.(PMAND(L)+PMAND(L-1))/2.) GOTO 6 ELSE IF(PP.GT.(PMAND(L)+PMAND(L+1))/2. .AND. . PP.LE.(PMAND(L)+PMAND(L-1))/2.) GOTO 6 ENDIF 5 CONTINUE CALL SABORT('OBOGRAM - UNABLE TO LOCATE DATA LEVEL') 6 CONTINUE IF(KK.EQ.1) THEN STAT(KX,L,LB,IAR,1,1) = STAT(KX,L,LB,IAR,1,1) + 1. STAT(KX,L,LB,IAR,1,2) = STAT(KX,L,LB,IAR,1,2) + BUFF12 STAT(KX,L,LB,IAR,1,3) = STAT(KX,L,LB,IAR,1,3) + BUFF12**2 ELSE STAT(KX,L,LB,IAR,1,1) = STAT(KX,L,LB,IAR,1,1) + 1. STAT(KX,L,LB,IAR,1,2) = STAT(KX,L,LB,IAR,1,2) + BUFF11 STAT(KX,L,LB,IAR,1,3) = STAT(KX,L,LB,IAR,1,3) + BUFF11**2 STAT(KX,L,LB,IAR,2,1) = STAT(KX,L,LB,IAR,2,1) + 1. STAT(KX,L,LB,IAR,2,2) = STAT(KX,L,LB,IAR,2,2) + BUFF12 STAT(KX,L,LB,IAR,2,3) = STAT(KX,L,LB,IAR,2,3) + BUFF12**2 ENDIF ENDIF 10 CONTINUE C COMPUTE THE AVERAGES FOR MEAN AND RMS'S (VECTOR SUMS FOR WINDS) C --------------------------------------------------------------- DO 20 KX=1,299 KK = MAX(KX/100,1) DO 20 L=1,15 DO 20 LB=1,6 DO 20 IAR=1,2 DO 20 KKX=1,KK N = STAT(KX,L,LB,IAR,KKX,1) IF(N.NE.0) THEN BIAS = STAT(KX,L,LB,IAR,KKX,2)/N STAT(KX,L,LB,IAR,KKX,2) = BIAS RMS = SQRT(STAT(KX,L,LB,IAR,KKX,3)/N) STAT(KX,L,LB,IAR,KKX,3) = RMS IF(KKX.EQ.2) THEN STAT(KX,L,LB,IAR,1,2) = SQRT(STAT(KX,L,LB,IAR,1,2)**2+BIAS**2) STAT(KX,L,LB,IAR,1,3) = SQRT(STAT(KX,L,LB,IAR,1,3)**2+ RMS**2) ENDIF ENDIF 20 CONTINUE C PRINT THE TOTAL, MEAN, AND RMS FOR ACCEPTS AND REJECTS C ------------------------------------------------------ DO 35 KX=1,300 DO 35 LB=1,6 IF(QOBS(KX,LB)) THEN KK = KX/100 WRITE(81,*)'------ZONE ',LB,' KX ',KX,' VAR ',KK,'------' DO 30 IAR=1,2 WRITE(LASC,2000) AR(IAR),(STAT(KX,L,LB,IAR,1,1),L=1,15) WRITE(LASC,2001) AR(IAR),(STAT(KX,L,LB,IAR,1,2),L=1,15) WRITE(LASC,2002) AR(IAR),(STAT(KX,L,LB,IAR,1,3),L=1,15) WRITE(LASC,2003) 30 CONTINUE ENDIF 35 CONTINUE C PRINT THE TOTAL REJECTED OB BY TYPE ZONE LEVEL C ---------------------------------------------- WRITE(LASC,*) WRITE(LASC,*) ' NUMBER OF OBS REJECTED BY OB TYPE, AND ZONE' WRITE(LASC,*) DO 45 KX=1,300 IF(QOBS(KX,1).OR.QOBS(KX,2).OR.QOBS(KX,3).OR.QOBS(KX,4).OR. . QOBS(KX,5).OR.QOBS(KX,6)) THEN C DIMENSION STAT(300,15,6,2,2,3),IOBS(300),RZONE(6) C-CRA RZONE = 0 DO LB=1,6 RZONE(LB) = 0 ENDDO C DO 40 LB=1,6 DO 40 L=1,15 RZONE(LB) = STAT(KX,L,LB,2,1,1) + RZONE(LB) 40 CONTINUE WRITE(LASC,2004)KX,(RZONE(LB),LB=1,6) ENDIF 45 CONTINUE C WRITE A BINARY VERSION OF THE OBOGRAM C ------------------------------------- I = 0 DO 55 KX=1,300 DO 50 LB=1,6 IF(QOBS(KX,LB)) THEN I = I+1 IOBS(I) = KX GOTO 55 ENDIF 50 CONTINUE 55 CONTINUE WRITE(LBIN) NMCDATE, . I,(IOBS(II),((((STAT(IOBS(II),LEV,LB,IAR,1,IST), . LEV=1,15), . IST=1,3), . IAR=1,2), . LB=1,6), . II=1,I) C FORMAT STATEMENTS FOR THE OBOGRAM C --------------------------------- 2000 FORMAT('TOTAL ',A3,15F6.1) 2001 FORMAT('MEAN ',A3,15F6.1) 2002 FORMAT('RMS ',A3,15F6.1) 2003 FORMAT(' ') 2004 FORMAT(I5,6F8.1) RETURN END C$$$ SUBPROGRAM DOCUMENTATION BLOCK C C SUBPROGRAM: OUTPUT C PRGMMR: WOOLLEN ORG: NMC22 DATE: 90-11-06 C C ABSTRACT: EXAMINES OUTPUTS FROM OI QUALITY CONTROL INTERPOLATION C AND MAKES DECISION WHETHER AN OB PASSES OR FAILS THE CHECKS. C ALSO PROVIDES FOR A DETAILED PRINT TRACING ALL DECISIONS MADE C IN REGARD TO SPECIFIED OBSERVATIONS IN THE QC PROCESS. C C PROGRAM HISTORY LOG: C 90-11-06 J. WOOLLEN C C USAGE: CALL OUTPUT C C INPUT ARGUMENTS: NONE C C OUTPUT ARGUMENTS: IN COMMON /OBMARK/ C C SUBPROGRAMS CALLED: C LIBRARY: C COMMON - ABORT C C REMARKS: C C C ATTRIBUTES: C LANGUAGE: CRAY FORTRAN C MACHINE: CRAY C C$$$ SUBROUTINE OUTPUT C-CRA TASK COMMON/HVECT/NXXYY,MCHK( 64 ,2),MOBS( 64 , 10 , 4 ,2), COMMON/HVECT/NXXYY,NYYZZ,MCHK( 64 ,2),MOBS( 64 , 10 , 4 ,2) 1, . MDIM( 64 , 4 ),AA( 64 , 4 ,2), . AE( 64 , 4 ,2),WT( 64 , 10 , 4 ,2), . FE( 64 ,3) COMMON /OBS/ NREP,NLEV,SID( 200000 ),INOB( 200000 ),IMAP(360,181) C-CRA COMMON /TOLLER / ITER,TB,TZ,TSAT,TWND,TJET C-CRA COMMON /OBPRT / NPT,SIDP(10),PRP(10),KXP(10) COMMON /TOLLER / ITER,ITERADD,TB,TZ,TSAT,TWND,TJET COMMON /OBPRT / NPT,NPTADD,SIDP(10),PRP(10),KXP(10) COMMON /TOSSUNIT/ ITUNIT DIMENSION IOMF(3),D( 4 ),DRES( 4 ),NDIM( 4 ),STR(2) DIMENSION DU( 4 ),DV( 4 ),TWT( 4 ),TWU( 4 ),TWV( 4 ) CHARACTER*8 CKEEP,CHECK(5) CHARACTER*3 TAG,FOC CHARACTER*1 VAR2(2),VAR3(3) LOGICAL KEEP DATA CHECK/'Z MATRIX','U MATRIX','V MATRIX','VERTICAL','MULTIVAR'/ DATA VAR2 / 'Z' , 'W' / DATA VAR3 / 'Z' , 'U' , 'V' / DATA IOMF / 2 , 1 , 2 / CHARACTER*1 PLUMIN DATA PI180/.017453/ C---------------------------------------------------------------------- C---------------------------------------------------------------------- C-CRA D = 0. C-CRA DU = 0. C-CRA DV = 0. C-CRA TWT = 0. C-CRA TWU = 0. C-CRA TWV = 0. C DIMENSION IOMF(3),D( 4 ),DRES( 4 ),NDIM( 4 ),STR(2) C DIMENSION DU( 4 ),DV( 4 ),TWT( 4 ),TWU( 4 ),TWV( 4 ) DO IJ=1, 4 D(IJ) = 0. DU(IJ) = 0. DV(IJ) = 0. TWT(IJ) = 0. TWU(IJ) = 0. TWV(IJ) = 0. ENDDO C C LOOP OVER EACH OB CHECKED IN THIS GROUP C --------------------------------------- DO 50 N=1,NXXYY I = MCHK(N,1) J = MCHK(N,2) II = FHD(I) PRES = FPR(I) OBER = FOE(I) CKSTA = SID(II) CLON = FLN(II) CLAT = FLT(II) KX = FKX(II) TAG = FOC(II) C LOOK FOR A KEEP FLAG C -------------------- IQM = FQM(I) IF(IQM.EQ.0) THEN KEEP = .TRUE. CKEEP = ' KEEP' ELSE KEEP = .FALSE. CKEEP = ' ' ENDIF C LOOP OVER THE VARIABLE(S) IN EACH OB C ------------------------------------ DO 45 K=J,2*J-1 KK = K-J+1 OMF = FOB(I,IOMF(K)) AMFSUM = 0. SUMWTS = 0. QC = 0. C LOOP OVER THE OUTCOME OF EACH TYPE OF CHECK (MICRO DECISION MAKER) C ------------------------------------------------------------------ DO 20 L=1, 4 AER = AE(N,L,KK) FER = FE(N,K) AMF = AA(N,L,KK) DRES(L) = OMF-AMF D(L) = DRES(L)/SQRT(OBER**2+AER**2) AWT = 1.-(AER/FER)**2 TWT(L) = MAX(.01,AWT) QC = QC + TWT(L)*D(L)**2 SUMWTS = SUMWTS + TWT(L) IF(L.LT.4) AMFSUM = AMFSUM + TWT(L)*AMF 20 CONTINUE AMFSUM = AMFSUM/(SUMWTS-TWT(4)) QC = SQRT(QC/SUMWTS) C-CRA TWT = TWT/SUMWTS DO IJ=1, 4 TWT(IJ) = TWT(IJ)/SUMWTS ENDDO C FOR WINDS CHECK THE VECTOR SUM OF AGAINST THE REFERENCE C ------------------------------------------------------- IF(K.EQ.1) THEN OMA = OMF - AMFSUM ELSE IF(K.EQ.2) THEN C-CRA DU = D C-CRA TWU = TWT DO IJ=1, 4 DU(IJ) = D(IJ) TWU(IJ) = TWT(IJ) ENDDO QCU = QC OMAU = OMF - AMFSUM OMFU = OMF GOTO 45 ELSE IF(K.EQ.3) THEN C-CRA DV = D C-CRA TWV = TWT DO IJ=1, 4 DV(IJ) = D(IJ) TWV(IJ) = TWT(IJ) ENDDO QCV = QC OMAV = OMF - AMFSUM OMFV = OMF C-CRA D = SQRT(DU**2+DV**2) C-CRA TWT = SQRT(TWU**2+TWV**2) DO IJ=1, 4 D(IJ) = SQRT(DU(IJ)**2+DV(IJ)**2) TWT(IJ) = SQRT(TWU(IJ)**2+TWV(IJ)**2) ENDDO QC = SQRT(QCU**2+QCV**2) OMA = SQRT(OMAU**2+OMAV**2) OMF = SQRT(OMFU**2+OMFV**2) ELSE PRINT*,'OUTPUT - K<>1,2,3' CALL SABORT('OUTPUT - K<>1,2,OR,3') ENDIF ENDIF C ADJUST THE TOLERANCE BY LOCATION AND VARIABLE C --------------------------------------------- IF(K.EQ.1) THEN SINLAT = (SIN(ABS(CLAT)*PI180)) TT = TB*TZ*SINLAT TMIN = 2.0 ELSE GNP = MIN(ABS(CLAT/30.),1.) IF(PRES.LE.400. .AND. PRES.GE.150.) THEN TT = TB*TJET*GNP ELSE TT = TB*TWND*GNP ENDIF TMIN = 3. ENDIF TT = MAX(TMIN,TT) IF(TAG.EQ.'SMI') TT = TT/SQRT(2.) C COMPUTE THE MULTIVARIATE REDUCTION RATIO AND ADJUST THE SUM C ----------------------------------------------------------- IF(OMF.NE.0.) THEN QCPC = MIN(MAX(ABS(OMA/OMF),.1),1.) TTT = TT/QCPC ELSE TTT = TT ENDIF C COMPUTE THE NORMALIZED DEVIATION FROM THE FIRST GUESS C ----------------------------------------------------- SDOMF = SQRT(OMF**2/(OBER**2+FER**2)) C MAKE A MESO DECISION FOR THIS ITERATION C --------------------------------------- IF(QC.GT.TTT) THEN KEEP = KEEP .AND. QC.LT.25. IF(.NOT.KEEP) THEN PACK = WOI(I,-1) PLUMIN = '-' ELSE PACK = WOI(I,1) PLUMIN = '+' ENDIF WRITE(ITUNIT,100)CKSTA,KX,CLON,CLAT,PRES,VAR2(J),OMF,SDOMF, . (D(L),TWT(L),L=1,4),QC,TTT,CKEEP,PLUMIN ELSE PACK = WOI(I,1) PLUMIN = '+' ENDIF C THIS WRITES DETAILED INFO ABOUT AN OB C ------------------------------------- DO 41 NN=1,NPT LU = 70 + NN pdif = abs(pres-prp(nn)) IF(CKSTA.EQ.SIDP(NN).AND.pdif.lt..1.AND.KX.EQ.KXP(NN)) THEN WRITE(LU,100)CKSTA,KX,CLON,CLAT,PRES,VAR2(J),OMF,SDOMF, . (D(L),TWT(L),L=1,4),QC,TTT,CKEEP,PLUMIN DO 44 KK=J,J*2-1 WRITE(LU,*) 'ITERATION=',ITER,' CHECK ON ',VAR3(KK) IF(KK.EQ.2) THEN SDOMF = SQRT(OMFU**2/(OBER**2+FER**2)) IF(TT.EQ.TTT) TTTT = TT IF(TT.NE.TTT) TTTT = TT/ABS(OMAU/OMFU) WRITE(LU,100)CKSTA,KX,CLON,CLAT,PRES,VAR3(KK),OMFU,SDOMF, . (DU(L),TWU(L),L=1,4),QCU,TTTT,CKEEP ELSE IF(KK.EQ.3) THEN SDOMF = SQRT(OMFV**2/(OBER**2+FER**2)) IF(TT.EQ.TTT) TTTT = TT IF(TT.NE.TTT) TTTT = TT/ABS(OMAV/OMFV) WRITE(LU,100)CKSTA,KX,CLON,CLAT,PRES,VAR3(KK),OMFV,SDOMF, . (DV(L),TWV(L),L=1,4),QCV,TTTT,CKEEP ENDIF DO 44 L=1, 4 WRITE(LU,*)CHECK(L),' FOR CHECK ON ',VAR3(KK),' DATUM ' SOINC = 0. DO 43 IOB=1,MDIM(N,L) ILEV = MOBS(N,IOB,L,1) KIV = MOBS(N,IOB,L,2) KI2 = MAX(KIV-1,1) IREP = FHD(ILEV) X = FLN(IREP) Y = FLT(IREP) OBID = SID(IREP) KIX = FKX(IREP) PPP = FPR(ILEV) OMF0 = FOB(ILEV,IOMF(KIV)) OIWT = WT(N,IOB,L,KI2) OINC = OIWT*OMF0 OMA0 = AA(N,L,KI2) AER0 = AE(N,L,KI2) WRITE(LU,102)SID(IREP),KIX,X,Y,PPP,VAR3(KIV),OMF0,OIWT,OINC, . OMA0,AER0 SOINC = SOINC+OINC 43 CONTINUE WRITE(LU,*) 'INTERPOLATED VALUE=',SOINC 44 CONTINUE ENDIF 41 CONTINUE 45 CONTINUE 50 CONTINUE 100 FORMAT(A6,I4,3F8.2,A4,2F8.2,4(' *',F6.1,' (',F3.2,')' ), . ' *',2F6.2,A6,A1) 101 FORMAT(A6,I4,3F8.2,A4,2F8.2,4(' *',F6.1,2X,I2),' *',2F6.2,A6,A1) 102 FORMAT(A6,I4,3F8.2,A4,6F8.2) RETURN END C$$$ SUBPROGRAM DOCUMENTATION BLOCK C . . . . C SUBPROGRAM: PILNLNP PRESSURE INTERPOLATION LINEAR IN LOG P C PRGMMR: WOOLLEN ORG: W/NMC22 DATE: 92-07-29 C C ABSTRACT: FUNCTION WHICH GIVEN AN PROFILE OF DESCENDING PRESSURES C AND A PROFILE OF QUANTITIES VALID AT THOSE PRESSURES, RETURNS C AN INTERPOLATED QUANTITY VALID AT A GIVEN ARBITRARY PRESSURE. C C PROGRAM HISTORY LOG: C 86-03-21 G. DIMEGO C 88-11-24 D. DEAVEN RECODED FOR CYBER 205 C C USAGE: X = PILNLNP(P,PARAY,QARAY,KMAX) C C INPUT ARGUMENTS: C P - PRESSURE TO INTERPOLATE TO C PARAY - GIVEN PRESSURE PROFILE C QARAY - QUANTITIES VALID FOR PRESSURE PROFILE C KMAX - LENGTH OF PROFILE C C FUNCTION RETURN VALUE C PILNLNP - INTERPOLATED QUANTITY C C ATTRIBUTES: C LANGUAGE: CDC FORTRAN200 C MACHINE: CYBER C C$$$ FUNCTION PILNLNP(P,PARAY,QARAY,KMAX) DIMENSION PARAY(KMAX),QARAY(KMAX) DO 10 LA=1,KMAX IF(P.GE.PARAY(LA)) GOTO 11 10 CONTINUE 11 IF(LA.GT.KMAX) THEN LA = KMAX LB = KMAX ELSE IF(LA.EQ.1) THEN LA = 1 LB = 1 ELSE LB = LA-1 ENDIF PA = PARAY(LA) PB = PARAY(LB) IF(PA.NE.PB) THEN WK = LOG(P/PB) / LOG(PA/PB) ELSE WK = 0. ENDIF PILNLNP = QARAY(LB) + (QARAY(LA)-QARAY(LB)) * WK RETURN END C$$$ SUBPROGRAM DOCUMENTATION BLOCK C C SUBPROGRAM: QCLOOP C PRGMMR: WOOLLEN ORG: NMC22 DATE: 90-11-06 C C ABSTRACT: ITERATES THE QC DECSION MAKING PROCESS UNTIL ALL DECISIONS C IN TWO CONSEQUETIVE ITERATIONS ARE THE SAME. AFTER FOUR ITERATIONS C ONLY REPORTS THAT FAIL THE CHECKS ARE CHECKED FURTHER. C C PROGRAM HISTORY LOG: C 90-11-06 J. WOOLLEN C C USAGE: CALL QCLOOP C C INPUT ARGUMENTS: NONE C C OUTPUT ARGUMENTS: NONE C C SUBPROGRAMS CALLED: C UNIQUE: - QCPTS C C LIBRARY: C COMMON - ABORT C C REMARKS: C C C ATTRIBUTES: C LANGUAGE: CRAY FORTRAN C MACHINE: CRAY C C$$$ SUBROUTINE QCLOOP(QCPTS,MAXTRY) COMMON /FPLIST / NFL,INFL( 300000 ) COMMON /CKLIST / NDD,INDD( 300000 ) C-CRA COMMON /TOLLER / ITER,T,TZ,TSAT,TWND,TJET COMMON /TOLLER / ITER,ITERADD,T,TZ,TSAT,TWND,TJET COMMON /TOSSUNIT/ ITUNIT NAMELIST /TOLLS / TB,TZ,TSAT,TWND,TJET LOGICAL DONE C---------------------------------------------------------------------- C----------------------------------------------------------------------- IF(NDD.EQ.0) THEN PRINT*,'QCLOOP - NDD = 0' RETURN ENDIF ITRY = 1 ITER = 1 IFLIP = 0 C SET UP THE TOLERANCES C --------------------- TB = 4. TZ = 5.0 TSAT = 5.0 TWND = 4.5 TJET = 5.5 REWIND (20) READ(20,TOLLS,END=1) PRINT*,'READ NAMELIST TOLLS' 1 WRITE(6,TOLLS) WRITE(60,TOLLS) TZ = TZ/TB TSAT = TSAT/TB TWND = TWND/TB TJET = TJET/TB T = TB C DO THE QCOI CHECKS ON LIST IN COMMON /CKLIST/ C --------------------------------------------- 10 REWIND(ITUNIT) CALL QCPTS C SET OB ERRS TO REFLECT LAST ROUND OF CHECKS - RECORD FLIPS C ---------------------------------------------------------- NBB = 0 NCK = 0 NFG = 0 NFB = 0 NFL = 0 DO 20 N=1,NDD ILEV = INDD(N) OBER = FOE(ILEV) IF(OBER.NE.0) THEN IF(FOI(ILEV).EQ.1) THEN PACK = WOE(ILEV,ABS(OBER)) IF(OBER.LT.0) THEN NFG = NFG + 1 NFL = NFL + 1 INFL(NFL) = ILEV ENDIF ELSE IF(FOI(ILEV).EQ.-1) THEN PACK = WOE(ILEV,-ABS(OBER)) IF(OBER.GT.0) THEN NFB = NFB + 1 NFL = NFL + 1 INFL(NFL) = ILEV ENDIF NBB = NBB + 1 ENDIF NCK = NCK + 1 ENDIF 20 CONTINUE C CHECK STATUS AND WRITE ITERATION REPORT C --------------------------------------- DONE = NFL.EQ.0 .OR. ITRY.EQ.MAXTRY WRITE(6 ,30)ITRY,NFG,NFB,NFL,NBB,NCK WRITE(60,30)ITRY,NFG,NFB,NFL,NBB,NCK 30 FORMAT('ITERATION ',I2,' : NFLIP=',3I6,' NBAD=',I6,' NTOT=',I6) C FINISHED YET? C ------------- IF(.NOT.DONE) THEN ITRY = ITRY+1 IF(ITRY.LT.MAXTRY) THEN ITER = ITRY GOTO 10 ELSE ITER = ITRY DO 40 N=1,NFL PACK = WOE(INFL(N),-ABS(FOE(INFL(N)))) 40 CONTINUE GOTO 10 ENDIF ENDIF C FINAL REPORT AND AWAY C --------------------- NBAD = NBB WRITE(6 ,*)' TOTAL BAD OBS=',NBAD WRITE(60,*)' TOTAL BAD OBS=',NBAD ITRY = 0 ITER = 0 RETURN END C$$$ SUBPROGRAM DOCUMENTATION BLOCK C . . . . C SUBPROGRAM: QCOI THE QC ANALYSIS KERNAL C PRGMMR: DEAVEN ORG: W/NMC22 DATE: 88-02-05 C PRGMMR: WOOLLEN ORG: W/NMC22 DATE: 90-03-24 C C ABSTRACT: COMPUTES AND SETS UP THE COEFFICIENT MATRIX AND RIGHT C HAND SIDES. SOLVES THE MATRIX PROBLEMS. APPLIES THE COMPUTED C WEIGHTS TO THE RESIDUALS TO OBTAIN THE ANALYSIS INCREMENTS AND C ANALYSIS ERRORS. C C USAGE: C INPUT ARGUMENTS: IN COMMON /HVECT/ C NXXYY - NUMBER OF OBS IN PACKAGE C XX - LONGITUDES OF OBS IN PACKAGE C YY - LATITUDES OF OBS IN PACKAGE C PP - PRESSURES OF OBS IN PACKAGE C II - INDEXES OF OBS IN PACKAGE C LL - REPORT LEVEL OF OBS IN PACKAGE C FESD - FORECAST ERRORS AT OB LOCATIONS C MKOBS - LIST OF OB INDEXES USED IN EACH INTERPOLATION C NDIM - NUMBER OF OBS USED IN EACH INTERPOLATION C JVMAX - NUMBER OF INTERPOLATING OBS FOR EACH CHECK IN C EACH PACKAGE C C OUTPUT ARGUMENTS: IN COMMON /HVECT/ C AA - INTERPOLATED VALUE OF OBS IN PACKAGE C AESD - ANALYSIS ERRORS AT OB LOCATIONS C C SUBPROGRAMS CALLED: C UNIQUE: - COORS,DRCTSL C C LIBRARY: C COMMON - ABORT C C EXIT STATES: C COND = 9 - ERROR DETECTED (SEE CONSOLE OUTPUT) C C REMARKS: C C ATTRIBUTES: C LANGUAGE: CRAY FORTRAN C MACHINE: CRAY C C$$$ C----------------------------------------------------------------------- SUBROUTINE QCOI PARAMETER( NQC = 4 , . NMX = 64 , . NN = 10 *( 10 +1)/2 , . N = 10 , . LVEC = 64 * 10 , . LMAT = 64 *NN ) COMMON /CORP/ CHLP(1500),CVC,FFLAT(91),IFA(N,N) C-CRA TASK COMMON/HVECT/NXXYY,MCHK( 64 ,2),MOBS( 64 , 10 , 4 ,2), COMMON/HVECT/NXXYY,NYYZZ,MCHK( 64 ,2),MOBS( 64 , 10 , 4 ,2) 1, . MDIM( 64 , 4 ),AA( 64 , 4 ,2), . AE( 64 , 4 ,2),WT( 64 , 10 , 4 ,2), . FE( 64 ,3) DIMENSION CHLS(LVEC) , XXCP(LVEC) , XXOB(LVEC) , C1(LMAT) , . FLAT(LVEC) , YYCP(LVEC) , YYOB(LVEC) , C2(LMAT) , . TZFA(LVEC) , KKCP(LVEC) , KKOB(LVEC) , DD(LMAT) , . NGRP(LVEC) , PPCP(LVEC) , PPOB(LVEC) , XI(LMAT) , . NIOB(LVEC) , OEOB(LVEC) , YI(LMAT) , . OBOB(LVEC) , KI(LMAT) , . PI(LMAT) , . XJ(LMAT) , . YJ(LMAT) , . KJ(LMAT) , . PJ(LMAT) , . CH(LMAT) , . FL(LMAT) , . TZ(LMAT) , . NG(LMAT) , . NI(LMAT) , . NJ(LMAT) DIMENSION FA(NMX,NN), RA(NMX,N,2), DP(NMX,2), IVAR(3) DATA IVAR/2,1,2/ C----------------------------------------------------------------------- C----------------------------------------------------------------------- IF(NXXYY.EQ.0) RETURN C-CRA AA = 0. C-CRA AE = 0. C-CRA WT = 0. C . MDIM( 64 , 4 ),AA( 64 , 4 ,2), C . AE( 64 , 4 ,2),WT( 64 , 10 , 4 ,2), DO I=1, 64 * 4 *2 AA(I,1,1) = 0. AE(I,1,1) = 0. ENDDO DO I=1, 64 * 10 * 4 *2 WT(I,1,1,1) = 0. ENDDO C LOOP THRU EACH TYPE OF QC CHECK C ------------------------------- DO 100 IQ=1,NQC IVMAX = 0 C CLEAR THE MATRIX ARRAYS C ----------------------- C-CRA FA = 0. C-CRA RA = 0. C-CRA DP = 0. C DIMENSION FA(NMX,NN), RA(NMX,N,2), DP(NMX,2), IVAR(3) DO I=1,NMX*NN FA(I,1) = 0. ENDDO DO I=1,NMX*N*2 RA(I,1,1) = 0. ENDDO DO I=1,NMX*2 DP(I,1) = 0. ENDDO C SET UP VECTORS WITH INFORMATION ABOUT EACH OB TO BE CHECKED C ----------------------------------------------------------- DO 15 IGRP=1,NXXYY NDIM = MDIM(IGRP,IQ) IF(NDIM.EQ.0) GOTO 15 ILEV = MCHK(IGRP,1) IREP = FHD(ILEV) TZF = FTZ(ILEV) XXC = FLN(IREP) YYC = FLT(IREP) PPC = FPR(ILEV) KKC = MCHK(IGRP,2) CHL = CHLP(INT(MIN(PPC,1500.))) FLA = FFLAT(INT(ABS(YYC)+1.51)) PPC = LOG(PPC) DO 10 IOB=1,NDIM JLEV = MOBS(IGRP,IOB,IQ,1) JTYP = MOBS(IGRP,IOB,IQ,2) JREP = FHD(JLEV) IVMAX = IVMAX+1 CHLS(IVMAX) = CHL FLAT(IVMAX) = FLA TZFA(IVMAX) = TZF XXCP(IVMAX) = XXC YYCP(IVMAX) = YYC KKCP(IVMAX) = KKC PPCP(IVMAX) = PPC XXOB(IVMAX) = FLN(JREP) YYOB(IVMAX) = FLT(JREP) PPOB(IVMAX) = LOG(FPR(JLEV)) KKOB(IVMAX) = JTYP OEOB(IVMAX) = FOE(JLEV) OBOB(IVMAX) = FOB(JLEV,IVAR(JTYP)) NGRP(IVMAX) = IGRP NIOB(IVMAX) = IOB 10 CONTINUE 15 CONTINUE IF(IVMAX.EQ.0) GOTO 61 C COMPUTE THE RIGHT HAND SIDES C ---------------------------- CALL COOR2(IVMAX,CHLS,FLAT,TZFA,XXOB,YYOB,KKOB,XXCP,YYCP,KKCP, . DD,C1,C2) DO 20 I=1,IVMAX IGRP = NGRP(I) IOB = NIOB(I) VCOR = 1./(1.+CVC*(PPOB(I)-PPCP(I))**2) RA(IGRP,IOB,1) = C1(I)*VCOR RA(IGRP,IOB,2) = C2(I)*VCOR 20 CONTINUE C COMPUTE THE OFF-DIAGONAL MATRIX TERMS C ------------------------------------- K = 0 L = 0 DO 40 IGRP=1,NXXYY NDIM = MDIM(IGRP,IQ) DO 30 IOB = 1,NDIM-1 DO 30 JOB = IOB+1,NDIM I = IOB+K J = JOB+K L = L+1 XI(L) = XXOB(I) YI(L) = YYOB(I) KI(L) = KKOB(I) PI(L) = PPOB(I) XJ(L) = XXOB(J) YJ(L) = YYOB(J) KJ(L) = KKOB(J) PJ(L) = PPOB(J) CH(L) = CHLS(I) FL(L) = FLAT(I) TZ(L) = TZFA(I) NG(L) = NGRP(I) NI(L) = NIOB(I) NJ(L) = NIOB(J) 30 CONTINUE K = K+NDIM 40 CONTINUE CALL COOR1(L,CH,FL,TZ,XI,YI,KI,XJ,YJ,KJ,DD,C1,C2) C STORE OFF-DIAGONAL MATRIX TERMS IN THE SOLVER ARRAY C --------------------------------------------------- DO 50 I=1,L IGRP = NG(I) IOB = IFA(NI(I),NJ(I)) FA(IGRP,IOB) = C1(I)/(1.+CVC*(PI(I)-PJ(I))**2) 50 CONTINUE C COMPUTE AND STORE THE ON-DIAGONAL MATRIX TERMS C ---------------------------------------------- DO 55 I=1,IVMAX IGRP = NGRP(I) IOB = IFA(NIOB(I),NIOB(I)) FA(IGRP,IOB) = 1. + (OEOB(I)/FE(IGRP,KKOB(I)))**2 55 CONTINUE C SOLVE THE LINEAR INTERPOLATION EQUATIONS C ---------------------------------------- CALL DRCTSL(FA,RA,DP,MDIM(1,IQ),NXXYY) C COMPUTE AND STORE ANALYSIS INCREMENTS AND ERRORS C ------------------------------------------------ DO 60 I=1,IVMAX IGRP = NGRP(I) IOB = NIOB(I) FERAT = FE(IGRP,KKCP(I))/FE(IGRP,KKOB(I)) WT(IGRP,IOB,IQ,1) = RA(IGRP,IOB,1) * FERAT WT(IGRP,IOB,IQ,2) = RA(IGRP,IOB,2) * FERAT AA(IGRP,IQ,1) = AA(IGRP,IQ,1) + OBOB(I)*WT(IGRP,IOB,IQ,1) AA(IGRP,IQ,2) = AA(IGRP,IQ,2) + OBOB(I)*WT(IGRP,IOB,IQ,2) 60 CONTINUE 61 DO 65 I=1,NXXYY AE(I,IQ,1) = FE(I,MCHK(I,2)) * SQRT(1.-DP(I,1)) AE(I,IQ,2) = FE(I,MCHK(I,2)) * SQRT(1.-DP(I,2)) 65 CONTINUE C END OF LOOP OVER EACH QC CHECK C ------------------------------ 100 CONTINUE RETURN END C$$$ SUBPROGRAM DOCUMENTATION BLOCK C C SUBPROGRAM: QCPTS C PRGMMR: WOOLLEN ORG: NMC22 DATE: 90-11-06 C C ABSTRACT: PUTS TOGETHER PACKAGES OF OBSERVATIONS TO BE CHECKED C IN A GIVEN ITERATION AND INVOKES AUTOTASKING OF THE DATA C SELECTION, OI, AND DECISION MAKING PROCEDURES. C C PROGRAM HISTORY LOG: C 90-11-06 J. WOOLLEN C C USAGE: C INPUT ARGUMENTS: IN COMMON /CKLIST/ C NDD - NUMBER OF OBSERVATIONS TO CHECK C INDD - ARRAY OF OBSERVATION INDEXES TO CHECK C C OUTPUT ARGUMENTS: NONE C C SUBPROGRAMS CALLED: C UNIQUE: - SELDAT,QCOI,OUTPUT C C LIBRARY: C COMMON - ABORT C C REMARKS: C C C ATTRIBUTES: C LANGUAGE: CRAY FORTRAN C MACHINE: CRAY C C$$$ SUBROUTINE QCPTS COMMON /CKLIST/ NDD,INDD( 300000 ) DIMENSION IND( 5000 ),NXY( 5000 ) DATA NMX / 64 / C---------------------------------------------------------------------- C---------------------------------------------------------------------- C ASSEMBLE LISTS OF OBS FOR MULTI PROCESSING C ------------------------------------------ IND(1) = 1 IPK = 1 C-CRA NXY = 0 C DIMENSION IND( 5000 ),NXY( 5000 ) DO I=1, 5000 NXY(I) = 0 ENDDO DO 5 I=1,NDD ILEV = INDD(I) IF(NXY(IPK)+1.GT.NMX) THEN IF(IPK+1.GT. 5000 ) CALL SABORT('QCPTS - IPK > MAXPAK') IPK = IPK+1 IND(IPK) = I NXY(IPK) = 1 ELSE NXY(IPK) = NXY(IPK)+1 ENDIF 5 CONTINUE C MULTI-PROCESS QC OB LISTS C ------------------------- CMIC$ DOALL AUTOSCOPE DO 10 I=1,IPK CALL SELDAT(IND(I),NXY(I)) CALL QCOI CALL OUTPUT 10 CONTINUE RETURN END C$$$ SUBPROGRAM DOCUMENTATION BLOCK C . . . . C SUBPROGRAM: RDOBER ASSIGN OBSERVATIONAL ERRORS C PRGMMR: DEAVEN ORG: W/NMC22 DATE: 88-02-10 C C ABSTRACT: INITIALIZES ARRAY 'QIKE' WITH OBSERVATION ERRORS FOR C ALL KINDS OF OBSERVATION TYPES LISTED IN THE OBSERVATION ERROR C FILE. C C PROGRAM HISTORY LOG: C 86-03-21 G. DIMEGO C 88-11-24 D. DEAVEN RECODED FOR CYBER 205 C 92-07-29 J. WOOLLEN UPDATED FOR OIQC C C USAGE: CALL RDOBER C C INPUT ARGUMENTS: C LU - LOGICAL UNIT FOR OB ERR INPUT C C ATTRIBUTES: C LANGUAGE: CDC FORTRAN200 C MACHINE: CYBER C C$$$ SUBROUTINE RDOBER(LU) COMMON /OBERRS/ QIKE(21,300),PMAND(21),P50(21),CHL(21) COMMON /OBDESC/ DESC(300) COMMON /GETOE/ USEOE CHARACTER*80 CKX,DESC DIMENSION RAWERR(21) LOGICAL USEOE DATA RAWERR / 7.0 , 7.5 , 8.5 , 11.0 , 13.7 , 15.2 , 16.0 , . 16.1 , 17.5 , 19.1 , 20.5 , 24.0 , 24.0 , 24.0 , . 24.0 , 24.0 , 24.0 , 24.0 , 24.0 , 24.0 , 24.0 / DATA DUMMY/0./ C---------------------------------------------------------------------- C---------------------------------------------------------------------- C DEFAULT VALUES FOR OB DESCRIPTIONS AND ERRORS C --------------------------------------------- DO 1 I=1,300 DESC(I) = ' ' DO 1 L=1,21 QIKE(L,I) = 0. IF(I.EQ.120) QIKE(L,I) = RAWERR(L) 1 CONTINUE USEOE = .FALSE. C READ OB ERROR VALUES FROM OB ERROR FILE C --------------------------------------- 10 READ(LU,'(A80)',END=20) CKX IF(CKX(1:1).EQ.' ') GOTO 10 READ(CKX,'(I3)',ERR=15) KX IF(KX.GT.0 .AND. KX.LT.300) THEN READ(LU,'(7F8.2)',END=20) (QIKE(L,KX),L=1,21) DESC(KX) = CKX USEOE = .TRUE. GOTO 10 ENDIF 15 PRINT*,'******************************************************' PRINT*,'******* UNRECOGNISED RECORD IN OBERRS FILE *********' PRINT*, CKX PRINT*,'******************************************************' GOTO 10 20 RETURN END C$$$ SUBPROGRAM DOCUMENTATION BLOCK C C SUBPROGRAM: SATPTS C PRGMMR: WOOLLEN ORG: NMC22 DATE: 90-11-06 C C ABSTRACT: PUTS TOGETHER PACKAGES OF OBSERVATIONS TO BE CHECKED C IN A GIVEN ITERATION AND INVOKES AUTOTASKING OF THE DATA C SELECTION, OI, AND DECISION MAKING PROCEDURES. C C PROGRAM HISTORY LOG: C 90-11-06 J. WOOLLEN C C USAGE: C INPUT ARGUMENTS: IN COMMON /CKLIST/ C NDD - NUMBER OF OBSERVATIONS TO CHECK C INDD - ARRAY OF OBSERVATION INDEXES TO CHECK C C OUTPUT ARGUMENTS: NONE C C SUBPROGRAMS CALLED: C UNIQUE: - SELDAT,QCOI,SATPUT C C LIBRARY: C COMMON - ABORT C C REMARKS: C C C ATTRIBUTES: C LANGUAGE: CRAY FORTRAN C MACHINE: CRAY C C$$$ SUBROUTINE SATPTS COMMON /CKLIST/ NDD,INDD( 300000 ) DIMENSION IND( 5000 ),NXY( 5000 ) DATA NMX / 64 / C---------------------------------------------------------------------- C---------------------------------------------------------------------- C ASSEMBLE LISTS OF OBS FOR MULTI PROCESSING C ------------------------------------------ IND(1) = 1 C-CRA NXY = 0 C DIMENSION IND( 5000 ),NXY( 5000 ) DO I=1, 5000 NXY(I) = 0 ENDDO IPK = 1 NPK = 0 LREP = 0 DO 5 I=1,NDD ILEV = INDD(I) IREP = FHD(ILEV) IF(IREP.NE.LREP) THEN LREP = IREP KLEV = FNL(IREP) IF(NXY(IPK)+KLEV.GT.NMX) THEN IF(IPK+1.GT. 5000 ) CALL SABORT('QCPTS - IPK > MAXPAK') IPK = IPK+1 IND(IPK) = I NXY(IPK) = KLEV ELSE NXY(IPK) = NXY(IPK) + KLEV ENDIF ENDIF 5 CONTINUE C MULTI-PROCESS QC OB LISTS C ------------------------- CMIC$ DOALL AUTOSCOPE DO 10 I=1,IPK CALL SELDAT(IND(I),NXY(I)) CALL QCOI CALL SATPUT 10 CONTINUE RETURN END C$$$ SUBPROGRAM DOCUMENTATION BLOCK C C SUBPROGRAM: SATPUT C PRGMMR: WOOLLEN ORG: NMC22 DATE: 90-11-06 C C ABSTRACT: EXAMINES OUTPUTS FROM OI QUALITY CONTROL INTERPOLATION C AND MAKES DECISION WHETHER A SATELLITE PASSES OR FAILS THE CHECKS. C DUE TO SPECIFIC CHARACTERISTIC DIFFERENCES OF SATELLITE C RETRIEVALS FROM CONVENTIONAL 'IN-SITU' OBSERVATIONS, A SPECIAL C ALGORITHM IS EMPLOYED IN JUDGING THE OUTCOMES OF THE QC CHECKS. C ALSO PROVIDES FOR A DETAILED PRINT TRACING ALL DECISIONS MADE C IN REGARD TO SPECIFIED OBSERVATIONS IN THE QC PROCESS. C C PROGRAM HISTORY LOG: C 90-11-06 J. WOOLLEN C C USAGE: C INPUT ARGUMENTS: IN COMMON /HVECT/ C NXXYY - NUMBER OF OBS IN PACKAGE C XX - LONGITUDES OF OBS IN PACKAGE C YY - LATITUDES OF OBS IN PACKAGE C PP - PRESSURES OF OBS IN PACKAGE C II - INDEXES OF OBS IN PACKAGE C LL - REPORT LEVEL OF OBS IN PACKAGE C FESD - FORECAST ERRORS AT OB LOCATIONS C AA - INTERPOLATED VALUE OF OBS IN PACKAGE C AESD - ANALYSIS ERRORS AT OB LOCATIONS C MKOBS - LIST OF OB INDEXES USED IN EACH INTERPOLATION C NDIM - NUMBER OF OBS USED IN EACH INTERPOLATION C JVMAX - NUMBER OF INTERPOLATING OBS FOR EACH CHECK IN C EACH PACKAGE C C OUTPUT FILES: C FORT.60 - TOSS.OUT - DETAILED TOSSLIST C FORT.71 - OBPRT1 - OB DETAILED PRINTOUTS FOR QC DIAGNOSIS C FORT.72 - OBPRT2 - OB DETAILED PRINTOUTS FOR QC DIAGNOSIS C FORT.73 - OBPRT3 - OB DETAILED PRINTOUTS FOR QC DIAGNOSIS C FORT.74 - OBPRT4 - OB DETAILED PRINTOUTS FOR QC DIAGNOSIS C FORT.75 - OBPRT5 - OB DETAILED PRINTOUTS FOR QC DIAGNOSIS C FORT.76 - OBPRT6 - OB DETAILED PRINTOUTS FOR QC DIAGNOSIS C FORT.77 - OBPRT7 - OB DETAILED PRINTOUTS FOR QC DIAGNOSIS C FORT.78 - OBPRT8 - OB DETAILED PRINTOUTS FOR QC DIAGNOSIS C FORT.79 - OBPRT9 - OB DETAILED PRINTOUTS FOR QC DIAGNOSIS C FORT.80 - OBPRT10 - OB DETAILED PRINTOUTS FOR QC DIAGNOSIS C C OUTPUT ARGUMENTS: IN COMMON /OBMARK/ C OBIND - CONTAINS QC OUTCOME FOR EACH REPORT AND LEVEL: C -2 = BELOW MODEL SURFACE (SET IN SUB PSFILE) C -1 = FAILED OI CHECKS C 0 = NO DATA IN REPORT C 1 = PASSED OI CHECKS C C SUBPROGRAMS CALLED: C LIBRARY: C COMMON - ABORT C C EXIT STATES: C COND = 135 - ERROR DETECTED (SEE CONSOLE OUTPUT) C C REMARKS: C C C ATTRIBUTES: C LANGUAGE: CRAY FORTRAN C MACHINE: CRAY C C$$$ SUBROUTINE SATPUT C-CRA TASK COMMON/HVECT/NXXYY,MCHK( 64 ,2),MOBS( 64 , 10 , 4 ,2), COMMON/HVECT/NXXYY,NYYZZ,MCHK( 64 ,2),MOBS( 64 , 10 , 4 ,2) 1, . MDIM( 64 , 4 ),AA( 64 , 4 ,2), . AE( 64 , 4 ,2),WTS( 64 , 10 , 4 ,2), . FE( 64 ,3) COMMON /OBS/ NREP,NLEV,SID( 200000 ),INOB( 200000 ),IMAP(360,181) COMMON /CORP / CHLP(1500),CVC,FFLAT(91),IFA( 10 , 10 ) COMMON /SATCOR / EPSCOR(20) C-CRA COMMON /TOLLER / ITER,TB,TZ,TSAT,TWND,TJET COMMON /TOLLER / ITER,ITERADD,TB,TZ,TSAT,TWND,TJET COMMON /TOSSUNIT/ ITUNIT DIMENSION OMA1(50),OMA2(50),EPV1(50),EPV2(50),QCT(50),QCS(50) DIMENSION ASD(50,2) DIMENSION RAF(50,2) DIMENSION OMA(50,2) DIMENSION LEV(50,2) DIMENSION OSD(50,2) DIMENSION FSD(50,2) DIMENSION PRS(50,2) DIMENSION NSEG(2) CHARACTER*8 SEG(2) DATA SEG/'1000-100','099 - 00'/ c DATA EPSCOR/ .99, .88,. 61, .36, .18, .05,-.05,-.12,-.15,-.14, c . -.11,-.08,-.05,-.03,-.02,-.02,-.01,-.01, 0., 0./ LOGICAL LAST DATA PI180/.017453/, MAXNSEG/50/ C---------------------------------------------------------------------- C---------------------------------------------------------------------- LLEV = MCHK(1,1) LREP = FHD(LLEV) C-CRA NSEG = 0 C DIMENSION NSEG(2),SEG(2) DO I=1,2 NSEG(I) = 0 ENDDO C LOOP OVER EACH SAT LEVEL IN THIS GROUP C -------------------------------------- DO 50 I=1,NXXYY LAST = I.EQ.NXXYY ILEV = MCHK(I,1) IREP = FHD(ILEV) C ITEMIZE COMPLETE PROFILES INTO SEGMENTS C --------------------------------------- 10 IF(IREP.EQ.LREP) THEN PRES = FPR(ILEV) OBER = FOE(ILEV) OMF = FOB(ILEV,2) IF(PRES.LE.2000. .AND. PRES.GE.100.) THEN K = 1 ELSE IF(PRES.LT.100. .AND. PRES.GT.0.) THEN K = 2 ENDIF NSEG(K) = NSEG(K) + 1 IF(NSEG(K).GT.MAXNSEG) 1 CALL SABORT('SATPUT - MAX NSEG EXCEEDED') J = NSEG(K) C STORE THE RESULTS OF THE CHECKS C ------------------------------- ASD(J,K) = 0. OMA(J,K) = 0. SWT = 0. DO 20 L=1,3 AER = AE(I,L,1) FER = FE(I,1) AMF = AA(I,L,1) WT = MAX(.01,1.-(AER/FER)**2) ASD(J,K) = ASD(J,K) + WT*AER**2 OMA(J,K) = OMA(J,K) + (OMF-AMF)*WT SWT = SWT + WT 20 CONTINUE ASD(J,K) = SQRT(ASD(J,K)/SWT) OMA(J,K) = OMA(J,K)/SWT RAF(J,K) = ASD(J,K)/FER LEV(J,K) = ILEV FSD(J,K) = FER OSD(J,K) = OBER PRS(J,K) = PRES IF(.NOT.LAST) GOTO 50 ENDIF C A REPORT IS ITEMIZED AND READY TO JUDGE C --------------------------------------- STAID = SID(LREP) CLAT = FLT(LREP) CLON = FLN(LREP) KX = FKX(LREP) LREP = IREP C COMPUTE THE TOLERANCE VALUE C --------------------------- SINLAT = SIN(ABS(CLAT)*PI180) TT = TB*TSAT*SINLAT TMIN = 2.0 TT = MAX(TMIN,TT) C SUM THE RESULTS OVER THE PROFILE SEGMENTS - MAKE AND STORE DECISION C ------------------------------------------------------------------- DO 35 K=1,2 IF (NSEG(K).EQ.0) GOTO 35 QC1 = 0. QC2 = 0. C-CRA QCT = 0. C-CRA QCS = 0. C-CRA OMA1 = 0. C-CRA OMA2 = 0. C-CRA EPV1 = 0. C-CRA EPV2 = 0. C DIMENSION OMA1(50),OMA2(50),EPV1(50),EPV2(50),QCT(50),QCS(50) DO L=1,50 QCT(L) = 0. QCS(L) = 0. OMA1(L) = 0. OMA2(L) = 0. EPV1(L) = 0. EPV2(L) = 0. ENDDO TADJ = 0. DO 30 J=1,NSEG(K) TADJ = TADJ + MAX(1.,(OSD(J,K)/FSD(J,K))**2) DT = OMA(J,K)**2 EPV = (ASD(J,K)**2+OSD(J,K)**2) OMA1(J) = DT EPV1(J) = EPV QCT(J) = SQRT(DT/EPV) QC1 = MAX(QC1,QCT(J)) IF(J.GT.1) THEN RPLN = LOG(PRS(J-1,K)/PRS(J,K)) PL=20.*RPLN + 1.001 IF(PL.LT.20.) THEN KP=PL XP=PL-KP OBCOR = (1.-XP)*EPSCOR(KP) + XP*EPSCOR(KP+1) ELSE OBCOR = 0. ENDIF COVO = OSD(J-1,K)*OSD(J,K)*OBCOR FCCOR = 1./(1.+CVC*RPLN**2) COVA = ASD(J-1,K)*ASD(J,K)*FCCOR STAB = (OMA(J-1,K)-OMA(J,K))**2 ESTAB = EPV+EPV0 - 2.*(COVO+COVA) OMA2(J-1) = STAB EPV2(J-1) = ESTAB QCS(J-1) = SQRT(STAB/ESTAB) QC2 = MAX(QC2,QCS(J-1)) ENDIF EPV0 = EPV 30 CONTINUE TADJ = SQRT(TADJ/NSEG(K)) TTT = TT/TADJ QC = MAX(QC1,QC2) IQC = SIGN(1.,TTT-QC) IF(IQC.EQ.-1) THEN WRITE(ITUNIT,31)STAID,KX,CLON,CLAT,SEG(K),QC1,QC2,QC,TTT 31 FORMAT(A6,I4,2F8.2,A8,2x,4F8.2) ENDIF DO 32 J=1,NSEG(K) KIQC = IQC IF(CLAT.LT.0. .OR. PRS(J,K).LE.100.) THEN QC = MAX(QCT(J),QCS(J)) KIQC = SIGN(1.,TTT-QC) ENDIF QCFLG = WOI(LEV(J,K),KIQC) 32 CONTINUE 35 CONTINUE C WRAP BACK TO START THE NEXT PROFILE C ----------------------------------- IF(.NOT.LAST) THEN C-CRA NSEG = 0 C DIMENSION NSEG(2),SEG(2) DO L=1,2 NSEG(L) = 0 ENDDO GOTO 10 ENDIF C END OF LOOP OVER LEVELS C ----------------------- 50 CONTINUE RETURN END C----------------------------------------------------------------------- C----------------------------------------------------------------------- SUBROUTINE SEARCH PARAMETER (MAXLV=255) PARAMETER (MTOTH=10 ) PARAMETER (MTOTV=5 ) COMMON /OBS/ NREP,NLEV,SID( 200000 ),INOB( 200000 ),IMAP(360,181) DIMENSION ITW(MAXLV,MTOTH,2),IPF(MAXLV,MTOTV,2) DATA DUMMY /0/ C AUTOTASKING LOOP FOR DATA SEARCH C -------------------------------- CMIC$ DOALL AUTOSCOPE PRIVATE(ITW,IPF,ILV,NLV) DO NPTR=1,NREP CALL SEARCHS(NPTR,ITW,IPF) CMIC$ GUARD 0 ILV = FLV(NPTR)-1 NLV = FNL(NPTR) DO K=1,2 DO J=1,MAX(MTOTH,MTOTV) DO L=1,NLV IF(J.LE.MTOTH) PACK = WTW(ILV+L,J,K,ITW(L,J,K)) IF(J.LE.MTOTV) PACK = WPF(ILV+L,J,K,IPF(L,J,K)) ENDDO ENDDO ENDDO CMIC$ END GUARD 0 ENDDO RETURN END C----------------------------------------------------------------------- C----------------------------------------------------------------------- SUBROUTINE SEARCHS(NPTR,ITW,IPF) PARAMETER (MAXLV=255) PARAMETER (MTOTH=10 ) PARAMETER (MTOTV=5 ) DIMENSION ITW(MAXLV,MTOTH,2),IPF(MAXLV,MTOTV,2) COMMON /OBS / NREP,NLEV,SID( 200000 ),INOB( 200000 ),IMAP(360,181) COMMON /CORP/ CHLP(1500),CVC,FFLAT(91),IFA( 10 , 10 ) CHARACTER*8 SID CHARACTER*3 TAG,FOC DIMENSION NIND( 200000 ),RLAT( 200000 ),RLON( 200000 ), . RDIS( 200000 ),RDIR( 200000 ), . DIST( 50000 ),DIRN( 50000 ), . ICOR( 50000 ),IDIR( 50000 ), . LIND( 50000 ),LTYP( 50000 ),PRES( 50000 ), . INL(MAXLV),INI(MAXLV), . SCAT( 50000 ,4,4),CAT( 50000 ,4,4),NONE(4,4), . IPIC( 10 ),CATCOR(101) LOGICAL SCAT,CAT,NONE,SAME,MORE DATA CATCOR /21*5.,30*4.,20*3.,20*2.,10*1./ DATA PI180 /.0174532/, RADE /6371./, RSCAN /1000./ DATA MAXCYL / 50000 / , MAXREP / 200000 / C---------------------------------------------------------------------- C---------------------------------------------------------------------- C INITIALIZE THE RESULTS FOR THE ZERO CASE C ---------------------------------------- C-CRA ITW = 0 C-CRA IPF = 0 C DIMENSION ITW(MAXLV,MTOTH,2),IPF(MAXLV,MTOTV,2) DO I=1,MAXLV*MTOTH*2 ITW(I,1,1) = 0 IPF(I,1,1) = 0 ENDDO IF(FQL(NPTR).LE.0) RETURN C DO THE DATA SEARCH IF SO INSTRUCTED C ----------------------------------- BLON = FLN(NPTR) BLAT = FLT(NPTR) BINQ = FQL(NPTR) NBKX = FKX(NPTR) TAG = FOC(NPTR) C MAKE SPECIAL CASE ADJUSTMENTS (RECCOS,SATEMP) C --------------------------------------------- IF(MOD(NBKX,100).EQ.32) NC = 3 IF(MOD(NBKX,100).NE.32) NC = 8 IF(TAG.EQ.'SAT') BINQ = BINQ - 1. C OUTLINE THE SEARCH AREA C ----------------------- DELLAT = RSCAN/(PI180*RADE) COSLAT = COS(BLAT*PI180) DELLON = MIN(DELLAT/COSLAT,180.0) X1 = BLON-DELLON X2 = BLON+DELLON Y1 = BLAT-DELLAT Y2 = BLAT+DELLAT IF(X1.LT. 0.) X1 = X1 + 360. IF(X2.GE.360.) X2 = X2 - 360. IF(Y1.LT.-90.) Y1 = -90. IF(Y2.GT. 90.) Y2 = 90. C ITEMIZE REPORTS IN THE SEACH CYLINDER C ------------------------------------- NPRF = 0 NOBS = 0 I1 = X1 + 1. I2 = X2 + 1. J1 = Y1 + 91. J2 = Y2 + 91. DO 10 K=1,2 IA = I1 IB = I2 IF(I1.GT.I2) THEN IF(K.EQ.1) IB = 360 IF(K.EQ.2) IA = 1 ELSE IF(K.EQ.2) GOTO 11 ENDIF DO 10 J=J1,J2 DO 10 N=IMAP(IA,J),IMAP(IB,J) NN = INOB(N) QLA = FQL(NN) IF(QLA.GE.0. .AND. QLA.LE.BINQ) THEN NPRF = NPRF+1 NIND(NPRF) = NN RLON(NPRF) = FLN(NN) RLAT(NPRF) = FLT(NN) ENDIF 10 CONTINUE 11 IF(NPRF.GT.MAXREP) CALL SABORT('SEARCH - NPRF > MAXREP') C COMPUTE DISTANCE AND DIRECTION FROM THE CYLINDER CENTER C ------------------------------------------------------ CALL CHDIST(BLON,BLAT,RLON,RLAT,RDIS,RDIR,NPRF) C STORE ALL LEVELS OF EACH REPORT LOCATED IN THE CYLINDER C ------------------------------------------------------- DO 21 N=1,NPRF SAME = SID(NPTR)(1:NC) .EQ. SID(NIND(N))(1:NC) IF(RDIS(N).LE.RSCAN .AND. .NOT.SAME) THEN ILEV = FLV(NIND(N)) KLEV = ILEV+FNL(NIND(N))-1 KX = FKX(NIND(N)) DO 20 JLEV=ILEV,KLEV IF(FOI(JLEV).GT.0) THEN NOBS = NOBS+1 LIND(NOBS) = JLEV PRES(NOBS) = LOG(FPR(JLEV)) DIST(NOBS) = (RDIS(N)/RADE)**2 DIRN(NOBS) = AMOD(RDIR(N)/90.,4.) + 1.01 LTYP(NOBS) = KX/100 ENDIF 20 CONTINUE IF(NOBS.GT.MAXCYL) CALL SABORT('SEARCH - NOBS>MAXCYL') ENDIF 21 CONTINUE C LOOP OVER THE OBSERVATION LEVELS TO BE CHECKED C ---------------------------------------------- ILEV = FLV(NPTR) NLV = FNL(NPTR) DO 50 L=1,NLV IL = ILEV-1+L C STORE DATA FOR VERTICAL CHECK (UPA ONLY) C ---------------------------------------- IF(TAG.EQ.'UPA') THEN IF(L.EQ.1) CALL SRTPRS(ILEV,NLV,INL,INI) DO J=1,MTOTV IF(INL(L)-J.GT.0 ) IPF(L,J,1) = INI(INL(L)-J) IF(INL(L)+J.LE.NLV) IPF(L,J,2) = INI(INL(L)+J) ENDDO ENDIF IF(NOBS.EQ.0) GOTO 50 C USE 3D ZZ CORRELATIONS AND DIRECTION ANGLES TO SELECT HORZ DATA C --------------------------------------------------------------- BPRS = FPR(IL) IF(BPRS.LE.0 .OR. BPRS.GT.2000) GOTO 50 IPRS = BPRS BPRS = LOG(BPRS) DO 24 N=1,NOBS CHRZ = EXP(-CHLP(IPRS)*(DIST(N))) CVRT = 1./(1.+ CVC*(BPRS-PRES(N))**2) IC = CHRZ*CVRT*100.+1. ICOR(N) = CATCOR(IC) IDIR(N) = DIRN(N) 24 CONTINUE DO 25 J=1,4 DO 25 I=1,4 DO 25 N=1,NOBS CAT(N,I,J) = ICOR(N).EQ.J .AND. IDIR(N).EQ.I 25 CONTINUE C LOOP OVER DIFFERENT CHECKS FOR THE SAME OB C ------------------------------------------ DO 45 IQ=1,2 C HORIZONTAL CHECK USING HEIGHT DATA IF(IQ.EQ.1) THEN DO 26 J=1,4 DO 26 I=1,4 DO 26 N=1,NOBS SCAT(N,I,J) = CAT(N,I,J) .AND. LTYP(N).EQ.1 26 CONTINUE ENDIF C HORIZONTAL CHECK USING WIND DATA IF(IQ.EQ.2) THEN DO 27 J=1,4 DO 27 I=1,4 DO 27 N=1,NOBS SCAT(N,I,J) = CAT(N,I,J) .AND. LTYP(N).EQ.2 27 CONTINUE ENDIF C SET THE NONE ARRAY TO FALSE C --------------------------- DO 28 J=1,4 DO 28 I=1,4 NONE(I,J) = .FALSE. 28 CONTINUE C SELECTION LOOPS GO UNTIL THEY PICK UP ENOUGH OR ALL POSSIBLE DATA C ----------------------------------------------------------------- NTOT = 0 DO 35 J=1,4 30 MORE = .FALSE. DO 34 I=1,4 IF(NTOT.EQ.MTOTH) GOTO 36 IF(NONE(I,J) ) GOTO 34 DO 32 N=1,NOBS IF(SCAT(N,I,J)) THEN NTOT = NTOT+1 IPIC(NTOT) = N SCAT(N,I,J) = .FALSE. MORE = .TRUE. GOTO 34 ENDIF 32 CONTINUE NONE(I,J) = .TRUE. 34 CONTINUE IF(MORE) GOTO 30 35 CONTINUE 36 CONTINUE C STORE THE DATA INDEXES SELECTED FOR HORIZONTAL CHECKS C ----------------------------------------------------- DO 40 N=1,NTOT ITW(L,N,IQ) = LIND(IPIC(N)) 40 CONTINUE 45 CONTINUE 50 CONTINUE RETURN END C----------------------------------------------------------------------- C----------------------------------------------------------------------- SUBROUTINE SRTPRS(ILEV,NLV,INL,INI) C-CRA DIMENSION INL(NLV),INI(NLV),PRS(NLV),SRT(NLV) DIMENSION INL(NLV),INI(NLV) PARAMETER(NLVPM=1000) DIMENSION PRS(NLVPM),SRT(NLVPM) C IF(NLV.GT.NLVPM) THEN PRINT *,'NLV.GT.NLVPM. NLV=',NLV CALL SABORT('NLV.GT.NLVPM ') ENDIF C----------------------------------------------------------------------- C----------------------------------------------------------------------- DO L=1,NLV ILV = ILEV+L-1 INI(L) = ILV PRS(L) = FPR(ILV) SRT(L) = PRS(L) ENDDO C SORT INDEXES IN DESCENDING ORDER OF PRESSURE C -------------------------------------------- DO I=1 ,NLV DO J=I+1,NLV IF(SRT(I).LT.SRT(J)) THEN STMP = SRT(I) SRT(I) = SRT(J) SRT(J) = STMP ITMP = INI(I) INI(I) = INI(J) INI(J) = ITMP ENDIF ENDDO ENDDO C MAKE A SET OF LOOK BACK POINTERS C -------------------------------- DO 10 I=1,NLV INL(I) = I IF(PRS(I).NE.SRT(I)) THEN DO J=1,NLV IF(PRS(I).EQ.SRT(J)) THEN INL(I) = J GOTO 10 ENDIF ENDDO CALL SABORT('SRTPRS - SORTED PRESSURE NO FOUND') ENDIF 10 ENDDO RETURN END C---------------------------------------------------------------------- C---------------------------------------------------------------------- SUBROUTINE SELDAT(IND,NXY) C-CRA TASK COMMON/HVECT/NXXYY,MCHK( 64 ,2),MOBS( 64 , 10 , 4 ,2), COMMON/HVECT/NXXYY,NYYZZ,MCHK( 64 ,2),MOBS( 64 , 10 , 4 ,2) 1, . MDIM( 64 , 4 ),AA( 64 , 4 ,2), . AE( 64 , 4 ,2),WT( 64 , 10 , 4 ,2), . FE( 64 ,3) COMMON /CKLIST/ NDD,INDD( 300000 ) DIMENSION MAXPIC(4),MAXN(4) DATA MAXPIC/10,10,5,5/ , MAXN/10,10,1,1/ , NMX/ 64 / C---------------------------------------------------------------------- C---------------------------------------------------------------------- C-CRA MDIM = 0 C . MDIM( 64 , 4 ),AA( 64 , 4 ,2), DO I=1, 64 * 4 MDIM(I,1) = 0 ENDDO M = 0 DO 50 NPTR=IND,IND+NXY-1 I = INDD(NPTR) J = FHD(I) KX = FKX(J)/100 DO 45 K=KX,KX M = M + 1 MCHK(M,1) = I MCHK(M,2) = K FE (M,1) = FFE(I,1) FE (M,2) = FFE(I,2) FE (M,3) = FE (M,2) C PICK DATA FROM T AND W GROUPS C ----------------------------- DO 20 IG=1,2 NPIC = 0 DO 15 IP=1,MAXPIC(IG) IF(NPIC.EQ.MAXN(IG)) GOTO 20 IPIC = FTW(I,IP,IG) IF(IPIC.EQ.0) GOTO 20 IF(FOE(IPIC).GT.0) THEN NPIC = NPIC+1 DO 10 IQ=IG,IG*2-1 MDIM(M,IQ) = NPIC MOBS(M,NPIC,IQ,1) = IPIC MOBS(M,NPIC,IQ,2) = IQ 10 CONTINUE ENDIF 15 CONTINUE 20 CONTINUE C PICK DATA FROM SAME PROFILE C --------------------------- NPIC = 0 DO 30 IG=1,2 DO 25 IP=1,5 IPIC = FPF(I,IP,IG) IF(IPIC.EQ.0) GOTO 30 IF(FOE(IPIC).GT.0) THEN DO 24 KK=K,K*2-1 NPIC = NPIC+1 MDIM(M,4) = NPIC MOBS(M,NPIC,4,1) = IPIC MOBS(M,NPIC,4,2) = KK 24 CONTINUE GOTO 30 ENDIF 25 CONTINUE 30 CONTINUE 45 CONTINUE 50 CONTINUE C CHECK THE NUMBER OF OBS TO STORE C -------------------------------- IF(M.GT.NMX) CALL SABORT('SELDAT - M GT NMX') C STORE THE NUMBER OF OBS TO CHECK C -------------------------------- NXXYY = M RETURN END C$$$ SUBPROGRAM DOCUMENTATION BLOCK C . . . . C SUBPROGRAM: SETCON STORES CONSTANTS C PRGMMR: DEAVEN ORG: W/NMC22 DATE: 88-02-10 C C ABSTRACT: INITIALIZES CONSTANTS AND PARAMETERS USING DATA STATEMENTS. C ALSO READS A FILE OF OBSERVATION ERRORS FROM UNIT 15. C C PROGRAM HISTORY LOG: C 86-03-21 G. DIMEGO C 88-11-24 D. DEAVEN CODED FOR CYBER 205 C 92-07-29 J. WOOLLEN UPDATED VERSION FOR OIQC C C USAGE: CALL SETCON C C ATTRIBUTES: C LANGUAGE: CDC FORTRAN200 C MACHINE: CYBER C C$$$ BLOCK DATA COMMON /SATCOR / EPSCOR(20) DATA EPSCOR/ .99, .88,. 61, .36, .18, .05,-.05,-.12,-.15,-.14, . -.11,-.08,-.05,-.03,-.02,-.02,-.01,-.01, 0., 0./ COMMON /CORP / CHLP(1500),CVC,FFLAT(91),IFA( 10 , 10 ) COMMON /FCRATS/ RATION(21),RATIOS(21),TROPUV(21) COMMON /OBERRS/ QIKE(21,300),PMAND(21),P50(21),CHL(21) DATA RATION / .7 , .7 , .7 , .7 , .7 , .7 , .7 , . .7 , .7 , .7 , .6 , .5 , .5 , .5 , . .5 , .5 , .5 , .5 , .5 , .5 , .5 / DATA RATIOS / .7 , .7 , .7 , .7 , .7 , .7 , .7 , . .7 , .7 , .7 , .6 , .5 , .5 , .5 , . .5 , .5 , .5 , .5 , .5 , .5 , .5 / DATA TROPUV / 4.0 , 4.0 , 4.5 , 6.5 , 8.0 , 8.0 , 9.0 , . 10.0 , 13.0 , 13.0 , 11.0 , 8.0 , 8.0 , 8.0 , . 8.0 , 8.0 , 8.0 , 8.0 , 8.0 , 8.0 , 8.0 / DATA CHL / 120. , 100. , 80. , 80. , 80. , 80. , 80. , . 80. , 80. , 60. , 60. , 40. , 40. , 40. , . 40. , 40. , 40. , 40. , 40. , 40. , 40. / DATA FFLAT / 11*0 , .035 , .06 , .12 , .18 , .26 , .35 , . .45 , .54 , .63 , .72 , .80 , .86 , .92 , . .97 , 66*1 / DATA PMAND / 1000 , 850 , 700 , 500 , 400 , 300 , 250 , . 200 , 150 , 100 , 70 , 50 , 30 , 20 , . 10 , 7 , 5 , 3 , 2 , 1 , .4 / DATA P50 / 1000 , 950 , 900 , 850 , 800 , 750 , 700 , . 650 , 600 , 550 , 500 , 450 , 400 , 350 , . 300 , 250 , 200 , 150 , 100 , 50 , .4 / DATA CVC / 5.0 / END SUBROUTINE SETCON COMMON /CORP / CHLP(1500),CVC,FFLAT(91),IFA( 10 , 10 ) COMMON /FCRATS/ RATION(21),RATIOS(21),TROPUV(21) COMMON /OBERRS/ QIKE(21,300),PMAND(21),P50(21),CHL(21) COMMON /OBHEIR/ QUALOB(300) C-CRA COMMON /OBPRT / NPT,IDSTA(10),PSTA(10),KXSTA(10) COMMON /OBPRT / NPT,NPTADD,IDSTA(10),PSTA(10),KXSTA(10) CHARACTER*8 IDSTA c DATA RATION / .7 , .7 , .7 , .7 , .7 , .7 , .7 , c . .7 , .7 , .7 , .6 , .5 , .5 , .5 , c . .5 , .5 , .5 , .5 , .5 , .5 , .5 / c DATA RATIOS / .7 , .7 , .7 , .7 , .7 , .7 , .7 , c . .7 , .7 , .7 , .6 , .5 , .5 , .5 , c . .5 , .5 , .5 , .5 , .5 , .5 , .5 / c DATA TROPUV / 4.0 , 4.0 , 4.5 , 6.5 , 8.0 , 8.0 , 9.0 , c . 10.0 , 13.0 , 13.0 , 11.0 , 8.0 , 8.0 , 8.0 , c . 8.0 , 8.0 , 8.0 , 8.0 , 8.0 , 8.0 , 8.0 / c DATA CHL / 120. , 100. , 80. , 80. , 80. , 80. , 80. , c . 80. , 80. , 60. , 60. , 40. , 40. , 40. , c . 40. , 40. , 40. , 40. , 40. , 40. , 40. / c DATA FFLAT / 11*0 , .035 , .06 , .12 , .18 , .26 , .35 , c . .45 , .54 , .63 , .72 , .80 , .86 , .92 , c . .97 , 66*1 / c DATA PMAND / 1000 , 850 , 700 , 500 , 400 , 300 , 250 , c . 200 , 150 , 100 , 70 , 50 , 30 , 20 , c . 10 , 7 , 5 , 3 , 2 , 1 , .4 / c DATA P50 / 1000 , 950 , 900 , 850 , 800 , 750 , 700 , c . 650 , 600 , 550 , 500 , 450 , 400 , 350 , c . 300 , 250 , 200 , 150 , 100 , 50 , .4 / c c DATA CVC / 5.0 / C------------------------------------------------------------------- C------------------------------------------------------------------- C STORE A SYMMETRIC MATRIX MAP C ---------------------------- LLL = 1 DO 10 I = 1, 10 DO 10 J = 1,I IFA(I,J) = LLL IFA(J,I) = LLL LLL = LLL + 1 10 CONTINUE C MAKE A TABLE OF LENGTH SCALES BY MILIBARS C ----------------------------------------- DO 20 I=1,1500 CHLP(I) = PILNLNP(FLOAT(I),PMAND,CHL,21) 20 CONTINUE C READ IN SOME POINTS TO PRINT C ---------------------------- DO 30 I=1,10 READ(18,'(A6,I5,I4)',END=35) IDSTA(I),IP,KXSTA(I) PSTA(I) = IP PRINT'(A6,I5,I4)',IDSTA(I),IP,KXSTA(I) 30 CONTINUE 35 NPT = I-1 C SET UP SOME DATA PRIORITIES C --------------------------- DO 5 I=1,300 5 QUALOB(I) = 2. QUALOB(120) = 1. QUALOB(132) = 1. QUALOB(220) = 1. QUALOB(221) = 1. QUALOB(230) = 1. QUALOB(231) = 1. QUALOB(232) = 1. QUALOB(233) = 1. QUALOB(283) = 3. QUALOB(161) = 3. QUALOB(162) = 4. QUALOB(163) = 5. QUALOB(166) = 3. QUALOB(167) = 4. QUALOB(168) = 5. QUALOB(171) = 3. QUALOB(172) = 4. QUALOB(173) = 5. QUALOB(176) = 3. QUALOB(177) = 4. QUALOB(178) = 5. C SO LONG FOR NOW C --------------- RETURN END C$$$ SUBPROGRAM DOCUMENTATION BLOCK C C SUBPROGRAM: SSMIPTS C PRGMMR: WOOLLEN ORG: NMC22 DATE: 90-11-06 C C ABSTRACT: PUTS TOGETHER PACKAGES OF SSMI OBSERVATION LOCATIONS, C PERFORMS A MULTIVARIATE SURFACE WIND ANALYSIS, AND ASSIGNS C ANALYSED DIRECTION TO THE SSMI WIND SPEED OBSERVATION TO C PRODUCE AN SSMI WIND VECTOR TO BE CHECKED BY OIQC AND PASSED C ON THE OBJECTIVE ANALLYSIS PROGRAM. C C PROGRAM HISTORY LOG: C 90-11-06 J. WOOLLEN C C USAGE: CALL SSMIPTS C C OUTPUT ARGUMENTS: NONE C C SUBPROGRAMS CALLED: C UNIQUE: - SSMISRCH,OIQC,SSMIPUT C C LIBRARY: C COMMON - ABORT C C REMARKS: C C C ATTRIBUTES: C LANGUAGE: CRAY FORTRAN C MACHINE: CRAY C C$$$ SUBROUTINE SSMIPTS COMMON /CKLIST/ NDD,INDD( 300000 ) DIMENSION IND( 5000 ),NXY( 5000 ) DATA NMX / 64 / C---------------------------------------------------------------------- C---------------------------------------------------------------------- C ASSEMBLE LISTS OF OBS FOR MULTI PROCESSING C ------------------------------------------ IND(1) = 1 IPK = 1 C-CRA NXY = 0 C DIMENSION IND( 5000 ),NXY( 5000 ) DO I=1, 5000 NXY(I) = 0 ENDDO DO 5 I=1,NDD ILEV = INDD(I) IF(NXY(IPK)+1.GT.NMX) THEN IF(IPK+1.GT. 5000 ) CALL SABORT('SSMIPTS - IPK > MAXPAK') IPK = IPK+1 IND(IPK) = I NXY(IPK) = 1 ELSE NXY(IPK) = NXY(IPK)+1 ENDIF 5 CONTINUE C MULTI-PROCESS QC OB LISTS C ------------------------- CMIC$ DOALL AUTOSCOPE DO 10 I=1,IPK CALL SSMISRCH(IND(I),NXY(I)) CALL QCOI CMIC$ GUARD CALL SSMIPUT CMIC$ ENDGUARD 10 CONTINUE RETURN END C$$$ SUBPROGRAM DOCUMENTATION BLOCK C C SUBPROGRAM: SSMIPUT C PRGMMR: WOOLLEN ORG: NMC22 DATE: 92-10-06 C C ABSTRACT: USES THE OUTPUT OF A MULTIVARIATE INTERPOLATION OF C NEARBY OBSERVATIONS TO ASSIGN DIRECTIONS TO SSMI WIND SPEED DATA. C ALSO PROVIDES FOR A DETAILED PRINT SHOWING OBSERVATIONS USED IN C THE INTERPOLATION ALONG WITH THE WEIGHTS AND INCREMENTS RESULTING. C C PROGRAM HISTORY LOG: C 90-11-06 J. WOOLLEN C C USAGE: CALL SSMIPUT C C OUTPUT FILES: C FORT.71 - FORT.80 - DETAILED PRINTOUTS FOR ANALYSIS C C SUBPROGRAMS CALLED: C LIBRARY: C COMMON - ABORT C C REMARKS: C C C ATTRIBUTES: C LANGUAGE: CRAY FORTRAN C MACHINE: CRAY C C$$$ SUBROUTINE SSMIPUT C-CRA TASK COMMON/HVECT/NXXYY,MCHK( 64 ,2),MOBS( 64 , 10 , 4 ,2), COMMON/HVECT/NXXYY,NYYZZ,MCHK( 64 ,2),MOBS( 64 , 10 , 4 ,2) 1, . MDIM( 64 , 4 ),AA( 64 , 4 ,2), . AE( 64 , 4 ,2),WT( 64 , 10 , 4 ,2), . FE( 64 ,3) COMMON /OBS/ NREP,NLEV,SID( 200000 ),INOB( 200000 ),IMAP(360,181) C-CRA COMMON /OBPRT / NPT,SIDP(10),PRP(10),KXP(10) COMMON /OBPRT / NPT,NPTADD,SIDP(10),PRP(10),KXP(10) DIMENSION IOMF(3) CHARACTER*1 VAR2(2),VAR3(3) DATA VAR2 / 'Z' , 'W' / DATA VAR3 / 'Z' , 'U' , 'V' / DATA IOMF / 2 , 1 , 2 / CHARACTER*1 PLUMIN DATA PI180/.017453/ C---------------------------------------------------------------------- C---------------------------------------------------------------------- C LOOP OVER EACH OB CHECKED IN THIS GROUP C --------------------------------------- DO 50 N=1,NXXYY ILEV = MCHK(N,1) IREP = FHD(ILEV) PRES = FPR(ILEV) OBER = FOE(ILEV) CKSTA = SID(IREP) CLON = FLN(IREP) CLAT = FLT(IREP) KX = FKX(IREP) GESU = FMW(IREP,1) GESV = FMW(IREP,2) C COMPUTE THE ANALYSED WIND DIRECTION AT THIS OB C ---------------------------------------------- ANAU = GESU + AA(N,1,1) ANAV = GESV + AA(N,1,2) IF(ANAU.GT.0.) THEN ANADIR = ATAN(ANAV/ANAU) ELSE IF(ANAU.LT.0.) THEN ANADIR = ATAN(ANAV/ANAU) + 180.*PI180 ELSE IF(ANAU.EQ.0.) THEN ANADIR = SIGN(90.*PI180,ANAV) ENDIF C COMBINE THE ANA DIRECTION WITH THE SSMI SPEED C --------------------------------------------- SSMISPD = ABS(FOB(ILEV,2)) SSMIU = SSMISPD*COS(ANADIR) SSMIV = SSMISPD*SIN(ANADIR) C COMPUTE THE OBSERVED INCREMENT AND STORE IN THE PROPER PLACE C ------------------------------------------------------------ PACK = WOB(ILEV,1,SSMIU-GESU) PACK = WOB(ILEV,2,SSMIV-GESV) 50 CONTINUE 100 FORMAT(A6,I4,3F8.2,A4,2F8.2,4(' *',F6.1,' (',F3.2,')' ), . ' *',2F6.2,A6,A1) 101 FORMAT(A6,I4,3F8.2,A4,2F8.2,4(' *',F6.1,2X,I2),' *',2F6.2,A6,A1) 102 FORMAT(A6,I4,3F8.2,A4,6F8.2) RETURN END C----------------------------------------------------------------------- C C----------------------------------------------------------------------- SUBROUTINE SSMISRCH(IND,NXY) COMMON /CKLIST/ NDD,INDD( 300000 ) COMMON /CORP / CHLP(1500),CVC,FFLAT(91),IFA( 10 , 10 ) COMMON /OBS/ NREP,NLEV,SID( 200000 ),INOB( 200000 ),IMAP(360,181) C-CRA TASK COMMON/HVECT/NXXYY,MCHK( 64 ,2),MOBS( 64 , 10 , 4 ,2), COMMON/HVECT/NXXYY,NYYZZ,MCHK( 64 ,2),MOBS( 64 , 10 , 4 ,2) 1, . MDIM( 64 , 4 ),AA( 64 , 4 ,2), . AE( 64 , 4 ,2),WT( 64 , 10 , 4 ,2), . FE( 64 ,3) CHARACTER*3 TAG,FOC DIMENSION NIND( 200000 ),RLAT( 200000 ),RLON( 200000 ), . DIST( 200000 ),DIRN( 200000 ), . LIND( 50000 ),LTYP( 50000 ),PRES( 50000 ), . ICOR( 50000 ),IDIR( 50000 ), . IPIC( 10 ),CATCOR(101) LOGICAL CAT( 50000 ,4,4),NONE(4,4),MORE,GO DATA CATCOR /21*5.,30*4.,20*3.,20*2.,10*1./ DATA PI180 /.0174532/, RADE /6371./, RSCAN /1000./ DATA MAXCYL / 50000 / , MAXREP / 200000 / DATA MTOTH /10/ , MTOTV /5/ C---------------------------------------------------------------------- C---------------------------------------------------------------------- NXXYY = 0 C-CRA MDIM = 0 C . MDIM( 64 , 4 ),AA( 64 , 4 ,2), DO I=1, 64 * 4 MDIM(I,1) = 0 ENDDO C LOOP OVER PROFILES TO BE CHECKED C -------------------------------- DO 60 NPT=IND,IND+NXY-1 MPTR = INDD(NPT) NPTR = FHD(MPTR) BLON = FLN(NPTR) BLAT = FLT(NPTR) BINQ = FQL(NPTR)-1 TAG = FOC(NPTR) IF(TAG.NE.'SMI') CALL SABORT('SSMISRCH-SEARCHING NON SSMI DATA') C OUTLINE THE SEARCH AREA C ----------------------- DELLAT = RSCAN/(PI180*RADE) COSLAT = COS(BLAT*PI180) DELLON = MIN(DELLAT/COSLAT,180.0) X1 = BLON-DELLON X2 = BLON+DELLON Y1 = BLAT-DELLAT Y2 = BLAT+DELLAT IF(X1.LT. 0.) X1 = X1 + 360. IF(X2.GE.360.) X2 = X2 - 360. IF(Y1.LT.-90.) Y1 = -90. IF(Y2.GT. 90.) Y2 = 90. C ITEMIZE REPORTS IN THE SEACH CYLINDER C ------------------------------------- NPRF = 0 NOBS = 0 I1 = X1 + 1. I2 = X2 + 1. J1 = Y1 + 91. J2 = Y2 + 91. DO 10 K=1,2 IA = I1 IB = I2 IF(I1.GT.I2) THEN IF(K.EQ.1) IB = 360 IF(K.EQ.2) IA = 1 ELSE IF(K.EQ.2) GOTO 11 ENDIF DO 10 J=J1,J2 DO 10 N=IMAP(IA,J),IMAP(IB,J) NN = INOB(N) QLA = FQL(NN) IF(QLA.GE.0. .AND. QLA.LE.BINQ) THEN NPRF = NPRF+1 NIND(NPRF) = NN RLON(NPRF) = FLN(NN) RLAT(NPRF) = FLT(NN) ENDIF 10 CONTINUE 11 IF(NPRF.GT.MAXREP) CALL SABORT('SSMISRCH - NPRF > MAXREP') C COMPUTE DISTANCE AND DIRECTION FROM THE CYLINDER CENTER C ------------------------------------------------------- CALL CHDIST(BLON,BLAT,RLON,RLAT,DIST,DIRN,NPRF) C STORE ALL LEVELS OF EACH REPORT LOCATED IN THE CYLINDER C ------------------------------------------------------- DO 21 N=1,NPRF IF(DIST(N).LE.RSCAN .AND. SID(NPTR).NE.SID(NIND(N))) THEN ILEV = FLV(NIND(N)) KLEV = ILEV+FNL(NIND(N))-1 KX = FKX(NIND(N)) DO 20 JLEV=ILEV,KLEV IF(FOI(JLEV).GT.0) THEN NOBS = NOBS+1 LIND(NOBS) = JLEV PRES(NOBS) = LOG(FPR(JLEV)) DIST(NOBS) = (DIST(N)/RADE)**2 DIRN(NOBS) = AMOD(DIRN(N)/90.,4.) + 1.01 LTYP(NOBS) = KX/100 ENDIF 20 CONTINUE IF(NOBS.GT.MAXCYL) CALL SABORT('SSMISRCH - NOBS>MAXCYL') ENDIF 21 CONTINUE C LOOP OVER THE OBSERVATION LEVELS TO BE CHECKED C ---------------------------------------------- ILEV = FLV(NPTR) KLEV = ILEV+FNL(NPTR)-1 DO 50 IL=ILEV,KLEV BPRS = FPR(IL) IPRS = BPRS BPRS = LOG(BPRS) C INITIALIZE THE SLECTION ARRAYS C ------------------------------ NXXYY = NXXYY + 1 MCHK(NXXYY,1) = IL MCHK(NXXYY,2) = 2 FE (NXXYY,1) = FFE(IL,1) FE (NXXYY,2) = FFE(IL,2) FE (NXXYY,3) = FFE(IL,2) IF(NOBS.EQ.0) GOTO 50 C COMPUTE ZZ HRZ AND VRT CORRELATIONS FOR EACH OB AND CATEGORIZE C -------------------------------------------------------------- DO 24 N=1,NOBS CHRZ = EXP(-CHLP(IPRS)*(DIST(N))) CVRT = 1./(1.+ CVC*(BPRS-PRES(N))**2) IC = CHRZ*CVRT*100.+1. ICOR(N) = CATCOR(IC) IDIR(N) = DIRN(N) 24 CONTINUE DO 25 J=1,4 DO 25 I=1,4 DO 25 N=1,NOBS CAT(N,I,J) = ICOR(N).EQ.J .AND. IDIR(N).EQ.I 25 CONTINUE C PICK THE CLOSEST 10 (MULTIVARIATE) HORIZONTAL OBSERVATIONS C ---------------------------------------------------------- DO 28 J=1,4 DO 28 I=1,4 NONE(I,J) = .FALSE. 28 CONTINUE NTOT = 0 NPIC = 0 DO 35 J=1,4 30 MORE = .FALSE. DO 34 I=1,4 IF(NPIC.GE.MTOTH) GOTO 36 IF(NONE(I,J) ) GOTO 34 DO 32 N=1,NOBS IF(CAT(N,I,J)) THEN NTOT = NTOT+1 IPIC(NTOT) = N NPIC = NPIC+LTYP(N) CAT(N,I,J) = .FALSE. MORE = .TRUE. GOTO 34 ENDIF 32 CONTINUE NONE(I,J) = .TRUE. 34 CONTINUE IF(MORE) GOTO 30 35 CONTINUE 36 IF(NPIC.GT.MTOTH) NTOT = NTOT-1 C STORE THE DATA INDEXES SELECTED FOR HORIZONTAL CHECKS C ----------------------------------------------------- NPIC = 0 DO 40 N=1,NTOT NN = IPIC(N) KK = LTYP(NN) DO 40 K=KK,KK*2-1 NPIC = NPIC+1 MDIM(NXXYY,1) = MDIM(NXXYY,1) + 1 MOBS(NXXYY,NPIC,1,1) = LIND(NN) MOBS(NXXYY,NPIC,1,2) = K 40 CONTINUE 50 CONTINUE 60 CONTINUE RETURN END C$$$ SUBPROGRAM DOCUMENTATION BLOCK C C SUBPROGRAM: STOERR C PRGMMR: WOOLLEN ORG: NMC22 DATE: 90-11-06 C C ABSTRACT: INTERPOLATES TO AND STORES FORECAST AND OBSERVATION ERRORS C WITH EACH OBSERVATION IN THE DATASET IN COMMON /DUMMY/. ALSO C MAKES A LIST OF OBS BY 2D LOCATION IN COMMON /TGRID/ IN ORDER C TO EXPEDITE VARIOUS DATA SEARCH PROCEDURES. SOME CHECKING ON THE C PLAUSIBILITY OF OB COORDINATES IS PERFORMED AS WELL. C C PROGRAM HISTORY LOG: C 90-11-06 J. WOOLLEN C C USAGE: C INPUT ARGUMENTS: NONE C OUTPUT ARGUMENTS IN COMMON /TGRID/ C NREP - NUMBER OF REPORT LOCATIONS C XQC - ARRAY OF REPORT LONGITUDES C YQC - ARRAY OF REPORT LATITUDES C IQC - ARRAY OF REPORT INDEXES INTO THE BUFFS ARRAY C TQC - ARRAY OF REPORT TYPES (HEIGHT OR WIND) C FERR - ARRAY OF FORECAST ERRORS AT REPORT LOCATIONS C C SUBPROGRAMS CALLED: C LIBRARY: C COMMON - ABORT C C EXIT STATES: C COND = 11 - PROBLEM DETECTED (SEE CONSOLE OUTPUT) C C REMARKS: C C C ATTRIBUTES: C LANGUAGE: CRAY FORTRAN C MACHINE: CRAY C C$$$ C----------------------------------------------------------------------- SUBROUTINE STOERR COMMON /OBS/ NREP,NLEV,SID( 200000 ),INOB( 200000 ),IMAP(360,181) COMMON /OBERRS/ QIKE(21,300),PMAND(21),P50(21),CHL(21) COMMON /CORP / CHLP(1500),CVC,FFLAT(91),IFA( 10 , 10 ) COMMON /GETOE / USEOE CHARACTER*3 FOC,CLASS LOGICAL USEOE DATA XMSG/32700./ C---------------------------------------------------------------------- C----------------------------------------------------------------------- C GO THROUGH ALL THE STORED DATA C ------------------------------ DO 30 IREP=1,NREP ILEV = FLV(IREP) CLASS = FOC(IREP) STID = SID(IREP) DDLAT = FLT(IREP) DDLON = FLN(IREP) KLEV = FNL(IREP) KX = FKX(IREP) KKX = KX/100 IF(KLEV.EQ.0) GOTO 30 C CHECK FOR UNNATURAL INFORMATION IN THE HEADER C --------------------------------------------- IF( . ( KX.LT.120 .OR. KX.GT.299 ) .OR. . ( KLEV.LE.0 ) .OR. . ( DDLAT.LT.-90. .OR. DDLAT.GT.90. ) .OR. . ( DDLON.LT.-180. .OR. DDLON.GT.360.) . ) THEN PRINT*,'FROM STOERR---SKIP THIS MESS---FROM STOERR' PRINT'(''STATION= '',A8)',STID PRINT*,'LAT = ',DDLAT PRINT*,'LON = ',DDLON PRINT*,'KX = ',KX PRINT*,'NLEVS = ',KLEV PRINT*,'FROM STOERR---SKIP THIS MESS---FROM STOERR' GOTO 30 ENDIF C LOOP THRU THE REPORT LEVELS C --------------------------- DO 20 KL=1,KLEV L = ILEV+KL-1 C CHECK TO SEE IF A SURFACE OBSERVATION IS PRESENT C ------------------------------------------------ IF(FOE(L).EQ.0.) THEN GOTO 20 ENDIF C GET FORECAST ERRORS FOR THIS REPORT C ----------------------------------- DELP = FPR(L)*.053 DDPRS = FPR(L) DDPRA = DDPRS - DELP DDPRB = DDPRS + DELP ZERA = GETZER(DDLAT,DDLON,DDPRA) ZERB = GETZER(DDLAT,DDLON,DDPRB) WERR = GETWER(DDLAT,DDLON,DDPRS) C FIND THE OB ERROR FOR THIS REPORT C --------------------------------- IF(KKX.EQ.1) THEN OBERA = PILNLNP(DDPRA,PMAND,QIKE(1,KX),21) OBERB = PILNLNP(DDPRB,PMAND,QIKE(1,KX),21) FPE = 1./(1.+CVC*LOG(DDPRA/DDPRB)**2) OBER = SQRT(OBERA**2 + OBERB**2 - 2*OBERA*OBERB*FPE) ELSE IF(KKX.EQ.2) THEN OBER = PILNLNP(DDPRS,P50,QIKE(1,KX),21) ENDIF C STORE THE ERRORS IN THEIR RIGHTFUL PLACE C ---------------------------------------- IF(ZERA.EQ.0 .OR. ZERB.EQ.0. .OR. WERR.EQ.0.) THEN PRINT*,'STOERR - FORECAST ERROR IS ZERO' CALL SABORT('STOERR - FORECAST ERROR IS ZERO') ENDIF FPE = 1./(1.+CVC*LOG(DDPRA/DDPRB)**2) TER = SQRT(ZERA**2+ZERB**2-2*ZERA*ZERB*FPE) WER = WERR FPE1 = 1./(1.+CVC*LOG(DDPRA/DDPRS)**2) FPE2 = 1./(1.+CVC*LOG(DDPRB/DDPRS)**2) TZF = MAX((ZERA*FPE1 - ZERB*FPE2),0.) / TER PACK = WFE(L,1,TER) PACK = WFE(L,2,WER) PACK = WTZ(L ,TZF) C SPECIAL TREATMENT FOR SPECIAL CASES C ----------------------------------- IF(CLASS.EQ.'SFC' .AND. KL.EQ.2) THEN PACK = WFE(L,1,(ZERA+ZERB)/2.) OBER = (OBERA+OBERB)/2. PACK = WTZ(L,1.) ENDIF IF(CLASS.NE.'SAT' .AND. DDPRS.GT.700.) THEN UPERR = 1.5*DDPRS/700. PACK = WFE(L,1,FFE(L,1)*UPERR) ENDIF C STORE AN OIQC OBSERVATION ERROR MAYBE C ------------------------------------- IF(USEOE) PACK = WOE(L,OBER) C END OF THE LOOPS C ---------------- 20 CONTINUE 30 CONTINUE PRINT*,' IN STORERR - NREP=',NREP RETURN END C$$$ SUBPROGRAM DOCUMENTATION BLOCK C C SUBPROGRAM: TOSSLIST C PRGMMR: WOOLLEN ORG: NMC22 DATE: 90-11-06 C C ABSTRACT: WRITES A LIST OF REJECTED OBSERVATIONS PLUS A SUMMARY C C PROGRAM HISTORY LOG: C 90-11-06 J. WOOLLEN C C USAGE: C INPUT ARGUMENTS: NONE C LO - LOGICAL UNIT FOR TOSSLIST OUTPUT C C OUTPUT ARGUMENTS: NONE C C OUTPUT FILES: C FORT.LO - TOSSLIST - OFFICIAL OPERATIONAL TOSSLIST C C SUBPROGRAMS CALLED: C LIBRARY: C COMMON - ABORT C C REMARKS: C C C ATTRIBUTES: C LANGUAGE: CRAY FORTRAN C MACHINE: CRAY C C$$$ C----------------------------------------------------------------------- SUBROUTINE TOSSLIST(LO) COMMON /OBS/ NREP,NLEV,SID( 200000 ),INOB( 200000 ),IMAP(360,181) CHARACTER*3 TAG,FOC DIMENSION IBAD(30,3),ITOT(30),ISFZ(30),PZ(30),PU(30),PT(30) LOGICAL SFZ,SAT DATA XMSG/32767./ C---------------------------------------------------------------------- C----------------------------------------------------------------------- C START UP C -------- CRA IBAD = 0 CRA ITOT = 0 CRA ISFZ = 0 CRA PZ = 0 CRA PU = 0 CRA PT = 0 C DIMENSION IBAD(30,3),ITOT(30),ISFZ(30),PZ(30),PU(30),PT(30) DO I=1,30*3 IBAD(I,1) = 0 ENDDO DO I=1,30 ITOT(I) = 0 ISFZ(I) = 0 PZ (I) = 0 PU (I) = 0 PT (I) = 0 ENDDO C COUNT BAD DATA IN THE BIG BUFFER - WRITE THE TOSSLIST C ----------------------------------------------------- N = 0 DO 10 IREP=1,NREP X = FLN(IREP) Y = FLT(IREP) STAID = SID(IREP) KX = FKX(IREP) TAG = FOC(IREP) ILEV = FLV(IREP) KLEV = ILEV+FNL(IREP)-1 DO 10 JLEV=ILEV,KLEV SFZ = TAG.EQ.'SFC' .AND. JLEV.EQ.KLEV SAT = TAG.EQ.'SAT' IF(FOI(JLEV).NE.0) THEN IF(SFZ) THEN ISFZ(KX/10) = ISFZ(KX/10) + 1 ELSE ITOT(KX/10) = ITOT(KX/10) + 1 ENDIF ENDIF IF(FOI(JLEV).LT.0) THEN IP = FPR(JLEV) IQ = FQM(JLEV) OBER = FOE(JLEV) BUFF11 = FOB(JLEV,1) BUFF12 = FOB(JLEV,2) IF(FOI(JLEV).EQ.-1) THEN IF(SFZ) THEN N = N+1 WRITE(LO,2002) N,STAID,X,Y,IP,KX,IQ,OBER,BUFF11,BUFF12 IBAD(KX/10,3) = IBAD(KX/10,3) + 1 ELSE IF(.NOT.SAT) THEN N = N+1 WRITE(LO,2000) N,STAID,X,Y,IP,KX,IQ,OBER,BUFF11,BUFF12 ENDIF IBAD(KX/10,1) = IBAD(KX/10,1) + 1 ENDIF ELSE IF(FOI(JLEV).EQ.-2) THEN N = N+1 WRITE(LO,2001) N,STAID,X,Y,IP,KX,IQ,OBER,BUFF11,BUFF12 IBAD(KX/10,2) = IBAD(KX/10,2) + 1 ENDIF ENDIF 10 CONTINUE C COMPUTE GROSS PERCENTAGE OF TOSSES C ---------------------------------- DO 20 I=10,29 IF(ISFZ(I).GT.0) THEN PZ(I) = FLOAT(IBAD(I,3))/FLOAT(ISFZ(I)) ENDIF IF(ITOT(I).GT.0) THEN PU(I) = FLOAT(IBAD(I,2))/FLOAT(ITOT(I)) PT(I) = FLOAT(IBAD(I,1))/FLOAT(ITOT(I)) ENDIF 20 CONTINUE C WRITE A SHORT SUMMARY OF TOSSES C ------------------------------- WRITE( 6,2003)(I,I=100,290,10),(IBAD(I,3),I=10,29),(PZ(I),I=10,29) WRITE( 6,2004)(I,I=100,290,10),(IBAD(I,2),I=10,29),(PU(I),I=10,29) WRITE( 6,2005)(I,I=100,290,10),(IBAD(I,1),I=10,29),(PT(I),I=10,29) WRITE(60,2003)(I,I=100,290,10),(IBAD(I,3),I=10,29),(PZ(I),I=10,29) WRITE(60,2004)(I,I=100,290,10),(IBAD(I,2),I=10,29),(PU(I),I=10,29) WRITE(60,2005)(I,I=100,290,10),(IBAD(I,1),I=10,29),(PT(I),I=10,29) C FORMAT STATEMENTS FOR THE TOSS LIST C ----------------------------------- 2000 FORMAT(' ',I5,1X,'J',3X,A8,2X,4X,2F10.2,I9,2X,2I5,F8.1,F13.1, . 5X,F13.1,10X,F6.1) 2001 FORMAT(' ',I5,1X,'P',3X,A8,2X,4X,2F10.2,I9,2X,2I5,F8.1,F13.1, . 5X,F13.1,10X,F6.1) 2002 FORMAT(' ',I5,1X,'Z',3X,A8,2X,4X,2F10.2,I9,2X,2I5,F8.1,F13.1, . 5X,F13.1,10X,F6.1) 2003 FORMAT(/' SURFACE Z TOSSES '/20I5/20I5/20F5.2) 2004 FORMAT(/' BELOW GROUND TOSSES'/20I5/20I5/20F5.2) 2005 FORMAT(/' QC OI TOSSES '/20I5/20I5/20F5.2) RETURN END CCBREAKCC C$$$ SUBPROGRAM DOCUMENTATION BLOCK C C SUBPROGRAM: VSOLVE C PRGMMR: WOOLLEN ORG: NMC22 DATE: 90-11-06 C C ABSTRACT: CHOLESKY TYPE SOLUTION FOR ARRAYS OF POSITIVE DEFINITE C SYMMETRIC MATRIXES. C C PROGRAM HISTORY LOG: C 90-11-06 J. WOOLLEN C C USAGE: C INPUT ARGUMENTS: C A - ARRAY OF SYMMETRIC MATRIXES C B - ARRAY OF RIGHT HAND SIDE VECTORS C NDIM - ARRAY OF MATRIX RANKS C BAD - BAD MATRIX INDICATOR C NFT - NUMBER OF RIGHT HAND SIDES PER MATRIX C NS - NUMBER OF MATRIXES TO SOLVE C MAXDIM - LARGEST RANK MATRIX IN STACK C C OUTPUT ARGUMENTS: C B - CONTAINS THE SOLUTIONS UPON RETURN C C SUBPROGRAMS CALLED: NONE C C EXIT STATES: NONE C C REMARKS: TRIANGULAR REPRESENTATION LISTS ELEMENTS TOWARDS THE C DIAGONAL. ALGORITHM FROM H.CARUS. C C ATTRIBUTES: C LANGUAGE: CRAY FORTRAN C MACHINE: CRAY C C$$$ C----------------------------------------------------------------------- SUBROUTINE VSOLVE (A,B,NDIM,BAD,NS,MAXDIM) PARAMETER( NMX = 64 , . NN = 10 *( 10 +1)/2 , . N = 10 ) DIMENSION A(NMX,NN),B(NMX,N,2),NDIM(NMX),BAD(NMX) LOGICAL BAD DIMENSION T(NMX) DATA CNUM/1.E-15/, NFT/2/ C---------------------------------------------------------------------- IX (I,J) = I*(I-1)/2 + J C---------------------------------------------------------------------- DO 1 M=1,NS 1 BAD(M) = .FALSE. C DECOMPOSE THE MATRIXES C ---------------------- DO 10 I=1,MAXDIM DO 10 J=1,I DO 2 M=1,NS 2 T(M) = A(M,IX(I,J)) DO 3 K=1,J-1 DO 3 M=1,NS 3 T(M) = T(M) - A(M,IX(I,K)) * A(M,IX(J,K)) IF(I.GT.J) THEN DO 4 M=1,NS 4 A(M,IX(I,J)) = T(M) * A(M,IX(J,J)) ELSE DO 5 M=1,NS IF(T(M).LT.CNUM .AND. NDIM(M).GE.I) BAD(M) = .TRUE. IF(T(M).LE.0) T(M) = 1. 5 CONTINUE DO 6 M=1,NS 6 A(M,IX(I,I)) = 1./SQRT(T(M)) ENDIF 10 CONTINUE C SOLVE FOR ALL RIGHT HAND SIDES C ------------------------------ DO 40 NF=1,NFT C FORWARD SUBSTITUTION C -------------------- DO 20 I=1,MAXDIM DO 11 M=1,NS 11 T(M) = B(M,I,NF) DO 12 J=1,I-1 DO 12 M=1,NS 12 T(M) = T(M) - A(M,IX(I,J)) * B(M,J,NF) DO 20 M=1,NS 20 B(M,I,NF) = T(M) * A(M,IX(I,I)) C BACKWARD SUBSTITUTION C --------------------- DO 30 I=MAXDIM,1,-1 DO 21 M=1,NS 21 T(M) = B(M,I,NF) IF(I.NE.MAXDIM) THEN DO 22 J=I+1,MAXDIM DO 22 M=1,NS 22 T(M) = T(M) - A(M,IX(J,I)) * B(M,J,NF) ENDIF DO 30 M=1,NS 30 B(M,I,NF) = T(M) * A(M,IX(I,I)) 40 CONTINUE RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C .=================================================================== C | DESCRIPTION | UNPK | PACK | BITS | * | DIMENSIONS C |=============*======*======*======*===*============================ C | LEVEL LINK | FLV | WLV | WORD | I | ?MAXREP C |-------------+------+------+------+---+---------------------------- C | OB CLASS | FOC | WOC | C3 | C | ?MAXREP C |-------------+------+------+------+---+---------------------------- C | LATITUDE | FLT | WLT | WORD | R | ?MAXREP C |-------------+------+------+------+---+---------------------------- C | LONGITUDE | FLN | WLN | WORD | R | ?MAXREP C |-------------+------+------+------+---+---------------------------- C | OB TIME | FTM | WTM | WORD | R | ?MAXREP C |-------------+------+------+------+---+---------------------------- C | OB TYPE | FKX | WKX | WORD | I | ?MAXREP C |-------------+------+------+------+---+---------------------------- C | N OF LEVELS | FNL | WNL | WORD | I | ?MAXREP C |-------------+------+------+------+---+---------------------------- C | USE REP FLG | FQL | WQL | WORD | R | ?MAXREP C |-------------+------+------+------+---+---------------------------- C | PREPDA REC | FPRC | WPRC | WORD | I | ?MAXREP C |-------------+------+------+------+---+---------------------------- C | PREPDA FLD | FPFD | WPFD | WORD | I | ?MAXREP C |-------------+------+------+------+---+---------------------------- C | PREPDA LEV | FPLV | WPLV | WORD | I | ?MAXLEV C |-------------+------+------+------+---+---------------------------- C | GES SFC U,V | FMW | WMW | WORD | R | ?MAXREP,2 C |-------------+------+------+------+---+---------------------------- C | HEADER LINK | FHD | WHD | WORD | I | ?MAXLEV C |-------------+------+------+------+---+---------------------------- C | PRESSURE | FPR | WPR | WORD | R | ?MAXLEV C |-------------+------+------+------+---+---------------------------- C | PREDA FLAG | FQM | WQM | WORD | I | ?MAXLEV C |-------------+------+------+------+---+---------------------------- C | OIQC FLAG | FOI | WOI | WORD | I | ?MAXLEV C |-------------+------+------+------+---+---------------------------- C | OB INCS | FOB | WOB | WORD | R | ?MAXLEV,2 C |-------------+------+------+------+---+---------------------------- C | OB ERROR | FOE | WOE | WORD | R | ?MAXLEV C |-------------+------+------+------+---+---------------------------- C | T/W FC ERR | FFE | WFE | WORD | R | ?MAXLEV,2 C |-------------+------+------+------+---+---------------------------- C | T/Z FACTOR | FTZ | WTZ | WORD | R | ?MAXLEV C |-------------+------+------+------+---+---------------------------- C | T/W OBS | FTW | WTW | WORD | I | ?MAXLEV,10,2 C |-------------+------+------+------+---+---------------------------- C | PROF A/B | FPF | WPF | WORD | I | ?MAXLEV,5,2 C `=============^======^======^======^===^============================ C C----------------------------------------------------------------------- C----------------------------------------------------------------------- FUNCTION FLV(I) PARAMETER (IM = 200000 ) COMMON /STOR1/IAR(IM) INTEGER IAR FLV = IAR(I) RETURN END C----------------------------------------------------------------------- FUNCTION WLV(I,V) PARAMETER (IM = 200000 ) COMMON /STOR1/IAR(IM) INTEGER IAR,V IAR(I) = V WLV = 0 RETURN END C----------------------------------------------------------------------- FUNCTION FOC(I) PARAMETER (IM = 200000 ) COMMON /STOR2/IAR(IM) CHARACTER*3 IAR,FOC FOC = IAR(I) RETURN END C----------------------------------------------------------------------- FUNCTION WOC(I,V) PARAMETER (IM = 200000 ) COMMON /STOR2/IAR(IM) CHARACTER*3 IAR,V IAR(I) = V WOC = 0 RETURN END C----------------------------------------------------------------------- FUNCTION FLT(I) PARAMETER (IM = 200000 ) COMMON /STOR3/IAR(IM) REAL IAR FLT = IAR(I) RETURN END C----------------------------------------------------------------------- FUNCTION WLT(I,V) PARAMETER (IM = 200000 ) COMMON /STOR3/IAR(IM) REAL IAR,V IAR(I) = V WLT = 0 RETURN END C----------------------------------------------------------------------- FUNCTION FLN(I) PARAMETER (IM = 200000 ) COMMON /STOR4/IAR(IM) REAL IAR FLN = IAR(I) RETURN END C----------------------------------------------------------------------- FUNCTION WLN(I,V) PARAMETER (IM = 200000 ) COMMON /STOR4/IAR(IM) REAL IAR,V IAR(I) = V WLN = 0 RETURN END C----------------------------------------------------------------------- FUNCTION FTM(I) PARAMETER (IM = 200000 ) COMMON /STOR5/IAR(IM) REAL IAR FTM = IAR(I) RETURN END C----------------------------------------------------------------------- FUNCTION WTM(I,V) PARAMETER (IM = 200000 ) COMMON /STOR5/IAR(IM) REAL IAR,V IAR(I) = V WTM = 0 RETURN END C----------------------------------------------------------------------- FUNCTION FKX(I) PARAMETER (IM = 200000 ) COMMON /STOR6/IAR(IM) INTEGER IAR FKX = IAR(I) RETURN END C----------------------------------------------------------------------- FUNCTION WKX(I,V) PARAMETER (IM = 200000 ) COMMON /STOR6/IAR(IM) INTEGER IAR,V IAR(I) = V WKX = 0 RETURN END C----------------------------------------------------------------------- FUNCTION FNL(I) PARAMETER (IM = 200000 ) COMMON /STOR7/IAR(IM) INTEGER IAR FNL = IAR(I) RETURN END C----------------------------------------------------------------------- FUNCTION WNL(I,V) PARAMETER (IM = 200000 ) COMMON /STOR7/IAR(IM) INTEGER IAR,V IAR(I) = V WNL = 0 RETURN END C----------------------------------------------------------------------- FUNCTION FQL(I) PARAMETER (IM = 200000 ) COMMON /STOR8/IAR(IM) REAL IAR FQL = IAR(I) RETURN END C----------------------------------------------------------------------- FUNCTION WQL(I,V) PARAMETER (IM = 200000 ) COMMON /STOR8/IAR(IM) REAL IAR,V IAR(I) = V WQL = 0 RETURN END C----------------------------------------------------------------------- FUNCTION FPRC(I) PARAMETER (IM = 200000 ) COMMON /STOR9/IAR(IM) INTEGER IAR FPRC = IAR(I) RETURN END C----------------------------------------------------------------------- FUNCTION WPRC(I,V) PARAMETER (IM = 200000 ) COMMON /STOR9/IAR(IM) INTEGER IAR,V IAR(I) = V WPRC = 0 RETURN END C----------------------------------------------------------------------- FUNCTION FPFD(I) PARAMETER (IM = 200000 ) COMMON /STOR10/IAR(IM) INTEGER IAR FPFD = IAR(I) RETURN END C----------------------------------------------------------------------- FUNCTION WPFD(I,V) PARAMETER (IM = 200000 ) COMMON /STOR10/IAR(IM) INTEGER IAR,V IAR(I) = V WPFD = 0 RETURN END C----------------------------------------------------------------------- FUNCTION FPLV(I) PARAMETER (IM = 300000 ) COMMON /STOR11/IAR(IM) INTEGER IAR FPLV = IAR(I) RETURN END C----------------------------------------------------------------------- FUNCTION WPLV(I,V) PARAMETER (IM = 300000 ) COMMON /STOR11/IAR(IM) INTEGER IAR,V IAR(I) = V WPLV = 0 RETURN END C----------------------------------------------------------------------- FUNCTION FMW(I,J) PARAMETER (IM = 200000 ) PARAMETER (JM = 2) COMMON /STOR12/IAR(IM,JM) REAL IAR FMW = IAR(I,J) RETURN END C----------------------------------------------------------------------- FUNCTION WMW(I,J,V) PARAMETER (IM = 200000 ) PARAMETER (JM = 2) COMMON /STOR12/IAR(IM,JM) REAL IAR,V IAR(I,J) = V WMW = 0 RETURN END C----------------------------------------------------------------------- FUNCTION FHD(I) PARAMETER (IM = 300000 ) COMMON /STOR13/IAR(IM) INTEGER IAR FHD = IAR(I) RETURN END C----------------------------------------------------------------------- FUNCTION WHD(I,V) PARAMETER (IM = 300000 ) COMMON /STOR13/IAR(IM) INTEGER IAR,V IAR(I) = V WHD = 0 RETURN END C----------------------------------------------------------------------- FUNCTION FPR(I) PARAMETER (IM = 300000 ) COMMON /STOR14/IAR(IM) REAL IAR FPR = IAR(I) RETURN END C----------------------------------------------------------------------- FUNCTION WPR(I,V) PARAMETER (IM = 300000 ) COMMON /STOR14/IAR(IM) REAL IAR,V IAR(I) = V WPR = 0 RETURN END C----------------------------------------------------------------------- FUNCTION FQM(I) PARAMETER (IM = 300000 ) COMMON /STOR15/IAR(IM) INTEGER IAR FQM = IAR(I) RETURN END C----------------------------------------------------------------------- FUNCTION WQM(I,V) PARAMETER (IM = 300000 ) COMMON /STOR15/IAR(IM) INTEGER IAR,V IAR(I) = V WQM = 0 RETURN END C----------------------------------------------------------------------- FUNCTION FOI(I) PARAMETER (IM = 300000 ) COMMON /STOR16/IAR(IM) INTEGER IAR FOI = IAR(I) RETURN END C----------------------------------------------------------------------- FUNCTION WOI(I,V) PARAMETER (IM = 300000 ) COMMON /STOR16/IAR(IM) INTEGER IAR,V IAR(I) = V WOI = 0 RETURN END C----------------------------------------------------------------------- FUNCTION FOB(I,J) PARAMETER (IM = 300000 ) PARAMETER (JM = 2) COMMON /STOR17/IAR(IM,JM) REAL IAR FOB = IAR(I,J) RETURN END C----------------------------------------------------------------------- FUNCTION WOB(I,J,V) PARAMETER (IM = 300000 ) PARAMETER (JM = 2) COMMON /STOR17/IAR(IM,JM) REAL IAR,V IAR(I,J) = V WOB = 0 RETURN END C----------------------------------------------------------------------- FUNCTION FOE(I) PARAMETER (IM = 300000 ) COMMON /STOR18/IAR(IM) REAL IAR FOE = IAR(I) RETURN END C----------------------------------------------------------------------- FUNCTION WOE(I,V) PARAMETER (IM = 300000 ) COMMON /STOR18/IAR(IM) REAL IAR,V IAR(I) = V WOE = 0 RETURN END C----------------------------------------------------------------------- FUNCTION FFE(I,J) PARAMETER (IM = 300000 ) PARAMETER (JM = 2) COMMON /STOR19/IAR(IM,JM) REAL IAR FFE = IAR(I,J) RETURN END C----------------------------------------------------------------------- FUNCTION WFE(I,J,V) PARAMETER (IM = 300000 ) PARAMETER (JM = 2) COMMON /STOR19/IAR(IM,JM) REAL IAR,V IAR(I,J) = V WFE = 0 RETURN END C----------------------------------------------------------------------- FUNCTION FTZ(I) PARAMETER (IM = 300000 ) COMMON /STOR20/IAR(IM) REAL IAR FTZ = IAR(I) RETURN END C----------------------------------------------------------------------- FUNCTION WTZ(I,V) PARAMETER (IM = 300000 ) COMMON /STOR20/IAR(IM) REAL IAR,V IAR(I) = V WTZ = 0 RETURN END C----------------------------------------------------------------------- FUNCTION FTW(I,J,K) PARAMETER (IM = 300000 ) PARAMETER (JM = 10) PARAMETER (KM = 2) COMMON /STOR21/IAR(IM,JM,KM) INTEGER IAR FTW = IAR(I,J,K) RETURN END C----------------------------------------------------------------------- FUNCTION WTW(I,J,K,V) PARAMETER (IM = 300000 ) PARAMETER (JM = 10) PARAMETER (KM = 2) COMMON /STOR21/IAR(IM,JM,KM) INTEGER IAR,V IAR(I,J,K) = V WTW = 0 RETURN END C----------------------------------------------------------------------- FUNCTION FPF(I,J,K) PARAMETER (IM = 300000 ) PARAMETER (JM = 5) PARAMETER (KM = 2) COMMON /STOR22/IAR(IM,JM,KM) INTEGER IAR FPF = IAR(I,J,K) RETURN END C----------------------------------------------------------------------- FUNCTION WPF(I,J,K,V) PARAMETER (IM = 300000 ) PARAMETER (JM = 5) PARAMETER (KM = 2) COMMON /STOR22/IAR(IM,JM,KM) INTEGER IAR,V IAR(I,J,K) = V WPF = 0 RETURN END C$$$ SUBPROGRAM DOCUMENTATION BLOCK C . . . . C SUBPROGRAM: STORE STORES OBSERVED DATA IN BUFFER C PRGMMR: DEAVEN ORG: W/NMC22 DATE: 88-02-10 C C ABSTRACT: STORES ALL DATA AND CORRESPONDING INDEX FOR LATER USE C IN DATA COLLECT AND SEARCH ROUTINES. C C PROGRAM HISTORY LOG: C 86-03-21 G. DIMEGO C 88-11-24 D. DEAVEN CODED FOR CYBER 205 C 90-11-06 J. WOOLLEN MODIFIED COMMON /BUFFS/ FOR QCOI C 92-07-28 J. WOOLLEN MODIFIED FOR STORSTAR VERSION C C USAGE: CALL STORE(LU) C C INPUT ARGUMENTS: C LU - LOGICAL UNIT FOR INCREMENT FILE C C ATTRIBUTES: C LANGUAGE: CDC FORTRAN200 C MACHINE: CYBER C C$$$ SUBROUTINE STORE(LUDAT,LUBFI) COMMON /OBS/ NREP,NLEV,SID( 200000 ),INOB( 200000 ),IMAP(360,181) COMMON /OBHEIR/ QUALOB(300) CHARACTER*80 SUCCESS(3),WARNING(3) CHARACTER*80 HDSTR,OBSTR,FCSTR,OESTR,QMSTR CHARACTER*8 SUBSET CHARACTER*3 TAG DIMENSION HDR(5),OBS(7,255),FCS(6,255),QMS(5,255),OES(5,255) LOGICAL QX(300),STUCK,PCHK DATA SUCCESS/'**************************************************', . '**************** STORE SUCCESSFUL ****************', . '**************************************************'/ DATA WARNING/'**************************************************', . '!! WARNING: NOT ALL OBS CAN BE STORED :WARNING !!', . '**************************************************'/ DATA HDSTR /'SID XOB YOB DHR TYP '/ DATA OBSTR /'POB QOB TOB ZOB UOB VOB CAT '/ DATA FCSTR /'PFC QFC TFC ZFC UFC VFC '/ DATA QMSTR /'PQM QQM TQM ZQM WQM '/ DATA OESTR /'POE QOE TOE ZOE WOE '/ DATA MAXREP / 200000 / DATA MAXLEV / 300000 / DATA ROG / 29.261 / DATA QX / 300*.TRUE./ C---------------------------------------------------------------------- PA(POB) = POB*(1.-.053) PB(POB) = POB*(1.+.053) C---------------------------------------------------------------------- C INITIALIZE C ---------- STUCK = .FALSE. NREC = 0 NREP = 0 NLEV = 0 MAXNLEV = 0 C OPEN THE DATA FILE AND CHECK THE DATE C ------------------------------------- PRINT* PRINT*,'****** READING DATA IN STORE ****** unit=',LUBFI CALL UFBREW(LUBFI,MSGS) CALL OPENBF(LUBFI,'IN',LUBFI) READ(LUDAT,'(8X,I8)',END=900,ERR=900) JDATE PRINT*,'****** READING DATA IN STORE ******' PRINT* C LOOP THROUGH THE INPUT MESSAGES - READ THE NEXT SUBSET C ------------------------------------------------------ CMIC$ DO ALL AUTOSCOPE CMIC$.SHARED (JDATE,NREC,NREP,NLEV,MAXNLEV,NREP0,NLEV0,QUALOB) CMIC$.SHARED (LUBFI,HDSTR,OBSTR,FCSTR,QMSTR,STUCK,SID,ROG) CMIC$.PRIVATE (SUBSET,IDATE,IREC,IFLD,HDR,OBS,FCS,QMS) CMIC$.PRIVATE (PCHK,ILEV,KLEV,KX,TAG,PACK,LEV) CMIC$.PRIVATE (POB,TOB,ZOB,UOB,VOB) CMIC$.PRIVATE (PFC,TFC,ZFC,UFC,VFC) CMIC$.PRIVATE (IQM,OB1,OB2,IOI,OER) CMIC$.PRIVATE (PQM,TQM,ZQM,WQQ) DO 50 KREC=1,MSGS CMIC$ GUARD 0 IF(IREADMG(LUBFI,SUBSET,IDATE).NE.0) call sabort('store - readmg') IF(IDATE.NE.JDATE) CALL SABORT('STORE - DATES DONT MATCH ') NREC = NREC+1 IREC = NREC IFLD = 0 CMIC$ ENDGUARD 0 11 DO WHILE(IREADSB(LUBFI).EQ.0) IFLD = IFLD+1 CALL UFBINT(LUBFI,HDR,5, 1,KLEV,HDSTR) CALL UFBINT(LUBFI,OBS,7,255,KLEV,OBSTR) CALL UFBINT(LUBFI,FCS,6,255,KLEV,FCSTR) CALL UFBINT(LUBFI,QMS,5,255,KLEV,QMSTR) C DECIDE WHETHER OIQC SHOULD PROCESS - SET A TAG C ---------------------------------------------- KX = HDR(5) IF(QX(KX)) THEN IF(KX.GE.180 .AND. KX.LE.199) THEN TAG = 'SFC' ELSEIF(KX.GE.160 .AND. KX.LE.179) THEN TAG = 'SAT' ELSEIF(KX.EQ.283) THEN TAG = 'SMI' ELSE TAG = 'UPA' ENDIF ELSE GOTO 11 ENDIF C SAVE THE CURRENT COUNTERS IN CASE OF DISASTER OR ARRAY OVERFLOW C --------------------------------------------------------------- CMIC$ GUARD 1 IF(.NOT.STUCK) THEN NREP0 = NREP NLEV0 = NLEV ENDIF C STORE THE HEADER C ---------------- NREP = NREP+1 ILEV = 0 IF(NREP.LE.MAXREP .AND. .NOT.STUCK) THEN SID(NREP) = HDR(1) PACK = WLN(NREP,HDR(2)) PACK = WLT(NREP,HDR(3)) PACK = WTM(NREP,HDR(4)) PACK = WKX(NREP,KX) PACK = WNL(NREP,ILEV) PACK = WLV(NREP,NLEV+1) PACK = WOC(NREP,TAG) PACK = WQL(NREP,QUALOB(KX)) PACK = WPRC(NREP,IREC) PACK = WPFD(NREP,IFLD) ELSE STUCK = .TRUE. ENDIF C STORE THE OB LEVELS C ------------------- DO 20 LEV=1,KLEV POB = OBS(1,LEV) TOB = OBS(3,LEV) ZOB = OBS(4,LEV) UOB = OBS(5,LEV) VOB = OBS(6,LEV) PFC = FCS(1,LEV) TFC = FCS(3,LEV) ZFC = FCS(4,LEV) UFC = FCS(5,LEV) VFC = FCS(6,LEV) PQM = QMS(1,LEV) TQM = QMS(3,LEV) ZQM = QMS(4,LEV) WQQ = QMS(5,LEV) IF(POB.LE.0 .OR. POB.GT.2000) PQM = 8 IF(PQM.EQ.8) TQM = MAX(PQM,TQM) IF(PQM.EQ.8) ZQM = MAX(PQM,ZQM) IF(PQM.EQ.8) WQQ = MAX(PQM,WQQ) IF(PQM.GT.8) GOTO 20 PCHK = TAG.EQ.'SFC' IQM = 0 OB1 = 0 OB2 = 0 IOI = 0 OER = 0 IF(TQM.LT.8) THEN IQM = TQM OB1 = TOB-TFC OB2 = ROG*OB1*LOG(PB(POB)/PA(POB)) OER = 1 IOI = 1 ELSEIF(TQM.EQ.8) THEN IOI = -2 ELSEIF(WQQ.LT.8) THEN IQM = WQQ OB1 = UOB-UFC OB2 = VOB-VFC OER = 1 IOI = 1 ELSEIF(WQQ.EQ.8) THEN IOI = -2 ELSEIF(.NOT.PCHK) THEN GOTO 20 ENDIF 15 NLEV = NLEV+1 ILEV = ILEV+1 IF(NLEV.LE.MAXLEV .AND. .NOT.STUCK) THEN IF(TAG.EQ.'SMI') THEN PACK = WMW(NREP,1,UFC) PACK = WMW(NREP,2,VFC) OB1 = 0. OB2 = SQRT(UOB**2+VOB**2) ENDIF IF(ABS(OB1).GT.3000.) OB1 = 3000. IF(ABS(OB2).GT.3000.) OB2 = 3000. PACK = WNL (NREP,ILEV) PACK = WHD (NLEV,NREP) PACK = WPLV(NLEV,LEV) PACK = WPR (NLEV,POB) PACK = WOB (NLEV,1,OB1) PACK = WOB (NLEV,2,OB2) PACK = WQM (NLEV,IQM) PACK = WOE (NLEV,OER) PACK = WOI (NLEV,IOI) ELSE STUCK = .TRUE. ENDIF IF(PCHK) THEN PCHK = .FALSE. IQM = MAX(PQM,ZQM) IF(IQM.LT.8) THEN IOI = 1 OB1 = 0 OB2 = ZOB-ZFC OER = 1 ELSEIF(IQM.EQ.8) THEN IOI = -2 OB1 = 0 OB2 = 0 OER = 0 ELSE IOI = 0 OB1 = 0 OB2 = 0 OER = 0 ENDIF GOTO 15 ENDIF 20 CONTINUE C KEEP TRACK OF MAX LEVELS C ------------------------ IF(ILEV.EQ.0) NREP = NREP-1 MAXNLEV = MAX(ILEV,MAXNLEV) CMIC$ ENDGUARD 1 C END OF READ LOOPS AND PARALLEL SECTION C -------------------------------------- ENDDO 50 ENDDO CALL CLOSBF(LUBFI) C ALL DONE STORING - RECORD THE OUTCOME AND EXIT C ---------------------------------------------- DO I=1,2 IF(I.EQ.1) LO = 6 IF(I.EQ.2) LO = 60 IF(NREP.LE.MAXREP .AND. NLEV.LE.MAXLEV) THEN WRITE(LO,'(A80)') WRITE(LO,'(A80)') SUCCESS WRITE(LO,101 ) MSGS WRITE(LO,102 ) NREP,NREP WRITE(LO,103 ) NLEV,NLEV WRITE(LO,104 ) MAXNLEV WRITE(LO,'(A80)') SUCCESS WRITE(LO,'(A80)') ELSE WRITE(LO,'(A80)') WRITE(LO,'(A80)') WARNING WRITE(LO,101 ) MSGS WRITE(LO,102 ) NREP,NREP0 WRITE(LO,103 ) NLEV,NLEV0 WRITE(LO,104 ) MAXNLEV WRITE(LO,'(A80)') WARNING WRITE(LO,'(A80)') NREP = NREP0 NLEV = NLEV0 ENDIF ENDDO RETURN 101 FORMAT('TOTAL RECORDS: ',I6) 102 FORMAT('TOTAL REPORTS: ',I6,' --- STORED: ',I6) 103 FORMAT('TOTAL OBS : ',I6,' --- STORED: ',I6) 104 FORMAT('MAX OBS/PRFLE: ',I6) 900 CALL SABORT('STORE - BAD OR MISSING NMCDATE FILE ') END C$$$ SUBPROGRAM DOCUMENTATION BLOCK C C SUBPROGRAM: PREPQM C PRGMMR: WOOLLEN ORG: NMC22 DATE: 90-11-06 C C ABSTRACT: REWRITES THE BUFRDA FILE WITH UPDATED QUALITY MARKS C C PROGRAM HISTORY LOG: C 90-11-06 J. WOOLLEN C C USAGE: C INPUT ARGUMENTS: NONE C LPIN - LOGICAL UNIT FOR PREPDA INPUT C LPOUT - LOGICAL UNIT FOR PREPQM OUTPUT C C INPUT FILES: C FORT.LPIN - PREPQM - PREPDA WITH UPDATED QUALITY MARKS C C OUTPUT ARGUMENTS: NONE C C OUTPUT FILES: C FORT.LPOUT - PREPQM - PREPDA WITH UPDATED QUALITY MARKS C C SUBPROGRAMS CALLED: C LIBRARY: C COMMON - ABORT C C REMARKS: C C C ATTRIBUTES: C LANGUAGE: CRAY FORTRAN C MACHINE: CRAY C C$$$ SUBROUTINE PREPQM(LUBIN,LUBOT) PARAMETER (MAXM=4000,MAXS=500) COMMON /OBS/ NREP,NLEV,SID( 200000 ),INOB( 200000 ),IMAP(360,181) c CHARACTER*30 OSTR(3),ESTR(3) CHARACTER*31 OSTR(3),ESTR(3) CHARACTER*8 SUBSET CHARACTER*3 TAG,FOC DIMENSION OBS(10,255),EVN(10,255),JREP(MAXM,MAXS) LOGICAL SFC,SMI DATA OSTR /'POB NUL','TOB QOB','UOB VOB'/ DATA ESTR /'POB PQM PPC PRC ', . 'TOB TQM TPC TRC QOB QQM QPC QRC', . 'UOB VOB WQM WPC WRC '/ DATA BMISS /10E10/ C----------------------------------------------------------------------- C----------------------------------------------------------------------- PRINT*,'******WRITE THE BUFRQM FILE*******' NREC = 0 C-CRA JREP = 0 C DIMENSION OBS(10,255),EVN(10,255),JREP(MAXM,MAXS) DO I=1,MAXM*MAXS JREP(I,1) = 0 ENDDO C SORT OUT THE FILE POSITIONS OF EACH OB IN STORAGE C ------------------------------------------------- DO IREP=1,NREP JREC = FPRC(IREP) JFLD = FPFD(IREP) IF(JREC.GT.MAXM) CALL SABORT('PREPQM - JREC GT MAXM') IF(JFLD.GT.MAXS) CALL SABORT('PREPQM - JFLD GT MAXS') JREP(JREC,JFLD) = IREP ENDDO C INITIALIZE THE INPUT AND OUTPUT FILES C ------------------------------------- CALL UFBREW(LUBIN,MSGS) CALL OPENBF(LUBIN,'IN',LUBIN) CALL OPENBF(LUBOT,'OUT',LUBIN) CALL UFBQCD(LUBIN,'OIQC',QCD) C READ MESSAGES, SUBSETS, DATA, EVENTS - TRY TO KEEP IT ALL STRAIGHT C ------------------------------------------------------------------ CMIC$ DO ALL AUTOSCOPE CMIC$.SHARED (LUBIN,LUBOT,NREC,QCD,ESTR,OSTR,JREP,BMISS) CMIC$.PRIVATE (SUBSET,IDATE,IREC,IFLD,OBS,EVN,UOB,VOB) CMIC$.PRIVATE (IREP,ILEV,JLEV,KLEV,JUPD,LEV,LEO) CMIC$.PRIVATE (TAG,KKX,SFC,SMI,QMA,QMB,OBM,L) CMIC$.PRIVATE (sid1,sid2) DO KREC=1,MSGS CMIC$ GUARD 0 IF(IREADMG(LUBIN,SUBSET,IDATE).NE.0)CALL SABORT('PREPQM - READMG') CALL OPENMB(LUBOT,SUBSET,IDATE) NREC = NREC+1 IREC = NREC IFLD = 0 CMIC$ ENDGUARD 0 DO WHILE(IREADSB(LUBIN).EQ.0) CALL UFBCPY(LUBIN,LUBOT) IFLD = IFLD+1 C FIND THIS OB IN STORAGE OR JUST WRITE IT BACK OUT C ------------------------------------------------- IF(JREP(IREC,IFLD).GT.0.AND.IREC.LE.MAXM.AND.IFLD.LE.MAXS) THEN IREP = JREP(IREC,IFLD) ILEV = FLV (IREP) KLEV = ILEV+NINT(FNL(IREP))-1 CALL UFBINT(LUBIN,SID1,1,1,LEV,'SID') SID2 = SID(IREP) IF(SID1.NE.SID2) CALL SABORT('SID1<>SID2') if(irec.ne.fprc(irep)) call sabort('irec.ne.fprc') if(ifld.ne.fpfd(irep)) call sabort('ifld.ne.fpfd') TAG = FOC(IREP) KKX = FKX(IREP) SFC = TAG.EQ.'SFC' SMI = TAG.EQ.'SMI' .AND. KKX/100.EQ.2 C CHECK ALL THE LEVELS OF THIS OB FOR EVENTS C ------------------------------------------ DO L=ILEV,KLEV JLEV = FPLV(L) QMB = FQM(L) OBM = FOI(L) QMA = QMB IF(OBM.EQ.-1 .AND. QMB.LE.3) QMA = QMB+4 IF(OBM.EQ. 1 .AND. QMB.GE.4) QMA = MIN(NINT(QMB)-4,3) IF(QMA.NE.QMB .OR. SMI) THEN C-CRA EVN = BMISS C DIMENSION OBS(10,255),EVN(10,255),JREP(MAXM,MAXS) DO I=1,10*255 EVN(I,1)= BMISS ENDDO IF( (SFC .AND. L.NE.ILEV)) JUPD = 1 IF(.NOT.(SFC .AND. L.NE.ILEV)) JUPD = KKX/100+1 CALL UFBINT(LUBIN,OBS,10,255,LEV,OSTR(JUPD)) IF(JUPD.EQ.1) THEN IF(OBS(1,JLEV).LT.BMISS) THEN EVN(1,JLEV) = OBS(1,JLEV) EVN(2,JLEV) = QMA EVN(3,JLEV) = QCD EVN(4,JLEV) = 0 CALL UFBINT(LUBOT,EVN,10,LEV,LEO,ESTR(JUPD)) ELSE CALL SABORT('PREPQM - MISSING PRESSURE') ENDIF ELSEIF(JUPD.EQ.2) THEN IF(OBS(1,JLEV).LT.BMISS) THEN EVN(1,JLEV) = OBS(1,JLEV) EVN(2,JLEV) = QMA EVN(3,JLEV) = QCD EVN(4,JLEV) = 0 ELSE print'(2a8,i8)',sid1,sid2,kkx CALL SABORT('PREPQM - MISSING TEMPERATURE') ENDIF IF(OBS(2,JLEV).LT.BMISS) THEN EVN(5,JLEV) = OBS(2,JLEV) EVN(6,JLEV) = QMA EVN(7,JLEV) = QCD EVN(8,JLEV) = 0 ENDIF CALL UFBINT(LUBOT,EVN,10,LEV,LEO,ESTR(JUPD)) ELSEIF(JUPD.EQ.3) THEN IF(TAG.NE.'SMI') THEN UOB = OBS(1,JLEV) VOB = OBS(2,JLEV) ELSE UOB = FMW(IREP,1) + FOB(ILEV,1) VOB = FMW(IREP,2) + FOB(ILEV,2) ENDIF IF(UOB.LT.BMISS .AND. VOB.LT.BMISS) THEN EVN(1,JLEV) = UOB EVN(2,JLEV) = VOB EVN(3,JLEV) = QMA EVN(4,JLEV) = QCD EVN(5,JLEV) = 0 CALL UFBINT(LUBOT,EVN,10,LEV,LEO,ESTR(JUPD)) ELSE print'(2a8,i8)',sid1,sid2,kkx CALL SABORT('PREPQM - MISSING WIND COMPONENT(S)') ENDIF ELSE CALL SABORT('PREPQM - JUPD OUT OF RANGE') ENDIF ENDIF ENDDO ENDIF CALL WRITSB(LUBOT) ENDDO CALL CLOSMG(LUBOT) ENDDO C HERE AT THE END C --------------- PRINT*,'******PREPQM PROCESSED ',MSGS,' BUFRDA RECORDS******' CALL CLOSBF(LUBIN) CALL CLOSBF(LUBOT) RETURN END SUBROUTINE SABORT(STRING) CHARACTER*(*) STRING WRITE(*,*) 'ABORT:',STRING STOP 8 END