c
c------------------------------------------------------------------
c Begin file /home/sting/foster/timegcm/mksrc/acd2d.f
c------------------------------------------------------------------
c
      subroutine acd2d(nspec,ispecno,acdname,hist,luhist,gcmlat,jmx,
     +   htlb,speclb,rholb,ispecpsi,iplt,ieqnx,iframe)
c
c Get acd 2d model history and find lower boundary conditions
c for selected species.
c
c On input:
c   nspec = number of desired species
c   ispecno(nspec) = numbers corresponding to desired species
c   acdname(nspec) = char*8 field names
c   hist = char*32 mss path name of desired history
c   luhist = fortran logical unit for history
c   gcmlat(jmx) = tgcm latitude grid
c   htlb = desired height for lower boundaries
c   rholb(jmx) = total density (rho) at gcm lat grid in no. density
c               (from cira88)
c
c On output:
c   speclb(jmx,nspec) = acd 2d values of species at htlb height
c                       gcmlat latitude grid
c
c Note: this routine distilled from RAJA processor sol9.f (see ~/acd2d)
c
      dimension ispecno(nspec),speclb(jmx,nspec),gcmlat(jmx),rholb(jmx),
     +  ispecpsi(nspec)
      character*8 acdname(nspec)
      character*16 lochist
      character*32 hist
      character*80 errmsg
      parameter(mxspec=63, lmax=35, niz=86, ndays=26)
      parameter(acdlat1=-85.,dlat=5.)
      parameter(acdalt1=0.,dalt=1.)
      parameter(iday1=1087, delday=15)
      dimension acdlat(lmax),acdalt(niz),acddata(mxspec,lmax,niz),
     +          ispecall(mxspec)
      data ilog/1/
c
c Use acd day 1462 (1/2) for solstice, 1177 (3/23) for equinox:
c
      igetday = 1462
      if (ieqnx.gt.0) igetday = 1177
c
c Check some input values:
c (luhist cannot be 2 because of conflict with gmeta lu)
c
      if (luhist.lt.1.or.luhist.eq.5.or.luhist.eq.6.or.luhist.eq.2) then
        write(6,"('>>> acd2d: bad luhist=',i3)") luhist
        stop 'luhist'
      endif
      if (nspec.lt.1.or.nspec.gt.mxspec) then
        write(6,"('>>> acd2d: bad nspec=',i5)") nspec
        stop 'nspec'
      endif
      do i=1,nspec
        if (ispecno(i).lt.1.or.ispecno(i).gt.mxspec) then
          write(6,"('>>> acd2d: bad ispecno(',i2,') = ',i4)")
     +      i,ispecno(i)
          stop 'ispecno'
        endif
      enddo
c
c Set up acd grid:
      do j=1,lmax
        acdlat(j) = acdlat1 + (j-1)*dlat
      enddo
      do k=1,niz
        acdalt(k) = acdalt1 + (k-1)*dalt
      enddo
      do i=1,mxspec
        ispecall(i) = i
      enddo
      if (htlb.lt.acdalt(1).or.htlb.gt.acdalt(niz)) then
        write(6,"('>>> acd2d: bad htlb=',e12.4,' acdalt(1)=',
     +    e12.4,' acdalt(',i3,')=',e12.4)") htlb,acdalt(1),
     +    niz,acdalt(niz)
        stop 'htlb'
      endif
c
c The acd model vertical grid is 0 -> 85 km by 1 km
c Lower boundary for timegcm  (htlb) is 62 km, so 
c there is no need to interpolate. ixhtlb is index
c to htlb in acddata:  
c
      do k=2,niz
        if (htlb.ge.acdalt(k-1).and.htlb.le.acdalt(k)) then
          ixhtlb = k-1
          goto 202
        endif
      enddo
      write(6,"('>>> acd2d: could not find acdalt index to ',
     +  'htlb=',e12.4,' acdalt=',/(6f9.3))") htlb,acdalt
      stop 'acdalt'
 202  continue
c
c Get acd history volume:
      call tail(hist,lochist)
      write(6,"(' ')")
      write(6,"('acd2d acquiring acd history file ',a,
     +   ' to local disk file ',a)") hist(1:lenstr(hist)),lochist
      call msread(ier,lochist,hist,' ',' ')
      if (ier.eq.3) then
        write(6,"('acd2d: acd history file ',a,' was found on',
     +    ' the disk, and not read from mss')") hist
      elseif (ier.ne.0) then
        call mserror(errmsg)
        write(6,"('>>> acd2d: error acquiring acd history ',a,
     +    ' errmsg=',a)") hist,errmsg
        stop 'acdhist'
      else
        write(6,"('acd2d: successful mss read of acd history ',a)")
     +    hist
      endif
c
c Assign it to luhist:
      call assgn(lochist,luhist)
