C***********************************************************************
C *****                 PROGRAM    C H E C K                       ****
C **********************************************************************
      program access
c-----------------------------------------------------------------------------
      parameter (idm=1437,jdm=1437,kdm=16,ms=30,msd=30,maxp=10000)! whole Atl.
      parameter (ipj=idm+jdm)
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,kdm),v(idm,jdm,kdm),ubaro(idm,jdm)
     .              ,vbaro(idm,jdm),dp(idm,jdm,kdm)
     .              ,saln(idm,jdm,kdm),temp(idm,jdm,kdm)  
      common/work/work(idm,jdm)
      logical outln
      character util(idm*jdm+14)*2,flnm*60,text*8,preambl(5)*79
      character ch3*3,ch2*2,ch5*5,ch12*12
      common/geog/jref,reflon,equat,gridsz,lalolb
c --- 'thref' = reference specific volume (cm^3/g)
      data thref/1./,g/980.6/,onecm/.01/
c
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
      data jref,reflon,equat,gridsz/1,-98.,1071.0,.08/
c
c --- latitude as a function of grid index
      xlat(i)=2.*atan(exp(gridsz*(equat-i)/57.29578))-1.5707963268
c
c--- longitude as a function of grid index
      xlon(j)=reflon+gridsz*(j-jref)
c
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
c     iunc=0
      flnm='0001,0001.02.222'
      flnm(11:12)=ch2
      flnm(14:16)=ch3
      write (*,'('' now opening data file '',a)') flnm
c
      ni=14
      open (unit=ni,file=flnm,
     .status='old',form='formatted')
c
c read topography file
c
      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
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
c --- read variables   
c
c  Fields dumped once:
c     ubavg
c     vbavg
c     pbavg
c     montg
c     thmix
c
c  read ubaro and unpack
      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)
c  read vbaro and unpack
      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)
c  read additional fields not needed for visualization
      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
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
      do k=1,kdm                     ! number of layers
      write(6,*)'reading k=',k
c read ubaclin
      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(a(1,1),idm,ii,jj,util,lgth)
c create total velocity u
      do i=1,ii
       do j=1,jj
        u(i,j,k)=a(i,j)+ubaro(i,j)
        if(zbot(i,j).lt.0.1) u(i,j,k)=99999.
       enddo
      enddo
c read vbaclin
      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(a(1,1),idm,ii,jj,util,lgth)
c create total velocity v
      do i=1,ii
       do j=1,jj
        v(i,j,k)=a(i,j)+vbaro(i,j)
        if(zbot(i,j).lt.0.1) v(i,j,k)=99999.
       enddo
      enddo
c read layer thickness dp
      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(a(1,1),idm,ii,jj,util,lgth)
      do i=1,ii
       do j=1,jj
        dp(i,j,k)=a(i,j)
       enddo
      enddo
c
      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(a(1,1),idm,ii,jj,util,lgth)
      do i=1,ii
       do j=1,jj
        temp(i,j,k)=a(i,j)
        if(zbot(i,j).lt.0.1) temp(i,j,k)=99999.
      enddo
      enddo
c
      read (ni,104,end=6) text,nstep,time,thet,ii,jj,
     .     lgth,(util(l),l=1,lgth)
c      write(*,100) text(1:8),nstep,time,thet,ii,jj,lgth,(util(l),l=1,14)
c      call unpakk(a(1,1),idm,ii,jj,util,lgth)
c      do i=1,ii
c       do j=1,jj
c        saln(i,j,k)=a(i,j)
c        if(zbot(i,j).lt.0.1) saln(i,j,k)=99999.
c       enddo
c      enddo
c
      enddo
c


c     write out u,v,dp in your favorite format here
c
      stop '(normal)'
6     stop '(not enough data to read)'
      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

    
      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

