c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine getgcmz(vol,nvol,istart,slev,fields,mtms,nlon,nlat,
     +  nf,iden,lu,iprhdr,ivol,spval,ln,ier)
c
c Define fields array from gcm hist at selected zp or interpolate to ht:
c   if slev <= zptop, then assumed to be zp, otherwise is selected height
c nf = number of fields (last dim of fields)
c      if nf=17, then fields array returned is 16 fields from hist + N2
c      if nf=18, then add fof2 on to fields in 18th slot
c iden and other parameters are as in getgcm()
c
      include "tgcmhdr.h"
      include "getgcm.h"
      character*(*) vol
      dimension mtms(3),fields(nlon,nlat,nf),vol(nvol),gcmht(imx,kmx),
     +  ln(nf),htscol(kmx),xnecol(kmx)
      character*128 flnm,flnm1
      logical isopen,exists
      data ivc/1/
      save ivc
c
c Check some inputs:     
c
      ier = 0
      if (lu.le.0.or.lu.gt.99.or.lu.eq.5.or.lu.eq.6) then
        write(6,"('>>> getgcmz: bad lu = ',i5)") lu
        ier = 2
        return
      endif
      if (mtms(2).lt.0.or.mtms(2).gt.23) then
        write(6,"('>>> getgcmz: bad model hour = ',i3,
     +    ' (mtms(2) = ',i5)") it,mtms(2)
        ier = 3
        return
      endif
      if (mtms(3).lt.0.or.mtms(3).gt.59) then
        write(6,"('>>> getgcmz: bad model minute = ',i3,
     +    ' (mtms(3) = ',i5)") it,mtms(3)
        ier = 4
        return
      endif
      if (nvol.le.0) then
        write(6,"('>>> getgcmz: need history volumes: nvol=',i3)")
     +    nvol
        ier = 5
        return
      endif
      if (iden.gt.2.or.iden.lt.0) then
        write(6,"('>>> getgcmz: bad iden=',i3,' (will default to ',
     +    ' no.densities (iden=1))')") iden
        iden = 1
      endif
      if (nlon.ne.imx.or.nlat.ne.jmx) then
        write(6,"('>>> warning getgcmz: nlon,nlat=',2i4,
     +    ' imx,jmx=',2i4,/' (these values should be the ',
     +    'same and are not)')") nlon,nlat,imx,jmx
      endif
      if (nf.gt.nflds+2) then
        write(6,"('>>> warning getgcmz: nf too big = ',i3,' max=',
     +    i3,' for 16 hist fields, plus N2, plus (optionally) FOF2')")
     +    nf,nflds+2
      endif
c     if (nf.eq.nflds+1) then
c       write(6,"('getgcmz will get ',i2,' fields: ',i2,
c    +    ' from hist + N2')") nf,nflds
c     elseif (nf.eq.nflds+2) then
c       write(6,"('getgcmz will get ',i2,' fields: ',i2,
c    +    ' from hist + N2 + FOF2')") nf,nflds
c     endif
c
c Define pressure column
c
      do k=1,kmx
        gcmzp(k) = zp1 + (k-1)*dzp
      enddo 