c
c Read through history to desired day:
c     parameter(iday1=1087, delday=15, igetday=1462) (ndays=26)
c
      do id=1,ndays
        call gdata(acddata,iday,luhist,ispecall,mxspec)
        write(6,"('acd2d: id=',i3,' searching for day=',i5,
     +    ' read day ',i5)") id,igetday,iday
        if (iday.eq.igetday) then
          write(6,"('acd2d: found day ',i5,' on acd hist ',a)")
     +      igetday,hist
          goto 100
        endif
      enddo
      write(6,"('>>> acd2d: could not find day ',i5,' on acd history ',
     +  a)") igetday,hist
      stop 'acdhist'
 100  continue
c
c Interpolate to desired latitude grid:
c
      do j=1,jmx
        if (gcmlat(j).lt.acdlat(1)) then
          j1 = 1
          j2 = 2
          goto 201
        elseif (gcmlat(j).gt.acdlat(lmax)) then
          j1 = lmax-1
          j2 = lmax
          goto 201
        else
          do jj=2,lmax
            if (gcmlat(j).ge.acdlat(jj-1).and.
     +          gcmlat(j).le.acdlat(jj))then
              j1 = jj-1
              j2 = jj
              goto 201
            endif
          enddo
        endif
        write(6,"('>>> acd2d: could not find acdlat index to ',
     +    ' gcmlat(j)=',f10.3,' acdlat=',/(6f9.3))") gcmlat(j),acdlat
        stop 'acdlat'
 201    continue
c
c Do interpolation to tgcm latitude: 
c
        do ip=1,nspec  
          f1 = (gcmlat(j)-acdlat(j1)) / (acdlat(j2)-acdlat(j1))
          speclb(j,ip) = f1*acddata(ispecno(ip),j2,ixhtlb)+
     +                 (1.-f1)*acddata(ispecno(ip),j1,ixhtlb)
        enddo
c
c End latitude loop:
      enddo
c
c Some species on acd histories are in mass mix ratios (ispecpsi(ip)=1),
c   and some are in number density (ispecpsi(ip)=0) -- 
c   We need to convert the number density species to mass mixing ratios 
c   for the data volume (use rho from zatmos (rholb(jmx) is in no.dens)):
c
      do ip=1,nspec
        if (ispecpsi(ip).le.0) speclb(:,ip) = speclb(:,ip) / rholb(:) 
      enddo
c
c Contour and line plot acd 2d data, if desired:
c
      if (iplt.gt.0) then
        call pltacd(acddata,mxspec,nspec,lmax,niz,acdalt,acdlat,
     +              acdname,htlb,ixhtlb,ieqnx,ispecno,ilog,cpspval,
     +              iframe)
      endif
c
      return
      end
c23456789012345678901234567890123456789012345678901234567890123456789012 
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c Routine to read acd history tape (from ~/acd2d/sol9.f)
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine gdata(data,day,unit,spcno,nspec)
          parameter(niz=86,lmax=35,nbcon=38,ncm=24)
          real data(nspec,lmax,niz),tol(lmax)
          integer day,unit,iz,spcno(nspec)
          integer retcod
          real dummy(lmax,niz)
c
      do 120 ic=1,nbcon
        call arread(dummy,niz,lmax,unit,day,retcod)
        if (retcod .ne. 0) stop 'hmu'
        do 110 nd=1,nspec
          if (spcno(nd) .eq. ic) then
            do 100 lm=1,lmax
            do 100 iz=1,niz
               data(nd,lm,iz)=dummy(lm,iz)
  100       continue
c           print*,'ic,spcno(nd),data(nd,lmax,niz)'
c           print*,ic,spcno(nd),data(nd,lmax,niz)
          endif
  110   continue
  120 continue
c
      do 220 nc=1,ncm
        call arread(dummy,niz,lmax,unit,day,retcod)
        if (retcod .ne. 0) stop 'nmu'
        do 210 nd=1,nspec
          if (spcno(nd) .eq. (nc+nbcon)) then
            do 200 lm=1,lmax
            do 200 iz=1,niz
              data(nd,lm,iz)=dummy(lm,iz)
  200       continue
          endif
  210   continue
  220 continue
c
      call arread(dummy,niz,lmax,unit,day,retcod)
      if (retcod .ne. 0) stop 'temp'
      if(spcno(nspec).eq.63 ) then
        do 300 lm=1,lmax
        do 300 iz=1,niz
          data(nspec,lm,iz)=dummy(lm,iz)
  300   continue
      endif
c
      call arread(tol,1,lmax,unit,day,retcod)
      if (retcod .ne. 0) stop 'tol'
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine arread( data, np, nl, iunit, iday, retcod )
      integer     np, nl, iunit, retcod, iday
      integer     j, k, nitems, iparm, id2
      real        data(nl,np), dd(86)
c     read(iunit)  iday, iparm, id2, nitems,
      read(iunit,end= 100)  iday, iparm, id2, nitems,
     -             ( dd(k) , (data(j,k),j = 1,nl), k = 1,np )

      if( np*nl .ne. nitems ) then
         retcod = -1
      else
         retcod = 0
      end if
c!    write(8,600) iday
      return
  100 write(8,800) 
      retcod = 9
      return
  600 format(1x,' arread iday=',i4)    
800   format(1x,' end of file')    
      end