w3fp10.f

Go to the documentation of this file.

C> @file
C> @brief Printer contour subroutine.
C> @author Ralph Jones @date 1989-09-08
 
C> Prints a two-dimensional grid of any shape, with
C> contouring, if desired. Grid values are scaled according to
C> to constants specified by the programer, rounded, and printed
C> as 4,3, or 2 digit integers with sign, the sign marking the
C> grid position of the printed number. If contouring is requested,
C> bessel's interpolation formula is used to optain the contour lines.
C> Contours are indicated by alphabetic characters ranging from a to
C> h or numeric characters from 0 to 9. Contour origin and interval
C> are specified by the programmer in terms of printed values.
C>
C> ### Program History Log
C> Date | Programmer | Comments
C> -----|------------|---------
C> 1989-09-08 | Ralph Jones | Initial
C> 1992-05-02 | Ralph Jones | Convert to cray cft77 fortran, add save.
C>
C> @param[in] RDATA Real array of grid data to be printed.
C> @param[in] KTBL Integer array with shape of array.
C> @param[in] CNST Real array of four elements, used in
C> scaling for printing and contouring.
C> @param[in] TITLE Is a array of 132 characters or less of
C> hollerith data, 1st char. must be blank.
C> printed at bottom of the map.
C> @param[in] KRECT 1 if grid is rectangular, 0 otherwise.
C> @param[in] KCONTR 1 for contouring , 0 otherwise.
C> @param[in] LINEV 0 is for 6 lines per vertical inch,
C> non-zero 8 lines per vertical inch.
C> @param[in] IWIDTH Number of characters in print line,
C> 132 is standard printer.
C>
C> Return conditions: Normal subroutine return, unless number of rows is
C> greater than 200, prints error message and exits.
C>
C> @note Special version of w3fp05(), 1st point is upper left hand
C> corner. Written on request of peter chase because some
C> grib fields can start with the upper left hand corner
C> as the 1st point of a grid.
C>
C> @author Ralph Jones @date 1989-09-08
      SUBROUTINE w3fp10(RDATA,KTBL,CNST,TITLE,KRECT,KCONTR,
     & LINEV,IWIDTH)
C
       REAL            CNST(4)
       REAL            RDATA(*)
       REAL            RWA(28)
       REAL            RWB(28)
       REAL            RWC(28)
       REAL            RWD(28)
       REAL            VDJA(29)
       REAL            VDJB(28)
       REAL            VDJC(28)
C
       INTEGER         TITLE(33)
       INTEGER         KRLOC(200)
       INTEGER         KTBL(*)
       INTEGER         OUTPUT
       INTEGER         PAGNL
       INTEGER         PAGNR
       INTEGER         PAGN3
       INTEGER         PCCNT
       INTEGER         PCFST
       INTEGER         PGCNT
       INTEGER         PGCNTA
       INTEGER         PGFST
       INTEGER         PGFSTA
       INTEGER         PGMAX
C
       LOGICAL         DONE
       LOGICAL         LCNTR
       LOGICAL         RECT
C
       CHARACTER*1     KALFA(16)
       CHARACTER*1     KALPH(20)
       CHARACTER*1     KHASTR
       CHARACTER*1     KHBLNK
       CHARACTER*1     KHDOLR
       CHARACTER*1     KHMNS
       CHARACTER*1     KHPLUS
       CHARACTER*1     KHRSTR
       CHARACTER*1     KHTBL(10)
       CHARACTER*1     KLINE(126)
       CHARACTER*1     KLINES(132)
       CHARACTER*1     KNUMB(20)
C
       equivalence(crmx,vdja(29))
       equivalence(kline(1),klines(8))
       equivalence(vdjc(1),rwa(1))
C
C      ... THE VALUE CRMX IS MACHINE DEPENDENT, IT SHOULD BE
C      ... SET TO A VALUE A LITTLE LESS THAN THE LARGEST POSITIVE
C      ... FLOATING  POINT NUMBER FOR THE COMPUTER.
C
       SAVE
