! updated for bigger arrays ! _MAXM_ number of bufr message writen out ! _MAXS_ in each bufr message, number of subsection ! MAXCYL number of obs (same gross type) ! _MAXREP_ number of reports ! _MAXLEV_ number of level data ! _MAXPK_ something to with multiprocessing and packing, multiple of 64 ! #define _MAXPAK_ 12800 #define _MAXPAK_ 51200 #define _MAXM_ 15000 #define _MAXS_ 500 #define _MAXCYL_ 200000 #define _MAXREP_ 800000 ! want _MAXREP_ < _MAXLEV_ #define _MAXLEV_ 2000000 !$$$ MAIN PROGRAM DOCUMENTATION BLOCK ! . . . . ! MAIN PROGRAM: OIQCTEMP OPTIMAL INTERPOLATION QUALITY CONTROL ! PRGMMR: WOOLLEN ORG: W/NMC2 DATE: 90-11-06 ! ! ABSTRACT: QUALITY CONTROLS OBSERVATIONS TO BE USED IN THE GLOBAL ! ANALYSIS-FORECAST SYSTEM. OBSERVATIONAL DATASET GENERATED BY ! THE GESTEMP PREPROCESSOR IS INPUT TO OIQC WHICH THEN PERFORMS ! A "COMPLEX" QUALITY CONTROL (FROM GANDIN,1985) EMPLOYING ! MULTIVARIATE OI TO CHECK ALL OBSERVATIONS AGAINST NEARBY ! NEIGHBORS. SEVERAL ITERATIONS OF THE BASIC SCHEME ARE REQUIRED ! TO COMPLETE THE PROCESS. DURING EACH ITERATION ALL OBSERVATIONS ! ARE SUBJETED TO FIVE INTERPOLATION CHECKS. CHECK 1 INTERPOLATES ! COMPARISON VALUES FROM NEARBY TEMPERATURE DATA. CHECK 2 AND CHECK ! 3 INTERPOLATE COMPARISON VALUES FROM NEARBY ZONAL AND LONGITUDINAL ! WIND COMPONENTS RESPECTIVELY. CHECK 4 IS AN INTERPOLATION ! FROM THE NEAREST OBS IN THE SAME PROFILE AS THE ONE BEING CHECKED. ! A WEIGHTED COMBINATION OF THE OUTCOME OF THESE FOUR CHECKS IS ! PRODUCED ON THE BASIS OF FACTORS WHICH DETERMINE WHICH ! OF THE INTERPOLATIONS (IF ANY) MIGHT BE EXPECTED TO BE MORE ! RELEVANT THAN THE OTHERS; FOR EXAMPLE THE NUMBER AND LOCATION ! IN SPACE AND TIME OF OBSERVATIONS USED IN A GIVEN CHECK ! AS WELL AS RELATIVE QUALITY OF THE OBSERVATIONS INVOLVED IS USED TO ! INDICATE WHICH OF THE OUTCOMES SHOULD BE RELIED MORE HEAVILY UPON ! TO MAKE A DECISION IN A PARTICULAR CASE. CHECK 5 INTERPOLATES ! A VALUE FROM THE MULTIVARIATE SET OF OBS USED IN CHECKS 1-3, ! WHICH IS USED TO MEASURE HOW THE ANALYSIS WOULD DRAW (OR NOT DRAW) ! TO THE OB BEING CHECKED. THE TOLERANCE ALLOWED FOR DEVIATION FROM ! THE EXPECTED OUTCOME OF THE INTERPOLATION CHECKS IS PROPORTIONAL ! TO A MEASURE OF THE DRAW COMPUTED USING THE RESULT OF CHECK 5. ! DURING EACH ITERATION EACH OBSEVATION EITHER PASSES OR FAILS THE ! COMPLEX OF CHECKS JUST DESCRIBED. A "PASS" MEANS THAT OB MAY BE ! USED FOR CHECKING OBS DURING THE SUBSEQUENT ITERATION. A "FAIL" ! INDICATES THE OPPOSITE. THE SYSTEM IS ITERATED UNTIL IDENTICAL ! RESULTS IN TERMS OF OBS THAT PASS AND FAIL ARE OBTAINED IN TWO ! CONSEQUETIVE ITERATIONS UP TO FOUR COMPLETE ITERATIONS AT WHICH ! POINT A FINAL ARBITRATION PROCEDURE IS INVOKED TO RESOLVE A (SMALL) ! NUMBER OF AMBIGUOUS CASES WHICH ARE PREVENTING CONVERGENCE OF THE ! SYSTEM. ! ! PROGRAM HISTORY LOG: ! 90-11-06 J. WOOLLEN ORIGINAL VERSION FOR IMPLEMENTATION ! 90-04-23 J. WOOLLEN VERSION FOR IMPLEMENTATION W/ SPECTRAL OI ! 13-08-23 W. EBISUZAKI ADDED CHECK NOBS TOO BIG ! CHANGED TO SYMBOLIC VARIBLES ! 14-07-05 J. WOOLEN OLD: ABORT IF MISSING OBS ! NEW: UPDATE QM AND CONTINUE, INCREASE _MAXLEV_ ! ! USAGE: ! INPUT FILES: ! FORT.11 - PREPDA - OBSERVATION FILE ! FORT.12 - PREPQC - OBSERVATION INCREMENT FILE FROM GEST ! FORT.14 - OBERRS - OBSERVATION ERROR FILE ! FORT.16 - NMCDATE - NMC DATE FILE ! FORT.18 - OBPRT.IN - LIST OF IPOINTS FOR DETAILED PRINT ! FORT.20 - TOLLS.IN - TOLERANCE LEVELS (FOR EXPERIMENTATION) ! ! OUTPUT FILES: ! FORT.60 - OBCNT.OUT - SUMMARY OF OIQC RUN ! FORT.61 - TOSS.SFZ - DETAILED TOSSLIST FOR SURFACE HEIGHTS ! FORT.62 - TOSS.OUT - DETAILED TOSSLIST FOR TEMPS AND WINDS ! FORT.63 - TOSS.SAT - DETAILED TOSSLIST FOR SAT TEMPS ! FORT.64 - TOSS.SMI - DETAILED TOSSLIST FOR SSMI WIND SPEEDS ! FORT.65 - TOSSLIST - OFFICIAL OPERATIONAL TOSSLIST ! FORT.66 - PSFILE.OUT - LIST OF OBS TOSSED FOR BEING UNDERGROUND ! FORT.70 - PREPQM - OBSERVATIONAL DATA FOR OI INGEST ! FORT.71 - 80 - RESERVED FOR DIAGNOSTIC TESTING ! FORT.81 - OBOGRAM.OUT - SUMMARY OF QC RESULTS (ASCII) ! FORT.82 - OBOGRAM.BIN - SUMMARY OF QC RESULTS (BINARY) ! ! SUBPROGRAMS CALLED: (*AUTOTASKED) ! UNIQUE: - SETCON,RDOER,STORE,STOERR,PSFILE,MAKEAMAP, ! SEARCH,CHDIST,QCLOOP,QCPTS,SATPTS, ! *SELDAT,*QCOI,*COORS,*DRCTSL,*VSOLVE,*OUTPUT,*SATPUT, ! OBOGRAM,TOSSLIST,PREPQM ! ! LIBRARY: ! COMMON - W3LIB ! ! EXIT STATES: ! COND = 0 - SUCCESSFUL RUN ! ! IF ANY DISASTORIUS PROBLEM SITUATIONS ARE DETECTED THE PROGRAM ! WILL ABORT WITH A MESSAGE DETAILING THE CALAMITY. ! ! ! REMARKS: ! ! ATTRIBUTES: ! LANGUAGE: CRAY FORTRAN ! MACHINE: CRAY ! !$$$ !----------------------------------------------------------------------- !TREE GETFCER !TREE PILNLNP !TREE SATPTS !TREE SATPUT !----------------------------------------------------------------------- PROGRAM OIQCTEMP COMMON /OBS/ NREP,NLEV,SID(_MAXREP_),INOB(_MAXREP_),IMAP(360,181) COMMON /CKLIST / NDD,INDD(_MAXLEV_) COMMON /TOSSUNIT/ ITUNIT CHARACTER*3 FOC EXTERNAL QCPTS,SATPTS DATA LUDAT /11/ DATA LUBFI /14/ DATA LUQCE /17/ DATA LUBFQ /70/ !----------------------------------------------------------------------- ! UNIFIED QC ANALYSIS SETUP PROCEDURE IS AS FOLLOWS: ! ------------------------------------------------- ! 0) W3LOG - WRITE AN OPERATIONS START MESSAGE ! 1) SETCON - SET CONSTANTS AND PARAMETERS ! 2) STORE - READ IN FERR (OB/INCREMENT) FILE AND STORE IT ! 3) RDOBER - READ IN THE OBSERVATION ERRORS ! 4) STOERR - STORE THE OB AND FC ERRORS WITH EACH REPORT LEVEL ! 5) MAKEAMAP - MAKE AN OB POINTER LIST SORTED BY LAT-LON !----------------------------------------------------------------------- ! CALL W3LOG('$S','$M') CALL SETCON CALL STORE(LUDAT,LUBFI) CALL RDOBER (LUQCE) CALL STOERR CALL MAKEAMAP ! OIQC CHECKS OF SURFACE HEIGHT/PRESSURE DATA ! ------------------------------------------- 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 ! OIQC CHECKS OF UPA AND SFC TEMP AND WIND DATA ! ------------------------------------------- NDD = 0 NDDALL = 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 NDDALL = NDDALL + FNL(IREP) 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' WRITE(6,*)' NDDALL ',NDDALL WRITE(60,*)' NDDALL ',NDDALL ITUNIT = 62 WRITE(6 ,*) ' >>call search' WRITE(60,*) ' >>call search' CALL SEARCH WRITE(6 ,*) ' >>call qcloop' WRITE(60,*) ' >>call qcloop' CALL QCLOOP(QCPTS,4) ENDIF ! OIQC CHECKS OF SATEM DATA ! ------------------------- WRITE(6 ,*)' >>nrep=',nrep WRITE(60,*)' >>nrep=',nrep 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 ! DIRECTION ASSIGNMENT AND OIQC CHECKS FOR SSMI DATA ! -------------------------------------------------- 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 ! COMPUTE STATISTICS, WRITE LISTS, UPDATE QUALITY MARKS ! ----------------------------------------------------- CALL OBOGRAM (11,81,82) CALL TOSSLIST (65) CALL PREPQM (LUBFI,LUBFQ) ! WRITE AN OPERATIONS END MESSAGE ! ------------------------------- ! CALL W3LOG('$E') STOP END !$$$ SUBPROGRAM DOCUMENTATION BLOCK ! ! SUBPROGRAM: CHDIST ! PRGMMR: WOOLLEN ORG: NMC22 DATE: 90-11-06 ! ! ABSTRACT: ! COMPUTES CHORD LENGTH DISTANCE FROM A FIXED ! POINT TO AN ARRAY OF POINTS USING THE FORMULA: ! S**2/2 = 1 - COS(Y1-Y2) + COS(Y1)*COS(Y2)*(1-COS(X1-X2)). ! ALSO THE DIRECTION OF EACH POINT IS COMPUTED WITH RESPECT TO ! THE FIXED POINT AND RETURNED AS WELL. ! ! PROGRAM HISTORY LOG: ! 90-11-06 J. WOOLLEN ! ! USAGE: ! INPUT ARGUMENTS: ! X1 - X-COORDINATE (LONGITUDE) OF FIXED POINT ! Y1 - Y-COORDINATE (LATITUDE ) OF FIXED POINT ! X2 - X-COORDINATES (LONGITUDE) FOR SET OF POINTS ! Y2 - Y-COORDINATES (LATITUDE ) FOR SET OF POINTS ! NP - NUMBER OF OUTPUTS REQUESTED ! ! OUTPUT ARGUMENTS: ! DIST - CHORD LENGTH DISTANCES FROM FIXED POINT (KM) ! DIRN - DIRECTION FROM FIXED POINT (DEG) ! ! SUBPROGRAMS CALLED: NONE ! ! EXIT STATES: NONE ! ! REMARKS: ! ! ! ATTRIBUTES: ! LANGUAGE: CRAY FORTRAN ! MACHINE: CRAY ! !$$$ SUBROUTINE CHDIST(X1,Y1,X2,Y2,DIST,DIRN,NP) DIMENSION X2(NP),Y2(NP),DIST(NP),DIRN(NP) DATA PI180/.0174532 /,RADE/6371./ !---------------------------------------------------------------------- !---------------------------------------------------------------------- IF(NP.EQ.0) RETURN ! COMPUTE THE DISTANCE ! -------------------- 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 ! COMPUTE DIRECTIONS ! ------------------ 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 !$$$ SUBPROGRAM DOCUMENTATION BLOCK ! ! SUBPROGRAM: COORS ! PRGMMR: WOOLLEN ORG: NMC22 DATE: 90-11-06 ! ! ABSTRACT: COMPUTES CHORD LENGTH DISTANCE FOR A MATRIX OF ! LOCATIONS <(X,Y)I,(X,Y)J> USING THE NORMAL ANGLE FORMULA: ! S**2/2 = 1 - COS(Y1-Y2) + COS(Y1)*COS(Y2)*(1-COS(X1-X2)) ! DERIVATIVES OF DISTANCE WITH RESPECT TO LAT AND LON ARE ALSO ! COMPUTED AND COMBINED WITH APPROPRIATE CORRELATION FUNCTIONS ! AND DERIVATIVES WITH RESPECT TO DISTANCE TO FORM THE MULTI- ! VARIATE CORRELATIONS FOR THE HEIGHT WIND ANALYSIS. ! ! PROGRAM HISTORY LOG: ! 90-11-06 J. WOOLLEN ! ! USAGE: ! INPUT ARGUMENTS: ! NP - NUMBER OF OUTPUTS REQUESTED ! CH - VECTOR OF CONSTANTS FOR THE GAUSSIAN LENGTH SCALE ! X1 - X-COORDINATES (LONGITUDE) FOR POINT 1 ! Y1 - Y-COORDINATES (LATITUDE ) FOR POINT 1 ! K1 - VARIABLE TYPE AT POINT 1 (I.E. 1-Z , 2-U , 3-V) ! X2 - X-COORDINATES (LONGITUDE) FOR POINT 2 ! Y2 - Y-COORDINATES (LATITUDE ) FOR POINT 2 ! K2 - VARIABLE TYPE AT POINT 2 (I.E. 1-Z , 2-U , 3-V) ! FL - DECOUPLING FACTOR ! TZ - HEIGHT-TEMPERATURE FACTOR ! ! OUTPUT ARGUMENTS: ! D - DISTANCES BETWEEN ARRAYS OF POINTS IN RADIANS ! C1 - CORRELATION COEFFICIENTS AS SPECIFIED ! C2 - CORRELATION COEFFICIENTS AS SPECIFIED ! ! SUBPROGRAMS CALLED: NONE ! ! EXIT STATES: NONE ! ! REMARKS: THIS PROGRAM COMPUTES A GAUSSIAN CORRELATION WITH ! LENGTH SCALE GIVEN BY THE CHTWV INPUT ARGUMENT. ! ! ! ATTRIBUTES: ! LANGUAGE: CRAY FORTRAN ! MACHINE: CRAY ! !$$$ 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/ !---------------------------------------------------------------------- !---------------------------------------------------------------------- IF(NP.EQ.0) RETURN ! LOOP OVER SET OF INPUT POINTS ! ----------------------------- DO 20 I=1,NP C1(I) = 0. C2(I) = 0. ! DETERMINE CORRELATION TYPE ! -------------------------- 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 ! COMPUTE THE MATRIX OF SINES AND COSINES ! --------------------------------------- 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) ! COMPUTE THE NORMAL ANGLE ! ------------------------ S = 1.0-COSDY+COSY1*COSY2*(1.0-COSDX) S = SQRT(2.*S) ! CALCULATE THE FIRST PARTIALS WITH RESPECT TO X1 AND Y1 ! ------------------------------------------------------ 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 ! CALCULATE THE MIXED PARTIAL DERIVATIVES ! --------------------------------------- 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) ! COMPUTE THE VARIOUS CORRELATIONS ! -------------------------------- 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/ !---------------------------------------------------------------------- !---------------------------------------------------------------------- IF(NP.EQ.0) RETURN ! LOOP OVER SET OF INPUT POINTS ! ----------------------------- DO 20 I=1,NP C1(I) = 0. C2(I) = 0. ! DETERMINE CORRELATION TYPE ! -------------------------- 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 ! COMPUTE THE MATRIX OF SINES AND COSINES ! --------------------------------------- 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) ! COMPUTE THE NORMAL ANGLE ! ------------------------ S = 1.0-COSDY+COSY1*COSY2*(1.0-COSDX) S = SQRT(2.*S) ! CALCULATE THE FIRST PARTIALS WITH RESPECT TO X1 AND Y1 ! ------------------------------------------------------ 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 ! CALCULATE THE MIXED PARTIAL DERIVATIVES ! --------------------------------------- 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) ! COMPUTE THE VARIOUS CORRELATIONS ! -------------------------------- 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 !$$$ SUBPROGRAM DOCUMENTATION BLOCK ! ! SUBPROGRAM: DRCTSL ! PRGMMR: WOOLLEN ORG: NMC22 DATE: 90-11-06 ! ! ABSTRACT: DRIVER FOR CHOLESKY TYPE LINEAR EQUATION SOLVER. ! ! PROGRAM HISTORY LOG: ! 90-11-06 J. WOOLLEN ! ! USAGE: ! INPUT ARGUMENTS: ! FAALL - ARRAY OF SYMMETRIC MATRIXES ! RAALL - ARRAY OF MATRIX RIGHT HAND SIDES ! NDIM - ARRAY OF MATRIX RANKS ! MAXDIM - MAXIMUM RANK MATRIX IN MATRIX ARRAY ! NXXYY - NUMBER OF MATRIXES IN MATRIX ARRAY ! NFT - NUMBER OF RIGHT HAND SIDE VECTORS PER MATRIX ! OUTPUT ARGUMENTS: ! RAALL - ARRAY OF MATRIX SOLUTIONS ! DOTPRD - ARRAY OF DOT PRODUCTS OF RHS VECTORS WITH MATRIX ! SOLUTIONS ! ! SUBPROGRAMS CALLED: ! UNIQUE: - VSOLVE ! ! REMARKS: ILL CONDITIONED OR NON POSITIVE DEFINITE MATRIXES ARE ! IDENTIFIED BY DOT PRODUCTS GT 1 OR LT 0 OR BY A MESSAGE ! FROM VSOLVE. FIVE ATTEMPTS ARE MADE TO SOLVE BAD ONES, ! BY RIDGE REGRESSION, AFTER WHICH A NULL SOLUTION IS RETURNED. ! ! ATTRIBUTES: ! LANGUAGE: CRAY FORTRAN ! MACHINE: CRAY ! !$$$ !----------------------------------------------------------------------- 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/ !---------------------------------------------------------------------- !---------------------------------------------------------------------- ! LOOP FOR POSSIBILITY OF BAD MATRIXS ! ----------------------------------- DO 50 ITRY=1,6 NBAD = 0 ! FIND THE LARGEST MATRIX IN THE STACK ! ------------------------------------ MAXDIM = 0 DO 5 I=1,NXXYY MAXDIM = MAX(MAXDIM,NDIM(I)) 5 CONTINUE IF(MAXDIM.EQ.0) RETURN ! INITIALIZE THE WORKING ARRAYS ! ----------------------------- 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 ! CALL THE MATRIX SOLVER ! ---------------------- CALL VSOLVE (A,B,NDIM,BAD,NXXYY,MAXDIM) ! MAKE THE DOT PRODUCTS OF SOLUTIONS WITH RIGHT HAND SIDES ! -------------------------------------------------------- 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 ! NORMALIZE THE WEIGHTS IF THEY ARE TOO BIG ! ----------------------------------------- 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 ! CHECK FOR BAD ONES - DO IT AGAIN IF NECESSARY ! --------------------------------------------- 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 ! COPY SOLUTIONS INTO OUTPUT ARRAY - ZERO BAD ONES OUT ! ---------------------------------------------------- 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 !$$$ SUBPROGRAM DOCUMENTATION BLOCK ! . . . . ! SUBPROGRAM: GETFCER GET FORECAST ERRORS ! PRGMMR: WOOLLEN ORG: NMC22 DATE: 92-07-29 ! ! ABSTRACT: FUNCTION WHICH RETURNS A FORECAST ERROR IN THE FORCAST ! HEIGHT OR WIND FIELD FOR A GIVEN LAT, LON, PRESSURE. ENTRY POINTS ! GETZER AND GETWER RETURN ERROR VALUES; GETFCER RETURNS 0. ! ! PROGRAM HISTORY LOG: ! 92-07-29 J. WOOLLEN ! ! USAGE: CALL GETZER(SLAT,SLON,PRES) - (RETURNS HEIGHT ERROR) ! CALL GETWER(SLAT,SLON,PRES) - (RETURNS WIND ERROR) ! ! INPUT ARGUMENTS: ! SLAT - LATITUDE ! SLON - LONGITUDE ! PRES - PRESSURE ! ! FUNCTION RETURN VALUE: ! GETZER - FORECAST HEIGHT ERROR ! GETWER - FORECAST WIND ERROR ! ! SUBPROGRAMS CALLED: ! ! EXIT STATES: ! ! REMARKS: ! ! ! ATTRIBUTES: ! LANGUAGE: CRAY FORTRAN ! MACHINE: CRAY ! !----------------------------------------------------------------------- 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 / !--------------------------------------------------------------------- !--------------------------------------------------------------------- GETFCER = 0. RETURN ENTRY GETZER(SLAT,SLON,PRES) WIND = .FALSE. GOTO 10 ENTRY GETWER(SLAT,SLON,PRES) WIND = .TRUE. ! COMPUTE A HEIGHT ERROR PROPORTIONAL TO RAWINSONDE ERROR ! ------------------------------------------------------- 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 ! COMPUTE THE GEOSTROPHIC WIND ERROR FOR GETWER ! --------------------------------------------- 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 !----------------------------------------------------------------------- ! !----------------------------------------------------------------------- SUBROUTINE MAKEAMAP COMMON /OBS/ NREP,NLEV,SID(_MAXREP_),INOB(_MAXREP_),IMAP(360,181) COMMON /CKLIST/ NDD,INDD(_MAXLEV_) !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! MAKE SURE NONE OF THE LONGITUDES EQUAL 360. ! ------------------------------------------- DO 1 N=1,NREP IF(FLN(N).EQ.360.) PACK = WLN(N,0.) 1 CONTINUE ! LONGITUDE SORT ! -------------- NN = 0 DO 10 NX=1,360 DO 10 N=1,NREP IX = FLN(N)+1. IF(IX.EQ.NX) THEN NN = NN+1 if (NN.GT._MAXLEV_) CALL SABORT('NN > _MAXLEV_') INDD(NN) = N ENDIF 10 CONTINUE ! LATITUDE SORT ! ------------- 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 if (NOB.GT._MAXLEV_) CALL SABORT('NOB > _MAXREP_') INOB(NOB) = INDD(N) ENDIF 20 CONTINUE ! CHECK TO MAKE SURE ALL REPS ARE ACCOUNTED FOR ! --------------------------------------------- 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 ! INITIALIZE IMAP ARRAY ! --------------------- IMAP = 0 ! 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 ! FILL GAPS IN IMAP FOR CONTINUITY (BACKWARDS) ! -------------------------------------------- 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 !$$$ SUBPROGRAM DOCUMENTATION BLOCK ! ! SUBPROGRAM: OBOGRAM ! PRGMMR: WOOLLEN ORG: NMC22 DATE: 90-11-06 ! ! ABSTRACT: EXAMINES THE RESULTS OF FINAL QC DECISIONS AND MAKES ! TABLES OF STATISTICS FOR ACCEPTED AND REJECTED REPORTS BY ! REPORT TYPE AND LEVEL. ! ! PROGRAM HISTORY LOG: ! 90-11-06 J. WOOLLEN ! ! USAGE: CALL OBOGRAM(LASC,LBIN) ! ! INPUT ARGUMENTS: ! LASC - LOGICAL UNIT FOR ASCII OUTPUT ! LBIN - LOGICAL UNIT FOR BINARY OUTPUT ! ! OUTPUT ARGUMENTS: NONE ! ! OUTPUT FILES: ! FORT.LASC - OBOGRAM.OUT - SUMMARY OF QC RESULTS (ASCII) ! FORT.LBIN - OBOGRAM.BIN - SUMMARY OF QC RESULTS (BINARY) ! ! SUBPROGRAMS CALLED: ! LIBRARY: ! COMMON - ABORT ! ! REMARKS: ! ! ! ATTRIBUTES: ! LANGUAGE: CRAY FORTRAN ! MACHINE: CRAY ! !$$$ !----------------------------------------------------------------------- SUBROUTINE OBOGRAM(LDAT,LASC,LBIN) COMMON /OBS/ NREP,NLEV,SID(_MAXREP_),INOB(_MAXREP_),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./ !---------------------------------------------------------------------- !----------------------------------------------------------------------- REWIND (LDAT) READ(LDAT,1) NMCDATE 1 FORMAT(8X,I8) WRITE(LASC,*) WRITE(LASC,*)' STATS FOR OBS CHECKED ',NMCDATE WRITE(LASC,*) ! START UP ! -------- ! LOGICAL QOBS(300,6) QOBS = .FALSE. ! DIMENSION STAT(300,15,6,2,2,3),IOBS(300),RZONE(6) STAT = 0 ! GO THROUGH THE BIG BUFFER TO COUNT OBS FOR OBOGRAM ! -------------------------------------------------- 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 ! GO THROUGH LEVELS FOR ACCEPT REJECT STATISTICS ! ---------------------------------------------- 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 ! COMPUTE THE AVERAGES FOR MEAN AND RMS'S (VECTOR SUMS FOR WINDS) ! --------------------------------------------------------------- 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 ! PRINT THE TOTAL, MEAN, AND RMS FOR ACCEPTS AND REJECTS ! ------------------------------------------------------ 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 ! PRINT THE TOTAL REJECTED OB BY TYPE ZONE LEVEL ! ---------------------------------------------- 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 ! DIMENSION STAT(300,15,6,2,2,3),IOBS(300),RZONE(6) !-CRA RZONE = 0 DO LB=1,6 RZONE(LB) = 0 ENDDO ! 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 ! WRITE A BINARY VERSION OF THE OBOGRAM ! ------------------------------------- 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) ! FORMAT STATEMENTS FOR THE OBOGRAM ! --------------------------------- 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 !$$$ SUBPROGRAM DOCUMENTATION BLOCK ! ! SUBPROGRAM: OUTPUT ! PRGMMR: WOOLLEN ORG: NMC22 DATE: 90-11-06 ! ! ABSTRACT: EXAMINES OUTPUTS FROM OI QUALITY CONTROL INTERPOLATION ! AND MAKES DECISION WHETHER AN OB PASSES OR FAILS THE CHECKS. ! ALSO PROVIDES FOR A DETAILED PRINT TRACING ALL DECISIONS MADE ! IN REGARD TO SPECIFIED OBSERVATIONS IN THE QC PROCESS. ! ! PROGRAM HISTORY LOG: ! 90-11-06 J. WOOLLEN ! ! USAGE: CALL OUTPUT ! ! INPUT ARGUMENTS: NONE ! ! OUTPUT ARGUMENTS: IN COMMON /OBMARK/ ! ! SUBPROGRAMS CALLED: ! LIBRARY: ! COMMON - ABORT ! ! REMARKS: ! ! ! ATTRIBUTES: ! LANGUAGE: CRAY FORTRAN ! MACHINE: CRAY ! !$$$ SUBROUTINE OUTPUT !-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(_MAXREP_),INOB(_MAXREP_),IMAP(360,181) !-CRA COMMON /TOLLER / ITER,TB,TZ,TSAT,TWND,TJET !-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/ !---------------------------------------------------------------------- !---------------------------------------------------------------------- !-CRA D = 0. !-CRA DU = 0. !-CRA DV = 0. !-CRA TWT = 0. !-CRA TWU = 0. !-CRA TWV = 0. ! DIMENSION IOMF(3),D( 4 ),DRES( 4 ),NDIM( 4 ),STR(2) ! 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 ! ! LOOP OVER EACH OB CHECKED IN THIS GROUP ! --------------------------------------- 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) ! LOOK FOR A KEEP FLAG ! -------------------- IQM = FQM(I) IF(IQM.EQ.0) THEN KEEP = .TRUE. CKEEP = ' KEEP' ELSE KEEP = .FALSE. CKEEP = ' ' ENDIF ! LOOP OVER THE VARIABLE(S) IN EACH OB ! ------------------------------------ DO 45 K=J,2*J-1 KK = K-J+1 OMF = FOB(I,IOMF(K)) AMFSUM = 0. SUMWTS = 0. QC = 0. ! LOOP OVER THE OUTCOME OF EACH TYPE OF CHECK (MICRO DECISION MAKER) ! ------------------------------------------------------------------ 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) !-CRA TWT = TWT/SUMWTS DO IJ=1, 4 TWT(IJ) = TWT(IJ)/SUMWTS ENDDO ! FOR WINDS CHECK THE VECTOR SUM OF AGAINST THE REFERENCE ! ------------------------------------------------------- IF(K.EQ.1) THEN OMA = OMF - AMFSUM ELSE IF(K.EQ.2) THEN !-CRA DU = D !-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 !-CRA DV = D !-CRA TWV = TWT DO IJ=1, 4 DV(IJ) = D(IJ) TWV(IJ) = TWT(IJ) ENDDO QCV = QC OMAV = OMF - AMFSUM OMFV = OMF !-CRA D = SQRT(DU**2+DV**2) !-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 ! ADJUST THE TOLERANCE BY LOCATION AND VARIABLE ! --------------------------------------------- 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.) ! COMPUTE THE MULTIVARIATE REDUCTION RATIO AND ADJUST THE SUM ! ----------------------------------------------------------- IF(OMF.NE.0.) THEN QCPC = MIN(MAX(ABS(OMA/OMF),.1),1.) TTT = TT/QCPC ELSE TTT = TT ENDIF ! COMPUTE THE NORMALIZED DEVIATION FROM THE FIRST GUESS ! ----------------------------------------------------- SDOMF = SQRT(OMF**2/(OBER**2+FER**2)) ! MAKE A MESO DECISION FOR THIS ITERATION ! --------------------------------------- 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 ! THIS WRITES DETAILED INFO ABOUT AN OB ! ------------------------------------- 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,')' ), 100 FORMAT(A6,I4,3F8.2,A4,2F8.2,4(' *',F6.1,' (',F5.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 !$$$ SUBPROGRAM DOCUMENTATION BLOCK ! . . . . ! SUBPROGRAM: PILNLNP PRESSURE INTERPOLATION LINEAR IN LOG P ! PRGMMR: WOOLLEN ORG: W/NMC22 DATE: 92-07-29 ! ! ABSTRACT: FUNCTION WHICH GIVEN AN PROFILE OF DESCENDING PRESSURES ! AND A PROFILE OF QUANTITIES VALID AT THOSE PRESSURES, RETURNS ! AN INTERPOLATED QUANTITY VALID AT A GIVEN ARBITRARY PRESSURE. ! ! PROGRAM HISTORY LOG: ! 86-03-21 G. DIMEGO ! 88-11-24 D. DEAVEN RECODED FOR CYBER 205 ! ! USAGE: X = PILNLNP(P,PARAY,QARAY,KMAX) ! ! INPUT ARGUMENTS: ! P - PRESSURE TO INTERPOLATE TO ! PARAY - GIVEN PRESSURE PROFILE ! QARAY - QUANTITIES VALID FOR PRESSURE PROFILE ! KMAX - LENGTH OF PROFILE ! ! FUNCTION RETURN VALUE ! PILNLNP - INTERPOLATED QUANTITY ! ! ATTRIBUTES: ! LANGUAGE: CDC FORTRAN200 ! MACHINE: CYBER ! !$$$ 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 !$$$ SUBPROGRAM DOCUMENTATION BLOCK ! ! SUBPROGRAM: QCLOOP ! PRGMMR: WOOLLEN ORG: NMC22 DATE: 90-11-06 ! ! ABSTRACT: ITERATES THE QC DECSION MAKING PROCESS UNTIL ALL DECISIONS ! IN TWO CONSEQUETIVE ITERATIONS ARE THE SAME. AFTER FOUR ITERATIONS ! ONLY REPORTS THAT FAIL THE CHECKS ARE CHECKED FURTHER. ! ! PROGRAM HISTORY LOG: ! 90-11-06 J. WOOLLEN ! ! USAGE: CALL QCLOOP ! ! INPUT ARGUMENTS: NONE ! ! OUTPUT ARGUMENTS: NONE ! ! SUBPROGRAMS CALLED: ! UNIQUE: - QCPTS ! ! LIBRARY: ! COMMON - ABORT ! ! REMARKS: ! ! ! ATTRIBUTES: ! LANGUAGE: CRAY FORTRAN ! MACHINE: CRAY ! !$$$ SUBROUTINE QCLOOP(QCPTS,MAXTRY) COMMON /FPLIST / NFL,INFL(_MAXLEV_) COMMON /CKLIST / NDD,INDD(_MAXLEV_) !-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 !---------------------------------------------------------------------- !----------------------------------------------------------------------- IF(NDD.EQ.0) THEN PRINT*,'QCLOOP - NDD = 0' RETURN ENDIF ITRY = 1 ITER = 1 IFLIP = 0 ! SET UP THE TOLERANCES ! --------------------- 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 ! DO THE QCOI CHECKS ON LIST IN COMMON /CKLIST/ ! --------------------------------------------- 10 REWIND(ITUNIT) CALL QCPTS ! SET OB ERRS TO REFLECT LAST ROUND OF CHECKS - RECORD FLIPS ! ---------------------------------------------------------- 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 ! CHECK STATUS AND WRITE ITERATION REPORT ! --------------------------------------- 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) ! FINISHED YET? ! ------------- 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 ! FINAL REPORT AND AWAY ! --------------------- NBAD = NBB WRITE(6 ,*)' TOTAL BAD OBS=',NBAD WRITE(60,*)' TOTAL BAD OBS=',NBAD ITRY = 0 ITER = 0 RETURN END !$$$ SUBPROGRAM DOCUMENTATION BLOCK ! . . . . ! SUBPROGRAM: QCOI THE QC ANALYSIS KERNAL ! PRGMMR: DEAVEN ORG: W/NMC22 DATE: 88-02-05 ! PRGMMR: WOOLLEN ORG: W/NMC22 DATE: 90-03-24 ! ! ABSTRACT: COMPUTES AND SETS UP THE COEFFICIENT MATRIX AND RIGHT ! HAND SIDES. SOLVES THE MATRIX PROBLEMS. APPLIES THE COMPUTED ! WEIGHTS TO THE RESIDUALS TO OBTAIN THE ANALYSIS INCREMENTS AND ! ANALYSIS ERRORS. ! ! USAGE: ! INPUT ARGUMENTS: IN COMMON /HVECT/ ! NXXYY - NUMBER OF OBS IN PACKAGE ! XX - LONGITUDES OF OBS IN PACKAGE ! YY - LATITUDES OF OBS IN PACKAGE ! PP - PRESSURES OF OBS IN PACKAGE ! II - INDEXES OF OBS IN PACKAGE ! LL - REPORT LEVEL OF OBS IN PACKAGE ! FESD - FORECAST ERRORS AT OB LOCATIONS ! MKOBS - LIST OF OB INDEXES USED IN EACH INTERPOLATION ! NDIM - NUMBER OF OBS USED IN EACH INTERPOLATION ! JVMAX - NUMBER OF INTERPOLATING OBS FOR EACH CHECK IN ! EACH PACKAGE ! ! OUTPUT ARGUMENTS: IN COMMON /HVECT/ ! AA - INTERPOLATED VALUE OF OBS IN PACKAGE ! AESD - ANALYSIS ERRORS AT OB LOCATIONS ! ! SUBPROGRAMS CALLED: ! UNIQUE: - COORS,DRCTSL ! ! LIBRARY: ! COMMON - ABORT ! ! EXIT STATES: ! COND = 9 - ERROR DETECTED (SEE CONSOLE OUTPUT) ! ! REMARKS: ! ! ATTRIBUTES: ! LANGUAGE: CRAY FORTRAN ! MACHINE: CRAY ! !$$$ !----------------------------------------------------------------------- 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) !-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/ !----------------------------------------------------------------------- !----------------------------------------------------------------------- IF(NXXYY.EQ.0) RETURN ! . MDIM( 64 , 4 ),AA( 64 , 4 ,2), ! . AE( 64 , 4 ,2),WT( 64 , 10 , 4 ,2), AA = 0. AE = 0. WT = 0. ! LOOP THRU EACH TYPE OF QC CHECK ! ------------------------------- DO 100 IQ=1,NQC IVMAX = 0 ! CLEAR THE MATRIX ARRAYS ! ----------------------- ! DIMENSION FA(NMX,NN), RA(NMX,N,2), DP(NMX,2), IVAR(3) FA = 0. RA = 0. DP = 0. ! SET UP VECTORS WITH INFORMATION ABOUT EACH OB TO BE CHECKED ! ----------------------------------------------------------- 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 ! COMPUTE THE RIGHT HAND SIDES ! ---------------------------- 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 ! COMPUTE THE OFF-DIAGONAL MATRIX TERMS ! ------------------------------------- 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) ! STORE OFF-DIAGONAL MATRIX TERMS IN THE SOLVER ARRAY ! --------------------------------------------------- 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 ! COMPUTE AND STORE THE ON-DIAGONAL MATRIX TERMS ! ---------------------------------------------- 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 ! SOLVE THE LINEAR INTERPOLATION EQUATIONS ! ---------------------------------------- CALL DRCTSL(FA,RA,DP,MDIM(1,IQ),NXXYY) ! COMPUTE AND STORE ANALYSIS INCREMENTS AND ERRORS ! ------------------------------------------------ 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 ! END OF LOOP OVER EACH QC CHECK ! ------------------------------ 100 CONTINUE RETURN END !$$$ SUBPROGRAM DOCUMENTATION BLOCK ! ! SUBPROGRAM: QCPTS ! PRGMMR: WOOLLEN ORG: NMC22 DATE: 90-11-06 ! ! ABSTRACT: PUTS TOGETHER PACKAGES OF OBSERVATIONS TO BE CHECKED ! IN A GIVEN ITERATION AND INVOKES AUTOTASKING OF THE DATA ! SELECTION, OI, AND DECISION MAKING PROCEDURES. ! ! PROGRAM HISTORY LOG: ! 90-11-06 J. WOOLLEN ! ! USAGE: ! INPUT ARGUMENTS: IN COMMON /CKLIST/ ! NDD - NUMBER OF OBSERVATIONS TO CHECK ! INDD - ARRAY OF OBSERVATION INDEXES TO CHECK ! ! OUTPUT ARGUMENTS: NONE ! ! SUBPROGRAMS CALLED: ! UNIQUE: - SELDAT,QCOI,OUTPUT ! ! LIBRARY: ! COMMON - ABORT ! ! REMARKS: ! ! ! ATTRIBUTES: ! LANGUAGE: CRAY FORTRAN ! MACHINE: CRAY ! !$$$ SUBROUTINE QCPTS COMMON /CKLIST/ NDD,INDD(_MAXLEV_ ) DIMENSION IND( _MAXPAK_ ),NXY( _MAXPAK_ ) DATA NMX / 64 / !---------------------------------------------------------------------- !---------------------------------------------------------------------- ! ASSEMBLE LISTS OF OBS FOR MULTI PROCESSING ! ------------------------------------------ IND(1) = 1 IPK = 1 ! DIMENSION IND( _MAXPAK_ ),NXY( _MAXPAK_ ) NXY = 0 DO 5 I=1,NDD ILEV = INDD(I) IF(NXY(IPK)+1.GT.NMX) THEN IF(IPK+1.GT. _MAXPAK_) CALL SABORT('QCPTS - IPK > MAXPAK') IPK = IPK+1 IND(IPK) = I NXY(IPK) = 1 ELSE NXY(IPK) = NXY(IPK)+1 ENDIF 5 CONTINUE ! MULTI-PROCESS QC OB LISTS ! ------------------------- !MIC$ DOALL AUTOSCOPE DO 10 I=1,IPK CALL SELDAT(IND(I),NXY(I)) CALL QCOI CALL OUTPUT 10 CONTINUE RETURN END !$$$ SUBPROGRAM DOCUMENTATION BLOCK ! . . . . ! SUBPROGRAM: RDOBER ASSIGN OBSERVATIONAL ERRORS ! PRGMMR: DEAVEN ORG: W/NMC22 DATE: 88-02-10 ! ! ABSTRACT: INITIALIZES ARRAY 'QIKE' WITH OBSERVATION ERRORS FOR ! ALL KINDS OF OBSERVATION TYPES LISTED IN THE OBSERVATION ERROR ! FILE. ! ! PROGRAM HISTORY LOG: ! 86-03-21 G. DIMEGO ! 88-11-24 D. DEAVEN RECODED FOR CYBER 205 ! 92-07-29 J. WOOLLEN UPDATED FOR OIQC ! ! USAGE: CALL RDOBER ! ! INPUT ARGUMENTS: ! LU - LOGICAL UNIT FOR OB ERR INPUT ! ! ATTRIBUTES: ! LANGUAGE: CDC FORTRAN200 ! MACHINE: CYBER ! !$$$ 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./ !---------------------------------------------------------------------- !---------------------------------------------------------------------- ! DEFAULT VALUES FOR OB DESCRIPTIONS AND ERRORS ! --------------------------------------------- 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. ! READ OB ERROR VALUES FROM OB ERROR FILE ! --------------------------------------- 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 !$$$ SUBPROGRAM DOCUMENTATION BLOCK ! ! SUBPROGRAM: SATPTS ! PRGMMR: WOOLLEN ORG: NMC22 DATE: 90-11-06 ! ! ABSTRACT: PUTS TOGETHER PACKAGES OF OBSERVATIONS TO BE CHECKED ! IN A GIVEN ITERATION AND INVOKES AUTOTASKING OF THE DATA ! SELECTION, OI, AND DECISION MAKING PROCEDURES. ! ! PROGRAM HISTORY LOG: ! 90-11-06 J. WOOLLEN ! ! USAGE: ! INPUT ARGUMENTS: IN COMMON /CKLIST/ ! NDD - NUMBER OF OBSERVATIONS TO CHECK ! INDD - ARRAY OF OBSERVATION INDEXES TO CHECK ! ! OUTPUT ARGUMENTS: NONE ! ! SUBPROGRAMS CALLED: ! UNIQUE: - SELDAT,QCOI,SATPUT ! ! LIBRARY: ! COMMON - ABORT ! ! REMARKS: ! ! ! ATTRIBUTES: ! LANGUAGE: CRAY FORTRAN ! MACHINE: CRAY ! !$$$ SUBROUTINE SATPTS COMMON /CKLIST/ NDD,INDD(_MAXLEV_) DIMENSION IND( _MAXPAK_ ),NXY( _MAXPAK_ ) DATA NMX / 64 / !---------------------------------------------------------------------- !---------------------------------------------------------------------- ! ASSEMBLE LISTS OF OBS FOR MULTI PROCESSING ! ------------------------------------------ IND(1) = 1 ! DIMENSION IND( _MAXPAK_ ),NXY( _MAXPAK_ ) NXY = 0 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. _MAXPAK_ ) CALL SABORT('QCPTS - IPK > MAXPAK') IPK = IPK+1 IND(IPK) = I NXY(IPK) = KLEV ELSE NXY(IPK) = NXY(IPK) + KLEV ENDIF ENDIF 5 CONTINUE ! MULTI-PROCESS QC OB LISTS ! ------------------------- !MIC$ DOALL AUTOSCOPE DO 10 I=1,IPK CALL SELDAT(IND(I),NXY(I)) CALL QCOI CALL SATPUT 10 CONTINUE RETURN END !$$$ SUBPROGRAM DOCUMENTATION BLOCK ! ! SUBPROGRAM: SATPUT ! PRGMMR: WOOLLEN ORG: NMC22 DATE: 90-11-06 ! ! ABSTRACT: EXAMINES OUTPUTS FROM OI QUALITY CONTROL INTERPOLATION ! AND MAKES DECISION WHETHER A SATELLITE PASSES OR FAILS THE CHECKS. ! DUE TO SPECIFIC CHARACTERISTIC DIFFERENCES OF SATELLITE ! RETRIEVALS FROM CONVENTIONAL 'IN-SITU' OBSERVATIONS, A SPECIAL ! ALGORITHM IS EMPLOYED IN JUDGING THE OUTCOMES OF THE QC CHECKS. ! ALSO PROVIDES FOR A DETAILED PRINT TRACING ALL DECISIONS MADE ! IN REGARD TO SPECIFIED OBSERVATIONS IN THE QC PROCESS. ! ! PROGRAM HISTORY LOG: ! 90-11-06 J. WOOLLEN ! ! USAGE: ! INPUT ARGUMENTS: IN COMMON /HVECT/ ! NXXYY - NUMBER OF OBS IN PACKAGE ! XX - LONGITUDES OF OBS IN PACKAGE ! YY - LATITUDES OF OBS IN PACKAGE ! PP - PRESSURES OF OBS IN PACKAGE ! II - INDEXES OF OBS IN PACKAGE ! LL - REPORT LEVEL OF OBS IN PACKAGE ! FESD - FORECAST ERRORS AT OB LOCATIONS ! AA - INTERPOLATED VALUE OF OBS IN PACKAGE ! AESD - ANALYSIS ERRORS AT OB LOCATIONS ! MKOBS - LIST OF OB INDEXES USED IN EACH INTERPOLATION ! NDIM - NUMBER OF OBS USED IN EACH INTERPOLATION ! JVMAX - NUMBER OF INTERPOLATING OBS FOR EACH CHECK IN ! EACH PACKAGE ! ! OUTPUT FILES: ! FORT.60 - TOSS.OUT - DETAILED TOSSLIST ! FORT.71 - OBPRT1 - OB DETAILED PRINTOUTS FOR QC DIAGNOSIS ! FORT.72 - OBPRT2 - OB DETAILED PRINTOUTS FOR QC DIAGNOSIS ! FORT.73 - OBPRT3 - OB DETAILED PRINTOUTS FOR QC DIAGNOSIS ! FORT.74 - OBPRT4 - OB DETAILED PRINTOUTS FOR QC DIAGNOSIS ! FORT.75 - OBPRT5 - OB DETAILED PRINTOUTS FOR QC DIAGNOSIS ! FORT.76 - OBPRT6 - OB DETAILED PRINTOUTS FOR QC DIAGNOSIS ! FORT.77 - OBPRT7 - OB DETAILED PRINTOUTS FOR QC DIAGNOSIS ! FORT.78 - OBPRT8 - OB DETAILED PRINTOUTS FOR QC DIAGNOSIS ! FORT.79 - OBPRT9 - OB DETAILED PRINTOUTS FOR QC DIAGNOSIS ! FORT.80 - OBPRT10 - OB DETAILED PRINTOUTS FOR QC DIAGNOSIS ! ! OUTPUT ARGUMENTS: IN COMMON /OBMARK/ ! OBIND - CONTAINS QC OUTCOME FOR EACH REPORT AND LEVEL: ! -2 = BELOW MODEL SURFACE (SET IN SUB PSFILE) ! -1 = FAILED OI CHECKS ! 0 = NO DATA IN REPORT ! 1 = PASSED OI CHECKS ! ! SUBPROGRAMS CALLED: ! LIBRARY: ! COMMON - ABORT ! ! EXIT STATES: ! COND = 135 - ERROR DETECTED (SEE CONSOLE OUTPUT) ! ! REMARKS: ! ! ! ATTRIBUTES: ! LANGUAGE: CRAY FORTRAN ! MACHINE: CRAY ! !$$$ SUBROUTINE SATPUT !-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(_MAXREP_),INOB(_MAXREP_),IMAP(360,181) COMMON /CORP / CHLP(1500),CVC,FFLAT(91),IFA( 10 , 10 ) COMMON /SATCOR / EPSCOR(20) !-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'/ ! DATA EPSCOR/ .99, .88,. 61, .36, .18, .05,-.05,-.12,-.15,-.14, ! . -.11,-.08,-.05,-.03,-.02,-.02,-.01,-.01, 0., 0./ LOGICAL LAST DATA PI180/.017453/, MAXNSEG/50/ !---------------------------------------------------------------------- !---------------------------------------------------------------------- LLEV = MCHK(1,1) LREP = FHD(LLEV) !-CRA NSEG = 0 ! DIMENSION NSEG(2),SEG(2) DO I=1,2 NSEG(I) = 0 ENDDO ! LOOP OVER EACH SAT LEVEL IN THIS GROUP ! -------------------------------------- DO 50 I=1,NXXYY LAST = I.EQ.NXXYY ILEV = MCHK(I,1) IREP = FHD(ILEV) ! ITEMIZE COMPLETE PROFILES INTO SEGMENTS ! --------------------------------------- 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) ! STORE THE RESULTS OF THE CHECKS ! ------------------------------- 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 ! A REPORT IS ITEMIZED AND READY TO JUDGE ! --------------------------------------- STAID = SID(LREP) CLAT = FLT(LREP) CLON = FLN(LREP) KX = FKX(LREP) LREP = IREP ! COMPUTE THE TOLERANCE VALUE ! --------------------------- SINLAT = SIN(ABS(CLAT)*PI180) TT = TB*TSAT*SINLAT TMIN = 2.0 TT = MAX(TMIN,TT) ! SUM THE RESULTS OVER THE PROFILE SEGMENTS - MAKE AND STORE DECISION ! ------------------------------------------------------------------- DO 35 K=1,2 IF (NSEG(K).EQ.0) GOTO 35 QC1 = 0. QC2 = 0. !-CRA QCT = 0. !-CRA QCS = 0. !-CRA OMA1 = 0. !-CRA OMA2 = 0. !-CRA EPV1 = 0. !-CRA EPV2 = 0. ! 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 ! WRAP BACK TO START THE NEXT PROFILE ! ----------------------------------- IF(.NOT.LAST) THEN !-CRA NSEG = 0 ! DIMENSION NSEG(2),SEG(2) DO L=1,2 NSEG(L) = 0 ENDDO GOTO 10 ENDIF ! END OF LOOP OVER LEVELS ! ----------------------- 50 CONTINUE RETURN END !----------------------------------------------------------------------- !----------------------------------------------------------------------- SUBROUTINE SEARCH PARAMETER (MAXLV=255) PARAMETER (MTOTH=10 ) PARAMETER (MTOTV=5 ) COMMON /OBS/ NREP,NLEV,SID(_MAXREP_),INOB(_MAXREP_),IMAP(360,181) DIMENSION ITW(MAXLV,MTOTH,2),IPF(MAXLV,MTOTV,2) DATA DUMMY /0/ ! AUTOTASKING LOOP FOR DATA SEARCH ! -------------------------------- !MIC$ DOALL AUTOSCOPE PRIVATE(ITW,IPF,ILV,NLV) !$OMP PARALLEL DO PRIVATE(ITW,IPF,ILV,NLV,K,J,L) DO NPTR=1,NREP CALL SEARCHS(NPTR,ITW,IPF) !$OMP CRITICAL !MIC$ 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 !$OMP END CRITICAL !MIC$ END GUARD 0 ENDDO !$OMP END PARALLEL DO RETURN END !----------------------------------------------------------------------- !----------------------------------------------------------------------- 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(_MAXREP_),INOB(_MAXREP_),IMAP(360,181) COMMON /CORP/ CHLP(1500),CVC,FFLAT(91),IFA( 10 , 10 ) CHARACTER*8 SID CHARACTER*3 TAG,FOC DIMENSION NIND(_MAXREP_),RLAT(_MAXREP_),RLON(_MAXREP_), . RDIS(_MAXREP_),RDIR(_MAXREP_), . DIST(_MAXCYL_),DIRN(_MAXCYL_), . ICOR(_MAXCYL_),IDIR(_MAXCYL_), . LIND(_MAXCYL_),LTYP(_MAXCYL_),PRES(_MAXCYL_), . INL(MAXLV),INI(MAXLV), . SCAT(_MAXCYL_ ,4,4),CAT(_MAXCYL_ ,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 / _MAXCYL_ / , MAXREP / _MAXREP_ / !---------------------------------------------------------------------- !---------------------------------------------------------------------- ! INITIALIZE THE RESULTS FOR THE ZERO CASE ! ---------------------------------------- ! DIMENSION ITW(MAXLV,MTOTH,2),IPF(MAXLV,MTOTV,2) ITW = 0 IPF = 0 IF(FQL(NPTR).LE.0) RETURN ! DO THE DATA SEARCH IF SO INSTRUCTED ! ----------------------------------- BLON = FLN(NPTR) BLAT = FLT(NPTR) BINQ = FQL(NPTR) NBKX = FKX(NPTR) TAG = FOC(NPTR) ! MAKE SPECIAL CASE ADJUSTMENTS (RECCOS,SATEMP) ! --------------------------------------------- IF(MOD(NBKX,100).EQ.32) NC = 3 IF(MOD(NBKX,100).NE.32) NC = 8 IF(TAG.EQ.'SAT') BINQ = BINQ - 1. ! OUTLINE THE SEARCH AREA ! ----------------------- 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. ! ITEMIZE REPORTS IN THE SEACH CYLINDER ! ------------------------------------- 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 IF(NPRF.GT.MAXREP) CALL SABORT('SEARCH - NPRF > MAXREP') NIND(NPRF) = NN RLON(NPRF) = FLN(NN) RLAT(NPRF) = FLT(NN) ENDIF 10 CONTINUE 11 IF(NPRF.GT.MAXREP) CALL SABORT('SEARCH - NPRF > MAXREP') ! COMPUTE DISTANCE AND DIRECTION FROM THE CYLINDER CENTER ! ------------------------------------------------------ CALL CHDIST(BLON,BLAT,RLON,RLAT,RDIS,RDIR,NPRF) ! STORE ALL LEVELS OF EACH REPORT LOCATED IN THE CYLINDER ! ------------------------------------------------------- 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 IF (NOBS.GT.MAXCYL) write(*,*) '1110 nobs=',NOBS, ' maxcyl=',maxcyl IF (NOBS.GT.MAXCYL) write(*,*) '1110 n=',n,1,nprf IF (NOBS.GT.MAXCYL) write(*,*) '1110 jlev=',jlev,ilev,klev IF (NOBS.GT.MAXCYL) CALL SABORT('SEARCH - NOBS>MAXCYL') 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) write(*,*) '1111 nobs=',NOBS, ' maxcyl=',maxcyl IF(NOBS.GT.MAXCYL) CALL SABORT('SEARCH - NOBS>MAXCYL') ENDIF 21 CONTINUE ! LOOP OVER THE OBSERVATION LEVELS TO BE CHECKED ! ---------------------------------------------- ILEV = FLV(NPTR) NLV = FNL(NPTR) DO 50 L=1,NLV IL = ILEV-1+L ! STORE DATA FOR VERTICAL CHECK (UPA ONLY) ! ---------------------------------------- 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 ! USE 3D ZZ CORRELATIONS AND DIRECTION ANGLES TO SELECT HORZ DATA ! --------------------------------------------------------------- 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 ! LOOP OVER DIFFERENT CHECKS FOR THE SAME OB ! ------------------------------------------ DO 45 IQ=1,2 ! 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 ! 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 ! SET THE NONE ARRAY TO FALSE ! --------------------------- DO 28 J=1,4 DO 28 I=1,4 NONE(I,J) = .FALSE. 28 CONTINUE ! SELECTION LOOPS GO UNTIL THEY PICK UP ENOUGH OR ALL POSSIBLE DATA ! ----------------------------------------------------------------- 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 ! STORE THE DATA INDEXES SELECTED FOR HORIZONTAL CHECKS ! ----------------------------------------------------- DO 40 N=1,NTOT ITW(L,N,IQ) = LIND(IPIC(N)) 40 CONTINUE 45 CONTINUE 50 CONTINUE RETURN END !----------------------------------------------------------------------- !----------------------------------------------------------------------- SUBROUTINE SRTPRS(ILEV,NLV,INL,INI) !-CRA DIMENSION INL(NLV),INI(NLV),PRS(NLV),SRT(NLV) DIMENSION INL(NLV),INI(NLV) PARAMETER(NLVPM=1000) DIMENSION PRS(NLVPM),SRT(NLVPM) ! IF(NLV.GT.NLVPM) THEN PRINT *,'NLV.GT.NLVPM. NLV=',NLV CALL SABORT('NLV.GT.NLVPM ') ENDIF !----------------------------------------------------------------------- !----------------------------------------------------------------------- DO L=1,NLV ILV = ILEV+L-1 INI(L) = ILV PRS(L) = FPR(ILV) SRT(L) = PRS(L) ENDDO ! SORT INDEXES IN DESCENDING ORDER OF PRESSURE ! -------------------------------------------- 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 ! MAKE A SET OF LOOK BACK POINTERS ! -------------------------------- 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 !---------------------------------------------------------------------- !---------------------------------------------------------------------- SUBROUTINE SELDAT(IND,NXY) !-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(_MAXLEV_) DIMENSION MAXPIC(4),MAXN(4) DATA MAXPIC/10,10,5,5/ , MAXN/10,10,1,1/ , NMX/ 64 / !---------------------------------------------------------------------- !---------------------------------------------------------------------- ! . MDIM( 64 , 4 ),AA( 64 , 4 ,2), MDIM = 0 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) ! PICK DATA FROM T AND W GROUPS ! ----------------------------- 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 ! PICK DATA FROM SAME PROFILE ! --------------------------- 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 ! CHECK THE NUMBER OF OBS TO STORE ! -------------------------------- IF(M.GT.NMX) CALL SABORT('SELDAT - M GT NMX') ! STORE THE NUMBER OF OBS TO CHECK ! -------------------------------- NXXYY = M RETURN END !$$$ SUBPROGRAM DOCUMENTATION BLOCK ! . . . . ! SUBPROGRAM: SETCON STORES CONSTANTS ! PRGMMR: DEAVEN ORG: W/NMC22 DATE: 88-02-10 ! ! ABSTRACT: INITIALIZES CONSTANTS AND PARAMETERS USING DATA STATEMENTS. ! ALSO READS A FILE OF OBSERVATION ERRORS FROM UNIT 15. ! ! PROGRAM HISTORY LOG: ! 86-03-21 G. DIMEGO ! 88-11-24 D. DEAVEN CODED FOR CYBER 205 ! 92-07-29 J. WOOLLEN UPDATED VERSION FOR OIQC ! ! USAGE: CALL SETCON ! ! ATTRIBUTES: ! LANGUAGE: CDC FORTRAN200 ! MACHINE: CYBER ! !$$$ 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) !-CRA COMMON /OBPRT / NPT,IDSTA(10),PSTA(10),KXSTA(10) COMMON /OBPRT / NPT,NPTADD,IDSTA(10),PSTA(10),KXSTA(10) CHARACTER*8 IDSTA ! 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 / !------------------------------------------------------------------- !------------------------------------------------------------------- ! STORE A SYMMETRIC MATRIX MAP ! ---------------------------- 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 ! MAKE A TABLE OF LENGTH SCALES BY MILIBARS ! ----------------------------------------- DO 20 I=1,1500 CHLP(I) = PILNLNP(FLOAT(I),PMAND,CHL,21) 20 CONTINUE ! READ IN SOME POINTS TO PRINT ! ---------------------------- 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 ! SET UP SOME DATA PRIORITIES ! --------------------------- 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. ! SO LONG FOR NOW ! --------------- RETURN END !$$$ SUBPROGRAM DOCUMENTATION BLOCK ! ! SUBPROGRAM: SSMIPTS ! PRGMMR: WOOLLEN ORG: NMC22 DATE: 90-11-06 ! ! ABSTRACT: PUTS TOGETHER PACKAGES OF SSMI OBSERVATION LOCATIONS, ! PERFORMS A MULTIVARIATE SURFACE WIND ANALYSIS, AND ASSIGNS ! ANALYSED DIRECTION TO THE SSMI WIND SPEED OBSERVATION TO ! PRODUCE AN SSMI WIND VECTOR TO BE CHECKED BY OIQC AND PASSED ! ON THE OBJECTIVE ANALLYSIS PROGRAM. ! ! PROGRAM HISTORY LOG: ! 90-11-06 J. WOOLLEN ! ! USAGE: CALL SSMIPTS ! ! OUTPUT ARGUMENTS: NONE ! ! SUBPROGRAMS CALLED: ! UNIQUE: - SSMISRCH,OIQC,SSMIPUT ! ! LIBRARY: ! COMMON - ABORT ! ! REMARKS: ! ! ! ATTRIBUTES: ! LANGUAGE: CRAY FORTRAN ! MACHINE: CRAY ! !$$$ SUBROUTINE SSMIPTS COMMON /CKLIST/ NDD,INDD(_MAXLEV_) DIMENSION IND( _MAXPAK_ ),NXY( _MAXPAK_ ) DATA NMX / 64 / !---------------------------------------------------------------------- !---------------------------------------------------------------------- ! ASSEMBLE LISTS OF OBS FOR MULTI PROCESSING ! ------------------------------------------ IND(1) = 1 IPK = 1 ! DIMENSION IND( _MAXPAK_ ),NXY( _MAXPAK_ ) NXY = 0 DO 5 I=1,NDD ILEV = INDD(I) IF(NXY(IPK)+1.GT.NMX) THEN IF(IPK+1.GT. _MAXPAK_ ) CALL SABORT('SSMIPTS - IPK > MAXPAK') IPK = IPK+1 IND(IPK) = I NXY(IPK) = 1 ELSE NXY(IPK) = NXY(IPK)+1 ENDIF 5 CONTINUE ! MULTI-PROCESS QC OB LISTS ! ------------------------- !MIC$ DOALL AUTOSCOPE DO 10 I=1,IPK CALL SSMISRCH(IND(I),NXY(I)) CALL QCOI !MIC$ GUARD CALL SSMIPUT !MIC$ ENDGUARD 10 CONTINUE RETURN END !$$$ SUBPROGRAM DOCUMENTATION BLOCK ! ! SUBPROGRAM: SSMIPUT ! PRGMMR: WOOLLEN ORG: NMC22 DATE: 92-10-06 ! ! ABSTRACT: USES THE OUTPUT OF A MULTIVARIATE INTERPOLATION OF ! NEARBY OBSERVATIONS TO ASSIGN DIRECTIONS TO SSMI WIND SPEED DATA. ! ALSO PROVIDES FOR A DETAILED PRINT SHOWING OBSERVATIONS USED IN ! THE INTERPOLATION ALONG WITH THE WEIGHTS AND INCREMENTS RESULTING. ! ! PROGRAM HISTORY LOG: ! 90-11-06 J. WOOLLEN ! ! USAGE: CALL SSMIPUT ! ! OUTPUT FILES: ! FORT.71 - FORT.80 - DETAILED PRINTOUTS FOR ANALYSIS ! ! SUBPROGRAMS CALLED: ! LIBRARY: ! COMMON - ABORT ! ! REMARKS: ! ! ! ATTRIBUTES: ! LANGUAGE: CRAY FORTRAN ! MACHINE: CRAY ! !$$$ SUBROUTINE SSMIPUT !-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(_MAXREP_),INOB(_MAXREP_),IMAP(360,181) !-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/ !---------------------------------------------------------------------- !---------------------------------------------------------------------- ! LOOP OVER EACH OB CHECKED IN THIS GROUP ! --------------------------------------- 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) ! COMPUTE THE ANALYSED WIND DIRECTION AT THIS OB ! ---------------------------------------------- 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 ! COMBINE THE ANA DIRECTION WITH THE SSMI SPEED ! --------------------------------------------- SSMISPD = ABS(FOB(ILEV,2)) SSMIU = SSMISPD*COS(ANADIR) SSMIV = SSMISPD*SIN(ANADIR) ! COMPUTE THE OBSERVED INCREMENT AND STORE IN THE PROPER PLACE ! ------------------------------------------------------------ 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,')' ), 100 FORMAT(A6,I4,3F8.2,A4,2F8.2,4(' *',F6.1,' (',F5.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 !----------------------------------------------------------------------- ! !----------------------------------------------------------------------- SUBROUTINE SSMISRCH(IND,NXY) COMMON /CKLIST/ NDD,INDD(_MAXLEV_) COMMON /CORP / CHLP(1500),CVC,FFLAT(91),IFA( 10 , 10 ) COMMON /OBS/ NREP,NLEV,SID(_MAXREP_),INOB(_MAXREP_),IMAP(360,181) !-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(_MAXREP_),RLAT(_MAXREP_),RLON(_MAXREP_), . DIST(_MAXCYL_),DIRN(_MAXCYL_), . LIND(_MAXCYL_),LTYP(_MAXCYL_),PRES(_MAXCYL_), . ICOR(_MAXCYL_),IDIR(_MAXCYL_), . IPIC( 10 ),CATCOR(101) LOGICAL CAT(_MAXCYL_ ,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 / _MAXCYL_ / , MAXREP / _MAXREP_ / DATA MTOTH /10/ , MTOTV /5/ !---------------------------------------------------------------------- !---------------------------------------------------------------------- NXXYY = 0 ! . MDIM( 64 , 4 ),AA( 64 , 4 ,2), MDIM = 0 ! LOOP OVER PROFILES TO BE CHECKED ! -------------------------------- 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') ! OUTLINE THE SEARCH AREA ! ----------------------- 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. ! ITEMIZE REPORTS IN THE SEACH CYLINDER ! ------------------------------------- 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') ! COMPUTE DISTANCE AND DIRECTION FROM THE CYLINDER CENTER ! ------------------------------------------------------- CALL CHDIST(BLON,BLAT,RLON,RLAT,DIST,DIRN,NPRF) ! STORE ALL LEVELS OF EACH REPORT LOCATED IN THE CYLINDER ! ------------------------------------------------------- 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 IF (NOBS.GT.MAXCYL) write(*,*) '1112: nobs=',nobs IF (NOBS.GT.MAXCYL) write(*,*) 'n=',n,1,nprf IF (NOBS.GT.MAXCYL) write(*,*) 'jlev=',jlev,ilev,klev IF (NOBS.GT.MAXCYL) CALL SABORT('SSMISRCH - NOBS>MAXCYL') 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 ! LOOP OVER THE OBSERVATION LEVELS TO BE CHECKED ! ---------------------------------------------- ILEV = FLV(NPTR) KLEV = ILEV+FNL(NPTR)-1 DO 50 IL=ILEV,KLEV BPRS = FPR(IL) IPRS = BPRS BPRS = LOG(BPRS) ! INITIALIZE THE SLECTION ARRAYS ! ------------------------------ 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 ! COMPUTE ZZ HRZ AND VRT CORRELATIONS FOR EACH OB AND CATEGORIZE ! -------------------------------------------------------------- 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 ! PICK THE CLOSEST 10 (MULTIVARIATE) HORIZONTAL OBSERVATIONS ! ---------------------------------------------------------- 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 ! STORE THE DATA INDEXES SELECTED FOR HORIZONTAL CHECKS ! ----------------------------------------------------- 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 !$$$ SUBPROGRAM DOCUMENTATION BLOCK ! ! SUBPROGRAM: STOERR ! PRGMMR: WOOLLEN ORG: NMC22 DATE: 90-11-06 ! ! ABSTRACT: INTERPOLATES TO AND STORES FORECAST AND OBSERVATION ERRORS ! WITH EACH OBSERVATION IN THE DATASET IN COMMON /DUMMY/. ALSO ! MAKES A LIST OF OBS BY 2D LOCATION IN COMMON /TGRID/ IN ORDER ! TO EXPEDITE VARIOUS DATA SEARCH PROCEDURES. SOME CHECKING ON THE ! PLAUSIBILITY OF OB COORDINATES IS PERFORMED AS WELL. ! ! PROGRAM HISTORY LOG: ! 90-11-06 J. WOOLLEN ! ! USAGE: ! INPUT ARGUMENTS: NONE ! OUTPUT ARGUMENTS IN COMMON /TGRID/ ! NREP - NUMBER OF REPORT LOCATIONS ! XQC - ARRAY OF REPORT LONGITUDES ! YQC - ARRAY OF REPORT LATITUDES ! IQC - ARRAY OF REPORT INDEXES INTO THE BUFFS ARRAY ! TQC - ARRAY OF REPORT TYPES (HEIGHT OR WIND) ! FERR - ARRAY OF FORECAST ERRORS AT REPORT LOCATIONS ! ! SUBPROGRAMS CALLED: ! LIBRARY: ! COMMON - ABORT ! ! EXIT STATES: ! COND = 11 - PROBLEM DETECTED (SEE CONSOLE OUTPUT) ! ! REMARKS: ! ! ! ATTRIBUTES: ! LANGUAGE: CRAY FORTRAN ! MACHINE: CRAY ! !$$$ !----------------------------------------------------------------------- SUBROUTINE STOERR COMMON /OBS/ NREP,NLEV,SID(_MAXREP_),INOB(_MAXREP_),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./ !---------------------------------------------------------------------- !----------------------------------------------------------------------- ! GO THROUGH ALL THE STORED DATA ! ------------------------------ 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 ! CHECK FOR UNNATURAL INFORMATION IN THE HEADER ! --------------------------------------------- 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 ! LOOP THRU THE REPORT LEVELS ! --------------------------- DO 20 KL=1,KLEV L = ILEV+KL-1 ! CHECK TO SEE IF A SURFACE OBSERVATION IS PRESENT ! ------------------------------------------------ IF(FOE(L).EQ.0.) THEN GOTO 20 ENDIF ! GET FORECAST ERRORS FOR THIS REPORT ! ----------------------------------- 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) ! FIND THE OB ERROR FOR THIS REPORT ! --------------------------------- 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 ! STORE THE ERRORS IN THEIR RIGHTFUL PLACE ! ---------------------------------------- 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) ! SPECIAL TREATMENT FOR SPECIAL CASES ! ----------------------------------- 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 ! STORE AN OIQC OBSERVATION ERROR MAYBE ! ------------------------------------- IF(USEOE) PACK = WOE(L,OBER) ! END OF THE LOOPS ! ---------------- 20 CONTINUE 30 CONTINUE PRINT*,' IN STORERR - NREP=',NREP RETURN END !$$$ SUBPROGRAM DOCUMENTATION BLOCK ! ! SUBPROGRAM: TOSSLIST ! PRGMMR: WOOLLEN ORG: NMC22 DATE: 90-11-06 ! ! ABSTRACT: WRITES A LIST OF REJECTED OBSERVATIONS PLUS A SUMMARY ! ! PROGRAM HISTORY LOG: ! 90-11-06 J. WOOLLEN ! ! USAGE: ! INPUT ARGUMENTS: NONE ! LO - LOGICAL UNIT FOR TOSSLIST OUTPUT ! ! OUTPUT ARGUMENTS: NONE ! ! OUTPUT FILES: ! FORT.LO - TOSSLIST - OFFICIAL OPERATIONAL TOSSLIST ! ! SUBPROGRAMS CALLED: ! LIBRARY: ! COMMON - ABORT ! ! REMARKS: ! ! ! ATTRIBUTES: ! LANGUAGE: CRAY FORTRAN ! MACHINE: CRAY ! !$$$ !----------------------------------------------------------------------- SUBROUTINE TOSSLIST(LO) COMMON /OBS/ NREP,NLEV,SID(_MAXREP_),INOB(_MAXREP_),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./ !---------------------------------------------------------------------- !----------------------------------------------------------------------- ! START UP ! -------- !RA IBAD = 0 !RA ITOT = 0 !RA ISFZ = 0 !RA PZ = 0 !RA PU = 0 !RA PT = 0 ! DIMENSION IBAD(30,3),ITOT(30),ISFZ(30),PZ(30),PU(30),PT(30) IBAD = 0 DO I=1,30 ITOT(I) = 0 ISFZ(I) = 0 PZ (I) = 0 PU (I) = 0 PT (I) = 0 ENDDO ! COUNT BAD DATA IN THE BIG BUFFER - WRITE THE TOSSLIST ! ----------------------------------------------------- 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 ! COMPUTE GROSS PERCENTAGE OF TOSSES ! ---------------------------------- 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 ! WRITE A SHORT SUMMARY OF TOSSES ! ------------------------------- 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) ! FORMAT STATEMENTS FOR THE TOSS LIST ! ----------------------------------- 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 !CBREAKCC !$$$ SUBPROGRAM DOCUMENTATION BLOCK ! ! SUBPROGRAM: VSOLVE ! PRGMMR: WOOLLEN ORG: NMC22 DATE: 90-11-06 ! ! ABSTRACT: CHOLESKY TYPE SOLUTION FOR ARRAYS OF POSITIVE DEFINITE ! SYMMETRIC MATRIXES. ! ! PROGRAM HISTORY LOG: ! 90-11-06 J. WOOLLEN ! ! USAGE: ! INPUT ARGUMENTS: ! A - ARRAY OF SYMMETRIC MATRIXES ! B - ARRAY OF RIGHT HAND SIDE VECTORS ! NDIM - ARRAY OF MATRIX RANKS ! BAD - BAD MATRIX INDICATOR ! NFT - NUMBER OF RIGHT HAND SIDES PER MATRIX ! NS - NUMBER OF MATRIXES TO SOLVE ! MAXDIM - LARGEST RANK MATRIX IN STACK ! ! OUTPUT ARGUMENTS: ! B - CONTAINS THE SOLUTIONS UPON RETURN ! ! SUBPROGRAMS CALLED: NONE ! ! EXIT STATES: NONE ! ! REMARKS: TRIANGULAR REPRESENTATION LISTS ELEMENTS TOWARDS THE ! DIAGONAL. ALGORITHM FROM H.CARUS. ! ! ATTRIBUTES: ! LANGUAGE: CRAY FORTRAN ! MACHINE: CRAY ! !$$$ !----------------------------------------------------------------------- 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/ !---------------------------------------------------------------------- IX (I,J) = I*(I-1)/2 + J !---------------------------------------------------------------------- DO 1 M=1,NS 1 BAD(M) = .FALSE. ! DECOMPOSE THE MATRIXES ! ---------------------- 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 ! SOLVE FOR ALL RIGHT HAND SIDES ! ------------------------------ DO 40 NF=1,NFT ! FORWARD SUBSTITUTION ! -------------------- 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)) ! BACKWARD SUBSTITUTION ! --------------------- 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 !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ! ! .=================================================================== ! | DESCRIPTION | UNPK | PACK | BITS | * | DIMENSIONS ! |=============*======*======*======*===*============================ ! | LEVEL LINK | FLV | WLV | WORD | I | ?MAXREP ! |-------------+------+------+------+---+---------------------------- ! | OB CLASS | FOC | WOC | C3 | C | ?MAXREP ! |-------------+------+------+------+---+---------------------------- ! | LATITUDE | FLT | WLT | WORD | R | ?MAXREP ! |-------------+------+------+------+---+---------------------------- ! | LONGITUDE | FLN | WLN | WORD | R | ?MAXREP ! |-------------+------+------+------+---+---------------------------- ! | OB TIME | FTM | WTM | WORD | R | ?MAXREP ! |-------------+------+------+------+---+---------------------------- ! | OB TYPE | FKX | WKX | WORD | I | ?MAXREP ! |-------------+------+------+------+---+---------------------------- ! | N OF LEVELS | FNL | WNL | WORD | I | ?MAXREP ! |-------------+------+------+------+---+---------------------------- ! | USE REP FLG | FQL | WQL | WORD | R | ?MAXREP ! |-------------+------+------+------+---+---------------------------- ! | PREPDA REC | FPRC | WPRC | WORD | I | ?MAXREP ! |-------------+------+------+------+---+---------------------------- ! | PREPDA FLD | FPFD | WPFD | WORD | I | ?MAXREP ! |-------------+------+------+------+---+---------------------------- ! | PREPDA LEV | FPLV | WPLV | WORD | I | ?MAXLEV ! |-------------+------+------+------+---+---------------------------- ! | GES SFC U,V | FMW | WMW | WORD | R | ?MAXREP,2 ! |-------------+------+------+------+---+---------------------------- ! | HEADER LINK | FHD | WHD | WORD | I | ?MAXLEV ! |-------------+------+------+------+---+---------------------------- ! | PRESSURE | FPR | WPR | WORD | R | ?MAXLEV ! |-------------+------+------+------+---+---------------------------- ! | PREDA FLAG | FQM | WQM | WORD | I | ?MAXLEV ! |-------------+------+------+------+---+---------------------------- ! | OIQC FLAG | FOI | WOI | WORD | I | ?MAXLEV ! |-------------+------+------+------+---+---------------------------- ! | OB INCS | FOB | WOB | WORD | R | ?MAXLEV,2 ! |-------------+------+------+------+---+---------------------------- ! | OB ERROR | FOE | WOE | WORD | R | ?MAXLEV ! |-------------+------+------+------+---+---------------------------- ! | T/W FC ERR | FFE | WFE | WORD | R | ?MAXLEV,2 ! |-------------+------+------+------+---+---------------------------- ! | T/Z FACTOR | FTZ | WTZ | WORD | R | ?MAXLEV ! |-------------+------+------+------+---+---------------------------- ! | T/W OBS | FTW | WTW | WORD | I | ?MAXLEV,10,2 ! |-------------+------+------+------+---+---------------------------- ! | PROF A/B | FPF | WPF | WORD | I | ?MAXLEV,5,2 ! `=============^======^======^======^===^============================ ! !----------------------------------------------------------------------- !----------------------------------------------------------------------- FUNCTION FLV(I) PARAMETER (IM = _MAXREP_) COMMON /STOR1/IAR(IM) INTEGER IAR FLV = IAR(I) RETURN END !----------------------------------------------------------------------- FUNCTION WLV(I,V) PARAMETER (IM = _MAXREP_) COMMON /STOR1/IAR(IM) INTEGER IAR,V IAR(I) = V WLV = 0 RETURN END !----------------------------------------------------------------------- FUNCTION FOC(I) PARAMETER (IM = _MAXREP_) COMMON /STOR2/IAR(IM) CHARACTER*3 IAR,FOC FOC = IAR(I) RETURN END !----------------------------------------------------------------------- FUNCTION WOC(I,V) PARAMETER (IM = _MAXREP_) COMMON /STOR2/IAR(IM) CHARACTER*3 IAR,V IAR(I) = V WOC = 0 RETURN END !----------------------------------------------------------------------- FUNCTION FLT(I) PARAMETER (IM = _MAXREP_) COMMON /STOR3/IAR(IM) REAL IAR FLT = IAR(I) RETURN END !----------------------------------------------------------------------- FUNCTION WLT(I,V) PARAMETER (IM = _MAXREP_) COMMON /STOR3/IAR(IM) REAL IAR,V IAR(I) = V WLT = 0 RETURN END !----------------------------------------------------------------------- FUNCTION FLN(I) PARAMETER (IM = _MAXREP_) COMMON /STOR4/IAR(IM) REAL IAR FLN = IAR(I) RETURN END !----------------------------------------------------------------------- FUNCTION WLN(I,V) PARAMETER (IM = _MAXREP_) COMMON /STOR4/IAR(IM) REAL IAR,V IAR(I) = V WLN = 0 RETURN END !----------------------------------------------------------------------- FUNCTION FTM(I) PARAMETER (IM = _MAXREP_) COMMON /STOR5/IAR(IM) REAL IAR FTM = IAR(I) RETURN END !----------------------------------------------------------------------- FUNCTION WTM(I,V) PARAMETER (IM = _MAXREP_) COMMON /STOR5/IAR(IM) REAL IAR,V IAR(I) = V WTM = 0 RETURN END !----------------------------------------------------------------------- FUNCTION FKX(I) PARAMETER (IM = _MAXREP_) COMMON /STOR6/IAR(IM) INTEGER IAR FKX = IAR(I) RETURN END !----------------------------------------------------------------------- FUNCTION WKX(I,V) PARAMETER (IM = _MAXREP_) COMMON /STOR6/IAR(IM) INTEGER IAR,V IAR(I) = V WKX = 0 RETURN END !----------------------------------------------------------------------- FUNCTION FNL(I) PARAMETER (IM = _MAXREP_) COMMON /STOR7/IAR(IM) INTEGER IAR FNL = IAR(I) RETURN END !----------------------------------------------------------------------- FUNCTION WNL(I,V) PARAMETER (IM = _MAXREP_) COMMON /STOR7/IAR(IM) INTEGER IAR,V IAR(I) = V WNL = 0 RETURN END !----------------------------------------------------------------------- FUNCTION FQL(I) PARAMETER (IM = _MAXREP_) COMMON /STOR8/IAR(IM) REAL IAR FQL = IAR(I) RETURN END !----------------------------------------------------------------------- FUNCTION WQL(I,V) PARAMETER (IM = _MAXREP_) COMMON /STOR8/IAR(IM) REAL IAR,V IAR(I) = V WQL = 0 RETURN END !----------------------------------------------------------------------- FUNCTION FPRC(I) PARAMETER (IM = _MAXREP_) COMMON /STOR9/IAR(IM) INTEGER IAR FPRC = IAR(I) RETURN END !----------------------------------------------------------------------- FUNCTION WPRC(I,V) PARAMETER (IM = _MAXREP_) COMMON /STOR9/IAR(IM) INTEGER IAR,V IAR(I) = V WPRC = 0 RETURN END !----------------------------------------------------------------------- FUNCTION FPFD(I) PARAMETER (IM = _MAXREP_) COMMON /STOR10/IAR(IM) INTEGER IAR FPFD = IAR(I) RETURN END !----------------------------------------------------------------------- FUNCTION WPFD(I,V) PARAMETER (IM = _MAXREP_) COMMON /STOR10/IAR(IM) INTEGER IAR,V IAR(I) = V WPFD = 0 RETURN END !----------------------------------------------------------------------- FUNCTION FPLV(I) PARAMETER (IM = _MAXLEV_) COMMON /STOR11/IAR(IM) INTEGER IAR FPLV = IAR(I) RETURN END !----------------------------------------------------------------------- FUNCTION WPLV(I,V) PARAMETER (IM = _MAXLEV_) COMMON /STOR11/IAR(IM) INTEGER IAR,V IAR(I) = V WPLV = 0 RETURN END !----------------------------------------------------------------------- FUNCTION FMW(I,J) PARAMETER (IM = _MAXREP_) PARAMETER (JM = 2) COMMON /STOR12/IAR(IM,JM) REAL IAR FMW = IAR(I,J) RETURN END !----------------------------------------------------------------------- FUNCTION WMW(I,J,V) PARAMETER (IM = _MAXREP_) PARAMETER (JM = 2) COMMON /STOR12/IAR(IM,JM) REAL IAR,V IAR(I,J) = V WMW = 0 RETURN END !----------------------------------------------------------------------- FUNCTION FHD(I) PARAMETER (IM = _MAXLEV_) COMMON /STOR13/IAR(IM) INTEGER IAR FHD = IAR(I) RETURN END !----------------------------------------------------------------------- FUNCTION WHD(I,V) PARAMETER (IM = _MAXLEV_) COMMON /STOR13/IAR(IM) INTEGER IAR,V IAR(I) = V WHD = 0 RETURN END !----------------------------------------------------------------------- FUNCTION FPR(I) PARAMETER (IM = _MAXLEV_) COMMON /STOR14/IAR(IM) REAL IAR FPR = IAR(I) RETURN END !----------------------------------------------------------------------- FUNCTION WPR(I,V) PARAMETER (IM = _MAXLEV_) COMMON /STOR14/IAR(IM) REAL IAR,V IAR(I) = V WPR = 0 RETURN END !----------------------------------------------------------------------- FUNCTION FQM(I) PARAMETER (IM = _MAXLEV_) COMMON /STOR15/IAR(IM) INTEGER IAR FQM = IAR(I) RETURN END !----------------------------------------------------------------------- FUNCTION WQM(I,V) PARAMETER (IM = _MAXLEV_) COMMON /STOR15/IAR(IM) INTEGER IAR,V IAR(I) = V WQM = 0 RETURN END !----------------------------------------------------------------------- FUNCTION FOI(I) PARAMETER (IM = _MAXLEV_) COMMON /STOR16/IAR(IM) INTEGER IAR FOI = IAR(I) RETURN END !----------------------------------------------------------------------- FUNCTION WOI(I,V) PARAMETER (IM = _MAXLEV_) COMMON /STOR16/IAR(IM) INTEGER IAR,V IAR(I) = V WOI = 0 RETURN END !----------------------------------------------------------------------- FUNCTION FOB(I,J) PARAMETER (IM = _MAXLEV_) PARAMETER (JM = 2) COMMON /STOR17/IAR(IM,JM) REAL IAR FOB = IAR(I,J) RETURN END !----------------------------------------------------------------------- FUNCTION WOB(I,J,V) PARAMETER (IM = _MAXLEV_) PARAMETER (JM = 2) COMMON /STOR17/IAR(IM,JM) REAL IAR,V IAR(I,J) = V WOB = 0 RETURN END !----------------------------------------------------------------------- FUNCTION FOE(I) PARAMETER (IM = _MAXLEV_) COMMON /STOR18/IAR(IM) REAL IAR FOE = IAR(I) RETURN END !----------------------------------------------------------------------- FUNCTION WOE(I,V) PARAMETER (IM = _MAXLEV_) COMMON /STOR18/IAR(IM) REAL IAR,V IAR(I) = V WOE = 0 RETURN END !----------------------------------------------------------------------- FUNCTION FFE(I,J) PARAMETER (IM = _MAXLEV_) PARAMETER (JM = 2) COMMON /STOR19/IAR(IM,JM) REAL IAR FFE = IAR(I,J) RETURN END !----------------------------------------------------------------------- FUNCTION WFE(I,J,V) PARAMETER (IM = _MAXLEV_) PARAMETER (JM = 2) COMMON /STOR19/IAR(IM,JM) REAL IAR,V IAR(I,J) = V WFE = 0 RETURN END !----------------------------------------------------------------------- FUNCTION FTZ(I) PARAMETER (IM = _MAXLEV_) COMMON /STOR20/IAR(IM) REAL IAR FTZ = IAR(I) RETURN END !----------------------------------------------------------------------- FUNCTION WTZ(I,V) PARAMETER (IM = _MAXLEV_) COMMON /STOR20/IAR(IM) REAL IAR,V IAR(I) = V WTZ = 0 RETURN END !----------------------------------------------------------------------- FUNCTION FTW(I,J,K) PARAMETER (IM = _MAXLEV_) PARAMETER (JM = 10) PARAMETER (KM = 2) COMMON /STOR21/IAR(IM,JM,KM) INTEGER IAR FTW = IAR(I,J,K) RETURN END !----------------------------------------------------------------------- FUNCTION WTW(I,J,K,V) PARAMETER (IM = _MAXLEV_) PARAMETER (JM = 10) PARAMETER (KM = 2) COMMON /STOR21/IAR(IM,JM,KM) INTEGER IAR,V IAR(I,J,K) = V WTW = 0 RETURN END !----------------------------------------------------------------------- FUNCTION FPF(I,J,K) PARAMETER (IM = _MAXLEV_) PARAMETER (JM = 5) PARAMETER (KM = 2) COMMON /STOR22/IAR(IM,JM,KM) INTEGER IAR FPF = IAR(I,J,K) RETURN END !----------------------------------------------------------------------- FUNCTION WPF(I,J,K,V) PARAMETER (IM = _MAXLEV_) PARAMETER (JM = 5) PARAMETER (KM = 2) COMMON /STOR22/IAR(IM,JM,KM) INTEGER IAR,V IAR(I,J,K) = V WPF = 0 RETURN END !$$$ SUBPROGRAM DOCUMENTATION BLOCK ! . . . . ! SUBPROGRAM: STORE STORES OBSERVED DATA IN BUFFER ! PRGMMR: DEAVEN ORG: W/NMC22 DATE: 88-02-10 ! ! ABSTRACT: STORES ALL DATA AND CORRESPONDING INDEX FOR LATER USE ! IN DATA COLLECT AND SEARCH ROUTINES. ! ! PROGRAM HISTORY LOG: ! 86-03-21 G. DIMEGO ! 88-11-24 D. DEAVEN CODED FOR CYBER 205 ! 90-11-06 J. WOOLLEN MODIFIED COMMON /BUFFS/ FOR QCOI ! 92-07-28 J. WOOLLEN MODIFIED FOR STORSTAR VERSION ! ! USAGE: CALL STORE(LU) ! ! INPUT ARGUMENTS: ! LU - LOGICAL UNIT FOR INCREMENT FILE ! ! ATTRIBUTES: ! LANGUAGE: CDC FORTRAN200 ! MACHINE: CYBER ! !$$$ SUBROUTINE STORE(LUDAT,LUBFI) COMMON /OBS/ NREP,NLEV,SID(_MAXREP_),INOB(_MAXREP_),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 / _MAXREP_ / DATA MAXLEV / _MAXLEV_ / DATA ROG / 29.261 / DATA QX / 300*.TRUE./ !---------------------------------------------------------------------- PA(POB) = POB*(1.-.053) PB(POB) = POB*(1.+.053) !---------------------------------------------------------------------- ! INITIALIZE ! ---------- STUCK = .FALSE. NREC = 0 NREP = 0 NLEV = 0 MAXNLEV = 0 ! OPEN THE DATA FILE AND CHECK THE DATE ! ------------------------------------- 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* ! LOOP THROUGH THE INPUT MESSAGES - READ THE NEXT SUBSET ! ------------------------------------------------------ !MIC$ DO ALL AUTOSCOPE !MIC$.SHARED (JDATE,NREC,NREP,NLEV,MAXNLEV,NREP0,NLEV0,QUALOB) !MIC$.SHARED (LUBFI,HDSTR,OBSTR,FCSTR,QMSTR,STUCK,SID,ROG) !MIC$.PRIVATE (SUBSET,IDATE,IREC,IFLD,HDR,OBS,FCS,QMS) !MIC$.PRIVATE (PCHK,ILEV,KLEV,KX,TAG,PACK,LEV) !MIC$.PRIVATE (POB,TOB,ZOB,UOB,VOB) !MIC$.PRIVATE (PFC,TFC,ZFC,UFC,VFC) !MIC$.PRIVATE (IQM,OB1,OB2,IOI,OER) !MIC$.PRIVATE (PQM,TQM,ZQM,WQQ) DO 50 KREC=1,MSGS !MIC$ 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 !MIC$ 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) ! DECIDE WHETHER OIQC SHOULD PROCESS - SET A TAG ! ---------------------------------------------- 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 ! SAVE THE CURRENT COUNTERS IN CASE OF DISASTER OR ARRAY OVERFLOW ! --------------------------------------------------------------- !MIC$ GUARD 1 IF(.NOT.STUCK) THEN NREP0 = NREP NLEV0 = NLEV ENDIF ! STORE THE HEADER ! ---------------- 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 ! STORE THE OB LEVELS ! ------------------- 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 ! KEEP TRACK OF MAX LEVELS ! ------------------------ IF(ILEV.EQ.0) NREP = NREP-1 MAXNLEV = MAX(ILEV,MAXNLEV) !MIC$ ENDGUARD 1 ! END OF READ LOOPS AND PARALLEL SECTION ! -------------------------------------- ENDDO 50 ENDDO CALL CLOSBF(LUBFI) ! ALL DONE STORING - RECORD THE OUTCOME AND EXIT ! ---------------------------------------------- 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 !$$$ SUBPROGRAM DOCUMENTATION BLOCK ! ! SUBPROGRAM: PREPQM ! PRGMMR: WOOLLEN ORG: NMC22 DATE: 90-11-06 ! ! ABSTRACT: REWRITES THE BUFRDA FILE WITH UPDATED QUALITY MARKS ! ! PROGRAM HISTORY LOG: ! 90-11-06 J. WOOLLEN ! ! USAGE: ! INPUT ARGUMENTS: NONE ! LPIN - LOGICAL UNIT FOR PREPDA INPUT ! LPOUT - LOGICAL UNIT FOR PREPQM OUTPUT ! ! INPUT FILES: ! FORT.LPIN - PREPQM - PREPDA WITH UPDATED QUALITY MARKS ! ! OUTPUT ARGUMENTS: NONE ! ! OUTPUT FILES: ! FORT.LPOUT - PREPQM - PREPDA WITH UPDATED QUALITY MARKS ! ! SUBPROGRAMS CALLED: ! LIBRARY: ! COMMON - ABORT ! ! REMARKS: ! ! ! ATTRIBUTES: ! LANGUAGE: CRAY FORTRAN ! MACHINE: CRAY ! !$$$ SUBROUTINE PREPQM(LUBIN,LUBOT) PARAMETER (MAXM= _MAXM_ , MAXS= _MAXS_ ) COMMON /OBS/ NREP,NLEV,SID(_MAXREP_),INOB(_MAXREP_),IMAP(360,181) ! 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/ !----------------------------------------------------------------------- !----------------------------------------------------------------------- PRINT*,'******WRITE THE BUFRQM FILE*******' NREC = 0 ! DIMENSION OBS(10,255),EVN(10,255),JREP(MAXM,MAXS) JREP = 0 ! SORT OUT THE FILE POSITIONS OF EACH OB IN STORAGE ! ------------------------------------------------- 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 ! INITIALIZE THE INPUT AND OUTPUT FILES ! ------------------------------------- CALL UFBREW(LUBIN,MSGS) CALL OPENBF(LUBIN,'IN',LUBIN) CALL OPENBF(LUBOT,'OUT',LUBIN) CALL UFBQCD(LUBIN,'OIQC',QCD) ! READ MESSAGES, SUBSETS, DATA, EVENTS - TRY TO KEEP IT ALL STRAIGHT ! ------------------------------------------------------------------ !MIC$ DO ALL AUTOSCOPE !MIC$.SHARED (LUBIN,LUBOT,NREC,QCD,ESTR,OSTR,JREP,BMISS) !MIC$.PRIVATE (SUBSET,IDATE,IREC,IFLD,OBS,EVN,UOB,VOB) !MIC$.PRIVATE (IREP,ILEV,JLEV,KLEV,JUPD,LEV,LEO) !MIC$.PRIVATE (TAG,KKX,SFC,SMI,QMA,QMB,OBM,L) !MIC$.PRIVATE (sid1,sid2) DO KREC=1,MSGS !MIC$ 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 !MIC$ ENDGUARD 0 DO WHILE(IREADSB(LUBIN).EQ.0) CALL UFBCPY(LUBIN,LUBOT) IFLD = IFLD+1 ! FIND THIS OB IN STORAGE OR JUST WRITE IT BACK OUT ! ------------------------------------------------- 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 ! CHECK ALL THE LEVELS OF THIS OB FOR EVENTS ! ------------------------------------------ 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 ! DIMENSION OBS(10,255),EVN(10,255),JREP(MAXM,MAXS) EVN = BMISS 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 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)) 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 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 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 CALL SABORT('PREPQM - JUPD OUT OF RANGE') ENDIF ENDIF ENDDO ENDIF CALL WRITSB(LUBOT) ENDDO CALL CLOSMG(LUBOT) ENDDO ! HERE AT THE END ! --------------- 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