c  $Header: /usr/people/flittner/radtran/tomrad/pv2.2/src/RCS/matscn.f,v 2.22 2000/08/30 16:38:09 flittner Exp $
      program matscn
c     author:  c. seftor
c     date:    12-mar-1992
c     last modified 9/22/93  .... zia ahmad
c     purpose: to include the effect of molecular anisotropy in
c              rayleigh phase function
c
c     purpose: used in conjunction with the mateer scan code to
c              generate radiance tables for a given ozone profile
c              and for given scan and azimuth angles (modified from
c              profil).
c
c     last modified 03/07/95... dave flittner
c     purpose: compute spherical correction to the out going beam
c     when integrating after iteration process.  Set with lsphout
c     logical variable.
c
c     last modified 03/08/95... dave flittner
c     purpose: store Z1=I1/(-3/8*muo*sqrt(1-muo^2)*sqrt(1-mu^2)) and
c     Z2=I2/(3/32*(1-muo^2)*(1-mu^2)/mu) so can be used as input into
c     interpolation tables.
c
c     last modified 03/14/95...dave flittner
c     purpose: set pressure scale height used in gravity correction
c     to rayleigh scattering od.  Create new variable pscaleforg and
c     pass in common block consts.  Also use logical switch lgcorrect
c     to impliment the gravity correction to the rayleigh scattering
c     optical depth.
c
c     last modified 03/22/96...dave flittner
c     purpose: use parameter statements for the scan angles, solar
c     zenith angles and azimuth angles.
c
c     last modified 05/10/96...dave flittner
c     purpose: allow for absorption by 2 gases
c
c     last modified 05/22/96...dave flittner
c     purpose: print-out table for downward flux at sfc.
c
c     last modified 09/03/96...dave flittner
c     purpose: allow computation of the transmitted radiance at the surface.
c     This is activated by the logical switch ldown set to TRUE.
c
c     last modified 09/05/96...dave flittner
c     purpose: set switchs in subroutine rdenv
c
c     last modified 09/20/96...dave flittner
c     purpose: fix call to rdenv.  Add ipsudo=2 option for computting
c     single scattering in spherical atmosphere.  Resultant radiance is
c     sum of spherical single scatter and flat multiple scatter.
c     In addition, clean-up the declaration of the variables.
c
c     last modified 8/30/97...dave flittner
c     purpose: include actinic flux info in the flux.out file.
c
c     last modified 10/13/97...dave flittner
c     purpose: change output to iter file so to pass the refdir instead of
c     exp(-tautot), which was only for the nadir direction.
c
c     last modified: July, 1998...def
c     purpose: correct actinic flux notes and mod. for elevated sfc in the
c     input profile.  This was used in the flux intercomparison.
c
c     last modified: Sep. 2, 1998...def
c     purpose: use include for all common blocks.  Also pass ipsudo to rdenv
c     through the "in" common block.  Add routine to read the command line.
c
c     last mod. Sep. 6, 1998...def
c     purpose: combine several commons into switches.cmn. Delete ctol routine.
c
c     last mod. Sep. 25, 1998...Ed Celarier
c     purpose: add some comments and fix index when storing flux values.
c
c     last mod. Sep. 29, 1998...eac
c     purpose: mod. output statements for version given to Batluck as version 2.9.
c
c     last mod. May 3, 1999...def
c     purpose: add option to compute addition absorbing gas profile using the
c     internal, implicit neutral density profile and a constant mixing ratio.
c     Currently it is set to do only o2 or o4 absorption and is controled by
c     the lo2abs and lo4abs logicals read in the ENV file.
c
c     last mod. May 6, 1999...def
c     purpose: move call to relayr outside the wavelength loop and pass
c     pshold and hhold through the common block contrl instead of command line.
c
c     last mod. May 10, 1999...def
c     purpose: delete recl from open of units 23, 33.
c
c last mod. Nov 2, 1999...def
c purpose: use .not. instead of .eqv. in logical testing.
c
c last mod.: Mar. 14, 2000...def
c purpose: Change sb and sbp to single elements.
c
c last mod.: May 16, 2000...def
c purpose: add switch to use o3 weighted temperatures with the user
c input profile.
c
c***************************************************************************
c     
      implicit none
      include 'parameter.inc'
