GIStemp STEP345_toSBBXgrid

Overview

We will be looking at the toSBBXgrid.f program. This FORTRAN program is compiled and run by the top level script do_comb_step3.sh and has no wrapper script. Yet another variation on a coding theme… I’ve chosen to leave the program listings as one large block rather than interlieve comments. This is so that experienced programmers can just “see the code” and read it. It also means that “newbies” can just read the commentary and skip the code.

It will only frustrate the folks who sort of remember the programming language and would really like the description interlieved so they would not need to keep hopping back and forth. The code is short enough that that ought not to be much of a problem.

The FORTRAN toSBBXgrid.f

Here is the listing. The explanation will follow below, after the === bar.

C*********************************************************************
C *** PROGRAM READS NCAR FILES
C *** Input files:  31-36 NEWUPD.1-NEWUPD.6
C ***
C *** Output files: 10    subbox.data
C ***               11    box.data
C ***
C *** This program grids the station data
C*********************************************************************
C****
C**** This program interpolates the given station data or their
C**** ANOMALIES with respect to 1951-1980 to a prescribed grid.
C****
C**** Input files:    units 31,32,...,30+INM (see NCARSURF TODISK)
C****
C****  Record 1: M1,INFOI(2),...,INFOI(8),M1L,TITLEI  header record
C****  Record 2: IDATA(M1-->M1L),LT,LN,ID,HT,NAME,M2,M2L  station 1
C****  Record 3: IDATA(M2-->M2L),LT,LN,ID,HT,NAME,M3,M3L  station 2
C****            etc.
C****  Set ITRIM=0 if the input files are not trimmed (i.e. M1L=
C****  M2L=...= INFOI(4) , INFOI(9) and trailing M2L,... missing)
C****
C**** Output files:     units 10 and 11 (regional means)
C****
C****  10: Record 1: INFO(1),...,INFO(8),TITLEO   header record
C****      Record 2: AVG(1 --> MONM),LTS,LTN,LNW,LNE  grid point 1
C****      Record 3: AVG(1 --> MONM),LTS,LTN,LNW,LNE  grid point 2
C****            etc.
C****  11: Record 1: INFOI(1),...,INFOI(8),TITLER     header record
C****      Record 2: AVGR(1-->M0),WTR(1-->M0),NG,LTS,LTN,LNW,LNE box 1
C****      Record 3: AVGR(1-->M0),WTR(1-->M0),NG,LTS,LTN,LNW,LNE box 2
C****            etc.
C****  AVG(1-->MONM) is a full time series, starting at January
C****  of year IYRBEG and ending at December of year IYREND.
C****  IDATA(1) and AVGR(1) refer to Jan of year IYRBG0 which may
C****  be less than IYRBEG, M0=MONM0 is the max. length of an input
C****  time series, and WTR(M) is the area of the part of the region
C****  that contained valid data for month M (in square meters).
C****
C****  NG is the total number of non-missing data on that record.
C****  LTS,LTN is the latitude of the southern,northern edge of the
C****  (sub)box, LNW,LNE the longitude of the western,eastern edge
C****  in hundredths of degrees (all 4-byte integers).
C****
C****  INFO(1),...,INFO(8)  are 4-byte integers,
C****  TITLEO is an 80-byte character string,
C****  AVG,AVGR,WTR are 4-byte reals.
C****      INFO 1 = 1 (indicates that time series are not trimmed)
C****           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 = MONM  (length of each time record)
C****           5 = MONM+4 (size of data record length)
C****           6 = IYRBEG (first year of each time record)
C****           7 = flag for missing data
C****           8 = flag for precipitation trace
C****  INFO(I)=INFOI(I) for I=2,3,7,8
C****  In the output file missing data are flagged by
C****  the real number  XBAD = FLOAT( INFO(7) )
C****
C**** The spatial averaging is done as follows:
C**** Stations within RCRIT km of the grid point P contribute
C**** to the mean at P with weight  1.- d/1200,  (d = distance
C**** between station and grid point 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 (individually
C**** for each month). If that common period is less than 20(NCRIT)
C**** years, the station is disregarded. To decrease that chance,
C**** stations are combined successively in order of the length of
C**** their time record. A final shift then reverses the mean shift
C**** OR (to get anomalies) causes the 1951-1980 mean to become
C**** zero for each month.
C****
C**** Regional means are computed similarly except that the weight
C**** of a grid box with valid data is set to its area.
C**** Separate programs were written to combine regional data in
C**** the same way, but using the weights saved on unit 11.
C****
C?*** Input parameters (# of input files, time period)
      PARAMETER (INM=6,KM0=12,ITRIM=1,KQM=20)
C?*** Out put parameters (output time period, base period)
      PARAMETER (IYBASE=1951,LYBASE=1980)
C?*** Grid dimensions
      PARAMETER (NRM=80,NCM=100)
C?*** Earth radius and limits
      PARAMETER (REARTH=6375.,NCRIT=20,RMAX=999.,RMIN=-999.)
C?*** Work array sizes
      PARAMETER (NSTAM=6000,nStaYR=300000)
      CHARACTER*80 TITLEI,TITLEO
C**** Input arrays
      DIMENSION INFO(8),INFOI(8+ITRIM),ITRL(14+ITRIM)
      integer, allocatable :: idata(:)
      real*4, allocatable :: rdata(:)
C**** Station dependent arrays
      DIMENSION I1S(NSTAM+1),MFS(NSTAM),IORD(NSTAM),ISOFI(NSTAM),
     *   LEN(NSTAM),WTI(NSTAM),ID(NSTAM),ID0(NSTAM),nuseid(nstam)
      REAL*4 CSLONS(NSTAM),CSLATS(NSTAM),SNLONS(NSTAM),SNLATS(NSTAM)
      REAL*4 LONS(NSTAM),LATS(NSTAM) ! lon/lat in degrees
C**** Output-related arrays (dim of output array = MONM < MONM0)
      real*4, allocatable :: avg(:,:),dnew(:),wt(:),avgr(:),wtr(:)
      DIMENSION WTM(12),IORDR(NCM),LENC(NCM)
      COMMON/DIAG/NSTCMB,NSTMNS
C?*** SKIPR and output grid dependent arrays (regions and centers)
      LOGICAL SKIPR/.false./,SKIP(INM,NRM)
      INTEGER INFOR(8)
      REAL*4 PI180,SCALE(KQM)/3*.1,9*0.,.1,7*0./,
     *  CSLONC(NCM,NRM),CSLATC(NCM,NRM),SNLONC(NCM,NRM),SNLATC(NCM,NRM)
      COMMON/GRIDC/PI180,XS(NRM),XN(NRM),XE(NRM),XW(NRM),SKIP,
     *  CSLONC,CSLATC,SNLONC,SNLATC,AREA(NCM,NRM),LATLON(4,NCM+1,NRM)
      COMMON/LIMIT/XBAD,NOVRLP,BIAS(12)
      PI180=DATAN(1.D0)/45.D0
      TOMETR=4.*180.*PI180 * REARTH**2/8000.
      RCRIT=1200.
      IF(IARGC().GT.1) THEN
        CALL GETARG(2,TITLEI)
        READ(TITLEI,*) ICRIT
        RCRIT=ICRIT
      END IF
      RBYRC=REARTH/RCRIT
CSRG  RBYRC0=9004.32/RCRIT
      CSCRIT=COS(RCRIT/REARTH)
      NOVRLP=NCRIT
      IYRBEG=1880
      IF(IARGC().GT.0) THEN
        CALL GETARG(1,TITLEI)
        READ(TITLEI,*) IYRBEG
      END IF
C****
C**** Read and use the header record of an input file
C****
      READ (31) INFOI,TITLEI
      KQ=INFOI(2)
      SCL=1.
      IF(KQ.LE.KQM) SCL=SCALE(KQ)
      IF(SCL.EQ.0.) 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(INFOI(3).EQ.6) KM=12
      IF(INFOI(3).EQ.7) KM=4
      MSIZE = km*nStaYR
      allocate ( rdata(msize),idata(msize) )
      ML=INFOI(4)
      allocate ( avg(ML,ncm),dnew(ML),wt(ML),avgr(ML),wtr(ML) )
      IYREND = ML/KM + INFOI(6) - 1
      MONM0=ML
      NYRSIN=INFOI(4)/KM
      IF(IYRBEG.LT.INFOI(6)) IYRBEG=INFOI(6)
      MONM=(IYREND-IYRBEG+1)*KM
      MFOUT=1+(IYRBEG-INFOI(6))*KM
      NFB=1+IYBASE-INFOI(6)
      NLB=1+LYBASE-INFOI(6)
      MBAD=INFOI(7)
      LAST=INFOI(7)
      XBAD=MBAD
      TRACE=INFOI(8)
C**** Create and write out the header record of the output file
      DO 10 I=1,8
      INFOR(I)=INFOI(I)
   10 INFO(I)=INFOI(I)
      INFO(1)=1
      INFOR(1)=1
      INFO(4)=MONM
      INFO(5)=INFO(4)+7
      INFOR(5)=2*INFOI(4)+1+4
      INFO(6)=IYRBEG
      TITLEO=TITLEI
      TITLEO(20:39)='ANOM (C)  CR     KM '
      WRITE(TITLEO(33:36),'(I4)') ICRIT
      IF(KQ.EQ.2) TITLEO(25:28)='(MM)'
C?    IF(KQ.EQ.3) TITLEO(31:44)='(MB)          '
C?    IF(KQ.EQ.13) TITLEO(32:50)='(PERCENT)          '
C?    IF(IYBASE.GT.0) TITLEO(61:69)='ANOMALIES'
      WRITE (TITLEO(46:57),'(I4,A)') IYRBEG,'-present'
      WRITE(6,'(1X,A80)') TITLEI,TITLEO
      WRITE(6,'('' INFO BOX   '',8I10)') INFOR
      WRITE(6,'('' INFO SUBBOX'',8I10)') INFO
      WRITE (10) INFO,TITLEO
      IF(.NOT.SKIPR) WRITE(11) INFOR,TITLEO(1:39),TITLEI(40:80)
C**** Find grid-dependent arrays (GRIDC-COMMON)
      CALL GRIDEA (NRM,NCM, INM, RBYRC)
C**** Convert areas to units of square-meters
      DO 20 NR=1,NRM
      DO 20 NC=1,NCM
   20 AREA(NC,NR)=AREA(NC,NR)*TOMETR
C****
C**** Loop over NRM large regions
C****
      DO 300 NR=1,NRM
C**** Collect the station data needed for region NR
      IS=0
      I1S(1)=1
      DO 90 IN=1,INM
      IF(SKIP(IN,NR)) GO TO 90
      REWIND 30+IN
      READ (30+IN) INFOI
      MF=INFOI(1)
      IF(ITRIM.GT.0) ML=INFOI(8+ITRIM)
   50 IF(MF.GE.LAST) GO TO 90
      IF(I1S(IS+1)+ML-MF.GT.MSIZE) STOP 'ERROR: NStaYR TOO SMALL'
      CALL SREAD (30+IN,ML+1-MF,IDATA(I1S(IS+1)),ITRL,14+ITRIM)
      XLAT=.1*ITRL(1)
      XLON=.1*ITRL(2)
      MFS(IS+1)=MF
      LENGTH=ML-MF+1
      MF=ITRL(14)
      IF(ITRIM.GT.0) ML=ITRL(14+ITRIM)
      IF(XLAT.LE.XS(NR).OR.XLAT.GE.XN(NR)) GO TO 50
      XLONR=XLON
      IF(XLON.GT.XE(NR)) XLONR=XLON-360.
      IF(XLON.LT.XW(NR)) XLONR=XLON+360.
      IF(XLONR.LE.XW(NR).OR.XLONR.GE.XE(NR)) GO TO 50
      IS=IS+1
      nuseid(is)=0
      IF(IS.GT.NSTAM) STOP 'ERROR: NStaM TOO SMALL'
      I1S(IS+1)=I1S(IS)+LENGTH
      LATS(IS)=XLAT
      LONS(IS)=XLON
      CSLATS(IS)=COS(XLAT*PI180)
      SNLATS(IS)=SIN(XLAT*PI180)
      CSLONS(IS)=COS(XLON*PI180)
      SNLONS(IS)=SIN(XLON*PI180)
      ID(IS)=ITRL(3)
      GO TO 50
   90 CONTINUE
      NSTA=IS
      MSIZE0=I1S(NSTA+1)
      IF(NSTA.EQ.0) THEN
         NGOOD=0
         DO 100 M=1,MONM0
         WTR(M)=0.
  100    AVGR(M)=XBAD
         DO 110 N=1,NCM
  110    CALL SWRITE(10,MONM,AVGR,LATLON(1,NC,NR),XBAD)
         GO TO 290
      END IF
C**** Convert data to real numbers if necessary (KQ<100)
      IF(KQ.LT.100) THEN
         DO 115 N=1,I1S(NSTA+1)-1
  115    RDATA(N)=IDATA(N)
      END IF
C**** Order the NSTA stations according to length of time record
      DO 120 IS=1,NSTA
C 120 LEN(IS)=I1S(IS+1)-I1S(IS)
      LEN(IS)=0
      DO 120 M=I1S(IS),I1S(IS+1)-1
C?LIM IF(RDATA(M) >RMAX or  <RMIN) RDATA(M)=XBAD
      IF(RDATA(M).EQ.XBAD) GO TO 120
      LEN(IS)=LEN(IS)+1
C?*** Change data if necessary (e.g. trace flag for precip)
C?PRC IF(RDATA(M).EQ.TRACE) RDATA(M)=0.
C?LIM IF(RDATA(M).GT.RMAX or <RMIN) RDATA(M)=RMAX or RMIN
  120 CONTINUE
      CALL SORT (IORD,NSTA,LEN)
C****
C**** Find the time series for all NCM centers in region NR
C****
      DO 200 NC=1,NCM
C**** the subbox edges
      Esou=.01*latlon(1,nc,nr)
      Enor=.01*latlon(2,nc,nr)
      Ewst=.01*latlon(3,nc,nr)
      Eest=.01*latlon(4,nc,nr)
C**** Loop over all stations in memory
      IS0=0
      DO 130 N=1,NSTA
      IS=IORD(N)
C**** Find distance between center and station
      CSDBYR=SNLATS(IS)*SNLATC(NC,NR)+CSLATS(IS)*CSLATC(NC,NR)*
     *  (CSLONS(IS)*CSLONC(NC,NR)+SNLONS(IS)*SNLONC(NC,NR))
CSRG  The next line is replaced by 'IF D/RC>1' below in the orig code
      IF(CSDBYR.LE.CSCRIT) THEN
C****   Keep stations in the current subbox even if D>Rcrit
        dbyrc=1
        if(lats(is).gt.Enor) go to 130
        if(lats(is).lt.Esou) go to 130
        if(abs(lats(is)).lt.89..and.lons(is).gt.Eest) go to 130
        if(abs(lats(is)).lt.89..and.lons(is).lt.Ewst) go to 130
      else
        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))