C> The value CRMX is machine dependent, it should be
C> set to a value a little less than the largest positive
C> floating point number for the computer.
       DATA  crmx  /10.e70/
       DATA  kalfa /'A',' ','B',' ','C',' ','D',' ','E',' ','F',
     &              ' ','G',' ','H',' '/
       DATA  khastr/'*'/
       DATA  khblnk/' '/
       DATA  khdolr/'$'/
       DATA  khmns /'-'/
       DATA  khplus/'+'/
       DATA  khrstr/'1'/
       DATA  khtbl /'0','1','2','3','4','5','6','7','8','9'/
 
C> LIMNRW is limit on number of rows allowed and is dimension of KRLOC
       DATA  limnrw/200/
       DATA  knumb /'0',' ','1',' ','2',' ','3',' ','4',' ',
     &              '5',' ','6',' ','7',' ','8',' ','9',' '/
       DATA  output/6/
       DATA  r5    /.2/
       DATA  r50   /.02/
C
 8000  FORMAT (1h0,10x,44herror from w3fp10 ... number of rows in your,
     &         9h array = ,i4,24h which exceeds limit of ,i4)
 8100  FORMAT ( 1ht)
 8200  FORMAT ( 1hs)
 8300  FORMAT ( 1h ,/,1h ,/,1h )
 8400  FORMAT ( 1h ,/,1h )
 8500  FORMAT ( 132a1)
 8600  FORMAT ( 33a4)
C
C           COMPUTE VALUES FOR PRINTER WIDTH
C
         IF (iwidth.GE.132.OR.iwidth.LE.0) pgmax = 25
         IF (iwidth.GE.1.AND.iwidth.LE.22) pgmax = 3
         IF (iwidth.GT.22.AND.iwidth.LT.132) pgmax = (iwidth-7) / 5
         lw            = (pgmax * 5 + 7) / 4
         pagn3         =  pgmax + 3
         vdja(pagn3+1) =  crmx
         mxpg          =  pgmax * 5 + 7
C
         IF (linev .NE. 0) THEN
C
C     ...OTHERWISE LINEV IS NON-ZERO, SO 8 LINES/INCH IS DESIRED...
C
           linate = 1
           r4     = 0.250
           r32    = 0.03125
           con2   = 10.0
           nbtwn  = 3
C
         ELSE
C
           linate = 2
           r4     = 0.33333333
           r32    = 1.0 / 18.0
           con2   = 6.0
           nbtwn  = 2
         ENDIF
C
           pgcnta = 0
           pgfsta = 0
           rect   = .false.
           done   = .false.
           kz     = 0
           kza    = 1000
           a      = cnst(1)
           kca    = 2 * (1 - krect)
C
C        TO SET NO. OF DIGITS TO BE PRINTED
C        WHICH IS A FUNCTION OF THE TENS POSITION IN KCONTR
C
         nodig = iabs(kcontr/10)
         nodig = 3 - nodig
C
C        WHERE C(NODIG) + 1 IS NO. OF DIGITS TO BE PRINTED
C
         IF (nodig.LT.1 .OR. nodig.GT.3) nodig = 3
C
C        ANY OUT-OF-RANGE WILL GET 4 DIGITS
C
         lcntr = .false.
         nconq = iabs(mod(kcontr,10))
         IF (nconq .EQ. 0) GO TO 400
         IF (nconq .LE. 2) GO TO 300
C
C        OTHERWISE RESET NCONQ
C
         nconq = 0
         GO TO 400
C
 300     CONTINUE
         lcntr = .true.
C
C        WITH NCONQ = 1 FOR LETTERS,AND = 2 FOR NUMBERS IN CONTOUR BANDS
C
 400     CONTINUE
         IF (nconq .NE. 2) THEN
C
C        OTHERWISE SET AS LETTERS
C
           kcow = 16
           DO 500  j = 1,kcow
             kalph(j) = kalfa(j)
 500       CONTINUE
C
         ELSE
C
           kcow = 20
           DO 700  j = 1,kcow
             kalph(j) = knumb(j)
 700       CONTINUE
C
         ENDIF
C
           radj = 4 * kcow
           kd   = 1
C
C *** SET UP TABLE OF INDICES CORRESPONDING TO FIRST ITEM IN EACH ROW
C ***     THIS IS KRLOC
C *** PICK OUT SIZE AND ROW NUMBER OF LARGEST ROW (KCMX AND KCLMX)
C *** KZA LEFT-JUSTIFIES MAP IF ALL ROWS HAVE COMMON MINIMAL OFFSET
C
         IF (ktbl(1 ).NE.(-1)) THEN
