ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c            SBDART version 1.4
c            
c            Created on Thu Jul 13 16:43:35 PDT 2000
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c sbdart


c
C============================================================================
c
C PURPOSE: COMPUTE PLANE-PARALLEL RADIATIVE TRANSFER IN THE EARTH'S ATMOSPHERE
C
C INPUT/OUTPUT VARIABLES   described in rt.doc
C
c============================================================================
c
c
c The maximum number of radiation streams is set by the parameter NSTRMS
c while the maximum number of atmospheric layers is set by the parameter MXLYR
c To increase the number of streams or layers beyond the current limits, these 
c parameters must be changed in each of the subroutines in which they appear
c (i.e., do a global replace, changing nstrms=40 to nstrms=??, or changing
c mxlyr=50 to mxlry=??).  Remember to recompile after the changes are saved.
c
C  SUBROUTINE CALLS
c
c absgas          transfer gas absorption data
c absint          compute path integrals for gaseous absorption
c atms            load standard atmospheres (P,T,WH,WO)
c chkin           check input parameters
c chkprn          diagnostic print out
c cloudpar        compute cloud rt properties (single scatter albedo...)
c depthscl        compute radiative quantities for DISORT input
c disort          solve multi-stream rt equation
c getmom          compute phase functions
c kdistr          compute k-distribution fits
c locate          search on a 1-D table 
c modatm          modify atmosphere according to user inputs
c modmix          modify constituent mix of atmosphere
c modrsc          modify Rayleigh scattering strength for sensitivity studies
c openckf         open diagnostic output files (when specified)
c rayleigh        compute rayleigh scattering optical depth
c stdaer          setup standard boundary layer aerosols
c saturate        set relative humidity of cloud layer
c setfilt         select sensor filter function
c stdout0         write output file header
c stdout1         write information at each frequency
c stdout2         write frequency integrated information
c suralb          load surface albedo
c tauaero         compute aerosol optical depth
c taugas          compute absorption due to molecular species
c usraer          setup user defined boundary layer aerosols.
c zensun          compute solar zenith and azimuth given day and time
c zlayer          find the computational layer with altitude z
c
c logical units            files
c
c    11           INPUT
c
c    12           pmom direct access file
c                 (stays open throughout run)
c
c    13           usrcloud.dat, atms.dat, filter.dat, albedo.dat,
c                 aerosol.dat, aeroden.dat
c                 (each file is closed after read)
c
c    20-35        diagnostic output rtinfo.?? (?? = 00 --> 15)
c 
C=======================================================================


      program sbdart,57
      parameter (lunc=20)
      parameter (lun=11)
      parameter (nstrms=40)
      parameter (mxly=50)
      parameter (maxulv=mxly+1)
      parameter (naerz=5, ncldz=5)
      parameter (re=6371.2)
      parameter (nk=3)
      real*8 zip
      parameter (zip=-1.D0)         

      dimension z(mxly),p(mxly),t(mxly),wh(mxly),wo(mxly)

      dimension taucld(mxly),gcld(mxly),wcld(mxly)
      dimension dtaug(mxly),wreal(mxly)
      dimension dtaur(mxly),taur(mxly)
      dimension dtaua(mxly),waer(mxly),gaer(mxly),dtauab(mxly)

      dimension dtcs(mxly),dtls(mxly)
      dimension dtcv(mxly),dtlv(mxly)
      dimension tauk(mxly,nk),dtauk(mxly,nk),wtk(mxly,nk),
     &          work(mxly,nk)

      dimension umu(nstrms),utau(mxly)

      parameter (nrstrm=nstrms)

      dimension uzen(nrstrm),phi(nrstrm)
      dimension sc(4)

      parameter (mxq=63)
      
      dimension uu(mxq,mxly)
      parameter (nta=9)

      parameter (naerw=47)
      parameter (naerwm=naerw*nstrms)
      
      dimension wlbaer(naerw),qbaer(naerw),wbaer(naerw),gbaer(naerw)
      dimension pmaer(naerwm)

      dimension airmass(mxly),gasabs(nta,mxly),trnsgas(nta,mxly)

      dimension fdird(maxulv),flxdn(maxulv),flxup(maxulv),
     &   dfdt(maxulv),uavg(maxulv),uur(nrstrm,maxulv,nrstrm),
     &   u0u(nrstrm,maxulv),albmed(nstrms),trnmed(nstrms),hl(0:nstrms)
      dimension zout(2)
      dimension zcloud(ncldz),tcloud(ncldz),lcld(ncldz)
      dimension zaer(naerz),taerst(naerz),jaer(naerz),laer(naerz)
      dimension zbaer(mxly),dbaer(mxly)
      dimension pmom(0:nstrms,mxly),temper(0:mxly)

      parameter (npa=200)        ! phase function arrays for high 
      dimension sphase(npa)      ! angular resolution scattering
      dimension sangle(npa)      ! 

      logical prnt(7),deltam,lamber,onlyfl,azmavg,usrtau,usrang,plank

      real lwp(ncldz),nre(ncldz)

      logical prnton
      character*127 header

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

      data pi/3.1415926536/,wvnmlo/0.0/,wvnmhi/0.0/,accur/0.0/

      data  nz, idatm,isat,nf,iaer,isalb,ickp,iout,imom
     &    / 33,   4,   0,  2,   0,    0,   0,  10,   0 /
      data wlinf,   wlsup,   wlinc
     &     / .55,     .55,     0./
      data   vis,  rhaer, albcon,   sza, csza, rhcld, amix
     &     / 23., zip,      0.,    0.,  zip,   zip, 0. /

      data  tbaer  /zip/
      data  wlbaer /naerw*zip/
      data  qbaer  /naerw*zip/
      data  wbaer  /naerw*zip/
      data  gbaer  /naerw*zip/
      data  dbaer  /mxly*zip/
      data  zbaer  /mxly*zip/
      data  abaer  /zip/
      data  pmaer  /naerwm*zip/

      data sc/1.,0.,0.,0./
      
      data iday,alat,alon,time      /0,-64.767,-64.067,16./
      data solfac                   /1.0/

      data sclh2o,  pbar,    uw,   uo3, o3trp,  ztrp
     &     /  zip,   zip,   zip,   zip,   zip,   0./
      data    xn2,   xo2,  xco2,  xch4,  xn2o
     &     /  zip,   zip,   zip,   zip,   zip/
      data    xco,  xno2,  xso2,  xnh3,   xno,  xhno3
     &     /  zip,   zip,   zip,   zip,   zip,  zip/

      data xo4/1.0/

      data nosct/0/
      data xrsc/1.0/
      data zpres/zip/
