GIStemp STEP345_zonav

Overview

For now, this will be just a listing of the programs. Exposition will have to wait until I’ve gotten through some of the earilier bits.

We will be looking at the script zonav and the two FORTRAN programs zonav.f and annzon.f in this section. 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.

A word or two about execution order. In most cases we find that there is a wrapper script around a FORTRAN program that does the compile, execution, deletion. This case is different. In this case the script zonav does the compile and the execution of two programs.

First, the script: zonav

Begin script listing

#!/bin/ksh

fortran_compile=$FC
if [[ $FC = '' ]]
then echo "set an environment variable FC to the fortran_compile_command like f90"
     echo "or do all compilation first and comment the compilation lines"
     exit
fi

label=$1 ; rad=1200

if [[ ! -s BX.Ts.${label}.$rad ]]
then echo "BX.Ts.${label}.$rad not found"
     exit
fi

${fortran_compile} zonav.f -o zonav.exe

ln BX.Ts.${label}.$rad fort.11
zonav.exe > zonav.Ts.${label}.log
cp fort.10 work_files/ZON.Ts.${label}.1200.step1
rm fort.11     zonav.exe

${fortran_compile} annzon.f -o annzon.exe
annzon.exe  > annzon.Ts.${label}.log   ; rm annzon.exe
mv fort.96 ZonAnn.Ts.${label}.txt
mv fort.97 GLB.Ts.${label}.txt
mv fort.98 NH.Ts.${label}.txt
mv fort.99 SH.Ts.${label}.txt
mv fort.12 work_files/ZON.Ts.${label}.1200
mv fort.11 work_files/ANNZON.Ts.${label}.1200
rm fort.10
mkdir results 2> /dev/null
mv *txt results/.

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

End script listing

So what does it do?

Notice that it compiles and runs both zonav.f and annzon.f so this one is a two-fer. It opens with the usual compiler environment variable, error if missing and exit, then follows with a check for the input file and linking its input file name to fort.11 while after the run, fort.10 gets copied to work_files (so I presume that’s the output file).

After the usual delete the binary, we do it all over again for annzon.f and finish up with some housekeeping.

That’s it. The whole thing.

So after we run the script what then?

The FORTRAN zonav.f

Here is the listing. The explanation will follow below at some future date, after the === bar.

    C*********************************************************************
