#! /bin/csh -f # Script to drive program to read Copenhagen magnetometer data # received on floppies (see rdcphmag for the version which reads # the tapes): # (1) Change the read and write file names in the next two lines. # (2) Modify the write format to suit your needs. Look for # occurances of the string "WRITE (IOUN" # (3) Various diagnostics are written to stdout (fort.6), so it # may be best to execute each job by entering, for example, # # rdcphdisk > gdf.out onintr cleanup cat << EOF > in$$ MY940415.GDF my940415.asc EOF cat << EOF > src$$.f PROGRAM RDCPH C Sample program showing use of RDCMAG to read Copenhagen C magnetometer data CHARACTER*80 INAM , ONAM PARAMETER (INUN=11 , MSGUN=6 , RMISS=32767., + LLRB=432 , LWD=4 , LLRW=LLRB/LWD , + LBH=32 , LAH=40 , NDS=180) C Declarations for unpacking a logical record C LLRB = length of logical record in bytes C LWD = No bytes in a word C LLRW = length of logical record in words C LBH = length of binary header in bytes C LAH = length of ASCII header in bytes C NDS = number of data samples in a logical record C LBH = (expected) length of binary header in bytes C LAH = (expected) length of ASCII header in bytes C NDS = (expected) number of data samples in a logical record DIMENSION LREC(LLRW) , IVALS(NDS) , VALS(NDS) CHARACTER*1 IB(LLRB) EQUIVALENCE (LREC,IB) CHARACTER*3 STA CHARACTER*1 PAR CHARACTER*36 STR36 C Declarations for counting for each station, nlrecs and saving C first and last begin times PARAMETER (MXNSTA=100) CHARACTER*3 KSTA(MXNSTA) DIMENSION KNLR(MXNSTA) , KBD1(MXNSTA) , KBT1(MXNSTA), + KBDN(MXNSTA) , KBTN(MXNSTA) DATA KNTS /0/ C Get the input and output file names and open READ (*,'(A,/,A)',END=1000) INAM , ONAM OPEN (INUN,FILE=INAM,ACCESS='DIRECT',RECL=LLRB,STATUS='OLD') OPEN (IOUN,FILE=ONAM) WRITE (MSGUN,'('' RDCPHDISK reading: '',A,/, + '' writing: '',A)') INAM,ONAM NLRS = 0 NLRW = 0 NLR = 0 100 NLR = NLR + 1 C READ (INUN,REC=NLR,END=1000) LREC READ (INUN,REC=NLR,ERR=1000) LREC C Unpack (most of) the lrec C STA = IAGA Station code C PAR = Element "H,D,Z,X,Y,E,N" where H...Y have standard C meanings; E is local magnetic east and perpendicular C to H; and N is an integer identifying a component of C a system defined in the binary header. C DNP = Geographic North Pole distance in degrees. C GGLON = Geographic Longitude (East) in degrees. C COLATI = Invariant colatitude in degrees (not always supplied). C IBEGD = Date of start of record YYMMDD. C IBEGT = Time of start of record HHMMSS. C IDELT = Sample interval in seconds C IBASE = Tabular base in nT or Degrees (for element D). C IEC = Element code, where: 1 => D(.1deg), 2 => I(.1deg), C 3 => H(nT), 4 => F(nT), 5 => X(nT), 6 => Y(nT), C 7 => Z(nT), 8 => E(nT), 9 => H1(nT) perp to H2, C 10 => H2(nT) positive east magnetic latitude. C IRF = Record flag, where: 0 => okay, 1 => all data missing, C 2 => erroneous data, 4 => data rec manually corrected, C 9 => not a regular data record. C RL = Reference level, units match VALS(see IEC), missing=32767. C UN = Uncertainty or max daily drift, units match VALS, C missing=32767. C VALS = Logical record data values; a time series of PAR values, C starting at IBEGD,IBEGT with IDELT sample interval. C Compare anticipated pointer values with those unpacked: CALL GINT (IB( 1),1,LLRD) CALL GINT (IB( 3),1,LBHRD) CALL GINT (IB( 5),1,LAHRD) CALL GINT (IB(11),1,NDSRD) IF (LLRD .NE.LLRB .OR. LBHRD.NE.LBH .OR. LAHRD.NE.LAH .OR. + NDSRD.NE.NDS) THEN NWARN = NWARN + 1 IF (NWARN .LT. 3) + WRITE(MSGUN,'('' RDCDISK ignoring inconsistant header fields:'' + ,/,'' Read Expected'' + ,/,'' logical record length ='',2I8 + ,/,'' binary header length ='',2I8 + ,/,'' ASCII header length ='',2I8 + ,/,'' number of data samples ='',2I8)') + LLRD,LLRB , LBHRD,LBH , LAHRD,LAH , NDSRD,NDS ENDIF C Unpack the rest of the header CALL GINT (IB( 9),1,IDELT) CALL GINT (IB(13),1,IRL ) CALL GINT (IB(15),1,IUN ) IRF = ICHAR (IB(25)) ISC = ICHAR (IB(26)) IEC = ICHAR (IB(29)) CALL GSTR (IB(33), 3,STA) CALL GSTR (IB(36), 1,PAR) CALL GSTR (IB(37),36,STR36) READ (STR36,'(3F6.3,3I6)',ERR=150) DNP,GGLON,COLATI,IBEGD,IBEGT, + IBASE GO TO 160 150 NLRS = NLRS + 1 IF (NLRS .LT. 10) +WRITE (MSGUN,'(''RDCPH: trouble reformating: '',A,'' in lrec'', + I8,'' - SKIPPING THIS REC'')') STR36,NLR GO TO 100 C Unpack the data values into integers 160 CALL GINT (IB(73),NDS,IVALS) C Convert 2-s complement 16 bit integers to negatives, then scale C them as indicated by the 'SC' field SCLFCT = 1. IF (ISC .GT. 0 .AND. ISC .LT. 9) THEN SCLFCT = 2.**(3-ISC) ELSEIF (ISC .GT. 8) THEN SCLFCT = 10.**(10-ISC) ENDIF IF (IRL .GT. 32767) IRL = IRL - 65536 IF (IUN .GT. 32767) IUN = IUN - 65536 RL = REAL (IRL) * SCLFCT UN = REAL (IUN) * SCLFCT IF (ICHAR (IB(32)) .EQ. 0) THEN RL = RMISS UN = RMISS ENDIF DO 200 I=1,NDS IF (IVALS(I) .GT. 32767) IVALS(I) = IVALS(I) - 65536 VALS(I) = REAL (IVALS(I)) * SCLFCT 200 CONTINUE C Print the first few lrecs IF (PAR.EQ.'H' .AND. STR36(25:30) .EQ. 0) WRITE (MSGUN,'('' '',/, +'' LOGICAL RECORD:'',I8,/, +'' Station code: '',A3,/, +'' Element: '',A1,/, +'' North pole distance: '',A6,'' thousandths degrees'',/, +'' Geographic longitude: '',A6,'' thousandths degrees'',/, +'' Invariant colatitude: '',A6,'' thousandths degrees'',/, +'' Start date: '',A6,'' (yymmdd)'',/, +'' Start time: '',A6,'' (hhmmss)'',/, +'' Sample interval:'',I8,'' sec'',/, +''Extended element code:'',I8, + '' (1or2 => units=0.1 degrees, 3-10 => units=nT)'',/, +'' Record flag:'',I8,/, +'' D Tabular base: '',A6,'' nT or deg'',/, +'' Reference Level:'',F8.2,'' nT or deg (missing=32767.)'',/, +'' Uncertainty:'',F8.2,'' nT or deg (missing=32767.)'',/, +'' Data values:'',/,(10F8.2))') +NLR,STA,PAR,STR36(1:6),STR36(7:12), + STR36(13:18),STR36(19:24),STR36(25:30), + IDELT,IEC,IRF,STR36(31:36),RL,UN,VALS C +NLR,STA,PAR,DNP,GGLON, C +COLATI,IBEGD,IBEGT,IDELT,IEC,IRF,IBASE,RL,UN,VALS C Write this lrec to IOUN as a header record and integer records C (which must be interpreted by multiplying by SCLFCT as derived C from ISC; see the above example) NLRW = NLRW + 1 WRITE (IOUN,'(6I7,1X,A3,1X,A1,1X,A36,I7,/,(10I7))') IDELT,IRL, + IUN,IRF,ISC,IEC,STA,PAR,STR36, NDS,IVALS C Accumulate statistics for final summary print MATCHD = 0 DO 300 I=1,KNTS IF (STA .EQ. KSTA(I)) THEN MATCHD = 1 KNLR(I) = KNLR(I) + 1 KBDN(I) = IBEGD KBTN(I) = IBEGT ENDIF 300 CONTINUE IF (MATCHD .EQ. 0) THEN KNTS = KNTS + 1 IF (KNTS .GT. MXNSTA) STOP 'mxnsta' KSTA(KNTS) = STA KBD1(KNTS) = IBEGD KBT1(KNTS) = IBEGT KNLR(KNTS) = 1 ENDIF GO TO 100 1000 NLR = NLR - 1 WRITE (MSGUN,'('' '',/, + '' --- FIRST --- --- LAST ---'',/, + '' STA # lrecs IBEGD IBEGT IBEGD IBEGT'',/, + (1X,A3,5I8) )') + (KSTA(I),KNLR(I),KBD1(I),KBT1(I),KBDN(I),KBTN(I),I=1,KNTS) WRITE (MSGUN,'('' '',I8'' TOTAL written '',/, + '' '',I8'' total read '')') NLRW,NLR IF (NWARN .GT. 0) WRITE (MSGUN,'('' '',/,''There were '',I8, + '' records with inconsistent header fields'')') NWARN IF (NLRS .GT. 0) WRITE (MSGUN,'('' '',/,''There were '',I8, + '' records skipped due to reformat problems'')') NLRS END SUBROUTINE GINT (ISTR,NVALS,IVALS) C Assign adjacent bytes as 16-bit two's complement integers C ISTR = Input array of packed 2-byte integers C NVALS = Number of 2-byte integers to assign C IVALS = Assigned values (word aligned). CHARACTER*1 ISTR(*) DIMENSION IVALS(*) J = -1 DO 100 I=1,NVALS J = J + 2 IVALS(I) = ICHAR (ISTR(J))*256 + ICHAR (ISTR(J+1)) C IVALS(I) = ICHAR (ISTR(J+1))*256 + ICHAR (ISTR(J)) (byte reversed) IF (IVALS(I) .GT. 32767) IVALS(I) = IVALS(I) - 65536 100 CONTINUE RETURN END SUBROUTINE GSTR (ISTR,LSTR,STR) C Assign adjacent bytes to character variable CHARACTER*1 ISTR(*) CHARACTER*(*) STR STR = ' ' DO 100 I=1,LSTR 100 STR(I:I) = ISTR(I) RETURN END EOF f77 src$$.f -o exe$$ exe$$ < in$$ cleanup: rm -f *{$$}*