w3fp04.f

Go to the documentation of this file.

C> @file
C> @brief Print array of data points at lat/lon points.
C> @author J. Horodeck @date 1980-01-15
 
C> Given an array of meteorological data and corresponding
C> latitude/longitude position for each data point, these data
C> values are printed at their approximate latitude/longitude
C> positions on a polar stereographic projection.
C>
C> PROGRAM HISTORY LOG:
C> - J. Horodeck 1980-01-15
C> - Ralph Jones 1985-07-31 Change to cdc fortran 200
C> - Ralph Jones 1990-08-15 Change to cray cft77 fortran
C>
C> @param[in] IFLD Real or integer fullword array of data points.
C> @param[in] ALAT Real array of latitude positions (>0 for nh,
C> <0 for sh) for the data to be plotted.
C> @param[in] ALON Real array of longitudes (west of greenwich)
C> @param[in] TITLE Integer size 10 alphanumeric array of
C> characters for title to be written on printout.
C> @param[in] IDIM Integer number of data values to plot (size of
C> arrays ifld, alat and alon).
C> @param[in] CMIL Real left side of grid - minimum coarse mesh
C> i coordinate (minimum value of 1.0).
C> @param[in] CMIR Real right side of grid - maximum coarse mesh
C> i coordinate (maximum value of 65.0).
C> @param[in] CMJB Real bottom of grid - minimum coarse mesh
C> j coordinate (minimum value of 1.0).
C> @param[in] CMJT Real top of grid - maximum coarse mesh j
C> coordinate (maximum value of 65.0).
C> @param[in] INUM Integer three digit number for the following:
C> - Hundreds digit = type of data
C>  - 1 = Fixed point
C>  - 2 = Floating point
C>  - 3 = Alphanumeric
C> - Tens digit = hemispheric reference
C>  - 0 = Northern hemisphere
C>  - 1 = Southern hemisphere
C> - Units digit = number of characters to plot
C> - Minimum = 1 character
C> - Maximum = 4 characters
C> @param[in] XFAC Real map scale factor (desired map scale = xfac
C> * 1:30,000,000 (standard nmc 65x65 grid scale))
C> @param[out] IERR Integer return code.
C>
C> @note Because this code could produce considerable output
C> the subset of the nmc 65x65 grid which can be printed is a
C> function of the map scale factor, e.g. for xfac=5 the maximum
C> range of i and j is 27.0, for xfac=2 the range is 64.0.
C>
C> @author J. Horodeck @date 1980-01-15
      SUBROUTINE w3fp04(IFLD,ALAT,ALON,TITLE,IDIM,CMIL,CMIR,
     &                   CMJB,CMJT,INUM,XFAC,IERR)
C
       REAL            ALAT(IDIM), ALON(IDIM)
C
       INTEGER         IFLD(IDIM), TITLE(10)
       INTEGER         LINE(24), IL(17), IR(17), IJU(20), IJL(20)
C
       LOGICAL         A
C
C
       CHARACTER*1     KH(120,77), MEAN(4), KB, KM, KP, LC,
     &                 kk(5,77,20), DATA(4), l1, l2, l3,
     &                 ipole(4), kn, ks
       CHARACTER*4     L24, L116
       CHARACTER*8     IFMTT
       CHARACTER*24    FMT1
       CHARACTER*28    FMT2
       CHARACTER*24    FMT4
C
       equivalence(mean(1),imean), (DATA(1),lfld)
       equivalence(rfield,ifield)
       equivalence(ifmtt,ifmt)
C
       DATA   jjaa /116/
       DATA   jjbb / 77/
       DATA   fmt1 /"(6X,    ('+',I , X),//) "/
       DATA   fmt2 /"(' +',I3,1X,    A1,' +',I3) "/
       DATA   fmt4 /"(//, 6X,    ('+',I , X))"/
       DATA   kb   /' '/, km/'-'/, kp/'+'/, lc/'X'/
       DATA   l1/'1'/, l2/'2'/, l3/'3'/, l24/'  24'/, l116/' 116'/
       DATA   ipole/'P','O','L','E'/,  kn/'N'/, ks/'S'/