CSRG    IF(CSDBYR.LT.1.) DBYRC=RBYRC0*SQRT(1.-CSDBYR)
CSRG    IF(DBYRC.GT.1.) GO TO 130
      end if
      IS0=IS0+1
      WTI(IS0)=1.-DBYRC
      ISOFI(IS0)=IS
  130 CONTINUE
C**** Set the weight of the stations outside Rcrit to min(insiders)
      wmin=1.
      do is=1,is0
        if(wti(is).gt.0..and.wti(is).lt.wmin) wmin=wti(is)
      end do
      do is=1,is0
        if(wti(is).le.0.) wti(is)=wmin ! *fac (0<fac<1) ???
      end do
C****
C**** Combine the station data
C****
      NSTCMB=0
      NSTMNS=0
      DO 150 M=1,MONM0
      WT(M)=0.
  150 AVG(M,NC)=XBAD
      IF(IS0.EQ.0) THEN
      LENC(NC)=0
         CALL SWRITE(10,MONM,AVG(1,NC),LATLON(1,NC,NR),XBAD)
         WRITE(6,'('' NO STATIONS FOR CENTER '',2I6)') NR,NC
         GO TO 200
      END IF
C**** Start with the station with the longest time record
      IS=ISOFI(1)
      IOFF=MFS(IS)-I1S(IS)
      NSTMNS=LEN(IS)
      NSTCMB=1
      ID0(1)=ID(IS)
      nuseid(is)=nuseid(is)+1
