C=============================================================================C C International Comprehensive Ocean-Atmosphere Data Set (ICOADS) 14 Nov 2012 C C Filename:level: qctrf1.f:01A Fortran 77 program C C Function: LMR/IMMA quality control and trimming flags Author: S.Lubker C C=============================================================================C C Software Revision Information (previous version: 22 Aug 2012, qctrf.f:01J): C Relative humidity passed to (instead of computed in) {flgrpt}. C-----------------------------------------------------------------------3456789 C Software documentation for the (invariant) user-interface routines {flgrpt, C qctrf}. C C Functionality: Calling {qctrf} will compute flags for the quality control C (Attm1) and trimming (Attm2) attachments of LMR. The output flags are C stored as coded values in integer arrays QCFLG and TRFLG, respectively; C input data are passed via the remaining arguments (please refer to {rdlmr6} C for the definition of all data structures). C C Similarly, data in IMMA format can be flagged by calling {flgrpt} from C within {rdimma}. In this case, the output flags are stored as floating- C point values in FTRUE locations ND through QCZ (omitting 2-degree box C number, and adding 10- and 1-degree box numbers; with input data passed C via the remaining arguments (refer to {rdimma}). C C External library required: {gsbytes.f}. External data files required C (available from: http://icoads.noaa.gov/software/other/other/ ): C ZKT2: Decadal Summary Untrimmed Limits (DSUL; ref. COADS Release 1 C supp. C. These are the original trimming limits for periods: C 1854-1909, 1910-49, 1950-79. Data after 1979 (before 1854) C are trimmed using the 1950-79 (1854-1909) limits. Beginning C with Release 2.5 relative humidity before 1910 is trimmed C using the 1910-49 limits. C LNDLM2_32: 32-bit version of 1-deg box landmask and 5-deg box QC limits C used for NCDC-QC (ref. COADS Release 1, supp. J). This is C a reformatting of data received 18 November 1981 from NCDC. C Machine dependencies: Assumes 32-bit word. C=============================================================================C C WARNING: Code beyond this point should not require any modification. C C=============================================================================C C-----------------------------------------------------------------------3456789 SUBROUTINE FLGRPT(FTRUE,FMISS,FERR,FUNITS +,B10,YR,MO,DY,HR,LON,LAT,DCK,SID,PT,D,W,VV,WW,W1,SLP +,AT,WBT,DPT,SST,N,NH,CL,H,CM,CH,WD,WP,WH,SD,SP,SH,A,PPP +,B1,ND,ZNC,LZ,QCZ,RH) C RDIMMA INTERFACE TO QCTRF IMPLICIT INTEGER(A-E,G-Z) DIMENSION FTRUE(*),FUNITS(*) +,QCFLG(15),TRFLG(23) LOGICAL B10XY ITRUE(I)=NINT(FTRUE(I)/FUNITS(I)) C IF (FTRUE(YR).EQ.FMISS .OR. FTRUE(YR).EQ.FERR +.OR. FTRUE(MO).EQ.FMISS .OR. FTRUE(MO).EQ.FERR +.OR. FTRUE(LON).EQ.FMISS .OR. FTRUE(LON).EQ.FERR +.OR. FTRUE(LAT).EQ.FMISS .OR. FTRUE(LAT).EQ.FERR) RETURN IF (FTRUE(DCK).EQ.FMISS .OR. FTRUE(DCK).EQ.FERR +.OR. FTRUE(SID).EQ.FMISS .OR. FTRUE(SID).EQ.FERR) RETURN IF (FTRUE(B10).EQ.FMISS) THEN IF (.NOT.B10XY(ITRUE(LON),ITRUE(LAT),BOX10)) STOP 'B10XY' FTRUE(B10)=BOX10 ENDIF FTRUE(B1)=MSQ(ITRUE(LON),ITRUE(LAT),MSQ1,MSQ5) CALL QCTRF(FTRUE,FMISS,FERR,FUNITS,QCFLG,TRFLG +,B10,YR,MO,DY,HR,LON,LAT,DCK,SID,PT,D,W,VV,WW,W1,SLP +,AT,WBT,DPT,SST,N,NH,CL,H,CM,CH,WD,WP,WH,SD,SP,SH,A,PPP,RH) DO 180 I=1,14 FTRUE(ZNC-1+I)=QCFLG(I) 180 CONTINUE DO 190 I=2,8 FTRUE(ND-2+I)=TRFLG(I) 190 CONTINUE IF (NINT(FTRUE(ND)).EQ.0)FTRUE(ND)=FMISS QCE=LZ-1 FTRUE(QCE)=TRFLG( 9)/2*2*2*2*2*2 + +TRFLG(10)/2*2*2*2*2 + +TRFLG(11)/2*2*2*2 + +TRFLG(12)/2*2*2 + +TRFLG(13)/2*2 + +TRFLG(14)/2 IF (NINT(FTRUE(QCE)).EQ.0)FTRUE(QCE)=FMISS FTRUE(LZ)=TRFLG(18) IF (NINT(FTRUE(LZ)).EQ.0)FTRUE(LZ)=FMISS FTRUE(QCZ)=TRFLG(19)*2*2*2*2 + +TRFLG(20)*2*2*2 + +TRFLG(21)*2*2 + +TRFLG(22)*2 + +TRFLG(23) IF (NINT(FTRUE(QCZ)).EQ.0)FTRUE(QCZ)=FMISS END C-----------------------------------------------------------------------3456789 SUBROUTINE QCTRF(FTRUE,FMISS,FERR,FUNITS,QCFLG,TRFLG +,B10,YR,MO,DY,HR,LON,LAT,DCK,SID,PT,D,W,VV,WW,W1,SLP +,AT,WBT,DPT,SST,N,NH,CL,H,CM,CH,WD,WP,WH,SD,SP,SH,A,PPP,RH) C I-COADS QUALITY CONTROL AND TRIMMING FLAGS IMPLICIT INTEGER(A-E,G-Z) INTEGER*8 QLTY,ATT1,ATT2 DIMENSION FTRUE(*),FUNITS(*),QCFLG(*),TRFLG(*) COMMON /QC/MSQ1,MSQ5,INVNQC(43),INVNF(14,11),MISS,ILL,LNDLCK(100) +,M4PT8(4,12,5),P4PT8(4,12,5),M5PT8(4,12,5),P5PT8(4,12,5) +/QC0/SHIPF,WINDF,VISF,PRSWXF,PSTWXF,PRESSF,DRYF,WETF,DEWF,SEAF +,CLOUDF,SEAWVF,SWLWVF,PTENDF,VALU(10) COMMON /ATT/AL(15),AD(255,15) +,SUP(255),SUPLEN,ERRNUM(51),ERRLEN(51),ERR(15,51),NERR DATA BOX/0/ SAVE ITRUE(A)=IFUNC(FTRUE(A),FUNITS(A),FMISS,FERR) C MISS=NINT(FMISS) ILL=NINT(FERR) C BOX10=ITRUE(B10) IF (BOX10.NE.BOX) THEN CALL LNDLMT(BOX10) BOX=BOX10 ENDIF B1=MSQ(ITRUE(LON),ITRUE(LAT),MSQ1,MSQ5) CALL DUPQC(QLTY,ATT1 +,ITRUE(SID),ITRUE(YR),ITRUE(MO),ITRUE(DY),ITRUE(HR),ITRUE(LAT) +,ITRUE(D),ITRUE(W),ITRUE(VV),ITRUE(WW),ITRUE(W1),ITRUE(SLP) +,ITRUE(AT),ITRUE(WBT),ITRUE(DPT),ITRUE(SST),ITRUE(N),ITRUE(NH) +,ITRUE(CL),ITRUE(H),ITRUE(CM),ITRUE(CH),ITRUE(WD),ITRUE(WP) +,ITRUE(WH),ITRUE(SD),ITRUE(SP),ITRUE(SH),ITRUE(A),ITRUE(PPP) +,RPT) AD(1,1)=SHIPF AD(2,1)=WINDF AD(3,1)=VISF AD(4,1)=PRSWXF AD(5,1)=PSTWXF AD(6,1)=PRESSF AD(7,1)=DRYF AD(8,1)=WETF AD(9,1)=DEWF AD(10,1)=SEAF AD(11,1)=CLOUDF AD(12,1)=SEAWVF AD(13,1)=SWLWVF AD(14,1)=PTENDF AD(15,1)=QLTY DO 1 I=1,15 QCFLG(I)=AD(I,1) 1 CONTINUE C CALL TRIM(ATT2 +,ITRUE(B10),ITRUE(YR),ITRUE(MO),ITRUE(HR),ITRUE(LON),ITRUE(LAT) +,ITRUE(DCK),ITRUE(SID),ITRUE(PT) +,FTRUE(D),FTRUE(W),FTRUE(SLP),FTRUE(AT),FTRUE(DPT),FTRUE(SST) +,FTRUE(RH) +,MISS,ILL) DO 2 I=1,23 TRFLG(I)=AD(I,2) 2 CONTINUE END C-----------------------------------------------------------------------3456789 FUNCTION IFUNC(FTRUE,FUNITS,FMISS,FERR) C CONVERT FTRUE TO ITRUE IMPLICIT INTEGER(A-E,G-Z) C IF (FTRUE.EQ.FMISS) THEN IFUNC=NINT(FMISS) ELSE IF (FTRUE.EQ.FERR) THEN IFUNC=NINT(FERR) ELSE IFUNC=NINT(FTRUE/FUNITS) ENDIF END C-----------------------------------------------------------------------3456789 LOGICAL FUNCTION B10XY(X,Y,B10) C CONVERT LON AND LAT IN HUNDREDTHS TO B10 C IMPLICIT INTEGER(A-E,G-Z) C B10XY=.FALSE. IF(X.LT.0.OR.X.GT.35999.OR.ABS(Y).GT.9000)RETURN C IF(X.EQ.0.OR.X.GT.18000)THEN C=35-MOD(36000-X,36000)/1000 ELSE C=X/1000 ENDIF C R=8-SIGN(MIN(ABS(Y),8999)/1000,Y) IF(Y.LT.0)R=R+1 C B10=R*36+MOD(C+36-3,36)+1 B10XY=.TRUE. END C-----------------------------------------------------------------------3456789 BLOCK DATA BDATT C COMMON BLOCK DATA STATEMENTS IMPLICIT INTEGER(A-E,G-Z) COMMON /ATT/AL(15),AD(255,15) +,SUP(255),SUPLEN,ERRNUM(51),ERRLEN(51),ERR(15,51),NERR C DATA AL/15*0/,SUPLEN/0/,NERR/0/ END C$B25-------------------------------------------------------------------3456789 FUNCTION B25(B2,B0,B26) C BOX2 WITHIN BOX10 1-25 IMPLICIT INTEGER(A-Z) IF(B2.GT.1.AND.B2.LT.16202)THEN B25=MOD((B2-2)/180,5)*5+MOD(B2-2,5)+1 ELSE IF(B2.EQ.1)THEN B25=B0 ELSE B25=B26 ENDIF END C$B2QXY-----------------------------------------------------------------3456789 LOGICAL FUNCTION B2QXY(Q,X,Y,B2) IMPLICIT INTEGER(A-E,G-Z) C B2QXY=.FALSE. IF(Q.LT.1.OR.Q.GT.4.OR.X.LT.0.OR.X.GT.35999.OR.ABS(Y).GT.9000) +RETURN C IF(ABS(Y).LT.9000)THEN IF(MOD(Q,2).EQ.0)THEN C=MIN(X/200,89) ELSE C=179-MIN(MOD(36000-X,36000)/200,89) ENDIF IF(Q/3.EQ.0)THEN R=89-(9000+Y)/200 ELSE R=(9000-Y)/200 ENDIF B2=2+R*180+C C ELSE IF(Y.EQ.9000)THEN B2=1 C ELSE B2=16202 ENDIF C B2QXY=.TRUE. END C$GETDSUL---------------------------------------------------------------3456789 SUBROUTINE GETDSUL(B10,MO,B2,YR,CODED,FTRUE,FMISS) C READ DECADAL SUMMARY UNTRIMMED LIMITS IMPLICIT INTEGER(A-E,G-Z) DIMENSION BUF(12,3,25) COMMON /DSUL1/FTRUEL(23),FTRUEU(23) +,FUNITS(23),FBASE(23),BITS(23),OFFSET(23) DIMENSION CODED(5+3*6),FTRUE(5+3*6) DATA UNIT/7/ DATA OLD/0/ SAVE LOGICAL OPEN C INQUIRE(UNIT,OPENED=OPEN) IF(.NOT.OPEN) +OPEN(UNIT,ACCESS='DIRECT',FORM='UNFORMATTED',RECL=4*12*3*25 +,FILE='ZKT2') C NEW=REC(B10,MO,B2) IF (NEW.NE.OLD) THEN READ(UNIT,REC=NEW)BUF OLD=NEW ENDIF CALL RPTGET(BUF(1,PERIOD(YR),B25(B2,1,1)) +,CODED,FTRUE,FMISS,FUNITS,FBASE,BITS,OFFSET,23,0,5) C C PRINT 200,(FTRUE(I),I=1,5),((FTRUE(5+(J-1)*3+I),J=1,6),I=1,3) 200 FORMAT(/' BOX10 ',F4.0,' MONTH ',F3.0,' BOX2 ',F6.0 +,' PERIOD ',F5.0,' CHECKSUM ',F6.0/ +2X,7X,'S',7X,'A',7X,'U',7X,'V',7X,'P',7X,'R'/ +1X,'L',5F8.2,F8.1/ +1X,'G',5F8.2,F8.1/ +1X,'U',5F8.2,F8.1) END C$PERIOD----------------------------------------------------------------3456789 FUNCTION PERIOD(YR) C DSUL PERIOD 1-3 IMPLICIT INTEGER(A-Z) PERIOD=1+YR/1910+YR/1950 END C$REC-------------------------------------------------------------------3456789 FUNCTION REC(BOX,MO,B2) C DSUL DIRECT ACCESS RECORD NUMBER IMPLICIT INTEGER(A-Z) C IF(B2.EQ.1)THEN B10=1 ELSE IF(B2.EQ.16202)THEN B10=648 ELSE B10=BOX ENDIF C REC=(B10-1)*12+MO IF(B10.EQ.1)THEN REC=REC+MO IF(B2.EQ.1)REC=REC-1 ELSE REC=REC+12 IF(B10.EQ.648)THEN REC=REC+MO IF(B2.NE.16202)REC=REC-1 ENDIF ENDIF END C$BDDSUL1---------------------------------------------------------------3456789 BLOCK DATA BDDSUL1 IMPLICIT INTEGER(A-E,G-Z) C COMMON /DSUL1/FTRUEL(23),FTRUEU(23) +,FUNITS(23),FBASE(23),BITS(23),OFFSET(23) C DATA FTRUEL/3*1.,1800.,0. +,3*-5.,3*-88.,6*-102.2,3*870.,3*0./ C DATA FTRUEU/648.,12.,16202.,2054.,4094. +,3*40.,3*58.,6*102.2,3*1074.6,3*100./ C DATA FUNITS/5*1. +,15*.01,3*.1/ C DATA FBASE/3*0.,1799.,0. +,3*-501.,3*-8801.,6*-10221.,3*86999.,3*-1./ C DATA BITS/10,4,14,8,12 +,18*16/ C DATA OFFSET/ 16, 26, 30, 44, 52 +, 64, 80, 96,112,128,144,160,176,192,208 +,224,240,256,272,288,304,320,336/ END C$GETRPT----------------------------------------------------------------3456789 SUBROUTINE RPTGET(RPT,CODED,FTRUE,FMISS +,FUNITS,FBASE,BITS,OFFSET,NUMBER,RPTID,INDXCK) C UNPACK REPORT AND CONVERT CODED TO TRUE VALUES C IMPLICIT INTEGER(A-E,G-Z) C CHARACTER*(*) RPT DIMENSION CODED(*),FTRUE(*),FUNITS(*),FBASE(*),BITS(*),OFFSET(*) C c IF(MOD(ICHAR(RPT(2:2)),16).NE.RPTID)STOP 'RPTID' C CODED(INDXCK)=0 DO 190 I=1,NUMBER IF(I.EQ.INDXCK)GOTO 190 CALL GBYTE(RPT,CODED(I),OFFSET(I),BITS(I)) IF(CODED(I).EQ.0)THEN FTRUE(I)=FMISS ELSE FTRUE(I)=(CODED(I)+FBASE(I))*FUNITS(I) CODED(INDXCK)=CODED(INDXCK)+CODED(I) ENDIF 190 CONTINUE I=INDXCK FTRUE(I)=MOD(CODED(I),2**BITS(I)-1) CALL GBYTE(RPT,CODED(I),OFFSET(I),BITS(I)) IF(FTRUE(I).NE.CODED(I))STOP 'CHECKSUM' END C$MSQ-------------------------------------------------------------------3456789 INTEGER FUNCTION MSQ(LON,LAT,MSQ1,MSQ5) IMPLICIT INTEGER(A-Z) UNITS(B)=MOD(B/100,10) IF (ABS(LAT).NE.9000) THEN MSQ1=UNITS(ABS(LAT))*10 + +UNITS(MIN(ABS(LON-LON/18000*36000),17999)) ELSE MSQ1=99 ENDIF MSQ5=MSQ1/50*2+MOD(MSQ1,10)/5+1 MSQ=MSQ1 END C$NIGHT-----------------------------------------------------------------3456789 INTEGER FUNCTION NIGHT(MO,HR,LON,B2,MISS,ILL) C 0=MISSING, 1=NIGHT, 2=DAY IMPLICIT INTEGER(A-Z) COMMON /DT/DT(90,12) IF (HR.EQ.MISS .OR. HR.EQ.ILL) THEN NIGHT=0 ELSE TABLE=DT((MIN(MAX(B2,2),16201)-2)/180+1,MO) IF (ABS(MOD(HR*.01+LON*.01/15.,24.)-12.)*10..GT.TABLE + .OR. TABLE.EQ.0) THEN NIGHT=1 ELSE NIGHT=2 ENDIF ENDIF END C$BDDT------------------------------------------------------------------3456789 BLOCK DATA BDDT C -----HALF LENGTH OF THE DURATION OF DAYLIGHT IN TENTHS C (ODD LATITUDES FROM 89 TO -89 BY MONTH) IMPLICIT INTEGER(A-Z) COMMON /DT/DT(90,12) DATA ((DT(I,J),J=1,12),I= 1,10) +/ 0, 0, 0,120,120,120,120,120,120, 0, 0, 0 +, 0, 0, 28,120,120,120,120,120,120, 0, 0, 0 +, 0, 0, 42,120,120,120,120,120, 88, 0, 0, 0 +, 0, 0, 48,120,120,120,120,120, 79, 0, 0, 0 +, 0, 0, 51,120,120,120,120,120, 74, 14, 0, 0 +, 0, 0, 52,100,120,120,120,120, 72, 27, 0, 0 +, 0, 0, 54, 91,120,120,120,120, 70, 33, 0, 0 +, 0, 20, 54, 86,120,120,120,107, 68, 38, 0, 0 +, 0, 27, 55, 82,120,120,120, 97, 67, 41, 0, 0 +, 0, 32, 56, 79,114,120,120, 91, 66, 43, 11, 0/ DATA ((DT(I,J),J=1,12),I=11,20) +/ 0, 35, 56, 77,102,120,120, 87, 66, 45, 20, 0 +, 16, 38, 57, 75, 96,120,106, 84, 65, 46, 26, 0 +, 23, 40, 57, 74, 91,105, 99, 82, 65, 48, 30, 15 +, 27, 42, 57, 73, 88, 98, 94, 80, 64, 49, 33, 22 +, 30, 43, 57, 72, 85, 94, 90, 78, 64, 50, 36, 26 +, 33, 45, 58, 71, 83, 90, 87, 77, 64, 50, 38, 30 +, 36, 46, 58, 70, 81, 88, 85, 75, 63, 51, 40, 32 +, 38, 47, 58, 69, 79, 85, 83, 74, 63, 52, 41, 35 +, 39, 48, 58, 69, 78, 83, 81, 73, 63, 52, 43, 37 +, 41, 49, 58, 68, 77, 81, 79, 72, 63, 53, 44, 39/ DATA ((DT(I,J),J=1,12),I=21,30) +/ 42, 50, 58, 67, 75, 80, 78, 71, 63, 53, 45, 40 +, 44, 50, 58, 67, 74, 78, 77, 70, 62, 54, 46, 42 +, 45, 51, 59, 66, 73, 77, 76, 70, 62, 54, 47, 43 +, 46, 52, 59, 66, 72, 76, 74, 69, 62, 55, 48, 44 +, 47, 52, 59, 66, 71, 75, 73, 68, 62, 55, 49, 45 +, 48, 53, 59, 65, 71, 74, 72, 68, 62, 55, 50, 46 +, 49, 53, 59, 65, 70, 73, 72, 67, 62, 56, 50, 47 +, 50, 54, 59, 64, 69, 72, 71, 67, 62, 56, 51, 48 +, 50, 54, 59, 64, 69, 71, 70, 66, 61, 56, 52, 49 +, 51, 55, 59, 64, 68, 70, 69, 66, 61, 57, 52, 50/ DATA ((DT(I,J),J=1,12),I=31,40) +/ 52, 55, 59, 64, 67, 69, 68, 65, 61, 57, 53, 51 +, 52, 55, 59, 63, 67, 68, 68, 65, 61, 57, 54, 52 +, 53, 56, 59, 63, 66, 68, 67, 64, 61, 57, 54, 52 +, 54, 56, 59, 63, 66, 67, 66, 64, 61, 58, 55, 53 +, 54, 57, 59, 62, 65, 66, 66, 64, 61, 58, 55, 54 +, 55, 57, 59, 62, 64, 66, 65, 63, 61, 58, 56, 54 +, 55, 57, 60, 62, 64, 65, 65, 63, 61, 58, 56, 55 +, 56, 58, 60, 62, 63, 64, 64, 63, 61, 58, 57, 56 +, 57, 58, 60, 61, 63, 64, 63, 62, 61, 59, 57, 56 +, 57, 58, 60, 61, 63, 63, 63, 62, 60, 59, 58, 57/ DATA ((DT(I,J),J=1,12),I=41,50) +/ 58, 59, 60, 61, 62, 63, 62, 62, 60, 59, 58, 57 +, 58, 59, 60, 61, 62, 62, 62, 61, 60, 59, 58, 58 +, 59, 59, 60, 61, 61, 61, 61, 61, 60, 60, 59, 59 +, 59, 60, 60, 60, 61, 61, 61, 61, 60, 60, 59, 59 +, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60 +, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60, 60 +, 61, 60, 60, 60, 59, 59, 59, 59, 60, 60, 61, 61 +, 61, 61, 60, 59, 59, 59, 59, 59, 60, 60, 61, 61 +, 62, 61, 60, 59, 58, 58, 58, 59, 60, 61, 62, 62 +, 62, 61, 60, 59, 58, 57, 58, 58, 60, 61, 62, 63/ DATA ((DT(I,J),J=1,12),I=51,60) +/ 63, 62, 60, 59, 57, 57, 57, 58, 60, 61, 62, 63 +, 63, 62, 60, 59, 57, 56, 57, 58, 59, 61, 63, 64 +, 64, 62, 60, 58, 57, 56, 56, 57, 59, 62, 63, 64 +, 65, 63, 60, 58, 56, 55, 55, 57, 59, 62, 64, 65 +, 65, 63, 61, 58, 56, 54, 55, 57, 59, 62, 64, 66 +, 66, 63, 61, 58, 55, 54, 54, 56, 59, 62, 65, 66 +, 66, 64, 61, 57, 54, 53, 54, 56, 59, 62, 65, 67 +, 67, 64, 61, 57, 54, 52, 53, 56, 59, 63, 66, 68 +, 68, 65, 61, 57, 53, 52, 52, 55, 59, 63, 66, 68 +, 68, 65, 61, 56, 53, 51, 52, 55, 59, 63, 67, 69/ DATA ((DT(I,J),J=1,12),I=61,70) +/ 69, 65, 61, 56, 52, 50, 51, 54, 59, 63, 68, 70 +, 70, 66, 61, 56, 51, 49, 50, 54, 59, 64, 68, 71 +, 70, 66, 61, 56, 51, 48, 49, 53, 58, 64, 69, 72 +, 71, 67, 61, 55, 50, 47, 48, 53, 58, 64, 70, 73 +, 72, 67, 61, 55, 49, 46, 48, 52, 58, 65, 70, 74 +, 73, 68, 61, 54, 49, 45, 47, 52, 58, 65, 71, 75 +, 74, 68, 61, 54, 48, 44, 46, 51, 58, 65, 72, 76 +, 75, 69, 61, 54, 47, 43, 44, 50, 58, 66, 73, 77 +, 76, 70, 62, 53, 46, 42, 43, 50, 58, 66, 74, 78 +, 78, 70, 62, 53, 45, 40, 42, 49, 57, 67, 75, 80/ DATA ((DT(I,J),J=1,12),I=71,80) +/ 79, 71, 62, 52, 43, 39, 41, 48, 57, 67, 76, 81 +, 81, 72, 62, 51, 42, 37, 39, 47, 57, 68, 77, 83 +, 82, 73, 62, 51, 41, 35, 37, 46, 57, 68, 79, 85 +, 84, 74, 62, 50, 39, 32, 35, 45, 57, 69, 80, 88 +, 87, 75, 62, 49, 37, 30, 33, 43, 56, 70, 82, 90 +, 90, 77, 63, 48, 35, 26, 30, 42, 56, 70, 84, 94 +, 93, 78, 63, 47, 32, 22, 26, 40, 56, 71, 87, 98 +, 97, 80, 63, 46, 29, 15, 21, 38, 55, 72, 90,105 +,104, 82, 63, 45, 24, 0, 14, 36, 55, 74, 94,120 +,120, 85, 64, 43, 18, 0, 0, 33, 54, 75,100,120/ DATA ((DT(I,J),J=1,12),I=81,90) +/120, 88, 64, 41, 6, 0, 0, 29, 54, 77,109,120 +,120, 93, 65, 38, 0, 0, 0, 23, 53, 79,120,120 +,120,100, 66, 34, 0, 0, 0, 13, 52, 82,120,120 +,120,120, 66, 29, 0, 0, 0, 0, 50, 87,120,120 +,120,120, 68, 20, 0, 0, 0, 0, 48, 93,120,120 +,120,120, 69, 0, 0, 0, 0, 0, 46,106,120,120 +,120,120, 72, 0, 0, 0, 0, 0, 41,120,120,120 +,120,120, 78, 0, 0, 0, 0, 0, 32,120,120,120 +,120,120, 92, 0, 0, 0, 0, 0, 0,120,120,120 +,120,120,120, 0, 0, 0, 0, 0, 0,120,120,120/ END C$QB10------------------------------------------------------------------3456789 FUNCTION QB10(B10) IMPLICIT INTEGER(A-E,G-Z) QB10=-1 IF(B10.LT.1.OR.B10.GT.648)RETURN QB10=2+(B10-1)/324*2-MOD(B10-1+3,36)/18 END C$RDPT------------------------------------------------------------------3456789 FUNCTION RDPT(SLP,AT,WBT,FMISS,FERR) C COMPUTE DEW POINT TEMPERATURE RDPT=FMISS IF(SLP.EQ.FMISS .OR. SLP.EQ.FERR)THEN P=1015. ELSE P=SLP ENDIF IF(AT.EQ.FMISS .OR. AT.EQ.FERR)RETURN IF(WBT.EQ.FMISS .OR. WBT.EQ.FERR)RETURN IF(AT.LT.WBT)RETURN ACON=7.5 BCON=237.3 ESW=6.1078*10.**(WBT*ACON/(WBT+BCON)) E=ESW-(.00066*P)*(((.00115*WBT)+1)*(AT-WBT)) IF(E.LT.0.)RETURN CCON=ALOG10(E/6.1078) RDPT=BCON*CCON/(ACON-CCON) END C$RELHUM----------------------------------------------------------------3456789 FUNCTION RELHUM(AT,DPT,MISS,ILL) IMPLICIT REAL(A-Z) IF ( AT.EQ.MISS .OR. AT.EQ.ILL +.OR. NINT(AT*10.).LT.-880 +.OR. NINT(AT*10.).GT. 580 +.OR. DPT.EQ.MISS .OR. DPT.EQ.ILL +.OR. NINT((AT-DPT)*10.).LT. -5 +.OR. NINT((AT-DPT)*10.).GT.700)THEN RELHUM=MISS ELSE RELHUM=MIN(HUM(AT,MIN(AT,DPT)),100.) ENDIF END C$SCALAR----------------------------------------------------------------3456789 FUNCTION SCALAR(D,W,FMISS,FLAG) C SCALAR WIND INTEGER FLAG,R,A,B,J,K,L,M,N,Q,S PARAMETER (R=1 ,A=2 ,B=3 ,J=4 ,K=5 ,L=6 ,M=7 ,N=8 ,Q=9 ,S=10) PARAMETER (C=3.14159265359/180.) SAVE C IF(NINT(D).GE.1 .AND. NINT(D).LE.360)THEN I=0 ELSE IF(NINT(D).EQ.361)THEN I=4 ELSE IF(NINT(D).EQ.362)THEN I=8 ELSE I=12 ENDIF C IF(NINT(W*10.).EQ.0)THEN I=I+1 ELSE IF(NINT(W*10.).GE.1 .AND. NINT(W*10.).LE.31)THEN I=I+2 ELSE IF(NINT(W*10.).GE.32 .AND. NINT(W*10.).LE.1022)THEN I=I+3 ELSE I=I+4 ENDIF C IF(I.EQ.1 .AND. FLAG.EQ.M)THEN GOTO 801 ELSE IF(I.EQ.8 .AND. (FLAG.EQ.A .OR. FLAG.EQ.J))THEN SCALAR=0. ELSE SCALAR=W ENDIF C GOTO(1,2,2,3 ,2,4,4,1 ,4,2,5,3 ,1,1,3,6)I 1 IF(FLAG.EQ.A .OR. FLAG.EQ.J .OR. FLAG.EQ.M)GOTO 150 GOTO 800 2 IF(FLAG.EQ.R .OR. FLAG.EQ.J)GOTO 150 GOTO 800 3 IF(FLAG.EQ.M)GOTO 150 GOTO 800 4 IF(FLAG.EQ.A .OR. FLAG.EQ.J)GOTO 150 GOTO 800 5 IF(FLAG.EQ.J)GOTO 150 GOTO 800 6 IF(FLAG.EQ.M .OR. FLAG.EQ.S)GOTO 150 GOTO 800 C 150 IF(FLAG.LE.J) +GOTO(8,9,9,7 ,8,7,9,8 ,8,7,7,7 ,8,7,7,7)I C 7 U=FMISS V=FMISS RETURN C 8 U=0. V=0. RETURN C 9 ANG=MIN(D,360.)*C U=-W*SIN(ANG) V=-W*COS(ANG) RETURN C ENTRY UW() C VECTOR WIND EASTWARD COMPONENT UW=U RETURN C ENTRY VW() C VECTOR WIND NORTHWARD COMPONENT VW=V RETURN C 800 PRINT *,' ERROR IN ROUTINE SCALAR D=',D,', W=',W,', FLAG=',FLAG 801 SCALAR=FMISS U=FMISS V=FMISS END C$TRIM------------------------------------------------------------------3456789 SUBROUTINE TRIM(ATT2,B10,YR,MO,HR,LON,LAT,DCK,SID,PT +,D,W,SLP,AT,DPT,SST,RH,MISS,ILL) C TRIMMING FLAGS ATTACHMENT IMPLICIT INTEGER(A-E,G-Z) INTEGER*8 ATT2 COMMON /ATT/AL(15),AD(255,15) +,SUP(255),SUPLEN,ERRNUM(51),ERRLEN(51),ERR(15,51),NERR PARAMETER (NFDSUL=5+3*6) DIMENSION CDSUL(NFDSUL),FDSUL(NFDSUL) DATA FDSUL/NFDSUL*-99999./ REAL D,W,SLP,AT,DPT,SST REAL FF,SCALAR,UW,VW,RH,RELHUM LOGICAL B2QXY LOGICAL NCDC LOGICAL MEDS SAVE FF()=SCALAR(D,W,MISS*1.,AD( 2,1)) C! RH()=RELHUM(AT,DPT,MISS*1.,ILL*1.) NCDC(B1,B2)=AD(B1,1).GE.B2 .AND. AD(B1,1).NE.10 MEDS(B1,B2)=SUPLEN.GE.B1 .AND.(SUP(B1).EQ.B2 + .OR. B2.EQ.ICHAR('3') .AND. SUP(B1).EQ.ICHAR('4')) C IF (.NOT.B2QXY(QB10(B10),LON,LAT,B2)) STOP 'B2QXY' C IF ( NINT(FDSUL(3)).NE.B2 +.OR. NINT(FDSUL(2)).NE.MO +.OR. PERIOD(NINT(FDSUL(4))).NE.PERIOD(YR)) THEN CALL GETDSUL(B10,MO,B2,YR,CDSUL,FDSUL,MISS*1.) IF ( NINT(FDSUL(3)).NE.B2 + .OR. NINT(FDSUL(2)).NE.MO + .OR. PERIOD(NINT(FDSUL(4))).NE.PERIOD(YR)) STOP 'GETDSUL' ENDIF C C MISCELLANEOUS FIELDS AD( 1,2)=B2 AD( 2,2)=NIGHT(MO,HR,LON,B2,MISS,ILL) C C COMPOSITE QC FLAGS (NCDC/MEDS) DO 190 I=9,23 AD( I,2)=0 190 CONTINUE IF (NCDC( 1,7)) AD( 9,2)=1 IF (NCDC(10,8)) AD(10,2)=1 IF (NCDC( 7,8)) AD(11,2)=1 IF (NINT(FF()).EQ.MISS) THEN IF (NINT(W).NE.MISS .AND. NINT(W).NE.ILL) + AD(12,2)=1 ENDIF IF (NCDC( 6,8)) AD(13,2)=1 IF (NCDC( 8,8) +.OR.NCDC( 9,8)) AD(14,2)=1 IF (DCK.EQ.714) THEN C CALL GETSUP(AL,AD,SUP,SUPLEN) IF (MEDS( 41,ICHAR('3'))) AD( 9,2)=AD( 9,2)+2 IF (MEDS( 73,ICHAR('3'))) AD(10,2)=AD(10,2)+2 IF (MEDS( 83,ICHAR('3'))) AD(11,2)=AD(11,2)+2 IF (MEDS(103,ICHAR('3')) + .OR.MEDS(113,ICHAR('3'))) AD(12,2)=AD(12,2)+2 IF (MEDS( 93,ICHAR('3'))) AD(13,2)=AD(13,2)+2 IF (MEDS(123,ICHAR('3'))) AD(14,2)=AD(14,2)+2 ENDIF C C TRIMMING FLAGS AD( 3,2)=TRMFLG(1,SST ,CDSUL(6),FDSUL(6),MISS,ILL) AD( 4,2)=TRMFLG(2,AT ,CDSUL(6),FDSUL(6),MISS,ILL) AD( 5,2)=TRMFLG(3,UW(),CDSUL(6),FDSUL(6),MISS,ILL) AD( 6,2)=TRMFLG(4,VW(),CDSUL(6),FDSUL(6),MISS,ILL) AD( 7,2)=TRMFLG(5,SLP ,CDSUL(6),FDSUL(6),MISS,ILL) C! AD( 8,2)=TRMFLG(6,RH(),CDSUL(6),FDSUL(6),MISS,ILL) AD( 8,2)=TRMFLG(6,RH,CDSUL(6),FDSUL(6),MISS,ILL) IF (DCK.EQ.714) THEN IF (AD( 3,2).EQ.12 .AND. MEDS( 73,ICHAR('1'))) AD( 3,2)=11 IF (AD( 7,2).EQ.12 .AND. MEDS( 93,ICHAR('1'))) AD( 7,2)=11 ENDIF C C COMPOSITE QC FLAGS (NCDC-ONLY) IF (NCDC( 4,2)) AD(15,2)=1 IF (NCDC(11,2)) AD(16,2)=1 IF (NCDC(12,2)) AD(17,2)=1 IF (NCDC( 4,4)) AD(15,2)=2 IF (NCDC(11,4)) AD(16,2)=2 IF (NCDC(12,4)) AD(17,2)=2 IF (NCDC( 4,7)) AD(15,2)=3 IF (NCDC(11,7)) AD(16,2)=3 IF (NCDC(12,7)) AD(17,2)=3 C C LANDLOCK FLAG IF (CDSUL(6).EQ.65534) AD(18,2)=1 C C EXCLUSION FLAGS IF (PT.EQ.13 +.OR. PT.EQ.14 .OR. PT.EQ.16 +.OR. PT.EQ.6 .AND. DCK.GE.876 .AND. DCK.LE.883 +.AND. (HR.EQ.MISS .OR. HR.EQ.ILL +.OR. MOD(HR,100).EQ.50 .OR. MOD((HR+50)/100,3).NE.0)) THEN AD(19,2)=1 AD(20,2)=1 AD(21,2)=1 AD(22,2)=1 AD(23,2)=1 ELSE IF (PT.EQ.7 .OR. SID.EQ.70 .OR. SID.EQ.71) THEN AD(21,2)=1 ENDIF C CALL PBYTE (ATT2,AD( 1,2), 0,14) CALL PBYTE (ATT2,AD( 2,2),14, 2) CALL PBYTES(ATT2,AD( 3,2),16, 4,0,6) CALL PBYTES(ATT2,AD( 9,2),40, 2,0,9) CALL PBYTES(ATT2,AD(18,2),58, 1,0,6) END C$TRMFLG----------------------------------------------------------------3456789 INTEGER FUNCTION TRMFLG(J,FTRUE,CDSUL2,FDSUL2,MISS,ILL) C TRIMMING FLAG IMPLICIT INTEGER(A-E,G-Z) DIMENSION CDSUL2(3,6),FDSUL2(3,6) DIMENSION FTRUEL(6),FTRUEU(6) DATA (FTRUEL(I),FTRUEU(I),I=1,2)/-5.,40. ,-88.,58./ REAL LOWER(6),UPPER(6) DATA (LOWER(I),UPPER(I),I=1,6) +/-3.,40. ,-50.,50. ,-50.,50. ,-50.,50. ,920.,1060. ,0.,100./ C IF(NINT(FTRUE).EQ.MISS)THEN TRMFLG=15 ELSE IF(NINT(FTRUE).EQ.ILL +.OR. J.LE.2.AND.(FTRUE.LT.FTRUEL(J).OR.FTRUE.GT.FTRUEU(J)))THEN TRMFLG=14 ELSE IF(CDSUL2(1,J).EQ.65534)THEN TRMFLG=13 ELSE IF(CDSUL2(1,J).EQ.0)THEN TRMFLG=12 ELSE IF(FTRUE.GT.FDSUL2(2,J)+(FDSUL2(3,J)-FDSUL2(2,J))*4.5/3.5 +.OR. FTRUE.GT.UPPER(J))THEN TRMFLG=7 ELSE IF(FTRUE.LT.FDSUL2(2,J)-(FDSUL2(2,J)-FDSUL2(1,J))*4.5/3.5 +.OR. FTRUE.LT.LOWER(J))THEN TRMFLG=6 ELSE IF(FTRUE.GT.FDSUL2(3,J))THEN TRMFLG=5 ELSE IF(FTRUE.LT.FDSUL2(1,J))THEN TRMFLG=4 ELSE IF(FTRUE.GT.FDSUL2(2,J)+(FDSUL2(3,J)-FDSUL2(2,J))*.8)THEN TRMFLG=3 ELSE IF(FTRUE.LT.FDSUL2(2,J)-(FDSUL2(2,J)-FDSUL2(1,J))*.8)THEN TRMFLG=2 ELSE TRMFLG=1 ENDIF END C$PBYTE-----------------------------------------------------------------3456789 SUBROUTINE PBYTE(P,U,Q,B) IMPLICIT INTEGER(A-Z) IF(U.LT.0 .OR. U.GT.2**B-1)STOP 'SBYTE' CALL SBYTE(P,U,Q,B) END C$PBYTES----------------------------------------------------------------3456789 SUBROUTINE PBYTES(P,U,Q,B,S,N) IMPLICIT INTEGER(A-Z) DIMENSION U(*) DO 1 I=1,N 1 IF(U(I).LT.0 .OR. U(I).GT.2**B-1)STOP 'SBYTES' CALL SBYTES(P,U,Q,B,S,N) END cQC.17#############################################NOAA Lubker 03/08/18 subroutine DUPQC(qltycd,attm1,xsid,xyr,xmnth,xday,xhr,xy,xwddir c***** quality control flags and quality code +,xwdspd,xvis,xpreswx,xpastwx,xpress,xdryblb,xwetblb,xdewpt +,xseatmp,xn,xnh,xcl,xh,xcm,xch,xwvdir,xwvper,xwvhgt,xswldir +,xswlper,xswlhgt,xa,xppp,lmr) c c -----------revision history--------------------------------------- c level author date description c ===== ====== ========== ==================== c .11 sl 85/04/11. qcfix.2 routines pg5 pg7 and conditions c for calling pg9 (w/ wddir = MISS, ILL) c substituted in; c routine PRINVN prints QC level. c .12 sl 92/03/11. a and ppp dummy arguments; c hr in hundredths; c call ATTMS removed. c .14 sl 94/06/13. y in hundredths. c .16 sl 02/04/16. buffer in replaced with direct access c read. c .17 sl 03/08/18. check for leap century added; check for c wdspd = MISS added for if ILL = MISS. c c ------------------------------------------------------------------ c c*********************************************************************** c c input- c sid thru swlhgt c coded+base c lmr c packed c c output- c qltycd c coded c c attm1 c quality control flags packed (56 bits) c 0=missing, 1=r, 2=a, 3=b, 4=j, 5=k, 6=l, 7=m, 8=n, 9=q, 10=s c c common- c lndlck(100) c msq1 0-msq1 99, 1=msq1 landlocked, 0=msq1 not landlocked c c m4pt8(4,12,5),p4pt8(4,12,5),m5pt8(4,12,5),p5pt8(4,12,5) c row 1=msq5 1 column 1=Jan level 1=pressure array m4pt8=m-4.8sd c row 2=msq5 2 column 2=Feb level 2=dry bulb array p4pt8=m+4.8sd c row 3=msq5 3 column 3=Mar level 3=wet bulb array m5pt8=m-5.8sd c row 4=msq5 4 column 4=Apr level 4=dew point array p5pt8=m+5.8sd c etc level 5=sea temp c c*********************************************************************** implicit integer(a-z) integer*8 qltycd,attm1 common /qc/msq1,msq5,invnqc(43),invnf(14,11),MISS,ILL,lndlck(100) +,m4pt8(4,12,5),p4pt8(4,12,5),m5pt8(4,12,5),p5pt8(4,12,5) +/qc0/shipf,windf,visf,prswxf,pstwxf,pressf,dryf,wetf,dewf,seaf +,cloudf,seawvf,swlwvf,ptendf,valu(10) dimension lmr(*) C data valu/0,1,1,2,2,2,3,3,3,3/,invnqc/43*0/,invnf/154*0/ sid =xsid yr =xyr mnth =xmnth day =xday hr =xhr y =xy wddir =xwddir wdspd =xwdspd vis =xvis preswx=xpreswx pastwx=xpastwx press =xpress dryblb=xdryblb wetblb=xwetblb dewpt =xdewpt seatmp=xseatmp n =xn nh =xnh cl =xcl h =xh cm =xcm ch =xch wvdir =xwvdir wvper =xwvper wvhgt =xwvhgt swldir=xswldir swlper=xswlper swlhgt=xswlhgt a =xa ppp =xppp shipf=1 windf=1 visf=1 prswxf=1 pstwxf=1 pressf=1 dryf=1 wetf=1 dewf=1 seaf=1 cloudf=1 seawvf=1 swlwvf=1 ptendf=1 c*********************************************************************** c***** ship position flag(m) if(lndlck(msq1+1).eq.1)goto 1006 if(day.lt.1.or.day.gt.31)goto 1006 if(hr.lt.0.or.hr.gt.2399)goto 1006 if(mnth.ne.4.and.mnth.ne.6.and.mnth.ne.9.and.mnth.ne.11)goto 1004 if(day.gt.30)goto 1006 goto 1008 1004 if(day.le.28.or.mnth.ne.2)goto 1008 if(mod(yr,4).ne.0.or.mod(yr,100).eq.0.and.mod(yr,400).ne.0) +goto 1006 if(day.le.29)goto 1008 1006 shipf=7 c*********************************************************************** c***** a, ppp, and illegal fields from attachments 1008 continue c1008 call ATTMS(sid,wddir,wdspd,vis,preswx,pastwx,press,dryblb c +,wetblb,dewpt,seatmp,n,nh,cl,h,cm,ch,wvdir,wvper,wvhgt,swldir c +,swlper,swlhgt,a,ppp,lmr) c*********************************************************************** c***** pressure flag(s,q,m,l) if(press.ne.MISS)goto 1010 pressf=10 goto 1060 1010 if(press.ne.ILL)goto 1020 pressf=7 goto 1060 1020 if(press.ge.9200.and.press.le.10600)goto 1030 pressf=9 goto 1060 1030 if(m5pt8(msq5,mnth,1).eq.MISS)goto 1060 if(press.ge.m5pt8(msq5,mnth,1).and.press.le.p5pt8(msq5,mnth,1)) +goto 1040 if(iabs(y).lt.3000.and.(wdspd.gt.154.and.wdspd.le.1029))goto 1050 pressf=9 goto 1060 1040 if(press.ge.m4pt8(msq5,mnth,1).and.press.le.p4pt8(msq5,mnth,1)) +goto 1060 1050 pressf=6 c*********************************************************************** c***** pressure tendency flag(s,m,k) 1060 if(a.ne.MISS.or.ppp.ne.MISS)goto 1070 ptendf=10 goto 1090 1070 if((a.ne.MISS.and.a.ne.ILL.and.ppp.ne.MISS.and.ppp.ne.ILL).and. +((ppp.eq.0.and.(a.eq.0.or.a.eq.4.or.a.eq.5)).or.(ppp.ne.0.and.a +.ne.4)))goto 1080 ptendf=7 goto 1090 1080 if(ppp.gt.150)ptendf=5 c*********************************************************************** c***** clouds flag(s,n,j,b), present weather flag(j,b) 1090 if(n.ne.MISS.or.nh.ne.MISS.or.cl.ne.MISS.or.h.ne.MISS.or.cm.ne. +MISS.or.ch.ne.MISS)goto 2010 cloudf=10 goto 2160 c 2010 if(n.ne.MISS.and.n.ne.ILL)goto 2050 if(nh.ne.8)goto 2020 cloudf=3 if(preswx.eq.43.or.preswx.eq.45.or.preswx.eq.47.or.preswx.eq.49) +prswxf=3 goto 2160 2020 if(nh.ne.9)goto 2030 cloudf=3 if(preswx.eq.42.or.preswx.eq.44.or.preswx.eq.46.or.preswx.eq.48) +prswxf=3 goto 2160 2030 if(nh.ne.MISS.and.nh.ne.ILL)goto 2040 cloudf=8 goto 2160 2040 cloudf=4 goto 2160 c 2050 if(n.ne.8)goto 2070 if(nh.ne.9)goto 2060 cloudf=3 goto 2160 2060 if(preswx.eq.43.or.preswx.eq.45.or.preswx.eq.47.or.preswx.eq.49) +prswxf=3 goto 2160 c 2070 if(n.ge.nh.or.nh.eq.MISS.or.nh.eq.ILL)goto 2090 if(nh.ne.8.and.nh.ne.9)goto 2080 cloudf=3 goto 2160 2080 cloudf=4 goto 2160 c 2090 if(n.ne.9)goto 2140 if(nh.ne.8)goto 2100 cloudf=3 goto 2160 2100 if(preswx.lt.0.or.preswx.gt.3)goto 2110 if(vis.lt.97.or.vis.gt.99)prswxf=4 cloudf=4 goto 2160 2110 if(nh.ne.9)goto 2120 if((preswx.lt.13.or.preswx.gt.29).and.preswx.ne.40)goto 2114 prswxf=4 goto 2117 2114 if(preswx.eq.42.or.preswx.eq.44.or.preswx.eq.46.or.preswx.eq.48) +prswxf=3 2117 if(h.ne.0.and.h.ne.10.and.h.ne.MISS)cloudf=3 goto 2160 2120 if((preswx.ge.13.and.preswx.le.29).or.preswx.eq.40.or.preswx +.eq.42.or.preswx.eq.44.or.preswx.eq.46.or.preswx.eq.48) +goto 2130 cloudf=3 goto 2160 2130 cloudf=4 prswxf=4 goto 2160 c 2140 if(n.ne.0.or.nh.ne.0)goto 2150 if((preswx.ge.49.and.preswx.le.99).or.preswx.eq.43.or.preswx.eq.45 +.or.preswx.eq.47)prswxf=4 if(h.ne.9)cloudf=3 goto 2160 c 2150 if(preswx.eq.43.or.preswx.eq.45.or.preswx.eq.47.or.preswx.eq.49) +prswxf=3 c*********************************************************************** c***** present weather flag(s,m,l), visibility flag(s,m) 2160 if(preswx.ne.MISS)goto 3010 prswxf=10 goto 3030 3010 if(preswx.ne.ILL)goto 3020 prswxf=7 goto 3030 3020 if(iabs(y).lt.2000.and.((preswx.ge.66.and.preswx.le.79).or.(preswx +.ge.83.and.preswx.le.88).or.(preswx.ge.36.and.preswx.le.39).or. +(preswx.ge.22.and.preswx.le.24).or.preswx.eq.26.or.preswx.eq.48 +.or.preswx.eq.49.or.preswx.eq.56.or.preswx.eq.57.or.preswx.eq.93 +.or.preswx.eq.94))prswxf=6 c 3030 if(vis.ne.MISS)goto 3040 visf=10 goto 3050 3040 if(vis.ne.ILL)goto 3050 visf=7 c*********************************************************************** c***** present weather flag(j,b), wind flag(j) 3050 if(prswxf.ge.4.or.visf.ne.1)goto 4110 c if(vis.gt.93)goto 4020 if(preswx.gt.3.and.preswx.ne.10)goto 4010 prswxf=4 goto 4110 4010 if((wdspd.ge.0.and.wdspd.lt.46).and.((preswx.ge.30.and.preswx.le. +39).or.preswx.eq.7))windf=4 goto 4110 c 4020 if(vis.gt.94)goto 4040 if(preswx.gt.3)goto 4030 prswxf=4 goto 4110 4030 if((wdspd.ge.0.and.wdspd.lt.46).and.((preswx.ge.30.and.preswx.le. +39).or.preswx.eq.7))windf=4 goto 4110 c 4040 if(vis.gt.96)goto 4080 if(preswx.lt.42.or.preswx.gt.49)goto 4050 prswxf=3 goto 4110 4050 if(preswx.gt.3)goto 4060 prswxf=4 goto 4110 4060 if(wdspd.lt.0.or.wdspd.ge.46)goto 4110 if(preswx.lt.33.or.preswx.gt.35)goto 4070 prswxf=4 goto 4110 4070 if(preswx.eq.36.or.preswx.eq.37)windf=4 goto 4110 c 4080 if(vis.gt.98)goto 4100 if((preswx.lt.42.or.preswx.gt.49).and.(preswx.lt.52.or.preswx.gt. +55).and.(preswx.lt.62.or.preswx.gt.65).and.(preswx.lt.72.or. +preswx.gt.75).and.preswx.ne.59.and.preswx.ne.69) +goto 4090 prswxf=3 goto 4110 4090 if((wdspd.ge.0.and.wdspd.lt.46).and.(preswx.ge.33.and.preswx.le.35 +))prswxf=4 goto 4110 c 4100 if((preswx.ge.30.and.preswx.le.99).or.(preswx.ge.4.and.preswx.le. +10))prswxf=4 c*********************************************************************** c***** dry bulb flag(s,q,n,m,l), wet bulb flag(s,q,n,m,l,b), dew point c flag(s,q,n,m,l,b) 4110 if(dryblb.ne.MISS)goto 5010 dryf=10 goto 5040 5010 if(dryblb.ne.ILL)goto 5020 dryf=7 goto 5040 5020 if(dryblb.ge.-491.and.dryblb.le.500)goto 5025 dryf=9 wetf=9 dewf=9 goto 5150 5025 if(m5pt8(msq5,mnth,2).eq.MISS)goto 5040 if(dryblb.ge.m5pt8(msq5,mnth,2).and.dryblb.le.p5pt8(msq5,mnth,2)) +goto 5030 dryf=9 wetf=9 dewf=9 goto 5150 5030 if(dryblb.ge.m4pt8(msq5,mnth,2).and.dryblb.le.p4pt8(msq5,mnth,2)) +goto 5040 dryf=6 wetf=6 dewf=6 goto 5150 c 5040 if(wetblb.ne.MISS)goto 5050 wetf=10 goto 5080 5050 if(wetblb.ne.ILL)goto 5060 wetf=7 goto 5080 5060 if(m5pt8(msq5,mnth,3).eq.MISS)goto 5080 if(wetblb.ge.m5pt8(msq5,mnth,3).and.wetblb.le.p5pt8(msq5,mnth,3)) +goto 5070 wetf=9 dewf=9 goto 5170 5070 if(wetblb.ge.m4pt8(msq5,mnth,3).and.wetblb.le.p4pt8(msq5,mnth,3)) +goto 5080 wetf=6 dewf=6 goto 5170 c 5080 if(dewpt.ne.MISS)goto 5090 dewf=10 goto 5120 5090 if(dewpt.ne.ILL)goto 5100 dewf=7 goto 5120 5100 if(m5pt8(msq5,mnth,4).eq.MISS)goto 5120 if(dewpt.ge.m5pt8(msq5,mnth,4).and.dewpt.le.p5pt8(msq5,mnth,4)) +goto 5110 dewf=9 goto 5190 5110 if(dewpt.ge.m4pt8(msq5,mnth,4).and.dewpt.le.p4pt8(msq5,mnth,4)) +goto 5120 dewf=6 goto 5190 c 5120 if(dryf.ne.1)goto 5140 if(dewf.ne.1)goto 5130 if(dryblb+5.ge.dewpt)goto 5125 dryf=8 dewf=8 if(wetf.eq.1)wetf=8 goto 5190 5125 if(dryblb.ge.dewpt)goto 5130 dewf=3 if(wetf.ne.1)goto 5190 wetf=3 goto 5133 5130 if(wetf.ne.1.)goto 5190 5133 if(dryblb+5.ge.wetblb)goto 5135 dryf=8 wetf=8 if(dewf.le.3)dewf=8 goto 5190 5135 if(dryblb.ge.wetblb)goto 5145 wetf=3 if(dewf.gt.3)goto 5190 dewf=3 goto 5147 5140 if(wetf.ne.1)goto 5190 5145 if(dewf.gt.3)goto 5190 5147 if(wetblb+5.ge.dewpt)goto 5149 wetf=8 dewf=8 goto 5190 5149 if(wetblb.ge.dewpt)goto 5190 wetf=3 dewf=3 goto 5190 c 5150 if(wetblb.ne.MISS)goto 5160 wetf=10 goto 5170 5160 if(wetblb.eq.ILL)wetf=7 5170 if(dewpt.ne.MISS)goto 5180 dewf=10 goto 5190 5180 if(dewpt.eq.ILL)dewf=7 c*********************************************************************** c***** sea temperature flag(s,q,m,l) 5190 if(seatmp.ne.MISS)goto 6010 seaf=10 goto 6040 6010 if(seatmp.ne.ILL)goto 6020 seaf=7 goto 6040 6020 if(seatmp.ge.-28.and.seatmp.le.500)goto 6025 seaf=9 goto 6040 6025 if(m5pt8(msq5,mnth,5).eq.MISS)goto 6040 if(seatmp.ge.m5pt8(msq5,mnth,5).and.seatmp.le.p5pt8(msq5,mnth,5)) +goto 6030 seaf=9 goto 6040 6030 if(seatmp.lt.m4pt8(msq5,mnth,5).or.seatmp.gt.p4pt8(msq5,mnth,5)) +seaf=6 c*********************************************************************** c***** present weather flag(j,b), dry bulb flag(j) 6040 if(dryf.ne.1.or.prswxf.ne.1)goto 7060 c if(dryblb.ge.-22)goto 7020 if((preswx.lt.58.or.preswx.gt.65).and.(preswx.lt.50.or.preswx.gt. +55).and.(preswx.lt.83.or.preswx.gt.84))goto 7010 prswxf=3 goto 7060 7010 if((preswx.ge.89.and.preswx.le.99).or.(preswx.ge.80.and.preswx +.le.82))dryf=4 goto 7060 c 7020 if(dryblb.lt.22)goto 7060 if((preswx.lt.48.or.preswx.gt.49).and.(preswx.lt.56.or.preswx +.gt.57).and.(preswx.lt.66.or.preswx.gt.67))goto 7030 prswxf=3 goto 7060 c 7030 if(dryblb.le.80)goto 7060 if((preswx.lt.83.or.preswx.gt.86).and.(preswx.lt.68.or.preswx +.gt.69))goto 7040 prswxf=3 goto 7060 7040 if(preswx.lt.70.or.preswx.gt.75)goto 7050 dryf=4 goto 7060 7050 if(preswx.ge.36.and.preswx.le.39)prswxf=4 c*********************************************************************** c***** past weather flag(s,m,j) 7060 if(pastwx.ne.MISS)goto 7070 pstwxf=10 goto 7090 7070 if(pastwx.ne.ILL)goto 7080 pstwxf=7 goto 7090 7080 if(dryf.eq.1.and.dryblb.gt.140.and.pastwx.eq.7)pstwxf=4 c*********************************************************************** c***** wind flag(s,q,m,j,a), sea waves flag(a) 7090 if(wddir.ne.MISS.or.wdspd.ne.MISS)goto 8010 windf=10 goto 8130 8010 if(wdspd.ne.ILL.or.wdspd.eq.MISS)goto 8020 windf=7 goto 8130 8020 if(wdspd.le.1029.or.wdspd.eq.MISS)goto 8030 windf=9 goto 8130 c 8030 if(wddir.lt.1.or.wddir.gt.360)goto 8060 if(wdspd.ne.0)goto 8050 if(wvhgt.le.0.or.wvhgt.gt.99)goto 8040 windf=7 goto 8130 8040 if(windf.eq.1)windf=2 goto 8130 8050 if(wdspd.eq.MISS)windf=7 goto 8130 c 8060 if(wddir.ne.361)goto 8070 if(wdspd.eq.0.or.windf.ne.1)goto 8130 windf=2 if((wdspd.gt.31.and.wdspd.le.1029) +.and.(wvhgt.gt.0.and.wvhgt.le.99))seawvf=2 goto 8130 c 8070 if(wddir.ne.362)goto 8100 if(wdspd.ne.MISS)goto 8080 windf=7 goto 8130 8080 if(wdspd.le.31)goto 8090 windf=4 goto 8130 8090 if(wdspd.eq.0.and.windf.eq.1)windf=2 goto 8130 c 8100 if(wddir.ne.MISS.or.wdspd.gt.31)goto 8120 if(windf.eq.1)windf=2 goto 8130 c 8120 windf=7 c*********************************************************************** c***** sea waves flag(s,q,n,m,j,b), wind flag(j) 8130 if(wvdir.ne.MISS.or.yr.lt.1968 +.or.wddir.eq.MISS.or.wddir.eq.ILL.or.wddir.eq.361 +.or.(wvper.eq.MISS.and.wvhgt.eq.MISS))goto 9005 if(wddir.gt.361)goto 9000 wvdir=36 goto 9005 9000 wvdir=38 c 9005 if(wvdir.ne.MISS.or.wvper.ne.MISS.or.wvhgt.ne.MISS)goto 9010 seawvf=10 goto 9240 9010 if(wvhgt.le.70.or.wvhgt.gt.99)goto 9020 seawvf=9 goto 9240 9020 if(yr.ge.1968)goto 9040 call pre68(.true.,wvdir,wvper,wvhgt,seawvf +,yr,mnth,seatmp,seaf,wdspd,windf) goto 9240 c 9040 if((wvdir.lt.1.or.(wvdir.gt.36.and.wvdir.ne.38)) +.and.seawvf.ne.2)goto 9210 if(wvhgt.ne.MISS.and.wvhgt.ne.ILL)goto 9050 seawvf=7 goto 9240 9050 if(wvhgt.eq.0)goto 9180 if(wvper.ne.0.and.wvper.ne.MISS.and.wvper.ne.ILL)goto 9060 if(wvper.ne.0)seawvf=3 goto 9130 9060 if(wvhgt.gt.5)goto 9070 if(wvper.gt.9+2*wvhgt)seawvf=3 goto 9130 9070 if(wvhgt.gt.12)goto 9080 if(wvper.gt.21)seawvf=3 goto 9130 9080 if(wvhgt.gt.20)goto 9090 if(wvper.lt.6)seawvf=3 goto 9130 9090 if(wvhgt.gt.29)goto 9100 if(wvper.lt.8)seawvf=3 goto 9130 9100 if(wvhgt.gt.40)goto 9110 if(wvper.lt.10)seawvf=3 goto 9130 9110 if(wvhgt.gt.54)goto 9120 if(wvper.lt.12)seawvf=3 goto 9130 9120 if(wvper.lt.14)seawvf=3 9130 if(windf.ge.4.or.wdspd.eq.MISS)goto 9240 if(seatmp.lt.20.or.seaf.ne.1)goto 9160 if(wvhgt.ne.1)goto 9140 if(wdspd.gt.108)goto 9170 goto 9240 9140 if(wvhgt.ne.2)goto 9150 if(wdspd.gt.170)goto 9170 goto 9240 9150 if(wvhgt.ne.3)goto 9160 if(wdspd.gt.242)goto 9170 goto 9240 9160 if(wvhgt.le.9)goto 9240 if(wdspd.lt.21)goto 9170 if(wvhgt.le.11)goto 9240 if(wdspd.lt.57)goto 9170 if(wvhgt.le.15)goto 9240 if(wdspd.lt.113)goto 9170 if(wvhgt.le.24)goto 9240 if(wdspd.lt.175)goto 9170 if(wvhgt.le.36)goto 9240 if(wdspd.lt.247)goto 9170 goto 9240 9170 seawvf=4 windf=4 goto 9240 9180 if(wvdir.ne.38.and.wvper.ne.MISS)goto 9190 seawvf=3 goto 9240 9190 if(wvper.ge.0.and.wvper.le.9)goto 9200 seawvf=8 goto 9240 9200 if(seaf.ne.1.or.windf.ge.4.or.seatmp.lt.20.or.wdspd.eq.MISS +.or.wdspd.le.67)goto 9240 seawvf=4 windf=4 goto 9240 c 9210 if((wvdir.eq.0.or.wvdir.eq.MISS).and.wvhgt.eq.0)goto 9220 seawvf=7 goto 9240 9220 if((wvper.ge.0.and.wvper.le.5).or.wvper.eq.MISS)goto 9230 seawvf=8 goto 9240 9230 if(wvdir.eq.0.and.wvper.eq.0)goto 9240 seawvf=3 c*********************************************************************** c***** swell waves flag(s,q,n,m,j,b) 9240 if(swldir.ne.MISS.or.swlper.ne.MISS.or.swlhgt.ne.MISS)goto 10010 swlwvf=10 goto 10180 10010 if(swlhgt.le.70.or.swlhgt.gt.99)goto 10020 swlwvf=9 goto 10180 10020 if((swldir.ne.11.and.swldir.ne.38).or.swlper.ne.11.or.swlhgt.ne. +11)goto 10030 swlwvf=8 goto 10180 10030 if(yr.ge.1968)goto 10040 call pre68(.false.,swldir,swlper,swlhgt,swlwvf +,yr,mnth,seatmp,seaf,wdspd,windf) goto 10180 c 10040 if(swldir.lt.1.or.(swldir.gt.36.and.swldir.ne.38))goto 10130 if(swlhgt.ne.MISS.and.swlhgt.ne.ILL)goto 10050 swlwvf=7 goto 10180 10050 if(swlhgt.eq.0)goto 10100 if(swlper.ne.0.and.swlper.ne.MISS.and.swlper.ne.ILL)goto 10060 if(swlper.ne.0)swlwvf=3 goto 10180 10060 if(swlhgt.ne.1)goto 10070 if(swlper.ge.12)swlwvf=3 goto 10180 10070 if(swlhgt.ne.2)goto 10080 if(swlper.ge.13)swlwvf=3 goto 10180 10080 if(swlhgt.le.12)goto 10180 if(swlper.le.5)goto 10090 if(swlhgt.le.20)goto 10180 if(swlper.le.7)goto 10090 if(swlhgt.le.29)goto 10180 if(swlper.le.9)goto 10090 if(swlhgt.le.40)goto 10180 if(swlper.le.11)goto 10090 if(swlhgt.le.54)goto 10180 if(swlper.le.13)goto 10090 goto 10180 10090 swlwvf=4 goto 10180 10100 if(swldir.ne.38.and.swlper.ne.MISS)goto 10110 swlwvf=3 goto 10180 10110 if(swlper.ge.0.and.swlper.le.8)goto 10180 if(swlper.ge.10.and.swlper.le.14)goto 10120 swlwvf=8 goto 10180 10120 swlwvf=3 goto 10180 c 10130 if(swldir.ne.0.and.swldir.ne.MISS)goto 10150 if(swldir.ne.0)goto 10140 if(swlhgt.ne.0.and.swlhgt.ne.MISS)goto 10150 goto 10160 10140 if(swlhgt.eq.0)goto 10160 10150 swlwvf=7 goto 10180 10160 if((swlper.ge.0.and.swlper.le.5).or.(swlper.ge.10.and.swlper.le.14 +).or.swlper.eq.MISS)goto 10170 swlwvf=8 goto 10180 10170 if(swldir.ne.0.or.swlper.ne.0.or.swlhgt.ne.0)swlwvf=3 c 10180 qltycd=valu(shipf)+valu(windf)+valu(visf)+valu(prswxf) ++valu(pstwxf)+valu(pressf)+valu(dryf)+valu(wetf)+valu(dewf) ++valu(seaf)+valu(cloudf)+valu(seawvf)+valu(swlwvf)+valu(ptendf) ++1 invnqc(qltycd)=invnqc(qltycd)+1 invnf(1,shipf+1)=invnf(1,shipf+1)+1 invnf(2,windf+1)=invnf(2,windf+1)+1 invnf(3,visf+1)=invnf(3,visf+1)+1 invnf(4,prswxf+1)=invnf(4,prswxf+1)+1 invnf(5,pstwxf+1)=invnf(5,pstwxf+1)+1 invnf(6,pressf+1)=invnf(6,pressf+1)+1 invnf(7,dryf+1)=invnf(7,dryf+1)+1 invnf(8,wetf+1)=invnf(8,wetf+1)+1 invnf(9,dewf+1)=invnf(9,dewf+1)+1 invnf(10,seaf+1)=invnf(10,seaf+1)+1 invnf(11,cloudf+1)=invnf(11,cloudf+1)+1 invnf(12,seawvf+1)=invnf(12,seawvf+1)+1 invnf(13,swlwvf+1)=invnf(13,swlwvf+1)+1 invnf(14,ptendf+1)=invnf(14,ptendf+1)+1 c call sbytes(attm1,shipf,0,4,0,14) attm1=attm1*16+shipf attm1=attm1*16+windf attm1=attm1*16+visf attm1=attm1*16+prswxf attm1=attm1*16+pstwxf attm1=attm1*16+pressf attm1=attm1*16+dryf attm1=attm1*16+wetf attm1=attm1*16+dewf attm1=attm1*16+seaf attm1=attm1*16+cloudf attm1=attm1*16+seawvf attm1=attm1*16+swlwvf attm1=attm1*16+ptendf attm1=attm1*16*16 + +attm1/16/16/16/16/16/16/16/16/16/16/16/16/16/16 end c####################################################################### subroutine pre68(seawv,dir,per,hgt,wvf c***** sea waves and swell waves prior to 1968 +,yr,mnth,seatmp,seaf,wdspd,windf) implicit integer(a-z) logical seawv common /qc/msq1,msq5,invnqc(43),invnf(14,11),MISS,ILL,lndlck(100) +,m4pt8(4,12,5),p4pt8(4,12,5),m5pt8(4,12,5),p5pt8(4,12,5) c c***** waves flag(n,m,j,b), wind flag(j) if((dir.lt.1.or.dir.gt.38).and.wvf.ne.2)goto 1170 if(hgt.ne.MISS.and.hgt.ne.ILL)goto 1010 wvf=7 return 1010 if(hgt.eq.0)goto 1140 if(per.ne.0.and.per.ne.MISS.and.per.ne.ILL)goto 1020 if(per.ne.0)wvf=3 goto 1090 1020 if(hgt.gt.5)goto 1030 if(per.gt.9+2*hgt)wvf=3 goto 1090 1030 if(hgt.gt.12)goto 1040 if(per.gt.21)wvf=3 goto 1090 1040 if(hgt.gt.20)goto 1050 if(per.lt.6)wvf=3 goto 1090 1050 if(hgt.gt.29)goto 1060 if(per.lt.8)wvf=3 goto 1090 1060 if(hgt.gt.40)goto 1070 if(per.lt.10)wvf=3 goto 1090 1070 if(hgt.gt.54)goto 1080 if(per.lt.12)wvf=3 goto 1090 1080 if(per.lt.14)wvf=3 1090 if(.not.seawv.or.windf.ge.4.or.yr.lt.1963.or.(yr.eq.1963.and.mnth +.lt.7).or.wdspd.eq.MISS)return if(seatmp.lt.20.or.seaf.ne.1)goto 1120 if(hgt.ne.1)goto 1100 if(wdspd.gt.108)goto 1130 return 1100 if(hgt.ne.2)goto 1110 if(wdspd.gt.170)goto 1130 return 1110 if(hgt.ne.3)goto 1120 if(wdspd.gt.242)goto 1130 return 1120 if(hgt.le.9)return if(wdspd.lt.21)goto 1130 if(hgt.le.11)return if(wdspd.lt.57)goto 1130 if(hgt.le.15)return if(wdspd.lt.113)goto 1130 if(hgt.le.24)return if(wdspd.lt.175)goto 1130 if(hgt.le.36)return if(wdspd.lt.247)goto 1130 return 1130 wvf=4 windf=4 return 1140 if(dir.lt.37.and.per.ne.MISS)goto 1150 wvf=3 return 1150 if(per.ge.0.and.per.le.9)goto 1165 if(per.eq.21.or.per.eq.22)goto 1160 wvf=8 return 1160 wvf=3 1165 if(.not.seawv.or.windf.ge.4.or.seaf.ne.1.or.yr.lt.1963 +.or.(yr.eq.1963.and.mnth.lt.7).or.seatmp.lt.20 +.or.wdspd.le.67.or.wdspd.eq.MISS)return wvf=4 windf=4 return c 1170 if(dir.ne.0.and.dir.ne.MISS)goto 1190 if(dir.ne.0)goto 1180 if(hgt.ne.0.and.hgt.ne.MISS)goto 1190 goto 1200 1180 if(hgt.eq.0)goto 1200 1190 wvf=7 return 1200 if((per.ge.0.and.per.le.5).or.per.eq.21.or.per.eq.22 +.or.per.eq.MISS)goto 1210 wvf=8 return 1210 if(dir.ne.0.or.per.ne.0.or.hgt.ne.0)wvf=3 end c####################################################################### subroutine LNDLMT(box10) c***** landlock and limits for box10 c*********************************************************************** c c common- c lndlck(100) c msq1 0-msq1 99, 1=msq1 landlocked, 0=msq1 not landlocked c c m4pt8(4,12,5),p4pt8(4,12,5),m5pt8(4,12,5),p5pt8(4,12,5) c row 1=msq5 1 column 1=Jan level 1=pressure array m4pt8=m-4.8sd c row 2=msq5 2 column 2=Feb level 2=dry bulb array p4pt8=m+4.8sd c row 3=msq5 3 column 3=Mar level 3=wet bulb array m5pt8=m-5.8sd c row 4=msq5 4 column 4=Apr level 4=dew point array p5pt8=m+5.8sd c etc level 5=sea temp c c local- c buffer(725) c box10 (32-bit integer) c m(4,12,5),sd(4,12,5),n(4,12,5) c row 1=msq5 1 column 1=Jan level 1=pressure array m=mean (32-bit integer to hundredths) c row 2=msq5 2 column 2=Feb level 2=dry bulb array sd=standard deviation (32-bit integer to hundredths) c row 3=msq5 3 column 3=Mar level 3=wet bulb array n=number (32-bit integer) c row 4=msq5 4 column 4=Apr level 4=dew point c etc level 5=sea temp c landlock (packed into 4 32-bit words) c msq1 0-msq1 99, 1=msq1 landlocked, 0=msq1 not landlocked c c*********************************************************************** implicit integer(a-z) common /qc/msq1,msq5,invnqc(43),invnf(14,11),MISS,ILL,lndlck(100) +,m4pt8(4,12,5),p4pt8(4,12,5),m5pt8(4,12,5),p5pt8(4,12,5) dimension buffer(725),m(4,12,5),sd(4,12,5),n(4,12,5) equivalence (buffer(2),m),(buffer(242),sd),(buffer(482),n) logical open c***** read inquire(9,opened=open) if(.not.open) +open(9,access='direct',form='unformatted',recl=4*725 +,file='LNDLM2_32',convert='big_endian') read(9,rec=box10)buffer if(buffer(1).ne.box10)goto 800 c***** unpack landlock C call gbytes(buffer(722),lndlck,0,1,0,100) k=0 do j=1,100 lndlck(j)=ibits(buffer(722+k/32),31-mod(k,32),1) k=k+1 enddo c***** compute limits do 410 lev=1,5 do 410 col=1,12 do 410 row=1,4 if(n(row,col,lev).lt.25)goto 400 p5pt8(row,col,lev)=(m(row,col,lev)+5.8*sd(row,col,lev))/10. + +sign(.5,(m(row,col,lev)+5.8*sd(row,col,lev))/10.) m5pt8(row,col,lev)=(m(row,col,lev)-5.8*sd(row,col,lev))/10. + +sign(.5,(m(row,col,lev)-5.8*sd(row,col,lev))/10.) p4pt8(row,col,lev)=(m(row,col,lev)+4.8*sd(row,col,lev))/10. + +sign(.5,(m(row,col,lev)+4.8*sd(row,col,lev))/10.) m4pt8(row,col,lev)=(m(row,col,lev)-4.8*sd(row,col,lev))/10. + +sign(.5,(m(row,col,lev)-4.8*sd(row,col,lev))/10.) goto 410 400 p5pt8(row,col,lev)=MISS m5pt8(row,col,lev)=MISS p4pt8(row,col,lev)=MISS m4pt8(row,col,lev)=MISS 410 continue return c***** error 800 print 810,box10,buffer(1) 810 format(' subroutine lndlmt box10 ',i3 +/' landlock/limits file positioned at box10 ',i3) stop end C Software for calculation of moisture variables in COADS 14 Jun 1991 C C ============================================================================ C C The following is excerpted from "Comprehensive Ocean-Atmosphere Data Set; C Release 1," pg. A18. C ---------------------------------------------------------------------------- C C 4.4 Moisture Variables C C The derived moisture variables (Q, R, and QS) are computed using the C FORTRAN functions that are given in [10] and referenced as follows: C Q = SSH(P,A - DP) C R = HUM(A,A - DP) C QS = SSH(P,S) C Inside SSH the mixing ratio is approximated by function WMR. The method C of computing vapor pressure differs in the untrimmed and trimmed C summaries. Function ESLO was used in the untrimmed summaries. C Unfortunately, ESLO is unreliable at physically unrealistic conditions, C although tests have demonstrated that, at least, no R exceeded 100%. C Function ES was used instead in the trimmed summaries. These algorithms C were chosen because of their accuracy and computational efficiency. For C more detailed information including the original source of these C techniques see [10]. C C [10] Schlatter, T. W., and D. V. Baker, 1981: Algorithms for thermodynamic C calculations. NOAA/ERL PROFS Program Office, Boulder, CO, 34 pp. C C ============================================================================ C C The following text and code is extracted from a later version of [10], C which differs for these routines only in that some comment lines have C been updated. In addition, the original code used functions ESAT and C ESW in HUM and WMR, respectively. The COADS (trimmed) implementation C substituted ES for ESAT or ESW. C ---------------------------------------------------------------------------- C C These algorithms were collected, edited, commented, and tested by Thomas W. C Schlatter and Donald V. Baker from August to October 1981 in the PROFS Program C Office, NOAA Environmental Research Laboratories, Boulder, Colorado. Where C possible, credit has been given to the original author of the algorithm and a C reference provided. C FUNCTION ESLO(T) C INCLUDE 'LIB_DEV:[GUDOC]EDFVAXBOX.FOR/LIST' C Baker, Schlatter 17-MAY-1982 Original version. C THIS FUNCTION RETURNS THE SATURATION VAPOR PRESSURE OVER LIQUID C WATER ESLO (MILLIBARS) GIVEN THE TEMPERATURE T (CELSIUS). THE C FORMULA IS DUE TO LOWE, PAUL R.,1977: AN APPROXIMATING POLYNOMIAL C FOR THE COMPUTATION OF SATURATION VAPOR PRESSURE, JOURNAL OF APPLIED C METEOROLOGY, VOL 16, NO. 1 (JANUARY), PP. 100-103. C THE POLYNOMIAL COEFFICIENTS ARE A0 THROUGH A6. DATA A0,A1,A2,A3,A4,A5,A6 1 /6.107799961, 4.436518521E-01, 1.428945805E-02, 2 2.650648471E-04, 3.031240396E-06, 2.034080948E-08, 3 6.136820929E-11/ ES = A0+T*(A1+T*(A2+T*(A3+T*(A4+T*(A5+A6*T))))) IF (ES.LT.0.) ES = 0. ESLO = ES RETURN END FUNCTION ES(T) C THIS FUNCTION RETURNS THE SATURATION VAPOR PRESSURE ES (MB) OVER C LIQUID WATER GIVEN THE TEMPERATURE T (CELSIUS). THE FORMULA APPEARS C IN BOLTON, DAVID, 1980: "THE COMPUTATION OF EQUIVALENT POTENTIAL C TEMPERATURE," MONTHLY WEATHER REVIEW, VOL. 108, NO. 7 (JULY), C P. 1047, EQ.(10). THE QUOTED ACCURACY IS 0.3% OR BETTER FOR C -35 < T < 35C. C INCLUDE 'LIB_DEV:[GUDOC]EDFVAXBOX.FOR/LIST' C Baker, Schlatter 17-MAY-1982 Original version. C ES0 = SATURATION VAPOR PRESSURE OVER LIQUID WATER AT 0C DATA ES0/6.1121/ ES = ES0*EXP(17.67*T/(T+243.5)) RETURN END FUNCTION HUM(T,TD) C INCLUDE 'LIB_DEV:[GUDOC]EDFVAXBOX.FOR/LIST' C G.S. Stipanuk 1973 Original version. C Reference Stipanuk paper entitled: C "ALGORITHMS FOR GENERATING A SKEW-T, LOG P C DIAGRAM AND COMPUTING SELECTED METEOROLOGICAL C QUANTITIES." C ATMOSPHERIC SCIENCES LABORATORY C U.S. ARMY ELECTRONICS COMMAND C WHITE SANDS MISSILE RANGE, NEW MEXICO 88002 C 33 PAGES C Baker, Schlatter 17-MAY-1982 C THIS FUNCTION RETURNS RELATIVE HUMIDITY (%) GIVEN THE C TEMPERATURE T AND DEW POINT TD (CELSIUS). AS CALCULATED HERE, C RELATIVE HUMIDITY IS THE RATIO OF THE ACTUAL VAPOR PRESSURE TO C THE SATURATION VAPOR PRESSURE. HUM= 100.*(ES(TD)/ES(T)) RETURN END FUNCTION SSH(P,T) C INCLUDE 'LIB_DEV:[GUDOC]EDFVAXBOX.FOR/LIST' C Baker, Schlatter 17-MAY-1982 Original version. C THIS FUNCTION RETURNS SATURATION SPECIFIC HUMIDITY SSH (GRAMS OF C WATER VAPOR PER KILOGRAM OF MOIST AIR) GIVEN THE PRESSURE P C (MILLIBARS) AND THE TEMPERATURE T (CELSIUS). THE EQUATION IS GIVEN C IN STANDARD METEOROLOGICAL TEXTS. IF T IS DEW POINT (CELSIUS), THEN C SSH RETURNS THE ACTUAL SPECIFIC HUMIDITY. C COMPUTE THE DIMENSIONLESS MIXING RATIO. W = .001*WMR(P,T) C COMPUTE THE DIMENSIONLESS SATURATION SPECIFIC HUMIDITY. Q = W/(1.+W) SSH = 1000.*Q RETURN END FUNCTION WMR(P,T) C THIS FUNCTION APPROXIMATES THE MIXING RATIO WMR (GRAMS OF WATER C VAPOR PER KILOGRAM OF DRY AIR) GIVEN THE PRESSURE P (MB) AND THE C TEMPERATURE T (CELSIUS). THE FORMULA USED IS GIVEN ON P. 302 OF THE C SMITHSONIAN METEOROLOGICAL TABLES BY ROLAND LIST (6TH EDITION). C INCLUDE 'LIB_DEV:[GUDOC]EDFVAXBOX.FOR/LIST' C Baker, Schlatter 17-MAY-1982 Original version. C EPS = RATIO OF THE MEAN MOLECULAR WEIGHT OF WATER (18.016 G/MOLE) C TO THAT OF DRY AIR (28.966 G/MOLE) DATA EPS/0.62197/ C THE NEXT TWO LINES CONTAIN A FORMULA BY HERMAN WOBUS FOR THE C CORRECTION FACTOR WFW FOR THE DEPARTURE OF THE MIXTURE OF AIR C AND WATER VAPOR FROM THE IDEAL GAS LAW. THE FORMULA FITS VALUES C IN TABLE 89, P. 340 OF THE SMITHSONIAN METEOROLOGICAL TABLES, C BUT ONLY FOR TEMPERATURES AND PRESSURES NORMALLY ENCOUNTERED IN C IN THE ATMOSPHERE. X = 0.02*(T-12.5+7500./P) WFW = 1.+4.5E-06*P+1.4E-03*X*X FWESW = WFW*ES(T) R = EPS*FWESW/(P-FWESW) C CONVERT R FROM A DIMENSIONLESS RATIO TO GRAMS/KILOGRAM. WMR = 1000.*R RETURN END c-----------------------------------------------------------------------3456789 block data bdqc0 implicit integer(a-z) common /qc/msq1,msq5,invnqc(43),invnf(14,11),MISS,ILL,lndlck(100) +,m4pt8(4,12,5),p4pt8(4,12,5),m5pt8(4,12,5),p5pt8(4,12,5) +/qc0/shipf,windf,visf,prswxf,pstwxf,pressf,dryf,wetf,dewf,seaf +,cloudf,seawvf,swlwvf,ptendf,valu(10) data valu/0,1,1,2,2,2,3,3,3,3/,invnqc/43*0/,invnf/154*0/ DATA MISS/-999999/,ILL/-999999/ end