C
C *** ONE-DIMENSIONAL FORM
C
         ktf  = 3
         kza  = 0
         imin = ktbl(2)
         jmin = ktbl(3)
         jmax = ktbl(3) + ktbl(1) - 1
         nrws = ktbl(1)
         IF (nrws .GT. limnrw) THEN
           WRITE (output,8000) nrws , limnrw
           RETURN
         ENDIF
         kc   = 1
C
         DO 1000 j = 1,nrws
           krloc(j) = kd
           IF (ktbl(kc+4) + ktbl(kc+3).LE.kz ) GO TO 900
             kclmx = j
             imax  = ktbl(kc+4) + ktbl(kc+3)
             kz    = imax
             kcmx  = krloc(j) + ktbl(kc+4)
 900       CONTINUE
             kd    = kd + ktbl(kc+4)
             kc    = kc + kca
 1000    CONTINUE
C
         ELSE
C
C *** TWO-DIMENSIONAL FORM
C *** THE TWO-DIMENSIONAL FORM IS COMPILER-DEPENDENT
C *** IT DEPENDS ON THE TWO-DIMENSIONAL ARRAY BEING STORED COLUMN-WISE
C *** THAT IS, WITH THE FIRST INDEX VARYING THE FASTEST
C
           imin = ktbl(6)
           jmin = ktbl(7)
           nrws = ktbl(5)
           IF (nrws .GT. limnrw) THEN
             WRITE (output,8000) nrws , limnrw
             RETURN
           ENDIF
C
         jmax  = ktbl(7) + ktbl(5) -1
         kc    = 1
         DO 1500 j = 1,nrws
           krloc(j) = ktbl(2) * (ktbl(4)-nrws+j-1) + ktbl(kc+7) + 1
           IF (ktbl(kc+7) + ktbl(kc+8).LE.kz) GO TO 1400
             imax  = ktbl(kc+7) + ktbl(kc+8)
             kz    = imax
             kcmx  = krloc(j) + ktbl(kc+8)
             kclmx = j
 1400      CONTINUE
           IF (ktbl(kc+7).LT.kza) kza = ktbl(kc+7)
             kc = kc + kca
 1500    CONTINUE
           imax = imax - kza
           ktf  = 7
         ENDIF
C
           pagnl = 0
           pagnr = pgmax
           IF (.NOT.lcntr) GO TO 1700
           adc   = (cnst(1) - cnst(4)) / cnst(3) + radj
           bc    =  cnst(2) / cnst(3)
C
C *** PRINT I-LABELS ACROSS TOP OF MAP
C
 1700    CONTINUE
C
C ***    WHICH PREPARES CDC512 PRINTER FOR 8 LINES PER INCH
C
         IF (linate.EQ.1) WRITE (output,8100)
C
C     ...WHICH PREPARES PRINTER FOR 6 LINES PER INCH
C
         IF (linate.EQ.2) WRITE (output,8200)
           klines(1) = khrstr
           kbr       = 1
           GO TO 6900
C
 1800    CONTINUE
         IF (.NOT.lcntr) GO TO 2000
C
C *** INITIALIZE CONTOUR WORKING AREA
C
           DO 1900 j = 1,pagn3
             rwc(j) = crmx
             rwd(j) = crmx
 1900      CONTINUE
C
C *** SET UP CONTOUR DATA AND PAGE LIMITERS FOR FIRST TWO ROWS
C
 2000    CONTINUE
           kra  = 1
           kc   = ktf + 1
           kbr  = 2
           GO TO 5900
C
 2100    CONTINUE
           kra = 2
           kc  = kc + kca
           kbr = 3
           GO TO 5900
C
 2200    CONTINUE
           kr   = 0
C
C *** TEST IF THIS IS LAST PAGE
C
           IF (imax.GT.pgmax-1) GO TO 2300
           lmr  = imax * 5 + 2
           done = .true.
C
C *** DO LEFT J-LABELS
C
 2300    CONTINUE
            jcurr = jmin