C**** First update of full data and weight arrays
      WMAX=WTI(1)
      DO 155 K=1,KM
      WTM(K)=WTI(1)
  155 BIAS(K)=0.
      DO 160 M=I1S(IS),I1S(IS+1)-1
      AVG(M+IOFF,NC)=RDATA(M)
      IF(RDATA(M).LT.XBAD) WT(M+IOFF)=WTI(1)
  160 CONTINUE
C**** Add in the remaining stations
      DO 190 I=2,IS0
      IS=ISOFI(I)
      IOFF=MFS(IS)-I1S(IS)
C**** Extend the new data into a full series
      DO 170 M=1,MONM0
  170 DNEW(M)=XBAD
      DO 180 M=I1S(IS),I1S(IS+1)-1
  180 DNEW(M+IOFF)=RDATA(M)
      NF1=1+(MFS(IS)-1)/KM
      NL1=1+(I1S(IS+1)-2+IOFF)/KM
C**** Shift new data, then combine them with current mean
C**** Save the average shifts in the array BIAS
      CALL CMBINE (AVG(1,NC),WT, DNEW,NF1,NL1,WTI(I),WTM, KM,
     *  ID(IS),NSM)
      NSTMNS=NSTMNS+NSM
      IF(NSM.EQ.0) GO TO 190
      NSTCMB=NSTCMB+1
      ID0(NSTCMB)=ID(IS)
      nuseid(is)=nuseid(is)+1
      IF(WMAX.LT.WTI(I)) WMAX=WTI(I)
  190 CONTINUE
