GIStemp STEP2_PApars

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.

Advertisements

About E.M.Smith

A technical managerial sort interested in things from Stonehenge to computer science. My present "hot buttons' are the mythology of Climate Change and ancient metrology; but things change...
This entry was posted in GISStemp Technical and Source Code and tagged , , , , , . Bookmark the permalink.