c
c array indices are: first index is latitude, going from north to south
c                    second index is longitude, from west to east
c
c fields are read in the following order:
c ubaro (southward barotropic velocity in cm/s)
c vbaro (eastward barotropic velocity in cm/s)
c pbaro (barotropic pressure)
c montg (mixed layer Montgomery potential) 
c thmix (mixed layer density in sigm_theta units)
c the sea surface height in cm/s is calculated as montg+thref*pbaro/g
c do k=1,16
c u     (southward baroclinic velocity in cm/s)
c v     (eastward baroclinic velocity in cm/s)
c dp    (layer thickness in pressure units, dp/g is in cm)
c tem   (layer temperature, degrees Celsius)
c saln  (layer salinity in parts per thousand)
c enddo
c other cycle for mass exchange between layers
c
c For a smaller region, set istart, jstart, ie, je
c read and unpakk each variable,
c pakk and write variables for the smaller region
c-----------------------------------------------------------------------------
c select a subregion and put the index limits into istart, etc
c is1, ie1, js1, je1 can correspond to an already extracted data set
c
      parameter (is1=1,ie1=1437,js1=1,je1=1437)
      parameter (istart=???,ie=???,jstart=???,je=???)
      parameter (idm=ie-istart+1,jdm=je-jstart+1,idm1=ie1-is1+1,
     .           jdm1=je1-js1+1,kdm=16,ms=30,msd=30,maxp=5000)
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-----------------------------------------------------------------------------
      parameter (ifull=1437,jfull=1437)
c
      common/fields/u(idm,jdm,2*kdm),v(idm,jdm,2*kdm),dp(idm,jdm,kdm),
     .              temp(idm,jdm,2*kdm),saln(idm,jdm,2*kdm),theta(kdm),
     .              p(idm,jdm,kdm+1),thmix(idm,jdm),srfht(idm,jdm),
     .              ubaro(idm,jdm),vbaro(idm,jdm),pbaro(idm,jdm),
     .              depths(idm,jdm),scu(idm),scv(idm),scp(idm),
     .              scu2(idm),scv2(idm),scp2(idm),montg(idm,jdm),
                    ,subduc(idm,jdm,kdm),trac(idm,jdm,2*kdm)
      real wrkfull(ifull,jfull),wrk1(idm1,jdm1)
      character util(ifull*jfull+14)*2,flnm*60,text*8,preambl(5)*79
      logical initl
c
c --- 'reflon' = longitude of grid point 'jref'
c --- 'gridsz' = grid size in degrees longitude
c --- 'equat'  = i index giving location of equator in grid
c --- 'thref' = reference specific volume (cm^3/g)
c --- 'g' = gravity (cm/s^2)
c
      data jref,reflon,equat,gridsz,lalolb/1,-98.,1071.,.08,5/  
      data thref/1./,g/980.6/
      data radian/57.29578/,pi/3.14159265/
c
c --- input and output file
      data ni/15/,no/16/
c
      common/geog/jref,reflon,equat,gridsz,lalolb
c
c replace by relative values if needed
c
      equat=equat-istart+1
      jref=jref-jstart+1
c
c --- latitude as a function of grid index
c      xlat(i)=2.*atan(exp(gridsz*(equat-i)/radian))-pi/2.
c --- grid index xi as a function of latitude 
c      xi(alat)=equat-alog(tan((2.*alat/radian+pi)/4.))*radian/gridsz
c
c --- longitude as a function of grid index
c      ylon(j)=reflon+gridsz*(j-jref)
c --- grid index yj as a function of longitude
c      xj(alon)==jref+(alon-reflon)/gridsz
c
 104  format (3x,a8,i8,2f8.0,2i4,i6/(1x,64a2))
 114  format (3x,a8,i8,f8.1,f8.1,2i4,i6/(1x,64a2))
 100  format (' ::',a8,i8,f8.1,f8.5,2i4,i6,1x,14a2)
 110  format (' ::',a8,i8,f8.1,f8.1,2i4,i6,1x,14a2)
 120  format (5(a79/),3i8/(40a2))
c
c --- acquire basin depths for the area
c
      open (unit=9,file='depthdamee.1437x1437'
     .     ,status='old',readonly)
      read (9,120,end=6) preambl,i1,j1,lgth,(util(l),l=1,lgth)
      write(6,*)'i1,j1,lgth',i1,j1,lgth
      call unpakk(wrkfull,i1,i1,j1,util,lgth)
      do i=1,idm
      do j=1,jdm
      depths(i,j)=wrkfull(i+istart-1,j+jstart-1)
      enddo
      enddo
      close (unit=9)