C**** Set BIAS=time average over the base period if IYBASE > 0
      IF(NFB.GT.0) CALL TAVG(AVG(1,NC),KM,NYRSIN, NFB,NLB, NR,NC,0.)
C**** Subtract BIAS, then scale and write the result to disk
      LENC(NC)=0
      M=0
      DO 195 IY=1,MONM0/KM
      DO 195 K=1,KM
      M=M+1
      IF(AVG(M,NC).EQ.XBAD) GO TO 195
      AVG(M,NC)=SCL*(AVG(M,NC)-BIAS(K))
      LENC(NC)=LENC(NC)+1
  195 CONTINUE
      CALL SWRITE(10,MONM,AVG(MFOUT,NC),LATLON(1,NC,NR),RCRIT*(1.-WMAX))
CW    WRITE(6,'('' CENTER '',3I6,'' STATIONS USED'')') NR,NC,IS0
      IF(NSTCMB.EQ.0.OR.NSTCMB.GT.8) GO TO 200
CW    WRITE(6,'(19I7)') (ID(ISOFI(I)),I=1,IS0)
CW    WRITE(6,'(19I7)') (LEN0(ISOFI(I)),I=1,IS0)
      LATC=.5*(LATLON(1,NC,NR)+LATLON(2,NC,NR))
      LONC=.5*(LATLON(3,NC,NR)+LATLON(4,NC,NR))
      WRITE(6,'('' LAT,LON,STN-MNTHS,STNS,IDS'',3I6,I8,2X,8I10)')
     *  LATC,LONC,NSTMNS,NSTCMB,(ID0(I),I=1,NSTCMB)
