c
      subroutine rdlsd(ngcmlon,gcmlon,ngcmlat,gcmlat,p0)
c
c Read a ccm2proc lsd file, return ccm fields on tgcm horizontal grid in 
c   fccm (pointer pfccm in common).
c Also fix units to match those coming from tgcm processor. This includes
c   omega to vertical velocity (see omega2vel).
c p0 appropriate for tgcm/t21 or tgcm/maccm2 is returned
c
      include "ccm.h"
c
c ccmzp1 and ccmdzp are hardwired because use of minihdr is not accurate
c enough. Also, bottom mb level has been fixed at 500mb (see PRESSLE in
c ccm2proc inputs, so the ln(p0/p) calculation comes out -20.64+ rather 
c than -20.5). These levels are same for t21 and maccm. 
c
      parameter(ccmzp1=-20.5, ccmdzp=0.5)
      real gcmlon(ngcmlon),gcmlat(ngcmlat)
      real fgcmglb(ngcmlon,ngcmlat)
      character*8 fname,fnameprev
      pointer (pdatrec,datrec(1)), (pihead,ihead(1)),
     +        (pfccmglb,fccmglb(1))
      integer minihdr(50)
      character*80 assgncmd,dskfile
      data ncalls/0/
      data fnameprev/'        '/, fname/'        '/
      save nltot,nccmtms,lenrec
c
      ncalls = ncalls+1
      if (ncalls.eq.1) then
        write(6,"(' ')")
        call clearstr(dskfile)
        call tail(ccmlsd,dskfile)	! disk file name is tail of mss path 
        call rdmss(ccmlsd,dskfile)
        call clearstr(assgncmd)
        write(assgncmd,"('assign -a ',a,' -s sbin fort.',i2)")
     +    dskfile(1:lenstr(dskfile)),luccm
        call assign(assgncmd(1:lenstr(assgncmd)))
c
c Get record length and read master header:
c
        read(luccm) lenrec
        rewind(luccm)
        call alloc(pihead,lenrec)
        read(luccm) (ihead(i),i=1,lenrec) ! master hdr (read full record)
c       write(6,"('First 11 words of master header:',/11i6)")
c    +    (ihead(i),i=1,11)
        nccmfld = ihead(3)		! number of fields
        nccmlev = ihead(4)		! number of levels each field
        nltot = ihead(5)		! *total* number of levels on file
        nccmtms = ihead(7)		! number of times on lsd file
        isrftype = ihead(9)		! surface type flag
        write(6,"('rdlsd read master header:'/'  lenrec=',i5,
     +    ' nccmfld=',i2,' nccmlev=',i2,' nltot=',i3,' nccmtms=',i3,
     +    ' isrftype=',i2)") lenrec,nccmfld,nccmlev,nltot,nccmtms,
     +    isrftype 
        call hpdeallc(pihead,ier,1)
      endif	! ncalls == 1
c
c Allocate space for a data record (local):
c (pdatrec is local, pccmlev_mb and pccmlev_zp are in common) 
c
      call alloc(pdatrec,lenrec)
c
c Allocate space for ccm level arrays (common):
c
      call alloc(pccmlev_mb,nccmlev)
      call alloc(pccmlev_zp,nccmlev)
c
c Read data records until EOF:
c
      nrec = 0
      ip = 1
      k = 0
      ixz2 = 0
      ixomega = 0
      do ilev = 1,nltot
        read(luccm,end=20) minihdr,(datrec(i),i=1,lenrec-50)
        write(fname,"(a8)") minihdr(7)
c
c Allocate space and intialize if reading first data record:
c
        if (nrec.eq.0.and.ilev.eq.1) then
          nccmlon = minihdr(14)
          nccmlat = minihdr(15)
          ccmday = float(minihdr(12)) + float(minihdr(13))/(24.*3600.)
c
c Allocate ccm field array (will be on gcm grid):
c fccm(ngcmlon,nccmlev,ngcmlat,nccmfld) 
c (pointer pfccm is in common, pfccmglb is local)
c
          if (ilev.eq.1) then
            call alloc(pfccm,ngcmlon*nccmlev*ngcmlat*nccmfld)
            call alloc(pfccmglb,nccmlon*nccmlat)
          endif
          fnameprev = fname
          ccmflab8(ip) = fname
c
c Get lon and lat arrays (pccmlon and pccmlat are in common):
c Note ccmlon(1) must be hardwired -- use of float(minihdr(16))/1000.
c   causes problems in grid transformation sub ccm2gcm (isrchgt, called
c   from tabulat returns 0). 
c         ccmlon(i) = float(minihdr(16))/1000.+dlon*float(i-1)
c
          if (ilev.eq.1) then
            call alloc(pccmlon,nccmlon)
            call alloc(pccmlat,nccmlat)
          endif
          dlon = 360./float(nccmlon)
          ccmlon(1) = -180.
          do i=1,nccmlon
            ccmlon(i) = ccmlon(1)+dlon*float(i-1)
          enddo 
c
c ccm latitudes are hard-wired in glat_t21 and glat_maccm according to 
c t21 or maccm. Use number of latitudes to determine T21 vs MACCM:
c
          if (nccmlat.eq.nlat_t21.and.nccmlon.eq.nlon_t21) then
            do j=1,nccmlat
              ccmlat(j) = glat_t21(j)
            enddo
            p0 = 5.412154e-4
            ccmres = 'CCM2-T21        '
          elseif (nccmlat.eq.nlat_maccm.and.nccmlon.eq.nlon_maccm) then
            do j=1,nccmlat
              ccmlat(j) = glat_maccm(j)
            enddo
            p0 = 4.15872e-4
            ccmres = 'CCM2-MACCM      '
          elseif (nccmlat.eq.nlat_t42.and.nccmlon.eq.nlon_t42) then
            do j=1,nccmlat
              ccmlat(j) = glat_maccm(j)	! lats same for ccm3t42 and maccm2
            enddo
            p0 = 5.412154e-4	! interface is same for ccm2t21 and ccm3t42
            ccmres = 'CCM3-T42        '
          else
            write(6,"('>>> unknown latitude resolution: nccmlat=',
     +        ' nlat_t21=',i2,' nlat_maccm=',i2)") nccmlat,
     +        nlat_t21,nlat_maccm
            stop 'nlat'
          endif 