c
c --- should be saved to a separate file than the data
c
c --- open file
c 
      write(6,*)'enter year and day'
      read(5,*)iyr,iday
      flnm='0001,0001.15.003'
      write(flnm(11:12),'(i2.2)')iyr
      write(flnm(14:16),'(i3.3)')iday
      write(6,*)'flnm',flnm
      flnm='/h/ncar/worley/'//flnm
      open (unit=ni,name=flnm,status='old',form='formatted')
 13   continue
c      
c --- now unpack barotropic velocity field
c
      read (ni,104,end=6) text,nstep,time,thet,i,j,
     .     lgth,(util(l),l=1,lgth)
      write (*,100) text(1:8),nstep,time,thet,i,j,lgth,(util(l),l=1,14)
      write(*,*)'before unpakk ubaro ii i',ii,i
      call unpakk(wrk1,idm1,i,j,util,lgth)
      do i=1,idm
      do j=1,jdm
      ubaro(i,j)=wrk1(i+istart-is1,j+jstart-js1)
      enddo
      enddo
c 
c --- pakk new file and write to unit no
c 
      call pakk(ubaro,idm,idm,jdm,util,length)
      write (no,104) textnstep,time,thet,idm,jdm,
     .     lgth,(util(l),l=1,lgth)
c
      read (ni,104,end=6) text,nstep,time,thet,i,j,
     .     lgth,(util(l),l=1,lgth)
      write (*,100) text(1:8),nstep,time,thet,i,j,lgth,(util(l),l=1,14)
      call unpakk(wrk1,idm1,i,j,util,lgth)
      do i=1,idm
      do j=1,jdm
      vbaro(i,j)=wrk1(i+istart-is1,j+jstart-js1)
      enddo
      enddo
c 
c --- pakk new file and write to unit no
c 
      call pakk(vbaro,idm,idm,jdm,util,length)
      write (no,104) textnstep,time,thet,idm,jdm,
     .     lgth,(util(l),l=1,lgth)
c
      read (ni,104,end=6) text,nstep,time,thet,i,j,
     .     lgth,(util(l),l=1,lgth)
      write (*,100) text(1:8),nstep,time,thet,i,j,lgth,(util(l),l=1,14)
      call unpakk(wrk1,idm1,i,j,util,lgth)
      do i=1,idm
      do j=1,jdm
      pbaro(i,j)=wrk1(i+istart-is1,j+jstart-js1)
      enddo
      enddo
c 
c --- pakk new file and write to unit no
c 
      call pakk(pbaro,idm,idm,jdm,util,length)
      write (no,104) textnstep,time,thet,idm,jdm,
     .     lgth,(util(l),l=1,lgth)
c
      read (ni,104,end=6) text,nstep,time,thet,i,j,
     .     lgth,(util(l),l=1,lgth)
      write (*,100) text(1:8),nstep,time,thet,i,j,lgth,(util(l),l=1,14)
      call unpakk(wrk1,idm1,i,j,util,lgth)
      do i=1,idm
      do j=1,jdm
      montg(i,j)=wrk1(i+istart-is1,j+jstart-js1)
      enddo
      enddo
c 
c --- pakk new file and write to unit no
c 
      call pakk(montg,idm,idm,jdm,util,length)
      write (no,104) textnstep,time,thet,idm,jdm,
     .     lgth,(util(l),l=1,lgth)
c
      read (ni,104,end=6) text,nstep,time,thet,i,j,
     .     lgth,(util(l),l=1,lgth)
      write (*,100) text(1:8),nstep,time,thet,i,j,lgth,(util(l),l=1,14)
      call unpakk(wrk1,idm1,i,j,util,lgth)
      do i=1,idm
      do j=1,jdm
      thmix(i,j)=wrk1(i+istart-is1,j+jstart-js1)
      enddo
      enddo
c 
c --- pakk new file and write to unit no
c 
      call pakk(thmix,idm,idm,jdm,util,length)
      write (no,104) textnstep,time,thet,idm,jdm,
     .     lgth,(util(l),l=1,lgth)
c
      do 14 k=1,kdm
      read (ni,104,end=6) text,nstep,time,thet,i,j,
     .     lgth,(util(l),l=1,lgth)
      write (*,100) text(1:8),nstep,time,thet,i,j,lgth,(util(l),l=1,14)
      call unpakk(wrk1,idm1,i,j,util,lgth)
      do i=1,idm
      do j=1,jdm
      u(i,j,k)=wrk1(i+istart-is1,j+jstart-js1)
      enddo
      enddo
c 
c --- pakk new file and write to unit no
c 
      call pakk(u(1,1,k),idm,idm,jdm,util,length)
      write (no,104) textnstep,time,thet,idm,jdm,
     .     lgth,(util(l),l=1,lgth)