CW    WRITE(6,'(19I7)') (NINT(1.E6*WTI(I)),I=1,IS0)
CW    WRITE(6,'(1X,12I4,5X,12I4)')(NINT(10.*AVG(MFOUT-1+M,NC)),M=1,MONM)
  200 CONTINUE
      nuse=0
      do n=1,nstam
      if(nuseid(n).gt.0) then
          write(77,*) 'used station ',id(n),nuseid(n),' times'
          nuse=nuse+1
      end if
      end do
      write(77,*) nuse,' stations used for region', nr
C****
C**** Find the regionally averaged time series
C****
      IF(SKIPR) GO TO 300
      CALL SORT (IORDR,NCM,LENC)
      NC=IORDR(1)
      DO 205 K=1,KM
      WTM(K)=AREA(NC,NR)
  205 BIAS(K)=0.
      DO 210 M=1,MONM0
      WTR(M)=0.
      IF(AVG(M,NC).LT.XBAD) WTR(M)=AREA(NC,NR)
  210 AVGR(M)=AVG(M,NC)
C**** Add in the series of the remaining centers in region NR
      DO 220 N=2,NCM
      NC=IORDR(N)
      IF(LENC(N).EQ.0) GO TO 230
  220 CALL CMBINE (AVGR,WTR,AVG(1,NC),1,MONM0/KM,AREA(NC,NR),WTM,KM,
     *  NC,NSM)
  230 CONTINUE
C**** Set BIAS=time average over the base period if IYBASE > 0
      IF(NFB.GT.0) CALL TAVG (AVGR,KM,NYRSIN, NFB,NLB, NR,0, 0.)
      NGOOD=0
      M=0
      DO 240 IY=1,MONM0/KM
      DO 240 K=1,KM
      M=M+1
      IF(AVGR(M).EQ.XBAD) GO TO 240
      AVGR(M)=AVGR(M)-BIAS(K)
      NGOOD=NGOOD+1
  240 CONTINUE