c
c If slev (selected level) <= gcmzp(kmx), then it is assumed to be
c   a pressure level (zp), and ixzp index is found. Otherwise
c   (slev > gcmzp(kmx), it assumed to be a selected height in which
c   case ixzp is left zero, and height interpolation is done to
c   return fields() at constant height (or spval if slev is above
c   model height at current grid point.
c
      ixzp = 0
      if (slev.le.gcmzp(kmx)) then
        if (slev.lt.gcmzp(1).or.slev.gt.gcmzp(kmx)) then
          write(6,"('>>> getgcmz: bad slev (assumed zp)=',f10.3,/
     +      ' Selected pressure slev must be ',f5.1,' < slev < ',
     +      f5.1)") slev,gcmzp(1),gcmzp(kmx)
          ier = 8
          return
        endif
        do k=1,kmx
          if (slev.lt.gcmzp(k)+0.5*dzp.and.slev.ge.gcmzp(k)-0.5*dzp) 
     +      ixzp=k
        enddo 
        if (ixzp.eq.0) then
          write(6,"('>>> getgcmz: could not find zp index to slev=',
     +      f10.2)") slev
          ier = 9
          return
        endif
      endif
c
c Set starting volume for search:
      if (istart.gt.0) then
        iv1 = 1
      else
        iv1 = ivc
      endif
c
c Volume search loop:
      do 300 iv=iv1,nvol
        ivc = iv
        inquire(lu,opened=isopen)
        if (istart.gt.0.or..not.isopen) then
          call getvol(vol(iv),lu,isopen,iier)
          if (iier.eq.2) then
            write(6,"('>>> getgcmz: getvol could not find vol ',a,
     +        ' on mss -- ',/12x,'skipping to next vol')") vol(iv)
            goto 300
          elseif (iier.ne.0) then
            write(6,"('>>> getgcmz: error ',i3,' from getvol')") iier
            ier = 6
            return
          endif 
          rewind(lu)
        endif
c
c Now we should be positioned for history search:
c
        it = 0
        lenvol = lenstr(vol(iv))
 100    continue			! search for history
        read(lu,end=900) iter,nday,nhr,nmin,label,
     +    date,output,start,stp,hist,sav,step,mag,difhor,iuivi,
     +    sdtide,ipower,aurora,dispos,data,source,sourct,dtide,
     +    dum,rdate,naurp,hp,cp,byimf
        it = it+1
        if (mtms(1).eq.nday.and.mtms(2).eq.nhr.and.
     +    mtms(3).eq.nmin) goto 200
        write(6,"('Getgcmz: Searching for ',i2,':',i2,':',i2,
     +    ' Found ',i2,':',i2,':',i2,' vol ',a,' it ',i3)")
     +    (mtms(i),i=1,3),nday,nhr,nmin,vol(iv)(1:lenvol),it
c
c Read past unwanted history:
c
        read(lu) dummy	! summary
        do j=1,nlat
          read(lu) dummy
        enddo
        goto 100
c
c Found desired history -- process lat slices:
c
 200    continue
        write(6,"(' ')")
        write(6,"('Getgcmz found history ',i2,':',i2,':',i2,' on vol ',
     +    a,' at it=',i3)") (mtms(i),i=1,3),vol(iv)(1:lenvol),it
        ivol = iv
        read(lu) dummy	! summary
        do j=1,nlat
          buffer in(lu,1)(frd(1),frd(nwlat))
          if (unit(lu)) 10,11,11
 11       write(6,"(' >>> Getgcmz: io problem buffering in j=',
     +      i3,' lu=',i3,' unit(lu)=',e10.3)") j,lu,unit(lu)
          ier = 7
          return
 10       continue
          call proclat(j,iden)
c
c Define field array:
c Do interpolation to height surface if necessary, as described above:
c
          if (ixzp.le.0) then
            do i=1,nlon
              gcmht(i,:) = f(i,:,ixz)
            enddo
          endif
          do ip=1,nflds
            if (ixzp.gt.0) then
              do i=1,nlon
                fields(i,j,ip) = f(i,ixzp,ip)
              enddo
            else
              call interpht(fields(1,j,ip),imx,f(1,1,ip),imxp1,kmx,
     +          gcmht,slev,ln(ip),spval)
            endif
          enddo
c
c N2 is in nflds+1 (17) slot:
          if (ixzp.gt.0) then
            do i=1,nlon
              fields(i,j,nflds+1) = gcmn2(i,ixzp)
            enddo
          else
            call interpht(fields(1,j,nflds+1),imx,gcmn2,imxp1,kmx,
     +        gcmht,slev,ln(nflds+1),spval)
          endif
c
c fof2 is optional at end of fields array:
c (hmf2 is stubbed in -- can add later if nf=nflds+3=19)
c
          if (nf.eq.nflds+2) then
            do i=1,nlon
              xnecol(:) = f(i,:,ixne) 
              htscol(:) = f(i,:,ixz) 
              call fof2int(xnecol,htscol,kmx,hmf2out,fof2out,0)
              fields(i,j,nflds+2) = fof2out
c             fields(i,j,nflds+3) = hmf2out
            enddo
          endif
        enddo	! j-loop
c
c Print summary of header if desired:
c
        if (iprhdr.gt.0) call printhdr(vol(iv),mtms(1))
c
c Save a link of the volume to /usr/tmp/TIGCM for possible future
c reference, and return:
c
        call mkflnm(vol(iv),flnm)
        flnm1 = "/usr/tmp/TIGCM/"//flnm(1:lenstr(flnm))
        inquire(file=flnm1,exist=exists)
        if (.not.exists) istat = link(flnm,"/usr/tmp/TIGCM")
        return
c
c Got end of file:
c
 900    write(6,"(' ')")
        write(6,"('Getgcmz: EOF encountered on unit ',i3,' it=',i3)") 
     +    lu,it
        write(6,"('        History ',3i4,' not found on vol ',a,/8x,
     +    '(closing unit',i3,')')") (mtms(i),i=1,3),
     +    vol(iv)(1:lenvol),lu
        close(lu)
 300  enddo
      write(6,"('Getgcmz: could not find history ',3i3,' The following',
     +  ' volumes were searched:')") (mtms(i),i=1,3)
      do i=1,nvol
        write(6,"(a)") vol(i)(1:lenstr(vol(i)))
      enddo
      ier = 1
      ivol = -1
      return
      end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine interpht(fout,imx,fin,imxp1,kmx,gcmht,sht,ln,spv)
c
c Return fout(imx) as field interpolated from fin(imxp1,kmx) to
c   selected height sht. Gcm heights are in gcmht(imx,kmx).
c   Return spv in fout(i) if sht > gcmht(i,maxht)
c This routine called at each j:
c             call interpht(fields(1,j,ip),imx,f(1,1,ip),imxp1,kmx,
c    +          gcmht,ht,slev,spval)
c
      dimension fout(imx),fin(imxp1,kmx),gcmht(imx,kmx)
c
c Lon loop:
      do i=1,imx
c
c Bracket sht:      
        do k=1,kmx-1
          if (sht.ge.gcmht(i,k).and.sht.le.gcmht(i,k+1)) then
            k0 = k
            k1 = k+1
            goto 100
          endif
        enddo
        fout(i) = spv
        goto 101
 100    continue
c
c Do linear or log interpolation:
c
        if (ln.le.0) then
          fout(i) = fin(i,k0) + (sht-gcmht(i,k0))*(fin(i,k1)-fin(i,k0))/
     +                           (gcmht(i,k1)-gcmht(i,k0))
        else
          exparg = (alog(fin(i,k1) / fin(i,k0)) / 
     +             (gcmht(i,k1) - gcmht(i,k0))) * (sht - gcmht(i,k0)) 
          fout(i) = fin(i,k0) * exp(exparg)
        endif
 101    continue
      enddo	! i=1,imx
      return
      end