! 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 _MAXCYL_ 800000
! #define _MAXREP_ 800000
#define _MAXREP_ 1200000
! 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_
! 15-03-17 W. EBISUZAKI ADDED OpenMP SEARCH LOOP
!
! 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
write(*,*) 'LIMITS MAXLEV(NN):', nn, _MAXLEV_
! 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
write(*,*) 'LIMITS MAXREP (NOB):', NOB, _MAXREP_
! 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
! write(*,*) 'LIMITS MAXPAK(IPK,2):', IPK, _MAXPAK_
! 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
! write(*,*) 'LIMITS MAXPAK(IPK,QCPTS)', ipk , _MAXPAK_
! 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
! write(*,*) 'LIMITS MAXPAK(IPK): ', IPK, _MAXPAK_
! 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