c=======================================================================















c file:                  cloudpar.f
c
c external routines:     cloudpar, cloudqwp
c
c internal routines:     clddata (block data)
c
c internal common block: cloudblk
c
c a direct access file "pmom.dat" is required by cloudqwp
C=======================================================================


      SUBROUTINE cloudpar(wl,re,qc,wc,gc) 3
c
c input:
c   wl      wavelength in microns
c   re      effective cloud drop radius (microns)
c output:
c   qc      mean extinction efficiency (unitless)
c   wc      single scatter albedo      (unitless)
c   gc      asymmetry factor           (unitless)
c     

      parameter (mxwv=400, mre=14)

c the parameters contained in this subroutine were computed using a mie
c scattering code provided by Paul Stackhouse.  The radius integrations
c were performed over a modified gamma distribution
c f(x)=(x/xm)^(p-1)*exp(-x/xm) where p=7.

c NOTE:

c re < 0 selects scattering parameters for cirrus cloud ice particles
c with a effective radius of 106um.  This data is stored in array
c elements (*,14) and is valid in the the wavelength range 0.29-333.33 um. 
c When re<0 the absolute value of re is used as a multiplier of the single
c scattering co-albedo hence re=-2 doubles the amount of absorption per 
c extinction event. 

c

      parameter (wlmin=0.29, wlmax=333.33)
      parameter (eps=.000001)
      common /cloudblk/ qq(mxwv,mre),ww(mxwv,mre),gg(mxwv,mre)
                      
      data icall,wmin,wmax,wstep/0,0.,0.,0./

      if (icall .eq. 0) then
        wmin=log(wlmin)
        wmax=log(wlmax)
        wstep=(wmax-wmin)/(mxwv-1)
        icall=1
      endif

      fw=1+(log(wl)-wmin)/wstep
      fw=min(max(fw,1.),float(mxwv)-eps)
      iw=int(fw)
      fw=fw-iw

      if(re .lt. 0.) then 
        ir=14
        qc=qq(iw  ,ir  )*(1.-fw) + qq(iw+1,ir  )*fw
        wc=ww(iw  ,ir  )*(1.-fw) + ww(iw+1,ir  )*fw
        gc=gg(iw  ,ir  )*(1.-fw) + gg(iw+1,ir  )*fw
        wc=1.-abs(re)*(1.-wc)
        wc=max(wc,eps)
      else
        fr=1.+((log(re))/log(2.)-1.)*2
        fr=min(max(fr,1.),float(mre)-eps)
        ir=int(fr)
        fr=fr-ir
        qc=qq(iw  ,ir  )*(1.-fw)*(1.-fr) + qq(iw+1,ir  )*fw*(1.-fr)+
     &     qq(iw  ,ir+1)*(1.-fw)*fr      + qq(iw+1,ir+1)*fw*fr

        wc=ww(iw  ,ir  )*(1.-fw)*(1.-fr) + ww(iw+1,ir  )*fw*(1.-fr)+
     &     ww(iw  ,ir+1)*(1.-fw)*fr      + ww(iw+1,ir+1)*fw*fr

        gc=gg(iw  ,ir  )*(1.-fw)*(1.-fr) + gg(iw+1,ir  )*fw*(1.-fr)+
     &     gg(iw  ,ir+1)*(1.-fw)*fr      + gg(iw+1,ir+1)*fw*fr

      endif
      return
      end