c             
      data  zcloud,  tcloud,     nre,     lcld
     &     /ncldz*0, ncldz*0.,  ncldz*8., ncldz*0/
      data zaer,      taerst,   jaer,    laer
     &     /naerz*0, naerz*0., naerz*0, naerz*0/
      data prnt/ 7*.false./
      data  deltam,  lamber, onlyfl, azmavg, usrang, usrtau
     &   /  .true.,   .true., .true., .true., .false., .false./
      
      data ntau,numu,nzen,nphi/4*0/
      data umu/nstrms*0./ 
      data uzen/nrstrm*zip/
      data phi/nrstrm*zip/
      data nothrm,ibcnd,phi0,fisot,temis/-1,0,0.,0.,0./
      data btemp,ttemp /zip,zip/
      data zout/0.,100./ 
      data kdist/3/
      data negprn/0/

      data zgrid1,zgrid2,ngrid          /0.,30.,mxly/
      data krhclr /0/

      data ipth /1/

      data maxcmu/nstrms/
      data maxumu,maxphi/2*nrstrm/
      data nstr/4/

      data lwp/ncldz*0./

      data sphase/npa*zip/
      data sangle/npa*zip/

c     wavenumber resolution in CM-1
      data wnres/20./

      data ibeam/0/
c-----------------------------------------------------------------------

      namelist /input/

     &     idatm, amix, isat, wlinf, wlsup, wlinc, sza, csza, solfac,
     &     nf, iday, time, alat, alon, zpres, pbar, sclh2o, uw, uo3,
     &     o3trp, ztrp, xrsc, xn2, xo2, xco2, xch4, xn2o, xco, xno2,
     &     xso2, xnh3, xno, xhno3, xo4, isalb, albcon, sc, zcloud,
     &     tcloud, lwp, nre, rhcld, krhclr, jaer, zaer, taerst, iaer,
     &     vis, rhaer, tbaer, wlbaer, qbaer, abaer, wbaer, gbaer, pmaer,
     &     zbaer, dbaer, nothrm, nosct, kdist, zgrid1, zgrid2, ngrid,
     &     ickp, zout, iout, deltam, lamber, ibcnd, phi0, prnt, ipth,
     &     fisot, temis, nstr, nzen, uzen, nphi, phi, imom, negprn,
     &     sphase, sangle, ttemp, btemp, ibeam


      namelist /dinput/
     &     deltam,lamber,ibcnd,phi0,prnt,ipth,fisot,temis,nstr,
     &     nzen,uzen,nphi,phi,ttemp,btemp,imom