C
 1001  FORMAT('1',16x,'PANEL #',i2,' OF ',i2,4x,10a8,/,/)
 1003  FORMAT(6x,116a1)
 2001  FORMAT(///,20x,'UPPER LEFT CORNER--LAT =',f6.2,' LON =',f7.2,'W'
     &        , 3x,'UPPER RIGHT CORNER--LAT =',f6.2,' LON =',f7.2,'W')
 2002  FORMAT(20x,'LOWER LEFT CORNER--LAT =',f6.2,' LON =',f7.2,'W'
     &        , 3x,'LOWER RIGHT CORNER--LAT =',f6.2,' LON =',f7.2,'W')
 2003  FORMAT(/,/, 16x, 'PANEL #', i2, ' OF ', i2, 4x, 10a8)
 9001  FORMAT(/,5x,'CMIL = ',f8.1,' CMIR = ',f8.1,' HIGH AND LOW'
     &         ,' VALUES REVERSED......RETURN......')
 9002  FORMAT(/,5x,'CMJB = ',f8.1,' CMJT = ',f8.1,' HIGH AND LOW'
     &         ,' VALUES REVERSED......RETURN......')
 9003  FORMAT(/,5x,f8.1,' IS ILLEGAL VALUE FOR LOW I. IT IS NOW 1.0')
 9004  FORMAT(/,5x,f8.1,' IS ILLEGAL VALUE FOR HIGH I. IT IS NOW 65.0')
 9005  FORMAT(/,5x,f8.1,' IS ILLEGAL VALUE FOR LOW J. IT IS NOW 1.0')
 9006  FORMAT(/,5x,f8.1,' IS ILLEGAL VALUE FOR HIGH J. IT IS NOW 65.0')
 9007  FORMAT(/,5x,'REQUESTED NUMBER OF CHARACTERS TO PLOT(',i2,' )IS'
     &         ,' NOT ALLOWED. FOUR(4) IS MAXIMUM. THATS ALL YOU GET')
 9008  FORMAT(/,5x,'REQUESTED SUBSET OF 65X65 GRID CANNOT CURRENTLY '
     &         ,'BE PLOTTED WITH MAP SCALE FACTOR',f5.1,/5x,'IF PLOT '
     &         ,'IS NECESSARY, CONTACT JOHN M. HORODECK,ESQ. NMC/DD'
     &         ,'/SEB FOR ASSISTANCE')
 9009  FORMAT(/,5x,i4,' IS INVALID HEMISPHERIC REFERENCE'
     &         , '......RETURN......')
 9010  FORMAT(/,5x,'HUNDREDS DIGIT OF INUM(INUM =',i4,') IS'
     &         , ' INVALID......RETURN......')
C
C    TEST I,J VALUES FOR RANGE AND ORDER
C
         IF  (cmir.GT.cmil)  GO TO 1
          ierr = 1
           print 9001, cmil, cmir
           RETURN
    1    CONTINUE
         IF  (cmjt.GT.cmjb)  GO TO 2
           ierr = 1
           print 9002, cmjb, cmjt
           RETURN
    2    CONTINUE
         IF  (cmil.GE.1.0)  GO TO 3
           print 9003, cmil
           cmil = 1.0
    3    CONTINUE
         IF  (cmir.LE.65.0)  GO TO 4
           print 9004, cmir
           cmir = 65.0
    4    CONTINUE
         IF  (cmjb.GE.1.0)  GO TO 5
           print 9005, cmjb
           cmjb = 1.0
    5    CONTINUE
         IF  (cmjt.LE.65.0)  GO TO 6
           print 9006, cmjt
           cmjt = 65.0
    6    CONTINUE
C
C    CALCULATE VARIOUS LIMITS
C
         lnum = mod(inum,10)
         nref = (mod(inum,100))/10
C
C    TEST FOR INCORRECT ARGUMENTS PASSED
C
         IF  (lnum.LE.4)  GO TO 7
           print 9007, lnum
           lnum = 4
    7    CONTINUE
         IF  (nref.LE.1)  GO TO 8
           ierr = 1
           print 9009, nref
           RETURN
    8    CONTINUE
         IF  ((inum/100).LE.3)  GO TO 81
           ierr = 1
           print 9010, inum
           RETURN
   81    CONTINUE