C
C   Local Variables
C
      integer*4 nwavl,jjprint,iter_lun,iuntasc,iuntbin,j
      integer*4 iuntcontrb
      integer*4 num,nc,ncp1,lmax,ijk,k,jlev,iter_num,ialbedo
      integer*4 nprttabaz,iter_pos,iaz,i
      real*8 p(101),t(101),tl(101),z(202)
      real*8 wavel(max_wave), tmpcf(2),cofx(4,487),h(487),ps(487)
      real*8 fracin,xpnot,thenot,f1,f2
      real*8 fsp(max_az,max_scan)
C
C   Setting up virtual memory for output parameters
C
      real*8 i0_out(max_scan, max_sza, max_wave,max_az)
      real*8 Z1_out(max_scan, max_sza, max_wave,max_az)
      real*8 Z2_out(max_scan, max_sza, max_wave,max_az)
      real*8 t_out(max_scan, max_sza, max_wave,max_az)
      real*8 sb_out(max_wave)
C new ones needed for the flux.out tables
      real*8 f0atten_out(max_sza, max_wave),gg_out(max_sza, max_wave)
      real*8 ggp_out(max_sza, max_wave),sbp_out(max_wave)
      real*8 f0_out(max_wave),tau_tots(max_wave),tau_tota(max_wave)
C
C   Functions
C
      integer char2int
      logical get_coef
      external char2int,get_coef
C
C   Common
C
      include 'switches.cmn'
      include 'kmat.cmn'
      include 'log.cmn'
      include 'inchr.cmn'
      include 'emm.cmn'
      include 'consts.cmn'
      include 'thkns.cmn'
      include 'atmos.cmn'
      include 'input.cmn'
      include 'in.cmn'
      include 'contrl.cmn'
      include 'out.cmn'
      include 'buff1.cmn'
      include 'buff2.cmn'
      include 'buff3.cmn'
      include 'buff4.cmn'
      include 'depolt.cmn'
c
      include 'cgas2.cmn'
      include 'czfunc.cmn'
      include 'czsing.cmn'
      include 'crefdir.cmn'
c
      include 'cfilenames.cmn'
      include 'contrb.cmn'

c======================================================================
c INITIALIZATION SECTION
c
c get the environment variables
 
      call rdenv

c read command line options

      call rdcmdln
c
      jjprint = 0
 
c initialize constants, etc.

      call iniclz
      call readin
 
c output profile file

      open (unit=23,file=outprffn,status='unknown')
 
c check to see need to open debuging file outerr

      jjprint = jprint(1) + jprint(2) + jprint(3) + jprint(4)
     1 + jprint(5) + jprint(6) + jprint(7) + jprint(8)
      if(jjprint.ne.0)then
        open(unit=33,file=outerrfn,status='unknown')
      endif
 
c if requested, open file for summary of the results

      if(jprint(9).eq.1)then
        open (unit=9,file=sumryfn,status='unknown')
      endif

c if requested, open file for iterations

      if (write_iter_file) then
         call get_lun(iter_lun)
         open (iter_lun, file=iterfn, status='unknown',
     *         form='unformatted')
      end if 

c if we're going to output the fluxes, open the file

      if(lprtflx)then
        call get_lun(iuntasc)
        open(unit=iuntasc,file=outflxfnasc,status='unknown')
        write(iuntasc,8000)
        write(6,*)'Writing flux info to ascii file:',outflxfnasc
      endif