c-----------------------------------------------------------------------
c*****
      iprint(kkk,nnn)=mod((kkk/2**nnn),2)
      prnton(kkk,nnn)=iprint(kkk,nnn) .eq. 1
      lunit(kkk,nnn)=iprint(kkk,nnn)*(lunc+nnn)
c*****
      
      open(unit=11,file='INPUT',status='old',iostat=ios)
      if(ios.eq.0) then
        read(11,input,end=1)
        read(11,dinput,end=2)
        goto 2
 1      continue
        stop 'error: namelist block $INPUT not found'
 2      continue
      else
        write(*,input)
        stop
      endif
      zipp=zip

c assign user radiance angles

      if(nphi.gt.1) then
        if(numset(zipp,phi,nrstrm).ne.2)
     &       write(*,'(a)') 'Error in MAIN -- ' //
     &       'must specify exactly 2 values of phi when nphi is set'
        if(nphi.gt.nrstrm)
     &       write(*,'(a)') 'Error in Main -- ' //
     &       'specified nphi larger than nrstrm'
        p1=min(phi(1),phi(2))
        p2=max(phi(1),phi(2))
        do i=1,nphi
          phi(i)=p1+(i-1)*(p2-p1)/float(nphi-1)
        enddo
      else
        nphi=numset(zipp,phi,nrstrm)
      endif

      if(nzen.ne.0) then
        if(numset(zipp,uzen,nrstrm).ne.2)
     &       write(*,'(a)') 'Error in MAIN -- ' //
     &       'must specify exactly 2 values of uzen when nzen is set'
        if(nzen.gt.nrstrm)
     &       write(*,'(a)') 'Error in Main -- ' //
     &       'specified nzen larger than nrstrm'
        z1=min(uzen(1),uzen(2))
        z2=max(uzen(1),uzen(2))
        do i=1,nzen
          uzen(i)=z1+(i-1)*(z2-z1)/float(nzen-1)
        enddo
      else
        nzen=numset(zipp,uzen,nrstrm)
      endif

c check input parameters

      call chkin(isat,nf,wlinf,wlsup,solfac,isalb,albcon,
     &     sc,idatm,zpres,pbar,zcloud,tcloud,nre,jaer,zaer,
     &     taerst,iaer,zout,iout,prnt,lwp,nphi,nzen)

      if(iout.le.0) stop

      call openckf(ickp,lunc)