C
         lnump1 = lnum + 1
         i1 = (cmil-1.0)*xfac + 1.0
         i2 = (cmir-1.0)*xfac + 1.0
         j1 = (cmjb-1.0)*xfac + 1.0
         j2 = (cmjt-1.0)*xfac + 1.0
C
C    WILL THIS PLOT BE TOO LARGE?
C
         IF  (((i2-i1).LT.139).AND.((j2-j1).LT.139))  GO TO 9
           ierr = 1
           print 9008, xfac
           RETURN
    9    CONTINUE
C
         offi   = i1 - 1
         offj   = j1 - 1
         jja    = (i2-i1)*5 + 1
         jjb    = (j2-j1)*4 + 1
         jjam1  = jja - 1
         jjbbm1 = jjbb - 1
         jjaam1 = jjaa - 1
         jjaapn = jjaa + lnum
         ibegin = lnump1 + 1
         ipage  = (jjam1/jjaa) + 1
         jpage  = (jjb/jjbb) + 1
         xmesh  = 381.0/xfac
         xip    = 32.0*xfac + 1.0
         xjp    = 32.0*xfac + 1.0
         iixip  = (xip-offi)*5 - 4
         jjxjp  = (xjp-offj)*4 - 3
C
C    PLOT DATA ONE PANEL AT A TIME IN SECTIONS
C
         DO  150  nx=1,ipage
         a = .false.
C
C    SET LIMITS OF I TO BE PRINTED
C
           il(nx) = i1 + (23*(nx-1))
           IF  (nx.NE.ipage)  ir(nx) = i1 + (23*nx)
           IF  (nx.EQ.ipage)  ir(nx) = i2
           imax   = ir(nx) - offi
           imin   = il(nx) - offi
           m      = 0
C
C    FILL ARRAY WITH VALUES OF I TO BE PRINTED AT TOP OF PAGE
C
           DO  10  i = imin,imax
             m       = m + 1
             line(m) = i
   10      CONTINUE
C
C    CALCULATE WIDTH OF PANEL IN INTEGERS AND
C    CHARACTERS FROM WHICH DETERMINE FORMAT
C    FIELD COUNT AND CONVERT BINARY TO ASCII
C
C    PRINT TOP LINE OF I
C
           la  = (imax-imin) + 1
           mmm = (la*5) - 4
           IF  (la.EQ.24)  GO TO 13
             CALL w3ai15(la,ifmt,1,4,kp)
             fmt1(5:8)   = ifmtt(1:4)
             fmt4(9:12)  = ifmtt(1:4)
             CALL w3ai15(mmm,ifmt,1,4,kp)
             fmt2(13:16) = ifmtt(1:4)
             GO TO 16
   13      CONTINUE
             fmt1(5:8)   = l24
             fmt2(13:16) = l116
             fmt4(9:12)  = l24
   16      CONTINUE
           IF  (la.LT.100) GO TO 19
             fmt1(15:15) = l3
             fmt1(17:17) = l1
             fmt4(19:19) = l3
             fmt4(21:21) = l1
             GO TO 22
   19      CONTINUE
             fmt1(15:15) = l2
             fmt1(17:17) = l2
             fmt4(19:19) = l2
             fmt4(21:21) = l2
   22      CONTINUE
           print 1001, nx, ipage, title
           WRITE(6,fmt1) (line(n), n=1,la)
C
C    PREPARE TO PRINT SECTIONS OF EACH PANEL
C
           DO  140  jnx=1,jpage
C
C    SET LIMITS OF J TO BE PRINTED
C
             iju(jnx) = j2 - (19*(jnx-1))
             IF  (jnx.NE.jpage) ijl(jnx) = j2 - (19*jnx)
             IF  (jnx.EQ.jpage) ijl(jnx) = j1
             jmax = iju(jnx) - offj
             jmin = ijl(jnx) - offj
             ju   = jjb - (4*jmax-3)
             jl   = jjb - (4*jmin-3)
             nnn  = (jmax-jmin)*4 + 1