c open file for contribution function
      if(lcontrb)then
        call get_lun(iuntcontrb)
        open(unit=iuntcontrb,file=outcontrbfn,status='unknown')
        write(6,*)'Writing contribution function info to ascii',
     &' file:',outcontrbfn
      endif
       

c write the header for the profile file

      call prthead
c
      write (23,6400) prfnam(1:8)
      write (23,6450)
      write (23,6475) pres
      write (23,6490) nthet
      write (23,6492) nscan
      write (23,6493) iazmth
      write (23,6494) nalb
      write (23,6496) wavel_start, wavel_stop
      write (23,6497)
      write (23,6498) (jprint(j),j=1,10)
      write (23,6500) (j,j=1,11)
      write (23,6600) xprf
      write (23,6700) tmpprf
      write (23,6725) (j,j=1,nthet)
      write (23,6750) (theta(j),j=1,nthet)
      write (23,6751) (scan(j),j=1,nscan)
      write (23,6752) (azmth(j),j=1,iazmth)
      write (23,6753) (alb(j),j=1,nalb)

c======================================================================
c MAIN LOOP
c

      tmp0   = 273.
      pnot   = pres
      iter_pos = 1
      num = 0
c
c move layering routine to outside the wavelength loop
            call relayr(pnot,xpnot)
c
c loop over all wavelengths

      do while (get_coef(wavel_start, wavel_stop, wavelen, alpha0,
     *        tmpcf(1), tmpcf(2), beta, rhon) .and. num .lt. max_wave)
c            if (alpha0 .eq. 0.0) alpha0 = 1.0d-10
cc for testing
c            tmpcf(1)=0.0d0
c            tmpcf(2)=0.0d0
            num = num + 1
            wavel(num) = wavelen
            if (jjprint.ne.0) write (33,7902) wavel(num)
            lambda = wavelen
            if(mod(num,100).eq.1) write (6,6200) wavelen

c get the coes for the other gases (if any)

            if(ngas.gt.1)then
              call getgascoe(wavelen)
            endif

c set the max # of iterations

            do while (iter_pos .lt. num_iter .and.
     *                wave_iter(iter_pos+1) .lt. wavelen)
              iter_pos = iter_pos + 1
            end do

            itmax  = iter(iter_pos)
c
c compute the index of refraction profile for this wavelength
c
            call crfndx(wavelen)

c           new statements

            if (jjprint .ne. 0) write(33,2121) wavelen
            call iniclz1
            call emmat(imu,imuz,iazmth,thnot,thta,azmth)

c           end of new statements

            call opthik(p,t,tl,cofx,tmp0,tmpcf)
            call dtaus(cofx,nc,xpnot,ncp1)
c store the total scattering and absorption optical depths
            tau_tota(num)=xpnot*alfaef
            tau_tots(num)=tautot-tau_tota(num)
c

c change the order of the subroutines, placing the call to slant after
c dtaus instead of after reflex. !def

            lmax=82
            layer=6
            call slant(h,ps,cofx,fracin,lmax)
            if(lsphout)then
              call exponesph(h,ps,cofx,lmax,nc,ncp1,imu)
            else
              do j=1,imu
                refdir(j)=0.0d0
                if(.not.ldown)refdir(j)=dexp(-tautot/emu(j))
              enddo
              call expone(nc,ncp1,imu)
            endif

c compute the integration factors for the upward flux at the top
c of the atmosphere

            call cek67(nc,ncp1)

c compute the reflectance and transmittance of the atmosphere 
c when illuminated from below by a diffuse source (sbar).

            call reflex(itmax,nek,ncp1)

