c
      subroutine rdlsd(ngcmlon,gcmlon,ngcmlat,gcmlat,p0)
c
c Read a ccm2proc lsd file, return ccm fields on tgcm horizontal grid in 
c   ccmf (pointer pccmf 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 fgcm(ngcmlon-1,ngcmlat),fgcmp(ngcmlon,ngcmlat)
      character*8 fname,fnameprev
      pointer (pdatrec,datrec(1)), (pihead,ihead(1)),
     +        (pfccm,fccm(1))
      integer minihdr(50)
      character*80 assgncmd,dskfile
c
      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 fields
c
c Allocate space for a data record, and ccmlev arrays:
c (pdatrec is local, pccmlev_mb and pccmlev_zp are in common) 
c
      call alloc(pdatrec,lenrec)
      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
 10   continue      
      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) 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 ccmf(ngcmlon,ngcmlat,nccmlev,nccmfld) 
c (pointer pccmf is in common, pfccm is local)
c
        call alloc(pccmf,ngcmlon*ngcmlat*nccmlev*nccmfld)
        call alloc(pfccm,nccmlon*nccmlat)
        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
        call alloc(pccmlon,nccmlon)
        call alloc(pccmlat,nccmlat)
        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) then
          do j=1,nccmlat
            ccmlat(j) = glat_t21(j)
          enddo
          p0 = 5.412154e-4
          ccmres = 'CCM2-T21        '
        elseif (nccmlat.eq.nlat_maccm) then
          do j=1,nccmlat
            ccmlat(j) = glat_maccm(j)
          enddo
          p0 = 4.15872e-4
          ccmres = 'CCM2-MACCM      '
        else
          write(6,"('>>> unknown latitude resolution: nccmlat=',
     +      ' nlat_t21=',i2,' nlat_maccm=',i2)") nccmlat,
     +      nlat_t21,nlat_maccm
          stop 'nlat'
        endif 
        write(6,"('rdlsd: nccmfld=',i2,' nccmlev=',i2,' nccmlon=',
     +    i2,' nccmlat=',i2,' ccmday=',f9.4,' ccmres=',a)")
     +    nccmfld,nccmlev,nccmlon,nccmlat,ccmday,
     +    ccmres(1:lenstr(ccmres))
      endif	! first record
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 (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
      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 ccmf, 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
c
c Define fccm(nccmlon,nccmlat) for current level and current field:
c
      do j=1,nccmlat
        do i=1,nccmlon
          fccm((j-1)*nccmlon+i) = datrec((j-1)*nccmlon+i)
        enddo
      enddo
      call fminmax(fccm,nccmlat*nccmlon,fmin,fmax,1.e36)
c
c Bilinearly interpolate ccm field to tgcm geographic grid, and
c   define ccmf(ngcmlon,ngcmlat,nccmlev,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
      call ccm2gcm(ngcmlon-1,ngcmlat,nccmlon,nccmlat,1,ccmlon,ccmlat,
     +  gcmlon,gcmlat,fccm,fgcm,ivec)
      do i=1,ngcmlon-1
        fgcmp(i,:) = fgcm(i,:)
      enddo
      fgcmp(ngcmlon,:) = fgcmp(1,:) ! longitude wrap-around
      scale = 1.
      if (fname.eq.'Z2      ') then
        scale = 1.e-3	! m -> km
      endif
c   ccmf(ngcmlon,ngcmlat,nccmlev,nccmfld)
      do j=1,ngcmlat
        do i=1,ngcmlon
          ccmf((ip-1)*nccmlev*ngcmlat*ngcmlon+(k-1)*ngcmlat*ngcmlon+
     +      (j-1)*ngcmlon+i) = fgcmp(i,j)*scale
        enddo
      enddo
c
c Report to stdout and return for next data record:
c
      if (k.eq.nccmlev)
     +  write(6,"('  completed read of ccm field ',a)") ccmflab8(ip)
      nrec = nrec+1
      fnameprev = fname
      goto 10
 20   continue
      write(6,"('EOF at nrec=',i4)") nrec 
c
c Convert ccm omega to vertical velocity (m/s):
c
      if (ixz2.eq.0) then
        write(6,"(/'>>> Need ccm2 field Z2 for OMEGA conversion')")
        stop 'Z2'
      endif
      call omega2vel(ccmf((ixomega-1)*nccmlev*ngcmlat*ngcmlon+1),
     +  ccmf((ixz2-1)*nccmlev*ngcmlat*ngcmlon+1),ccmlev_mb,
     +    ngcmlat,ngcmlon,nccmlev)

c   ccmf(ngcmlon,ngcmlat,nccmlev,nccmfld)
c     do kk=1,nccmlev
c       ix = (ixz2-1)*nccmlev*ngcmlon*ngcmlat+(kk-1)*ngcmlon*ngcmlat+1
c       aveht = fglbm(ccmf(ix),ngcmlon,ngcmlat,gcmlat,5.,5.,1.e36)
c       write(6,"('rdlsd: kk=',i2,' ixz2=',i2,' ix=',i7,' ccm aveht=',
c    +    f10.2)") kk,ixz2,ix,aveht
c     enddo

c
c Release local space:
c
      call hpdeallc(pihead,ier,1)
      call hpdeallc(pdatrec,ier,1)
      call hpdeallc(pfccm,ier,1)
      return
      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,nlat,nlev) = ccm2 OMEGA (Pa/s)
c   hts(nlon,nlat,nlev) = 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,nlat,nlev),hts(nlon,nlat,nlev),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,j,k+1) - hts(i,j,k-1)
          enddo
          pzps(1) = 4.*hts(i,j,2) - 3.*hts(i,j,1) - hts(i,j,3)
          pzps(nlev) = 3.*hts(i,j,nlev) - 4.*hts(i,j,nlev-1) +
     +      hts(i,j,nlev-2)
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