C *** PROGRAM READS BXdata
C *** Input files:  11    box.data (BX1977.T1200)
C ***
C *** Output files: 10    zonal.means (ZON1977.T1200)
C ***
C*********************************************************************
C****
C**** This program combines the given gridded data (anomalies)
C**** to produce AVERAGES over various LATITUDE BELTS.
C****
C**** Input file:    unit 11   (=output of job NCARSURF GRIDDING)
C****
C****  11: Record 1: INFOI(1),...,INFOI(8),TITLEI   header record
C****      Record 2: AR(1-->MONM0),WTR(1-->MONM0),NG  grid point 1
C****      Record 3: AR(1-->MONM0),WTR(1-->MONM0),NG  grid point 2
C****            etc.
C****
C**** Output files:     unit 10
C****
C****  10: Record 1: INFO(1),...,INFO(8),TITLEO,TXT    header record
C****      Record 2: DATA(1-->MONM),WT(1-->MONM),TITLE1     belt 1
C****      Record 3: DATA(1-->MONM),WT(1-->MONM),TITLE2     belt 2
C****            etc.
C****  DATA(1-->MONM) is a full time series, starting at January
C****  of year IYRBEG and ending at December of year IYREND.
C****  WT is proportional to the area containing valid data.
C****  AR(1) refers to Jan of year IYRBG0 which may
C****  be less than IYRBEG, MONM0 is the length of an input time
C****  series, and WTR(M) is the area of the part of the region
C****  that contained valid data for month M.
C****  NG is the total number of non-missing data on that record.
C****  TITLE1,... describe the latitude belts.
C****
C****  INFO(1),...,INFO(8)  are 4-byte integers,
C****  All TITLEs and TXT are 80-byte character strings,
C****  all other entries in records 2,3,... are 4-byte reals.
C****      INFO 1 and 5 are irrelevant
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****           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) except perhaps for I=4 and I=6.
C****  In the output file missing data are flagged by
C****  the real number  XBAD = FLOAT( INFO(7) )
C****
C**** JBM zonal means are computed first, combining successively
C**** the appropriate regional data (AR with weight WTR). To remove
C**** the regional bias, the data of a new region are shifted
C**** so that the mean over the common period remains unchanged
C**** after its addition. If that common period is less than
C**** 20(NCRIT) years, the region is disregarded. To avoid that
C**** case as much as possible, regions are combined in order of
C**** the length of their time record. A final shift causes the
C**** 1951-1980 mean to become zero (for each month).
C****
C**** All other means (incl. hemispheric and global means) are
C**** computed from these zonal means using the same technique.
C**** NOTE: the weight of a zone may be smaller than its area
C**** since data-less parts are disregarded; this also causes the
C**** global mean to be different from the mean of the hemispheric
C**** means.
C****
C?*** Input parameters (# of input files, time period)
C?*** Out put parameters (output time period, base period)
      PARAMETER (IYBASE=1951,LYBASE=1980)
C?*** Grid dimensions and limit parameters
      PARAMETER (NRM=80,JBM=8,IBMM=16,NZS=3 ,NCRIT=20)
      CHARACTER*80 TITLEI,TITLEO,TITLEZ(JBM+NZS+3)/
     * JBM*'see subr GRIDx',NZS*'see subr GRIDx',
     * 'NORTHERN HEMISPHERE','SOUTHERN HEMISPHERE','GLOBAL MEANS'/
C?*** Special title for Jim's plots
      CHARACTER*40 TXT(2)/
     *  ' zones:  90->64.2->44.4->23.6->0->-23.6-',
     *  '>-44.4->-64.2->-90                      '/
C**** Input arrays
      DIMENSION INFOI(8),LENR(IBMM)
      real*4, allocatable :: ar(:,:),wtr(:,:)
C**** Output-related arrays (dim of output array = MONM < MONM0)
      DIMENSION INFO(8),IORD(JBM+IBMM),LENZ(JBM)
      real*4, allocatable :: AVG(:,:),WT(:,:),AVGG(:),WTG(:)
      COMMON/GRIDC/IBM(JBM),KZONE(JBM,NZS+3)
      COMMON/LIMIT/XBAD,NOVRLP,BIAS(12)
      NOVRLP=NCRIT
C****
C**** Read and use the header record of an input file
C****
      READ (11) INFOI,TITLEI
      KQ=INFOI(2)
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
      ML=INFOI(4)
      NYRSIN=INFOI(4)/KM
      IYRBEG=INFOI(6)
      MONM=ML ; IYREND = MONM/KM + IYRBEG - 1
      allocate ( ar(ml,IBMM),wtr(ml,IBMM), AVG(ml,JBM),WT(ml,JBM) )
      allocate ( AVGG(ml),WTG(ml) )
      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
   10 INFO(I)=INFOI(I)
      INFO(4)=MONM
      INFO(6)=IYRBEG
      TITLEO=TITLEI
      WRITE(6,'(1X,A80)') TITLEI,TITLEO
      WRITE(6,'(8I10)') INFOI
      WRITE(6,'(8I10)') INFO
      WRITE (10) INFO,TITLEO,TXT
C**** Find grid-dependent arrays (GRIDC-COMMON)
      CALL GRIDEA (NRM,TITLEZ,IBMM,NZS)
C****
C**** Find the means over the JBM basic latitude belts
C****
      DO 300 JB=1,JBM
C**** Collect the data for the IBM(JB) regions in the belt JB
      DO 30 IB=1,IBM(JB)
      CALL SREAD(11,ML   ,AR(1,IB),WTR(1,IB),LENR(IB))
   30 CONTINUE
C**** Check whether there are any data
      LNTOT=0
      DO 40 N=1,IBM(JB)
      LNTOT=LNTOT+LENR(N)
   40 CONTINUE
      IF(LNTOT.EQ.0) THEN
         WRITE(6,*) ' **** NO DATA FOR ZONE ',JB,TITLEZ(JB)
         GO TO 300
      END IF
C**** Order the IBM regions according to length of time record
      CALL SORT (IORD,IBM(JB),LENR)
      NR=IORD(1)
C**** Start with longest time series
      DO 210 M=1,MONM
      WT(M,JB)=WTR(M,NR)
  210 AVG(M,JB)=AR(M,NR)
C**** Add in the series of the remaining regions in belt JB
      DO 220 N=2,IBM(JB)
      NR=IORD(N)
      IF(LENR(N).EQ.0) GO TO 230
      CALL CMBINE (AVG(1,JB),WT(1,JB), AR(1,NR), 1,NYRSIN,
     *             WTR(1,NR),WTM, KM, JB)
  220 CONTINUE
  230 CONTINUE
C**** Set BIAS=time average over the base period if IYBASE > 0
      CALL TAVG(AVG(1,JB),KM,NYRSIN, NFB,NLB, JB, 0.)
      LENZ(JB)=0
      M=0
      DO 240 IY=1,NYRSIN
      DO 240 K=1,KM
      M=M+1
      IF(AVG(M,JB).EQ.XBAD) GO TO 240
      AVG(M,JB)=AVG(M,JB)-BIAS(K)
      LENZ(JB)=LENZ(JB)+1
  240 CONTINUE
      write(*,*) 'zonal mean',JB,TITLEZ(JB)
      WRITE(6,'(1X,12I4,5X,12I4)')(NINT(10.*AVG(M,JB)),M=1,MONM)
      write(*,*) 'weights'
      WRITE(6,'(1X,12I4,5X,12I4)')(NINT(WT(M,JB)/60000.+.499),M=1,MONM)
      CALL SWRITE(10,MONM,AVG(MFOUT,JB),WT(MFOUT,JB),TITLEZ(JB))
      WRITE(6,*) AVG(MFOUT,JB),WT(MFOUT,JB),TITLEZ(JB),' SAVED ON DISK'
  300 CONTINUE
C****
C**** Find the means for the remaining NZS+3 bands
C****
      CALL SORT (IORD,JBM,LENZ)
      DO 400 JZ=1,NZS+3
C**** Start with the longest time series
      IF(LENZ(1).EQ.0) THEN
         WRITE(6,*) ' **** NO DATA FOR ZONE ',JBM+JZ,TITLEZ(JBM+JZ)
         GO TO 400
      END IF
      DO 310 J1=1,JBM
      IF(KZONE(IORD(J1),JZ).GT.0) GO TO 320
  310 CONTINUE
  320 JB=IORD(J1)
      DO 330 M=1,MONM
      WTG(M)=WT(M,JB)
  330 AVGG(M)=AVG(M,JB)
C**** Add in the remaining latitude belts JB with KZONE(JB,JZ)=1
      DO 340 J=J1+1,JBM
      IF(KZONE(IORD(J),JZ).EQ.0) GO TO 340
      JB=IORD(J)
      IF(LENZ(J).EQ.0) GO TO 360
      CALL CMBINE (AVGG,WTG,AVG(1,JB), 1,NYRSIN, WT(1,JB),WTM,KM,JBM+JZ)
  340 CONTINUE
C**** Set BIAS=time average over the base period if IYBASE > 0
      IF(NFB.GT.0) CALL TAVG(AVGG,KM,NYRSIN, NFB,NLB, JBM+JZ, 0.)
      M=0
      DO 350 IY=1,NYRSIN
      DO 350 K=1,KM
      M=M+1
      IF(AVGG(M).NE.XBAD) AVGG(M)=AVGG(M)-BIAS(K)
  350 CONTINUE
  360 CONTINUE
      write(*,*) 'zone',JBM+JZ
      WRITE(6,'(1X,12I4,5X,12I4)')(NINT(10.*AVGG(M)),M=1,MONM)
      write(*,*) 'weights',JBM+JZ
      WRITE(6,'(1X,12I4,5X,12I4)')(NINT(WTG(M)/60000.+.49999),M=1,MONM)
      CALL SWRITE(10,MONM,AVGG(MFOUT),WTG(MFOUT),TITLEZ(JBM+JZ))
      WRITE(6,*) AVGG(MFOUT),WTG(MFOUT),TITLEZ(JBM+JZ),' SAVED ON DISK'
  400 CONTINUE
      STOP
      END

      SUBROUTINE GRIDEA (NRM0,TITLEZ, IBMM0,NZS0)
C****
C**** This output grid dependent routine sets the parameters
C**** IBM(J) (= number of regions in latitude belt J)  and
C**** KZONE(J,N) (=1 if belt J belongs to domain N) and the TITLES
C****
C**** Current order of boxes:   1-4   north , west->east ...
C****                            .      to
C****                          77-80  south , west->east
C**** Order of subboxes:        1-10  south , west->east ...
C****                            .      to
C****                         91-100  north , west->east
C****   Sergei's ordering       1-10  west , south->north
C****                            .      to
C****                         91-100  east , south->north
C****
      PARAMETER (NRM=80,ICM=10,JCM=10, NZS=3, NCM=ICM*JCM)
      CHARACTER*80 TITLEZ(*),TITLES(1+NZS)/
     *    '  LATITUDE BELT FROM  XX.X ? TO  XX.X ?',
     *    '  NORTHERN LATITUDES: 23.6 N TO  90.0 N',
     *    '       LOW LATITUDES: 23.6 S TO  23.6 N',
     *    '  SOUTHERN LATITUDES: 90.0 S TO  23.6 S'/
      COMMON/GRIDC/IBM(8),KZONE(8,NZS+3)
C****
C**** Sergej's equal area grid
C****
C**** Grid constants for latitude zones
      REAL BANDS(9)/90.,64.2,44.4,23.6,0.,-23.6,-44.4,-64.2,-90./
      INTEGER NUMJ(8)/4,8,12,16,16,12,8,4/
C**** Check grid dimensions
      IF(NRM0.NE.NRM.OR.NZS.NE.NZS0) THEN
         WRITE(6,'(4I9)') NRM0,NRM, NZS0,NZS
         STOP 'ERROR: INCONSISTENCY FOR NRM NCM NZS'
      END IF
C**** Loop over all ZONES
      DO 50 J=1,8
      IBM(J)=NUMJ(J)
      IF(IBM(J).LE.IBMM0) GO TO 20
      WRITE(6,*) 'IBMM SHOULD BE AT LEAST EQUAL TO',IBM(J),' J=',J
      STOP 'IBMM TOO SMALL'
C**** Define the large bands as combination of the basic zones
   20 KZONE(J,1)=0
      IF(J.LE.3) KZONE(J,1)=1
      KZONE(J,2)=0
      IF(J.EQ.4.OR.J.EQ.5) KZONE(J,2)=1
      KZONE(J,3)=0
      IF(J.GE.6) KZONE(J,3)=1
      KZONE(J,NZS+1)=0
      IF(J.LE.4) KZONE(J,NZS+1)=1
      KZONE(J,NZS+2)=1-KZONE(J,NZS+1)
      KZONE(J,NZS+3)=1
C**** Set the titles TITLEZ for the basic latitude bands
      TITLEZ(J)=TITLES(1)
      WRITE(TITLEZ(J)(23:26),'(F4.1)') ABS(BANDS(J+1))
      WRITE(TITLEZ(J)(28:28),'(A1)') 'N'
      IF(BANDS(J+1).LT.0.) WRITE(TITLEZ(J)(28:28),'(A1)') 'S'
      IF(BANDS(J+1).EQ.0.) WRITE(TITLEZ(J)(22:28),'(A7)') 'EQUATOR'
      WRITE(TITLEZ(J)(34:37),'(F4.1)') ABS(BANDS(J))
      WRITE(TITLEZ(J)(39:39),'(A1)') 'N'
      IF(BANDS(J).LT.0.) WRITE(TITLEZ(J)(39:39),'(A1)') 'S'
      IF(BANDS(J).EQ.0.) WRITE(TITLEZ(J)(33:39),'(A7)') 'EQUATOR'
   50 CONTINUE
C**** Set the titles for the special zones
      DO 60 J=1,NZS
   60 TITLEZ(8+J)=TITLES(1+J)
      RETURN
      END

      SUBROUTINE CMBINE (AVG,WT, DNEW,NF1,NL1,WT1,WTM, KM, ID)
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,*),WT1(KM,*),MISSNG(12)
C**** Loop over months or seasons if appropriate
      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