c
      read (ni,104,end=6) text,nstep,time,thet,i,j,
     .     lgth,(util(l),l=1,lgth)
      write (*,100) text(1:8),nstep,time,thet,i,j,lgth,(util(l),l=1,14)
      call unpakk(wrk1,idm1,i,j,util,lgth)
      do i=1,idm
      do j=1,jdm
      v(i,j,k)=wrk1(i+istart-is1,j+jstart-js1)
      enddo
      enddo
c 
c --- pakk new file and write to unit no
c 
      call pakk(v(1,1,k),idm,idm,jdm,util,length)
      write (no,104) textnstep,time,thet,idm,jdm,
     .     lgth,(util(l),l=1,lgth)
c
      read (ni,104,end=6) text,nstep,time,thet,i,j,
     .     lgth,(util(l),l=1,lgth)
      write (*,100) text(1:8),nstep,time,thet,i,j,lgth,(util(l),l=1,14)
      call unpakk(wrk1,idm1,i,j,util,lgth)
      do i=1,idm
      do j=1,jdm
      dp(i,j,k)=wrk1(i+istart-is1,j+jstart-js1)
      enddo
      enddo
c 
c --- pakk new file and write to unit no
c 
      call pakk(dp(1,1,k),idm,idm,jdm,util,length)
      write (no,104) textnstep,time,thet,idm,jdm,
     .     lgth,(util(l),l=1,lgth)
c
      read (ni,104,end=6) text,nstep,time,thet,i,j,
     .     lgth,(util(l),l=1,lgth)
      write (*,100) text(1:8),nstep,time,thet,i,j,lgth,(util(l),l=1,14)
      call unpakk(wrk1,idm1,i,j,util,lgth)
      do i=1,idm
      do j=1,jdm
      temp(i,j,k)=wrk1(i+istart-is1,j+jstart-js1)
      enddo
      enddo
c 
c --- pakk new file and write to unit no
c 
      call pakk(temp(1,1,k),idm,idm,jdm,util,length)
      write (no,104) textnstep,time,thet,idm,jdm,
     .     lgth,(util(l),l=1,lgth)
c
      read (ni,104,end=6) text,nstep,time,thet,i,j,
     .     lgth,(util(l),l=1,lgth)
      write (*,100) text(1:8),nstep,time,thet,i,j,lgth,(util(l),l=1,14)
      call unpakk(wrk1,idm1,i,j,util,lgth)
      do i=1,idm
      do j=1,jdm
      saln(i,j,k)=wrk1(i+istart-is1,j+jstart-js1)
      enddo
      enddo
c 
c --- pakk new file and write to unit no
c 
      call pakk(saln(1,1,k),idm,idm,jdm,util,length)
      write (no,104) textnstep,time,thet,idm,jdm,
     .     lgth,(util(l),l=1,lgth)
c
14    continue
c
      ii=idm
      jj=jdm
      ii1=ii-1
      jj1=jj-1
      do 40 j=1,jj1
      do 40 i=1,ii1
 40   srfht(i,j)=(montg(i,j)+thref*pbaro(i,j))/g
c
c     return
c
 6    continue
      stop '(e-o-f)'
      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.)
ccc      if (i1.gt.37) i1=i1+6
ccc      if (i1.gt.11) i1=i1+7
ccc      i1=i1+14
ccc      if (i2.gt.37) i2=i2+6
ccc      if (i2.gt.11) i2=i2+7
ccc      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)
ccc      if (i1.gt.96) i1=i1-6
ccc      if (i1.gt.64) i1=i1-7
ccc      i1=i1-14
ccc      if (i2.gt.96) i2=i2-6
ccc      if (i2.gt.64) i2=i2-7
ccc      i2=i2-14
c
 4    array(i,j)=scal*float(64*(i1-32)+(i2-32))+base
      if (lngth.ne.length) then
      write(6,*)'lngth',lngth,'length',length
      stop 'UNPACK'
      endif
      return
c
      ENTRY UNPAKK1(ARRAY,IDIM,II,JJ,COMPAC,LENGTH)
c
      do 19 l=1,14
 19     comp2(l)=compac(l)
ccc      read (comp14(1),101) base
ccc      read (comp14(2),101) scal
      lngth=14
      do 14 i=1,ii
        do 14 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)
ccc      if (i1.gt.96) i1=i1-6
ccc      if (i1.gt.64) i1=i1-7
ccc      i1=i1-14
ccc      if (i2.gt.96) i2=i2-6
ccc      if (i2.gt.64) i2=i2-7
ccc      i2=i2-14
c
ccc 14    array(i,j)=scal*float(64*(i1-32)+(i2-32))+base
  14   array(i,j)=0.
      if (lngth.ne.length) stop 'UNPACK1'
      return
      end