The programs PApars, PApars.f, flags.f t2fit.f and tr2.f
Here is the listing. The explanation will follow some months or weeks from now below, after the === bar.
This one is a bit complicated. First we have the wrapper script PApars, then the program PApars.f but it is complied with flags.f, t2fit.f, and tr2.f as subroutine libraries of a sort. See do_comb_step2.sh for the actual compile…
Once again we see the trick of using parameter zero, the execution name, as a way to avoid hard coding the program to compile (presuming you don’t rename your script…) Oh, and notice that f77 is called out explicitly but that doesn’t really matter since the compile lines are commented out.
I’ve put it all in one listing since it looks like these bits are all really one big program.
First up, the PApars script:
if [[ $# -lt 1 ]] then echo "Usage: $0 source (e.g. GHCN.CL) radius(km) overlap_cond(20)" exit; fi rad=1000 ; if [[ $# -gt 1 ]] ; then rad=$2 ; fi lap=20 ; if [[ $# -gt 2 ]] ; then lap=$3 ; fi i="./ANN.dTs.$1" echo "inputfiles: $i.1-6 rural neighborhood radius:$rad km overlap_cond:$lap" if [[ ! -s $i.2 ]] then echo "input files $1.1-6 missing"; exit; fi n=1 while [[ $n -le 6 ]] do cp ${i}.$n fort.3$n (( n=$n + 1 )) done # f77 $0.f tr2.f t2fit.f ./$0.exe $rad $lap > $0.$1.$rad.$lap.log # f77 flags.f ; a.out >> $0.$1.$rad.$lap.log ./flags.exe >> $0.$1.$rad.$lap.log ; rm -f fort.78 ## Output files: echo "The following files were created:" fl=$0'.list' echo 'CCdStationID slope-l slope-r knee Yknee slope Ymid RMS RMSl 3-rur+urb ext.range flag' > $fl cat fort.2 >> $fl ; echo ${fl} echo " and for diagnostic purposes:" mv fort.66 $0.statn.log.${1}.$rad.$lap mv fort.77 $0.statn.use.${1}.$rad.$lap rm fort.2 mv fort.79 $0.noadj.stations.list ; echo $0.noadj.stations.list echo $0.$1.$rad.$lap.log # Clean-up rm -f fort.3[0-6] a.out $0.o t2fit.o tr2.o echo $0.statn.use.${1}.$rad.$lap echo $0.statn.log.${1}.$rad.$lap
Next is the FORTRAN program PApars.f listing:
C********************************************************************* C *** unit# C *** Input files: 31-36 ANN.dTs.GHCN.CL.1 ... ANN.dTs.GHCN.CL.6 C *** C *** Output file: 78 list of ID's of Urban stations with C *** homogenization info (text file) C *** Header line is added in subsequent step C********************************************************************* C**** C**** This program combines for each urban station the rural stations C**** within R=1000km and writes out parameters for broken line C**** approximations to the difference of urban and combined rural C**** annual anomaly time series. C**** C**** Input files: units 31,32,...,36 C**** Record 1: I1,INFO(2),...,INFO(8),I1L,TITLE header record C**** Record 2: IDATA(I1-->I1L),LT,LN,ID,HT,NAME,I2,I2L station 1 C**** Record 3: IDATA(I2-->I2L),LT,LN,ID,HT,NAME,I3,I3L station 2 C**** etc. NAME(31:31)=brightnessIndex 1=dark->3=bright C**** etc. NAME(32:32)=pop.flag R/S/U rur/sm.town/urban C**** etc. NAME(34:36)=country code C**** C**** IDATA(1) refers to year IYRBEG=INFO(6) C**** IYRM=INFO(4) is the max. length of an input time series, C**** C**** INFO(1),...,INFO(8) are 4-byte integers, C**** TITLE is an 80-byte character string, C**** INFO 2 = KQ (quantity flag, see below) C**** 3 = MAVG (time avg flag: 1 - 4 DJF - SON, 5 ANN, C**** 6 MONTHLY, 7 SEAS, 8 - 19 JAN - DEC ) C**** 4 = IYRM (length of each time record) C**** 5 = IYRM+4 (size of data record length) C**** 6 = IYRBEG (first year of full time record) C**** 7 = flag for missing data C**** 8 = flag for precipitation trace C**** Combining of rural stations C**** =========================== C**** Stations within Rngbr km of the urban center U contribute C**** to the mean at U with weight 1.- d/Rngbr (d = distance C**** between rural and urban station in km). To remove the station C**** bias, station data are shifted before combining them with the C**** current mean. The shift is such that the means over the time C**** period they have in common remains unchanged. If that common C**** period is less than 20(NCRIT) years, the station is disregarded. C**** To decrease that chance, stations are combined successively in C**** order of the length of their time record. C**** The homogeneity adjustment parameters C**** ===================================== C**** To minimize the impact of the natural local variability, only C**** that part of the combined rural record is actually used that is C**** supported by at least 3 stations, i.e. heads and tails of the C**** record that are based on only 1 or 2 stations are dropped. The C**** difference between that truncated combination and the non-rural C**** record is found and the best linear fit and best fit by a broken C**** line (with a variable "knee") to that difference series are found. C**** The parameters defining those 2 approximations are tabulated. C**** C**** Note: No attempt is made to find the longterm trends for urban C**** and rural combination separately; using the difference only C**** minimizes the impact of short term regional events that C**** affect both rural and urban stations, hence cancel out. C?*** Input parameters (# of input files, time period) PARAMETER (INM=6,KM0=1,IYRM0=KM0*(3000-1880+1),ITRIM=1) C?*** Earth radius, lower overlap limit, min.coverage of approx.range PARAMETER (REARTH=6375.,NCRIT=20,NRURM=3,XCRIT=2./3.) C?*** Work array sizes PARAMETER (NSTAM=8000,MSIZE=700000) CHARACTER*80 TITLEI,TITLEO,NAME*36,fn*7 C**** Input arrays DIMENSION IDATA(MSIZE),INFO(8+ITRIM),ITRL(14+ITRIM) REAL RDATA(MSIZE) EQUIVALENCE (ITRL(5),NAME),(IDATA,RDATA) C**** Station dependent arrays DIMENSION I1SU(NSTAM),ILSU(NSTAM),MFSU(NSTAM),IDU(NSTAM), * I1SR(NSTAM),ILSR(NSTAM),MFSR(NSTAM),IDR(NSTAM), * LEN(NSTAM),WTI(NSTAM),IORD(NSTAM),ISOFI(NSTAM) ! rural * ,nuseid(nstam),lenis(nstam) ! for diagn. purposes only real*8 r(nstam) REAL*4 CSLONU(NSTAM),CSLATU(NSTAM),SNLONU(NSTAM),SNLATU(NSTAM) REAL*4 CSLONR(NSTAM),CSLATR(NSTAM),SNLONR(NSTAM),SNLATR(NSTAM) CHARACTER*3 CC(NSTAM) C**** Output-related arrays (dim of output array = IYRM < IYRM0) DIMENSION AVG(IYRM0),DNEW(IYRM0),WT(IYRM0),URB(IYRM0),IWT(IYRM0) COMMON/LIMIT/XBAD,NLAP COMMON/FITCOM/W(900),X(900),F(900),CF(900),DF(900),ZFP(20),ZR(20) 1 ,FPAR(20),DELTAP,DFSTOP,X0,TMEAN,RMEAN,TFMEAN,RMSFIT 2 ,YR(900),TS(900),TSFIT(900),RMSP(20),KPFIT(20) 3 ,KPAR,NXY,NCP,NWCYCL,NITERS,NFIBNO,LISTIT,IDW REAL*8 W,X,F,CF,DF,ZFP,ZR,FPAR,DELTAP,DFSTOP,X0,TMEAN,RMEAN,TM3 REAL*8 TFMEAN,RMSFIT,YR,TS,TSFIT,RMSP PI180=DATAN(1.D0)/45.D0 RngbrF=1000. IF(IARGC().GT.0) THEN CALL GETARG(1,TITLEI) READ(TITLEI,*) IRngbr RngbrF=IRngbr END IF RngbrH=RngbrF/2 RBYRCF=REARTH/RngbrF RBYRCH=REARTH/RngbrH CSCRIF=COS(RngbrF/REARTH) CSCRIH=COS(RngbrH/REARTH) NLAP=NCRIT IF(IARGC().GT.1) THEN CALL GETARG(2,TITLEI) READ(TITLEI,*) NLAP END IF C**** Open the input files do in=31,30+inm fn='fort.xx' write(fn(6:7),'(i2)') in open(in,file=fn,form='unformatted') end do C**** Output text file open(78,file='fort.78',form='formatted') ! table of adj parameters C**** Diagnostic output files in addition to log on standard output open(66,file='fort.66',form='formatted') ! combination info open(77,file='fort.77',form='formatted') ! station usage stats open(79,file='fort.79',form='formatted') ! isolated urban stations C**** Read and interpret the header record READ (31) INFO,TITLEI KQ=INFO(2) IF(KQ.NE.1) THEN WRITE(6,*) ' PROGRAM NOT READY FOR QUANTITY ',KQ STOP 'INAPPROPRIATE QUANTITY' END IF C KM = the number of time frames per year KM=1 IF(INFO(3).EQ.6) KM=12 IF(INFO(3).EQ.7) KM=4 IF(KM.NE.KM0) STOP 'ERROR: CHANGE KM0' ML=INFO(4) iyrm=INFO(4) NYRSIN=INFO(4) IF(IYRM0.LT.INFO(4)) THEN WRITE(6,'('' INCREASE IYRM0 TO AT LEAST'',I5)') INFO(4) STOP 'ERROR: IYRM0 NOT OK' END IF IYOFF=INFO(6)-1 MBAD=INFO(7) LAST=INFO(7) XBAD=MBAD C**** TRACE=INFO(8) ! not needed here C**** Collect the station data needed ISU=0 ISR=0 I1Snow=1 DO 90 IN=1,INM REWIND 30+IN READ (30+IN) INFO MF=INFO(1) IF(ITRIM.GT.0) ML=INFO(8+ITRIM) 50 IF(MF.GE.LAST) GO TO 90 IF(I1Snow+ML-MF.GT.MSIZE) STOP 'ERROR: MSIZE TOO SMALL' CALL SREAD (30+IN,ML+1-MF,IDATA(I1Snow),ITRL,14+ITRIM) MFCUR=MF LENGTH=ML-MF+1 MF=ITRL(14) IF(ITRIM.GT.0) ML=ITRL(14+ITRIM) XLAT=.1*ITRL(1) XLON=.1*ITRL(2) IF(NAME(31:32).EQ.' R'.or.NAME(31:31).EQ.'1') THEN ISR=ISR+1 ! count rural stations IF(ISR.GT.NSTAM) STOP 'ERROR: NSTAM TOO SMALL' MFSR(ISR)=MFCUR I1SR(ISR)=I1Snow ILSR(ISR)=I1Snow+LENGTH-1 CSLATR(ISR)=COS(XLAT*PI180) SNLATR(ISR)=SIN(XLAT*PI180) CSLONR(ISR)=COS(XLON*PI180) SNLONR(ISR)=SIN(XLON*PI180) IDR(ISR)=ITRL(3) ELSE ISU=ISU+1 ! count non-rural stations IF(ISU.GT.NSTAM) STOP 'ERROR: NSTAM TOO SMALL' MFSU(ISU)=MFCUR I1SU(ISU)=I1Snow ILSU(ISU)=I1Snow+LENGTH-1 CSLATU(ISU)=COS(XLAT*PI180) SNLATU(ISU)=SIN(XLAT*PI180) CSLONU(ISU)=COS(XLON*PI180) SNLONU(ISU)=SIN(XLON*PI180) IDU(ISU)=ITRL(3) CC(ISU)=NAME(34:36) END IF I1Snow=I1Snow+LENGTH GO TO 50 90 CONTINUE NSTAu=ISU ! total number of bright/urban or dim/sm.town stations NSTAr=ISR ! total number of dark/rural stations LDTOT=I1Snow-1 ! total length of IDATA used write(*,*) 'number of rural/urban stations',NSTAr,NSTAu C**** Convert data to real numbers (ann avgs were multiplied by 10) DO 115 N=1,LDTOT IF(IDATA(N).EQ.MBAD) THEN RDATA(N)=XBAD ELSE RDATA(N)=.1*IDATA(N) ! .01C to .1C END IF 115 CONTINUE C**** Order the NSTAr rural stations according to length of time record DO 121 IS=1,NSTAr C 121 LEN(IS)=ILSR(IS)-I1SR(IS)+1 ! (counting gaps also) LEN(IS)=0 nuseid(is)=0 DO 120 M=I1SR(IS),ILSR(IS) IF(RDATA(M).EQ.XBAD) GO TO 120 LEN(IS)=LEN(IS)+1 120 CONTINUE write(*,*) 'rural station:',is,' id:',idr(is),' #ok',len(is) 121 continue do m=1,nstar r(m)=len(m)+(m-1d0)/dfloat(nstar) ! make "lengths" unique end do CALL SORT (IORD,NSTAr,R) do m=1,nstar len(m)=r(m) end do do 122 IS=1,NSTAr 122 write(*,*) 'rural station:',is,' id:',idr(iord(is)),' #ok',len(is) C**** Combine time series for rural stations around each urban station DO 200 NURB=1,NSTAu C** Loop over all rural stations in memory; try 1/2 the radius first CSCRIT=CSCRIH RBYRC=RBYRCH Rngbr=RngbrH 125 IS0=0 ! counter for rural neighbors DO 130 N=1,NSTAr ISR=IORD(N) IYU1=MFSU(NURB)+IYOFF-1 ! subtract 1 for a possible partial yr IYU2=IYU1+ILSU(NURB)-I1SU(NURB)+2 ! add 1 for partial year C** Find cosine of distance between center and station CSDBYR=SNLATR(ISR)*SNLATU(NURB)+CSLATR(ISR)*CSLATU(NURB)* * (CSLONR(ISR)*CSLONU(NURB)+SNLONR(ISR)*SNLONU(NURB)) X0=1950. IF(CSDBYR.LE.CSCRIT) GO TO 130 DBYRC=0. C** The arc is replaced by the smaller chord (to avoid using ACOS) IF(CSDBYR.LT.1.) DBYRC=RBYRC*SQRT(2.*(1.-CSDBYR)) IS0=IS0+1 WTI(IS0)=1.-DBYRC ISOFI(IS0)=ISR lenis(IS0)=len(n) 130 CONTINUE C**** C**** Combine the station data C**** DO 150 M=1,IYRM WT(M)=0. IWT(M)=0 URB(M)=XBAD 150 AVG(M)=XBAD IF(IS0.EQ.0) THEN IF(CSCRIT.EQ.CSCRIH) THEN write(*,*) 'Trying full radius',RngbrF CSCRIT=CSCRIF RBYRC=RBYRCF Rngbr=RngbrF GO TO 125 END IF WRITE(79,'('' NO RURAL NEIGHBORS FOR '',I9)') IDU(NURB) GO TO 200 END IF ioff=MFSU(NURB)-I1SU(NURB) DO 151 M=I1SU(NURB),ILSU(NURB) c if(M+IOFF.gt.IYRM0) stop 231 151 URB(M+IOFF)=RDATA(M) write(*,'(a,i9,a,i4,a,2i5,f9.0)') 'urb stnID:',idu(nurb),' # rur:' * ,is0,' ranges:',MFSU(NURB)+IYOFF,ILSU(NURB)+ioff+IYOFF,Rngbr C**** Start with the station with the longest time record IS=ISOFI(1) IOFF=MFSR(IS)-I1SR(IS) nuseid(is)=nuseid(is)+1 C**** First update of full data and weight arrays DO 160 M=I1SR(IS),ILSR(IS) c if(M+IOFF.gt.IYRM0) stop 244 AVG(M+IOFF)=RDATA(M) IF(RDATA(M).LT.XBAD) WT(M+IOFF)=WTI(1) IF(RDATA(M).LT.XBAD) IWT(M+IOFF)=1 160 CONTINUE write(*,'(a,i5,a,i4,i6,i10)') 'longest rur range:', * MFSR(IS)+IYOFF,'-',ILSR(IS)+ioff+IYOFF,LENis(1),idr(is) C**** Add in the remaining stations DO 190 I=2,IS0 IS=ISOFI(I) IOFF=MFSR(IS)-I1SR(IS) write(*,'(a,i5,a,i5,a,i4,i6,i10)')'add stn',i,' range:', * MFSR(IS)+IYOFF,'-',ILSR(IS)+ioff+IYOFF,LENis(i),idr(is) C**** Extend the new data into a full series DO 170 M=1,IYRM 170 DNEW(M)=XBAD DO 180 M=I1SR(IS),ILSR(IS) c if(M+IOFF.gt.IYRM0) stop 259 180 DNEW(M+IOFF)=RDATA(M) NF1=MFSR(IS) NL1=ILSR(IS)+IOFF C**** Shift new data, then combine them with current mean CALL CMBINE (AVG(1),WT,IWT, DNEW,NF1,NL1,WTI(I), * IDR(IS),NSM,NCOM) write(*,*) 'data added: ',nsm,' overlap:',ncom,' years' IF(NSM.EQ.0) GO TO 190 nuseid(is)=nuseid(is)+1 190 CONTINUE C**** Subtract urban station and call a curve fitting program IY1=1 191 NXY=0 N3=0 N3F=0 N3L=0 NXX=0 TMEAN=0 if(IY1.eq.1) * write(66,'(a,i9)')'year dTs-urban dTs-rural StnID=',idu(nurb) DO 195 IY=IY1,IYRM IF(AVG(IY).NE.XBAD.OR.URB(IY).NE.XBAD) NXX=NXX+1 IF(AVG(IY).EQ.XBAD.OR.URB(IY).EQ.XBAD) GO TO 194 if(IWT(IY).ge.NRURM) then N3L=IY N3=N3+1 if(N3F.eq.0) N3F=IY end if if(N3.le.0) go to 195 NXY=NXY+1 TS(NXY)=AVG(IY)-URB(IY) F(NXY)=TS(NXY) TMEAN=TMEAN+F(NXY) YR(NXY)=IY+IYOFF X(NXY)=YR(NXY)-X0 W(NXY)=1. if(IWT(IY).ge.NRURM) THEN NXY3=NXY TM3=TMEAN END IF 194 if(nxx.gt.0.and.IY1.eq.1) write(66,'(i4,2f10.2)') * IY+IYOFF,URB(IY),AVG(IY) 195 CONTINUE IF(N3.LT.NCRIT) THEN IF(RBYRC.NE.RBYRCF) THEN write(*,*) 'trying full radius',RngbrF RBYRC=RBYRCF CSCRIT=CSCRIF Rngbr=RngbrF GO TO 125 END IF WRITE(79,'(a3,i9.9,a13,i5,a15,i5,a50,a5)') CC(NURB),IDU(NURB), * ' good years:',N3,' total years:',N3L-N3F+1, * ' too little rural-neighbors-overlap - drop station',' 9999' GO TO 200 ELSE IF(FLOAT(N3).LT.XCRIT*(N3L-N3F+1.-.1)) THEN ! Not enough good years for the given range (<66%) ! The -0.1 is to prevent equality with potential uncertainty ! due to hardware rounding or compiler behaviour. ! Nick Barnes, Ravenbrook Limited, 2008-09-10 ! Try to save cases in which the gaps are in the early part: IY1=N3L-(N3-1)/XCRIT if(IY1 < N3F+1) IY1=N3F+1 ! avoid infinite loop WRITE(79,'(a3,i9.9,a17,i5,a1,i4)') * CC(NURB),IDU(NURB),' drop early years',1+IYOFF,'-',IY1-1+IYOFF GO TO 191 ELSE TMEAN=TM3/NXY3 NXY=NXY3 CALL GETFIT(15) ! 15=unit for debug output - currently not used C**** Find extended range IYXTND=NINT(N3/XCRIT)-(N3L-N3F+1) write(*,*) 'possible range increase',IYXTND,N3,N3L-N3F+1 n1x=N3F+IYOFF n2x=N3L+IYOFF IF(IYXTND.lt.0) stop 'impossible' if(IYXTND.gt.0) then LXEND=IYU2-(N3L+IYOFF) IF(IYXTND.le.LXEND) then n2x=n2x+LXEND else n1x=n1x-(IYXTND-LXEND) if(n1x.lt.IYU1) n1x=IYU1 n2x=IYU2 end if end if C**** Write out a table entry for the table of adjustment parameters write(78,'(a3,i9.9,2f9.3,i5,5f9.3,I5,a1,I4,i5,a1,i4)') CC(nurb), * idu(nurb),(fpar(i),i=1,2),nint(fpar(3)+X0),(fpar(i),i=4,6), * (rmsp(i),i=1,2),N3F+IYOFF,'-',N3L+IYOFF,N1X,'-',N2X END IF 200 CONTINUE nuse=0 do n=1,NSTAr if(nuseid(n).gt.0) then write(77,*) 'used station ',idr(n),nuseid(n),' times' nuse=nuse+1 end if end do write(77,*) nuse,' rural stations were used' STOP END SUBROUTINE CMBINE (AVG,WT,IWT, DNEW,NF1,NL1,WT1, ID,NSM,NCOM) C**** C**** Bias of new data is removed by subtracting the difference C**** over the common domain. Then the new data are averaged in. C**** COMMON/LIMIT/XBAD,NLAP DIMENSION AVG(*),DNEW(*),WT(*),IWT(*) C**** Loop over years NSM=0 C**** Find means over common domain to compute bias SUMN=0 SUM=0 NCOM=0 DO 10 N=NF1,NL1 IF(AVG(N).GE.XBAD.OR.DNEW(N).GE.XBAD) GO TO 10 NCOM=NCOM+1 SUM=SUM+AVG(N) SUMN=SUMN+DNEW(N) 10 CONTINUE IF(NCOM.LT.NLAP) RETURN BIAS=(SUM-SUMN)/FLOAT(NCOM) C**** Update period of valid data, averages and weights DO 20 N=NF1,NL1 IF(DNEW(N).GE.XBAD) GO TO 20 WTNEW=WT(N)+WT1 AVG(N)=(WT(N)*AVG(N)+WT1*(DNEW(N)+BIAS))/WTNEW WT(N)=WTNEW IWT(N)=IWT(N)+1 NSM=NSM+1 20 CONTINUE 50 CONTINUE RETURN END SUBROUTINE SORT (INDX,NDIM,RLNGTH) C**** Sorts INDEX and LNGTH such that rLNGTH becomes decreasing integer INDX(NDIM) real*8 RLNGTH(NDIM) DO 10 N=1,NDIM 10 INDX(N)=N DO 30 N=1,NDIM-1 C**** Find maximum of LNGTH(N),...LNGTH(NDIM) NLMAX=N DO 20 NN=N+1,NDIM IF(RLNGTH(NN).GT.RLNGTH(NLMAX)) NLMAX=NN 20 CONTINUE C**** Switch positions N and NLMAX RMAX=RLNGTH(NLMAX) RLNGTH(NLMAX)=RLNGTH(N) RLNGTH(N)=RMAX IMAX=INDX(NLMAX) INDX(NLMAX)=INDX(N) 30 INDX(N)=IMAX RETURN END SUBROUTINE SREAD (IN,LEN,IDATA,ITRL,NTRL) C**** Speed read routine for input records DIMENSION IDATA(LEN),ITRL(NTRL) READ (IN) IDATA,ITRL RETURN END
The program t2fit.f listing
subroutine getfit(nw) COMMON/FITCOM/W(900),X(900),F(900),CF(900),DF(900),ZFP(20),ZR(20) 1 ,FPAR(20),DELTAP,DFSTOP,X0,TMEAN,RMEAN,TFMEAN,RMSFIT 2 ,YR(900),TS(900),TSFIT(900),RMSP(20),KPFIT(20) 3 ,KPAR,NXY,NCP,NWCYCL,NITERS,NFIBNO,LISTIT,IDW REAL*8 W,X,F,CF,DF,ZFP,ZR,FPAR,DELTAP,DFSTOP,X0,TMEAN,RMEAN,TM3 REAL*8 TFMEAN,RMSFIT,YR,TS,TSFIT,RMSP nhalf=nxy/2 RMSmin=1.e20 do n=6,nxy-5 Xknee=x(n) call TREND2(x,F,nxy,Xknee,9999.,2,2, ! input * sl1,sl2,Yknee,RMS,sl,Y0,RMS0) ! output if(RMS.lt.RMSmin) then RMSmin=RMS xmin=Xknee+X0 fpar(1)=sl1 fpar(2)=sl2 fpar(3)=Xknee fpar(4)=Yknee fpar(5)=sl fpar(6)=Y0 rmsp(1)=RMS/nxy rmsp(2)=RMS0/nxy end if end do c write(nw,*) xmin,RMSmin/nxy return end
Followed by flags.f
parameter (lshort=7,slplim=1.,slpx=.5) character*105 title,cc*3 sumch=0. sumchp=0. nsta=0 nstap=0 nyrs=0 nyrsp=0 nshort=0 nsl1=0 nsl2=0 ndsl=0 nok=0 nokx=0 nswitch=0 nsw0=0 open(78,file='fort.78',form='formatted') open( 2,file='fort.2',form='formatted') 10 read(78,'(a)',end=100) title read(title,'(a3,i9,2f9.3,i5,5f9.3,i5,1x,i4,i5,1x,i4)',end=100) * cc,id,sl1,sl2,knee,yk,sl,ylin,rms,rms0,iy1,iy2,iy1e,iy2e nsta=nsta+1 sumch=sumch+sl*(iy2-iy1+1) nyrs=nyrs+(iy2-iy1+1) if(sl.lt.0.) then nstap=nstap+1 nyrsp=nyrsp+(iy2-iy1+1) sumchp=sumchp+sl*(iy2-iy1+1) end if C*** classify : iflag: +1 for short legs etc iflag=0 if(knee.lt.iy1+lshort.or.knee.gt.iy2-lshort) iflag=iflag+1 if(knee.lt.iy1+lshort.or.knee.gt.iy2-lshort) nshort=nshort+1 if(abs(sl1).gt.slplim) iflag=iflag+20 if(abs(sl1).gt.slplim) nsl1=nsl1+1 if(abs(sl2).gt.slplim) iflag=iflag+10 if(abs(sl2).gt.slplim) nsl2=nsl2+1 if(abs(sl2-sl1).gt.slplim) iflag=iflag+100 if(abs(sl2-sl1).gt.slplim) ndsl=ndsl+1 if(abs(sl2-sl1).gt.slpx) iflag=iflag+100 if(iflag.eq.0) nok=nok+1 if(iflag.eq.100) nokx=nokx+1 if(sl1*sl2.lt.0..and.abs(sl1).gt..2.and. * abs(sl2).gt..2) iflag=iflag+1000 if(iflag.ge.1000) nswitch=nswitch+1 if(iflag.eq.1000) nsw0=nsw0+1 write(title(101:105),'(i5)') iflag write(2,'(a)') title go to 10 100 write (*,*) 'all',nsta,-sumch/nsta,-10.*sumch/nyrs write (*,*) 'urb warm',nstap,-sumchp/nstap,-10.*sumchp/nyrsp write(*,*) '# short,sl1,sl2,dsl,ok',nshort,nsl1,nsl2,ndsl,nok,nokx write (*,*) 'switches: all , else ok ',nswitch,nsw0 stop end
The program tr2.f listing:
SUBROUTINE TREND2(xc,A,LEN,Xmid,BAD,MIN1,MIN2, SL1,SL2,Ymid,RMS, * SL,Y0,RMS0) C**** finds a fit using regression analysis by a line C**** with a break in slope at Xmid. Returned are the 2 slopes C**** SL1,SL2 provided we have at least MIN1,MIN2 data. C**** Linear regression data are also computed (for emergencies) REAL*8 xc(*),A(*) REAL*8 sx(2),sxx(2),sxa(2),sa,saa,denom,xnum1,xnum2 INTEGER kount(2) sl1=bad sl2=bad Ymid=bad sa=0. saa=0. do k=1,2 kount(k)=0 sx(k)=0. sxx(k)=0. sxa(k)=0. end do do 100 n=1,len if(a(n).eq.BAD) go to 100 x=xc(n)-Xmid sa=sa+a(n) saa=saa+a(n)**2 k=1 if(x.gt.0.) k=2 kount(k)=kount(k)+1 sx(k)=sx(k)+x sxx(k)=sxx(k)+x**2 sxa(k)=sxa(k)+x*a(n) 100 continue ntot=kount(1)+kount(2) denom=ntot*sxx(1)*sxx(2)-sxx(1)*sx(2)**2-sxx(2)*sx(1)**2 xnum1=sx(1)*(sx(2)*sxa(2)-sxx(2)*sa)+sxa(1)*(ntot*sxx(2)-sx(2)**2) xnum2=sx(2)*(sx(1)*sxa(1)-sxx(1)*sa)+sxa(2)*(ntot*sxx(1)-sx(1)**2) if(kount(1).lt.MIN1.or.kount(2).lt.MIN2) return sl1=xnum1/denom sl2=xnum2/denom Ymid=(sa-sl1*sx(1)-sl2*sx(2))/ntot RMS=ntot*Ymid**2+saa-2*Ymid*(sa-sl1*sx(1)-sl2*sx(2))+ * sl1*sl1*sxx(1)+sl2*sl2*sxx(2)-2*sl1*sxa(1)-2*sl2*sxa(2) C**** linear regression sx(1)=sx(1)+sx(2) sxx(1)=sxx(1)+sxx(2) sxa(1)=sxa(1)+sxa(2) sl=(ntot*sxa(1)-sa*sx(1))/(ntot*sxx(1)-sx(1)**2) Y0=(sa-sl*sx(1))/ntot RMS0=ntot*Y0**2+saa+sl*sl*sxx(1)-2*Y0*(sa-sl*sx(1))-2*sl*sxa(1) return end
=========================================================
The analysis of this program will have to wait a fair while as I’m still finishing another step. It looks like it’s a well done library of utilities so I’ll likely leave it to the end.