C?    WTMNEW=WTM+WT1
C?    BIAS(K)=(WTM*BIAS(K)+WT1*BIASK)/WTMNEW
C?    WTM=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(K,N)
      AVG(K,N)=(WT(K,N)*AVG(K,N)+WT1(K,N)*(DNEW(K,N)+BIASK))/WTNEW
      WT(K,N)=WTNEW
   20 CONTINUE
      MISSED=MISSED-1
      MISSNG(K)=0
   50 CONTINUE
      IF(MISSED.GT.0) WRITE(6,90) ID,MISSNG
   90 FORMAT(' UNUSED DATA - ID/SUBBOX',I8,12I2)
      RETURN
      END

      SUBROUTINE TAVG (AVG,KM,NYRS, NFB,NLB, JB, DEFLT)
C****
C**** Data are shifted by KM constants such that the KM averages
C**** over the base period (year NFB to NLB) are zero
C****
      COMMON/LIMIT/XBAD,NOVRLP,BIAS(12)
      DIMENSION AVG(KM,*),LEN(12)
      MISSED=KM
      DO 50 K=1,KM
      BIAS(K)=DEFLT
      SUM=0.
      M=0
      DO 10 N=NFB,NLB
      IF(AVG(K,N).GE.XBAD) GO TO 10
      M=M+1
      SUM=SUM+AVG(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(JB*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,JB'',3I9)') K,JB
      SUM=0.
      M=0
      DO 60 N=1,NYRS
      IF(AVG(K,N).GE.XBAD) GO TO 60
      M=M+1
      SUM=SUM+AVG(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,DATA,WT,NG)
C**** Speed read routine for input records
      DIMENSION DATA(LEN),WT(LEN)
      READ (IN) DATA,WT,NG
      RETURN
      END

      SUBROUTINE SWRITE (IOUT,NDIM,AR,WT,TITLE)
      DIMENSION AR(NDIM),WT(NDIM),TITLE(20)
      WRITE(IOUT) AR,WT,TITLE
      RETURN
      END

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

A bit longer than the others so far, and the style is different. Different intern? Given the “old school” style of ALL CAPS and the f77 style of “do loop” this is likely one of the oldest parts of the “system”.

The exposition of details will have to wait until I’m done with the earlier steps.

What about annzon?

The FORTRAN annzon.f

Here is the listing. The explanation will follow at some future time below, after the === bar.

C*********************************************************************
C *** program reads ZONAL monthly means and recomputes REGIONAL means
C *** as well as annual means.
C *** Input file:  10    zonal.means (ZON1977.T1200)
C ***
C *** Output files: 11    annual.zonal.means (ANNZON1977.T1200)
C ***               12    zonal.means (ZON1977.T1200) changed
C*********************************************************************
C****
C****
C**** Displays data and their annual means
C****
C**** Input file:    unit 10   (=output of job NCARSURF ZONAVG)
C****
C****  10: Record 1: INFO(1),...,INFO(8),TITLE       header record
C****      Record 2: DATA(1-->MONM),WT(1-->MONM),TITLE1     belt 1
C****      Record 3: DATA(1-->MONM),WT(1-->MONM),TITLE2     belt 2
C****            etc.
C****  DATA(1-->MONM) is a full time series, starting at January
C****  of year IYRBEG and ending at December of year IYREND.
C****  WT is proportional to the area containing valid data.
C****  TITLE1,... describe the latitude belts.
C****
C****  INFO(1),...,INFO(8)  are 4-byte integers,
C****  All TITLEs are 80-byte character strings,
C****  all other entries in records 2,3,... are 4-byte reals.
C****      INFO 1 and 5 are irrelevant
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****           6 = IYRBEG (first year of each time record)
C****           7 = flag for missing data
C****           8 = flag for precipitation trace
C****  Flag for missing data:  XBAD = FLOAT( INFO(7) )
C****
C?*** Input parameters (time period - use LINFO to find it)
      PARAMETER (IY1TAB=1880)
C?*** Grid dimensions and limit parameters
C?*** Alternate global (IAVGGH=1) and hemispheric means (IAVGGH=2)
      PARAMETER (IAVGGH=2, JBM=8,NZS=3, JZM=JBM+NZS+3, JZP=14, MONMIN=6)
      CHARACTER*80 TITLE,TITL2,TITLEZ(JZM),BLANK/' '/,TITLEO/
     * 'ANNUALLY AVERAGED (7 OR MORE MONTHS) TEMPERATURE ANOMALIES (C)'/
      CHARACTER*10 TIT(3)/'    GLOBAL','N.HEMISPH.','S.HEMISPH.'/
C**** Input arrays
      DIMENSION INFO(8),INFOO(8)
      real*4, allocatable :: DATA(:,:,:),WT(:,:,:) ! (KM,IYRSX,JZM)
      INTEGER IORD(JZP)/14,12,13, 9,10,11, 1,2,3,4,5,6,7,8/,IOUT(18)
C**** Output-related arrays (dim of output array = MONM < MONM0)
      INTEGER JZG(4,2)/9,10,10,11, 9,4,5,11/
      REAL, dimension(4) :: WTSP=(/3.,2.,2.,3./)
      REAL, allocatable :: ANN(:,:),ANNW(:,:) !  (IYRSX,JZM)
C****
C**** Read and use the header record of an input file
C****
      READ (10) INFO,TITLE,TITL2
      KQ=INFO(2)
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
      IYRBEG=INFO(6)
      MONM=INFO(4)
      IYRS=MONM/KM
      allocate ( data(km,iyrs,jzm),wt(km,iyrs,jzm) )
      allocate ( ann (   iyrs,jzm),annw( iyrs,jzm) )
      IYREND=INFO(4)/KM+IYRBEG-1
      MBAD=INFO(7)
      XBAD=MBAD
C**** Create and write out the header record of the output file
      WRITE(6,'(1X,A80)') TITLE
      WRITE(96,*) 'Annual Temperature Anomalies (.01 C) - ',TITLE(29:80)
      WRITE(97,'(A80)') TITLE
      WRITE(98,'(A80)') TITLE
      WRITE(99,'(A80)') TITLE
      WRITE(6,'(8I10)') INFO
      DO 10 I=1,8
   10 INFOO(I)=INFO(I)
      INFOO(3)=5
      INFOO(4)=INFO(4)/12
      WRITE(TITLEO(20:20),'(I1)') MONMIN
      WRITE(6,'(''1'',A80)') TITLEO
      WRITE(6,'(8I10)') INFOO
C****
C**** Collect the JZM zonal means (basic+special+hemi+global)
C****
      DO 30 JZ=1,JZM
      CALL SREAD(10,MONM,DATA(1,1,JZ),WT(1,1,JZ),TITLEZ(JZ))
   30 CONTINUE
C****
C**** Find the annual means
C****
      DO 100 JZ=1,JZM
      DO 100 IY=1,IYRS
      ANN(IY,JZ)=XBAD
      ANNW(IY,JZ)=0.
      ANNIY=0.
      ANNWIY=0.
      MON=0
      DO 50 M=1,KM
      IF(DATA(M,IY,JZ).EQ.XBAD) GO TO 50
      MON=MON+1
      ANNIY=ANNIY+DATA(M,IY,JZ)
      ANNWIY=ANNWIY+WT(M,IY,JZ)
   50 CONTINUE
      IF(MON.GE.MONMIN) ANN(IY,JZ)=ANNIY/MON
      ANNW(IY,JZ)=ANNWIY/12.
  100 CONTINUE
C****
C**** Alternate global mean (from North.lats,Equ.reg,South.lats)
C****
      IF(IAVGGH.EQ.0) GO TO 180
      IGLB=IAVGGH
      DO 120 IY=1,IYRS
      GLOB=0.
      ANN(IY,JZM)=XBAD
      DO 110 J=1,4
      IF(ANN(IY,JZG(J,IGLB)).EQ.XBAD) GO TO 120
      GLOB=GLOB+ANN(IY,JZG(J,IGLB))*WTSP(J)
  110 CONTINUE
      ANN(IY,JZM)=.1*GLOB
  120 CONTINUE
      DO 130 IY=1,IYRS
      DO 130 M=1,12
      DATA(M,IY,JZM)=XBAD
      GLOB=0.
      DO 125 J=1,4
      IF(DATA(M,IY,JZG(J,IGLB)).EQ.XBAD) GO TO 130
      GLOB=GLOB+DATA(M,IY,JZG(J,IGLB))*WTSP(J)
  125 CONTINUE
      DATA(M,IY,JZM)=.1*GLOB
  130 CONTINUE
C****
C**** Alternate hemispheric means (similar to alt. global mean)
C****
      IF(IAVGGH.EQ.1) GO TO 180
      DO 150 IHEM=1,2
      DO 140 IY=1,IYRS
      ANN(IY,IHEM+11)=XBAD
      IF(ANN(IY,IHEM+3).NE.XBAD.AND.ANN(IY,2*IHEM+7).NE.XBAD)
     *   ANN(IY,IHEM+11)=.4*ANN(IY,IHEM+3)+.6*ANN(IY,2*IHEM+7)
  140 CONTINUE
      DO 150 IY=1,IYRS
      DO 150 M=1,12
      DATA(M,IY,IHEM+11)=XBAD
      IF(DATA(M,IY,IHEM+3).NE.XBAD.AND.DATA(M,IY,2*IHEM+7).NE.XBAD)
     *  DATA(M,IY,IHEM+11)=.4*DATA(M,IY,IHEM+3)+.6*DATA(M,IY,2*IHEM+7)
  150 CONTINUE
C****
C**** Display the annual means
C****
  180 DO 190 JZ=1,JZP
  190 WRITE(6,'(I6,2X,A80)') JZ,TITLEZ(IORD(JZ))
      IYRSP=IYRS
      if(data(12,iyrs,14).gt.8000.) IYRSP=IYRS-1 ! skip incomplete year
      DO 200 IY=IY1TAB-IYRBEG+1,IYRSP
      IF( (IY+IYRBEG.gt.IY1TAB+5.and.MOD(IY+IYRBEG-2,20).EQ.0)
     * .or. IY.eq.IY1TAB-IYRBEG+1 ) THEN
        WRITE(96,*)
        WRITE(96,'(''                           24N   24S   90S   '',
     *        ''  64N   44N   24N   EQU   24S   44S   64S   90S'')')
        WRITE(96,'(''Year  Glob  NHem  SHem    -90N  -24N  -24S   '',
     *    '' -90N  -64N  -44N  -24N  -EQU  -24S  -44S  -64S Year'')')
      ENDIF
      IYR=IYRBEG+IY-1
      WRITE(96,'(I4,3(1X,I5),2X,3(1X,I5),2X,8(1X,I5),I5)') IYR,
     *  (NINT(100.*ANN(IY,IORD(JZ))),JZ=1,JZP),IYR
  200 CONTINUE
      WRITE(96,'(''Year  Glob  NHem  SHem     24N   24S   90S   '',
     *  ''  64N   44N   24N   EQU   24S   44S   64S   90S Year'')')
      WRITE(96,'(''                          -90N  -24N  -24S   '',
     *      '' -90N  -64N  -44N  -24N  -EQU  -24S  -44S  -64S'')')
      WRITE(96,*)
C****
C**** Display the monthly global and hemispheric means
C****
      DO 208 J=1,3
      WRITE(96+J,'(A10,'' Temperature Anomalies'',
     * '' in .01 C     base period: 1951-1980'')') TIT(J)
      DO 207 IY=IY1TAB-IYRBEG+1,IYRS
      IF( (IY+IYRBEG.gt.IY1TAB+5.and.MOD(IY+IYRBEG-2,20).EQ.0)
     *   .or. IY.eq.IY1TAB-IYRBEG+1  ) THEN
        WRITE(96+J,*)
        WRITE(96+J,'(''Year   Jan  Feb  Mar  Apr  May  Jun  Jul  Aug'',
     *''  Sep  Oct  Nov  Dec    J-D D-N    DJF  MAM  JJA  SON  Year'')')
      END IF
      awin=9999.
      if(iy.gt.1)
     *awin=data(12,IY-1,IORD(J))+data(1,IY,IORD(J))+data(2,IY,IORD(J))
      aspr=data(3,IY,IORD(J))+data(4,IY,IORD(J))+data(5,IY,IORD(J))
      asmr=data(6,IY,IORD(J))+data(7,IY,IORD(J))+data(8,IY,IORD(J))
      afl=data(9,IY,IORD(J))+data(10,IY,IORD(J))+data(11,IY,IORD(J))
      IOUT(15)=100*MBAD
      if(awin.lt.8000.) IOUT(15)=NINT(100.*awin/3.)
      IOUT(16)=100*MBAD
      if(aspr.lt.8000.) IOUT(16)=NINT(100.*aspr/3.)
      IOUT(17)=100*MBAD
      if(asmr.lt.8000.) IOUT(17)=NINT(100.*asmr/3.)
      IOUT(18)=100*MBAD
      if(afl.lt.8000.)  IOUT(18)=NINT(100.*afl/3.)
      IOUT(14)=100*MBAD
      ann2=awin+aspr+asmr+afl
      if(ann2.lt.8000.) IOUT(14)=NINT(100.*ann2/12.)
      IOUT(13)=100*MBAD
      ann1=ANN(IY,IORD(J))
      if(iy.eq.iyrs .and. data(12,IY,IORD(J)).gt.8000.) ann1=9999.
      if(ann1.lt.8000.) IOUT(13)=NINT(100.*ANN(IY,IORD(J)))
      DO 202 M=1,12
  202 IOUT(M)=NINT(100.*DATA(M,IY,IORD(J)))
      IYR=IYRBEG-1+IY
  207 WRITE(96+J,'(I4,1X,12I5,2X,I5,I4,2X,4I5,I6)') IYR,IOUT,IYR
      WRITE(96+J,'(''Year   Jan  Feb  Mar  Apr  May  Jun  Jul  Aug'',
     *''  Sep  Oct  Nov  Dec    J-D D-N    DJF  MAM  JJA  SON  Year'')')
  208 CONTINUE
C****
C**** Save annual means on disk
C****
      WRITE(11) INFOO,TITLEO,TITLE(29:80),BLANK(1:28)
      DO 210 JZ=1,JZM
      CALL SWRITE(11,MONM/12,ANN(1,JZ),ANNW(1,JZ),TITLEZ(JZ))
  210 CONTINUE
C****
C**** Resave monthly means on disk
C****
      REWIND 10
      WRITE(12) INFO,TITLE,TITL2
      DO 310 JZ=1,JZM
      CALL SWRITE(12,MONM,DATA(1,1,JZ),WT(1,1,JZ),TITLEZ(JZ))
  310 CONTINUE
      STOP
      END

      SUBROUTINE SREAD (IN,NDIM,AR,WT,TITLE)
      CHARACTER*80 TITLE
      DIMENSION AR(NDIM),WT(NDIM)
      READ (IN) AR,WT,TITLE
C     WRITE(6,*) 'in ',TITLE,ndim,(AR(n),n=ndim-11,ndim)
      RETURN
      END

      SUBROUTINE SWRITE (IOUT,NDIM,AR,WT,TITLE)
      CHARACTER*80 TITLE
      DIMENSION AR(NDIM),WT(NDIM)
      WRITE (IOUT) AR,WT,TITLE
C     WRITE(6,*) 'out',TITLE,ndim,(AR(n),n=ndim-11,ndim)
      RETURN
      END

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

Another old school program.

That’s it. The exposition of annzon will also have to wait until I’m done with earlier steps. Until then you can at least admire the source code ;-)

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.