c solar spectrum, atmosphric profile, satellite sensor
c filtering, surface albedo and aerosol parameters

      if(iday.ne.0) then
        call zensun(abs(iday),time,alat,alon,sza,sazm,solfac)
        if(iday.lt.0) then 
          write(*,'(i5,6f9.3)') abs(iday),time,alat,alon,sza,sazm,solfac
          stop
        endif
      elseif(csza.ge.0.) then 
        sza=acos(csza)*180/pi
      endif
      if(abs(sza-90).lt..01) sza=95.

      nosun=0
      if(sza.ge.90.) nosun=1 

      call atms(idatm,amix,nz,z,p,t,wh,wo)

      if(zgrid1.ne.0.) call zgrid(nz,z,p,t,wh,wo,zgrid1,zgrid2,ngrid)

      if(zpres.gt.zip)  then
        call locate(z,nz,zpres,j)
        fj=(zpres-z(j))/(        pbar=p(j)*(      endif

      call modatm(nz,sclh2o,uw,uo3,o3trp,ztrp,pbar,z,p,wh,wo)
      call modmix(xn2,xo2,xco2,xch4,xn2o,xco,xno2,xso2,xnh3,xno,xhno3)

c array dimensions

      nstr1=nstr

      mcldz=max(numset(0.,zcloud,ncldz),
     &          numset(0.,tcloud,ncldz),
     &          numset(0.,lwp,ncldz))

      call zlayer(nz,z,mcldz,zcloud,lcld)
      call zlayer(nz,z,numset(0.,taerst,naerz),zaer,laer)
      lu1=lunit(ickp,1)
      if(rhcld.ge.0) then
        if(krhclr.eq.1) then 
          call satcloud(nz,ncldz,lcld,t,rhcld,wh,lu1)
        else
          call saturate(nz,ncldz,lcld,z,t,rhcld,wh,lu1)
        endif
      endif

      call absint(uu,nz,z,p,t,wh,wo,lunit(ickp,0))
      if(xrsc.ne.1.0) call modrsc(xrsc)

      temper(0)=t(nz)
      do 5 j=1,nz
        temper(j)=t(nz+1-j)
 5    continue
      if(btemp.lt.0.) btemp=temper(nz)
      if(ttemp.lt.0.) ttemp=temper(0)

c find nearest computational levels to given output altitudes

      call nearest(z,nz,zout(1),nbot)
      call nearest(z,nz,zout(2),ntop)
c      write(*,'(4(10f10.3)/)') z
c      write(*,*) 'ntop nbot ',ntop,nbot
      ntop=nz-ntop+1
      nbot=nz-nbot+2

c      write(*,*) 'ntop nbot ',ntop,nbot

c set up for radiance calculation

      if(nzen.ne.0) then
        
        mpa=numset(zip,sangle,npa) ! number of angles at which aerosol
                                   ! phase function is specified
        
        if(mpa.gt.0.and.tcloud(1).ne.0.) then
          write(*,*) 'Error -- cloud optical depth cannot be ',
     &         ' treated in high precision radiance mode'
          stop
        endif
        if(mpa.gt.0) then
          if(numset(zip,sphase,npa).eq.0) then
            g=gbaer(1)
c            write(*,'(2a11)') 'sangle','sphase'
            do i=1,mpa
              am=cos(sangle(i)*pi/180)
              sphase(i)=(1.-g**2)/(1.+g**2-2.*g*am)**1.5
c              write(*,'(1p,2e11.3)') sangle(i),sphase(i)
            enddo
          endif
        endif

        usrang=.true.
        onlyfl=.false.
        numu=nzen

        do 20 j=1,numu
          umu(j)=min(1.,max(cos(uzen(numu+1-j)*pi/180.),-1.))
          if(umu(j).eq.0.) then
            if(j.eq.numu) then 
              umu(j)=-.0001
            else
              umu(j)=.0001
            endif
          endif
 20     continue

        if(nphi.ne.0) then
          azmavg=.false.
        else
          nphi=1
          phi(1)=0.
        endif
      endif
        
      if(prnton(ickp,0)) 
     &   call chkprn(lunc,iday,time,alat,alon,sza,solfac,pbar,sclh2o,
     &               uw,uo3,o3trp,tcloud,zcloud,taerst,zaer,laer,
     &               zout,ntop,nbot,nz,z,p,t,wh,wo,uu)
      

      call suralb(isalb,albcon,sc)

c load extinction, absorption, and asymmetry parameters for boundary
c layer aerosols, i.e., either rural, urban, oceanic or tropospheric aerosols
c and performs interpolation over relative humidity
      
      if(rhaer.lt.0.) then
        rh=relhum(      else
        rh=rhaer
      endif

c set vertical profile 

      nzbaer=numset(zipp,zbaer,mxly)
      ndbaer=numset(zipp,dbaer,mxly)
      if(ndbaer.gt.0) then
        if(nzbaer.eq.1) then
          write(*,*) 'Error -- only one value of zbaer set'
          stop
        endif
        if(nzbaer.ne.ndbaer.and.nzbaer.gt.1) then
          write(*,*) 'Error -- number of elements must match:'
          write(*,'(a,/,(1p10e11.3))') 'zbaer',(zbaer(i),i=1,nzbaer)
          write(*,'(a,/,(1p10e11.3))') 'dbaer',(dbaer(i),i=1,ndbaer)
          stop
        endif
        call uaerden(max(nzbaer,ndbaer),zbaer,dbaer)
      endif

      if(iaer.eq.5) then
        call usraer(wlbaer,qbaer,wbaer,gbaer,pmaer,abaer,q55,imoma)
        if(tbaer.eq.zip) tbaer=q55
      else
        call stdaer(iaer,rh)
      endif

      amu0=cos(sza*pi/180.0)

c load spectral response function 

      call setfilt(isat,wlinf,wlsup,wlinc,wl1,wl2,nwl)

      call stdout0(iout,nwl,nz)
      call zeroit(trnsgas,nta*mxly)

c beginning of wavelength looping

      do 100 il=0,nwl-1
        if(wlinc.gt.1) then
          wl=(nwl-1)/(float(il)/wl2+float(nwl-1-il)/wl1)
          dw1=wl-(nwl-1)/(float(il-1)/wl2+float(nwl-il)/wl1)
          dw2=(nwl-1)/(float(il+1)/wl2+float(nwl-2-il)/wl1)-wl
        elseif(wlinc.lt.0.) then
          wr=wl2/wl1
          wl=wl1*wr**(float(il)/(nwl-1))
          dw1=wl-wl1*wr**(float(il-1)/(nwl-1))
          dw2=wl1*wr**(float(il+1)/(nwl-1))-wl
        else
          wl=wl1+il*wlinc
          dw1=wlinc
          dw2=wlinc
        endif

        if(nothrm.lt.0) then 
          plank=wl.gt.2.
        else
          plank=nothrm.eq.0
        endif

        if(il.eq.0.and.il.ne.nwl-1)     dw1=0.
        if(il.eq.nwl-1.and.il.ne.0)     dw2=0.
        dw=.5*(dw1+dw2)

        wvnm=(10000./wl)
        dwvn=.0001*wvnm
        wvnmlo=wvnm-dwvn
        wvnmhi=wvnm+dwvn
        
c input flxin as solar energy (W/m2/um) or an arbitrary constant 

        ff=filter(wl)

        if(nf.eq.0) then
          if(amu0.ne.0.) flxin=pi/amu0
        else
          flxin=solfac*solirr(wl,nf)
        endif
        
        if(iout.eq.3) ff=ff*flxin

        if(nosct.eq.1.or.nosun.eq.1) then
           flxin=0.
           amu0=1.
        endif

        rsfc=min(salbedo(wl),1.)


c calculate cloud optical depth and scattering parameters


          call taucloud(nz,ncldz,wl,lcld,lwp,tcloud,
     &             nre,wcld,gcld,taucld,imom,pmom,lu1)

        !write(*,'(a/1p(10e11.3))') 'wl',wl
        !write(*,'(a/1p(10e11.3))') 'z',(        !write(*,'(a/1p(10e11.3))') 'wcld',(wcld(jj),jj=1,nz)
        !write(*,'(a/1p(10e11.3))') 'gcld',(gcld(jj),jj=1,nz)
        !write(*,'(a/1p(10e11.3))') 'taucld',(taucld(jj),jj=1,nz)
        !write(*,'(a,1p2e11.3)') 'wl,total ',wl,total(taucld,1,nz)

c boundary layer and stratospheric aerosols

        call tauaero(wl,nz,z,vis,tbaer,naerz,taerst,laer,jaer,
     &       waer,gaer,dtaua,dtauab,imoma)

        !write(*,'(a/1p(10e11.3))') 'wl,tbaer',wl,tbear
        !write(*,'(a/1p(10e11.3))') 'z',(        !write(*,'(a/1p(10e11.3))') 'dtaua',(dtaua(jj),jj=1,nz)
        !write(*,'(a,1p2e11.3)') 'wl,tauaer ',wl,total(dtaua,1,nz)
        
c calculate optical thickness of gaseous absorption 
c and rayleigh scattering (taur, dtaur) 

        call rayleigh(wl,uu,nz,taur,dtaur)


c TAUGAS returns -ln(molecular_transmission)

        if(xo4.ne.1.0) call o4cfac(xo4)

        call taugas(wl,uu,amu0,nz,z,dtcs,dtls,lunit(ickp,2))
        call taugas(wl,uu,1.,nz,z,dtcv,dtlv,0)
        call absgas(airmass,gasabs)

        if(iout.eq.2) then
          kdist=0
          tauraytot=total(dtaur,1,nz)
          tauaertot=total(dtaua,1,nz)
        endif

c find number of passes, set up esft parameters for this wavelenth
        
        if(kdist.eq.0.or.total(dtlv,1,nz).lt..01) then
          kd=0
          npass=1
        else
          kd=kdist
          call kdistr(nz,mxly,mxq,uu,dtauk,tauk,wtk,work)
          
          test=max(tauk(nz,1),tauk(nz,2),tauk(nz,3))
          if(test.lt.0.01) then
            npass=1
          else
            npass=nk
          endif            
c            write(*,*) 'test = ',test,'  npass = ',npass
        endif
        
c       loop through k-distribution terms, calculate total optical
c       thickness (gases, rayleigh and clouds)
        
        
        do 90 k=1,npass
          
c           write(*,*) '0. maxulv=',maxulv
          

c skip call to DISORT if regular output not requested

          if(.not. (iout.eq.2.or.iout.eq.3)) then 

 60         continue
            
            call depthscl(kd,nz,k,npass,nstr,wl,amu0,z,dtls,
     &           dtcs,dtcv,dtlv,dtaur,dtaua,dtauab,taucld,dtauk,
     &           wtk,wcld,waer,gcld,gaer,dtaug,wreal,pmom,wt,
     &           imom,imoma,lunit(ickp,3))

            if(npass.gt.1) ipth=0
            
c            write(*,'(f10.3,3i4,1p7e12.4)')
c     &       wl,kd,k,npass,total(dtaur,1,nz),total(dtaua,1,nz),
c     &       total(taucld,1,nz),total(dtcv,1,nz),total(dtlv,1,nz),
c     &       total(dtaug,1,nz),dtaug(nz)

            call disort(nz,dtaug,wreal,pmom,temper,wvnmlo,wvnmhi,
     &           usrtau,ntau,utau,nstr,usrang,numu,umu,nphi,phi,ibcnd,
     &           flxin,amu0,phi0,fisot,lamber,rsfc,hl,btemp,ttemp,
     &           temis,deltam,plank,onlyfl,accur,prnt,header,mxly,
     &           maxulv,maxumu,maxcmu,maxphi,fdird,flxdn,flxup,dfdt,
     &           uavg,uur,u0u,albmed,trnmed)

c    add in atmospheric single scattering contribution when aerosol
c    phase function (not moments) is specified.

c      call zeroit(uur,nstrms*maxulv*nstrms)

            if(ibeam.eq.2) call beamsct(nphi,phi,numu,umu,phi0,amu0,
     &           nz,dtaug,wreal,mpa,sangle,sphase,flxin,dtaur,iout,
     &           nbot,ntop,uur)
c      stop
            numu=nzen

c            
c    change nstr if a negative radiance is found within 5 layers
c    of the output layer.
c
            
            rmin=0.
            if(iout.eq.20) then
              call radrng(uur,nzen,nz+1,nphi,ntop,5,radmin,radmax)
            elseif(iout.eq.21) then 
              call radrng(uur,nzen,nz+1,nphi,nbot,5,radmin,radmax)
            elseif(iout.eq.22) then 
              call radrng(uur,nzen,nz+1,nphi,ntop,nz+1,radmin,radmax)
            elseif(iout.eq.23) then 
              call radrng(uur,nzen,nz+1,nphi,ntop,5,radmin,radmax)
              call radrng(uur,nzen,nz+1,nphi,nbot,5,radmn,radmx)
              radmin=min(radmin,radmn)
              radmax=max(radmax,radmx)
            endif
          
            if(negprn.eq.1.and.rmin.lt.0.) then 
              write(*,*) 'number of streams:',nstr
              write(*,*) 'min,max radiance: ',radmin,radmax
              do iz=1,nz+1,nz
                write(*,'(x,a,i4)') 'iz = ',iz
                do izen=1,nzen
                  write(*,'(f10.3,1p11e11.3)') uzen(izen),
     &                 (uur(izen,iz,iphi),iphi=1,min(nphi,11))
                enddo
              enddo
            endif

             if(radmin.lt.-1.e-8*abs(radmax)) then
               nstr=nstr-2
               if(nstr.lt.max(4,nstr1-2)) nstr=min(nstrms,nstr1+2)
               if(nstr.eq.nstr1) then
                 write(*,*)'%%% Could not repair negative radiance'
                 stop
               endif
               go to 60
             else
               nstr=nstr1
             endif
            
c        write(*,*) 'numu: ',numu
c        write(*,*) 'umu:  ',umu
c        write(*,*) 'nphi: ',nphi
c        write(*,*) 'phi:  ',phi
c        write(*,*) 'uur:'
c        do ii=1,nzen
c          write(*,'(1p10e11.3)') (uur(ii,nz+1,kk),kk=1,nphi)
c        enddo
            
            if(prnton(ickp,4)) then
              write(lunc+4,'(a,i5,e11.3)') 'k, wt ',k,wt
              write(lunc+4,'(3x,5a11)')
     &          'dfdt','dflxup','fdird','flxnd','uavg'
              write(lunc+4,'(i3,1p5e11.3)')
     &          (i,dfdt(i),flxup(i),fdird(i),flxdn(i),uavg(i),i=1,nz+1)
            endif
            
          endif

          call stdout1(nz,z,ntop,nbot,iout,wl,dw,wt,
     &         fdird,flxdn,flxup,ff,nphi,nzen,sazm,phi,uzen,uur,
     &         gasabs,trnsgas,tauraytot,tauaertot,taucldtot,k,npass)

 90     continue

 100  continue
      
      call stdout2(nz,iout,wl1,wl2,nphi,nzen,sazm,phi,uzen,z,
     &     airmass,trnsgas)
      
 1000 format(/a/(1p10e11.3))

      end