CW    WRITE(6,'(1X,12I4,5X,12I4)')(NINT(10.*AVGR(M)),M=1,MONM0)
CW    WRITE(6,'(1X,12I4,5X,12I4)')(NINT(WTR(M)/TOMETR),M=1,MONM0)
  290 IF(.NOT.SKIPR) WRITE(11) AVGR,WTR,NGOOD,(LATLON(I,NCM+1,NR),I=1,4)
  300 WRITE(6,'('' REGION'',3I9,'' STATIONS USED'')') NR,MSIZE0,NSTA
      STOP
      END

      SUBROUTINE GRIDEA (NRM1,NCM1, INM1, RBYRC)
C****
C**** This output grid dependent routine sets up the latitude and
C**** longitude limits for the Rcrit km hull of the large regions
C**** and computes the relevant grid quantities.
C****
C**** Current order of boxes:   1-4   north , west->east ...
C****                            .      to
C****                          76-80  south , west->east
C**** Order of subboxes:        1-10  south , west->east ...
C****                            .      to
C****                         91-100  north , west->east
C****
      PARAMETER (NRM=80,ICM=10,JCM=10, INM=6,ISHORT=0, NCM=ICM*JCM)
      LOGICAL SKIP(INM,NRM)
      REAL*4 PI180,SNLATJ(JCM),CSLATJ(JCM),LT100S(JCM),LT100N(JCM),
     *  CSLONC(NCM,NRM),CSLATC(NCM,NRM),SNLONC(NCM,NRM),SNLATC(NCM,NRM)
      COMMON/GRIDC/PI180,XS(NRM),XN(NRM),XE(NRM),XW(NRM),SKIP,
     *  CSLONC,CSLATC,SNLONC,SNLATC,AREA(NCM,NRM),LATLON(4,NCM+1,NRM)
C****
C**** Sergej's equal area grid
C****
C**** Grid constants for latitude zones
CCCCC REAL BANDS(9)/90.,64.2,44.4,23.6,0.,-23.6,-44.4,-64.2,-90./
C**** We need only the SINEs of the (northern) band edges
      REAL SNNEDG(9)/1.,.9,.7,.4,0.,-.4,-.7,-.9,-1./
C     Input data sets - 1: 90N-60N  2: 60N-30N ... 6: 60S-90S 7:short
      INTEGER INFRST(8)/1,1,2,2,3,4,5,5/,NUMJ(8)/4,8,12,16,16,12,8,4/
      INTEGER INLAST(8)/2,2,3,4,5,5,6,6/
C**** Check grid dimensions
      IF(NRM1.NE.NRM.OR.ICM*JCM.NE.NCM1.OR.INM1.NE.INM) THEN
         WRITE(6,'(6I9)') NRM1,NRM, ICM*JCM,NCM1, INM1,INM
         STOP 'ERROR: GRID SIZES INCONSISTENT'
      END IF
C**** Don't skip any short files if ISHORT=1
      DO 10 NR=1,NRM
      SKIP(INM,NR)=.FALSE.
      DO 10 IN=1,INM-ISHORT
   10 SKIP(IN,NR)=.TRUE.
C**** Find DDLAT to extend each box R km in north and south direction
      DDLAT=1./(PI180*RBYRC)
C**** Loop over all BOXES (large regions)
      NR=0
      DO 50 J=1,8
C**** Find the sin of the latitudes of the centers in band J
      SNN=SNNEDG(J)
      SNS=SNNEDG(J+1)
      DSLATJ=(SNN-SNS)/JCM
      DO 20 JC=1,JCM
      LT100S(JC)=NINT(100./PI180*ASIN(SNS+(JC-1)*DSLATJ))
      LT100N(JC)=NINT(100./PI180*ASIN(SNS+JC*DSLATJ))
      SNLATJ(JC)=SNS+(JC-.5)*DSLATJ
   20 CSLATJ(JC)=SQRT(1.-SNLATJ(JC)**2)
      DO 50 I=1,NUMJ(J)
      NR=NR+1