C
 2400    CONTINUE
           kr  = kr + 1
           kra = kr + 2
           kc  = kc + kca
           kta = mod(jcurr,10)
           ktb = mod(jcurr,100)/10
           ktc = mod(jcurr,1000)/100
           IF (kr .EQ. 1 .OR. (.NOT. lcntr)) GO TO 2500
           GO TO  2600
C
 2500    CONTINUE
           IF (linate.EQ.1) WRITE (output,8300)
           IF (linate.EQ.2) WRITE (output,8400)
C
 2600    CONTINUE
         klines(2) = khplus
         klines(1) = khblnk
         IF (jcurr.LT.0) klines(2) = khmns
         kta       = iabs(kta)
         ktb       = iabs(ktb)
         ktc       = iabs(ktc)
         IF (ktc .EQ. 0) GO TO 2700
         klines(3) = khtbl(ktc+1)
         klines(4) = khtbl(ktb+1)
         klines(5) = khtbl(kta+1)
         GO TO 2800
C
 2700    CONTINUE
         klines(3) = khtbl(ktb+1)
         klines(4) = khtbl(kta+1)
         klines(5) = khblnk
C
 2800    CONTINUE
         DO 2900  j = 6,mxpg
           klines(j) = khblnk
 2900    CONTINUE
         IF (.NOT.done) GO TO 3000
C
C *** DO RIGHT J-LABELS IF LAST PAGE OF MAP
C
           kline(lmr)   = klines(2)
           kline(lmr+1) = klines(3)
           kline(lmr+2) = klines(4)
           kline(lmr+3) = klines(5)
C
C *** FETCH AND CONVERT GRID VALUES TO A1 FORMAT FOR WHOLE LINE
C
 3000    CONTINUE
         krx = krloc(kr)
         klx = 5 * pgfst + 1
         IF (pgcnt.EQ.0) GO TO 4000
           DO 3800 kk = 1,pgcnt
             temp       = rdata(krx) * cnst(2) + a
             ktemp      = abs(temp) + 0.5
             kline(klx) = khplus
             IF (temp.LT.0.0) kline(klx) = khmns
             GO TO (3300,3200,3100),nodig
 3100      CONTINUE
             kta = mod(ktemp,10000)/1000
C
 3200      CONTINUE
             ktb = mod(ktemp,1000)/100
C
 3300      CONTINUE
             ktc = mod(ktemp,100)/10
             ktd = mod(ktemp,10)
             GO TO (3400,3500,3600),nodig
C
 3400      CONTINUE
             kline(klx+1) = khtbl(ktc+1)
             kline(klx+2) = khtbl(ktd+1)
             GO TO  3700
C
 3500      CONTINUE
             kline(klx+1) = khtbl(ktb+1)
             kline(klx+2) = khtbl(ktc+1)
             kline(klx+3) = khtbl(ktd+1)
             GO TO  3700
C
 3600      CONTINUE
             kline(klx+1) = khtbl(kta+1)
             kline(klx+2) = khtbl(ktb+1)
             kline(klx+3) = khtbl(ktc+1)
             kline(klx+4) = khtbl(ktd+1)
C
 3700    CONTINUE
           klx = klx + 5
           krx = krx + 1
 3800    CONTINUE
C
C *** FOLLOWING CHECKS FOR POLE POINT AND INSERTS PROPER CHARACTER.
C
         IF (jcurr.NE.0) GO TO 4000
         IF (imin.LT.(-25).OR.imin.GT.0) GO TO 4000
         kx = -imin
         IF (kx.LT.pgfst.AND.kx.GT.pgcnt+pgfst) GO TO 4000
         kx = 5 * kx
         IF (kline(kx+1).EQ.khmns) GO TO 3900
         kline(kx) = khdolr
         GO TO 4000
C
 3900    CONTINUE
         kline(kx+1) = khastr
C
C *** PRINT LINE OF MAP DATA
C
 4000    CONTINUE
           WRITE (output,8500) (klines(ii),ii=1,mxpg)
           krloc(kr) = krx
           jcurr     = jcurr + 1
C           JCURR     = JCURR + JRWMP
C
C *** TEST BOTTOM OF MAP
C
         IF (kr.EQ.nrws) GO TO 5700
C
C *** SET UP CONTOUR DATA AND PAGE LIMITERS FOR NEXT ROW
C
         kbr = 4
         GO TO 5900
