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

      subroutine saturate(nz,ncldz,lcld,z,t,rhcld,wh,lunc) 1,1
      parameter (tzero=273.15)
      dimension t(*),z(*),wh(*),lcld(*)
c
c purpose: modify the watervapor density inside clouds to have a relative
c          humidity of RHCLD.  reduce water vapor density outside
c          of cloud to maintain total water vapor path
c
c          NOTE: due to normalization procedure the final water vapor
c                density inside the cloud may not exactly match 
c                (rhcld)*(saturated water vapor density), 
c                but it should be close.
c
c INPUT:
c   nz        number of atmospheric levels
c   ncldz     size cloud layer arrays
c   lcld      cloud layer array
c   z         altitude (km) of atmospheric layers (z(nz)=surface)
c   t         temperature at altitude z
c   rhcld     relative humidity of cloud layer, used to specify
c             water vapor density.  
c                                   
c
c h2osat is the mass density (g/m3) of water vapor at 100% saturation 
c (source: handbook of chemistry and physics, h2o vapor pressure table d-112,
c assuming density related to pressure by ideal gas law)
c
            satden(a)=a*exp(18.916758-a*(14.845878+a*2.4918766))
c
c  adjust cloud layer, don't adjust clear layers
c
c  adjust cloud layer and adjust clear layer to compensate for increased wvp
c
      if(lunc.gt.0) write(lunc,'(a/1p(10e11.3))') 'wh_in  ',
     &           (wh(i)/satden(tzero/t(i)),i=1,nz)

      wvp=0.
      zbot=z(1)
      do 5 i=1,nz
        if(i.eq.1) then
          ztop=.5*(        elseif(i.eq.nz) then
          ztop=z(nz)
        else
          ztop=.5*(        endif
        wvp=wvp+.1*(ztop-zbot)*wh(i)
        zbot=ztop
 5    continue
      if(wvp.eq.0.) return

      do 20 i=1,ncldz
        call levrng(ncldz,lcld,i,lbot,ltop)
        if(lbot.ne.0) then
          do 10 j=ltop,lbot
            jj=nz-j+1
            wh(jj)=-rhcld*satden(tzero/t(jj))
 10       continue
        endif
 20   continue

      wvpclr=0.
      wvpcld=0.
      zbot=z(1)
      do 30 i=1,nz-1
        if(i.eq.1) then
          ztop=.5*(        elseif(i.eq.nz) then
          ztop=z(nz)
        else
          ztop=.5*(        endif
        if(wh(i).gt.0.) wvpclr=wvpclr+.1*(ztop-zbot)*wh(i)
        if(wh(i).lt.0.) wvpcld=wvpcld-.1*(ztop-zbot)*wh(i)
        zbot=ztop
 30   continue

      if(wvpcld.eq.0) return

      if(wvpclr.eq.0) then
        cldfac=wvp/wvpcld
        clrfac=0.
      else
        clrfac=(wvp-wvpcld)/wvpclr
        cldfac=1.
        if(clrfac.lt.0) then
          clrfac=0.
          cldfac=wvp/wvpcld
        endif
      endif

      do 40 i=1,nz
        if(wh(i).lt.0) then
          wh(i)=-cldfac*wh(i)
        else
          wh(i)=clrfac*wh(i)
        endif
 40   continue

      wvpn=0.
      do 50 i=1,nz-1
        dz=(        den1=wh(i)
        den2=wh(i+1)
        tst1=abs(den1-den2)
        tst2=min(den1,den2)
        if(tst1.le..001*den1 .or. tst2.eq.0.) then
          du=.5*dz*(den1+den2)
        else
          du=dz*(den1-den2)/log(den1/den2)
        endif
        wvpn=wvpn+.1*du
 50   continue

      do 60 i=1,nz
        wh(i)=wh(i)*wvp/wvpn
 60   continue

      if(lunc.gt.0) then
      
        write(lunc,'(4a10/4f10.3)') 'wvp','wvpn','wvpclr','wvpcld',
     &       wvp,wvpn,wvpclr,wvpcld
        write(lunc,'(a/1p(10e11.3))') 'rh_out  ',
     &           (wh(i)/satden(tzero/t(i)),i=1,nz)
        write(lunc,'(a/1p(10e11.3))') 'wh_out  ',
     &           (wh(i),i=1,nz)
      endif      
      return
      end