c loop through solar zenith angles, theta

            do 2000 ijk=1,nthet
               if(jprint(9).eq.1)write (6,6300) theta(ijk)
               thenot=theta(ijk)
               if (ipsudo.eq.0)then 
                 call frstz2(z,thenot,fs,f1,f2,ncp1)
               else if (ipsudo.eq.1) then
                 if(irefrac.eq.1)then
                   call firstzr(h,ps,z,thenot,fracin,lmax,ncp1)
                 else
                   call firstz(h,ps,z,thenot,fracin,lmax,ncp1)
                 endif
               else if (ipsudo.eq.2) then
                 do j=1,nscan
                    do k=1,iazmth
                       call firstzmuop(h,ps,zsp(1,k,j),theta(ijk),
     &fracin,lmax,ncp1,scan(j),azmth(k))
                       fsp(k,j)=fs
                    enddo
                 enddo
                 call firstz(h,ps,z,theta(ijk),fracin,lmax,ncp1)
                 do j=1,nscan
                    do k=1,iazmth
                       do jlev=1,202
                          zsp(jlev,k,j)=zsp(jlev,k,j)-z(jlev)
                       enddo
                    enddo
                 enddo
               endif
               call itrate(z,ncp1,itmax,ijk)
               
c Storing output parameters (j indexes scan angle, k indexes azimuth)

               do j = 1, imu
                  do k=1,iazmth
                     i0_out(j,ijk,num,k)=eizero(j,k)
                     Z1_out(j,ijk,num,k)=Z1func(j,k)
                     Z2_out(j,ijk,num,k)=Z2func(j,k)
                     if(ipsudo.eq.2)then
                       t_out(j,ijk,num,k)=(fsp(k,j)+gg)*
     &                                    (tnstr(j)+refdir(j))
                     else
                       t_out(j,ijk,num,k)=(fs+gg)*(tnstr(j)+refdir(j))
	             endif
                     if (write_iter_file)write (iter_lun) wavelen,
     *                   theta(ijk), scan(j),
     *                   i0_out(j,ijk,num,k)/pi,Z1_out(j,ijk,num,k)/pi,
     *                   Z2_out(j,ijk,num,k)/pi,t_out(j,ijk,num,k)/pi,
     *                   sb, itmax+2,
     *                   (tnstrz(iter_num,j),iter_num=1,itmax+2),
     *                   (e0za(iter_num,j,k), iter_num = 1, itmax+2),
c    *                   (e1za(iter_num,j,k), iter_num = 1, itmax+2),
c    *                   (e2za(iter_num,j,k), iter_num = 1, itmax+2),
     *                   (ggz(iter_num), iter_num = 1, itmax+2),
     *                   (sbz(iter_num), iter_num = 1, itmax+2),
     *                   fs, gg, refdir(j), iazmth, nalb,
     *                   ((eittl(j, iaz, ialbedo), iaz = 1, iazmth),
     *                   ialbedo = 1, nalb),
     *                   ((eittr(j, iaz, ialbedo), iaz = 1, iazmth),
     *                   ialbedo = 1, nalb),
     *                   (eitu(j, iaz), iaz = 1, iazmth),azmth(k)
                  enddo
               end do

c Optional diagnostic output

               if (jprint(8).eq.1) call tbprnt(thenot)


c Write the new fluxes to the flux file, if requested

               if(lprtflx)then
                 f0_out(num)=1.0
                 f0atten_out(ijk,num)=z(ncp1)
                 gg_out(ijk,num)=gg
                 ggp_out(ijk,num)=ggp
                 write(iuntasc,'(f8.2,f5.1,1PE9.2,1PE11.3E3,1P4E10.3)')
     &wavelen,theta(ijk),f0_out(num),f0atten_out(ijk,num),gg,
     &sb,ggp,sbp
               endif
               if(lcontrb)then
