c
c------------------------------------------------------------------
c Begin file /home/sting/foster/cira/cira.f
c------------------------------------------------------------------
c
      subroutine cira(month,ludat,nlat,nvert,clat,calt,cpmb,
     +  zmtn,zmht,zmun,ier)
c
c Get NASA CIRA data from mss for tn, un, ht:
c
c On input:
c   month = integer of desired month of data (1 < month < 12)
c   ludat = logical unit on which to acquire and read mss data file   
c   nlat = 17 = number of latitudes in cira grid
c   nvert = 36 = number of vertical points in cira grid
c   
c On output:
c   clat(nlat) = latitude grid points
c   calt(nvert) = altitude grid points
c   cpmb(nvert) = pressure at vertical grid (mb)
c   zmtn(nlat,nvert) = zonal mean neutral temperature at desired month
c   zmht(nlat,nvert) = zonal mean heights (km) at desired month
c   zmun(nlat,nvert) = zonal mean zonal wind (m/s) at desired month
c   ier = 0 if all ok, < 0 otherwise (if ier < 0, this routine
c         returns without defining other output)
c
c Required subroutines:
c   tail(), lenstr() (available in sting:~foster/lib)
c
c Note: this routine reads data from mss file /FOSTER/cira/cira88.dat
c
      parameter (rlat1=-80.,dlat=10.)
      dimension zmtn(nlat,nvert),zmht(nlat,nvert),
     +  zmun(nlat,nvert),cpmb(nvert),calt(nvert),clat(nlat)
c
      character*3 chmon12(12),chmon
      character*40 mssdat,msstail
      character*80 errmsg,rec80,asncmd
      character*1 dum
      data mssdat /'/FOSTER/cira/cira88.dat                 '/
      logical isthere,found
      data chmon12 /'JAN','FEB','MAR','APR','MAY','JUN','JUL','AUG',
     +              'SEP','OCT','NOV','DEC'/
c
c Check input:
c
      ier = 0
      if (ludat.lt.1.or.ludat.gt.99.or.ludat.eq.5.or.ludat.eq.6) then
        write(6,"('>>> cira: bad ludat=',i5)") ludat
        ier = -1
        return
      endif
      if (ludat.eq.2) then
        write(6,"('>>> cira warning: ludat should not be 2 because',
     +    ' of potential conflict with NCAR gks unit -- ',
     +    '(this is non-fatal to cira) -- continuing...')") 
      endif
      if (nlat.ne.17) then
        write(6,"('>>> cira: need nlat=17: nlat=',i3)") nlat
        ier = -1
        return
      endif
      if (nvert.ne.36) then
        write(6,"('>>> cira: need nvert=36: nvert=',i3)") nvert
        ier = -1
        return
      endif
      if (month.lt.1.or.month.gt.12) then
        write(6,"('>>> cira: bad month=',i4,' (1 < month < 12)')")
     +    month
        ier = -1
        return
      endif
      chmon = chmon12(month) 
c
c Set up latitude grid:
c
      do j=1,nlat
        clat(j) = rlat1 + (j-1)*dlat
      enddo
c
c Acquire and assign mss data file:
c (msread not performed if disk file already in working dir)
c
      lenmss = lenstr(mssdat)
      call tail(mssdat,msstail)
      lentail = lenstr(msstail)
      inquire(file=msstail(1:lentail),exist=isthere)
      write(6,"(' ')")
      if (.not.isthere) then
        write(6,"('cira: msread on file ',a)") mssdat(1:lenmss)
        call msread(ier,msstail(1:lenmss),mssdat(1:lenmss),' ',
     +              ' ')
        if (ier.ne.0) then
          call mserror(errmsg)
          write(6,"('>>> cira: error obtaining mss file ',a,
     +      /' err msg=',a)") mssdat(1:lenmss),errmsg
          ier = -1
          return
        endif
      endif 
      call asnunit(ludat,'-s text -a '//msstail(1:lentail),ieras)
      if (ieras.ne.0) then
        write(6,"('>>> cira: error ',i3,' assigning disk file ',a,
     +    ' to unit ',i3)") ieras,msstail(1:lentail),ludat
        ier = -1
        return
      else
        write(6,"('cira: assigned disk file ',a,' to unit ',i3)")
     +    msstail(1:lentail),ludat
      endif
c
c Read through file:
c There are 12 months for each of 3 fields (tn, ht, un in that order)
c
      ifield = 0
      found = .false.
      do ir=1,100000
        read(ludat,"(a)",end=900) rec80
c
c Found desired month:
        if (rec80(1:3).eq.chmon) then
          found = .true.
          write(6,"('Cira found data for month ',i3)") month
          write(6,"(5x,a)") rec80(1:50)
          ifield = ifield+1
c
c Skip 4 line header between month line and data:
          do i=1,4
            read(ludat,"(a)") dum
          enddo
c
c Read data:
c (note data on file starts at top alt and goes down, but we want
c  data from bottom up, so invert indices (kk) as we read):
c
          do k=1,nvert
            kk = nvert-k+1
            if (ifield.eq.1) then
              read(ludat,"(f5.1,f8.0,f6.1,16f7.1)") calt(kk),cpmb(kk),
     +          (zmtn(j,kk),j=1,nlat)
            elseif (ifield.eq.2) then
              read(ludat,"(f5.1,f8.0,f6.1,16f7.1)") calt(kk),cpmb(kk),
     +          (zmht(j,kk),j=1,nlat)
            elseif (ifield.eq.3) then
              read(ludat,"(f5.1,f8.0,f6.1,16f7.1)") calt(kk),cpmb(kk),
     +          (zmun(j,kk),j=1,nlat)
            endif
          enddo ! k=1,nvert
c
c First values apparently missing on data file:
c (Do simple extrapolation from 2nd and 3rd latitude values)
c
          zmtn(1,1) = 2.*zmtn(2,1)-zmtn(3,1)
          zmht(1,1) = 2.*zmht(2,1)-zmht(3,1)
          zmun(1,1) = 2.*zmun(2,1)-zmun(3,1)
        endif ! rec80 == chmon
      enddo ! read loop
      goto 901
c
c EOF on mss data file:
 900  write(6,"('Cira: EOF on unit ',i3)") ludat
      goto 902
c
c 100k records on file? (should not happen):
 901  write(6,"('>>> cira: end of read loop')")
c
 902  continue
      if (.not.found) then
        write(6,"('>>> cira: could not find month ',i3)") month
        ier = -1
      endif
      return
      end