C***********************************************************************
C *****                 PROGRAM    C H E C K                       ****
C **********************************************************************
C MODIFIED FROM ACCESS1 TO READ DAMEE OUTPUT FILE PROCESSED BY CVTMPP
C AND MAKE SIMPLE PLOT
C **********************************************************************
      program access

C ---------------------------------------------------------------------
C --- CHANGE PARAMETERS FOR EACH SUBDOMAIN
      
c-----------------------------------------------------------------------------
      parameter (idm=1437,jdm=1437,kdm=16,ms=30,msd=30,maxp=10000)! whole Atl.
      parameter (ipj=idm+jdm)
      parameter (nplots=6)

c-----------------------------------------------------------------------------
c --- information in common block gindex keeps do loops from running into land
c     common/gindex/ip(idm,jdm),iu(idm,jdm),iv(idm,jdm),iq(idm,jdm),
c    .ifp(jdm,ms),ilp(jdm,ms),isp(jdm),jfp(idm,ms),jlp(idm,ms),jsp(idm),
c    .ifq(jdm,ms),ilq(jdm,ms),isq(jdm),jfq(idm,ms),jlq(idm,ms),jsq(idm),
c    .ifu(jdm,ms),ilu(jdm,ms),isu(jdm),jfu(idm,ms),jlu(idm,ms),jsu(idm),
c    .ifv(jdm,ms),ilv(jdm,ms),isv(jdm),jfv(idm,ms),jlv(idm,ms),jsv(idm)
c
c     common/diags/ifd(idm+jdm,msd),jfd(idm+jdm,msd),ild(idm+jdm,msd),
c    . jld(idm+jdm,msd),nsec(idm+jdm)
c-----------------------------------------------------------------------------
c
c --- ms-1  = max. number of interruptions of any grid row or column by land
c --- msd-1 = max. number of interruptions of diagonal rows/columns by land
c --- maxp  = max. number of points outlining the ocean basin
c
      common/dimens/ii,jj,ii1,jj1
c-----------------------------------------------------------------------------
      common/fields/a(idm,jdm),b(idm,jdm)
      common/topog/zbot(idm,jdm)
      common/fields2/u(idm,jdm),v(idm,jdm),ubaro(idm,jdm),vbaro(idm,jdm)
     .               ,dp(idm,jdm)
      common/work/work(idm,jdm)
      dimension scalen(nplots),qlon(nplots),qhin(nplots),qqn(nplots)
      character*8 charvar(nplots)
      logical outln
      character util(idm*jdm+14)*2,flnm*60,text*8,preambl(5)*79
      character ch3*3,ch2*2,ch5*5,ch12*12
      character label*80,label2*40,label3*40
      real plot(jdm,idm),plot1(jdm,idm),plot2(jdm,idm)
c     logical lpltk(kdm)
c     logical lplt(nplots)
      integer lpltk(kdm),lplt(nplots),iptcon(nplots)
      common/geog/jref,reflon,equat,gridsz,lalolb
      xlat(i)=2.*atan(exp(gridsz*(equat-i)/57.29578))-1.5707963268
      data jref,reflon,equat,gridsz,lalolb/1,-98.,1071.,.08,10/
      
c --- 'thref' = reference specific volume (cm^3/g)
      data thref/1./,g/980.6/,onecm/.01/
      data charvar/
     .'u(cm/s)','v(cm/s)','dpm (m)','tem (C)','sal psu','velo  '/
      data scalen/  1.,   1.,98060.,  1.,  1., 1./
      data qlon/ -300.,-300.,    0., -2., 32., 0./
      data qhin/  300., 300., 5000., 32., 38.4, 0./
      data qqn/    50.,  50.,  100.,  2.,  .5, 0./
c 1= .true., 0=.false.
      data lplt/1,1,1,1,1,1/
      data lpltk/1,13*1,2*1/
      data iptcon/0,0,1,0,0,1/
c lpltk is modified in reading it below.
      data iptcol/1/


c --- parameters describing the geographic setting of the (total) model grid:
c --- 'reflon' = longitude of grid point 'jref'
c --- 'gridsz' = grid size in degrees longitude
c --- 'equat'  = i index giving location of equator in grid
CCC      data jref,reflon,equat,gridsz/1,-98.,1071.0,.08/
c --- latitude as a function of grid index
CCC      xlat(i)=2.*atan(exp(gridsz*(equat-i)/57.29578))-1.5707963268
c:w
c--- longitude as a function of grid index
CCC      xlon(j)=reflon+gridsz*(j-jref)
c --- name of file containing the actual data is read from   s t d i n.
      write(*,*)'enter year'
      read(*,*)iyr
      write(*,*)'enter day'
      read(*,*)iday
      write(ch2,'(i2.2)')iyr
      write(ch3,'(i3.3)')iday
      write(*,*)'enter iunc, 1 to uncompress'
      read(*,*)iunc
