C$$$ MAIN PROGRAM DOCUMENTATION BLOCK C C MAIN PROGRAM: PREPOBS_OIQCBUFR C PRGMMR: WOOLLEN/KEYSER ORG: NP20 DATE: 2013-06-19 C C ABSTRACT: QUALITY CONTROLS OBSERVATIONS TO BE USED IN THE GLOBAL C ANALYSIS-FORECAST SYSTEM. OBSERVATIONAL DATASET GENERATED BY C THE GESTEMP PREPROCESSOR IS INPUT TO OIQC WHICH THEN PERFORMS C A "COMPLEX" QUALITY CONTROL (FROM GANDIN,1985) EMPLOYING C MULTIVARIATE OI TO CHECK ALL OBSERVATIONS AGAINST NEARBY C NEIGHBORS. SEVERAL ITERATIONS OF THE BASIC SCHEME ARE REQUIRED C TO COMPLETE THE PROCESS. DURING EACH ITERATION ALL OBSERVATIONS C ARE SUBJECTED TO FIVE INTERPOLATION CHECKS. CHECK #1 INTERPOLATES C COMPARISON VALUES FROM NEARBY TEMPERATURE DATA. CHECK #2 AND CHECK C #3 INTERPOLATE COMPARISON VALUES FROM NEARBY ZONAL AND LONGITUDINAL C WIND COMPONENTS RESPECTIVELY. CHECK #4 IS AN INTERPOLATION C FROM THE NEAREST OBS IN THE SAME PROFILE AS THE ONE BEING CHECKED. C A WEIGHTED COMBINATION OF THE OUTCOME OF THESE FOUR CHECKS IS C PRODUCED ON THE BASIS OF FACTORS WHICH DETERMINE WHICH C OF THE INTERPOLATIONS (IF ANY) MIGHT BE EXPECTED TO BE MORE C RELEVANT THAN THE OTHERS; FOR EXAMPLE THE NUMBER AND LOCATION C IN SPACE AND TIME OF OBSERVATIONS USED IN A GIVEN CHECK C AS WELL AS RELATIVE QUALITY OF THE OBSERVATIONS INVOLVED IS USED TO C INDICATE WHICH OF THE OUTCOMES SHOULD BE RELIED MORE HEAVILY UPON C TO MAKE A DECISION IN A PARTICULAR CASE. CHECK #5 INTERPOLATES C A VALUE FROM THE MULTIVARIATE SET OF OBS USED IN CHECKS #1-#3, C WHICH IS USED TO MEASURE HOW THE ANALYSIS WOULD DRAW (OR NOT DRAW) C TO THE OB BEING CHECKED. THE TOLERANCE ALLOWED FOR DEVIATION FROM C THE EXPECTED OUTCOME OF THE INTERPOLATION CHECKS IS PROPORTIONAL C TO A MEASURE OF THE DRAW COMPUTED USING THE RESULT OF CHECK #5. C DURING EACH ITERATION EACH OBSERVATION EITHER PASSES OR FAILS THE C COMPLEX OF CHECKS JUST DESCRIBED. A "PASS" MEANS THAT OB MAY BE C USED FOR CHECKING OBS DURING THE SUBSEQUENT ITERATION. A "FAIL" C INDICATES THE OPPOSITE. THE SYSTEM IS ITERATED UNTIL IDENTICAL C RESULTS IN TERMS OF OBS THAT PASS AND FAIL ARE OBTAINED IN TWO C CONSECUTIVE ITERATIONS UP TO FOUR COMPLETE ITERATIONS AT WHICH C POINT A FINAL ARBITRATION PROCEDURE IS INVOKED TO RESOLVE A (SMALL) C NUMBER OF AMBIGUOUS CASES WHICH ARE PREVENTING CONVERGENCE OF THE C SYSTEM. C C PROGRAM HISTORY LOG: C 1990-11-06 J. WOOLLEN ORIGINAL VERSION FOR IMPLEMENTATION C 1990-04-23 J. WOOLLEN VERSION FOR IMPLEMENTATION W/ SPECTRAL OI C 1998-09-25 J. WOOLLEN Y2K VERSION & RETURN CODE 4 FOR TOO MUCH DATA C 1999-09-27 J. WOOLLEN IBM-SP MPI/OMP VERSION IMPLEMENTED C 1999-12-21 J. WOOLLEN FIX FOR RECEIVE PROBLEM WITH TOO MUCH DATA C 2004-09-09 D. KEYSER CHANGE CALL FROM BUFRLIB ROUTINE BORT TO C W3LIB ROUTINE ERREXIT; BUMPED UP MAXREP FROM 200K TO 300K C (AND ALL ASSOC. ARRAYS), BUMPED UP MAXLEV FROM 300K TO 400K C (AND ALL ASSOC. ARRAYS) - NEEDED TO HANDLE METAR RPTS IN R. C TYPE 187; IMPROVED DOCBLOCKS C 2006-06-06 D. KEYSER BUMPED UP MAXREP FROM 300K TO 600K (AND ALL C ASSOC. ARRAYS), BUMPED UP MAXLEV FROM 400K TO 800K (AND ALL C ASSOC. ARRAYS) - NEEDED TO HANDLE WINDSAT RPTS IN R. TYPE C 289 C 2007-04-23 D. KEYSER BUMPED UP MAXREP FROM 600K TO 900K (AND ALL C ASSOC. ARRAYS), BUMPED UP MAXLEV FROM 800K TO 1100K (AND ALL C ASSOC. ARRAYS) - NEEDED TO HANDLE WINDSAT RPTS IN R. TYPE C 289 C 2007-04-23 J. WOOLLEN FIX IN SUBR NDDSPLIT TO AVOID LOAD BALANCING C ERROR MADE IN UNUSUAL SITUATIONS WITH SMALL AMOUNTS OF DATA C 2012-08-31 J. WOOLLEN FIX IN SUBR SEARCHS AND FCN LTX TO AVOID C ARRAY OVERFLOW IN IMAP FOR CERTAIN LAT/LON VALUES C 2012-11-20 J. WOOLLEN INITIAL PORT TO WCOSS: REMOVE ALL THREADING C SO CODE BECOMES STRICTLY MPI FOR SIMPLER AND BETTER LOAD C BALANCING C 2013-03-26 J. WOOLLEN FIX IN SUBR. PREPQM ADDRESSING ALIGNMENT C PROBLEM WHEN THE LAST BUFR MESSAGE CONTAINS NO DATA WHICH IS C OIQC PROCESSED C 2013-03-26 D. KEYSER FINAL CHANGES TO RUN ON WCOSS C 2013-06-19 J. WOOLLEN FIX SUBR. SYNCNDD TO PREVENT HANGUP IN C RECEIVE BY A PROCESSOR WHEN THERE IS VERY LITTLE DATA OF A C CERTAIN TYPE C 2019-06-20 S. Melchior Explicit declartion of bmiss in subr. STORE C C USAGE: C INPUT FILES: C UNIT 11 - NCEP PRODUCTION DATE FILE C UNIT 14 - OBSERVATION (PREPBUFR) PRIOR TO UPDATED QUALITY MARKS C UNIT 17 - OBSERVATION ERROR TABLE FILE SPECIFIC TO THIS PROGRAM C (THIS IS NOT THE ANALYSIS OBSERVATION ERROR FILE - C THOSE DATA ARE READ FROM THE PREPBUFR FILE) C UNIT 18 - LIST OF IPOINTS FOR DETAILED PRINT IN UNITS 71-80 C UNIT 20 - TOLERANCE LEVELS (FOR EXPERIMENTATION) C C OUTPUT FILES: C UNIT 06 - STANDARD OUTPUT PRINT C UNIT 61 - DETAILED TOSSLIST FOR SURFACE HEIGHTS C UNIT 62 - DETAILED TOSSLIST FOR TEMPS AND WINDS C UNIT 63 - DETAILED TOSSLIST FOR SAT TEMPS C UNIT 64 - DETAILED TOSSLIST FOR SSMI WIND SPEEDS C UNIT 65 - OFFICIAL OPERATIONAL TOSSLIST C UNIT 70 - OBSERVATION (PREPBUFR) FILE WITH UPDATED QUALITY MARKS C UNIT 71 - DETAILED PRINTOUTS FOR QC DIAGNOSIS (IPOINT IN UNIT 18) C UNIT 72 - DETAILED PRINTOUTS FOR QC DIAGNOSIS (IPOINT IN UNIT 18) C UNIT 73 - DETAILED PRINTOUTS FOR QC DIAGNOSIS (IPOINT IN UNIT 18) C UNIT 74 - DETAILED PRINTOUTS FOR QC DIAGNOSIS (IPOINT IN UNIT 18) C UNIT 75 - DETAILED PRINTOUTS FOR QC DIAGNOSIS (IPOINT IN UNIT 18) C UNIT 76 - DETAILED PRINTOUTS FOR QC DIAGNOSIS (IPOINT IN UNIT 18) C UNIT 77 - DETAILED PRINTOUTS FOR QC DIAGNOSIS (IPOINT IN UNIT 18) C UNIT 78 - DETAILED PRINTOUTS FOR QC DIAGNOSIS (IPOINT IN UNIT 18) C UNIT 79 - DETAILED PRINTOUTS FOR QC DIAGNOSIS (IPOINT IN UNIT 18) C UNIT 80 - DETAILED PRINTOUTS FOR QC DIAGNOSIS (IPOINT IN UNIT 18) C UNIT 81 - SUMMARY OF QC RESULTS (ASCII) C UNIT 82 - SUMMARY OF QC RESULTS (BINARY) C C SUBPROGRAMS CALLED: (*AUTOTASKED) C UNIQUE: - SETCON RDOBER STORE STOERR *COOR2 C MAKEAMAP SEARCH CHDIST QCLOOP QCPTS C SATPTS *SELDAT *QCOI *COOR1 *DRCTSL C *VSOLVE *OUTPUT *SATPUT OBOGRAM TOSSLIST C PREPQM SYNCSTOR MPISPLIT SSMIPUT SRTPRS C SEARCHS SSMIPTS SYNCNDD C LIBRARY: C MPI - MPI_BCST MPI_RECV MPI_SEND MPI_INIT C MPI_COMM_RANK MPI_COMM_SIZE C MPI_FINALIZE C MASS - VSCOS C W3NCO - W3TAGB W3TAGE ERREXIT C BUFRLIB - DATEBF CLOSBF UFBCPY OPENMB IREADMM C OPENBF UFBQCD UFBINT WRITSB UFBCNT C CLOSMG COPYMG ICOPYSB UFBMEM IREADSB C READMM C C EXIT STATES: C COND = 0 - SUCCESSFUL RUN C COND = 4 - TOO MUCH DATA, NOT ALL WAS CHECKED - NORMAL C TERMINATION C COND = 99 - DISASTOROUS PROBLEM SITUATION DETECTED - ABNORMAL C TERMINATION C C C REMARKS: NONE. C C ATTRIBUTES: C LANGUAGE: FORTRAN 90 C MACHINE: NCEP WCOSS C C$$$ C----------------------------------------------------------------------- CTREE GETFCER CTREE PILNLNP CTREE SATPTS CTREE SATPUT C----------------------------------------------------------------------- PROGRAM PREPOBS_OIQCBUFR PARAMETER (MAXREP=1000000) PARAMETER (MAXLEV=2000000) COMMON /OIQCDATA/ FLV (MAXREP ) , FLT (MAXREP ) , .FLN (MAXREP ) , FTM (MAXREP ) , FKX (MAXREP ) , .FNL (MAXREP ) , FMS (MAXREP ) , FMW (MAXREP,2 ) , .FQL (MAXREP ) , FHD (MAXLEV ) , FPRC(MAXLEV ) , .FPFD(MAXLEV ) , FPLV(MAXLEV ) , FPR (MAXLEV ) , .FQM (MAXLEV ) , FOI (MAXLEV ) , FOB (MAXLEV,2 ) , .FOE (MAXLEV ) , FFE (MAXLEV,2 ) , FTZ (MAXLEV ) , .FTW (MAXLEV,10,2) , FPF (MAXLEV,5,2 ) , FOC (MAXREP ) CHARACTER*3 FOC COMMON /OBS/ NREP,NLEV,SID(MAXREP),INOB(MAXREP),IMAP(360,181) COMMON /CKLIST / NDD,INDD(MAXLEV) COMMON /TOSSUNIT/ ITUNIT COMMON /RETCODE/ IRETCODE REAL*8 SID EXTERNAL QCPTS,SATPTS DATA LUDAT /11/ DATA LUBFI /14/ DATA LUQCE /17/ DATA LUBFQ /70/ DATA ITERS /6/ include "mpif.h" C----------------------------------------------------------------------- IRETCODE = 0 C----------------------------------------------------------------------- C UNIFIED QC ANALYSIS SETUP PROCEDURE IS AS FOLLOWS: C ------------------------------------------------- C -1) IRETCODE - INITIALLY 0, MAY BE CHANGED TO 4 IF TOO MUCH DATA C 0) W3TAGB - WRITE AN OPERATIONS START MESSAGE C 1) SETCON - SET CONSTANTS AND PARAMETERS C 2) STORE - READ IN FERR (OB/INCREMENT) FILE AND STORE IT C 3) RDOBER - READ IN THE OBSERVATION ERRORS C 4) STOERR - STORE THE OB AND FC ERRORS WITH EACH REPORT LEVEL C 5) MAKEAMAP - MAKE AN OB POINTER LIST SORTED BY LAT-LON C----------------------------------------------------------------------- call mpi_init(ierr) call mpi_comm_rank(MPI_COMM_WORLD,myid,ierr) call mpi_comm_size(MPI_COMM_WORLD,nprc,ierr) if(myid.eq.0) then print* print*,'nprocs = ',nprc print* CALL W3TAGB('PREPOBS_OIQCBUFR',2019,0171,0088,'NP22') print* print* print*, '--> WELCOME TO PREPOBS_OIQCBUFR - VERSION 06-20-2019' print* print* endif CALL SETCON CALL STORE(LUDAT,LUBFI) CALL RDOBER (LUQCE) CALL STOERR CALL MAKEAMAP C OIQC CHECKS OF SURFACE HEIGHT/PRESSURE DATA C ------------------------------------------- NDD = 0 DO IREP=1,NREP USEFLG = ABS(FQL(IREP)) FQL(IREP) = -USEFLG IF(FOC(IREP).EQ.'SFC') THEN FNL(IREP) = 1 ILEV = FLV(IREP)+1 IF(FOE(ILEV).GT.0.) THEN NDD = NDD+1 INDD(NDD) = ILEV FQL(IREP) = USEFLG FLV(IREP) = ILEV ENDIF ENDIF ENDDO IF(NDD.GT.0) THEN if(myid.eq.0) then WRITE(6 ,*)' CHECKING SURFACE HEIGHT DATA ',NDD,' OBS' endif ITUNIT = 61 CALL SEARCH CALL QCLOOP(QCPTS,ITERS) DO N=1,NDD IREP = FHD(INDD(N)) FLV(IREP) = INDD(N)-1 ENDDO ENDIF C OIQC CHECKS OF UPA AND SFC TEMP AND WIND DATA C ------------------------------------------- NDD = 0 DO IREP=1,NREP USEFLG = ABS(FQL(IREP)) FQL(IREP) = -USEFLG IF(FOC(IREP).EQ.'UPA' .OR. FOC(IREP).EQ.'SFC') THEN FQL(IREP) = USEFLG ILEV = FLV(IREP) KLEV = ILEV+FNL(IREP)-1 DO JLEV=ILEV,KLEV IF(FOE(JLEV).GT.0.) THEN NDD = NDD+1 INDD(NDD) = JLEV ENDIF ENDDO ENDIF ENDDO IF(NDD.GT.0) THEN if(myid.eq.0) then WRITE(6 ,*)' CHECKING TEMP AND WIND DATA ',NDD,' OBS' endif ITUNIT = 62 CALL SEARCH CALL QCLOOP(QCPTS,ITERS) ENDIF C OIQC CHECKS OF SATEM DATA C ------------------------- NDD = 0 DO IREP=1,NREP USEFLG = ABS(FQL(IREP)) IF(FOC(IREP).EQ.'SAT') THEN FQL(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 FQL(IREP) = -USEFLG ENDIF ENDDO IF(NDD.GT.0) THEN if(myid.eq.0) then WRITE(6 ,*)' CHECKING SAT TEMP DATA ',NDD,' OBS' endif ITUNIT = 63 CALL SEARCH CALL QCLOOP(SATPTS,1) ENDIF C DIRECTION ASSIGNMENT AND OIQC CHECKS FOR SSMI DATA C -------------------------------------------------- NDD = 0 DO IREP=1,NREP USEFLG = ABS(FQL(IREP)) IF(FOC(IREP).EQ.'SMI') THEN ILEV = FLV(IREP) IF(FOE(ILEV).GT.0.) THEN NDD = NDD+1 INDD(NDD) = ILEV FQL(IREP) = USEFLG ENDIF ELSE IF(FOC(IREP).EQ.'SAT') THEN FQL(IREP) = -USEFLG ELSE FQL(IREP) = 0. ENDIF ENDDO IF(NDD.GT.0) THEN if(myid.eq.0) then WRITE(6 ,*)' CHECKING SSMI DATA ',NDD,' OBS' endif ITUNIT = 64 CALL SSMIPTS call syncndd(fob(1,1),nlev,mpi_real) call syncndd(fob(1,2),nlev,mpi_real) CALL SEARCH CALL QCLOOP(QCPTS,ITERS) ENDIF C COMPUTE STATISTICS, WRITE LISTS, UPDATE QUALITY MARKS C ----------------------------------------------------- call mpi_barrier(mpi_comm_world,ier) if(myid.eq.0) then CALL OBOGRAM (11,81,82) CALL TOSSLIST (65) CALL PREPQM (LUBFI,LUBFQ) CALL W3TAGE('PREPOBS_OIQCBUFR') endif C TIE IT TOGETHER AND END UP C -------------------------- call mpi_finalize(mret) IF(IRETCODE.NE.0) THEN CALL ERREXIT(IRETCODE) ELSE STOP ENDIF END C$$$ SUBPROGRAM DOCUMENTATION BLOCK C C SUBPROGRAM: CHDIST C PRGMMR: WOOLLEN ORG: NP20 DATE: 1990-11-06 C C ABSTRACT: COMPUTES CHORD LENGTH DISTANCE FROM A FIXED POINT TO AN C ARRAY OF POINTS USING THE FORMULA: C S**2/2 = 1 - COS(Y1-Y2) + COS(Y1)*COS(Y2)*(1-COS(X1-X2)). C ALSO THE DIRECTION OF EACH POINT IS COMPUTED WITH RESPECT TO THE C FIXED POINT AND RETURNED AS WELL. C C PROGRAM HISTORY LOG: C 1990-11-06 J. WOOLLEN - ORIGINAL AUTHOR C C USAGE: CALL CHDIST(X1,Y1,X2,Y2,DIST,DIRN,NIND,NP) C INPUT ARGUMENT LIST: C X1 - X-COORDINATE (LONGITUDE) OF FIXED POINT C Y1 - Y-COORDINATE (LATITUDE ) OF FIXED POINT C X2 - X-COORDINATES (LONGITUDE) FOR SET OF POINTS C Y2 - Y-COORDINATES (LATITUDE ) FOR SET OF POINTS C NP - NUMBER OF OUTPUTS REQUESTED C C OUTPUT ARGUMENT LIST: C DIST - CHORD LENGTH DISTANCES FROM FIXED POINT (KM) C DIRN - DIRECTION FROM FIXED POINT (DEG) C C REMARKS: NONE. C C C ATTRIBUTES: C LANGUAGE: FORTRAN 90 C MACHINE: NCEP WCOSS C C$$$ SUBROUTINE CHDIST(X1,Y1,X2,Y2,DIST,DIRN,NIND,NP) DIMENSION X2(NP),Y2(NP),DIST(NP),DIRN(NP),NIND(NP) dimension cosy2(np),cosdx(np),cosdy(np) DATA PI180/.0174532 /,RADE/6371./ SAVE PI180,RADE C---------------------------------------------------------------------- !itab(x) = nint(mod(x+360.,360.)*100.) C---------------------------------------------------------------------- IF(NP.EQ.0) RETURN C COMPUTE THE DISTANCE C -------------------- COSY1 = cos(y1*pi180) do i=1,np cosy2(i) = cos(y2(i)*pi180) cosdx(i) = cos((x1-x2(i))*pi180) cosdy(i) = cos((y1-y2(i))*pi180) enddo !call vscos(cosy2,cosy2,np) !call vscos(cosdx,cosdx,np) !call vscos(cosdy,cosdy,np) DO I=1,NP S = 1.0-COSDY(I)+COSY1*COSY2(I)*(1.0-COSDX(I)) S = SQRT(2.*S) IF(S.LE..002) S = 0. DIST(I) = S*RADE ENDDO C COMPUTE DIRECTIONS C ------------------ DO I=1,NP DX = (X2(I)-X1)*COSY1 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. ENDDO RETURN END C$$$ SUBPROGRAM DOCUMENTATION BLOCK C C SUBPROGRAM: COOR1 C PRGMMR: WOOLLEN ORG: NP20 DATE: 1990-11-06 C C ABSTRACT: COMPUTES CHORD LENGTH DISTANCE FOR A MATRIX OF LOCATIONS C <(X,Y)I,(X,Y)J> USING THE NORMAL ANGLE FORMULA: C S**2/2 = 1 - COS(Y1-Y2) + COS(Y1)*COS(Y2)*(1-COS(X1-X2)) C DERIVATIVES OF DISTANCE WITH RESPECT TO LAT AND LON ARE ALSO C COMPUTED AND COMBINED WITH APPROPRIATE CORRELATION FUNCTIONS AND C DERIVATIVES WITH RESPECT TO DISTANCE TO FORM THE MULTI-VARIATE C CORRELATIONS FOR THE HEIGHT WIND ANALYSIS. C C PROGRAM HISTORY LOG: C 1990-11-06 J. WOOLLEN - ORINGIAL AUTHOR C C USAGE: CALL COOR1(NP,CH,FL,TZ,X1,Y1,K1,X2,Y2,K2,D,C1,C2) C INPUT ARGUMENT LIST: C NP - NUMBER OF OUTPUTS REQUESTED C CH - VECTOR OF CONSTANTS FOR THE GAUSSIAN LENGTH SCALE C X1 - X-COORDINATES (LONGITUDE) FOR POINT 1 C Y1 - Y-COORDINATES (LATITUDE ) FOR POINT 1 C K1 - VARIABLE TYPE AT POINT 1 (I.E. 1-Z , 2-U , 3-V) C X2 - X-COORDINATES (LONGITUDE) FOR POINT 2 C Y2 - Y-COORDINATES (LATITUDE ) FOR POINT 2 C K2 - VARIABLE TYPE AT POINT 2 (I.E. 1-Z , 2-U , 3-V) C FL - DECOUPLING FACTOR C TZ - HEIGHT-TEMPERATURE FACTOR C C OUTPUT ARGUMENT LIST: C D - DISTANCES BETWEEN ARRAYS OF POINTS IN RADIANS C C1 - CORRELATION COEFFICIENTS AS SPECIFIED C C2 - CORRELATION COEFFICIENTS AS SPECIFIED C C REMARKS: THIS PROGRAM COMPUTES A GAUSSIAN CORRELATION WITH C LENGTH SCALE GIVEN BY THE CHTWV INPUT ARGUMENT. C C C ATTRIBUTES: C LANGUAGE: FORTRAN 90 C MACHINE: NCEP WCOSS C C$$$ SUBROUTINE COOR1(NP,CH,FL,TZ,X1,Y1,K1,X2,Y2,K2,D,C1,C2) DIMENSION CH(NP),FL(NP),TZ(NP) DIMENSION X1(NP),Y1(NP),K1(NP) DIMENSION X2(NP),Y2(NP),K2(NP) DIMENSION D (NP),C1(NP),C2(NP) LOGICAL ZZ,ZU,ZV,UZ,UU,UV,VZ,VU,VV DATA PI180 /.0174532 /, SMIN/.002/ C---------------------------------------------------------------------- C---------------------------------------------------------------------- IF(NP.EQ.0) RETURN C LOOP OVER SET OF INPUT POINTS C ----------------------------- DO 20 I=1,NP C1(I) = 0. C2(I) = 0. C DETERMINE CORRELATION TYPE C -------------------------- ZZ = K1(I).EQ.1 .AND. K2(I).EQ.1 ZU = K1(I).EQ.1 .AND. K2(I).EQ.2 ZV = K1(I).EQ.1 .AND. K2(I).EQ.3 UZ = K1(I).EQ.2 .AND. K2(I).EQ.1 UU = K1(I).EQ.2 .AND. K2(I).EQ.2 UV = K1(I).EQ.2 .AND. K2(I).EQ.3 VZ = K1(I).EQ.3 .AND. K2(I).EQ.1 VU = K1(I).EQ.3 .AND. K2(I).EQ.2 VV = K1(I).EQ.3 .AND. K2(I).EQ.3 C COMPUTE THE MATRIX OF SINES AND COSINES C --------------------------------------- COSY1 = COS(Y1(I)*PI180) COSY2 = COS(Y2(I)*PI180) SINY1 = SIN(Y1(I)*PI180) SINY2 = SIN(Y2(I)*PI180) COSDX = COS((X1(I)-X2(I))*PI180) COSDY = COS((Y1(I)-Y2(I))*PI180) SINDX = SIN((X1(I)-X2(I))*PI180) SINDY = SIN((Y1(I)-Y2(I))*PI180) C COMPUTE THE NORMAL ANGLE C ------------------------ S = 1.0-COSDY+COSY1*COSY2*(1.0-COSDX) S = SQRT(2.*S) C CALCULATE THE FIRST PARTIALS WITH RESPECT TO X1 AND Y1 C ------------------------------------------------------ SI = MAX(S,SMIN) SI = 1./SI DY1 = SI*(SINDY-SINY1*COSY2*(1.0-COSDX)) DX1 = SI*(COSY1*COSY2*SINDX) DY2 = -SI*(SINDY+COSY1*SINY2*(1.0-COSDX)) DX2 = -DX1 C CALCULATE THE MIXED PARTIAL DERIVATIVES C --------------------------------------- DY1Y2= SI*(SINY1*SINY2*(1.0-COSDX)-COSDY-DY1*DY2) DY1X2= SI*(SINY1*COSY2*SINDX-DY1*DX2) DX1X2= -SI*(COSY1*COSY2*COSDX+DX1*DX2) DX1Y2= -SI*(COSY1*SINY2*SINDX+DX1*DY2) C COMPUTE THE VARIOUS CORRELATIONS C -------------------------------- RHO = EXP(-CH(I)*S*S) DRHO = -2.*CH(I)*S*RHO D2RHO = -2.*CH(I)*(RHO+S*DRHO) WW0 = 2.*CH(I) SRWW0 = SQRT(WW0) D(I) = S IF(S.LT.SMIN) THEN IF(ZZ) C1(I) = 1. IF(UU) C1(I) = 1. IF(VV) C1(I) = 1. D(I) = 0. ELSE IF(ZZ) THEN C1(I) = RHO ELSE IF(ZU) THEN C1(I) = FL(I)*TZ(I) * DRHO*DY2*SIGN(1.,-SINY2)/SRWW0 ELSE IF(ZV) THEN C1(I) = FL(I)*TZ(I) * DRHO*DX2*SIGN(1., SINY2)/(SRWW0*COSY2) ELSE IF(UZ) THEN C1(I) = FL(I)*TZ(I) * DRHO*DY1*SIGN(1.,-SINY1)/SRWW0 ELSE IF(UU) THEN C1(I) = DY1*DY2*D2RHO+DY1Y2*DRHO C1(I) = C1(I)*SIGN(1., SINY1*SINY2)/WW0 C1(I) = C1(I)*FL(I) + RHO*(1.-FL(I)) ELSE IF(UV) THEN C1(I) = DY1*DX2*D2RHO+DY1X2*DRHO C1(I) = FL(I) * C1(I)*SIGN(1.,-SINY1*SINY2)/(WW0*COSY2) ELSE IF(VZ) THEN C1(I) = FL(I)*TZ(I) * DRHO*DX1*SIGN(1., SINY1)/(SRWW0*COSY1) ELSE IF(VU) THEN C1(I) = DX1*DY2*D2RHO+DX1Y2*DRHO C1(I) = FL(I) * C1(I)*SIGN(1.,-SINY1*SINY2)/(WW0*COSY1) ELSE IF(VV) THEN C1(I) = DX1*DX2*D2RHO+DX1X2*DRHO C1(I) = C1(I)*SIGN(1., SINY1*SINY2)/(WW0*COSY1*COSY2) C1(I) = C1(I)*FL(I) + RHO*(1.-FL(I)) ENDIF 20 CONTINUE RETURN END C$$$ SUBPROGRAM DOCUMENTATION BLOCK C C SUBPROGRAM: COOR2 C PRGMMR: WOOLLEN ORG: NP20 DATE: 1990-11-06 C C ABSTRACT: COMPUTES CHORD LENGTH DISTANCE FOR A MATRIX OF LOCATIONS C <(X,Y)I,(X,Y)J> USING THE NORMAL ANGLE FORMULA: C S**2/2 = 1 - COS(Y1-Y2) + COS(Y1)*COS(Y2)*(1-COS(X1-X2)) C DERIVATIVES OF DISTANCE WITH RESPECT TO LAT AND LON ARE ALSO C COMPUTED AND COMBINED WITH APPROPRIATE CORRELATION FUNCTIONS AND C DERIVATIVES WITH RESPECT TO DISTANCE TO FORM THE MULTI-VARIATE C CORRELATIONS FOR THE HEIGHT WIND ANALYSIS. C C PROGRAM HISTORY LOG: C 1990-11-06 J. WOOLLEN - ORINGIAL AUTHOR C C USAGE: CALL COOR2(NP,CH,FL,TZ,X1,Y1,K1,X2,Y2,K2,D,C1,C2) C INPUT ARGUMENT LIST: C NP - NUMBER OF OUTPUTS REQUESTED C CH - VECTOR OF CONSTANTS FOR THE GAUSSIAN LENGTH SCALE C X1 - X-COORDINATES (LONGITUDE) FOR POINT 1 C Y1 - Y-COORDINATES (LATITUDE ) FOR POINT 1 C K1 - VARIABLE TYPE AT POINT 1 (I.E. 1-Z , 2-U , 3-V) C X2 - X-COORDINATES (LONGITUDE) FOR POINT 2 C Y2 - Y-COORDINATES (LATITUDE ) FOR POINT 2 C K2 - VARIABLE TYPE AT POINT 2 (I.E. 1-Z , 2-U , 3-V) C FL - DECOUPLING FACTOR C TZ - HEIGHT-TEMPERATURE FACTOR C C OUTPUT ARGUMENT LIST: C D - DISTANCES BETWEEN ARRAYS OF POINTS IN RADIANS C C1 - CORRELATION COEFFICIENTS AS SPECIFIED C C2 - CORRELATION COEFFICIENTS AS SPECIFIED C C REMARKS: THIS PROGRAM COMPUTES A GAUSSIAN CORRELATION WITH C LENGTH SCALE GIVEN BY THE CHTWV INPUT ARGUMENT. C C C ATTRIBUTES: C LANGUAGE: FORTRAN 90 C MACHINE: NCEP WCOSS C C$$$ SUBROUTINE COOR2(NP,CH,FL,TZ,X1,Y1,K1,X2,Y2,K2,D,C1,C2) DIMENSION CH(NP),FL(NP),TZ(NP) DIMENSION X1(NP),Y1(NP),K1(NP) DIMENSION X2(NP),Y2(NP),K2(NP) DIMENSION D (NP),C1(NP),C2(NP) LOGICAL ZZ,ZU,ZV,UZ,UU,UV,VZ,VU,VV DATA PI180 /.0174532 /, SMIN/.002/ C---------------------------------------------------------------------- C---------------------------------------------------------------------- IF(NP.EQ.0) RETURN C LOOP OVER SET OF INPUT POINTS C ----------------------------- DO 20 I=1,NP C1(I) = 0. C2(I) = 0. C DETERMINE CORRELATION TYPE C -------------------------- ZZ = K1(I).EQ.1 .AND. K2(I).EQ.1 ZU = K1(I).EQ.1 .AND. K2(I).EQ.2 ZV = K1(I).EQ.1 .AND. K2(I).EQ.3 UZ = K1(I).EQ.2 .AND. K2(I).EQ.1 UU = K1(I).EQ.2 .AND. K2(I).EQ.2 UV = K1(I).EQ.2 .AND. K2(I).EQ.3 VZ = K1(I).EQ.3 .AND. K2(I).EQ.1 VU = K1(I).EQ.3 .AND. K2(I).EQ.2 VV = K1(I).EQ.3 .AND. K2(I).EQ.3 C COMPUTE THE MATRIX OF SINES AND COSINES C --------------------------------------- COSY1 = COS(Y1(I)*PI180) COSY2 = COS(Y2(I)*PI180) SINY1 = SIN(Y1(I)*PI180) SINY2 = SIN(Y2(I)*PI180) COSDX = COS((X1(I)-X2(I))*PI180) COSDY = COS((Y1(I)-Y2(I))*PI180) SINDX = SIN((X1(I)-X2(I))*PI180) SINDY = SIN((Y1(I)-Y2(I))*PI180) C COMPUTE THE NORMAL ANGLE C ------------------------ S = 1.0-COSDY+COSY1*COSY2*(1.0-COSDX) S = SQRT(2.*S) C CALCULATE THE FIRST PARTIALS WITH RESPECT TO X1 AND Y1 C ------------------------------------------------------ SI = MAX(S,SMIN) SI = 1./SI DY1 = SI*(SINDY-SINY1*COSY2*(1.0-COSDX)) DX1 = SI*(COSY1*COSY2*SINDX) DY2 = -SI*(SINDY+COSY1*SINY2*(1.0-COSDX)) DX2 = -DX1 C CALCULATE THE MIXED PARTIAL DERIVATIVES C --------------------------------------- DY1Y2= SI*(SINY1*SINY2*(1.0-COSDX)-COSDY-DY1*DY2) DY1X2= SI*(SINY1*COSY2*SINDX-DY1*DX2) DX1X2= -SI*(COSY1*COSY2*COSDX+DX1*DX2) DX1Y2= -SI*(COSY1*SINY2*SINDX+DX1*DY2) C COMPUTE THE VARIOUS CORRELATIONS C -------------------------------- RHO = EXP(-CH(I)*S*S) DRHO = -2.*CH(I)*S*RHO D2RHO = -2.*CH(I)*(RHO+S*DRHO) WW0 = 2.*CH(I) SRWW0 = SQRT(WW0) D(I) = S IF(S.LT.SMIN) THEN IF(ZZ) C1(I) = 1. IF(ZZ) C2(I) = 1. IF(UU) C1(I) = 1. IF(VV) C1(I) = 1. IF(VU) C2(I) = 1. D(I) = 0. ELSE IF(ZZ) THEN C1(I) = RHO C2(I) = C1(I) ELSE IF(ZU) THEN C1(I) = FL(I)*TZ(I) * DRHO*DY2*SIGN(1.,-SINY2)/SRWW0 C2(I) = FL(I)*TZ(I) * DRHO*DX2*SIGN(1., SINY2)/(SRWW0*COSY2) ELSE IF(ZV) THEN C1(I) = FL(I)*TZ(I) * DRHO*DX2*SIGN(1., SINY2)/(SRWW0*COSY2) C2(I) = C1(I) ELSE IF(UZ) THEN C1(I) = FL(I)*TZ(I) * DRHO*DY1*SIGN(1.,-SINY1)/SRWW0 C2(I) = C1(I) ELSE IF(UU) THEN C1(I) = DY1*DY2*D2RHO+DY1Y2*DRHO C1(I) = C1(I)*SIGN(1., SINY1*SINY2)/WW0 C1(I) = C1(I)*FL(I) + RHO*(1.-FL(I)) C2(I) = DY1*DX2*D2RHO+DY1X2*DRHO C2(I) = FL(I) * C2(I)*SIGN(1.,-SINY1*SINY2)/(WW0*COSY2) ELSE IF(UV) THEN C1(I) = DY1*DX2*D2RHO+DY1X2*DRHO C1(I) = FL(I) * C1(I)*SIGN(1.,-SINY1*SINY2)/(WW0*COSY2) C2(I) = C1(I) ELSE IF(VZ) THEN C1(I) = FL(I)*TZ(I) * DRHO*DX1*SIGN(1., SINY1)/(SRWW0*COSY1) C2(I) = C1(I) ELSE IF(VU) THEN C1(I) = DX1*DY2*D2RHO+DX1Y2*DRHO C1(I) = FL(I) * C1(I)*SIGN(1.,-SINY1*SINY2)/(WW0*COSY1) C2(I) = DX1*DX2*D2RHO+DX1X2*DRHO C2(I) = C2(I)*SIGN(1., SINY1*SINY2)/(WW0*COSY1*COSY2) C2(I) = C2(I)*FL(I) + RHO*(1.-FL(I)) ELSE IF(VV) THEN C1(I) = DX1*DX2*D2RHO+DX1X2*DRHO C1(I) = C1(I)*SIGN(1., SINY1*SINY2)/(WW0*COSY1*COSY2) C1(I) = C1(I)*FL(I) + RHO*(1.-FL(I)) C2(I) = C1(I) ENDIF 20 CONTINUE RETURN END C$$$ SUBPROGRAM DOCUMENTATION BLOCK C C SUBPROGRAM: DRCTSL C PRGMMR: WOOLLEN ORG: NP20 DATE: 1990-11-06 C C ABSTRACT: DRIVER FOR CHOLESKY TYPE LINEAR EQUATION SOLVER. C C PROGRAM HISTORY LOG: C 1990-11-06 J. WOOLLEN - ORIGINAL AUTHOR C C USAGE: CALL DRCTSL(FA,RA,DP,NDIM,NXXYY) C INPUT ARGUMENT LIST: C FAALL - ARRAY OF SYMMETRIC MATRIXES C RAALL - ARRAY OF MATRIX RIGHT HAND SIDES C NDIM - ARRAY OF MATRIX RANKS C MAXDIM - MAXIMUM RANK MATRIX IN MATRIX ARRAY C NXXYY - NUMBER OF MATRIXES IN MATRIX ARRAY C NFT - NUMBER OF RIGHT HAND SIDE VECTORS PER MATRIX C OUTPUT ARGUMENT LIST: C RAALL - ARRAY OF MATRIX SOLUTIONS C DOTPRD - ARRAY OF DOT PRODUCTS OF RHS VECTORS WITH MATRIX C SOLUTIONS C OUTPUT FILES: C UNIT 06 - STANDARD OUTPUT PRINT C C REMARKS: ILL CONDITIONED OR NON POSITIVE DEFINITE MATRIXES ARE C IDENTIFIED BY DOT PRODUCTS GT 1 OR LT 0 OR BY A MESSAGE FROM C VSOLVE. FIVE ATTEMPTS ARE MADE TO SOLVE BAD ONES, BY RIDGE C REGRESSION, AFTER WHICH A NULL SOLUTION IS RETURNED. C C ATTRIBUTES: C LANGUAGE: FORTRAN 90 C MACHINE: NCEP WCOSS C C$$$ SUBROUTINE DRCTSL(FA,RA,DP,NDIM,NXXYY) PARAMETER( NMX = 64 , . NN = 10*(10+1)/2 , . N = 10 ) DIMENSION FA(NMX,NN), RA(NMX,N,2), DP(NMX,2), NDIM(NXXYY) DIMENSION A (NMX,NN), B (NMX,N,2) DIMENSION BAD(NMX), SMOOTH(6) LOGICAL BAD DATA SMOOTH /1.00,1.01,1.02,1.05,1.10,2.00/ C---------------------------------------------------------------------- C---------------------------------------------------------------------- C LOOP FOR POSSIBILITY OF BAD MATRIXS C ----------------------------------- DO 50 ITRY=1,6 NBAD = 0 C FIND THE LARGEST MATRIX IN THE STACK C ------------------------------------ MAXDIM = 0 DO 5 I=1,NXXYY MAXDIM = MAX(MAXDIM,NDIM(I)) 5 CONTINUE IF(MAXDIM.EQ.0) RETURN C INITIALIZE THE WORKING ARRAYS C ----------------------------- DO 10 J=1,MAXDIM*(MAXDIM+1)/2 DO 10 I=1,NXXYY A(I,J) = FA(I,J) 10 CONTINUE DO 11 J=1,MAXDIM DO 11 I=1,NXXYY B(I,J,1) = RA(I,J,1) B(I,J,2) = RA(I,J,2) 11 CONTINUE DO 12 I=1,NXXYY DP(I,1) = 0. DP(I,2) = 0. 12 CONTINUE DO 13 J=1,MAXDIM JJ = J*(J+1)/2 DO 13 I=1,NXXYY A(I,JJ) = FA(I,JJ)*SMOOTH(ITRY) 13 CONTINUE C CALL THE MATRIX SOLVER C ---------------------- CALL VSOLVE (A,B,NDIM,BAD,NXXYY,MAXDIM) C MAKE THE DOT PRODUCTS OF SOLUTIONS WITH RIGHT HAND SIDES C -------------------------------------------------------- DO 20 J=1,MAXDIM DO 20 I=1,NXXYY DP(I,1) = DP(I,1) + RA(I,J,1)*B(I,J,1) DP(I,2) = DP(I,2) + RA(I,J,2)*B(I,J,2) 20 CONTINUE C NORMALIZE THE WEIGHTS IF THEY ARE TOO BIG C ----------------------------------------- DO 25 I=1,NXXYY IF(DP(I,1).GT.1. .OR. DP(I,2).GT.1.) THEN DO 24 J=1,MAXDIM B(I,J,1) = B(I,J,1)/DP(I,1) B(I,J,2) = B(I,J,2)/DP(I,2) 24 CONTINUE DP(I,1) = 1. DP(I,2) = 1. ENDIF 25 CONTINUE C CHECK FOR BAD ONES - DO IT AGAIN IF NECESSARY C --------------------------------------------- DO 30 I=1,NXXYY BAD(I) = BAD(I) .OR. DP(I,1).LT.0. .OR. DP(I,2).LT.0. BAD(I) = BAD(I) .OR. DP(I,1).GT.1. .OR. DP(I,2).GT.1. 30 CONTINUE DO 40 I=1,NXXYY IF(BAD(I)) THEN DP(I,1) = 0. DP(I,2) = 0. DO 35 J=1,MAXDIM B(I,J,1) = 0. B(I,J,2) = 0. 35 CONTINUE NBAD = NBAD + 1 ENDIF 40 CONTINUE IF(NBAD.NE.0) PRINT*, 'NBAD=',NBAD,' ON TRY ',ITRY IF(NBAD.EQ.0) GOTO 55 50 CONTINUE C COPY SOLUTIONS INTO OUTPUT ARRAY - ZERO BAD ONES OUT C ---------------------------------------------------- 55 DO 60 J=1,MAXDIM DO 60 I=1,NXXYY RA(I,J,1) = B(I,J,1) RA(I,J,2) = B(I,J,2) 60 CONTINUE RETURN END C$$$ SUBPROGRAM DOCUMENTATION BLOCK C C SUBPROGRAM: GETFCER C PRGMMR: WOOLLEN ORG: NP20 DATE: 1992-07-29 C C ABSTRACT: FUNCTION WHICH RETURNS A FORECAST ERROR IN THE FORCAST C HEIGHT OR WIND FIELD FOR A GIVEN LAT, LON, PRESSURE. ENTRY POINTS C GETZER AND GETWER RETURN ERROR VALUES; GETFCER RETURNS 0. C C PROGRAM HISTORY LOG: C 1992-07-29 J. WOOLLEN - ORIGINAL AUTHOR C C USAGE: CALL GETZER(SLAT,SLON,PRES) - (RETURNS HEIGHT ERROR) C CALL GETWER(SLAT,SLON,PRES) - (RETURNS WIND ERROR) C INPUT ARGUMENT LIST: C SLAT - LATITUDE C SLON - LONGITUDE C PRES - PRESSURE C C FUNCTION RETURN VALUE: C GETZER - FORECAST HEIGHT ERROR C GETWER - FORECAST WIND ERROR C C REMARKS: NONE. C C C ATTRIBUTES: C LANGUAGE: FORTRAN 90 C MACHINE: NCEP WCOSS C C$$$ FUNCTION GETFCER(SLAT,SLON,PRES) COMMON /OBERRS/ QIKE(21,300),PMAND(21),P50(21),CHL(21) COMMON /CORP / CHLP(1500),CVC,FFLAT(91),IFA(10,10) COMMON /FCRATS/ RATION(21),RATIOS(21),TROPUV(21) LOGICAL WIND DATA PI180 / .0174532 / DATA OMEGA / .72921E-4 / DATA RADE / 6.371E6 / DATA GRAV / 9.8 / C--------------------------------------------------------------------- C--------------------------------------------------------------------- GETFCER = 0. RETURN ENTRY GETZER(SLAT,SLON,PRES) WIND = .FALSE. GOTO 10 ENTRY GETWER(SLAT,SLON,PRES) WIND = .TRUE. C COMPUTE A HEIGHT ERROR PROPORTIONAL TO RAWINSONDE ERROR C ------------------------------------------------------- 10 RAWZER = PILNLNP(PRES,PMAND,QIKE(1,120),21) IF(SLAT.LT.0.) THEN ZER = RAWZER/SQRT(PILNLNP(PRES,PMAND,RATIOS,21)) * 1.1 ELSE ZER = RAWZER/SQRT(PILNLNP(PRES,PMAND,RATION,21)) ENDIF IF(.NOT.WIND) THEN GETZER = ZER RETURN ENDIF C COMPUTE THE GEOSTROPHIC WIND ERROR FOR GETWER C --------------------------------------------- TRPWER = PILNLNP(PRES,PMAND,TROPUV,21) SRWW0 = SQRT(2.0*CHLP(MAX(IFIX(PRES),1))) IF(SLAT.EQ.0.) ALPHA = 0. IF(SLAT.NE.0.) ALPHA = GRAV/(2.*OMEGA*SIN(ABS(SLAT)*PI180)*RADE) IY = ABS(SLAT)+1.5 G = FFLAT(IY) GETWER = G*ZER*ALPHA*SRWW0 + (1-G)*TRPWER RETURN END C----------------------------------------------------------------------- C C----------------------------------------------------------------------- SUBROUTINE MAKEAMAP PARAMETER (MAXREP=1000000) PARAMETER (MAXLEV=2000000) COMMON /OIQCDATA/ FLV (MAXREP ) , FLT (MAXREP ) , .FLN (MAXREP ) , FTM (MAXREP ) , FKX (MAXREP ) , .FNL (MAXREP ) , FMS (MAXREP ) , FMW (MAXREP,2 ) , .FQL (MAXREP ) , FHD (MAXLEV ) , FPRC(MAXLEV ) , .FPFD(MAXLEV ) , FPLV(MAXLEV ) , FPR (MAXLEV ) , .FQM (MAXLEV ) , FOI (MAXLEV ) , FOB (MAXLEV,2 ) , .FOE (MAXLEV ) , FFE (MAXLEV,2 ) , FTZ (MAXLEV ) , .FTW (MAXLEV,10,2) , FPF (MAXLEV,5,2 ) , FOC (MAXREP ) CHARACTER*3 FOC COMMON /OBS/ NREP,NLEV,SID(MAXREP),INOB(MAXREP),IMAP(360,181) COMMON /CKLIST/ NDD,INDD(MAXLEV) REAL*8 SID C----------------------------------------------------------------------- C----------------------------------------------------------------------- C MAKE SURE NONE OF THE LONGITUDES EQUAL 360. C ------------------------------------------- DO 1 N=1,NREP IF(FLN(N).EQ.360.) FLN(N) = 0. 1 CONTINUE C LONGITUDE SORT C -------------- NN = 0 DO 10 NX=1,360 DO 10 N=1,NREP IX = FLN(N)+1. IF(IX.EQ.NX) THEN NN = NN+1 INDD(NN) = N ENDIF 10 CONTINUE C LATITUDE SORT C ------------- NOB = 0 DO 20 NY=1,181 DO 20 N=1,NN IY = FLT(INDD(N))+91. IF(IY.EQ.NY) THEN NOB = NOB+1 INOB(NOB) = INDD(N) ENDIF 20 CONTINUE C CHECK TO MAKE SURE ALL REPS ARE ACCOUNTED FOR C --------------------------------------------- IF(NREP.NE.NN .OR. NREP.NE.NOB) THEN PRINT* PRINT*,'*****************************************************' PRINT*,'** WARNING*WARNING*WARNING*WARNING*WARNING*WARNING **' PRINT*,'*****************************************************' PRINT*,'** IN MAKEAMAP ALL REPORTS HAVE NOT BEEN PROCESSED **' PRINT*,'*****************************************************' PRINT*,'** WARNING*WARNING*WARNING*WARNING*WARNING*WARNING **' PRINT*,'*****************************************************' PRINT* ENDIF C INITIALIZE IMAP ARRAY C --------------------- 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 C FILL GAPS IN IMAP FOR CONTINUITY (BACKWARDS) C -------------------------------------------- LASTN = NOB DO 40 NY=181,1,-1 DO 40 NX=360,1,-1 IF(IMAP(NX,NY).EQ.0) THEN IMAP(NX,NY) = LASTN ELSE LASTN = IMAP(NX,NY) ENDIF 40 CONTINUE RETURN END C$$$ SUBPROGRAM DOCUMENTATION BLOCK C C SUBPROGRAM: OBOGRAM C PRGMMR: WOOLLEN ORG: NP20 DATE: 1990-11-06 C C ABSTRACT: EXAMINES THE RESULTS OF FINAL QC DECISIONS AND MAKES TABLES C OF STATISTICS FOR ACCEPTED AND REJECTED REPORTS BY REPORT TYPE AND C LEVEL. C C PROGRAM HISTORY LOG: C 1990-11-06 J. WOOLLEN - ORIGINAL AUTHOR C C USAGE: CALL OBOGRAM(LDAT,LASC,LBIN) C INPUT ARGUMENT LIST: C LDAT - LOGICAL UNIT FOR NCEP PRODUCTION DATE FILE C LASC - LOGICAL UNIT FOR OUTPUT ASCII FILE C LBIN - LOGICAL UNIT FOR OUTPUT BINARY FILE C C INPUT FILES: C UNIT LDAT - NCEP PRODUCTION DATE FILE C C OUTPUT FILES: C UNIT 06 - STANDARD OUTPUT PRINT C UNIT LASC - SUMMARY OF QC RESULTS (ASCII) C UNIT LBIN - SUMMARY OF QC RESULTS (BINARY) C C REMARKS: NONE. C C ATTRIBUTES: C LANGUAGE: FORTRAN 90 C MACHINE: NCEP WCOSS C C$$$ SUBROUTINE OBOGRAM(LDAT,LASC,LBIN) PARAMETER (MAXREP=1000000) PARAMETER (MAXLEV=2000000) COMMON /OIQCDATA/ FLV (MAXREP ) , FLT (MAXREP ) , .FLN (MAXREP ) , FTM (MAXREP ) , FKX (MAXREP ) , .FNL (MAXREP ) , FMS (MAXREP ) , FMW (MAXREP,2 ) , .FQL (MAXREP ) , FHD (MAXLEV ) , FPRC(MAXLEV ) , .FPFD(MAXLEV ) , FPLV(MAXLEV ) , FPR (MAXLEV ) , .FQM (MAXLEV ) , FOI (MAXLEV ) , FOB (MAXLEV,2 ) , .FOE (MAXLEV ) , FFE (MAXLEV,2 ) , FTZ (MAXLEV ) , .FTW (MAXLEV,10,2) , FPF (MAXLEV,5,2 ) , FOC (MAXREP ) CHARACTER*3 FOC 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 DIMENSION STAT(300,15,6,2,2,3),IOBS(300),RZONE(6) LOGICAL QOBS(300,6) REAL*8 SID DATA AR/'ACC','REJ'/ DATA XMSG/32767./ C---------------------------------------------------------------------- C----------------------------------------------------------------------- REWIND (LDAT) READ(LDAT,1) NMCDATE 1 FORMAT(8X,I8) WRITE(LASC,*) WRITE(LASC,*)' STATS FOR OBS CHECKED ',NMCDATE WRITE(LASC,*) C START UP C -------- QOBS = .FALSE. STAT = 0 C GO THROUGH THE BIG BUFFER TO COUNT OBS FOR OBOGRAM C -------------------------------------------------- DO 10 IREP=1,NREP X = FLN(IREP) Y = FLT(IREP) KX = FKX(IREP) KK = KX/100 LB = MIN(IFIX((Y+90.)/30.+1.),6) TAG = FOC(IREP) IF(TAG.EQ.'SFC') FNL(IREP) = 2 ILEV = FLV(IREP) KLEV = ILEV+FNL(IREP)-1 C GO THROUGH LEVELS FOR ACCEPT REJECT STATISTICS C ---------------------------------------------- DO 10 JLEV=ILEV,KLEV IF(FOI(JLEV).NE.0) THEN IF(FOI(JLEV).EQ.1) THEN IAR = 1 ELSE IF(FOI(JLEV).EQ.-1) THEN IAR = 2 ELSE GOTO 10 ENDIF BUFF11 = FOB(JLEV,1) BUFF12 = FOB(JLEV,2) IF(TAG.EQ.'SFC' .AND. JLEV.EQ.KLEV) THEN KX = KX - 100 ELSE IF(KK.EQ.1) BUFF12 = BUFF11 ENDIF QOBS(KX,LB) = .TRUE. PP = FPR(JLEV) DO 5 L=1,15 IF(L.EQ.1) THEN IF(PP.GT.(PMAND(L)+PMAND(L+1))/2.) GOTO 6 ELSE IF(L.EQ.15) THEN IF(PP.LE.(PMAND(L)+PMAND(L-1))/2.) GOTO 6 ELSE IF(PP.GT.(PMAND(L)+PMAND(L+1))/2. .AND. . PP.LE.(PMAND(L)+PMAND(L-1))/2.) GOTO 6 ENDIF 5 CONTINUE PRINT *, 'OBOGRAM - UNABLE TO LOCATE DATA LEVEL' CALL W3TAGE('PREPOBS_OIQCBUFR') CALL ERREXIT(99) 6 CONTINUE IF(KK.EQ.1) THEN STAT(KX,L,LB,IAR,1,1) = STAT(KX,L,LB,IAR,1,1) + 1. STAT(KX,L,LB,IAR,1,2) = STAT(KX,L,LB,IAR,1,2) + BUFF12 STAT(KX,L,LB,IAR,1,3) = STAT(KX,L,LB,IAR,1,3) + BUFF12**2 ELSE STAT(KX,L,LB,IAR,1,1) = STAT(KX,L,LB,IAR,1,1) + 1. STAT(KX,L,LB,IAR,1,2) = STAT(KX,L,LB,IAR,1,2) + BUFF11 STAT(KX,L,LB,IAR,1,3) = STAT(KX,L,LB,IAR,1,3) + BUFF11**2 STAT(KX,L,LB,IAR,2,1) = STAT(KX,L,LB,IAR,2,1) + 1. STAT(KX,L,LB,IAR,2,2) = STAT(KX,L,LB,IAR,2,2) + BUFF12 STAT(KX,L,LB,IAR,2,3) = STAT(KX,L,LB,IAR,2,3) + BUFF12**2 ENDIF ENDIF 10 CONTINUE C COMPUTE THE AVERAGES FOR MEAN AND RMS'S (VECTOR SUMS FOR WINDS) C --------------------------------------------------------------- DO 20 KX=1,299 KK = MAX(KX/100,1) DO 20 L=1,15 DO 20 LB=1,6 DO 20 IAR=1,2 DO 20 KKX=1,KK N = STAT(KX,L,LB,IAR,KKX,1) IF(N.NE.0) THEN BIAS = STAT(KX,L,LB,IAR,KKX,2)/N STAT(KX,L,LB,IAR,KKX,2) = BIAS RMS = SQRT(STAT(KX,L,LB,IAR,KKX,3)/N) STAT(KX,L,LB,IAR,KKX,3) = RMS IF(KKX.EQ.2) THEN STAT(KX,L,LB,IAR,1,2) = SQRT(STAT(KX,L,LB,IAR,1,2)**2+BIAS**2) STAT(KX,L,LB,IAR,1,3) = SQRT(STAT(KX,L,LB,IAR,1,3)**2+ RMS**2) ENDIF ENDIF 20 CONTINUE C PRINT THE TOTAL, MEAN, AND RMS FOR ACCEPTS AND REJECTS C ------------------------------------------------------ DO 35 KX=1,300 DO 35 LB=1,6 IF(QOBS(KX,LB)) THEN KK = KX/100 WRITE(LASC,*)'------ZONE ',LB,' KX ',KX,' VAR ',KK,'------' DO 30 IAR=1,2 WRITE(LASC,2000) AR(IAR),(STAT(KX,L,LB,IAR,1,1),L=1,15) WRITE(LASC,2001) AR(IAR),(STAT(KX,L,LB,IAR,1,2),L=1,15) WRITE(LASC,2002) AR(IAR),(STAT(KX,L,LB,IAR,1,3),L=1,15) WRITE(LASC,2003) 30 CONTINUE ENDIF 35 CONTINUE C PRINT THE TOTAL REJECTED OB BY TYPE ZONE LEVEL C ---------------------------------------------- WRITE(LASC,*) WRITE(LASC,*) ' NUMBER OF OBS REJECTED BY OB TYPE, AND ZONE' WRITE(LASC,*) DO 45 KX=1,300 IF(QOBS(KX,1).OR.QOBS(KX,2).OR.QOBS(KX,3).OR.QOBS(KX,4).OR. . QOBS(KX,5).OR.QOBS(KX,6)) THEN RZONE = 0 DO 40 LB=1,6 DO 40 L=1,15 RZONE(LB) = STAT(KX,L,LB,2,1,1) + RZONE(LB) 40 CONTINUE WRITE(LASC,2004)KX,(RZONE(LB),LB=1,6) ENDIF 45 CONTINUE C WRITE A BINARY VERSION OF THE OBOGRAM C ------------------------------------- I = 0 DO 55 KX=1,300 DO 50 LB=1,6 IF(QOBS(KX,LB)) THEN I = I+1 IOBS(I) = KX GOTO 55 ENDIF 50 CONTINUE 55 CONTINUE WRITE(LBIN) NMCDATE, . I,(IOBS(II),((((STAT(IOBS(II),LEV,LB,IAR,1,IST), . LEV=1,15), . IST=1,3), . IAR=1,2), . LB=1,6), . II=1,I) C FORMAT STATEMENTS FOR THE OBOGRAM C --------------------------------- 2000 FORMAT('TOTAL ',A3,15F6.1) 2001 FORMAT('MEAN ',A3,15F6.1) 2002 FORMAT('RMS ',A3,15F6.1) 2003 FORMAT(' ') 2004 FORMAT(I5,6F8.1) RETURN END C$$$ SUBPROGRAM DOCUMENTATION BLOCK C C SUBPROGRAM: OUTPUT C PRGMMR: WOOLLEN ORG: NP20 DATE: 1990-11-06 C C ABSTRACT: EXAMINES OUTPUTS FROM OI QUALITY CONTROL INTERPOLATION AND C MAKES DECISION WHETHER AN OB PASSES OR FAILS THE CHECKS. ALSO C PROVIDES FOR A DETAILED PRINT TRACING ALL DECISIONS MADE IN REGARD C TO SPECIFIED OBSERVATIONS IN THE QC PROCESS. C C PROGRAM HISTORY LOG: C 1990-11-06 J. WOOLLEN - ORIGINAL AUTHOR C C USAGE: CALL OUTPUT C OUTPUT FILES: C UNIT 06 - STANDARD OUTPUT PRINT C UNIT 71 - DETAILED PRINTOUTS FOR QC DIAGNOSIS (IPOINT IN UNIT 18) C UNIT 72 - DETAILED PRINTOUTS FOR QC DIAGNOSIS (IPOINT IN UNIT 18) C UNIT 73 - DETAILED PRINTOUTS FOR QC DIAGNOSIS (IPOINT IN UNIT 18) C UNIT 74 - DETAILED PRINTOUTS FOR QC DIAGNOSIS (IPOINT IN UNIT 18) C UNIT 75 - DETAILED PRINTOUTS FOR QC DIAGNOSIS (IPOINT IN UNIT 18) C UNIT 76 - DETAILED PRINTOUTS FOR QC DIAGNOSIS (IPOINT IN UNIT 18) C UNIT 77 - DETAILED PRINTOUTS FOR QC DIAGNOSIS (IPOINT IN UNIT 18) C UNIT 78 - DETAILED PRINTOUTS FOR QC DIAGNOSIS (IPOINT IN UNIT 18) C UNIT 79 - DETAILED PRINTOUTS FOR QC DIAGNOSIS (IPOINT IN UNIT 18) C UNIT 80 - DETAILED PRINTOUTS FOR QC DIAGNOSIS (IPOINT IN UNIT 18) C UNIT ITUNIT - DETAILED TOSSLIST (SEE REMARKS) C C REMARKS: POSSIBLE VALUES FOR ITUNIT (PASSED IN VIA COMMMON C /TOSSUNIT/): C 61 - UNIT FOR FILE HOLDING DETAILED TOSSLIST FOR SURFACE HEIGHTS C 62 - UNIT FOR FILE HOLDING DETAILED TOSSLIST FOR TEMPS AND WINDS C 63 - UNIT FOR FILE HOLDING DETAILED TOSSLIST FOR SAT TEMPS C 64 - UNIT FOR FILE HOLDING DETAILED TOSSLIST FOR SSMI WIND SPEEDS C C ATTRIBUTES: C LANGUAGE: FORTRAN 90 C MACHINE: NCEP WCOSS C C$$$ SUBROUTINE OUTPUT PARAMETER (MAXREP=1000000) PARAMETER (MAXLEV=2000000) COMMON /OIQCDATA/ FLV (MAXREP ) , FLT (MAXREP ) , .FLN (MAXREP ) , FTM (MAXREP ) , FKX (MAXREP ) , .FNL (MAXREP ) , FMS (MAXREP ) , FMW (MAXREP,2 ) , .FQL (MAXREP ) , FHD (MAXLEV ) , FPRC(MAXLEV ) , .FPFD(MAXLEV ) , FPLV(MAXLEV ) , FPR (MAXLEV ) , .FQM (MAXLEV ) , FOI (MAXLEV ) , FOB (MAXLEV,2 ) , .FOE (MAXLEV ) , FFE (MAXLEV,2 ) , FTZ (MAXLEV ) , .FTW (MAXLEV,10,2) , FPF (MAXLEV,5,2 ) , FOC (MAXREP ) CHARACTER*3 FOC COMMON /HVECT / NXXYY,MCHK(64,2),MOBS(64,10,4,2), . 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) COMMON /TOLLER / ITER,TB,TZ,TSAT,TWND,TJET COMMON /OBPRT / NPT,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 CHARACTER*1 VAR2(2),VAR3(3) LOGICAL KEEP REAL*8 SID,CKSTA DATA CHECK/'Z MATRIX','U MATRIX','V MATRIX','VERTICAL','MULTIVAR'/ DATA VAR2 / 'Z' , 'W' / DATA VAR3 / 'Z' , 'U' , 'V' / DATA IOMF / 2 , 1 , 2 / CHARACTER*1 PLUMIN DATA PI180/.017453/ C---------------------------------------------------------------------- C---------------------------------------------------------------------- D = 0. DU = 0. DV = 0. TWT = 0. TWU = 0. TWV = 0. C LOOP OVER EACH OB CHECKED IN THIS GROUP C --------------------------------------- DO 50 N=1,NXXYY I = MCHK(N,1) J = MCHK(N,2) II = FHD(I) PRES = FPR(I) OBER = FOE(I) CKSTA = SID(II) CLON = FLN(II) CLAT = FLT(II) KX = FKX(II) TAG = FOC(II) C LOOK FOR A KEEP FLAG C -------------------- IQM = FQM(I) IF(IQM.EQ.0) THEN KEEP = .TRUE. CKEEP = ' KEEP' ELSE KEEP = .FALSE. CKEEP = ' ' ENDIF C LOOP OVER THE VARIABLE(S) IN EACH OB C ------------------------------------ DO 45 K=J,2*J-1 KK = K-J+1 OMF = FOB(I,IOMF(K)) AMFSUM = 0. SUMWTS = 0. QC = 0. C LOOP OVER THE OUTCOME OF EACH TYPE OF CHECK (MICRO DECISION MAKER) C ------------------------------------------------------------------ DO 20 L=1,4 AER = AE(N,L,KK) FER = FE(N,K) AMF = AA(N,L,KK) DRES(L) = OMF-AMF D(L) = DRES(L)/SQRT(OBER**2+AER**2) AWT = 1.-(AER/FER)**2 TWT(L) = MAX(.01,AWT) QC = QC + TWT(L)*D(L)**2 SUMWTS = SUMWTS + TWT(L) IF(L.LT.4) AMFSUM = AMFSUM + TWT(L)*AMF 20 CONTINUE AMFSUM = AMFSUM/(SUMWTS-TWT(4)) QC = SQRT(QC/SUMWTS) TWT = TWT/SUMWTS C FOR WINDS CHECK THE VECTOR SUM OF AGAINST THE REFERENCE C ------------------------------------------------------- IF(K.EQ.1) THEN OMA = OMF - AMFSUM ELSE IF(K.EQ.2) THEN DU = D TWU = TWT QCU = QC OMAU = OMF - AMFSUM OMFU = OMF GOTO 45 ELSE IF(K.EQ.3) THEN DV = D TWV = TWT QCV = QC OMAV = OMF - AMFSUM OMFV = OMF D = SQRT(DU**2+DV**2) TWT = SQRT(TWU**2+TWV**2) QC = SQRT(QCU**2+QCV**2) OMA = SQRT(OMAU**2+OMAV**2) OMF = SQRT(OMFU**2+OMFV**2) ELSE PRINT *, 'OUTPUT - K<>1,2,OR,3' CALL W3TAGE('PREPOBS_OIQCBUFR') CALL ERREXIT(99) ENDIF ENDIF C ADJUST THE TOLERANCE BY LOCATION AND VARIABLE C --------------------------------------------- IF(K.EQ.1) THEN SINLAT = (SIN(ABS(CLAT)*PI180)) TT = TB*TZ*SINLAT TMIN = 2.0 ELSE GNP = MIN(ABS(CLAT/30.),1.) IF(PRES.LE.400. .AND. PRES.GE.150.) THEN TT = TB*TJET*GNP ELSE TT = TB*TWND*GNP ENDIF TMIN = 3. ENDIF TT = MAX(TMIN,TT) IF(TAG.EQ.'SMI') TT = TT/SQRT(2.) C COMPUTE THE MULTIVARIATE REDUCTION RATIO AND ADJUST THE SUM C ----------------------------------------------------------- IF(OMF.NE.0.) THEN QCPC = MIN(MAX(ABS(OMA/OMF),.1),1.) TTT = TT/QCPC ELSE TTT = TT ENDIF C COMPUTE THE NORMALIZED DEVIATION FROM THE FIRST GUESS C ----------------------------------------------------- SDOMF = SQRT(OMF**2/(OBER**2+FER**2)) C MAKE A MESO DECISION FOR THIS ITERATION C --------------------------------------- IF(QC.GT.TTT) THEN KEEP = KEEP .AND. QC.LT.25. IF(.NOT.KEEP) THEN FOI(I) = -1 PLUMIN = '-' ELSE FOI(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 FOI(I) = 1 PLUMIN = '+' ENDIF C THIS WRITES DETAILED INFO ABOUT AN OB C ------------------------------------- DO 41 NN=1,NPT LU = 70 + NN pdif = abs(pres-prp(nn)) IF(CKSTA.EQ.SIDP(NN).AND.pdif.lt..1.AND.KX.EQ.KXP(NN)) THEN WRITE(LU,100)CKSTA,KX,CLON,CLAT,PRES,VAR2(J),OMF,SDOMF, . (D(L),TWT(L),L=1,4),QC,TTT,CKEEP,PLUMIN DO 44 KK=J,J*2-1 WRITE(LU,*) 'ITERATION=',ITER,' CHECK ON ',VAR3(KK) IF(KK.EQ.2) THEN SDOMF = SQRT(OMFU**2/(OBER**2+FER**2)) IF(TT.EQ.TTT) TTTT = TT IF(TT.NE.TTT) TTTT = TT/ABS(OMAU/OMFU) WRITE(LU,100)CKSTA,KX,CLON,CLAT,PRES,VAR3(KK),OMFU,SDOMF, . (DU(L),TWU(L),L=1,4),QCU,TTTT,CKEEP ELSE IF(KK.EQ.3) THEN SDOMF = SQRT(OMFV**2/(OBER**2+FER**2)) IF(TT.EQ.TTT) TTTT = TT IF(TT.NE.TTT) TTTT = TT/ABS(OMAV/OMFV) WRITE(LU,100)CKSTA,KX,CLON,CLAT,PRES,VAR3(KK),OMFV,SDOMF, . (DV(L),TWV(L),L=1,4),QCV,TTTT,CKEEP ENDIF DO 44 L=1,4 WRITE(LU,*)CHECK(L),' FOR CHECK ON ',VAR3(KK),' DATUM ' SOINC = 0. DO 43 IOB=1,MDIM(N,L) ILEV = MOBS(N,IOB,L,1) KIV = MOBS(N,IOB,L,2) KI2 = MAX(KIV-1,1) IREP = FHD(ILEV) X = FLN(IREP) Y = FLT(IREP) 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(A8,I4,2F7.2,F8.2,A4,2F8.2,4(' *',F6.1,' (',F5.2,')' ), . ' *',2F6.2,A6,A1) 101 FORMAT(A8,I4,2F7.2,F8.2,A4,2F8.2,4(' *',F6.1,2X,I2),' *',2F6.2,A6, . A1) 102 FORMAT(A8,I4,2F7.2,F8.2,A4,6F8.2) RETURN END C$$$ SUBPROGRAM DOCUMENTATION BLOCK C C SUBPROGRAM: PILNLNP C PRGMMR: WOOLLEN ORG: NP20 DATE: 1992-07-29 C C ABSTRACT: FUNCTION WHICH GIVEN AN PROFILE OF DESCENDING PRESSURES AND C A PROFILE OF QUANTITIES VALID AT THOSE PRESSURES, RETURNS AN C INTERPOLATED QUANTITY VALID AT A GIVEN ARBITRARY PRESSURE. C C PROGRAM HISTORY LOG: C 1986-03-21 G. DIMEGO - ORINIGAL AUTHOR C 1988-11-24 D. DEAVEN - RECODED FOR CYBER 205 C 1992-07-29 J. WOOLLEN - UPDATED FOR OIQC C C USAGE: CALL = PILNLNP(P,PARAY,QARAY,KMAX) C INPUT ARGUMENT LIST: C P - PRESSURE TO INTERPOLATE TO C PARAY - GIVEN PRESSURE PROFILE C QARAY - QUANTITIES VALID FOR PRESSURE PROFILE C KMAX - LENGTH OF PROFILE C C FUNCTION RETURN VALUE C PILNLNP - INTERPOLATED QUANTITY C C ATTRIBUTES: C LANGUAGE: FORTRAN 90 C MACHINE: NCEP WCOSS C C$$$ FUNCTION PILNLNP(P,PARAY,QARAY,KMAX) DIMENSION PARAY(KMAX),QARAY(KMAX) DO 10 LA=1,KMAX IF(P.GE.PARAY(LA)) GOTO 11 10 CONTINUE 11 IF(LA.GT.KMAX) THEN LA = KMAX LB = KMAX ELSE IF(LA.EQ.1) THEN LA = 1 LB = 1 ELSE LB = LA-1 ENDIF PA = PARAY(LA) PB = PARAY(LB) IF(PA.NE.PB) THEN WK = LOG(P/PB) / LOG(PA/PB) ELSE WK = 0. ENDIF PILNLNP = QARAY(LB) + (QARAY(LA)-QARAY(LB)) * WK RETURN END C$$$ SUBPROGRAM DOCUMENTATION BLOCK C C SUBPROGRAM: PREPQM C PRGMMR: WOOLLEN ORG: NP20 DATE: 2013-03-26 C C ABSTRACT: REWRITES THE PREPBUFR FILE WITH UPDATED QUALITY MARKS. C C PROGRAM HISTORY LOG: C 1990-11-06 J. WOOLLEN - ORIGINAL AUTHOR C 2013-03-26 J. WOOLLEN INITIAL PORT TO WCOSS: FIX ADDRESSING C ALIGNMENT PROBLEM WHEN THE LAST BUFR MESSAGE CONTAINS NO DATA C WHICH IS OIQC PROCESSED C C USAGE: CALL PREPQM(LUBIN,LUBOT) C INPUT ARGUMENT LIST: C LUBIN - LOGICAL UNIT FOR OBSERVATION (PREPBUFR) PRIOR TO C UPDATED QUALITY MARKS C LUBOT - LOGICAL UNIT FOR OBSERVATION (PREPBUFR) WITH UPDATED C QUALITY MARKS C C INPUT FILES: C UNIT LUBIN - OBSERVATION (PREPBUFR) PRIOR TO UPDATED QUALITY C MARKS C C OUTPUT FILES: C UNIT 06 - STANDARD OUTPUT PRINT C UNIT LUBOT - OBSERVATION (PREPBUFR) FILE WITH UPDATED QUALITY C MARKS C C REMARKS: NONE. C C C ATTRIBUTES: C LANGUAGE: FORTRAN 90 C MACHINE: NCEP WCOSS C C$$$ SUBROUTINE PREPQM(LUBIN,LUBOT) PARAMETER (MAXREP=1000000) PARAMETER (MAXLEV=2000000) COMMON /OIQCDATA/ FLV (MAXREP ) , FLT (MAXREP ) , .FLN (MAXREP ) , FTM (MAXREP ) , FKX (MAXREP ) , .FNL (MAXREP ) , FMS (MAXREP ) , FMW (MAXREP,2 ) , .FQL (MAXREP ) , FHD (MAXLEV ) , FPRC(MAXLEV ) , .FPFD(MAXLEV ) , FPLV(MAXLEV ) , FPR (MAXLEV ) , .FQM (MAXLEV ) , FOI (MAXLEV ) , FOB (MAXLEV,2 ) , .FOE (MAXLEV ) , FFE (MAXLEV,2 ) , FTZ (MAXLEV ) , .FTW (MAXLEV,10,2) , FPF (MAXLEV,5,2 ) , FOC (MAXREP ) CHARACTER*3 FOC COMMON /OBS/ NREP,NLEV,SID(MAXREP),INOB(MAXREP),IMAP(360,181) COMMON /DATAZZ/ KLEV,HDR(5), . OBS(8,255),QMS(6,255),PCS(6,255),RCS(6,255), . OBE(8,255),QME(6,255),OES(6,255),FES(6,255), . FCS(7,255) CHARACTER*30 OSTR,QSTR,PSTR,RSTR CHARACTER*8 SUBSET CHARACTER*3 TAG LOGICAL EVENT1,EVENTN,SMI,SFC REAL*8 OBS,QMS,OBE,QME,PCS,RCS,SID DATA OSTR /'POB TOB UOB VOB QOB'/ DATA QSTR /'PQM TQM WQM NUL QQM'/ DATA PSTR /'PPC TPC WPC NUL QPC'/ DATA RSTR /'PRC TRC WRC NUL QRC'/ DATA VMAX /10E10/ C----------------------------------------------------------------------- C----------------------------------------------------------------------- PRINT* PRINT*,'******WRITE THE PREPBUFR FILE*******' PRINT* C OPEN THE INPUT AND OUTPUT FILES C ------------------------------- C CALL OPENBF(LUBIN,'IN',LUBIN) CALL UFBQCD(LUBIN,'OIQC',QCD) CALL MAXOUT(25000) CALL OPENBF(LUBOT,'OUT',LUBIN) imsg = 1 irec = 1 ifld = 1 C LOOP OVER THE OBS CHECKED BY OIQC C --------------------------------- DO 100 ILEV=1,NLEV IREP = FHD(ILEV) OBM = FOI(ILEV) TAG = FOC(IREP) KX = FKX(IREP) QM = FQM(ILEV) SFC = TAG.EQ.'SFC' SMI = TAG.EQ.'SMI' JREC = FPRC(ILEV) JFLD = FPFD(ILEV) JLEV = FPLV(ILEV) IF(SFC.AND.FLV(IREP).NE.ILEV) THEN JUPD = 1 ELSE JUPD = KX/100+1 ENDIF C SEE IF THERE IS AN EVENT TO WRITE AND WHICH KIND IT WOULD BE C ------------------------------------------------------------ QMN = QM IF(OBM.EQ.-1 .AND. QM.GE.0 .AND. QM.LE.3) QMN = QM + 4 IF(OBM.EQ. 1 .AND. QM.GE.4 .AND. QM.LE.7) QMN = QM - 4 EVENT1 = QMN.NE.QM .OR. SMI EVENTN = EVENT1 .AND. IREC.EQ.JREC .AND. IFLD.EQ.JFLD EVENT1 = EVENT1 .AND. .NOT.EVENTN C----------------------------------------------------------------------- C DECIDE HOW TO PROCEED FOR THIS OB/EVENT: C IF THE FIRST OB, SYNCHRONIZE THE FILES; C IF A SUBSEQUENT OB: C FOR THE FIRST EVENT IN A REPORT, POSITION THE OUTPUT FILE C FOR SUBSEQUENT EVENTS IN A REPORT, SKIP THE POSITIONING C FOR NO EVENT AT ALL, CONTINUE THE LOOP C----------------------------------------------------------------------- IF(ILEV.EQ.1) GOTO 49 IF(EVENT1 ) GOTO 10 IF(EVENTN ) GOTO 50 IF(ILEV.LT.NLEV) CYCLE C FLUSH A COMPLETED SET OF EVENTS BEFORE POSITIONING THE OUTPUT FILE C ------------------------------------------------------------------ 10 CONTINUE CALL UFBINT(LUBOT,OBE,8,KLEV,IRET,OSTR) CALL UFBINT(LUBOT,QME,6,KLEV,IRET,QSTR) CALL UFBINT(LUBOT,PCS,6,KLEV,IRET,PSTR) CALL UFBINT(LUBOT,RCS,6,KLEV,IRET,RSTR) CALL WRITSB(LUBOT) C WIND AND COPY OBS UP TO THE NEXT EVENTED MESSAGE AND SUBSET C ----------------------------------------------------------- IF(IREC.LT.JREC.OR.ILEV.EQ.NLEV) THEN DO WHILE(ICOPYSB(LUBIN,LUBOT).EQ.0) ENDDO IF(IREADMM(imsg,SUBSET,IDATE).ne.0) THEN IF(ilev.lt.nlev) GOTO 900 goto 100 endif CALL UFBCNT(LUBIN,IREC,IFLD) DO WHILE(IREC.LT.JREC.OR.ILEV.EQ.NLEV) CALL CLOSMG(LUBOT) CALL COPYMG(LUBIN,LUBOT) IF(IREADMM(imsg ,SUBSET,IDATE).NE.0) GOTO 100 CALL UFBCNT(LUBIN,IREC,IFLD) ENDDO ENDIF DO WHILE(IFLD.LT.JFLD-1) CALL OPENMB(LUBOT,SUBSET,IDATE) IF(ICOPYSB(LUBIN,LUBOT).NE.0) GOTO 901 CALL UFBCNT(LUBIN,IREC,IFLD) ENDDO C READ THE NEXT SUBSET TO RECEIVE EVENTS C -------------------------------------- 49 CONTINUE IF(IREADSB(LUBIN).NE.0) then IF(IREADMM(imsg,SUBSET,IDATE).NE.0) GOTO 100 goto 49 endif CALL UFBINT(LUBIN,OBS,8,255,KLEV,OSTR) CALL UFBINT(LUBIN,QMS,6,255,KLEV,QSTR) CALL OPENMB(LUBOT,SUBSET,IDATE) CALL UFBCNT(LUBIN,IREC,IFLD) CALL UFBCPY(LUBIN,LUBOT) OBE = VMAX QME = VMAX PCS = VMAX RCS = VMAX 50 CONTINUE C ADD AN EVENT IF THE OIQC QM <> DATA QM C -------------------------------------- IF(QMS(JUPD,JLEV).NE.QMN) THEN IF(JUPD.LE.3) OBE(JUPD,JLEV) = OBS(JUPD,JLEV) IF(JUPD.EQ.3) OBE( 4,JLEV) = OBS( 4,JLEV) QME(JUPD,JLEV) = QMN PCS(JUPD,JLEV) = QCD RCS(JUPD,JLEV) = QMN ENDIF C LINK THE MOISTURE QM TO THE PRESSURE/TEMPERATURE QM C --------------------------------------------------- IF(QMN.GE.4 .AND. QMS(5,JLEV).LT.4) THEN OBE(5,JLEV) = OBS(5,JLEV) QME(5,JLEV) = QMS(5,JLEV)+4 PCS(5,JLEV) = QCD RCS(5,JLEV) = QMS(5,JLEV)+4 ENDIF C ADD AN EVENT CONTAINING ANALYSED SSMI WIND DIRECTION C ---------------------------------------------------- IF(SMI) THEN OBE(3,JLEV) = FMW(IREP,1) + FOB(ILEV,1) OBE(4,JLEV) = FMW(IREP,2) + FOB(ILEV,2) QME(3,JLEV) = QM PCS(3,JLEV) = QCD RCS(3,JLEV) = QM ENDIF C TRAP THE LAST OIQC RESULT IN ORDER TO COMPLETE THE OUTPUT FILE C -------------------------------------------------------------- 100 CONTINUE C CLOSE OUT THE FILES AND EXIT C ---------------------------- call ufbcnt(lubot,ireo,iflo) PRINT*,'PREPQM READ ',IREC,' WROTE ',IREO,' PREPBUFR RECORDS' CALL CLOSBF(LUBIN) CALL CLOSBF(LUBOT) RETURN 900 PRINT *, 'PREPQM - READMG ALIGNMENT PROBLEM' CALL W3TAGE('PREPOBS_OIQCBUFR') CALL ERREXIT(99) 901 PRINT *, 'PREPQM - COPYSB ALIGNMENT PROBLEM' CALL W3TAGE('PREPOBS_OIQCBUFR') CALL ERREXIT(99) 902 PRINT *, 'PREPQM - READMM ALIGNMENT PROBLEM' CALL W3TAGE('PREPOBS_OIQCBUFR') CALL ERREXIT(99) END C$$$ SUBPROGRAM DOCUMENTATION BLOCK C C SUBPROGRAM: QCLOOP C PRGMMR: WOOLLEN ORG: NP20 DATE: 1990-11-06 C C ABSTRACT: ITERATES THE QC DECSION MAKING PROCESS UNTIL ALL DECISIONS C IN TWO CONSEQUETIVE ITERATIONS ARE THE SAME. AFTER FOUR ITERATIONS C ONLY REPORTS THAT FAIL THE CHECKS ARE CHECKED FURTHER. C C PROGRAM HISTORY LOG: C 1990-11-06 J. WOOLLEN - ORIGINAL AUTHOR C C USAGE: CALL QCLOOP(QCPTS,MAXTRY) C INPUT ARGUMENT LIST: C QCPTS - NAME OF SUBROUTINE TO CALL C MAXTRY - MAXIMUM NUMBER OF TRIES C C INPUT FILES: C UNIT 20 - TOLERANCE LEVELS (FOR EXPERIMENTATION) C C OUTPUT FILES: C UNIT 06 - STANDARD OUTPUT PRINT C C REMARKS: NONE. C C ATTRIBUTES: C LANGUAGE: FORTRAN 90 C MACHINE: NCEP WCOSS C C$$$ SUBROUTINE QCLOOP(QCPTS,MAXTRY) PARAMETER (MAXREP=1000000) PARAMETER (MAXLEV=2000000) COMMON /OBS/ NREP,NLEV,SID(MAXREP),INOB(MAXREP),IMAP(360,181) COMMON /OIQCDATA/ FLV (MAXREP ) , FLT (MAXREP ) , .FLN (MAXREP ) , FTM (MAXREP ) , FKX (MAXREP ) , .FNL (MAXREP ) , FMS (MAXREP ) , FMW (MAXREP,2 ) , .FQL (MAXREP ) , FHD (MAXLEV ) , FPRC(MAXLEV ) , .FPFD(MAXLEV ) , FPLV(MAXLEV ) , FPR (MAXLEV ) , .FQM (MAXLEV ) , FOI (MAXLEV ) , FOB (MAXLEV,2 ) , .FOE (MAXLEV ) , FFE (MAXLEV,2 ) , FTZ (MAXLEV ) , .FTW (MAXLEV,10,2) , FPF (MAXLEV,5,2 ) , FOC (MAXREP ) CHARACTER*3 FOC COMMON /FPLIST / NFL,INFL(MAXLEV) COMMON /CKLIST / NDD,INDD(MAXLEV) COMMON /TOLLER / ITER,T,TZ,TSAT,TWND,TJET COMMON /TOSSUNIT/ ITUNIT NAMELIST /TOLLS / TB,TZ,TSAT,TWND,TJET LOGICAL DONE include "mpif.h" dimension istat(mpi_status_size) dimension ireq(100) C---------------------------------------------------------------------- call mpi_comm_rank(mpi_comm_world,myid,ierr) call mpi_comm_size(mpi_comm_world,nprc,ierr) mdtp = mpi_real mdes = 0 mcom = mpi_comm_world C----------------------------------------------------------------------- IF(NDD.EQ.0) THEN PRINT*,'QCLOOP - NDD = 0' RETURN ENDIF ITRY = 1 ITER = 1 IFLIP = 0 C SET UP THE TOLERANCES C --------------------- TB = 4. TZ = 5.0 TSAT = 5.0 TWND = 4.5 TJET = 5.5 REWIND (20) READ(20,TOLLS,END=1,err=1) 1 if(myid.eq.0) then C WRITE(6,TOLLS) endif TZ = TZ/TB TSAT = TSAT/TB TWND = TWND/TB TJET = TJET/TB T = TB C DO THE QCOI CHECKS ON LIST IN COMMON /CKLIST/ C --------------------------------------------- 10 continue C REWIND(ITUNIT) CALL QCPTS call syncndd(foi,nlev,mpi_real) C SET OB ERRS TO REFLECT LAST ROUND OF CHECKS - RECORD FLIPS C ---------------------------------------------------------- NBB = 0 NCK = 0 NFG = 0 NFB = 0 NFL = 0 DO 20 N=1,NDD ILEV = INDD(N) OBER = FOE(ILEV) IF(OBER.NE.0) THEN IF(FOI(ILEV).EQ.1) THEN FOE(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 FOE(ILEV) = -ABS(OBER) IF(OBER.GT.0) THEN NFB = NFB + 1 NFL = NFL + 1 INFL(NFL) = ILEV ENDIF NBB = NBB + 1 ENDIF NCK = NCK + 1 ENDIF 20 CONTINUE C CHECK STATUS AND WRITE ITERATION REPORT C --------------------------------------- DONE = NFL.EQ.0 .OR. ITRY.EQ.MAXTRY if(myid.eq.0) then WRITE(6 ,30)ITRY,NFG,NFB,NFL,NBB,NCK endif 30 FORMAT('ITERATION ',I2,' : NFLIP=',3I6,' NBAD=',I9,' NTOT=',I9) C FINISHED YET? C ------------- IF(.NOT.DONE) THEN ITRY = ITRY+1 IF(ITRY.LT.MAXTRY) THEN ITER = ITRY GOTO 10 ELSE ITER = ITRY DO 40 N=1,NFL FOE(INFL(N)) = -ABS(FOE(INFL(N))) 40 CONTINUE GOTO 10 ENDIF ENDIF C FINAL REPORT AND AWAY C --------------------- NBAD = NBB if(myid.eq.0) then WRITE(6 ,*)' TOTAL BAD OBS=',NBAD WRITE(6 ,*) endif ITRY = 0 ITER = 0 RETURN END C$$$ SUBPROGRAM DOCUMENTATION BLOCK C C SUBPROGRAM: QCOI C PRGMMR: WOOLLEN ORG: NP20 DATE: 1990-03-24 C C ABSTRACT: COMPUTES AND SETS UP THE COEFFICIENT MATRIX AND RIGHT HAND C SIDES. SOLVES THE MATRIX PROBLEMS. APPLIES THE COMPUTED WEIGHTS TO C THE RESIDUALS TO OBTAIN THE ANALYSIS INCREMENTS AND ANALYSIS C ERRORS. C C PROGRAM HISTORY LOG: C 1988-02-05 DEAVEN - ORIGINAL AUTHOR C 1990-03-24 J. WOOLLEN - UPDATED FOR OIQC C C USAGE: CALL QCOI C C REMARKS: INPUT INFORMATION IN COMMON /HVECT/: C NXXYY - NUMBER OF OBS IN PACKAGE C XX - LONGITUDES OF OBS IN PACKAGE C YY - LATITUDES OF OBS IN PACKAGE C PP - PRESSURES OF OBS IN PACKAGE C II - INDEXES OF OBS IN PACKAGE C LL - REPORT LEVEL OF OBS IN PACKAGE C FESD - FORECAST ERRORS AT OB LOCATIONS C MKOBS - LIST OF OB INDEXES USED IN EACH INTERPOLATION C NDIM - NUMBER OF OBS USED IN EACH INTERPOLATION C JVMAX - NUMBER OF INTERPOLATING OBS FOR EACH CHECK IN C EACH PACKAGE C C OUTPUT INFORMATION IN COMMON /HVECT/: C AA - INTERPOLATED VALUE OF OBS IN PACKAGE C AESD - ANALYSIS ERRORS AT OB LOCATIONS C C ATTRIBUTES: C LANGUAGE: FORTRAN 90 C MACHINE: NCEP WCOSS C C$$$ SUBROUTINE QCOI PARAMETER (MAXREP=1000000) PARAMETER (MAXLEV=2000000) COMMON /OIQCDATA/ FLV (MAXREP ) , FLT (MAXREP ) , .FLN (MAXREP ) , FTM (MAXREP ) , FKX (MAXREP ) , .FNL (MAXREP ) , FMS (MAXREP ) , FMW (MAXREP,2 ) , .FQL (MAXREP ) , FHD (MAXLEV ) , FPRC(MAXLEV ) , .FPFD(MAXLEV ) , FPLV(MAXLEV ) , FPR (MAXLEV ) , .FQM (MAXLEV ) , FOI (MAXLEV ) , FOB (MAXLEV,2 ) , .FOE (MAXLEV ) , FFE (MAXLEV,2 ) , FTZ (MAXLEV ) , .FTW (MAXLEV,10,2) , FPF (MAXLEV,5,2 ) , FOC (MAXREP ) CHARACTER*3 FOC 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) COMMON /HVECT / NXXYY,MCHK(64,2),MOBS(64,10,4,2), . MDIM(64,4),AA(64,4,2), . AE(64,4,2),WT(64,10,4,2), . FE(64,3) DIMENSION CHLS(LVEC) , XXCP(LVEC) , XXOB(LVEC) , C1(LMAT) , . FLAT(LVEC) , YYCP(LVEC) , YYOB(LVEC) , C2(LMAT) , . TZFA(LVEC) , KKCP(LVEC) , KKOB(LVEC) , DD(LMAT) , . NGRP(LVEC) , PPCP(LVEC) , PPOB(LVEC) , XI(LMAT) , . NIOB(LVEC) , OEOB(LVEC) , YI(LMAT) , . OBOB(LVEC) , KI(LMAT) , . PI(LMAT) , . XJ(LMAT) , . YJ(LMAT) , . KJ(LMAT) , . PJ(LMAT) , . CH(LMAT) , . FL(LMAT) , . TZ(LMAT) , . NG(LMAT) , . NI(LMAT) , . NJ(LMAT) DIMENSION FA(NMX,NN), RA(NMX,N,2), DP(NMX,2), IVAR(3) DATA IVAR/2,1,2/ C----------------------------------------------------------------------- C----------------------------------------------------------------------- IF(NXXYY.EQ.0) RETURN AA = 0. AE = 0. WT = 0. C LOOP THRU EACH TYPE OF QC CHECK C ------------------------------- DO 100 IQ=1,NQC IVMAX = 0 C CLEAR THE MATRIX ARRAYS C ----------------------- FA = 0. RA = 0. DP = 0. C SET UP VECTORS WITH INFORMATION ABOUT EACH OB TO BE CHECKED C ----------------------------------------------------------- DO 15 IGRP=1,NXXYY NDIM = MDIM(IGRP,IQ) IF(NDIM.EQ.0) GOTO 15 ILEV = MCHK(IGRP,1) IREP = FHD(ILEV) TZF = FTZ(ILEV) XXC = FLN(IREP) YYC = FLT(IREP) PPC = FPR(ILEV) KKC = MCHK(IGRP,2) CHL = CHLP(INT(MIN(PPC,1500.))) FLA = FFLAT(INT(ABS(YYC)+1.51)) PPC = LOG(PPC) DO 10 IOB=1,NDIM JLEV = MOBS(IGRP,IOB,IQ,1) JTYP = MOBS(IGRP,IOB,IQ,2) JREP = FHD(JLEV) IVMAX = IVMAX+1 CHLS(IVMAX) = CHL FLAT(IVMAX) = FLA TZFA(IVMAX) = TZF XXCP(IVMAX) = XXC YYCP(IVMAX) = YYC KKCP(IVMAX) = KKC PPCP(IVMAX) = PPC XXOB(IVMAX) = FLN(JREP) YYOB(IVMAX) = FLT(JREP) PPOB(IVMAX) = LOG(FPR(JLEV)) KKOB(IVMAX) = JTYP OEOB(IVMAX) = FOE(JLEV) OBOB(IVMAX) = FOB(JLEV,IVAR(JTYP)) NGRP(IVMAX) = IGRP NIOB(IVMAX) = IOB 10 CONTINUE 15 CONTINUE IF(IVMAX.EQ.0) GOTO 61 C COMPUTE THE RIGHT HAND SIDES C ---------------------------- CALL COOR2(IVMAX,CHLS,FLAT,TZFA,XXOB,YYOB,KKOB,XXCP,YYCP,KKCP, . DD,C1,C2) DO 20 I=1,IVMAX IGRP = NGRP(I) IOB = NIOB(I) VCOR = 1./(1.+CVC*(PPOB(I)-PPCP(I))**2) RA(IGRP,IOB,1) = C1(I)*VCOR RA(IGRP,IOB,2) = C2(I)*VCOR 20 CONTINUE C COMPUTE THE OFF-DIAGONAL MATRIX TERMS C ------------------------------------- K = 0 L = 0 DO 40 IGRP=1,NXXYY NDIM = MDIM(IGRP,IQ) DO 30 IOB = 1,NDIM-1 DO 30 JOB = IOB+1,NDIM I = IOB+K J = JOB+K L = L+1 XI(L) = XXOB(I) YI(L) = YYOB(I) KI(L) = KKOB(I) PI(L) = PPOB(I) XJ(L) = XXOB(J) YJ(L) = YYOB(J) KJ(L) = KKOB(J) PJ(L) = PPOB(J) CH(L) = CHLS(I) FL(L) = FLAT(I) TZ(L) = TZFA(I) NG(L) = NGRP(I) NI(L) = NIOB(I) NJ(L) = NIOB(J) 30 CONTINUE K = K+NDIM 40 CONTINUE CALL COOR1(L,CH,FL,TZ,XI,YI,KI,XJ,YJ,KJ,DD,C1,C2) C STORE OFF-DIAGONAL MATRIX TERMS IN THE SOLVER ARRAY C --------------------------------------------------- DO 50 I=1,L IGRP = NG(I) IOB = IFA(NI(I),NJ(I)) FA(IGRP,IOB) = C1(I)/(1.+CVC*(PI(I)-PJ(I))**2) 50 CONTINUE C COMPUTE AND STORE THE ON-DIAGONAL MATRIX TERMS C ---------------------------------------------- DO 55 I=1,IVMAX IGRP = NGRP(I) IOB = IFA(NIOB(I),NIOB(I)) FA(IGRP,IOB) = 1. + (OEOB(I)/FE(IGRP,KKOB(I)))**2 55 CONTINUE C SOLVE THE LINEAR INTERPOLATION EQUATIONS C ---------------------------------------- CALL DRCTSL(FA,RA,DP,MDIM(1,IQ),NXXYY) C COMPUTE AND STORE ANALYSIS INCREMENTS AND ERRORS C ------------------------------------------------ DO 60 I=1,IVMAX IGRP = NGRP(I) IOB = NIOB(I) FERAT = FE(IGRP,KKCP(I))/FE(IGRP,KKOB(I)) WT(IGRP,IOB,IQ,1) = RA(IGRP,IOB,1) * FERAT WT(IGRP,IOB,IQ,2) = RA(IGRP,IOB,2) * FERAT AA(IGRP,IQ,1) = AA(IGRP,IQ,1) + OBOB(I)*WT(IGRP,IOB,IQ,1) AA(IGRP,IQ,2) = AA(IGRP,IQ,2) + OBOB(I)*WT(IGRP,IOB,IQ,2) 60 CONTINUE 61 DO 65 I=1,NXXYY AE(I,IQ,1) = FE(I,MCHK(I,2)) * SQRT(1.-DP(I,1)) AE(I,IQ,2) = FE(I,MCHK(I,2)) * SQRT(1.-DP(I,2)) 65 CONTINUE C END OF LOOP OVER EACH QC CHECK C ------------------------------ 100 CONTINUE RETURN END C$$$ SUBPROGRAM DOCUMENTATION BLOCK C C SUBPROGRAM: QCPTS C PRGMMR: WOOLLEN ORG: NP20 DATE: 1990-11-06 C C ABSTRACT: PUTS TOGETHER PACKAGES OF OBSERVATIONS TO BE CHECKED IN A C GIVEN ITERATION AND INVOKES AUTOTASKING OF THE DATA SELECTION, OI, C AND DECISION MAKING PROCEDURES. C C PROGRAM HISTORY LOG: C 1990-11-06 J. WOOLLEN - ORIGINAL AUTHOR C C USAGE: CALL QCPTS C C REMARKS: INPUT INFORMATION IN COMMON /CKLIST/: C NDD - NUMBER OF OBSERVATIONS TO CHECK C INDD - ARRAY OF OBSERVATION INDEXES TO CHECK C C C ATTRIBUTES: C LANGUAGE: FORTRAN 90 C MACHINE: NCEP WCOSS C C$$$ SUBROUTINE QCPTS PARAMETER (MAXLEV=2000000) COMMON /HVECT / NXXYY,MCHK(64,2),MOBS(64,10,4,2), . MDIM(64,4),AA(64,4,2), . AE(64,4,2),WTS(64,10,4,2), . FE(64,3) common /nddspt/ nd1(0:100),nd2(0:100) COMMON /CKLIST/ NDD,INDD(MAXLEV) DATA NMX /64/ include "mpif.h" C----------------------------------------------------------------------- call mpi_comm_rank(MPI_COMM_WORLD,myid,ierr) c call nddsplit(myid,nd1,nd2) C----------------------------------------------------------------------- C MULTI-PROCESS QC OB LISTS C ------------------------- DO I=nd1(myid),nd2(myid),nmx CALL SELDAT(i,min(nmx,nd2(myid)-i+1)) CALL QCOI CALL OUTPUT ENDDO RETURN END C$$$ SUBPROGRAM DOCUMENTATION BLOCK C C SUBPROGRAM: RDOBER C PRGMMR: WOOLLEN ORG: NP20 DATE: 1992-07-29 C C ABSTRACT: INITIALIZES ARRAY 'QIKE' WITH OBSERVATION ERRORS FOR ALL C KINDS OF OBSERVATION TYPES LISTED IN THE OBSERVATION ERROR FILE. C C PROGRAM HISTORY LOG: C 1986-03-21 G. DIMEGO - ORIGINAL AUTHOR C 1988-11-24 D. DEAVEN - RECODED FOR CYBER 205 C 1992-07-29 J. WOOLLEN - UPDATED FOR OIQC C C USAGE: CALL RDOBER(LU) C INPUT ARGUMENT LIST: C LU - LOGICAL UNIT FOR OBSERVATION ERROR TABLE FILE SPECIFIC C TO THIS PROGRAM C C INPUT FILES: C UNIT LU - OBSERVATION ERROR TABLE FILE SPECIFIC TO THIS PROGRAM C (THIS IS NOT THE ANALYSIS OBSERVATION ERROR FILE - C THOSE DATA ARE READ FROM THE PREPBUFR FILE) C OUTPUT FILES: C UNIT 06 - STANDARD OUTPUT PRINT C C ATTRIBUTES: C LANGUAGE: FORTRAN 90 C MACHINE: NCEP WCOSS C C$$$ SUBROUTINE RDOBER(LU) COMMON /OBERRS/ QIKE(21,300),PMAND(21),P50(21),CHL(21) COMMON /OBDESC/ DESC(300) COMMON /GETOE/ USEOE CHARACTER*80 CKX,DESC DIMENSION RAWERR(21) LOGICAL USEOE DATA RAWERR / 7.0 , 7.5 , 8.5 , 11.0 , 13.7 , 15.2 , 16.0 , . 16.1 , 17.5 , 19.1 , 20.5 , 24.0 , 24.0 , 24.0 , . 24.0 , 24.0 , 24.0 , 24.0 , 24.0 , 24.0 , 24.0 / DATA DUMMY/0./ C---------------------------------------------------------------------- C---------------------------------------------------------------------- C DEFAULT VALUES FOR OB DESCRIPTIONS AND ERRORS C --------------------------------------------- DO 1 I=1,300 DESC(I) = ' ' DO 1 L=1,21 QIKE(L,I) = 0. IF(I.EQ.120) QIKE(L,I) = RAWERR(L) 1 CONTINUE USEOE = .FALSE. C READ OB ERROR VALUES FROM OB ERROR FILE C --------------------------------------- 10 READ(LU,'(A80)',END=20) CKX IF(CKX(1:1).EQ.' ') GOTO 10 READ(CKX,'(I3)',ERR=15) KX IF(KX.GT.0 .AND. KX.LT.300) THEN READ(LU,'(7F8.2)',END=20) (QIKE(L,KX),L=1,21) DESC(KX) = CKX USEOE = .TRUE. GOTO 10 ENDIF 15 PRINT*,'******************************************************' PRINT*,'******* UNRECOGNISED RECORD IN OBERRS FILE *********' PRINT*, CKX PRINT*,'******************************************************' GOTO 10 20 RETURN END C$$$ SUBPROGRAM DOCUMENTATION BLOCK C C SUBPROGRAM: SATPTS C PRGMMR: WOOLLEN ORG: NP20 DATE: 1990-11-06 C C ABSTRACT: PUTS TOGETHER PACKAGES OF OBSERVATIONS TO BE CHECKED IN A C GIVEN ITERATION AND INVOKES AUTOTASKING OF THE DATA SELECTION, OI, C AND DECISION MAKING PROCEDURES. C C PROGRAM HISTORY LOG: C 1990-11-06 J. WOOLLEN - ORIGINAL AUTHOR C C USAGE: CALL SATPTS C OUTPUT FILES: C UNIT 06 - STANDARD OUTPUT PRINT C C REMARKS: INPUT INFORMATION IN COMMON /CKLIST/: C NDD - NUMBER OF OBSERVATIONS TO CHECK C INDD - ARRAY OF OBSERVATION INDEXES TO CHECK C C C ATTRIBUTES: C LANGUAGE: FORTRAN 90 C MACHINE: NCEP WCOSS C C$$$ SUBROUTINE SATPTS PARAMETER (MAXREP=1000000) PARAMETER (MAXLEV=2000000) COMMON /OIQCDATA/ FLV (MAXREP ) , FLT (MAXREP ) , .FLN (MAXREP ) , FTM (MAXREP ) , FKX (MAXREP ) , .FNL (MAXREP ) , FMS (MAXREP ) , FMW (MAXREP,2 ) , .FQL (MAXREP ) , FHD (MAXLEV ) , FPRC(MAXLEV ) , .FPFD(MAXLEV ) , FPLV(MAXLEV ) , FPR (MAXLEV ) , .FQM (MAXLEV ) , FOI (MAXLEV ) , FOB (MAXLEV,2 ) , .FOE (MAXLEV ) , FFE (MAXLEV,2 ) , FTZ (MAXLEV ) , .FTW (MAXLEV,10,2) , FPF (MAXLEV,5,2 ) , FOC (MAXREP ) CHARACTER*3 FOC common /nddspt/ nd1(0:100),nd2(0:100) COMMON /CKLIST/ NDD,INDD(MAXLEV) DIMENSION IND(5000),NXY(5000) DATA NMX /64/ include "mpif.h" C----------------------------------------------------------------------- call mpi_comm_rank(MPI_COMM_WORLD,myid,ierr) c call nddsplit(myid,nd1,nd2) C----------------------------------------------------------------------- C ASSEMBLE LISTS OF OBS FOR MULTI PROCESSING C ------------------------------------------ IND(1) = 1 NXY = 0 IPK = 1 NPK = 0 LREP = 0 DO 5 I=nd1(myid),nd2(myid) ILEV = INDD(I) IREP = FHD(ILEV) IF(IREP.NE.LREP) THEN LREP = IREP KLEV = FNL(IREP) IF(NXY(IPK)+KLEV.GT.NMX) THEN IF(IPK+1.GT.5000) THEN PRINT *, 'QCPTS - IPK > MAXPAK' CALL W3TAGE('PREPOBS_OIQCBUFR') CALL ERREXIT(99) ENDIF IPK = IPK+1 IND(IPK) = I NXY(IPK) = KLEV ELSE NXY(IPK) = NXY(IPK) + KLEV ENDIF ENDIF 5 CONTINUE C MULTI-PROCESS QC OB LISTS C ------------------------- DO 10 I=1,IPK CALL SELDAT(IND(I),NXY(I)) CALL QCOI CALL SATPUT 10 CONTINUE RETURN END C$$$ SUBPROGRAM DOCUMENTATION BLOCK C C SUBPROGRAM: SATPUT C PRGMMR: WOOLLEN ORG: NP20 DATE: 1990-11-06 C C ABSTRACT: EXAMINES OUTPUTS FROM OI QUALITY CONTROL INTERPOLATION AND C MAKES DECISION WHETHER A SATELLITE PASSES OR FAILS THE CHECKS. DUE C TO SPECIFIC CHARACTERISTIC DIFFERENCES OF SATELLITE RETRIEVALS FROM C CONVENTIONAL 'IN-SITU' OBSERVATIONS, A SPECIAL ALGORITHM IS C EMPLOYED IN JUDGING THE OUTCOMES OF THE QC CHECKS. ALSO PROVIDES C FOR A DETAILED PRINT TRACING ALL DECISIONS MADE IN REGARD TO C SPECIFIED OBSERVATIONS IN THE QC PROCESS. C C PROGRAM HISTORY LOG: C 1990-11-06 J. WOOLLEN - ORIGINAL AUTHOR C C USAGE: CALL SATPUT C C OUTPUT FILES: C UNIT 06 - STANDARD OUTPUT PRINT C UNIT ITUNIT - DETAILED TOSSLIST (SEE REMARKS) C C REMARKS: POSSIBLE VALUES FOR ITUNIT (PASSED IN VIA COMMMON C /TOSSUNIT/): C 61 - UNIT FOR FILE HOLDING DETAILED TOSSLIST FOR SURFACE HEIGHTS C 62 - UNIT FOR FILE HOLDING DETAILED TOSSLIST FOR TEMPS AND WINDS C 63 - UNIT FOR FILE HOLDING DETAILED TOSSLIST FOR SAT TEMPS C 64 - UNIT FOR FILE HOLDING DETAILED TOSSLIST FOR SSMI WIND SPEEDS C C INPUT INFORMATION IN COMMON /HVECT/: C NXXYY - NUMBER OF OBS IN PACKAGE C XX - LONGITUDES OF OBS IN PACKAGE C YY - LATITUDES OF OBS IN PACKAGE C PP - PRESSURES OF OBS IN PACKAGE C II - INDEXES OF OBS IN PACKAGE C LL - REPORT LEVEL OF OBS IN PACKAGE C FESD - FORECAST ERRORS AT OB LOCATIONS C AA - INTERPOLATED VALUE OF OBS IN PACKAGE C AESD - ANALYSIS ERRORS AT OB LOCATIONS C MKOBS - LIST OF OB INDEXES USED IN EACH INTERPOLATION C NDIM - NUMBER OF OBS USED IN EACH INTERPOLATION C JVMAX - NUMBER OF INTERPOLATING OBS FOR EACH CHECK IN C EACH PACKAGE C C OUTPUT INFORMATION IN COMMON /OBMARK/: C OBIND - CONTAINS QC OUTCOME FOR EACH REPORT AND LEVEL: C -2 = BELOW MODEL SURFACE (SET IN SUB PSFILE) C -1 = FAILED OI CHECKS C 0 = NO DATA IN REPORT C 1 = PASSED OI CHECKS C C ATTRIBUTES: C LANGUAGE: FORTRAN 90 C MACHINE: NCEP WCOSS C C$$$ SUBROUTINE SATPUT PARAMETER (MAXREP=1000000) PARAMETER (MAXLEV=2000000) COMMON /OIQCDATA/ FLV (MAXREP ) , FLT (MAXREP ) , .FLN (MAXREP ) , FTM (MAXREP ) , FKX (MAXREP ) , .FNL (MAXREP ) , FMS (MAXREP ) , FMW (MAXREP,2 ) , .FQL (MAXREP ) , FHD (MAXLEV ) , FPRC(MAXLEV ) , .FPFD(MAXLEV ) , FPLV(MAXLEV ) , FPR (MAXLEV ) , .FQM (MAXLEV ) , FOI (MAXLEV ) , FOB (MAXLEV,2 ) , .FOE (MAXLEV ) , FFE (MAXLEV,2 ) , FTZ (MAXLEV ) , .FTW (MAXLEV,10,2) , FPF (MAXLEV,5,2 ) , FOC (MAXREP ) CHARACTER*3 FOC COMMON /HVECT / NXXYY,MCHK(64,2),MOBS(64,10,4,2), . 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) COMMON /TOLLER / ITER,TB,TZ,TSAT,TWND,TJET COMMON /TOSSUNIT/ ITUNIT CHARACTER*8 SEG(2) 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) REAL*8 SID,STAID DATA SEG/'1000-100','099 - 00'/ LOGICAL LAST DATA PI180/.017453/, MAXNSEG/50/ C---------------------------------------------------------------------- C---------------------------------------------------------------------- LLEV = MCHK(1,1) LREP = FHD(LLEV) NSEG = 0 C LOOP OVER EACH SAT LEVEL IN THIS GROUP C -------------------------------------- DO 50 I=1,NXXYY LAST = I.EQ.NXXYY ILEV = MCHK(I,1) IREP = FHD(ILEV) C ITEMIZE COMPLETE PROFILES INTO SEGMENTS C --------------------------------------- 10 IF(IREP.EQ.LREP) THEN PRES = FPR(ILEV) OBER = FOE(ILEV) OMF = FOB(ILEV,2) IF(PRES.LE.2000. .AND. PRES.GE.100.) THEN K = 1 ELSE IF(PRES.LT.100. .AND. PRES.GT.0.) THEN K = 2 ENDIF NSEG(K) = NSEG(K) + 1 IF(NSEG(K).GT.MAXNSEG) THEN PRINT *, 'SATPUT - MAX NSEG EXCEEDED' CALL W3TAGE('PREPOBS_OIQCBUFR') CALL ERREXIT(99) ENDIF J = NSEG(K) C STORE THE RESULTS OF THE CHECKS C ------------------------------- ASD(J,K) = 0. OMA(J,K) = 0. SWT = 0. DO 20 L=1,3 AER = AE(I,L,1) FER = FE(I,1) AMF = AA(I,L,1) WT = MAX(.01,1.-(AER/FER)**2) ASD(J,K) = ASD(J,K) + WT*AER**2 OMA(J,K) = OMA(J,K) + (OMF-AMF)*WT SWT = SWT + WT 20 CONTINUE ASD(J,K) = SQRT(ASD(J,K)/SWT) OMA(J,K) = OMA(J,K)/SWT RAF(J,K) = ASD(J,K)/FER LEV(J,K) = ILEV FSD(J,K) = FER OSD(J,K) = OBER PRS(J,K) = PRES IF(.NOT.LAST) GOTO 50 ENDIF C A REPORT IS ITEMIZED AND READY TO JUDGE C --------------------------------------- STAID = SID(LREP) CLAT = FLT(LREP) CLON = FLN(LREP) KX = FKX(LREP) LREP = IREP C COMPUTE THE TOLERANCE VALUE C --------------------------- SINLAT = SIN(ABS(CLAT)*PI180) TT = TB*TSAT*SINLAT TMIN = 2.0 TT = MAX(TMIN,TT) C SUM THE RESULTS OVER THE PROFILE SEGMENTS - MAKE AND STORE DECISION C ------------------------------------------------------------------- DO 35 K=1,2 IF (NSEG(K).EQ.0) GOTO 35 QC1 = 0. QC2 = 0. QCT = 0. QCS = 0. OMA1 = 0. OMA2 = 0. EPV1 = 0. EPV2 = 0. 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,A10,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 FOI(LEV(J,K)) = KIQC 32 CONTINUE 35 CONTINUE C WRAP BACK TO START THE NEXT PROFILE C ----------------------------------- IF(.NOT.LAST) THEN NSEG = 0 GOTO 10 ENDIF C END OF LOOP OVER LEVELS C ----------------------- 50 CONTINUE RETURN END C----------------------------------------------------------------------- C----------------------------------------------------------------------- SUBROUTINE SEARCH PARAMETER (MAXREP=1000000) PARAMETER (MAXLEV=2000000) COMMON /OIQCDATA/ FLV (MAXREP ) , FLT (MAXREP ) , .FLN (MAXREP ) , FTM (MAXREP ) , FKX (MAXREP ) , .FNL (MAXREP ) , FMS (MAXREP ) , FMW (MAXREP,2 ) , .FQL (MAXREP ) , FHD (MAXLEV ) , FPRC(MAXLEV ) , .FPFD(MAXLEV ) , FPLV(MAXLEV ) , FPR (MAXLEV ) , .FQM (MAXLEV ) , FOI (MAXLEV ) , FOB (MAXLEV,2 ) , .FOE (MAXLEV ) , FFE (MAXLEV,2 ) , FTZ (MAXLEV ) , .FTW (MAXLEV,10,2) , FPF (MAXLEV,5,2 ) , FOC (MAXREP ) CHARACTER*3 FOC REAL XTW(MAXLEV,10,2),XPF(MAXLEV,5,2) PARAMETER (MAXLV=255) PARAMETER (MTOTH=10 ) PARAMETER (MTOTV=5 ) common /nddspt/ nd1(0:100),nd2(0:100) COMMON /CKLIST/ NDD,INDD(MAXLEV) COMMON /OBS/ NREP,NLEV,SID(MAXREP),INOB(MAXREP),IMAP(360,181) REAL*8 SID include "mpif.h" C----------------------------------------------------------------------- call mpi_comm_rank(MPI_COMM_WORLD,myid,ierr) call mpi_comm_size(MPI_COMM_WORLD,nprc,ierr) C----------------------------------------------------------------------- ! ------------------------------------------------------- ! divide the obs to check into equal parts for the solver ! ------------------------------------------------------- npart=ndd/nprc do i=0,nprc-1 n1=(i*npart)+1 n2=n1+npart-1 if(i==nprc-1) n2=ndd nd1(i)=n1;nd2(i)=n2 enddo C ------------------------------------------------------------ C farm out iterations of the search loop to processors in turn C ------------------------------------------------------------ ftw = 0 fpf = 0 np1=fhd(indd(1)) np2=fhd(indd(ndd)) DO N=np1,np2 IF(MOD(N,NPRC)/=MYID) CYCLE IF(FQL(N)>0) CALL SEARCHS(n) ENDDO c ------------------------------------------------------------- c each process puts its result in a section of the global array c ------------------------------------------------------------- c this code combines (reduces) all sections to one for all c -------------------------------------------------------- nftw=maxlev*10*2;nfpf=maxlev*5*2;mcomm=mpi_comm_world call mpi_allreduce(ftw,xtw,nftw,mpi_real,mpi_sum,mcomm,ier1) call mpi_allreduce(fpf,xpf,nfpf,mpi_real,mpi_sum,mcomm,ier2) ftw=xtw; fpf=xpf ! copy the reduction back where it belongs RETURN END C----------------------------------------------------------------------- C----------------------------------------------------------------------- SUBROUTINE SEARCHS(NPTR) PARAMETER (MAXREP=1000000) PARAMETER (MAXLEV=2000000) COMMON /OIQCDATA/ FLV (MAXREP ) , FLT (MAXREP ) , .FLN (MAXREP ) , FTM (MAXREP ) , FKX (MAXREP ) , .FNL (MAXREP ) , FMS (MAXREP ) , FMW (MAXREP,2 ) , .FQL (MAXREP ) , FHD (MAXLEV ) , FPRC(MAXLEV ) , .FPFD(MAXLEV ) , FPLV(MAXLEV ) , FPR (MAXLEV ) , .FQM (MAXLEV ) , FOI (MAXLEV ) , FOB (MAXLEV,2 ) , .FOE (MAXLEV ) , FFE (MAXLEV,2 ) , FTZ (MAXLEV ) , .FTW (MAXLEV,10,2) , FPF (MAXLEV,5,2 ) , FOC (MAXREP ) CHARACTER*3 FOC COMMON /OBS / NREP,NLEV,SID(MAXREP),INOB(MAXREP),IMAP(360,181) COMMON /OBHEIR/ QUALOB(300),binqadj common /trigs/ llnd(-10:10),rcos(181),prln(0:1100) PARAMETER(IDMX=8) PARAMETER(JDMX=5) PARAMETER(KDMX=2) PARAMETER(PMAX=.5) PARAMETER(MXLV=255) PARAMETER(MTOTH=10) PARAMETER(MTOTV=5) CHARACTER*8 SID CHARACTER*3 TAG DIMENSION NIND(MAXREP),RLAT(MAXREP),RLON(MAXREP), . RDIS(MAXREP),RDIR(MAXREP),INL(MXLV),INI(MXLV), . PATX(IDMX,JDMX,KDMX,MXLV), . MATX(IDMX,JDMX,KDMX,MXLV),NTOT(KDMX) LOGICAL SAME DATA PI180 /.0174532/ DATA RADE /6371./ DATA RSCAN /1000./ data mmrep/MAXREP/ SAVE PI180,RADE,MMREP,RSCAN,SAME C---------------------------------------------------------------------- C---------------------------------------------------------------------- C COPY SOME DETAILS ABOUT EACH OBSERVATION TO CHECK C ------------------------------------------------- PATX = 9999 MATX = 0 RCDR = 1./( 360./FLOAT(IDMX)) RCDS = 1./(RSCAN/FLOAT(JDMX)) IPTR = FLV(NPTR) JPTR = FLV(NPTR)+FNL(NPTR)-1 NLVL = FNL(NPTR) BLON = FLN(NPTR) BLAT = FLT(NPTR) BINQ = FQL(NPTR)+binqadj NBKX = FKX(NPTR) TAG = FOC(NPTR) C MAKE SPECIAL CASE ADJUSTMENTS (RECCOS,SATEMP) C --------------------------------------------- IF(MOD(NBKX,100).EQ.32) NC = 3 IF(MOD(NBKX,100).NE.32) NC = 8 IF(TAG.EQ.'SAT') BINQ = BINQ - 1. C ITEMIZE REPORTS IN THE SEACH CYLINDER C ------------------------------------- NPRF = 0 NOBS = 0 lon = ifix(blon)+1 lat = ifix(blat)+91. lscan = nint(rscan*.01) do ltd=-lscan,+lscan ltt = lat+ltd if(ltt.ge.1.and.ltt.le.181) then lnd = min(ceiling(llnd(ltd)*rcos(ltt)),180) ln1 = lon-lnd ln2 = lon+lnd lt1 = ltx(ln1,ltt) lt2 = ltx(ln2,ltt) DO N=IMAP(ln1,lt1),IMAP(ln2,lt2) NN = INOB(N) QLA = FQL(NN) IF(QLA.GE.0. .AND. QLA.LE.BINQ .and. nprf+1.le.mmrep) THEN NPRF = NPRF+1 NIND(NPRF) = NN RLON(NPRF) = FLN(NN) RLAT(NPRF) = FLT(NN) ENDIF ENDDO endif enddo C COMPUTE DISTANCE AND DIRECTION FROM THE CYLINDER CENTER C ------------------------------------------------------ CALL CHDIST(BLON,BLAT,RLON,RLAT,RDIS,RDIR,NIND,NPRF) C LOOP OVER OBSERVATION PROFILES TO BE USED FOR CHECKING C ------------------------------------------------------ DO N=1,NPRF SAME = SID(NPTR)(1:NC).EQ.SID(NIND(N))(1:NC) IF(RDIS(N).LT.RSCAN.AND..NOT.SAME) THEN IDIR = IFIX(RDIR(N)*RCDR)+1 JDIS = IFIX(RDIS(N)*RCDS)+1 LTYP = IFIX(FKX(NIND(N))*.01) IIND = FLV(NIND(N)) JIND = FLV(NIND(N))+FNL(NIND(N))-1 DO LIND=IIND,JIND IF(FOI(LIND).GT.0) THEN DO LPTR=IPTR,JPTR plog1 = prln(nint(fpr(lptr))) plog2 = prln(nint(fpr(lind))) PDIF = ABS(plog1-plog2) LEVL = LPTR-IPTR+1 J = JDIS DO WHILE(PDIF.GT.PATX(IDIR,J,LTYP,LEVL).AND.J.LE.JDMX) J = J+1 if(j.gt.jdmx) exit ! DAK: 3/27/13 - prevents overflow ! in subscript 2 (J) of array ! PATX ENDDO if(j.le.jdmx) then MATX(IDIR,J,LTYP,LEVL) = LIND PATX(IDIR,J,LTYP,LEVL) = PDIF endif ENDDO ENDIF ENDDO ENDIF ENDDO C STORE THE DATA INDEXES SELECTED FOR THE HORIZONTAL AND VERTICAL CHECKS C ---------------------------------------------------------------------- DO LPTR=IPTR,JPTR LEVL = LPTR-IPTR+1 NTOT = 0 DO K=1,KDMX IF(NTOT(K).LT.MTOTH) THEN DO J=1,JDMX DO I=1,IDMX LIND = MATX(I,J,K,LEVL) IF(LIND.GT.0 .AND. NTOT(K).LT.MTOTH) THEN NTOT(K) = NTOT(K)+1 FTW(LPTR,NTOT(K),K) = LIND ENDIF ENDDO ENDDO ENDIF ENDDO IF(TAG.EQ.'UPA') THEN IF(LEVL.EQ.1) CALL SRTPRS(IPTR,NLVL,INL,INI) DO J=1,MTOTV IF(INL(LEVL)-J.GT.0 ) FPF(LPTR,J,1) = INI(INL(LEVL)-J) IF(INL(LEVL)+J.LE.NLVL) FPF(LPTR,J,2) = INI(INL(LEVL)+J) ENDDO ENDIF ENDDO C END OF INDEPENDENT LOOP ITERATIONS C ---------------------------------- 100 CONTINUE RETURN END C----------------------------------------------------------------------- C----------------------------------------------------------------------- SUBROUTINE SRTPRS(ILEV,NLV,INL,INI) PARAMETER (MAXREP=1000000) PARAMETER (MAXLEV=2000000) COMMON /OIQCDATA/ FLV (MAXREP ) , FLT (MAXREP ) , .FLN (MAXREP ) , FTM (MAXREP ) , FKX (MAXREP ) , .FNL (MAXREP ) , FMS (MAXREP ) , FMW (MAXREP,2 ) , .FQL (MAXREP ) , FHD (MAXLEV ) , FPRC(MAXLEV ) , .FPFD(MAXLEV ) , FPLV(MAXLEV ) , FPR (MAXLEV ) , .FQM (MAXLEV ) , FOI (MAXLEV ) , FOB (MAXLEV,2 ) , .FOE (MAXLEV ) , FFE (MAXLEV,2 ) , FTZ (MAXLEV ) , .FTW (MAXLEV,10,2) , FPF (MAXLEV,5,2 ) , FOC (MAXREP ) CHARACTER*3 FOC DIMENSION INL(NLV),INI(NLV),PRS(NLV),SRT(NLV) C----------------------------------------------------------------------- C----------------------------------------------------------------------- DO L=1,NLV ILV = ILEV+L-1 INI(L) = ILV PRS(L) = FPR(ILV) SRT(L) = PRS(L) ENDDO C SORT INDEXES IN DESCENDING ORDER OF PRESSURE C -------------------------------------------- DO I=1 ,NLV DO J=I+1,NLV IF(SRT(I).LT.SRT(J)) THEN STMP = SRT(I) SRT(I) = SRT(J) SRT(J) = STMP ITMP = INI(I) INI(I) = INI(J) INI(J) = ITMP ENDIF ENDDO ENDDO C MAKE A SET OF LOOK BACK POINTERS C -------------------------------- DO 10 I=1,NLV INL(I) = I IF(PRS(I).NE.SRT(I)) THEN DO J=1,NLV IF(PRS(I).EQ.SRT(J)) THEN INL(I) = J GOTO 10 ENDIF ENDDO PRINT *, 'SRTPRS - SORTED PRESSURE NO FOUND' CALL W3TAGE('PREPOBS_OIQCBUFR') CALL ERREXIT(99) ENDIF 10 ENDDO RETURN END C---------------------------------------------------------------------- C---------------------------------------------------------------------- SUBROUTINE SELDAT(IND,NXY) PARAMETER (MAXREP=1000000) PARAMETER (MAXLEV=2000000) COMMON /OIQCDATA/ FLV (MAXREP ) , FLT (MAXREP ) , .FLN (MAXREP ) , FTM (MAXREP ) , FKX (MAXREP ) , .FNL (MAXREP ) , FMS (MAXREP ) , FMW (MAXREP,2 ) , .FQL (MAXREP ) , FHD (MAXLEV ) , FPRC(MAXLEV ) , .FPFD(MAXLEV ) , FPLV(MAXLEV ) , FPR (MAXLEV ) , .FQM (MAXLEV ) , FOI (MAXLEV ) , FOB (MAXLEV,2 ) , .FOE (MAXLEV ) , FFE (MAXLEV,2 ) , FTZ (MAXLEV ) , .FTW (MAXLEV,10,2) , FPF (MAXLEV,5,2 ) , FOC (MAXREP ) CHARACTER*3 FOC COMMON /HVECT / NXXYY,MCHK(64,2),MOBS(64,10,4,2), . 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/ C---------------------------------------------------------------------- C---------------------------------------------------------------------- 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) C PICK DATA FROM T AND W GROUPS C ----------------------------- DO 20 IG=1,2 NPIC = 0 DO 15 IP=1,MAXPIC(IG) IF(NPIC.EQ.MAXN(IG)) GOTO 20 IPIC = FTW(I,IP,IG) IF(IPIC.EQ.0) GOTO 20 IF(FOE(IPIC).GT.0) THEN NPIC = NPIC+1 DO 10 IQ=IG,IG*2-1 MDIM(M,IQ) = NPIC MOBS(M,NPIC,IQ,1) = IPIC MOBS(M,NPIC,IQ,2) = IQ 10 CONTINUE ENDIF 15 CONTINUE 20 CONTINUE C PICK DATA FROM SAME PROFILE C --------------------------- NPIC = 0 DO 30 IG=1,2 DO 25 IP=1,5 IPIC = FPF(I,IP,IG) IF(IPIC.EQ.0) GOTO 30 IF(FOE(IPIC).GT.0) THEN DO 24 KK=K,K*2-1 NPIC = NPIC+1 MDIM(M,4) = NPIC MOBS(M,NPIC,4,1) = IPIC MOBS(M,NPIC,4,2) = KK 24 CONTINUE GOTO 30 ENDIF 25 CONTINUE 30 CONTINUE 45 CONTINUE 50 CONTINUE C CHECK THE NUMBER OF OBS TO STORE C -------------------------------- IF(M.GT.NMX) THEN PRINT *, 'SELDAT - M GT NMX' CALL W3TAGE('PREPOBS_OIQCBUFR') CALL ERREXIT(99) ENDIF C STORE THE NUMBER OF OBS TO CHECK C -------------------------------- NXXYY = M RETURN END C$$$ SUBPROGRAM DOCUMENTATION BLOCK C C SUBPROGRAM: SETCON C PRGMMR: WOOLLEN ORG: NP20 DATE: 1992-07-29 C C ABSTRACT: INITIALIZES CONSTANTS AND PARAMETERS USING DATA STATEMENTS. C ALSO READS A FILE OF OBSERVATION ERRORS FROM UNIT 15. C C PROGRAM HISTORY LOG: C 1986-03-21 G. DIMEGO - ORIGINAL AUTHOR C 1988-11-24 D. DEAVEN - CODED FOR CYBER 205 C 1992-07-29 J. WOOLLEN - UPDATED VERSION FOR OIQC C C USAGE: CALL SETCON C INPUT FILES: C UNIT 18 - LIST OF IPOINTS FOR DETAILED PRINT IN UNITS 71-80 C C OUTPUT FILES: C UNIT 06 - STANDARD OUTPUT PRINT C C ATTRIBUTES: C LANGUAGE: FORTRAN 90 C MACHINE: NCEP WCOSS C C$$$ 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),binqadj COMMON /OBPRT / NPT,IDSTA(10),PSTA(10),KXSTA(10) common /trigs/ llnd(-10:10),rcos(181),prln(0:1100) CHARACTER*8 IDSTA DATA PI180 /.0174532 / C------------------------------------------------------------------- C------------------------------------------------------------------- c compute trig lookup tables for searching data c --------------------------------------------- do i=-10,10 llnd(i) = sqrt(float(100-i*i)) enddo do i=1,181 rcos(i) = abs(1./cos((float(i)-90.01)*pi180)) enddo do i=1,1100 prln(i) = log(float(i)) enddo prln(0) = prln(1) C STORE A SYMMETRIC MATRIX MAP C ---------------------------- LLL = 1 DO 10 I = 1,10 DO 10 J = 1,I IFA(I,J) = LLL IFA(J,I) = LLL LLL = LLL + 1 10 CONTINUE C MAKE A TABLE OF LENGTH SCALES BY MILIBARS C ----------------------------------------- DO 20 I=1,1500 CHLP(I) = PILNLNP(FLOAT(I),PMAND,CHL,21) 20 CONTINUE C READ IN SOME POINTS TO PRINT C ---------------------------- DO 30 I=1,10 READ(18,'(A6,I5,I4)',END=35,err=35) IDSTA(I),IP,KXSTA(I) PSTA(I) = IP PRINT'(A6,I5,I4)',IDSTA(I),IP,KXSTA(I) 30 CONTINUE 35 NPT = I-1 C SET UP SOME DATA PRIORITIES C --------------------------- binqadj = 0 DO 5 I=1,300 5 QUALOB(I) = 2. QUALOB(120) = 1. QUALOB(132) = 1. QUALOB(220) = 1. QUALOB(221) = 1. QUALOB(230) = 1. QUALOB(231) = 1. QUALOB(232) = 1. QUALOB(233) = 1. QUALOB(283) = 3. QUALOB(161) = 3. QUALOB(162) = 4. QUALOB(163) = 5. QUALOB(166) = 3. QUALOB(167) = 4. QUALOB(168) = 5. QUALOB(171) = 3. QUALOB(172) = 4. QUALOB(173) = 5. QUALOB(176) = 3. QUALOB(177) = 4. QUALOB(178) = 5. C SO LONG FOR NOW C --------------- RETURN END C----------------------------------------------------------------------- C----------------------------------------------------------------------- BLOCKDATA SETCONS 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) COMMON /OBPRT / NPT,IDSTA(10),PSTA(10),KXSTA(10) COMMON /OBHEIR/ QUALOB(300),binqadj COMMON /SATCOR/ EPSCOR(20) 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 EPSCOR/ .99, .88,. 61, .36, .18, .05,-.05,-.12,-.15,-.14, . -.11,-.08,-.05,-.03,-.02,-.02,-.01,-.01, 0., 0./ DATA CVC / 5.0 / END C$$$ SUBPROGRAM DOCUMENTATION BLOCK C C SUBPROGRAM: SSMIPTS C PRGMMR: WOOLLEN ORG: NP20 DATE: 1990-11-06 C C ABSTRACT: PUTS TOGETHER PACKAGES OF SSMI OBSERVATION LOCATIONS, C PERFORMS A MULTIVARIATE SURFACE WIND ANALYSIS, AND ASSIGNS ANALYSED C DIRECTION TO THE SSMI WIND SPEED OBSERVATION TO PRODUCE AN SSMI C WIND VECTOR TO BE CHECKED BY OIQC AND PASSED ON THE OBJECTIVE C ANALLYSIS PROGRAM. C C PROGRAM HISTORY LOG: C 1990-11-06 J. WOOLLEN - ORIGNIAL AUTHOR C C USAGE: CALL SSMIPTS C C REMARKS: NONE. C C ATTRIBUTES: C LANGUAGE: FORTRAN 90 C MACHINE: NCEP WCOSS C C$$$ SUBROUTINE SSMIPTS PARAMETER (MAXLEV=2000000) common /nddspt/ nd1(0:100),nd2(0:100) COMMON /CKLIST/ NDD,INDD(MAXLEV) COMMON /OBHEIR/ QUALOB(300),binqadj DATA NMX /64/ include "mpif.h" C----------------------------------------------------------------------- call mpi_comm_rank(MPI_COMM_WORLD,myid,ierr) c call nddsplit(myid,nd1,nd2) C----------------------------------------------------------------------- binqadj = -1 C MULTI-PROCESS QC OB LISTS C ------------------------- CALL SEARCH DO I=nd1(myid),nd2(myid),nmx CALL SELDAT(i,min(nmx,nd2(myid)-i+1)) CALL QCOI CALL SSMIPUT ENDDO binqadj = 0 RETURN END C$$$ SUBPROGRAM DOCUMENTATION BLOCK C C SUBPROGRAM: SSMIPUT C PRGMMR: WOOLLEN ORG: NP20 DATE: 1992-10-06 C C ABSTRACT: USES THE OUTPUT OF A MULTIVARIATE INTERPOLATION OF NEARBY C OBSERVATIONS TO ASSIGN DIRECTIONS TO SSMI WIND SPEED DATA. ALSO C PROVIDES FOR A DETAILED PRINT SHOWING OBSERVATIONS USED IN THE C INTERPOLATION ALONG WITH THE WEIGHTS AND INCREMENTS RESULTING. C C PROGRAM HISTORY LOG: C 1990-11-06 J. WOOLLEN - ORIGINAL AUTHOR C C USAGE: CALL SSMIPUT C C REMARKS: NONE. C C ATTRIBUTES: C LANGUAGE: FORTRAN 90 C MACHINE: NCEP WCOSS C C$$$ SUBROUTINE SSMIPUT PARAMETER (MAXREP=1000000) PARAMETER (MAXLEV=2000000) COMMON /OIQCDATA/ FLV (MAXREP ) , FLT (MAXREP ) , .FLN (MAXREP ) , FTM (MAXREP ) , FKX (MAXREP ) , .FNL (MAXREP ) , FMS (MAXREP ) , FMW (MAXREP,2 ) , .FQL (MAXREP ) , FHD (MAXLEV ) , FPRC(MAXLEV ) , .FPFD(MAXLEV ) , FPLV(MAXLEV ) , FPR (MAXLEV ) , .FQM (MAXLEV ) , FOI (MAXLEV ) , FOB (MAXLEV,2 ) , .FOE (MAXLEV ) , FFE (MAXLEV,2 ) , FTZ (MAXLEV ) , .FTW (MAXLEV,10,2) , FPF (MAXLEV,5,2 ) , FOC (MAXREP ) CHARACTER*3 FOC COMMON /HVECT / NXXYY,MCHK(64,2),MOBS(64,10,4,2), . 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) COMMON /OBPRT / NPT,SIDP(10),PRP(10),KXP(10) DIMENSION IOMF(3) CHARACTER*1 VAR2(2),VAR3(3) REAL*8 SID,CKSTA DATA VAR2 / 'Z' , 'W' / DATA VAR3 / 'Z' , 'U' , 'V' / DATA IOMF / 2 , 1 , 2 / CHARACTER*1 PLUMIN DATA PI180/.017453/ C---------------------------------------------------------------------- C---------------------------------------------------------------------- C LOOP OVER EACH OB CHECKED IN THIS GROUP C --------------------------------------- DO 50 N=1,NXXYY ILEV = MCHK(N,1) IREP = FHD(ILEV) PRES = FPR(ILEV) OBER = FOE(ILEV) CKSTA = SID(IREP) CLON = FLN(IREP) CLAT = FLT(IREP) KX = FKX(IREP) GESU = FMW(IREP,1) GESV = FMW(IREP,2) C COMPUTE THE ANALYSED WIND DIRECTION AT THIS OB C ---------------------------------------------- ANAU = GESU + AA(N,1,1) ANAV = GESV + AA(N,1,2) IF(ANAU.GT.0.) THEN ANADIR = ATAN(ANAV/ANAU) ELSE IF(ANAU.LT.0.) THEN ANADIR = ATAN(ANAV/ANAU) + 180.*PI180 ELSE IF(ANAU.EQ.0.) THEN ANADIR = SIGN(90.*PI180,ANAV) ENDIF C COMBINE THE ANA DIRECTION WITH THE SSMI SPEED C --------------------------------------------- SSMISPD = ABS(FOB(ILEV,2)) SSMIU = SSMISPD*COS(ANADIR) SSMIV = SSMISPD*SIN(ANADIR) C COMPUTE THE OBSERVED INCREMENT AND STORE IN THE PROPER PLACE C ------------------------------------------------------------ FOB(ILEV,1) = SSMIU-GESU FOB(ILEV,2) = SSMIV-GESV 50 CONTINUE 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 C----------------------------------------------------------------------- C C----------------------------------------------------------------------- SUBROUTINE SSMISRCH(IND,NXY) PARAMETER (MAXREP=1000000) PARAMETER (MAXLEV=2000000) COMMON /OIQCDATA/ FLV (MAXREP ) , FLT (MAXREP ) , .FLN (MAXREP ) , FTM (MAXREP ) , FKX (MAXREP ) , .FNL (MAXREP ) , FMS (MAXREP ) , FMW (MAXREP,2 ) , .FQL (MAXREP ) , FHD (MAXLEV ) , FPRC(MAXLEV ) , .FPFD(MAXLEV ) , FPLV(MAXLEV ) , FPR (MAXLEV ) , .FQM (MAXLEV ) , FOI (MAXLEV ) , FOB (MAXLEV,2 ) , .FOE (MAXLEV ) , FFE (MAXLEV,2 ) , FTZ (MAXLEV ) , .FTW (MAXLEV,10,2) , FPF (MAXLEV,5,2 ) , FOC (MAXREP ) CHARACTER*3 FOC 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) COMMON /HVECT / NXXYY,MCHK(64,2),MOBS(64,10,4,2), . MDIM(64,4),AA(64,4,2), . AE(64,4,2),WT(64,10,4,2), . FE(64,3) CHARACTER*3 TAG DIMENSION NIND(MAXREP),RLAT(MAXREP),RLON(MAXREP), . DIST(MAXREP),DIRN(MAXREP), . LIND(50000),LTYP(50000),PRES(50000), . ICOR(50000),IDIR(50000), . IPIC(10),CATCOR(101) LOGICAL CAT(50000,4,4),NONE(4,4),MORE,GO REAL*8 SID DATA CATCOR /21*5.,30*4.,20*3.,20*2.,10*1./ DATA PI180 /.0174532/, RADE /6371./, RSCAN /1000./ DATA MAXCYL /50000/ DATA MTOTH /10/ , MTOTV /5/ C---------------------------------------------------------------------- C---------------------------------------------------------------------- NXXYY = 0 MDIM = 0 C LOOP OVER PROFILES TO BE CHECKED C -------------------------------- DO 60 NPT=IND,IND+NXY-1 MPTR = INDD(NPT) NPTR = FHD(MPTR) BLON = FLN(NPTR) BLAT = FLT(NPTR) BINQ = FQL(NPTR)-1 TAG = FOC(NPTR) IF(TAG.NE.'SMI') THEN PRINT *, 'SSMISRCH - SEARCHING NON SSMI DATA' CALL W3TAGE('PREPOBS_OIQCBUFR') CALL ERREXIT(99) ENDIF C OUTLINE THE SEARCH AREA C ----------------------- DELLAT = RSCAN/(PI180*RADE) COSLAT = COS(BLAT*PI180) DELLON = MIN(DELLAT/COSLAT,180.0) X1 = BLON-DELLON X2 = BLON+DELLON Y1 = BLAT-DELLAT Y2 = BLAT+DELLAT IF(X1.LT. 0.) X1 = X1 + 360. IF(X2.GE.360.) X2 = X2 - 360. IF(Y1.LT.-90.) Y1 = -90. IF(Y2.GT. 90.) Y2 = 90. C ITEMIZE REPORTS IN THE SEACH CYLINDER C ------------------------------------- NPRF = 0 NOBS = 0 I1 = X1 + 1. I2 = X2 + 1. J1 = Y1 + 91. J2 = Y2 + 91. DO 10 K=1,2 IA = I1 IB = I2 IF(I1.GT.I2) THEN IF(K.EQ.1) IB = 360 IF(K.EQ.2) IA = 1 ELSE IF(K.EQ.2) GOTO 11 ENDIF DO 10 J=J1,J2 DO 10 N=IMAP(IA,J),IMAP(IB,J) NN = INOB(N) QLA = FQL(NN) IF(QLA.GE.0. .AND. QLA.LE.BINQ) THEN NPRF = NPRF+1 NIND(NPRF) = NN RLON(NPRF) = FLN(NN) RLAT(NPRF) = FLT(NN) ENDIF 10 CONTINUE 11 IF(NPRF.GT.MAXREP) THEN PRINT *, 'SSMISRCH - NPRF > MAXREP' CALL W3TAGE('PREPOBS_OIQCBUFR') CALL ERREXIT(99) ENDIF C COMPUTE DISTANCE AND DIRECTION FROM THE CYLINDER CENTER C ------------------------------------------------------- CALL CHDIST(BLON,BLAT,RLON,RLAT,DIST,DIRN,NIND,NPRF) C STORE ALL LEVELS OF EACH REPORT LOCATED IN THE CYLINDER C ------------------------------------------------------- DO 21 N=1,NPRF IF(DIST(N).LE.RSCAN .AND. SID(NPTR).NE.SID(NIND(N))) THEN ILEV = FLV(NIND(N)) KLEV = ILEV+FNL(NIND(N))-1 KX = FKX(NIND(N)) DO 20 JLEV=ILEV,KLEV IF(FOI(JLEV).GT.0) THEN NOBS = NOBS+1 LIND(NOBS) = JLEV PRES(NOBS) = LOG(FPR(JLEV)) DIST(NOBS) = (DIST(N)/RADE)**2 DIRN(NOBS) = AMOD(DIRN(N)/90.,4.) + 1.01 LTYP(NOBS) = KX/100 ENDIF 20 CONTINUE IF(NOBS.GT.MAXCYL) THEN PRINT *, 'SSMISRCH - NOBS>MAXCYL' CALL W3TAGE('PREPOBS_OIQCBUFR') CALL ERREXIT(99) ENDIF ENDIF 21 CONTINUE C LOOP OVER THE OBSERVATION LEVELS TO BE CHECKED C ---------------------------------------------- ILEV = FLV(NPTR) KLEV = ILEV+FNL(NPTR)-1 DO 50 IL=ILEV,KLEV BPRS = FPR(IL) IPRS = BPRS BPRS = LOG(BPRS) C INITIALIZE THE SLECTION ARRAYS C ------------------------------ NXXYY = NXXYY + 1 MCHK(NXXYY,1) = IL MCHK(NXXYY,2) = 2 FE (NXXYY,1) = FFE(IL,1) FE (NXXYY,2) = FFE(IL,2) FE (NXXYY,3) = FFE(IL,2) IF(NOBS.EQ.0) GOTO 50 C COMPUTE ZZ HRZ AND VRT CORRELATIONS FOR EACH OB AND CATEGORIZE C -------------------------------------------------------------- DO 24 N=1,NOBS CHRZ = EXP(-CHLP(IPRS)*(DIST(N))) CVRT = 1./(1.+ CVC*(BPRS-PRES(N))**2) IC = CHRZ*CVRT*100.+1. ICOR(N) = CATCOR(IC) IDIR(N) = DIRN(N) 24 CONTINUE DO 25 J=1,4 DO 25 I=1,4 DO 25 N=1,NOBS CAT(N,I,J) = ICOR(N).EQ.J .AND. IDIR(N).EQ.I 25 CONTINUE C PICK THE CLOSEST 10 (MULTIVARIATE) HORIZONTAL OBSERVATIONS C ---------------------------------------------------------- DO 28 J=1,4 DO 28 I=1,4 NONE(I,J) = .FALSE. 28 CONTINUE NTOT = 0 NPIC = 0 DO 35 J=1,4 30 MORE = .FALSE. DO 34 I=1,4 IF(NPIC.GE.MTOTH) GOTO 36 IF(NONE(I,J) ) GOTO 34 DO 32 N=1,NOBS IF(CAT(N,I,J)) THEN NTOT = NTOT+1 IPIC(NTOT) = N NPIC = NPIC+LTYP(N) CAT(N,I,J) = .FALSE. MORE = .TRUE. GOTO 34 ENDIF 32 CONTINUE NONE(I,J) = .TRUE. 34 CONTINUE IF(MORE) GOTO 30 35 CONTINUE 36 IF(NPIC.GT.MTOTH) NTOT = NTOT-1 C STORE THE DATA INDEXES SELECTED FOR HORIZONTAL CHECKS C ----------------------------------------------------- NPIC = 0 DO 40 N=1,NTOT NN = IPIC(N) KK = LTYP(NN) DO 40 K=KK,KK*2-1 NPIC = NPIC+1 MDIM(NXXYY,1) = MDIM(NXXYY,1) + 1 MOBS(NXXYY,NPIC,1,1) = LIND(NN) MOBS(NXXYY,NPIC,1,2) = K 40 CONTINUE 50 CONTINUE 60 CONTINUE RETURN END C$$$ SUBPROGRAM DOCUMENTATION BLOCK C C SUBPROGRAM: STOERR C PRGMMR: WOOLLEN ORG: NP20 DATE: 1990-11-06 C C ABSTRACT: INTERPOLATES TO AND STORES FORECAST AND OBSERVATION ERRORS C WITH EACH OBSERVATION IN THE DATASET IN COMMON /DUMMY/. ALSO MAKES C A LIST OF OBS BY 2D LOCATION IN COMMON /TGRID/ IN ORDER TO EXPEDITE C VARIOUS DATA SEARCH PROCEDURES. SOME CHECKING ON THE PLAUSIBILITY C OF OB COORDINATES IS PERFORMED AS WELL. C C PROGRAM HISTORY LOG: C 1990-11-06 J. WOOLLEN - ORIGINAL AUTHOR C C USAGE: CALL STOERR C OUTPUT FILES: C UNIT 06 - STANDARD OUTPUT PRINT C C REMARKS: C OUTPUT INFORMATION IN COMMON /TGRID/: C NREP - NUMBER OF REPORT LOCATIONS C XQC - ARRAY OF REPORT LONGITUDES C YQC - ARRAY OF REPORT LATITUDES C IQC - ARRAY OF REPORT INDEXES INTO THE BUFFS ARRAY C TQC - ARRAY OF REPORT TYPES (HEIGHT OR WIND) C FERR - ARRAY OF FORECAST ERRORS AT REPORT LOCATIONS C C C ATTRIBUTES: C LANGUAGE: FORTRAN 90 C MACHINE: NCEP WCOSS C C$$$ SUBROUTINE STOERR PARAMETER (MAXREP=1000000) PARAMETER (MAXLEV=2000000) COMMON /OIQCDATA/ FLV (MAXREP ) , FLT (MAXREP ) , .FLN (MAXREP ) , FTM (MAXREP ) , FKX (MAXREP ) , .FNL (MAXREP ) , FMS (MAXREP ) , FMW (MAXREP,2 ) , .FQL (MAXREP ) , FHD (MAXLEV ) , FPRC(MAXLEV ) , .FPFD(MAXLEV ) , FPLV(MAXLEV ) , FPR (MAXLEV ) , .FQM (MAXLEV ) , FOI (MAXLEV ) , FOB (MAXLEV,2 ) , .FOE (MAXLEV ) , FFE (MAXLEV,2 ) , FTZ (MAXLEV ) , .FTW (MAXLEV,10,2) , FPF (MAXLEV,5,2 ) , FOC (MAXREP ) CHARACTER*3 FOC 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 CLASS LOGICAL USEOE REAL*8 SID,STID DATA XMSG/32700./ include "mpif.h" C---------------------------------------------------------------------- call mpi_comm_rank(MPI_COMM_WORLD,myid,ierr) C----------------------------------------------------------------------- C GO THROUGH ALL THE STORED DATA C ------------------------------ DO 30 IREP=1,NREP ILEV = FLV(IREP) CLASS = FOC(IREP) STID = SID(IREP) DDLAT = FLT(IREP) DDLON = FLN(IREP) KLEV = FNL(IREP) KX = FKX(IREP) KKX = KX/100 IF(KLEV.EQ.0) GOTO 30 C CHECK FOR UNNATURAL INFORMATION IN THE HEADER C --------------------------------------------- IF( . ( KX.LT.120 .OR. KX.GT.299 ) .OR. . ( KLEV.LE.0 ) .OR. . ( DDLAT.LT.-90. .OR. DDLAT.GT.90. ) .OR. . ( DDLON.LT.-180. .OR. DDLON.GT.360.) . ) THEN PRINT*,'FROM STOERR---SKIP THIS MESS---FROM STOERR' PRINT'(''STATION= '',A8)',STID PRINT*,'LAT = ',DDLAT PRINT*,'LON = ',DDLON PRINT*,'KX = ',KX PRINT*,'NLEVS = ',KLEV PRINT*,'FROM STOERR---SKIP THIS MESS---FROM STOERR' GOTO 30 ENDIF C LOOP THRU THE REPORT LEVELS C --------------------------- DO 20 KL=1,KLEV L = ILEV+KL-1 C CHECK TO SEE IF A SURFACE OBSERVATION IS PRESENT C ------------------------------------------------ IF(FOE(L).EQ.0.) THEN GOTO 20 ENDIF C GET FORECAST ERRORS FOR THIS REPORT C ----------------------------------- DELP = FPR(L)*.053 DDPRS = FPR(L) DDPRA = DDPRS - DELP DDPRB = DDPRS + DELP ZERA = GETZER(DDLAT,DDLON,DDPRA) ZERB = GETZER(DDLAT,DDLON,DDPRB) WERR = GETWER(DDLAT,DDLON,DDPRS) C FIND THE OB ERROR FOR THIS REPORT C --------------------------------- IF(KKX.EQ.1) THEN OBERA = PILNLNP(DDPRA,PMAND,QIKE(1,KX),21) OBERB = PILNLNP(DDPRB,PMAND,QIKE(1,KX),21) FPE = 1./(1.+CVC*LOG(DDPRA/DDPRB)**2) OBER = SQRT(OBERA**2 + OBERB**2 - 2*OBERA*OBERB*FPE) ELSE IF(KKX.EQ.2) THEN OBER = PILNLNP(DDPRS,P50,QIKE(1,KX),21) ENDIF C STORE THE ERRORS IN THEIR RIGHTFUL PLACE C ---------------------------------------- IF(ZERA.EQ.0 .OR. ZERB.EQ.0. .OR. WERR.EQ.0.) THEN print*,'zera=',zera print*,'zerb=',zerb print*,'werr=',werr print*,'ddprs,ddpra,ddprb',ddprs,ddpra,ddprb PRINT *,'STOERR - FORECAST ERROR IS ZERO' CALL W3TAGE('PREPOBS_OIQCBUFR') CALL ERREXIT(99) 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 FFE(L,1) = TER FFE(L,2) = WER FTZ(L) = TZF C SPECIAL TREATMENT FOR SPECIAL CASES C ----------------------------------- IF(CLASS.EQ.'SFC' .AND. KL.EQ.2) THEN FFE(L,1) = (ZERA+ZERB)/2. OBER = (OBERA+OBERB)/2. FTZ(L) = 1. ENDIF IF(CLASS.NE.'SAT' .AND. DDPRS.GT.700.) THEN UPERR = 1.5*DDPRS/700. FFE(L,1) = FFE(L,1)*UPERR ENDIF C STORE AN OIQC OBSERVATION ERROR MAYBE C ------------------------------------- IF(USEOE) FOE(L) = OBER C END OF THE LOOPS C ---------------- 20 CONTINUE 30 CONTINUE IF(MYID.EQ.0) PRINT*,' IN STORERR - NREP=',NREP IF(MYID.EQ.0) PRINT* RETURN END C$$$ SUBPROGRAM DOCUMENTATION BLOCK C C SUBPROGRAM: STORE C PRGMMR: WOOLLEN ORG: NP20 DATE: 1992-07-28 C C ABSTRACT: STORES ALL DATA AND CORRESPONDING INDEX FOR LATER USE IN C DATA COLLECT AND SEARCH ROUTINES. C C PROGRAM HISTORY LOG: C 1986-03-21 G. DIMEGO - ORIGINAL AUTHOR C 1988-11-24 D. DEAVEN - CODED FOR CYBER 205 C 1990-11-06 J. WOOLLEN - MODIFIED COMMON /BUFFS/ FOR QCOI C 1992-07-28 J. WOOLLEN - MODIFIED FOR STORSTAR VERSION C 2019-06-20 S. MELCHIOR - Explicit declaration of bmiss C C USAGE: CALL STORE(LUDAT,LUBFI) C INPUT ARGUMENT LIST: C LUDAT - LOGICAL UNIT FOR NCEP PRODCUTION DATE FILE C LUBFI - LOGICAL UNIT FOR OBSERVATION (PREPBUFR) PRIOR TO C UPDATED QUALITY MARKS C C INPUT FILES: C UNIT LUDAT - NCEP PRODUTION DATE FILE C UNIT LUBFI - OBSERVATION (PREPBUFR) PRIOR TO UPDATED QUALITY C MARKS C C OUTPUT FILES: C UNIT 06 - STANDARD OUTPUT PRINT C UNIT LUDAT - NCEP PRODUTION DATE FILE C C ATTRIBUTES: C LANGUAGE: FORTRAN 90 C MACHINE: NCEP WCOSS C C$$$ SUBROUTINE STORE(LUDAT,LUBFI) PARAMETER (MAXREP=1000000) PARAMETER (MAXLEV=2000000) COMMON /OIQCDATA/ FLV (MAXREP ) , FLT (MAXREP ) , .FLN (MAXREP ) , FTM (MAXREP ) , FKX (MAXREP ) , .FNL (MAXREP ) , FMS (MAXREP ) , FMW (MAXREP,2 ) , .FQL (MAXREP ) , FHD (MAXLEV ) , FPRC(MAXLEV ) , .FPFD(MAXLEV ) , FPLV(MAXLEV ) , FPR (MAXLEV ) , .FQM (MAXLEV ) , FOI (MAXLEV ) , FOB (MAXLEV,2 ) , .FOE (MAXLEV ) , FFE (MAXLEV,2 ) , FTZ (MAXLEV ) , .FTW (MAXLEV,10,2) , FPF (MAXLEV,5,2 ) , FOC (MAXREP ) CHARACTER*3 FOC COMMON /OBS/ NREP,NLEV,SID(MAXREP),INOB(MAXREP),IMAP(360,181) COMMON /OBHEIR/ QUALOB(300),binqadj COMMON /RETCODE/ IRETCODE 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 REAL*8 HDR,OBS,FCS,QMS,OES,SID,BMISS 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 VMAX / 10E10 / DATA ROG / 29.261 / DATA QX / 300*.TRUE./ include "mpif.h" C---------------------------------------------------------------------- call mpi_comm_rank(MPI_COMM_WORLD,myid,ierr) C---------------------------------------------------------------------- C INITIALIZE C ---------- STUCK = .FALSE. NREP0 = 0 NLEV0 = 0 MAXNLEV = 0 C CHECK THE DATA DATE - PROCEED ANYWAY IF NO NMCDATE FILE TO BE FOUND C ------------------------------------------------------------------- CALL DATEBF(LUBFI,IY,IM,ID,IH,IDATE) READ(LUDAT,'(8X,I8)',END=1,ERR=900) JDATE IF(IDATE.NE.JDATE) GOTO 901 1 IF(IDATE.NE.JDATE) THEN REWIND LUDAT WRITE(LUDAT,'(8X,I8)') IDATE ENDIF C GET READY TO READ BUFR C ---------------------- bmiss=10e10; call setbmiss(bmiss) call ufbmem(lubfi,0,nmsg,iunit) call mpisplit(myid,nmsg,m1,m2) imsg = m1 if(myid.eq.0) then PRINT* PRINT*,'****** READING DATA IN STORE ******' PRINT* endif C LOOP THROUGH THE INPUT MESSAGES - READ THE NEXT SUBSET C ------------------------------------------------------ DO M=M1,M2 CALL READMM(IMSG,SUBSET,IDATE,IRET) 10 do while(IREADSB(LUBFI).EQ.0) CALL UFBINT(LUBFI,HDR,5, 1,IRET,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) CALL UFBINT(LUBFI,OES,5,255,KLEV,OESTR) CALL UFBCNT(LUBFI,IREC,IFLD) C DECIDE WHETHER OIQC SHOULD PROCESS - SET A TAG C ---------------------------------------------- KX = HDR(5) IF(QX(KX)) THEN IF(KX.GE.180 .AND. KX.LE.199) THEN TAG = 'SFC' ELSEIF(KX.GE.160 .AND. KX.LE.179) THEN TAG = 'SAT' ELSEIF(KX.EQ.283) THEN TAG = 'SMI' ELSE TAG = 'UPA' ENDIF ELSE GOTO 10 ENDIF C STORE THE HEADER C ---------------- NREP0 = NREP0+1 ILEV = 0 IF(NREP0.LE.MAXREP .AND. .NOT.STUCK) THEN NREP = NREP0 NLEV = NLEV0 SID(NREP) = HDR(1) FLN(NREP) = HDR(2) FLT(NREP) = HDR(3) FTM(NREP) = HDR(4) FKX(NREP) = KX FNL(NREP) = ILEV FLV(NREP) = NLEV+1 FOC(NREP) = TAG FQL(NREP) = QUALOB(KX) ELSE STUCK = .TRUE. ENDIF C STORE THE OB LEVELS C ------------------- DO 20 LEV=1,KLEV POB = OBS(1,LEV) TOB = OBS(3,LEV) ZOB = OBS(4,LEV) UOB = OBS(5,LEV) VOB = OBS(6,LEV) CAT = OBS(7,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) POE = OES(1,LEV) TOE = OES(3,LEV) ZOE = OES(4,LEV) WER = OES(5,LEV) IF(PQM.GE.8) GOTO 20 IF(PQM.EQ.8) TQM = MAX(PQM,TQM) IF(PQM.EQ.8) ZQM = MAX(PQM,ZQM) IF(PQM.EQ.8) WQQ = MAX(PQM,WQQ) PCHK = TAG.EQ.'SFC' IQM = 0 OB1 = 0 OB2 = 0 IOI = 0 OER = 0 IF(TQM.LT.8) THEN IQM = TQM PRG = .053*POB PA = POB - PRG PB = POB + PRG OB1 = TOB-TFC OB2 = ROG*OB1*LOG(PB/PA) 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 CONTINUE ILEV = ILEV+1 NLEV0 = NLEV0+1 IF(NLEV0.LE.MAXLEV .AND. .NOT.STUCK) THEN NLEV = NLEV0 IF(TAG.EQ.'SMI') THEN FMW(NREP,1) = UFC FMW(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. FNL (NREP) = ILEV FHD (NLEV) = NREP FPRC (NLEV) = IREC FPFD (NLEV) = IFLD FPLV (NLEV) = LEV FPR (NLEV) = POB FOB (NLEV,1) = OB1 FOB (NLEV,2) = OB2 FQM (NLEV) = IQM FOE (NLEV) = OER FOI (NLEV) = IOI ELSE STUCK = .TRUE. ENDIF IF(PCHK) THEN PCHK = .FALSE. IQM = MAX(PQM,ZQM) IF(IQM.LT.8) THEN IOI = 1 OB1 = 0 OB2 = ZOB-ZFC OER = 1 ELSEIF(IQM.EQ.8) THEN IOI = -2 OB1 = 0 OB2 = 0 OER = 0 ELSE IOI = 0 OB1 = 0 OB2 = 0 OER = 0 ENDIF GOTO 15 ENDIF 20 CONTINUE C KEEP TRACK OF MAX LEVELS - GO READ ANOTHER REPORT C ------------------------------------------------- MAXNLEV = MAX(ILEV,MAXNLEV) IF(ILEV.EQ.0) THEN NREP0 = NREP0-1 IF(.NOT.STUCK) NREP = NREP-1 ENDIF ENDDO ENDDO c syncronize the input buffers c ---------------------------- call syncstor(nrep0,nlev0,maxnlev) C all procs do the following - set tags C ------------------------------------- DO I=1,NREP KX = FKX(I) 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 FOC(I) = TAG ENDDO C ALL DONE - RECORD THE OUTCOME C ----------------------------- IF(MYID.EQ.0) THEN IF(NREP.EQ.NREP0 .AND. NLEV.EQ.NLEV0) IRETCODE = 0 IF(NREP.NE.NREP0 .OR . NLEV.NE.NLEV0) IRETCODE = 4 WRITE(6,'(A80)') IF(IRETCODE.NE.0) WRITE(6,'(A80)') WARNING IF(IRETCODE.EQ.0) WRITE(6,'(A80)') SUCCESS WRITE(6,101 ) NMSG WRITE(6,102 ) NREP0,NREP WRITE(6,103 ) NLEV0,NLEV WRITE(6,104 ) MAXNLEV IF(IRETCODE.NE.0) WRITE(6,'(A80)') WARNING IF(IRETCODE.EQ.0) WRITE(6,'(A80)') SUCCESS WRITE(6,'(A80)') ENDIF 101 FORMAT('TOTAL RECORDS: ',I10) 102 FORMAT('TOTAL REPORTS: ',I10,' --- STORED: ',I10) 103 FORMAT('TOTAL OBS : ',I10,' --- STORED: ',I10) 104 FORMAT('MAX OBS/PRFLE: ',I10) 105 CONTINUE RETURN 900 PRINT *, 'STORE - BAD OR MISSING NCEP DATE FILE' CALL W3TAGE('PREPOBS_OIQCBUFR') CALL ERREXIT(99) 901 PRINT *, 'STORE - BUFR AND NCEP DATES DONT MATCH' PRINT902,IDATE,JDATE CALL W3TAGE('PREPOBS_OIQCBUFR') CALL ERREXIT(99) 902 FORMAT(' BUFR DATE= ',I10.8,' NCEP DATE= ',I10.8) END C----------------------------------------------------------------------- C----------------------------------------------------------------------- SUBROUTINE MPISPLIT(MYID,NSPLIT,NS1,NS2) include "mpif.h" C----------------------------------------------------------------------- call mpi_comm_size(MPI_COMM_WORLD,nprc,ierr) C----------------------------------------------------------------------- npart = nsplit/nprc + min(mod(nsplit,nprc),1) ns1 = npart*myid+1 ns2 = min(npart*(myid+1),nsplit) return end C----------------------------------------------------------------------- C----------------------------------------------------------------------- SUBROUTINE syncstor(nrep0,nlev0,maxnlev) PARAMETER (MAXREP=1000000) PARAMETER (MAXLEV=2000000) COMMON /OIQCDATA/ FLV (MAXREP ) , FLT (MAXREP ) , .FLN (MAXREP ) , FTM (MAXREP ) , FKX (MAXREP ) , .FNL (MAXREP ) , FMS (MAXREP ) , FMW (MAXREP,2 ) , .FQL (MAXREP ) , FHD (MAXLEV ) , FPRC(MAXLEV ) , .FPFD(MAXLEV ) , FPLV(MAXLEV ) , FPR (MAXLEV ) , .FQM (MAXLEV ) , FOI (MAXLEV ) , FOB (MAXLEV,2 ) , .FOE (MAXLEV ) , FFE (MAXLEV,2 ) , FTZ (MAXLEV ) , .FTW (MAXLEV,10,2) , FPF (MAXLEV,5,2 ) , FOC (MAXREP ) CHARACTER*3 FOC COMMON /OBS/ NREP,NLEV,SID(MAXREP),INOB(MAXREP),IMAP(360,181) real*8 sid dimension nreps(5) data mrep/MAXREP/,mlev/MAXLEV/ include "mpif.h" dimension istat(mpi_status_size) dimension dummy(MAXLEV) real*8 dummy C----------------------------------------------------------------------- call mpi_comm_rank(MPI_COMM_WORLD,myid,ierr) call mpi_comm_size(MPI_COMM_WORLD,nprc,ierr) mcom = mpi_comm_world mrea = mpi_real mre8 = mpi_real8 mint = mpi_integer mchr = mpi_character mdes = 0 C----------------------------------------------------------------------- nreps(1) = nrep nreps(2) = nlev nreps(3) = nrep0 nreps(4) = nlev0 nreps(5) = maxnlev if(myid.ne.mdes) then call mpi_send(nreps,5,mint,mdes,myid,mcom,ierr) call mpi_send(sid ,nrep,mre8,mdes,myid,mcom,ierr) call mpi_send(flv ,nrep,mrea,mdes,myid,mcom,ierr) call mpi_send(flt ,nrep,mrea,mdes,myid,mcom,ierr) call mpi_send(fln ,nrep,mrea,mdes,myid,mcom,ierr) call mpi_send(ftm ,nrep,mrea,mdes,myid,mcom,ierr) call mpi_send(fkx ,nrep,mrea,mdes,myid,mcom,ierr) call mpi_send(fnl ,nrep,mrea,mdes,myid,mcom,ierr) call mpi_send(fms ,nrep,mrea,mdes,myid,mcom,ierr) call mpi_send(fmw(1,1),nrep,mrea,mdes,myid,mcom,ierr) call mpi_send(fmw(1,2),nrep,mrea,mdes,myid,mcom,ierr) call mpi_send(fql ,nrep,mrea,mdes,myid,mcom,ierr) call mpi_send(fhd ,nlev,mrea,mdes,myid,mcom,ierr) call mpi_send(fprc ,nlev,mrea,mdes,myid,mcom,ierr) call mpi_send(fpfd ,nlev,mrea,mdes,myid,mcom,ierr) call mpi_send(fplv ,nlev,mrea,mdes,myid,mcom,ierr) call mpi_send(fpr ,nlev,mrea,mdes,myid,mcom,ierr) call mpi_send(fqm ,nlev,mrea,mdes,myid,mcom,ierr) call mpi_send(foi ,nlev,mrea,mdes,myid,mcom,ierr) call mpi_send(foe ,nlev,mrea,mdes,myid,mcom,ierr) call mpi_send(fob(1,1),nlev,mrea,mdes,myid,mcom,ierr) call mpi_send(fob(1,2),nlev,mrea,mdes,myid,mcom,ierr) else do i=1,nprc-1 call mpi_recv(nreps,5,mint,i,i,mcom,istat,ierr) jrep = nreps(1) jlev = nreps(2) nrep0 = nrep0+nreps(3) nlev0 = nlev0+nreps(4) maxnlev = max(maxnlev,nreps(5)) if(nrep+jrep.le.mrep .and. nlev+jlev.le.mlev) then k = nrep+1 l = nlev+1 call mpi_recv(sid (k ) ,jrep,mre8,i,i,mcom,istat,ierr) call mpi_recv(flv (k ) ,jrep,mrea,i,i,mcom,istat,ierr) call mpi_recv(flt (k ) ,jrep,mrea,i,i,mcom,istat,ierr) call mpi_recv(fln (k ) ,jrep,mrea,i,i,mcom,istat,ierr) call mpi_recv(ftm (k ) ,jrep,mrea,i,i,mcom,istat,ierr) call mpi_recv(fkx (k ) ,jrep,mrea,i,i,mcom,istat,ierr) call mpi_recv(fnl (k ) ,jrep,mrea,i,i,mcom,istat,ierr) call mpi_recv(fms (k ) ,jrep,mrea,i,i,mcom,istat,ierr) call mpi_recv(fmw (k,1) ,jrep,mrea,i,i,mcom,istat,ierr) call mpi_recv(fmw (k,2) ,jrep,mrea,i,i,mcom,istat,ierr) call mpi_recv(fql (k ) ,jrep,mrea,i,i,mcom,istat,ierr) call mpi_recv(fhd (l ) ,jlev,mrea,i,i,mcom,istat,ierr) call mpi_recv(fprc(l ) ,jlev,mrea,i,i,mcom,istat,ierr) call mpi_recv(fpfd(l ) ,jlev,mrea,i,i,mcom,istat,ierr) call mpi_recv(fplv(l ) ,jlev,mrea,i,i,mcom,istat,ierr) call mpi_recv(fpr (l ) ,jlev,mrea,i,i,mcom,istat,ierr) call mpi_recv(fqm (l ) ,jlev,mrea,i,i,mcom,istat,ierr) call mpi_recv(foi (l ) ,jlev,mrea,i,i,mcom,istat,ierr) call mpi_recv(foe (l ) ,jlev,mrea,i,i,mcom,istat,ierr) call mpi_recv(fob (l,1) ,jlev,mrea,i,i,mcom,istat,ierr) call mpi_recv(fob (l,2) ,jlev,mrea,i,i,mcom,istat,ierr) flv(k:nrep+jrep) = flv(k:nrep+jrep)+l-1 fhd(l:nlev+jlev) = fhd(l:nlev+jlev)+k-1 nrep = nrep+jrep nlev = nlev+jlev else call mpi_recv(dummy(1),jrep,mre8,i,i,mcom,istat,ierr) do n=1,20 call mpi_recv(dummy(1),jlev,mrea,i,i,mcom,istat,ierr) enddo endif enddo endif call mpi_bcast(nrep,1,mint,mdes,mcom,ierr) call mpi_bcast(nlev,1,mint,mdes,mcom,ierr) call mpi_bcast(sid ,nrep,mre8,mdes,mcom,ierr) call mpi_bcast(flv ,nrep,mrea,mdes,mcom,ierr) call mpi_bcast(flt ,nrep,mrea,mdes,mcom,ierr) call mpi_bcast(fln ,nrep,mrea,mdes,mcom,ierr) call mpi_bcast(ftm ,nrep,mrea,mdes,mcom,ierr) call mpi_bcast(fkx ,nrep,mrea,mdes,mcom,ierr) call mpi_bcast(fnl ,nrep,mrea,mdes,mcom,ierr) call mpi_bcast(fms ,nrep,mrea,mdes,mcom,ierr) call mpi_bcast(fmw(1,1),nrep,mrea,mdes,mcom,ierr) call mpi_bcast(fmw(1,2),nrep,mrea,mdes,mcom,ierr) call mpi_bcast(fql ,nrep,mrea,mdes,mcom,ierr) call mpi_bcast(fhd ,nlev,mrea,mdes,mcom,ierr) call mpi_bcast(fprc ,nlev,mrea,mdes,mcom,ierr) call mpi_bcast(fpfd ,nlev,mrea,mdes,mcom,ierr) call mpi_bcast(fplv ,nlev,mrea,mdes,mcom,ierr) call mpi_bcast(fpr ,nlev,mrea,mdes,mcom,ierr) call mpi_bcast(fqm ,nlev,mrea,mdes,mcom,ierr) call mpi_bcast(foi ,nlev,mrea,mdes,mcom,ierr) call mpi_bcast(foe ,nlev,mrea,mdes,mcom,ierr) call mpi_bcast(fob(1,1),nlev,mrea,mdes,mcom,ierr) call mpi_bcast(fob(1,2),nlev,mrea,mdes,mcom,ierr) RETURN END C$$$ SUBPROGRAM DOCUMENTATION BLOCK C C SUBPROGRAM: TOSSLIST C PRGMMR: WOOLLEN ORG: NP20 DATE: 1990-11-06 C C ABSTRACT: WRITES A LIST OF REJECTED OBSERVATIONS PLUS A SUMMARY. C C PROGRAM HISTORY LOG: C 1990-11-06 J. WOOLLEN - ORIGINAL AUTHOR C C USAGE: CALL TOSSLIST(LO) C INPUT ARGUMENT LIST: C LO - LOGICAL UNIT FOR TOSSLIST OUTPUT C C OUTPUT FILES: C UNIT 06 - STANDARD OUTPUT PRINT C UNIT LO - OFFICIAL OPERATIONAL TOSSLIST C C REMARKS: NONE. C C ATTRIBUTES: C LANGUAGE: FORTRAN 90 C MACHINE: NCEP WCOSS C C$$$ SUBROUTINE TOSSLIST(LO) PARAMETER (MAXREP=1000000) PARAMETER (MAXLEV=2000000) COMMON /OIQCDATA/ FLV (MAXREP ) , FLT (MAXREP ) , .FLN (MAXREP ) , FTM (MAXREP ) , FKX (MAXREP ) , .FNL (MAXREP ) , FMS (MAXREP ) , FMW (MAXREP,2 ) , .FQL (MAXREP ) , FHD (MAXLEV ) , FPRC(MAXLEV ) , .FPFD(MAXLEV ) , FPLV(MAXLEV ) , FPR (MAXLEV ) , .FQM (MAXLEV ) , FOI (MAXLEV ) , FOB (MAXLEV,2 ) , .FOE (MAXLEV ) , FFE (MAXLEV,2 ) , FTZ (MAXLEV ) , .FTW (MAXLEV,10,2) , FPF (MAXLEV,5,2 ) , FOC (MAXREP ) CHARACTER*3 FOC COMMON /OBS/ NREP,NLEV,SID(MAXREP),INOB(MAXREP),IMAP(360,181) CHARACTER*3 TAG DIMENSION IBAD(30,3),ITOT(30),ISFZ(30),PZ(30),PU(30),PT(30) LOGICAL SFZ,SAT REAL*8 SID,STAID DATA XMSG/32767./ C---------------------------------------------------------------------- C----------------------------------------------------------------------- C START UP C -------- IBAD = 0 ITOT = 0 ISFZ = 0 PZ = 0 PU = 0 PT = 0 C COUNT BAD DATA IN THE BIG BUFFER - WRITE THE TOSSLIST C ----------------------------------------------------- N = 0 DO 10 IREP=1,NREP X = FLN(IREP) Y = FLT(IREP) STAID = SID(IREP) KX = FKX(IREP) TAG = FOC(IREP) ILEV = FLV(IREP) KLEV = ILEV+FNL(IREP)-1 DO 10 JLEV=ILEV,KLEV SFZ = TAG.EQ.'SFC' .AND. JLEV.EQ.KLEV SAT = TAG.EQ.'SAT' IF(FOI(JLEV).NE.0) THEN IF(SFZ) THEN ISFZ(KX/10) = ISFZ(KX/10) + 1 ELSE ITOT(KX/10) = ITOT(KX/10) + 1 ENDIF ENDIF IF(FOI(JLEV).LT.0) THEN IP = FPR(JLEV) IQ = FQM(JLEV) OBER = FOE(JLEV) BUFF11 = FOB(JLEV,1) BUFF12 = FOB(JLEV,2) IF(FOI(JLEV).EQ.-1) THEN IF(SFZ) THEN N = N+1 WRITE(LO,2002) N,STAID,X,Y,IP,KX,IQ,OBER,BUFF11,BUFF12 IBAD(KX/10,3) = IBAD(KX/10,3) + 1 ELSE IF(.NOT.SAT) THEN N = N+1 WRITE(LO,2000) N,STAID,X,Y,IP,KX,IQ,OBER,BUFF11,BUFF12 ENDIF IBAD(KX/10,1) = IBAD(KX/10,1) + 1 ENDIF ELSE IF(FOI(JLEV).EQ.-2) THEN N = N+1 WRITE(LO,2001) N,STAID,X,Y,IP,KX,IQ,OBER,BUFF11,BUFF12 IBAD(KX/10,2) = IBAD(KX/10,2) + 1 ENDIF ENDIF 10 CONTINUE C COMPUTE GROSS PERCENTAGE OF TOSSES C ---------------------------------- DO 20 I=10,29 IF(ISFZ(I).GT.0) THEN PZ(I) = FLOAT(IBAD(I,3))/FLOAT(ISFZ(I)) ENDIF IF(ITOT(I).GT.0) THEN PU(I) = FLOAT(IBAD(I,2))/FLOAT(ITOT(I)) PT(I) = FLOAT(IBAD(I,1))/FLOAT(ITOT(I)) ENDIF 20 CONTINUE C WRITE A SHORT SUMMARY OF TOSSES C ------------------------------- WRITE( 6,2003)(I,I=100,290,10),(IBAD(I,3),I=10,29),(PZ(I),I=10,29) WRITE( 6,2004)(I,I=100,290,10),(IBAD(I,2),I=10,29),(PU(I),I=10,29) WRITE( 6,2005)(I,I=100,290,10),(IBAD(I,1),I=10,29),(PT(I),I=10,29) C FORMAT STATEMENTS FOR THE TOSS LIST C ----------------------------------- 2000 FORMAT(' ',I5,1X,'J',3X,A8,2X,4X,2F10.2,I9,2X,2I5,F8.1,F13.1, . 5X,F13.1,10X,F6.1) 2001 FORMAT(' ',I5,1X,'P',3X,A8,2X,4X,2F10.2,I9,2X,2I5,F8.1,F13.1, . 5X,F13.1,10X,F6.1) 2002 FORMAT(' ',I5,1X,'Z',3X,A8,2X,4X,2F10.2,I9,2X,2I5,F8.1,F13.1, . 5X,F13.1,10X,F6.1) 2003 FORMAT(/' SURFACE Z TOSSES '/20I5/20I5/20F5.2) 2004 FORMAT(/' BELOW GROUND TOSSES'/20I5/20I5/20F5.2) 2005 FORMAT(/' QC OI TOSSES '/20I5/20I5/20F5.2) RETURN END CCBREAKCC C$$$ SUBPROGRAM DOCUMENTATION BLOCK C C SUBPROGRAM: VSOLVE C PRGMMR: WOOLLEN ORG: NP20 DATE: 1990-11-06 C C ABSTRACT: CHOLESKY TYPE SOLUTION FOR ARRAYS OF POSITIVE DEFINITE C SYMMETRIC MATRIXES. C C PROGRAM HISTORY LOG: C 1990-11-06 J. WOOLLEN - ORIGINAL AUTHOR C C USAGE: CALL VSOLVE(A,B,NDIM,BAD,NS,MAXDIM) C INPUT ARGUMENT LIST: C A - ARRAY OF SYMMETRIC MATRIXES C B - ARRAY OF RIGHT HAND SIDE VECTORS C NDIM - ARRAY OF MATRIX RANKS C BAD - BAD MATRIX INDICATOR C NFT - NUMBER OF RIGHT HAND SIDES PER MATRIX C NS - NUMBER OF MATRIXES TO SOLVE C MAXDIM - LARGEST RANK MATRIX IN STACK C C OUTPUT ARGUMENT LIST: C B - CONTAINS THE SOLUTIONS UPON RETURN C C REMARKS: TRIANGULAR REPRESENTATION LISTS ELEMENTS TOWARDS THE C DIAGONAL. ALGORITHM FROM H.CARUS. C C ATTRIBUTES: C LANGUAGE: FORTRAN 90 C MACHINE: NCEP WCOSS C C$$$ SUBROUTINE VSOLVE (A,B,NDIM,BAD,NS,MAXDIM) PARAMETER( NMX = 64 , . NN = 10*(10+1)/2 , . N = 10 ) DIMENSION A(NMX,NN),B(NMX,N,2),NDIM(NMX),BAD(NMX) LOGICAL BAD DIMENSION T(NMX) DATA CNUM/1.E-15/, NFT/2/ C---------------------------------------------------------------------- IX (I,J) = I*(I-1)/2 + J C---------------------------------------------------------------------- DO 1 M=1,NS 1 BAD(M) = .FALSE. C DECOMPOSE THE MATRIXES C ---------------------- DO 10 I=1,MAXDIM DO 10 J=1,I DO 2 M=1,NS 2 T(M) = A(M,IX(I,J)) DO 3 K=1,J-1 DO 3 M=1,NS 3 T(M) = T(M) - A(M,IX(I,K)) * A(M,IX(J,K)) IF(I.GT.J) THEN DO 4 M=1,NS 4 A(M,IX(I,J)) = T(M) * A(M,IX(J,J)) ELSE DO 5 M=1,NS IF(T(M).LT.CNUM .AND. NDIM(M).GE.I) BAD(M) = .TRUE. IF(T(M).LE.0) T(M) = 1. 5 CONTINUE DO 6 M=1,NS 6 A(M,IX(I,I)) = 1./SQRT(T(M)) ENDIF 10 CONTINUE C SOLVE FOR ALL RIGHT HAND SIDES C ------------------------------ DO 40 NF=1,NFT C FORWARD SUBSTITUTION C -------------------- DO 20 I=1,MAXDIM DO 11 M=1,NS 11 T(M) = B(M,I,NF) DO 12 J=1,I-1 DO 12 M=1,NS 12 T(M) = T(M) - A(M,IX(I,J)) * B(M,J,NF) DO 20 M=1,NS 20 B(M,I,NF) = T(M) * A(M,IX(I,I)) C BACKWARD SUBSTITUTION C --------------------- DO 30 I=MAXDIM,1,-1 DO 21 M=1,NS 21 T(M) = B(M,I,NF) IF(I.NE.MAXDIM) THEN DO 22 J=I+1,MAXDIM DO 22 M=1,NS 22 T(M) = T(M) - A(M,IX(J,I)) * B(M,J,NF) ENDIF DO 30 M=1,NS 30 B(M,I,NF) = T(M) * A(M,IX(I,I)) 40 CONTINUE RETURN END C----------------------------------------------------------------------- C----------------------------------------------------------------------- SUBROUTINE syncndd(arr,n,mdtp) PARAMETER (MAXLEV=2000000) common /nddspt/ nd1(0:100),nd2(0:100) COMMON /CKLIST/ NDD,INDD(MAXLEV) dimension arr(n) include "mpif.h" dimension istat(mpi_status_size) dimension ireq(100) C----------------------------------------------------------------------- call mpi_comm_rank(MPI_COMM_WORLD,myid,ierr) call mpi_comm_size(MPI_COMM_WORLD,nprc,ierr) mcom = mpi_comm_world mdes = 0 C----------------------------------------------------------------------- if(myid.ne.mdes) then ib = indd(nd1(myid)) ie = indd(nd2(myid)) mlen = ie-ib+1 if(nd2(myid).gt.0) . call mpi_send(arr(ib),mlen,mdtp,mdes,myid,mcom,ierr) else do i=1,nprc-1 ib = indd(nd1(i)) ie = indd(nd2(i)) mlen = ie-ib+1 if(nd2(i).gt.0) . call mpi_recv(arr(ib),mlen,mdtp,i,i,mcom,istat,ierr) enddo endif call mpi_bcast (arr,n,mdtp,mdes,mcom,ierr) RETURN END C----------------------------------------------------------------------- C----------------------------------------------------------------------- SUBROUTINE NDDSPLIT(MYID,nd1,nd2) PARAMETER (MAXREP=1000000) PARAMETER (MAXLEV=2000000) COMMON /OIQCDATA/ FLV (MAXREP ) , FLT (MAXREP ) , .FLN (MAXREP ) , FTM (MAXREP ) , FKX (MAXREP ) , .FNL (MAXREP ) , FMS (MAXREP ) , FMW (MAXREP,2 ) , .FQL (MAXREP ) , FHD (MAXLEV ) , FPRC(MAXLEV ) , .FPFD(MAXLEV ) , FPLV(MAXLEV ) , FPR (MAXLEV ) , .FQM (MAXLEV ) , FOI (MAXLEV ) , FOB (MAXLEV,2 ) , .FOE (MAXLEV ) , FFE (MAXLEV,2 ) , FTZ (MAXLEV ) , .FTW (MAXLEV,10,2) , FPF (MAXLEV,5,2 ) , FOC (MAXREP ) CHARACTER*3 FOC COMMON /CKLIST/ NDD,INDD(MAXLEV) include "mpif.h" C----------------------------------------------------------------------- call mpi_comm_size(MPI_COMM_WORLD,nprc,ierr) C----------------------------------------------------------------------- C SUM THE REPORTS PER AMAP BOX C ---------------------------- JTOT = 0 DO N=1,NDD JTOT = JTOT+JNDD(N) ENDDO C COMPUTE THE range within the FREQUENCY DISTRIBUTION C --------------------------------------------------- npart = jtot/nprc + min(mod(jtot,nprc),1) ns1 = min(npart*myid,jtot) ns2 = min(npart*(myid+1),jtot) c print*,'ns',ns1,ns2 c pick the start and end points for this processor c ------------------------------------------------ jntt = 0 do n=1,ndd nd1 = n jntt = jntt+jndd(n) c print*,'n',n,jntt,ndd if(jntt.ge.ns1) then do m=min(n+1,ndd),ndd nd2 = m jntt = jntt+jndd(m) c print*,'m',m,jntt,ns2 if(jntt.ge.ns2) return enddo if(ndd.eq.1) return PRINT *, 'nddsplit inside alignment error' CALL W3TAGE('PREPOBS_OIQCBUFR') CALL ERREXIT(99) endif enddo c bad news if it gets this far! c ----------------------------- PRINT *, 'nddsplit outside alignment error' CALL W3TAGE('PREPOBS_OIQCBUFR') CALL ERREXIT(99) end C----------------------------------------------------------------------- C----------------------------------------------------------------------- function jndd(n) PARAMETER (MAXREP=1000000) PARAMETER (MAXLEV=2000000) COMMON /OIQCDATA/ FLV (MAXREP ) , FLT (MAXREP ) , .FLN (MAXREP ) , FTM (MAXREP ) , FKX (MAXREP ) , .FNL (MAXREP ) , FMS (MAXREP ) , FMW (MAXREP,2 ) , .FQL (MAXREP ) , FHD (MAXLEV ) , FPRC(MAXLEV ) , .FPFD(MAXLEV ) , FPLV(MAXLEV ) , FPR (MAXLEV ) , .FQM (MAXLEV ) , FOI (MAXLEV ) , FOB (MAXLEV,2 ) , .FOE (MAXLEV ) , FFE (MAXLEV,2 ) , FTZ (MAXLEV ) , .FTW (MAXLEV,10,2) , FPF (MAXLEV,5,2 ) , FOC (MAXREP ) CHARACTER*3 FOC COMMON /OBS/ NREP,NLEV,SID(MAXREP),INOB(MAXREP),IMAP(360,181) COMMON /CKLIST/ NDD,INDD(MAXLEV) REAL*8 SID DATA PI180 /.0174532/ save iold,jold C----------------------------------------------------------------------- !ltb(l) = min(max(1,l),181) !lnb(l) = min(max(1,l),360) C----------------------------------------------------------------------- if(n.eq.1) iold = 0 irep = fhd(indd(n)) if(iold.eq.irep) then jndd = jold return else iold = irep jndd = 1 endif lon = ifix(fln(irep))+1 lat = ifix(flt(irep))+91. ldn = 10 do ltt=lat-5,lat+5 if(ltt.ge.2.and.ltt.le.180) then ln1 = lon-ldn ln2 = lon+ldn lt1 = ltx(ln1,ltt) lt2 = ltx(ln2,ltt) jndd = jndd+min(abs(imap(ln1,lt1)-imap(ln2,lt2)),ldn*100) endif enddo jold = jndd return end C----------------------------------------------------------------------- C----------------------------------------------------------------------- function ltx(lon,lat) C----------------------------------------------------------------------- ltb(l) = min(max(1,l),181) C----------------------------------------------------------------------- if(lon.lt.1) then do while(lon<1);lon=360+lon;enddo ltx = ltb(lat-1) elseif(lon.gt.360) then do while(lon>360);lon=lon-360;enddo ltx = ltb(lat+1) else ltx = lat endif return end c----------------------------------------------------------------------------------- c----------------------------------------------------------------------------------- subroutine prttime(title) character(*) title data t0/0.00/ data t1/0.00/ save t0,t1 c----------------------------------------------------------------------------------- include "mpif.h" c----------------------------------------------------------------------------------- call mpi_comm_rank(MPI_COMM_WORLD,myid,ierr);if(myid/=0)return c----------------------------------------------------------------------------------- if(t0==0.00) t0=secnds(t0) print1,title,secnds(t0)-t1,secnds(t0) t1=secnds(t0) 1 format('time>>>',a20,3(2x,f8.2)) return end