c         write(6,"('nccmfld=',i2,' nccmlev=',i2,' nccmlon=',
c    +      i2,' nccmlat=',i2,' ccmday=',f9.4,' ',a)")
c    +      nccmfld,nccmlev,nccmlon,nccmlat,ccmday,
c    +      ccmres(1:lenstr(ccmres))
        endif	! nrec==0 and ilev==1
c
c Update level and field indices, define ccmlev if first field:
c       ccmlev_zp(k) = alog(p0*1.e-3/ccmlev_mb(k))	! ccmlev in zp
c
        k = k+1
        if (fname.ne.fnameprev) then
          ip = ip+1	! field index
          ccmflab8(ip) = fname
c
c ixz2 and ixomega are field indices to Z2 and OMEGA in fccm, for use
c by omega2vel after all fields have been read:
c
          if (fname.eq.'Z2      ') ixz2 = ip
          if (fname.eq.'OMEGA   ') ixomega = ip
          k = 1
        endif
        if (ip.eq.1) then
          ccmlev_mb(k) = float(minihdr(6))/1000.	! ccmlev in mb
          ccmlev_zp(k) = ccmzp1+(k-1)*ccmdzp		! ccmlev in zp
        endif
c
c Define fccmglb(nccmlon,nccmlat) for current level and current field:
c
        do j=1,nccmlat
          do i=1,nccmlon
            fccmglb((j-1)*nccmlon+i) = datrec((j-1)*nccmlon+i)
          enddo
        enddo
c
c Bilinearly interpolate ccm field to tgcm geographic grid, and
c   define fccm(ngcmlon,nccmlev,ngcmlat,nccmfld) at current k level 
c   and ip field. Convert Z2 from m to km (U,V already in m/s).
c
        ivec = 0
        if (fname.eq.'U       '.or.fname.eq.'V       ') ivec = 1
c
c 3/19/96: new interp routine (uses John Adams bicubic interp)
c
        call ccm2gcm(ngcmlon,ngcmlat,nccmlon,nccmlat,1,ccmlon,ccmlat,
     +    gcmlon,gcmlat,fccmglb,fgcmglb,ivec)
c
c Convert geopotential height from m to km:
c
        scale = 1.
        if (fname.eq.'Z2      ') scale = 1.e-3	! m -> km
c
c Fccm to contain ccm data at tgcm grid:
c
        do j=1,ngcmlat
          do i=1,ngcmlon
            fccm((ip-1)*ngcmlat*nccmlev*ngcmlon+(j-1)*nccmlev*ngcmlon+
     +        (k-1)*ngcmlon+i) = fgcmglb(i,j)*scale
          enddo
        enddo
c
c Report to stdout and return for next data record:
c
        if (k.eq.nccmlev.and.ncalls.eq.1)
     +    write(6,"('  completed read of ccm field ',a,' ccmday=',
     +      f12.5)") ccmflab8(ip),ccmday
        nrec = nrec+1
        fnameprev = fname
      enddo	! ilev=1,nltot
c
c Convert ccm omega to vertical velocity (m/s):
c   fccm(ngcmlon,nccmlev,ngcmlat,nccmfld)
c
      if (ixz2.eq.0) then
        write(6,"(/'>>> Need ccm2 field Z2 for OMEGA conversion')")
        stop 'Z2'
      endif
      call omega2vel(fccm((ixomega-1)*ngcmlat*nccmlev*ngcmlon+1),
     +  fccm((ixz2-1)*ngcmlat*nccmlev*ngcmlon+1),ccmlev_mb,
     +  ngcmlat,ngcmlon,nccmlev)
c
c Release local space:
c
      call hpdeallc(pfccmglb,ier,1)
      call hpdeallc(pdatrec,ier,1)
c
      return
 20   continue
      write(6,"('EOF at nrec=',i4)") nrec 
      stop 'lsdeof'
      end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine omega2vel(omega,hts,press,nlat,nlon,nlev)
c
c Convert ccm omega from pascals/sec to velocity m/s:
c Input: 
c   omega(nlon,nlev,nlat) = ccm2 OMEGA (Pa/s)
c   hts(nlon,nlev,nlat) = ccm2 geopotential heights (km)
c   press(nlev) = pressures (mb)
c Ouput:
c   omega has been changed to m/s (can now be plotted with tgcm proc's W)
c   (hts and press are unchanged)
c        
      real omega(nlon,nlev,nlat),hts(nlon,nlev,nlat),press(nlev)
      real pzps(nlev)	! scale height (local)
c
      do j=1,nlat
        do i=1,nlon
          do k=2,nlev-1
            pzps(k) = hts(i,k+1,j) - hts(i,k-1,j)
          enddo
          pzps(1) = 4.*hts(i,2,j) - 3.*hts(i,1,j) - hts(i,3,j)
          pzps(nlev) = 3.*hts(i,nlev,j) - 4.*hts(i,nlev-1,j) +
     +      hts(i,nlev-2,j)
c
c New omega: use pzps in m (is in km), and press in Pa (is in mb):
c (1 mb = 100 Pa)
c
          omega(i,:,j) = -omega(i,:,j)*pzps(:)*1000./(press(:)*100.)
        enddo
      enddo
      return
      end