C
C    FILL CHARACTER ARRAY WITH BLANKS AND PUT X MARKERS IN CORNERS
C    IF FIRST PANEL BLANK ENTIRE AREA,
C    OTHERWISE TRANSFER FIRST INUM I BYTES TO LARGE ARRAY
C    AND BLANK REMAINING ARRAY
C
             DO  37  j=1,jjbb
               IF  (nx.NE.1)  GO TO 31
                 DO  28  i=1,jjaapn
                   kh(i,j) = kb
   28            CONTINUE
                 GO TO 37
   31          CONTINUE
               DO  32  i=1,lnump1
                 kh(i,j) = kk(i,j,jnx)
   32          CONTINUE
               DO  34  i=ibegin,jjaapn
                 kh(i,j) = kb
   34          CONTINUE
   37        CONTINUE
             IF  (jnx.NE.1) GO TO 40
               kh(1,jjbb) = lc
               kh(mmm,jjbb) = lc
 200   CONTINUE
   40        CONTINUE
             IF  (jnx.NE.jpage)  GO TO 50
               kh(1,1) = lc
               kh(mmm,1) = lc
   50        CONTINUE
C
C    LOOP TO PUT DATA IN CHARACTER ARRAY
C
             DO  110  i=1,idim
C
C    TEST FOR BAD GEOGRAPHY
C
               IF  ((abs(alat(i)).GT.90.).OR.(alon(i).LT.0.0).OR.(alon
     a               (i).GT.360.0))  GO TO 90
C
C    CHANGE LAT,LON TO I,J
C
                 IF  (nref.EQ.0)  GO TO 51
                   CALL w3fb04(alat(i),alon(i),-xmesh,260.0,deli,delj)
                   GO TO 52
   51            CONTINUE
                   CALL w3fb04(alat(i),alon(i),xmesh,80.0,deli,delj)
   52            CONTINUE
                 xi = xip + deli
                 xj = xjp + delj
C
C    POSITION I,J COORDINATES IN CHARACTER ARRAY AND TEST
C    IF VALUES RETURNED ARE WITHIN LIMITS OF MAP AND WITHIN SECTIONS
C
                 ii = 1.0 + (xi-offi-0.9001)*5.0
                 jj = 1.0 + (xj-offj-0.8751)*4.0
                 iw = (jjaam1*(nx-1)) + 1
                 ix = (jjaam1*nx) + 1
                 iy = jjb - (jjbbm1*(jnx-1))
                 IF  (jnx.NE.jpage)  iz = jjb - (jjbbm1*jnx)
                 IF  (jnx.EQ.jpage)  iz = 1
                 IF  ((ii.LT.1).OR.(ii.GT.jja))  GO TO 100
                 IF  ((jj.LT.1).OR.(jj.GT.jjb))  GO TO 100
                 IF  ((ii.LT.iw).OR.(ii.GT.ix))  GO TO 100
                 IF  ((jj.GT.iy).OR.(jj.LT.iz))  GO TO 100
C
C    WRITE N+POLE IF IN THIS SECTION
C
                 IF  (.NOT.((iixip.GE.iw.AND.iixip.LE.ix).AND.
     a                 (jjxjp.LE.iy.AND.jjxjp.GE.iz)))  GO TO 56
                   iixxp = iixip - (jjaam1*(nx-1))
                   jjxxp = jjxjp - (iz-1)
                   IF  (nref.EQ.0)  kh(iixxp-1,jjxxp) = kn
                   IF  (nref.EQ.1)  kh(iixxp-1,jjxxp) = ks
                   kh(iixxp,jjxxp) = kp
                   DO  53  l=1,4
                     kh(iixxp+l,jjxxp) = ipole(l)
   53              CONTINUE
   56            CONTINUE
C
C    CONVERT CHARACTER ARRAY COORDINATES FROM
C    TOTAL MAP VALUES TO SECTION VALUES
C
                 ii = ii - (jjaam1*(nx-1))
                 IF  (jnx.NE.jpage)  jj = jj - (iz-1)