C
 4100    CONTINUE
         IF (.NOT.lcntr) GO TO 2400
C
C *** DO CONTOURING
C
         DO 4200 jj = 1,mxpg
             klines(jj) = khblnk
 4200    CONTINUE
C
C *** VERTICAL INTERPOLATIONS
C
         DO 4700 kk = 1,pagn3
           IF (rwb(kk).LT.crmx.AND.rwc(kk).LT.crmx) GO TO  4300
           vdjb(kk) = crmx
           vdjc(kk) = crmx
           GO TO 4600
C
 4300    CONTINUE
           IF (rwa(kk).LT.crmx.AND.rwd(kk).LT.crmx) GO TO  4400
           vdjc(kk) = 0.
           GO TO  4500
C
 4400    CONTINUE
           vdjc(kk) = r32*(rwa(kk)+rwd(kk)-rwb(kk)-rwc(kk))
C
 4500    CONTINUE
           vdjb(kk) = r4*(rwc(kk)-rwb(kk)-con2*vdjc(kk))
C
 4600    CONTINUE
           vdja(kk)=rwb(kk)
C
 4700    CONTINUE
C
C     ...DO 2 OR 3 ROWS OF CONTOURING BETWEEN GRID ROWS...
C
         DO 5600 ll = 1,nbtwn
           DO 4800 kk = 1,pagn3
             vdjb(kk) = vdjc(kk) + vdjb(kk)
             vdja(kk) = vdjb(kk) + vdja(kk)
 4800      CONTINUE
C
C     ...WHERE VDJA HAS THE INTERPOLATED VALUE FOR THIS INTER-ROW
C *** HORIZONTAL INTERPOLATIONS
C
           hdc = 0.0
           IF (vdja(1).GE.crmx) GO TO 4900
           hdc = r50*(vdja(4)+vdja(1)-vdja(2)-vdja(3))
C
 4900    CONTINUE
         kxb = 0
         DO 5200 kk = 1,pgmax
           IF (vdja(kk+1).GE.crmx) GO TO 5100
           hda = vdja(kk+1)
           IF (vdja(kk+2).GE.crmx) GO TO 5500
           IF (vdja(kk+3).GE.crmx) hdc = 0.0
           hdb = r5 * (vdja(kk+2) - vdja(kk+1) - 15.0 * hdc)
C
C *** COMPUTE AND STORE CONTOUR CHARACTERS, 5 PER POINT
C
           khda = hda
           kdb = iabs(mod(khda,kcow))
           kline(kxb+1) = kalph(kdb+1)
           DO 5000 jj = 2,5
             hdb  = hdb + hdc
             hda  = hda + hdb
             khda = hda
             kdb  = iabs(mod(khda,kcow))
             kxa  = kxb + jj
             kline(kxa) = kalph(kdb+1)
 5000    CONTINUE
         hdc = r50*(vdja(kk+4)+vdja(kk+1)-vdja(kk+2)-vdja(kk+3))
         IF (vdja(kk+4).GE.crmx) hdc = 0.0
C
 5100    CONTINUE
           kxb = kxb + 5
C
 5200    CONTINUE
C
 5300    CONTINUE
         WRITE (output,8500) (klines(ii),ii=1,mxpg)
         DO 5400 kk = 1,mxpg
           klines(kk) = khblnk
 5400    CONTINUE
         GO TO 5600
C
 5500    CONTINUE
           khda = hda
           kdb  = iabs(mod(khda,kcow))
           kline(kxb+1) = kalph(kdb+1)
           GO TO 5300
C
 5600    CONTINUE
         GO TO 2400
C
 5700    CONTINUE
         IF (linate.EQ.1)  WRITE (output,8300)
         IF (linate.EQ.2)  WRITE (output,8400)
         klines(1) = khblnk
C
C *** PRINT I-LABELS ACROSS BOTTOM OF PAGE
C
         kbr = 5
         GO TO 6900
C
 5800    CONTINUE
         IF (linate.EQ.1)  WRITE (output,8300)
         IF (linate.EQ.2)  WRITE (output,8400)
C
C *** PRINT TITLE
C
         WRITE (output,8600) (title(ii),ii=1,lw)
