|
| 1 | + program diffpcp |
| 2 | +! |
| 3 | +! BASED ON DIFFPCP FROM VERF_PRECIP AND MODIFIED |
| 4 | +! |
| 5 | +!$$$ MAIN PROGRAM DOCUMENTATION BLOCK |
| 6 | +! . . . . |
| 7 | +! MAIN PROGRAM: ADDPCP |
| 8 | +! |
| 9 | +! Programmer: Mallory Row |
| 10 | +! |
| 11 | +! ABSTRACT: program to read in two precip files, then subtracts the two |
| 12 | +! |
| 13 | +! nam_2007060712_000_012 ! Input #1 |
| 14 | +! nam_2007060712_012_024 ! Input #2 |
| 15 | +! nam_2007060712_000_030 ! output file |
| 16 | +! |
| 17 | +! ATTRIBUTES: |
| 18 | +! LANGUAGE: Intel FORTRAN 90 |
| 19 | +! |
| 20 | +!$$$ |
| 21 | + integer jpds(200),jgds(200),kpds1(200),kpds2(200),kgds(200) |
| 22 | + integer kpdso(200),kgdso(200) |
| 23 | + parameter(ji=5000*2000) |
| 24 | + logical*1 bitdiff(ji),bit1(ji),bit2(ji) |
| 25 | + real diff(ji),pcp1(ji), pcp2(ji) |
| 26 | + character*200 infile1,infile2,prefx,outfile |
| 27 | + character*18 datstr |
| 28 | + character*3 timeflag |
| 29 | + character*2 orflag |
| 30 | + INTEGER IDAT(8),JDAT(8) |
| 31 | + REAL RINC(5) |
| 32 | + |
| 33 | +! Read agrument |
| 34 | +nargs = iargc() ! iargc() - number of arguments |
| 35 | +if (nargs.lt.3) then |
| 36 | + write(*,*)'usage : diffpcp infile1 infile2 outfile' |
| 37 | + stop |
| 38 | +endif |
| 39 | + |
| 40 | +call getarg(1,infile1) |
| 41 | +call getarg(2,infile2) |
| 42 | +call getarg(3,outfile) |
| 43 | + |
| 44 | + |
| 45 | +! |
| 46 | +! Get precip. ETA and most others have pds(5)=61, but AVN and MRF |
| 47 | +! have pds(5)=59 (precip rates), and in those cases we should |
| 48 | +! add up the files and then divde the result by the number of files |
| 49 | +! that contributed to the sum (i.e. calculate the average precip rate |
| 50 | +! during the entire period |
| 51 | +! |
| 52 | + |
| 53 | + jpds=-1 |
| 54 | + jgds=-1 |
| 55 | + call baopenr(11,infile1,ierr) |
| 56 | + if (ierr .ne. 0) write(6,*)" failed to open ", infile1 |
| 57 | + call getgb(11,0,ji,0,jpds,jgds,kf1,kr,kpds1,kgds,bit1,pcp1,iret1) |
| 58 | + write(6,*) 'getgb ', infile1, 'kf1= ',kf1, ' iret=', iret1 |
| 59 | + call baopenr(12,infile2,ierr) |
| 60 | + call getgb(12,0,ji,0,jpds,jgds,kf2,kr,kpds2,kgds,bit2,pcp2,iret2) |
| 61 | + write(6,*) 'getgb ', infile2, 'kf2= ',kf2, ' iret=', iret2 |
| 62 | + |
| 63 | +! |
| 64 | +! Calculate the difference (subtract file1 from file2): |
| 65 | +! |
| 66 | + if (iret1.ne.0 .or. iret2.ne.0) & |
| 67 | + STOP 'iret1 and/or iret2 not zero!' |
| 68 | + |
| 69 | +! |
| 70 | +! Check to see if the two files have the same length: |
| 71 | + if (kf1.ne.kf2) & |
| 72 | + STOP 'File lengths (kf1,kf2) differ! STOP.' |
| 73 | + |
| 74 | +! |
| 75 | +! Check to see if the two time stamps are identical: |
| 76 | + if (kpds1(21).ne.kpds2(21) .or. kpds1( 8).ne.kpds2( 8) .or. & |
| 77 | + kpds1( 9).ne.kpds2( 9) .or. kpds1(10).ne.kpds2(10) .or. & |
| 78 | + kpds1(11).ne.kpds2(11) .or. kpds1(12).ne.kpds2(12) .or. & |
| 79 | + kpds1(13).ne.kpds2(13)) then |
| 80 | + write(6,30) & |
| 81 | + (kpds1(21)*100+kpds1(8))/100-1, mod(kpds1(8),100), & |
| 82 | + kpds1(9),kpds1(10),kpds1(11),kpds1(12),kpds1(13), & |
| 83 | + (kpds2(21)*100+kpds2(8))/100-1, mod(kpds2(8),100), & |
| 84 | + kpds2(9),kpds2(10),kpds2(11),kpds2(12),kpds2(13) |
| 85 | + 30 format('Time stamps differ! Date1=',7i2.2,' Date2=',7i2.2, ' STOP') |
| 86 | + stop |
| 87 | + endif |
| 88 | + |
| 89 | + write(6,*) 'Infile1 Time range 1 : kpds1(14)=', kpds1(14) |
| 90 | + write(6,*) 'Infile1 Time range 2 : kpds1(15)=', kpds1(15) |
| 91 | + write(6,*) 'Infile2 Time range 1 : kpds2(14)=', kpds2(14) |
| 92 | + write(6,*) 'Infile2 Time range 2 : kpds2(15)=', kpds2(15) |
| 93 | + |
| 94 | +! Check to see if the two 'time range 1' are identical: |
| 95 | + if (kpds1(14).ne.kpds2(14)) then |
| 96 | + write(6,*) 'Time range 1 differ: kpds1(14)=', kpds1(14), & |
| 97 | + ' kpds2(14)=', kpds2(14),' STOP' |
| 98 | + stop |
| 99 | + endif |
| 100 | + |
| 101 | +! Check 'time range 2' to make sure that kpds2(15) > kpds1(15): |
| 102 | + |
| 103 | + if (kpds1(15).ge.kpds2(15)) then |
| 104 | + write(6,*) 'Time range 2 problem: kpds1(15)=', kpds1(15), & |
| 105 | + ' kpds2(15)=', kpds2(15),' STOP' |
| 106 | + stop |
| 107 | + endif |
| 108 | + |
| 109 | + kpdso=kpds2 |
| 110 | + kgdso=kgds |
| 111 | +! |
| 112 | + kpdso(14) = kpds1(15) |
| 113 | + |
| 114 | + do 40 N=1,kf1 |
| 115 | + bitdiff(N)=bit2(N).and.bit1(N) |
| 116 | + if (bitdiff(N)) then |
| 117 | + diff(N)=pcp2(N)-pcp1(N) |
| 118 | + write(54,54) n, diff(n), pcp2(n), pcp1(n) |
| 119 | + 54 format(i8,2x,3(3x,f8.3)) |
| 120 | + else |
| 121 | + diff(N)=0. |
| 122 | + endif |
| 123 | + 40 continue |
| 124 | +! |
| 125 | + 999 continue |
| 126 | +! |
| 127 | +! set unit decimal scale factor |
| 128 | + kpdso(22) = 1 |
| 129 | + |
| 130 | +! set output time stamp |
| 131 | + WRITE(DATSTR,50) (KPDSO(21)-1)*100+KPDSO(8),KPDSO(9), & |
| 132 | + KPDSO(10),KPDSO(11),KPDSO(14),KPDSO(15) |
| 133 | + 50 FORMAT(I4.4,3I2.2,'_',I3.3,'_',I3.3) |
| 134 | + OUTFILE = outfile |
| 135 | + CALL BAOPEN(51,OUTFILE,ierr) |
| 136 | + call putgb(51,kf1,kpdso,kgdso,bitdiff,diff,iret) |
| 137 | + if (iret.eq.0) then |
| 138 | + write(6,*) 'PUTGB successful, iret=', iret, 'for ', outfile |
| 139 | + else |
| 140 | + write(6,*) 'PUTGB failed! iret=', iret, 'for ', outfile |
| 141 | + endif |
| 142 | + CALL BACLOSE(51,ierr) |
| 143 | + CALL W3TAGE('DIFFPCP ') |
| 144 | +! |
| 145 | + stop |
| 146 | + |
| 147 | + end |
0 commit comments