c-----------------------------------------------------------------------


      subroutine cloudqwp(wl,re,q,w,pmom) 2
c
c input:
c
c   wl      wavelength in microns
c
c   re      effective cloud drop radius (microns)
c
c output:
c
c   q       mean extinction efficiency (unitless)
c
c   w       single scatter albedo      (unitless)
c
c   pmom    Legendre expansion coeficients to the scattering phase 
c           phase function             (unitless)
c     
c
c   This routine reads a direct access file in the current working
c   directory.  Each data record in this file consists of
c   qbar,omega,pmom(1:nmom).  where qbar is the extinction efficiency,
c   omega is the single scattering albedo and pmom are the legendre
c   coefficients.  To maintain portability and reduce size of the file
c   each quantity is scaled to integer values. The physical variable
c   is recovered by inverting the scaling process using the realizable
c   range of each quantity.  Hence,
c
c        physical_value = pmin+(pmax-pmin)*float(i-imin)/(imax-imin), 
c
c    where, (pmin,pmax)=(0,5) for qbar, (0,1) for omega and (-1,1)
c    for pmom, and imin,imax is given in the first record.
c
c    The first record of the file contains an ascii string containing
c    the values: irec1, num_wl, num_re, imin, imax, written with format
c    '(3i10,2i15)', where irec1 is the locatation of the first data
c    record num_wl is the total number of wavelength points, num_re is
c    the number of effective radii, imin/imax is the minimum/maximum
c    integer values used in the scaling.  The second record contains the
c    values, wl_min,wl_max,re_min,re_max written with format '(4f15.4),
c    which provide the min and max physical values of wavelength and
c    effective radius.  The next several records contain ascii header
c    info describing the contents of the file. This routine ignors
c    all the ASCII data after the first two records.
c
c

cccc
c 
c the record length written by wri_pmom is currently help at 2+20=22
c  
cccc

      

      parameter (lunda=12)
      parameter (nstrms=40)
      parameter (nstrpm=20)
      parameter (lenrec=nstrpm+2)
      parameter (mone=-1.)
      dimension pm(lenrec)
      dimension pmom(*)

      integer m(lenrec)

      character*(lenrec*4) string

      real vmin(lenrec),vmax(lenrec)
      logical first
      
      data vmin/0.,0.,nstrpm*mone/
      data vmax/5.,1.,nstrpm*1./

      data first/.true./

      data istart,num_wl,num_re/0, 0, 0/

      data mmin,mmax /0,0/

      data wl_min, wl_max, re_min, re_max, wfac, rfac/
     &       0.,     0.,     0.,     0.,    0.,   0. /

      if(first) then 
        open(lunda,status='old', form='unformatted',
     &     access='direct',file='pmom.dat',recl=nstrpm+2)
        read(lunda,rec=1) string

        read(string,'(3i10,2i15)') istart,num_wl,num_re,mmin,mmax

        read(lunda,rec=2) string
        read(string,'(4f15.4)') wl_min,wl_max,re_min,re_max
        if(num_re.ne.1) rfac=(num_re-1)/log(re_max/re_min)
        if(num_wl.ne.1) wfac=(num_wl-1)/log(wl_max/wl_min)
c        write(*,*) 'mmin, mmax',mmin,mmax
        first=.false.
      endif

      fre=rfac*log(re/re_min)
      fwl=wfac*log(wl/wl_min)

      ire=int(fre)
      iwl=int(fwl)

      if(ire.lt.0.or.ire.gt.num_re-1) then
        write(*,*) 'error in cloudqwp --- re value out of bounds'
        write(*,*) 're,re_min,re_max: ',re,re_min,re_max
        stop
      endif

      if(iwl.lt.0.or.iwl.gt.num_wl-1) then
        write(*,*) 'error in cloudqwp --- wl value out of bounds'
        write(*,*) 'wl,wl_min,wl_max: ',wl,wl_min,wl_max
        stop
      endif

      wre=fre-ire
      wwl=fwl-iwl

      irec=istart+ire*num_wl+iwl

c      write(*,*) 'iwl,wwl = ',iwl,wwl
c      write(*,*) 'ire,wre = ',ire,wre

      f=(1.-wre)*(1.-wwl)
      if(f.ne.0.) then 
        read(lunda,rec=irec) (        do i=1,lenrec
          p=vmin(i)+(vmax(i)-vmin(i))*float(          pm(i)=f*p
c          if(i.eq.1) write(*,*) '1%', f,vmin(1),vmax(1),m(1),mmin,mmax,p
        enddo
      endif

      f=(1.-wre)*wwl
      if(f.ne.0.) then 
        read(lunda,rec=irec+1) (        do i=1,lenrec
          p=vmin(i)+(vmax(i)-vmin(i))*float(          pm(i)=f*p+pm(i)
c          if(i.eq.1) write(*,*) '2%', f,vmin(1),vmax(1),m(1),mmin,mmax,p
        enddo
      endif

      f=wre*(1.-wwl)
      if(f.ne.0.) then 
        read(lunda,rec=irec+num_wl) (        do i=1,lenrec
          p=vmin(i)+(vmax(i)-vmin(i))*float(          pm(i)=f*p+pm(i)
c          if(i.eq.1) write(*,*) '3%', f,vmin(1),vmax(1),m(1),mmin,mmax,p
        enddo
      endif

      f=wre*wwl
      if(f.ne.0.) then 
        read(lunda,rec=irec+num_wl+1) (        do i=1,lenrec
          p=vmin(i)+(vmax(i)-vmin(i))*float(          pm(i)=f*p+pm(i)
c          if(i.eq.1) write(*,*) '4%', f,vmin(1),vmax(1),m(1),mmin,mmax,p
        enddo
      endif

c NOTE: the zeroeth legendre moment (which is always equal to one)
c       is not returned.  The first legendre moment, pmom(1), is
c       equal to the asymmetry factor, g.

      q=pm(1)
      w=pm(2)

      do i=1,nstrpm
        pmom(i)=pm(2+i)
      enddo
c
c assume Henyey-Greenstein extrapolation for higher moments
c
      do i=nstrpm+1,nstrms
        pmom(i)=pmom(nstrpm)*(float(nstrpm)/float(i))**i
      enddo

      return
      end