c output the contribution function stuff
                 do j = 1, imu
                    do k=1,iazmth
                       write(iuntcontrb,*) wavelen,
     &theta(ijk), scan(j),azmth(k),totint(j,k)
                       do i=1,ncp1
                          if((i.ne.1).and.(i.ne.ncp1))then
                          write(iuntcontrb,*)pressure(i),
     &dJ_trans(i,j,k),
     &dtsp(i)/(pressure(i+1)-pressure(i-1))
                          else if(i.eq.1)then
                          write(iuntcontrb,*)pressure(i),
     &dJ_trans(i,j,k),
     &dtsp(i)/(pressure(i+1))
                          else if(i.eq.ncp1)then
                          write(iuntcontrb,*)pressure(i),
     &dJ_trans(i,j,k),
     &dtsp(i)/(pressure(i)-pressure(i-1))
                          endif
                       enddo
                    enddo
                 enddo
               endif

2000        continue

c End of the loop through the thetas.

         sb_out(num) = sb
         if(lprtflx)sbp_out(num)=sbp

      end do

c End of the loop through the wavelengths

c=========================================================================
c OUTPUT RESULTS THAT WEREN'T OUTPUT DURING THE MAIN LOOP
c
      nwavl = num
 
c write-out tomnval file in binary form
 
      if(jprint(10).eq.1)then
        open(unit=40, file=nvalfn,form='unformatted',
     *           status='unknown')
        if(lphiindep)then
          write (40) nscan, nthet, nwavl
          nprttabaz=1
          if(.not.lv7tabout)then
            write (40) (scan(j),j=1,nscan)
            write (40) (theta(j),j=1,nthet)
          endif
        else
          write (40) nscan, iazmth, nthet, nwavl
          if(.not.lv7tabout)then
            write (40) (scan(j),j=1,nscan)
            write (40) (theta(j),j=1,nthet)
            write (40) (azmth(j),j=1,iazmth)
          endif
          nprttabaz=iazmth
        endif
        write (40) (wavel(num), num = 1, nwavl)
        write (40) ((((i0_out(i, j, num,k), num=1,nwavl), j = 1, nthet),
     *                         i=1, nscan),k=1,nprttabaz)
        write (40) ((((Z1_out(i, j, num,k), num=1,nwavl), j = 1, nthet),
     *                         i=1, nscan),k=1,nprttabaz)
        write (40) ((((Z2_out(i, j, num,k), num=1,nwavl), j = 1, nthet),
     *                         i=1, nscan),k=1,nprttabaz)
        write (40) ((((t_out(i, j, num,k), num=1,nwavl), j = 1, nthet),
     *                         i=1, nscan),k=1,nprttabaz)
        write (40) (sb_out(num), num = 1, nwavl)
      endif

c write out binary flux files

      if(lprtflx)then
        call get_lun(iuntbin)
        close(iuntasc)
        write(6,*)'Writing flux info to binary file:',outflxfnbin
        open(unit=iuntbin,file=outflxfnbin,form='unformatted',
     &status='unknown')
        write(iuntbin) nthet, nwavl
        write(iuntbin) (theta(j),j=1,nthet)
        write(iuntbin) (wavel(num), num = 1, nwavl)
        write(iuntbin) ((f0atten_out(j,num), num=1,nwavl), j = 1, nthet)
        write(iuntbin) ((gg_out(j,num), num=1,nwavl), j = 1, nthet)
        write(iuntbin) (sb_out(num), num = 1, nwavl)
        write(iuntbin) ((ggp_out(j,num),num=1,nwavl), j = 1, nthet)
        write(iuntbin) (sbp_out(num), num = 1, nwavl)
        write(iuntbin) (tau_tots(num), num = 1, nwavl)
        write(iuntbin) (tau_tota(num), num = 1, nwavl)
        close(iuntbin)
      endif

c close all the rest of the files

      if (write_iter_file) close(iter_lun) 
      if(lcontrb) close(iuntcontrb)
      close (40)
      close (23)
      close (33)
      stop

 901   format(f4.1,1x,f5.1,1x,10(f7.2,1x))