C**** Complete definition of SKIP-array
      DO 30 IN=INFRST(J),INLAST(J)
   30 SKIP(IN,NR)=.FALSE.
      DLON=360./NUMJ(J)
      XEAST=-180.+I*DLON
      XWEST=XEAST-DLON
C**** Extend each box by half a box in each direction (>Rcrit km ?)
      XW(NR)=XWEST-.5*DLON
      IF(J.EQ.1.OR.J.EQ.8) XW(NR)=-180.
      XE(NR)=XEAST+.5*DLON
      IF(J.EQ.1.OR.J.EQ.8) XE(NR)=180.
      XN(NR)=ASIN(SNN)/PI180+DDLAT
      XS(NR)=ASIN(SNS)/PI180-DDLAT
C**** Find sine and cosine of the center latitudes and longitudes
      DLONC=DLON/ICM
      DO 40 IC=1,ICM
      LN100W=NINT(100.*(XWEST+(IC-1)*DLONC))
      LN100E=NINT(100.*(XWEST+IC*DLONC))
      RLONI=PI180*(XWEST+(IC-.5)*DLONC)
      SNLONI=SIN(RLONI)
      CSLONI=COS(RLONI)
      DO 40 JC=1,JCM
CSRG  ICJC=JC+JCM*(IC-1)
      ICJC=IC+ICM*(JC-1)
      LATLON(1,ICJC,NR)=LT100S(JC)
      LATLON(2,ICJC,NR)=LT100N(JC)
      LATLON(3,ICJC,NR)=LN100W
      LATLON(4,ICJC,NR)=LN100E
      AREA(ICJC,NR)=1.
      SNLONC(ICJC,NR)=SNLONI
      CSLONC(ICJC,NR)=CSLONI
      SNLATC(ICJC,NR)=SNLATJ(JC)
   40 CSLATC(ICJC,NR)=CSLATJ(JC)
      LATLON(1,NCM+1,NR)=LT100S(1)
      LATLON(2,NCM+1,NR)=LT100N(JCM)
      LATLON(3,NCM+1,NR)=NINT(100.*XWEST)
      LATLON(4,NCM+1,NR)=NINT(100.*(XWEST+DLON))
   50 CONTINUE
      RETURN
      END

      SUBROUTINE CMBINE (AVG,WT, DNEW,NF1,NL1,WT1,WTM, KM, ID,NSM)
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,NOVRLP,BIAS(12)
      DIMENSION AVG(KM,*),DNEW(KM,*),WT(KM,*),WTM(*),MISSNG(12)
C**** Loop over months or seasons if appropriate
      NSM=0
      MISSED=KM
      DO 50 K=1,KM
      MISSNG(K)=1
C**** Find means over common domain to compute bias
      SUMN=0
      SUM=0
      NCOM=0
      DO 10 N=NF1,NL1
      IF(AVG(K,N).GE.XBAD.OR.DNEW(K,N).GE.XBAD) GO TO 10
      NCOM=NCOM+1
      SUM=SUM+AVG(K,N)
      SUMN=SUMN+DNEW(K,N)
   10 CONTINUE
      IF(NCOM.LT.NOVRLP) GO TO 50
      BIASK=(SUM-SUMN)/FLOAT(NCOM)
C**** Find mean bias
      WTMNEW=WTM(K)+WT1
      BIAS(K)=(WTM(K)*BIAS(K)+WT1*BIASK)/WTMNEW
      WTM(K)=WTMNEW
C**** Update period of valid data, averages and weights
      DO 20 N=NF1,NL1
      IF(DNEW(K,N).GE.XBAD) GO TO 20
      WTNEW=WT(K,N)+WT1
      AVG(K,N)=(WT(K,N)*AVG(K,N)+WT1*(DNEW(K,N)+BIASK))/WTNEW
      WT(K,N)=WTNEW
      NSM=NSM+1
   20 CONTINUE
      MISSED=MISSED-1
      MISSNG(K)=0
   50 CONTINUE
CLOG? IF(MISSED.GT.0) WRITE(6,90) ID,WT1,MISSNG
   90 FORMAT(' UNUSED DATA - ID/SUBBOX,WT',I8,F5.2,12I2)
      RETURN
      END

      SUBROUTINE TAVG (DATA,KM,NYRS, NFB,NLB, NR,NC, DEFLT)