C
C ***    TEST END OF MAP
C
         IF (krloc(kclmx).EQ.kcmx) RETURN
C
C ***    ADJUST PAGE LINE BOUNDARIES
C
         IF (imax.GT.pgmax) imax = imax - pgmax
         imin  = ka
         pagnl = pagnl + pgmax
         pagnr = pagnr + pgmax
         GO TO   1700
C
C *** ROUTINE TO PRE-STORE ROWS FOR CONTOURING AND COMPUTE LINE LIMITERS
C
 5900    CONTINUE
         pgfst = pgfsta
         pgcnt = pgcnta
         IF (kra.GT.nrws) GO TO 6800
           krfst = ktbl(kc) - kza
           krcnt = ktbl(kc+1)
           kfx   = krloc(kra)
           IF (rect) GO TO  6100
           IF (krfst-pagnl.LE.(-1)) GO TO 6400
           pcfst  = krfst - pagnl + 1
           IF (pcfst.GE.pagn3) GO TO 6700
           pgfsta = pcfst-1
           pccnt  = min(pagnr-krfst+2,krcnt)
           IF (pgfsta.EQ.0) GO TO 6600
           pgcnta = min(pagnr-krfst,krcnt)
           IF (pgcnta.GT.0) GO TO 6000
           pgcnta = 0
           GO TO 6100
C
 6000    CONTINUE
           rect = krect.EQ.1.AND.pgcnta.LE.krcnt
C
 6100    CONTINUE
         IF (.NOT.lcntr) GO TO (1800,2100,2200,4100,5800) kbr
         DO 6200 kk = 1,pagn3
           rwa(kk) = rwb(kk)
           rwb(kk) = rwc(kk)
           rwc(kk) = rwd(kk)
           rwd(kk) = crmx
 6200    CONTINUE
C
         IF (pccnt.EQ.0) GO TO (1800,2100,2200,4100,5800) kbr
           kpc = pcfst + 1
         DO 6300 kk = 1,pccnt
           rwd(kpc) = rdata(kfx) * bc + adc
           kfx = kfx + 1
           kpc = kpc + 1
 6300    CONTINUE
         GO TO (1800,2100,2200,4100,5800) kbr
C
 6400    CONTINUE
           pcfst  = 0
           pgfsta = 0
           kfx    = kfx - 1
           pccnt  = krfst + krcnt - pagnl + 1
           IF (pccnt.LT.pagn3) GO TO 6500
           pccnt  = pagn3
           pgcnta = pgmax
           GO TO 6100
C
 6500    CONTINUE
           IF (pccnt.GT.0) GO TO 6600
           pgcnta = 0
           pccnt  = 0
           GO TO 6100
C
 6600    CONTINUE
           pgcnta = min(pgmax,krcnt+krfst-pagnl)
           GO TO 6100
C
 6700    CONTINUE
           pgcnta = 0
C
 6800    CONTINUE
           pccnt  = 0
           GO TO 6100
C
C *** ROUTINE TO PRINT I-LABELS
C
 6900    CONTINUE
         DO 7000 kk = 2,mxpg
           klines(kk) =  khblnk
 7000    CONTINUE
C
         kk  = 1
         ka  = imin
         lbl = min(imax,pgmax)
C
         DO 7300 jj = 1,lbl
           kline(kk) = khplus
           IF (ka.LT.0) kline(kk) = khmns
           kta = iabs(mod(ka,100)) / 10
           ktb = iabs(mod(ka,10))
           ktc = iabs(mod(ka,1000)) / 100
           IF (ktc .EQ. 0) GO TO 7100
           kline(kk+1) = khtbl(ktc+1)
           kline(kk+2) = khtbl(kta+1)
           kline(kk+3) = khtbl(ktb+1)
           GO TO  7200
C
 7100      CONTINUE
             kline(kk+1) = khtbl(kta+1)
             kline(kk+2) = khtbl(ktb+1)
C
 7200    CONTINUE
           kk = kk + 5
           ka = ka + 1
C
 7300    CONTINUE
C
         WRITE (output,8500) (klines(ii),ii=1,mxpg)
C
         GO TO (1800,2100,2200,4100,5800) kbr
C
 7400    CONTINUE
           RETURN
C
       END