c     iunc=0
      flnm='0001,0001.02.222'
      flnm(11:12)=ch2
      flnm(14:16)=ch3
      write (*,'('' now opening data file '',a)') flnm
      ni=14
      if(iunc.eq.1)call system('gunzip '//flnm)
      open (unit=ni,file=flnm,
     .status='old',form='formatted')

c read topography 
      open(9,file=depthdamee.1437x1437',
     .     form='formatted',status='old')
      read (9,120) preambl,ii,jj,length,(util(l),l=1,length)
 120  FORMAT (5(a79/),3i8/(40a2))
      close(9)
      call unpakk(zbot,idm,ii,jj,util,length)
c     call bigrid(zbot)

c --- --------------------------------------------------------------------
c --- plot
c --- --------------------------------------------------------------------

c --- plot choosen fields for checking
      call opngks
      call gsclip(0)


c     call cpsetr ('spv - special value',-.03125)
      call cpsetr ('spv - special value',99999.)
      call cpseti ('pai - parameter array index', -2)
      call cpseti ('clu',1)
      call cpseti ('clc',0)

c --- read preamble identifying geographic extent of grid fragment

c --- data are packed into character strings, 2 bytes per grid point,
c --- preceded by header information:
c --- 'text'  -- type of field
c --- 'nstep' -- model time step
c --- 'time'  -- model time (days)
c --- 'thet'  -- vertical coordinate
c ---            (zero for 2-d fields, [1-(alpha/thref)] for 3-d fields)
c --- 'i'     -- number of points in i direction
c --- 'j'     -- number of points in j direction
c --- 'lgth'  -- length of 2-byte character string to follow.

c --- unpack 2-dimensional fields
 104  format (3x,a8,i8,2f8.0,2i4,i8/(1x,64a2))
c 104  format (3x,a8,i8,2f8.0,2i4,i6/(64a2))
 100  format (' ::',a8,i8,f8.1,f8.5,2i4,i8,1x,14a2)

c --- topography mask
      do i=1,ii
      do j=1,jj
      plot1(j,i)=-.100
      if(zbot(ii+1-i,j).lt.0.1) plot1(j,i) =100.0 
      enddo
      enddo

c --- read variables
c  Fields dumped once:
c     ubavg
c     vbavg
c     pbavg
c     montg
c     thmix
c
c  Fields dumped at every layer:
c  First k loop
c     u
c     v
c     dp
c     temp
c     saln
c  Second k loop
c     subduc

c --- read variables   
      read (ni,104,end=6) text,nstep,time,thet,ii,jj,
     .     lgth,(util(l),l=1,lgth)
      write(*,100) text(1:8),nstep,time,thet,ii,jj,lgth,(util(l),l=1,14)
      call unpakk(ubaro(1,1),idm,ii,jj,util,lgth)
      read (ni,104,end=6) text,nstep,time,thet,ii,jj,
     .     lgth,(util(l),l=1,lgth)
      write(*,100) text(1:8),nstep,time,thet,ii,jj,lgth,(util(l),l=1,14)
      call unpakk(vbaro(1,1),idm,ii,jj,util,lgth)
      do nvar=3,5
      read (ni,104,end=6) text,nstep,time,thet,ii,jj,
     .     lgth,(util(l),l=1,lgth)
      write(*,100) text(1:8),nstep,time,thet,ii,jj,lgth,(util(l),l=1,14)
      enddo

      do kk=1,kdm
      write(6,*)'reading k=',kk
      read (ni,104,end=6) text,nstep,time,thet,ii,jj,
     .     lgth,(util(l),l=1,lgth)
      write(*,100) text(1:8),nstep,time,thet,ii,jj,lgth,(util(l),l=1,14)
      if(lpltk(kk).eq.1.and.(lplt(1).eq.1.or.lplt(6).eq.1))then
      call unpakk(a(1,1),idm,ii,jj,util,lgth)
      do i=1,ii
       do j=1,jj
        u(j,i)=a(ii+1-i,j)+ubaro(ii+1-i,j)
        if(zbot(ii+1-i,j).lt.0.1)u(j,i)=99999.
       enddo
      enddo
      call minmax(u,jdm,idm,99999.,umin,umax)
      write(6,*)'min max u layer',kk,umin,umax
      endif
      read (ni,104,end=6) text,nstep,time,thet,ii,jj,
     .     lgth,(util(l),l=1,lgth)
      write(*,100) text(1:8),nstep,time,thet,ii,jj,lgth,(util(l),l=1,14)
      if(lpltk(kk).eq.1.and.(lplt(2).eq.1.or.lplt(6).eq.1))then
      call unpakk(a(1,1),idm,ii,jj,util,lgth)
      do i=1,ii
       do j=1,jj
        v(j,i)=a(ii+1-i,j)+vbaro(ii+1-i,j)
        if(zbot(ii+1-i,j).lt.0.1)v(j,i)=99999.
       enddo
      enddo
      call minmax(v,jdm,idm,99999.,vmin,vmax)
      write(6,*)'min max v layer',kk,vmin,vmax
      endif
      read (ni,104,end=6) text,nstep,time,thet,ii,jj,
     .     lgth,(util(l),l=1,lgth)
      write(*,100) text(1:8),nstep,time,thet,ii,jj,lgth,(util(l),l=1,14)
      call unpakk(dp(1,1),idm,ii,jj,util,lgth)
      if(lpltk(kk).eq.1.and.(lplt(1).eq.1.or.lplt(2).eq.1.or.
     .   lplt(6).eq.1))then
      do i=1,ii
       do j=1,jj
        if(dp(ii+1-i,j).lt..5e5)then
        u(j,i)=99999.
        v(j,i)=99999.
        endif
       enddo
      enddo
      endif

      nvar=3
      if(lpltk(kk).eq.1.and.lplt(3).eq.1)then
      do i=1,ii
       do j=1,jj
        plot(j,i)=dp(ii+1-i,j)/scalen(nvar)
        if(zbot(ii+1-i,j).lt.0.1)plot(j,i)=99999.
c       if(dp(ii+1-i,j).lt..5e5)plot(j,i)=99999.
       enddo
      enddo
      call minmax(plot,jdm,idm,99999.,amin,amax)
      write(6,*)'min max ',charvar(nvar),' layer',kk,amin,amax
      nquad=-1        !  next plot full-frame
      offset=0.5
      outln=.true.
      ldash=0
      label=
     .'micom 0.08d  ???      lay         year='//ch2//' ,day='//ch3
      write(label(26:27),'(i2)')kk
      write(label(14:21),'(a8)')charvar(nvar)
      qlo=0.
      qhi=qhin(nvar)
      qq=qqn(nvar)
      if(amax.le.(qhi/10))qq=qq/2.
c     ldash=-366
      label2='min       max       contour int         '
      write(label2(5:9),'(f5.1)')amin
      write(label2(15:20),'(f6.1)')amax
      write(label2(36:40),'(f5.1)')qq
      if(iptcon(nvar).eq.1)then
      call hormp(plot1,jdm,idm,ii,jj,-.1,.1,.1,                 
     .                  ldash,label,label2,label3,offset,nquad,outln)
      write(6,*)'finished topography plot'
      call horplt(plot,jdm,idm,ii,jj,qlo,qhi,qq,
     .                  ldash,label,label2,label3,offset,nquad,outln,1)
      write(6,*)'finished dp plot layer',kk
      endif
      if(iptcol.eq.1)then
      ch5(1:3)=charvar(nvar)(1:3)
      write(ch5(4:5),'(i2.2)')kk
      ch12=ch5//flnm(10:16)
      do i=1,ii
       do j=1,jj
        plot2(j,ii+1-i)=plot(j,i)
       enddo
      enddo
      write(6,*)'doing ppm image'//ch5, 'layer', kk, 'qlo, qhi',qlo,qhi
      call imageppm(plot2,qlo,qhi,99999.,ch12)
      endif
      write(6,*)'finished variable plot'//ch5,kk

      endif

      do nvar=4,5
      read (ni,104,end=6) text,nstep,time,thet,ii,jj,
     .     lgth,(util(l),l=1,lgth)
      write(*,100) text(1:8),nstep,time,thet,ii,jj,lgth,(util(l),l=1,14)
      if(lpltk(kk).eq.1.and.(lplt(nvar).eq.1))then
      call unpakk(a(1,1),idm,ii,jj,util,lgth)
      do i=1,ii
       do j=1,jj
        plot(j,i)=a(ii+1-i,j)/scalen(nvar)
        if(zbot(ii+1-i,j).lt.0.1)plot(j,i)=99999.
        if(dp(ii+1-i,j).lt..5e5)plot(j,i)=99999.
       enddo
      enddo
      call minmax(plot,jdm,idm,99999.,amin,amax)
      write(6,*)'min max ',charvar(nvar),' layer',kk,amin,amax
      nquad=-1        !  next plot full-frame
      offset=0.5
      outln=.true.
      ldash=0
      label=
     .'micom 0.08d  ???      lay         year='//ch2//' ,day='//ch3
      write(label(26:27),'(i2)')kk
      write(label(14:21),'(a8)')charvar(nvar)
      qlo=qlon(nvar)
      qhi=qhin(nvar)
      qq=qqn(nvar)
c     ldash=-366
      label2='min       max       contour interval   '
      write(label2(5:9),'(f5.1)')amin
      write(label2(15:19),'(f5.1)')amax
      write(label2(38:40),'(f3.1)')qq
      if(iptcon(nvar).eq.1)then
      call hormp(plot1,jdm,idm,ii,jj,-.1,.1,.1,
     .                  ldash,label,label2,label3,offset,nquad,outln)
      write(6,*)'finished topography plot'
      call horplt(plot,jdm,idm,ii,jj,qlo,qhi,qq,
     .                  ldash,label,label2,label3,offset,nquad,outln,1)
      endif
c temperature
c     data bmin,bmax/0.,32./
c salinity
c     data bmin,bmax/32.,38.4/
c ssh
c     data bmin,bmax/-1.1,.72/
      ch5(1:3)=charvar(nvar)(1:3)
      write(ch5(4:5),'(i2.2)')kk
      ch12=ch5//flnm(10:16)
      if(iptcol.eq.1)then
      do i=1,ii
       do j=1,jj
        plot2(j,ii+1-i)=plot(j,i)
       enddo
      enddo
      write(6,*)'doing ppm image'//ch5, 'layer', kk, 'qlo, qhi',qlo,qhi
      call imageppm(plot2,qlo,qhi,99999.,ch12)
      endif
      write(6,*)'finished variable plot'//ch5,kk
      endif
      enddo

      if(lpltk(kk).eq.1.and.lplt(6).eq.1)then
      label='micom 0.08d  veloc   lay         year='//ch2//' ,day='//ch3
      write(label(26:27),'(i2)')kk
      write(6,'(a)')label
      nquad=-1        !  next plot full-frame
      offset=0.5
      outln=.true.
      ldash=0
c plot uv vector plot
     
      label2='                                 '                                      
      label3='                                 '
      call hormp(plot1,jdm,idm,ii,jj,-.1,.1,.1,
     .                  ldash,label,label2,label3,offset,nquad,outln)
      write(6,*)'ii,jj',ii,jj
      siz=.45
      xmid=.5
      ymid=.5
      call set(xmid-siz*amin1(1.,float(jj-1)/float(ii-1))
     .        ,xmid+siz*amin1(1.,float(jj-1)/float(ii-1))
     .        ,ymid-siz*amin1(1.,float(ii-1)/float(jj-1))
     .        ,ymid+siz*amin1(1.,float(ii-1)/float(jj-1))
     .        ,1.-offset,float(jj)-offset,1.-offset,float(ii)-offset,1)
      if(kk.le.14)q=.3
      if(kk.ge.15)q=.6
      cutoffv=0.
      cutoff=q*2.*cutoffv
      do i=2,ii,5
      do j=1,jj-1,5
      if(u(j,i).eq.99999..or.u(j,i-1).eq.99999..
     .   or.v(j,i).eq.99999.or.v(j+1,i).eq.99999.)go to 10
      qu=.5*q*(u(j,i)+u(j,i-1))
      qv=.5*q*(v(j,i)+v(j+1,i))
      call arrow(float(j)-qv,float(i)+qu
     .          ,float(j)+qv,float(i)-qu,cutoff)

10    continue
      enddo
      enddo
      qu=0.
      qv=q*50.
      xa=float(200)-qv 
      ya=float(200)+qu
      xb=float(200)+qv 
      yb=float(200)-qu
      call arrow(xa,ya,xb,yb,cutoff)
      call wtstr(xb+10.,ya,'50 cm/s',0,0,-1)

c     CALL FRSTPT(Xb,Yb)
c     q=1./max(.1,sqrt((xb-xa)**2+(yb-ya)**2))
c     CALL VECTOR(XB-(XB-XA-(YB-YA)*0.4)*Q,YB-(YB-YA+(XB-XA)*0.4)*Q)
c     CALL VECTOR(XB-(XB-XA+(YB-YA)*0.4)*Q,YB-(YB-YA-(XB-XA)*0.4)*Q)
      label2='                                 '                                      
      label3='                                 '
      call horplt(plot1,jdm,idm,ii,jj,qlo,qhi,qq,
     .                  ldash,label,label2,label3,offset,nquad,outln,0)
      endif
      enddo                                              !kk 
6     continue


      call clsgks
      call system ('mv gmeta gmeta.'//ch2//'.'//ch3)
      stop '(normal)'
      end

      SUBROUTINE ARROW(XAA,YAA,XBB,YBB,cutoff)
C
C --- DRAW ARROW POINTING FROM (XAA,YAA) TO (XBB,YBB)
C --- IF ARROW IS LONGER THAN 'CUTOFF', do not plot it.
C
c     DATA CUTOFF/2./,BARB/.2/,HEADSZ/2./
c     REDUC=.5*CUTOFF/AMAX1(CUTOFF,SQRT((XBB-XAA)**2+(YBB-YAA)**2))
      data headsz/1./
      reduc=.5
      XA=(.5+REDUC)*XAA+(.5-REDUC)*XBB
      XB=(.5+REDUC)*XBB+(.5-REDUC)*XAA
      YA=(.5+REDUC)*YAA+(.5-REDUC)*YBB
      YB=(.5+REDUC)*YBB+(.5-REDUC)*YAA
 2    Q=SQRT((XB-XA)**2+(YB-YA)**2)
      if(q.le.cutoff)return
      CALL FRSTPT(XA,YA)
      CALL VECTOR(XB,YB)
      Q=HEADSZ/Q
      CALL VECTOR(XB-(XB-XA-(YB-YA)*0.4)*Q,YB-(YB-YA+(XB-XA)*0.4)*Q)
      CALL VECTOR(XB-(XB-XA+(YB-YA)*0.4)*Q,YB-(YB-YA-(XB-XA)*0.4)*Q)
      CALL VECTOR(XB,YB)
      RETURN
      END


      SUBROUTINE PAKK(ARRAY,IDIM,II,JJ,COMPAC,LENGTH)
C
C --- Converts the contents of -ARRAY- into an ASCII character string which
C --- is stored in character*2 array -COMPAC-. COMPAC(1)...COMPAC(7) contain
C --- the base value, i.e., the minimum value encountered in -ARRAY-.
C --- COMPAC(8)...COMPAC(14) contain a scale factor by which the individual
C --- 6-bit integers encoded as ASCII character pairs in COMPAC(8),...
C --- COMPAC(LENGTH) must be multiplied before the base value is added
C --- during an UNPAKKing operation. Base value and scale factor are encoded
C --- in E14.7 format.
C
C --- The printable ASCII characters used to encode the integers include
C --- the numbers 0...9, upper- and lower-case letters A...Z, a...z, plus
C --- two additional characters '.' and '/' (total of 64).
C
C --- A packing operation fills (II*JJ+14) array elements in -COMPAC- which
C --- must be dimensioned accordingly in the calling program. The total
C --- number of occupied array elements is returned in -LENGTH-. In calls to
C --- UNPACK, -LENGTH- is treated as input variable.
C
      real array(idim,1)
      character*2 char,compac(1),comp2(14)
      character*14 comp14(2)
      equivalence (comp2,comp14)
      data nbits/12/
      base=1.e22
      do 1 i=1,ii
      do 1 j=1,jj
 1    base=amin1(base,array(i,j))
      scal=0.
      do 2 i=1,ii
      do 2 j=1,jj
 2    scal=amax1(scal,array(i,j)-base)
      scal=scal/float(2**nbits-1)
      i1=0
      i2=0
      length=14
      do 3 i=1,ii
      do 3 j=1,jj
      if (scal.eq.0.) go to 7
      numb=(array(i,j)-base)/scal+.5
      i1=numb/64
      i2=numb-64*i1
c
c --- Map 6-bit numbers onto character set consisting of numbers
c --- 0...9, letters A...Z, a...z, and the two characters '.' and '/'.
c --- (If mapping into the character range 32...95 -- which includes the
c --- characters !"#$%&'()*+,-./:;<=>?@[\]^_  -- is deemed safe, delete
c --- the next 6 lines.)
c      if (i1.gt.37) i1=i1+6
c      if (i1.gt.11) i1=i1+7
c      i1=i1+14
c      if (i2.gt.37) i2=i2+6
c      if (i2.gt.11) i2=i2+7
c      i2=i2+14
c
 7    length=length+1
      compac(length)(1:1)=char(i1+32)
      compac(length)(2:2)=char(i2+32)
 100  format (a2)
 3    continue
      write (comp14(1),101) base
      write (comp14(2),101) scal
 101  format (1pe14.7)
      do 8 l=1,14
 8    compac(l)=comp2(l)
c
      return
c
c
      ENTRY UNPAKK(ARRAY,IDIM,II,JJ,COMPAC,LENGTH)
c
      do 9 l=1,14
 9    comp2(l)=compac(l)
      read (comp14(1),101) base
      read (comp14(2),101) scal
      lngth=14
      do 4 i=1,ii
      do 4 j=1,jj
      lngth=lngth+1
      i1=ichar(compac(lngth)(1:1))
      i2=ichar(compac(lngth)(2:2))
c
c --- 6-bit numbers are mapped onto character set consisting of numbers
c --- 0...9, letters A...Z, a...z, and the two characters '.' and '/'.
c --- (If mapped into character range 32...95, delete next 6 lines)
c      if (i1.gt.96) i1=i1-6
c      if (i1.gt.64) i1=i1-7
c      i1=i1-14
c      if (i2.gt.96) i2=i2-6
c      if (i2.gt.64) i2=i2-7
c      i2=i2-14
c
 4    array(i,j)=scal*float(64*(i1-32)+(i2-32))+base
      if (lngth.ne.length) stop 'UNPACK'
      return
      end

      subroutine horplt(rotat,idm,jdm,ii1,jj1,qlo,qhi,qq,
     .                  ldash,label,label2,label3,offset,nquad,outln,
     .                  icont)
      character label*80,label2*40,label3*40,lab*4
      character xlab1*1
      dimension rotat(jdm,idm)
      common/geog/jref,reflon,equat,gridsz,lalolb
      logical outln,paint
      common /area1/ paint
c
      write(6,*)'horplt',idm,jdm,ii1,jj1,qlo,qhi,qq,ldash,
     .label,offset,nquad,outln	
      ii=min(ii1+1,idm)
      jj=min(jj1+1,jdm)
      write(6,*)'ii,jj',ii,jj
c
      if (outln) write (*,'('' horiz. map: ''a40)') label
      paint=outln
      if (nquad.lt.0) then
c     siz=.47
      siz=.45
      xmid=.5
      ymid=.5
      else
      siz=.22
      xmid=.5+.25*float(2*mod(nquad,2)-1)
      ymid=.49+.25*(1.-2.*float(nquad/2))
      end if
c     nquad=mod(nquad+1,4)
c
c --- rotate field 90 degrees
c
      call getset(pl,pr,pb,pt,hl,hr,hb,ht,ll)
      call set(xmid-siz*amin1(1.,float(jj-1)/float(ii-1))
     .        ,xmid+siz*amin1(1.,float(jj-1)/float(ii-1))
     .        ,ymid-siz*amin1(1.,float(ii-1)/float(jj-1))
     .        ,ymid+siz*amin1(1.,float(ii-1)/float(jj-1))
     .        ,1.-offset,float(jj)-offset,1.-offset,float(ii)-offset,1)
      write(6,*)'set',xmid-siz*amin1(1.,float(jj-1)/float(ii-1)),
     .xmid+siz*amin1(1.,float(jj-1)/float(ii-1)),
     .ymid-siz*amin1(1.,float(ii-1)/float(jj-1)),
     .ymid+siz*amin1(1.,float(ii-1)/float(jj-1)),
     .1.-offset,float(jj)-offset,1.-offset,float(ii)-offset,1
      if(icont.eq.1)then
      call conrec(rotat,jdm,jj1,ii1,qlo,qhi,qq,1,-1,ldash,0.)
      if (.999*qq.lt.qlo.and.qq.gt..999*qhi) go to 4
      endif
c
c --- print title
      chrsiz=.010*float(jj-1)/
     .       (2.*siz*amin1(1.,float(jj-1)/float(ii-1)))
      if (outln) call wtstr(float(jj/2),float(ii)-offset+1.6*chrsiz,
     .   label,1,0,0)
      if (outln) call wtstr(float(jj/2),float(-15)+offset-1.6*chrsiz,
     .   label2,1,0,0)
c
c --- add latitude labels
      write(xlab1,'(a1)')'-'
      do 1 lat=-30,70,lalolb
      y=float(ii)-xoflat(lat)
      if (y.le.1..or.y.ge.float(ii)) go to 1
      write (lab,'(i2)') iabs(lat)
      if (xmid.lt..5) then
      call wtstr(       1.-offset-chrsiz,y-offset,lab(1:2),0,0,+1)
      else
      call wtstr(float(jj)-offset+chrsiz,y-offset,lab(1:2),0,0,-1)
      call wtstr(1.-.1*chrsiz,y-offset,xlab1,0,0,-1)
      call wtstr(float(jj)+.1*chrsiz,y-offset,xlab1,0,0,+1)

      endif
 1    continue
      write(6,*)'end 1'
c
c --- add longitude labels
      do 3 lon=-100,10,lalolb
      x=(float(lon)-reflon)/gridsz+jref
      if (x.le.1..or.x.ge.float(jj)) go to 3
      write (lab,'(i4)') lon
      if      (lab(3:3).eq.' ') then
      call wtstr(x-offset,1.-offset-chrsiz,lab(4:4),0,0,0)
      else if (lab(2:2).eq.' ') then
      call wtstr(x-offset,1.-offset-chrsiz,lab(3:4),0,0,0)
      else if (lab(1:1).eq.' ') then
      call wtstr(x-offset,1.-offset-chrsiz,lab(2:4),0,0,0)
      else
      call wtstr(x-offset,1.-offset-chrsiz,lab(1:4),0,0,0)
      endif
      call wtstr(x-offset,1.,xlab1,0,90,-1)
      call wtstr(x-offset,float(ii),xlab1,0,90,1)
 3    continue
      write(6,*)'end 3'
      if (offset.ne.0.) then
       call set(xmid-siz*amin1(1.,float(jj-1)/float(ii-1))
     .         ,xmid+siz*amin1(1.,float(jj-1)/float(ii-1))
     .         ,ymid-siz*amin1(1.,float(ii-1)/float(jj-1))
     .         ,ymid+siz*amin1(1.,float(ii-1)/float(jj-1))
     .         ,1.,float(jj),1.,float(ii),1)
c      if (outln) call bordplt
       call set(xmid-siz*amin1(1.,float(jj-1)/float(ii-1))
     .         ,xmid+siz*amin1(1.,float(jj-1)/float(ii-1))
     .         ,ymid-siz*amin1(1.,float(ii-1)/float(jj-1))
     .         ,ymid+siz*amin1(1.,float(ii-1)/float(jj-1))
     .         ,1.-offset,float(jj)-offset,1.-offset,float(ii)-offset,1)
c      write(6,*)'end set bordplt'
      else				!  offset = 0
c      if (outln) call bordplt
c      write(6,*)'end bordplt'
      end if
      call perim(1,1,1,1)
      call set(pl,pr,pb,pt,hl,hr,hb,ht,ll)
      call frame

c
 4    return
      end


      subroutine hormp(rotat,idm,jdm,ii1,jj1,qlo,qhi,qq,
     .                  ldash,label,label2,label3,offset,nquad,outln)
      character label*80,label2*40,label3*40,lab*4
      character xlab1*1
      dimension rotat(jdm,idm)
      common/geog/jref,reflon,equat,gridsz,lalolb
      logical outln,paint
      common /area1/ paint
c
      ii=min(ii1+1,idm)
      jj=min(jj1+1,jdm)
c
      if (outln) write (*,'('' horiz. map: ''a40)') label
      paint=outln
      if (nquad.lt.0) then
c     siz=.47
      siz=.45
      xmid=.5
      ymid=.5
      else
      siz=.22
      xmid=.5+.25*float(2*mod(nquad,2)-1)
      ymid=.49+.25*(1.-2.*float(nquad/2))
      end if
c
c --- rotate field 90 degrees
c
      call cpgetr('CLL - CONTOUR-LINE LINE WIDTH',rline)
      write(6,*)'cll',rline
      call cpsetr('CLL - CONTOUR-LINE LINE WIDTH',2.)
      call getset(pl,pr,pb,pt,hl,hr,hb,ht,ll)
      call set(xmid-siz*amin1(1.,float(jj-1)/float(ii-1))
     .        ,xmid+siz*amin1(1.,float(jj-1)/float(ii-1))
     .        ,ymid-siz*amin1(1.,float(ii-1)/float(jj-1))
     .        ,ymid+siz*amin1(1.,float(ii-1)/float(jj-1))
     .        ,1.-offset,float(jj)-offset,1.-offset,float(ii)-offset,1)
      write(6,*)'set',xmid-siz*amin1(1.,float(jj-1)/float(ii-1)),
     .xmid+siz*amin1(1.,float(jj-1)/float(ii-1)),
     .ymid-siz*amin1(1.,float(ii-1)/float(jj-1)),
     .ymid+siz*amin1(1.,float(ii-1)/float(jj-1)),
     .1.-offset,float(jj)-offset,1.-offset,float(ii)-offset,1
      call conrec(rotat,jdm,jj1,ii1,qlo,qhi,qq,1,-1,ldash,2.)
      if (.999*qq.lt.qlo.and.qq.gt..999*qhi) go to 4
c
      call set(pl,pr,pb,pt,hl,hr,hb,ht,ll)
      call cpsetr('CLL - CONTOUR-LINE LINE WIDTH',0.)
 4    return
      end
c =====================================================================
      subroutine conrec(z,k,m,n,flo,fhi,finc,nset,nhilo,ndash,width)
c
c --- redirect   c o n r e c   calls  to   c o n p a c k
c
c --- declarations required for area fill:
c      parameter(lgthmp=1500000,lgthwk=20000)
      parameter(lgthmp=6000000,lgthwk=80000)
      common /area/ map(lgthmp),xwork(lgthwk),ywork(lgthwk)
      common /colopt/ ipalet,ibase(10)
      dimension iarea(64),igrup(64)
      logical label
      external shade,cpdrpl
c
      parameter (lgthr=30000,lgthi=50000)
      dimension rwork(lgthr),iwork(lgthi),z(k,1)
c     data spval/-.03125/
      data spval/99999./
c
      write (*,*) 'call to conrec => conpack'
c --- decide whether to label isolines
      label=.true.
      if (ipalet.gt.1) label=.false.
      write(6,*)'label',label
c
c --- activate special value feature
      call cpsetr('spv - special value',spval)
      call cpseti('pai - parameter array index',-2)
      call cpseti('clu - contour level use directive',0)
c
      if (nset.gt.0)
     .call cpseti('set - set call directive',0)
c
      call cpseti('cls - contour level selection',50)
      call cpsetr('cis - contour interval specifier',finc)
      call cpseti('rwc - real workspace for contours',500)
      if (label) then
        call cpseti('lis - label interval specifier',3)
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
        call cpgetr('CWM - character-width multiplier',cwm)
        CALL CPGETR('LLW - line label white space',rllw)
        call cpgetr('CLL - contour line width',rclls)
        write(6,*)' cwm llw cll',cwm,rllw,rclls
c --- switch from "default" to "regular" line labeling scheme
        call cpseti('llp - line label positioning',2)
        call cpsetr('rc1 - regular scheme constant 1',.1)
        call cpsetr('hlw - hi/lo label white space',.002)
        call cpsetr('llw - line label white space',.002)
        call cpsetr('cwm - character width multiplier',.9)
        call cpseti('llo - line label orientation',1)
        call cpsetc('ilt - information label text string',' ')
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      end if
c
      if (flo.eq.0. and. fhi.eq.0.) then
c --- set minimum > maximum to allow unconstrained contouring
        call cpsetr('cmn - contour minimum',1.)
        call cpsetr('cmx - contour maximum',0.)
      else
c --- pass on prescribed min/max values
        call cpsetr('cmn - contour minimum',flo)
        call cpsetr('cmx - contour maximum',fhi)
      end if
c
      call cprect(z,k,m,n,rwork,lgthr,iwork,lgthi)
      call cppkcl(z,rwork,iwork)
      call cpgeti('ncl - number of contour levels',num)
      if (num.le.0) return
c set line width for continent boundary
      do i=1,num
      call cpseti('pai - parameter array index',i)
      if(width.gt.0.)call cpsetr('CLL - CONTOUR-LINE LINE WIDTH',width)
      enddo
c
c --- set dash pattern (convert from 10 to 16 bit pattern)
      ndsh10=iabs(ndash)
      ndsh16=64*ndsh10+mod(ndsh10,64)
      if (ndash.ne.0) then
        do 1 i=1,num
          call cpseti('pai - parameter array index',i)
          if (ndash.gt.0) then
c --- dash all contour lines
            call cpseti('cld - contour line dash pattern',ndsh16)
          else
c --- dash negative contours only
            call cpgetr('clv - contour level value',valu)
            if (valu.lt.0.)
     .      call cpseti('cld - contour line dash pattern',ndsh16)
          end if
 1      continue
      end if
c
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      if (ipalet.gt.0) then		! add color to contour plots
      write (*,100) ibase(ipalet)+1,ibase(ipalet+1)
 100  format (' use color table entries',i4,' --',i4)
c
c --- initialize area fill routine
      call arinam(map,lgthmp)
      do 2 i=1,num
        call cpseti('pai - parameter array index',i)
        call cpseti('aia - area indicator above',0)
        call cpseti('aib - area indicator below',0)
 2    continue
c
      if (ipalet.eq.1) then
c --- find zero isoline
        valu1=0.
        izero=-1
        do 4 i=1,num
          valu2=valu1
          call cpseti('pai - parameter array index',i)
          call cpgetr('clv - contour level value',valu1)
          if (valu1.eq.0. .or. valu1.eq.-valu2) izero=i
 4      continue
      end if				!  ipalet = 1
c
      do 3 i=1,num
        call cpseti('pai - parameter array index',i)
        call cpgetr('clv - contour level value',valu1)
c
        if (ipalet.eq.1) then
c --- paint every 2nd gap between contour lines to create 'zebra' effect
          if (mod(i-izero+99,2).eq.0) go to 3
c --- use different colors for positive and negative areas
          indx=1
          if (valu1.lt.0.) indx=2
c
        else				!  ipalet > 1
c --- paint each gap with a different color
          indx=i
        end if				!  ipalet = or > 1
        indx=min(indx+ibase(ipalet),ibase(ipalet+1))
c
        call cpseti('aia - area indicator above',indx)
        write (*,102) valu1,indx
 102  format (' set area index above contour line ',1pe9.2,' to',i4)
        if (i.lt.num) then
          call cpseti('pai - parameter array index',i+1)
          call cpseti('aib - area indicator below',indx)
          call cpgetr('clv - contour level value',valu2)
          write (*,103) valu2,indx
 103  format (' set area index below contour line ',1pe9.2,' to',i4)
        end if
 3    continue
c
      call cpclam(z,rwork,iwork,map)
      if (label) call cplbam(z,rwork,iwork,map)
      call arscam(map,xwork,ywork,lgthwk,iarea,igrup,64,shade)
c
      end if				!  ipalet > 0
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c
      if (ipalet.le.0) then
        call cpcldr(z,rwork,iwork)
      else if (label) then
        call cpcldm(z,rwork,iwork,map,cpdrpl)
        call cplbdr(z,rwork,iwork)
      end if
c
      return
      end
c
c
      subroutine shade(xcs,ycs,ncs,iarea,igrup,ngrups)
c
c --- this routine does the actual coloring of contour line intervals
      dimension xcs(ncs),ycs(ncs),iarea(ngrups),igrup(ngrups)
ccc      write (*,'('' shade: group/area='',8i3))')
ccc     .   (igrup(n),iarea(n),n=1,ngrups)
      do 1 n=1,ngrups
c --- conpack defines areas formed by contour lines to be in group 3
      if (igrup(n).eq.3.and.iarea(n).gt.0) then
        call gsfaci(iarea(n))
        call gqcr(1,iarea(n),1,ierr,r,g,b)
ccc        write (*,100) 'colors set to',r,g,b,'(color index',iarea(n),
ccc     .   ')   ierr =',ierr
 100    format (a,3f6.3,5x,a,i4,a,i3)
        call gfa(ncs,xcs,ycs)
      end if
 1    continue
      return
      end

    
      function xoflat(lat)
c --- x as a function of latitude
      common/geog/jref,reflon,equat,gridsz
      data radian/57.2957795/
      xoflat=equat-radian*alog(tan(.5*(lat+90.)/radian))/gridsz
      return
      end

      subroutine minmax(a,imt,jmt,bad,amin,amax)
      dimension a(imt,jmt)
      amin=1.e30
      amax=-1.e30
      do i=1,imt
      do j=1,jmt
      if(a(i,j).lt.amin.and.a(i,j).ne.bad)then
      amin=a(i,j)
      imin=i
      jmin=j
      endif
      if(a(i,j).gt.amax.and.a(i,j).ne.bad)then
      amax=a(i,j)
      imax=i
      jmax=j
      endif
      enddo
      enddo
      write(6,*)' min ', amin, imin, jmin, ' max ', amax, imax, jmax
      return
      end

      subroutine imageppm(plot,qlo,qhi,bad,ch12)
      parameter (imt=1437,jmt=1437)
c note:change parameters in subroutines
      dimension plot(imt,jmt)
      character*80 fnam, cmmd
      character ch12*12
      common /itim/ifirst
      ifirst=1
      open(7,file='imageppm.'//ch12,form='formatted',status='unknown')
      call writch(plot,bad,qlo,qhi)
      close (7)
      cmmd='/usr/rcf/bin/ppmtogif < imageppm.'//ch12//' >'//ch12//'.gif'
      call system (cmmd)
      call system ('rm -f imageppm.'//ch12)
      return
      end 

      subroutine palette(b,ired,igreen,iblue,qlo,qhi)
      parameter (ncounts=256)
      common /itim/ifirst
      dimension ir(ncounts),ig(ncounts),ib(ncounts)
c temperature
c     data bmin,bmax/0.,32./
c salinity
c     data bmin,bmax/32.,38.4/
c ssh
c     data bmin,bmax/-1.1,.72/
      bmin=qlo
      bmax=qhi
      if(ifirst.eq.1)then
      open(1,file='palette.micom.dp',
     .form='formatted',status='old')
      do n=1,ncounts
      read(1,*)ir(n),ig(n),ib(n)
      enddo
      close(1)
      endif
      ifirst=0
      brange=bmax-bmin
      if(b.lt.bmin.or.b.gt.bmax)then
      ired=0
      igreen=0
      iblue=0
      return
      else
      icnt=(b-bmin)*256/brange
      ired=ir(icnt)*256/32
      igreen=ig(icnt)*256/32
      iblue=ib(icnt)*256/32
      endif
      return
      end

      subroutine writch(a,bad,qlo,qhi)
      parameter(imt=1437,jmt=1437)
      dimension a(imt,jmt)
      common/rbg/ired(imt,jmt),iblue(imt,jmt),igreen(imt,jmt)
      common /itim/ifirst
      do j=1,jmt
      do i=1,imt
      ired(i,j)=0
      igreen(i,j)=0
      iblue(i,j)=0
      if(a(i,j).ne.bad)
     .call palette(a(i,j),ired(i,j),igreen(i,j),iblue(i,j),qlo,qhi)
      enddo
      enddo
      call paths
c     call isobath
      write(7,'(a)')'P6'
      write(7,'(i4,1x,i4)')imt,jmt
      write(7,'(i3)')255
      write(7,'(6194907a1)')
     *((char(ired(i,j)),char(igreen(i,j)),char(iblue(i,j)),
     *i=1,imt),j=1,jmt)
      return
      end

      subroutine paths
      parameter (idm=1437,jdm=1437,imt=jdm,jmt=idm)
      common/grid/xlat(idm),xlong(jdm)
      common/rbg/ired(imt,jmt),iblue(imt,jmt),igreen(imt,jmt)
      common/topog/zbot(idm,jdm)
      dimension pathin(2,400,3),pathnin(2,400,3),pathsin(2,400,3)
      call lats
      open(2,file='path.mean'
     . ,form='formatted',status='old')
      open(3,file='path.north'
     . ,form='formatted',status='old')
      open(4,file='path.south'
     . ,form='formatted',status='old')

      do m=1,3
      ni=m+1
      read(ni,*)
      do n=1,400
      read(ni,*)pathin(1,n,m),pathin(2,n,m)
      enddo
      close(ni)
      enddo

      do m=1,3
      write(6,*)'starting long,lat',m,pathin(1,1,m),pathin(2,1,m)
      do n=1,400
      do i=1,idm
      if(pathin(1,n,m).gt.-75.)then
      if(xlat(i).gt.pathin(2,n,m).and.xlat(i+1).le.pathin(2,n,m))then
      j1=(pathin(1,n,m)-xlong(1))/(xlong(2)-xlong(1))+1 
      do i1=-1,1
      ired(j1,i+i1)=255
      iblue(j1,i+i1)=255
      igreen(j1,i+i1)=255
      enddo
      endif
      endif
      enddo
      enddo
      enddo
c plot 200 m isobath
      zbot1=200.
      niso=0
      do 50 i=1,idm
      do 50 j=1,jdm

      if(i.lt.idm)then
      if(zbot(i,j).eq.zbot1.or.
     .   (zbot(i,j).le.zbot1.and.zbot(i+1,j).ge.zbot1).or.
     .   (zbot(i,j).ge.zbot1.and.zbot(i+1,j).le.zbot1))then
      do 30 i1=0,1
      ired(j,i+i1)=125
      iblue(j,i+i1)=125
      igreen(j,i+i1)=125
30    continue
      endif
      endif
      if(j.lt.jdm)then
      if(zbot(i,j).eq.zbot1.or.
     .   (zbot(i,j).le.zbot1.and.zbot(i,j+1).ge.zbot1).or.
     .   (zbot(i,j).ge.zbot1.and.zbot(i,j+1).le.zbot1))then
      ired(j,i+1)=125
      iblue(j,i+1)=125
      igreen(j,i+1)=125
      endif
      endif

50    continue
      
      return
      end
      
      subroutine lats
      parameter (idm=1437,jdm=1437)
      common/grid/xlat(idm),xlong(jdm)
      data rad/57.29578/,pi/3.14159265/,grid/.08/,xpiv/1071./
      data rlong/-98.0/
      ALAT(DIST,GRIDD)=2.*ATAN(EXP(DIST*GRIDD/RAD))-PI/2.
c --- latitudes
      DO 4 I=1,IDM
      xlat(i)=alat(xpiv-i,grid)*rad
4     continue
      x1=alat(0.,grid)*rad
      if(abs(x1).gt.0.001)then
      write(6,*)'xlat at xpiv should be zero',
     .        x1,xpiv
      stop
      endif

c --- longitudes
      DO 5 j=1,jdm
      xlong(j)=rlong+(j-1)*grid
5     CONTINUE
      return
      end 