C
C    IF SPACE IS OCCUPIED SKIP THIS STATION
C
                 jnum = lnum + 1
                 DO  70  ik=1,jnum
                   in = ik - 1
                   IF  (kh(ii+in,jj).EQ.kb)  GO TO 60
                     GO TO 110
   60              CONTINUE
   70            CONTINUE
C
C    PLACE VALUE TO BE PLOTTED IN CHARACTER ARRAY
C
                 ifield = ifld(i)
C
C    TEST FOR TYPE OF DATA
C
                 IF  ((inum/100).EQ.3)  GO TO 82
                   IF  ((inum/100).EQ.1)  GO TO 73
                     jfld = rfield
                     GO TO 76
   73              CONTINUE
                   jfld = ifield
   76              CONTINUE
C
C    IF ORIGINALLY FIXED POINT OR HAS BEEN CONVERTED
C    FROM FLOATING POINT TO FIXED POINT
C
                   IF  ((jfld/10000).GE.1) jfld = mod(jfld,10000)
                   iiabs = iabs(jfld)
                   CALL w3ai15(iiabs,imean,1,lnum,kp)
                   IF  (jfld.LT.0)  kh(ii,jj) = km
                   IF  (jfld.GE.0)  kh(ii,jj) = kp
                   DO  79  ia=1,lnum
                     kh(ii+ia,jj) = mean(ia)
   79              CONTINUE
                   GO TO 110
   82            CONTINUE
C
C    FOR ALPHANUMERIC DATA
C
                 lfld = ifld(i)
                 kh(ii,jj) = kp
                 DO  85  iq=1,lnum
                   kh(ii+iq,jj) = DATA(iq)
   85            CONTINUE
   90          CONTINUE
  100          CONTINUE
  110        CONTINUE
             jjn = 0
C
C    PRINT JTH ROW AND VALUES OF J
C
             DO  130  j=ju,jl,4
               jn = nnn - (4*jjn)
               IF  (a)  GO TO 115
                 jx = (jjb-j)/4 + 1
                 WRITE(6,fmt2) jx, (kh(i,jn), i=1,mmm), jx
  115          CONTINUE
               jjn = jjn + 1
               IF  (jn.NE.1)  GO TO 118
C
C    SAVE LAST INUM BYTES OF I
C
                 DO  117  l=1,jjbb
                   DO  116  i=116,jjaapn
                     ia = i - 115
                     kk(ia,l,jnx) = kh(i,l)
  116              CONTINUE
  117            CONTINUE
                 a = .true.
                 GO TO 140
  118          CONTINUE
               DO  120  im=1,3
                 jn = jn - 1
                 print 1003, (kh(i,jn), i=1,mmm)
  120          CONTINUE
               a = .false.
  130        CONTINUE
  140      CONTINUE
           WRITE(6,fmt4) (line(n), n=1,la)
C
C    CALCULATE AND PRINT LAT/LON AT CORNERS
C
           al  = il(nx)
           ar  = ir(nx)
           xi1 = ((al-1.0)/xfac + 1.0) - 33.0
           xi2 = ((ar-1.0)/xfac + 1.0) - 33.0
           xj1 = cmjb - 33.0
           xj2 = cmjt - 33.0
           IF  (nref.EQ.0)  GO TO 142
             CALL w3fb05(xi1,xj1,-xmesh,260.0,alat1,alon1)
             CALL w3fb05(xi1,xj2,-xmesh,260.0,alat2,alon2)
             CALL w3fb05(xi2,xj2,-xmesh,260.0,alat3,alon3)
             CALL w3fb05(xi2,xj1,-xmesh,260.0,alat4,alon4)
             GO TO 144
  142      CONTINUE
             CALL w3fb05(xi1,xj1,xmesh,80.0,alat1,alon1)
             CALL w3fb05(xi1,xj2,xmesh,80.0,alat2,alon2)
             CALL w3fb05(xi2,xj2,xmesh,80.0,alat3,alon3)
             CALL w3fb05(xi2,xj1,xmesh,80.0,alat4,alon4)
  144      CONTINUE
           print 2001, alat2, alon2, alat3, alon3
           print 2002, alat1, alon1, alat4, alon4
           print 2003, nx, ipage, title
  150    CONTINUE
         ierr = 0
         RETURN
       END