C****
C**** TAVG computes the time averages (separately for each calendar
C**** month if KM=12) over the base period (year NFB to NLB) and
C**** saves them in BIAS. In case of no data, the average is set to
C**** DEFLT if NR=0 or computed over the whole period if NR>0.
C****
      COMMON/LIMIT/XBAD,NOVRLP,BIAS(12)
      DIMENSION DATA(KM,*),LEN(12)
      MISSED=KM
      DO 50 K=1,KM
      BIAS(K)=DEFLT
      SUM=0.
      M=0
      DO 10 N=NFB,NLB
      IF(DATA(K,N).GE.XBAD) GO TO 10
      M=M+1
      SUM=SUM+DATA(K,N)
   10 CONTINUE
      LEN(K)=M
      IF(M.EQ.0) GO TO 50
      BIAS(K)=SUM/FLOAT(M)
      MISSED=MISSED-1
   50 CONTINUE
      IF(NR*MISSED.EQ.0) RETURN
C**** If base period is data free, use bias with respect to whole series
      DO 100 K=1,KM
      IF(LEN(K).GT.0) GO TO 100
      WRITE(6,'(''0NO DATA IN BASE PERIOD - MONTH,NR,NC'',3I9)') K,NR,NC
      SUM=0.
      M=0
      DO 60 N=1,NYRS
      IF(DATA(K,N).GE.XBAD) GO TO 60
      M=M+1
      SUM=SUM+DATA(K,N)
   60 CONTINUE
      IF(M.EQ.0) GO TO 100
      BIAS(K)=SUM/FLOAT(M)
  100 CONTINUE
      RETURN
      END

      SUBROUTINE SORT (INDEX,NDIM,LNGTH)
C**** Sorts INDEX and LNGTH such that LNGTH becomes decreasing
      DIMENSION INDEX(NDIM),LNGTH(NDIM)
      DO 10 N=1,NDIM
   10 INDEX(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(LNGTH(NN).GT.LNGTH(NLMAX)) NLMAX=NN
   20 CONTINUE
C**** Switch positions N and NLMAX
      LMAX=LNGTH(NLMAX)
      LNGTH(NLMAX)=LNGTH(N)
      LNGTH(N)=LMAX
      IMAX=INDEX(NLMAX)
      INDEX(NLMAX)=INDEX(N)
   30 INDEX(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

      SUBROUTINE SWRITE (IOUT,NDIM,ARRAY,LATLON,DMIN)
      COMMON/DIAG/NSTNS,NSTMNS
      DIMENSION ARRAY(NDIM),LATLON(4)
      WRITE(IOUT) ARRAY,LATLON,NSTNS,NSTMNS,DMIN
      RETURN
      END

      SUBROUTINE LNFIT(DATA,NS,BAD, TREND, IGFRST,IGLAST,NGOOD)
      REAL*8 F1,F2,AX,A1,A2,A22
C**   Determines linear fit using regression analysis
C**   DATA(I) - time series to be fitted
C**   NS - length of the time series
C**   TREND = slope * total number of years
C**   IGFRST,IGLAST = first and last non-missing position
C**   NGOOD = number of good data (trend BAD if NGOOD<2)
      DIMENSION DATA(NS)
      IGFRST=0
      NGOOD=0
      NGSUM=0
      NGSUM2=0
      F1=0.
      F2=0.
      DO 10 I=1,NS
      IF(DATA(I).EQ.BAD) GO TO 10
      IF(IGFRST.EQ.0) IGFRST=I
      IGLAST=I
      NGOOD=NGOOD+1
      NGSUM=NGSUM+I
      NGSUM2=NGSUM2+I**2
      F1=F1+DATA(I)
      F2=F2+DATA(I)*I
   10 CONTINUE
      TREND=BAD
      IF(NGOOD.LT.2) RETURN
      A1=NGOOD
      A2=NGSUM
      A22=NGSUM2
      AX=A1*A22-A2**2
      TREND=NS*(F2*A1-F1*A2)/AX
C     A=(NGSUM2*F1-F2*NGSUM)/AX
      RETURN
      END

=========================================================

Almost painfully long, given the old school ALL CAPS style. The style of the “do loops” is more f77 (or even FORTRAN IV). This is probably one of the older parts of the “system”.

The analysis of this program will have to wait a bit as I’m still finishing an earlier step. It looks like it mostly just smears data over time and space to fill in missing chunks. I suspect that is fraught with “issues”.