c
c
c diagnostic main program
c
c        write(*,*)'---------'
c        do 1 t=.1,.5,.1
c          e1=expint(1,t)
c          ex=(e1+log(t)+0.5772156649)/t
c          write(*,'(1p7e15.7)') t,e1,ex
c   1    continue
c        write(*,*)'---------'
c        do 5 t=.5,1.,.1
c          e1=expint(1,t)
c          write(*,'(1p7e15.7)') t,e1
c   5    continue
c        write(*,*)'---------'
c        do 10 t=1.,8.
c          e1=expint(1,t)
c   10   continue
c        write(*,*)'---------'
c        end

c file:                  atms.f
c
c external routines:     atms,modatm,taucloud,saturate,satcloud
c
c internal routines      midsum,midwin,subsum,subwin,tropic,us62,useratm
c                        zgrid
c
c logical units:         unit 13, used to read user atmosphere from atms.dat

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

      subroutine taucloud(nz,ncldz,wl,lcld,lwp,tcloud, 1,6
     &              nre,wcld,gcld,taucld,imom,pmom,lunc)
c
c purpose:  compute cloud radiometric properties
c
c input: 
c
c  nz      number of atmospheric layers
c  ncldz   maximum number of cloud layers
c  wl      wavelength
c  lcld    cloud layer indecies
c  lwp     liquid water path of each cloud layer
c  tcloud  optical depth of each cloud layer
c  nre     effective radius within each cloud layer
c  lunc    logical unit for diagnostic print
c
c output:
c
c  wcld    cloud single scatter albedo at each vertical grid point
c  gcld    asymmetry factor at each vertical grid point
c  taucld  cloud optical thickness at each vertical grid point
c
c                            the effective radius of ice particles:
      parameter (reice=106.)
c                            the bulk density of ice (g/cm3):
      parameter (rhoice=.917) 

      parameter (mxly=50)
      parameter (nstrms=40)
      dimension pmcld(nstrms)
      dimension pmom(0:nstrms,*)

      real lwp(*),tcloud(*),nre(*),lwpth
      integer lcld(*),icnt(mxly)
      dimension wcld(*),gcld(*),taucld(*)
      dimension qcld(mxly),q550(mxly)

      save q550
     
      data ifirst/1/, ipo/1/
      
c      write(*,*) '.lcld1 ',lcld1      
c      write(*,*) '.lcld2 ',lcld2      

      do 5 j=1,nz
        qcld(j)=0.
        wcld(j)=0.
        gcld(j)=0.
        taucld(j)=0.
        icnt(j)=0
        do i=0,nstrms
          pmom(i,j)=0.
        enddo
 5    continue

      if(tcloud(1).gt.0.) then
        do 20 i=1,ncldz
          call levrng(ncldz,lcld,i,lbot,ltop)
          if(lbot.eq.0) goto 20
          if(tcloud(i).eq.0..and.lwp(i).eq.0.) goto 20
          do 10 j=ltop,lbot
            if(ltop.eq.lbot) then
              reff=nre(i)
              tcld=tcloud(i)
              lwpth=lwp(i)
            else
              wt=float(j-ltop)/(lbot-ltop)
              reff=nre(i+1)*(nre(i)/nre(i+1))**wt

c           if tcloud(i+1) is zero, spread the opacity uniformly over
c           the affected layers.  if tcloud(i+1) is nonzero then the
c           opacity is distrbuted such that it linearly increases
c           (or decreases) from the bottom to the top layer.  The rate
c           of increase is controled by tcloud(i+1) which sets the ratio
c           of the the opacity in the top layer to that of the bottom
c           layer.  The total opacity over all layers adds up to
c           tcloud(i)

              if(tcloud(i+1).eq.0.) then
                tcld=tcloud(i)/(lbot-ltop+1)
              else
                tcld=2*tcloud(i)/((lbot-ltop+1)*(1.+tcloud(i+1)))
                tcld=tcld+(lbot-j)*tcld*(tcloud(i+1)-1.)/(lbot-ltop)
              endif

              if(lwp(i+1).eq.0.) then
                lwpth=lwp(i)/(lbot-ltop+1)
              else
                lwpth=2*lwp(i)/((lbot-ltop+1)*(1.+lwp(i+1)))
                lwpth=lwpth+(lbot-j)*lwpth*(lwp(i+1)-1.)/(lbot-ltop)
              endif
            endif

            if(imom.ge.1) then 
              call cloudqwp(wl,reff,qc,wc,pmcld)
              gc=pmcld(1)
              do k=1,nstrms
                pmom(k,j)=pmcld(k)+pmom(k,j)
              enddo
            else
              call cloudpar(wl,reff,qc,wc,gc)
            endif
            qcld(j)=qc+qcld(j)
            wcld(j)=wc+wcld(j)
            gcld(j)=gc+gcld(j)
            icnt(j)=1+icnt(j)
            if(tcloud(i).ne.0.) then
              if(ifirst.eq.1) then
                if(imom.ge.1) then
                  call cloudqwp(0.55,reff,qc550,wc550,pmcld)
                else
                  call cloudpar(0.55,reff,qc550,wc550,gc550)
                endif
                q550(j)=qc550
              endif
              taucld(j)=tcld*qc/q550(j)+taucld(j)
            else
              if(lwpth.ne.0.) then
                if(reff.lt.0.) then
                  taucld(j)=.75*qc*lwpth/reice/rhoice+taucld(j)
                else
                  taucld(j)=.75*qc*lwpth/reff+taucld(j)
                endif
              endif
            endif
            
 10       continue
 20     continue
      elseif(tcloud(1).lt.0.) then 

c read cloud optical depth from file usrcloud

        open(unit=13,file='usrcloud',form='formatted',status='old')
        read(13,*) nzin
        if(nzin.ne.nz) then
          write(*,*) 'error in file usrcloud --- ' //
     &         'incorrect number of layers'
          stop
        endif

        do j=1,nz
          jm=nz-j+1
          read(13,*) zz,wlc,reff
          if(j.gt.1) taucld(jm+1)=taucld(jm+1)*(zz-zb)
          if(wlc.ne.0) then 
            lwpth=1000.*wlc        ! dz factor not included
            call cloudpar(wl,reff,qc,wc,gc)
            if(reff.lt.0.) then
              taucld(jm)=.75*qc*lwpth/rhoice/reice
            else
              taucld(jm)=.75*qc*lwpth/reff
            endif
            wcld(jm)=wc
            gcld(jm)=gc
          else
            taucld(jm)=0.
            wcld(jm)=0.
            gcld(jm)=0.
          endif     
          zb=zz
        enddo
        close(13)
      endif
      
      do 30 j=1,nz
        if(icnt(j).ne.0) then
          qcld(j)=qcld(j)/icnt(j)
          wcld(j)=wcld(j)/icnt(j)
          gcld(j)=gcld(j)/icnt(j)
          if(imom.ge.1) then 
            do i=1,nstrms
              pmom(i,j)=pmom(i,j)/icnt(j)
            enddo
          endif
        endif
        
        if (lunc.gt.0 .and. ifirst.eq.1 .and. taucld(j).ne.0.) then
          if(ipo.eq.1) write(lunc,'(8a11)')
     &         'wl','layer ','taucld','qc','ssa','g'
          ipo=0
          write (lunc,'(f11.5,i11,4(f11.5),i11)')
     &         wl,j,taucld(j),qcld(j),wcld(j),gcld(j),icnt(j)
        endif
 30   continue

      ifirst=0
      return
      end