2121        format('matscn wavelength', 1pe12.4)
6200  format(' matscn. now doing calc for wavelength ',f7.2)
6300  format(' matscn. now doing calc for sza ',f7.2)
6400  format(' ***** ',a8,' *****')
6450  format(//,' ***** input *****')
6475  format(//,' pres....',f4.2)
6490  format(' nthet...',i4)
6492  format(' nscan...',i4)
6493  format(' nazim...',i4)
6494  format(' nalb....',i4)
6496  format(' wavelength range...',f8.2,' to ',f8.2)
6497  format(' evalit',' evalrf',' expone',' itrate',' opthik',
     1 ' relayr','  slant',' tbprnt','  sumry','   nval')
6498  format(i5,9i7)
6500  format(/,' sbuv layr',t10,11i7)
6550  format(//,'dndw tables for layer',i3,/)
6600  format(' ozone amt',t12,11f7.2)
6700  format(' layer tmp',t12,11f7.2)
6725  format(/,' sza #',t10,10i7)
6750  format(' sza',t12,10f7.2)
6751  format(' scan angles',t12,15f7.2)
6752  format(' azm  angles',t12,8f7.2)
6753  format(' albedos    ',t12,8f7.2)
6775  format(/,' wavl #',t13,i7,11i12)
7000  format(' wavelen',t11,12f12.3)
7100  format(/,' alpha0',t12,12d12.5)
7200  format(' beta',t12,12d12.5)
7300  format(' tmpcof(1)',t12,12d12.5)
7400  format(' tmpcof(2)',t12,12d12.5)
7425  format(' niter',t12,i7,11i12)
7450  format(//,' ***** output *****')
7700  format(f7.3,t12,12d12.5)
7850  format('   sza')
7800  format(//,'n values for alb of ',f6.2)
7801  format('             scn of ',f6.2,' azm of ',f6.2,
     1//,'       wvl',t11,12f12.3)
7900  format(//,'dndw values for alb of ',f6.2)
7901  format('             scn of ',f6.2,' azm of ',f6.2,
     1//,'       wvl',t11,12f12.3)
7902  format(//,' ***** Output for wavelength ',f7.2,' *****',//)
8000  format(//,' EVERYTHING YOU NEED TO KNOW FOR DOWNWARD',
     &' FLUX AT SURFACE',//,
     &' Notation;',/,
     &'  mu0=cosine(solar zenith angle at the surface)',/,
     &'  f0=exoatmospheric solar flux through',
     &' a surface perpendicular to the direction',/,'     of',
     &' propagation.',/,
     &'  f0atten=the attenuated downward solar flux at the surface ',
     &'through a surface',/,'          perpendicular to the ',
     &'direction of propagation.  This does NOT',/,'          include',
     &' the cosine(solar zenith angle) factor.',/,
     &'  gg=downward diffuse flux at the surface for 0 surface',
     &' reflectance',/,
     &'  sb=reflectance of atmosphere when illuminated by upward',
     &' isotropic flux @ sfc.'/,
     &'  ggp=actinic flux due to downward diffuse at the surface for',
     &' 0 surface',/,9x,'reflectance',/,
     &'  sbp=reflectance of atmosphere (in terms of actinic flux) w',
     &'hen illuminated',/,9x,'by upward isotropic flux @ sfc.'//,
     &'  Total downward flux at sfc for reflectance rho =',/,
     &'     (f0atten*mu0+gg)/(1-rho*sb)',//,
     &'  Total actinic flux (4pi) at sfc for reflectance rho =',/,
     &'     f0atten + ggp + rho*(sbp+2)*(f0atten*mu0+gg)/(1-rho*sb)',//,
     &'  Total downward actinic flux at sfc for reflectance rho =',/,
     &'     f0atten + ggp + rho*sbp*(f0atten*mu0+gg)/(1-rho*sb)',//,
     &'wavelength Solar  f0',4x,'f0atten',4x,'gg',8x,'sb',8x,'ggp',7x,
     &'sbp',/,11x,'Zenith',/,11x,'[deg]',/,77('*'))
       end