        Options  /Extend_Source
      program uv_flux_overpass_bsza
c
c***********************************************************************
c
c     author
c       colin seftor
c
c     wornings:
c     only N7(ifm=1) or EP(ifm=3) data can be processed
c     only first 300000 rows in input file will be processed
c     ? still use noon sza in  call ery(
c     
c
c     date
c       23-jan-2003 use snow_bound=0.2 and br_sn>0
c       20-jan-2003 extended output: similar to version 1
c       20-jan-2003 Update snow correction algorithm       
c       19-jan-2003 TABLE.REFEP and TABLE.REFN7 touse with both satellites
c       16-jan-2003 to change to br_sza_av in call uvf(...)       
c       15-jan-2003 to read Sodankyla-toms timeseries input data
c       28-mar-2001
c
c     purpose
c       calculate Brewer-time surface uv flux for 305, 310, 325, and 380 nm
c       calculate ? noon-time erythemally-weighted flux
c       calculate erythemal dose (integral of erythemally-weighted flux
c                                 over the whole day)
c
c     method
c       performs table lookups on pre-computed flux values using
c       the TOMS Version 7 ozone and temperature profiles and
c       interpolation techniques 
c
c     variables
c       name       type     description
c       ----       ----     -----------
c       lat         r       latitude
c       long        r       longitude
c       pteran      r       surface pressure , [atm] (1atm=1013.25mb)
c       cofs       r(4,4)   lagrange coefficients
c       densol     r(4,7)   used in lagrange coef calc
c       denscn     r(4,3)   used in lagrange coef calc
c       latind     int(2)   latitude index (used in profile mixing fraction)
c       rad         r(3)    radiance
c                           1 = 360 nm
c                           2 = 340 nm for n7 ; 360 for EP
c                           3 = 380 nm for n7 ; 360 for EP
c       rad331      r       331 nm radiance
c       mler       real     minimum surface reflectivity
c       salb       real     surface albedo
c       ler        r(3)     lambert equivalent refl at terrain pressure
c                           1 = 360 nm
c                           2 = 340 nm for n7
c                           3 = 380 nm for n7
c       ler_p      r(2)     lambert equivalent refl at two pressure levels
c                           for 360 nm
c       res        real     331 nm residue
c       p1         r        angular factor to convert z1 and z2 to i1 and i2
c       pr         r        angular factor to convert z1 and z2 to i1 and i2
c       algflg     int      snow/ice flag (algflg >= 10)
c       surf       int      surface type
c       uv         r(4)     uv surface flux for 305, 310, 325, 380 nm
c       uvclear    r(4)     clear scene uv surface flux for 305, 310, 325, 380 nm
c       ramcor   r(4,3,2)   raman scattering correction 
c                           4 => 4 wavelengths 331, 360, 340, 380 nm
c                           3 => 3 pressures: 1.0, 0.7, 0.4 atm
c                           2 => 2 reflectivities: 15%, 80% ler
c       ssza       r(24)    solar zenith angles for each of the 24 hours
c       iloc       i(24)    table indices for sza at each of the 24 hours
c       wav        r(4)     wavelengths to calculate flux for
c       g10        r(5)     aerosol correction factors for 1.0 atm
c       g08        r(5)     aerosol correction factors for 1.0 atm
c       gcor       r(24)    aerosol corrections, max of 24 
c       szanode    r(5)     sza node points for g08 and g10
c
c    include files
c       const.inc           constant common block
c
c***********************************************************************
c
      include 'const.inc'
      real lat, long
      real densol(4,7), denscn(4,3), ramcor(4,3,2)
      real ramcor_n7(4,3,2),ramcor_ep(4,3,2)
      real p1, pr, mler, cofs(16), salb, ler(3), ler_p(2)
      real ssza(24), tcl4
      real prsnode(4), taunode(11), pres(3)
      real rad(3), rad331, totoz(2)
      real*8 rrad
      real*8 lsza,lsatza,lphi,lpteran,lrad,tcl,lmref,taucloud,cldtau
      real uv(4), uvclear(4)
      real eflux, eclear
      real edflux, edclear 
      real g08(5), g10(5), gcor(24), szanode(5)
      integer surf,algflg, latind(2), errflg, indprs, indtau
      integer iloc(24), indprf(2)
      integer seq
      real n331, n340, n360, n380

c     add variables for timeseries
      	integer ifm,nsite,isp, n_mbr
      	character*2100 header
      	logical success
	real  year, day, time,sat_toms
        real lat_toms, long_toms,pteran, pcloud, sza_toms,phi_toms, o3bst 
        real ref, ai_toms, soi, nwl3,nwl4,nwl5,nwl6
        real sn_br,sun_glint_ang,br305,brm305,br310,brm310,br324,brm324
        real br_cie,br_sza_min,br_sza_max,br_sza_av,dist
c
      data wav/305.088716,310.089969,324.093502,380.107892/
      data taunode/0.,1.,3.,5.,10.,20.,30.,40.,50.,70.,100./
      data prsnode/0.,-0.2231,-0.5108,-0.9163/
      data xzlog/0.0000000, 0.1438410, 0.3465740, 0.6931472,
     1           1.0728856, 1.4918724, 1.8551181, 2.2582959,
     2           2.6627451, 3.3552829/
      data xlog /0.0000000, 0.034668232, 0.14384104, 0.34657359,
     1           0.69314718, 1.0728856/
      data rrad /57.29578/
      data pres /1.0,0.7,0.4/
      data g10/0.055,0.070,0.090,0.120,0.20/
      data g08/0.055,0.070,0.090,0.120,0.20/
      data szanode/0,20.,40.,60.,80./
c     -- ** SATELLITE SPECIFIC INFORMAION ** --
c     -- N7
c	---          331    360    340   380nm
      data ramcor_n7/0.196,-0.853,-0.141,0.382,
     1               0.159,-0.677,-0.113,0.298,
     2               0.122,-0.501,-0.085,0.214,
     3               0.248,-0.958,-0.173,0.396,
     4               0.178,-0.670,-0.123,0.273,
     5               0.100,-0.371,-0.069,0.151/
c     -- EP
      data ramcor_ep/0.166,-0.704,0.,0.,
     1               0.135,-0.559,0.,0.,
     2               0.104,-0.414,0.,0.,
     3               0.210,-0.790,0.,0.,
     4               0.151,-0.553,0.,0.,
     5               0.084,-0.306,0.,0./
c
c     -- angle conversions
      rtod=180./3.14159
      dtor=3.14159/180.

c     -- calculate values needed in ler reflectivity table interpolation
      call dencal(densol,denscn)



c     -- open output file 
      open (25,file='sodankyla_ovp_v15_brsza.bound02',status='unknown')
      write (25,998)             ! write output header
      write (25,*) 'site yr  day  time isp    lat  long  ozone mler salb ler(360) ref(V7) 305_clear 310_clear 
     &324_clear, ery_clear, tcl,  305  310    324   ery algflg, ai(V7) res331 br_sza, br305,br310,br324,br_cie,
     & dist, snow'

c     -- open input file 
      open (24,
     *file='/misc/krotkov/OMI_UV/VALIDATION/FIN/sodankyla_all.input',
     *               ERR=900,status='old')
c     -- read input header
      do irec = 1,39      
      read (24,*)
      enddo

c	-- loop through the input data

	sat_year=0

      do irec = 1,300000

        read(24,*,err=901,end=800) nsite, year, day, time, isp, sat_toms, 
     1  lat_toms, long_toms,pteran, pcloud, sza_toms,phi_toms, o3bst, 
     2  ref_toms, ai_toms, soi, nwl3,nwl4,nwl5,nwl6,algflg,errflg,surf,
     3  sn_br,sun_glint_ang,br305,brm305,br310,brm310,br324,brm324,
     4  br_cie,br_sza_min,br_sza_max,br_sza_av,n_mbr,dist

c       write(6,*                ) nsite, year, day, time, isp, sat_toms, 
c    1  lat_toms, long_toms,pteran, pcloud, sza_toms,phi_toms, o3bst, 
c    2  ref_toms, ai_toms, soi, nwl3,nwl4,nwl5,nwl6,algflg,errflg,surf,
c    3  sn_br,sun_glint_ang,br305,brm305,br310,brm310,br324,brm324,
c    4  br_cie,br_sza_min,br_sza_max,br_sza_av,n_mbr,dist

c       read (24,*,end=800) iyr,iday,lat,long,pteran,algflg,sza,satza,
c    1   phi,ozone,n331,n340,n360,n380					! global input format

	iyr 	=int(year)
	iday	=int(day)
	lat	=lat_toms
	long	=long_toms
	pteran	=pteran
	algflg	=algflg
	sza	=sza_toms
	satza	=sat_toms
	phi	=phi_toms
	ozone	=o3bst

c       -- assign correct N values                    
	        IF(year.LE.1993) THEN
        ifm=1
	n331	=nwl3
	n340	=nwl4
	n360	=nwl5
	n380	=nwl6
c       write(6,*) 'FMI=',ifm,' : N7/TOMS'
        	ELSE
        ifm=3
	n331	=nwl5
	n340	=nwl6	! not used for EP  
	n360	=nwl6
	n380	=nwl6	! not used for EP  
c       write(6,*) 'FMI=',ifm,' : EP/TOMS'
        	ENDIF		!  end channel selection   

c       -- convert 360 and 331 nm nvalue to radiance
        rad(1) = 10.**(-.01*n360)
        rad(2) = 10.**(-.01*n340)	! equal 331 for EP
        rad(3) = 10.**(-.01*n380)	! equal 360 for EP
        rad331 = 10.**(-.01*n331)

c-- statements before 101 are needed only then starting a new year
c-- if the year is the same just skip everything until label 101

        IF(year.EQ.sat_year) GO TO 101
        sat_year=year
        write(6,*) 'read TABLE.REF for: ifm= ',ifm,'	sat_year= ',sat_year

c--      read appropriate TABLE.REF for EP or N7:
c--      for year >=1996 (EP/TOMS): ifm=3
c--      for year <=1993 (N7/TOMS): ifm=1

      call rdtables(ifm,success)
c       if(.NOT.success) goto 899

c--      choose appropriate Raman correction
c--      for year >=1996 (EP/TOMS): ifm=3
c--      for year <=1993 (N7/TOMS): ifm=1

	        IF(year.LE.1993) THEN
        DO i=1,4
                DO j=1,3
                	DO k=1,2
        ramcor(i,j,k)=ramcor_n7(i,j,k)
c       ramcor   r(4,3,2)   raman scattering correction
c                           4 => 4 wavelengths 331, 360, 340, 380 nm
c                           3 => 3 pressures: 1.0, 0.7, 0.4 atm
c                           2 => 2 reflectivities: 15%, 80% ler
                	ENDDO
                ENDDO
        ENDDO
        	ELSE
        DO i=1,4
                DO j=1,3
                	DO k=1,2
        ramcor(i,j,k)=ramcor_ep(i,j,k)
c       ramcor   r(4,3,2)   raman scattering correction
c                           4 => 4 wavelengths 331, 360, 340, 380 nm
c                           3 => 3 pressures: 1.0, 0.7, 0.4 atm
c                           2 => 2 reflectivities: 15%, 80% ler
                	ENDDO
                ENDDO
        ENDDO
        	ENDIF		!  end Raman correction    

 101    CONTINUE ! skip tables initialization if the year and satellite is the same, just different day

c       -- determine pressure index for both flux and nr tables
c       -- in log of pressure for flux table
        ptlog  = alog(pteran)
        call locate(prsnode,4,ptlog,indprs)
        fr_pr  = (ptlog-prsnode(indprs)) /
     1           (prsnode(indprs+1)-prsnode(indprs))
c       -- in pressure for nr table
        call locate(pres,3,pteran,ipres)
        pfrac = (pteran-pres(ipres)) / (pres(ipres+1)-pres(ipres))
c       -- day of year
        doy = float(iday)
c       -- sun-earth distance
        phase = 0.0172024238D0 * (dfloat(doy)-3.4)
        r  = 1.0 - 0.0167*dcos(phase)
        r2 = r*r
c       -- get minimum uv reflectivity (mler) from database
        call gmref(doy,lat,long,mler)  
c       -- get table indices, lagrange coefficients, and angular factors
        call tabpos(sza,satza,densol,denscn,indsol,indscn,iofset,cofs,
     1   p1,pr)
c       -- determine profile mixing fraction
        call pmfrac(lat,latind,prfrac)
c       -- get ozone profile indices
        call prfind(latind,pteran,ozone,indprf,totoz)
c       -- lambertian equivalent reflectivity calculation (ler)
c       -- ler(1) = 360 nm with Raman correction 
c       -- for nimbus we also have
c       -- ler(2) = 340 nm	! ler(2) = 360 nm for EP (without Raman)
c       -- ler(3) = 380 nm      ! ler(3) = 360 nm for EP (without Raman)
c       -- ************ set for ler(3)=380 nm for n7  *****************
c       nref = 3
c       -- ****** set for ler(1) = 360 nm for both n7  and ep   *******
        nref = 1    

        call reflec(nref,ipres,pfrac,iofset,cofs,phi,p1,pr,pteran,
     1   rad,ramcor,ler_p,ler,indram)
c       -- calculate 331 nm residue
        call residue(ipres,pfrac,indprf,iofset,cofs,phi,p1,pr,pteran,
     1   rad331,ramcor,indram,totoz,prfrac,ler_p,res)


c       -- set surface albedo to mler
        salb = mler

c       -- modify salb if there is snow or ice (steps 1,2,3):
c	-- [Krotkov et al, Opt. Eng.,41,12,p.3033,2002]

c---    Step I)  First, determine if there is snow on the ground:
c       if (algflg.ge.10) then  ! use TOMS climatology
        if (sn_br.gt.0) then	! Use measured snow depth at the site

c       Step II) Determine a snow albedo threshold (snowrefl) for a given spot on a given day:
c---    Currently the snowrefl is the mler value bounded from below by a constant value (snow_bound)
c---    This allows us to account for higher snow albedo
c---    for regions with permanent snow cover: Greenland, Antarctica, Arctic )
c---    The value 0.4 was selected as appropriate for snow covered urban/suburban populated ares i
c---    containing at least moderate densities of roads, houses, and trees.
c---    ! snow_bound can be adjusted for any specific site

c	snow_bound=0.4    ! default TOMS global value [Krotkov et al 2002]
	snow_bound=0.2    ! bounded snowrefl

          if (mler.lt.snow_bound) then
            snowrefl = snow_bound
          else       
            snowrefl = mler  
          endif	  

c       Step III) After snowrefl has been established, we estimate the actual snow albedo, by
c       	  comparing snowrefl with the measured LER at 360nm: ler(1)
c  algflg   ler(1)       	Use this reflectivity for salb   
c     ----   ----        	---------------------
c    =>10    >snowrefl+0.05   	snowrefl     (assumes this means clouds over snow)
c    =>10   <=snowrefl   	ler(1)  (so tau_cl will be zero; clear sky over snow)

          if (ler(1).lt.snowrefl+0.05) then
            salb = ler(1)       !!! make sure in this case tau_cl is also forced to zero!
          else   
            salb = snowrefl	! cloud over snow
          endif

        endif     ! ---  END modify salb if there is snow or ice


c-----        !!! ler(1) is hard coded in all calls to subroutines below

c       -- if 360 nm ler le mler, we have a clear scene. set tcl = 0
c       -- if mler >= 0.9 or terrain pressure < 0.75, assume we 
c       -- have a clear scene and set tcl = 0
c       -- otherwise calculate tcl 
c       if (ler(1).le.mler .or. salb.ge.0.9 .or. pteran.lt.0.75) then
        if (ler(1).le.salb .or. salb.ge.0.9 .or. pteran.lt.0.75) then
          tcl = 0.d0
        else
c         -- correct for cloud (using Celarier's tables and code)
          lsza    = sza
          lphi    = phi
          lsatza  = satza
          lpteran = pteran
          lmref   = salb
          lrad    = rad(1)	! USE 360nm nr
c         -- get cloud optical depth
          tcl=taucloud(lsza,lphi,lsatza,lpteran,lmref,lrad)
c         if (algflg.ge.10 .and. ler(1).gt.0.4) then
          if (sn_br.gt.0   .and. ler(1).gt.0.4) then
c           -- use linear interpolation for snow/ice cases
c	CHECK THAT RADHIGH.DAT and RADLOW.DAT correspond to 360nm -- nk 0119.03
            tcl = cldtau(lsza,lphi,lsatza,lpteran,lmref,lrad)
          else
c           -- use fancy lagrange stuff otherwise
            tcl = taucloud(lsza,lphi,lsatza,lpteran,lmref,lrad)
          endif
          if (tcl.lt.0.) tcl = 0.
        endif
c       -- determine cloud optical depth index
        tcl4 = tcl
        call locate(taunode,11,tcl4,indtau)
        fr_tau = (tcl-taunode(indtau)) /
     1           (taunode(indtau+1)-taunode(indtau))
c       -- get sza for all hours of the day, keep ones with sza <= 88.
c       -- also determine local noon (ie minimum) sza
        icnt = 1
        szamin = 90.
c       -- get sun solar zenith angles every hour
c       -- save those less than 88 degrees
c       -- get sun declination angle
        call sunang1(rrad,iday,sdec,tc)
        do i = 1,24
          hr   = float(i) - long/15. 
          if (hr.lt.0.) then
            hr = hr + 24.
          else if (hr.gt.0.) then
            hr = hr - 24.
          endif
c         -- get sun solar zenith angle for hour hr
          call sunang2(rrad,hr,sdec,tc,long,lat,szatmp,saz)
          if (szatmp.gt.0. .and. szatmp.le.88.) then
            ssza(icnt) = szatmp
            if (szatmp.lt.szamin) then
              szamin = szatmp
              imn    = icnt
            endif
            call locate(xzlog,10,alog(1./cos(dtor*szatmp)),iloc(icnt))
c           -- determine aerosol correction factor
            call locate(szanode,5,szatmp,indsza)
            if (indsza.lt.5) then
              szafrac = (sza - szanode(indsza)) /
     1                  (szanode(indsza+1) - szanode(indsza))
              gc10 = (1.-szafrac)*g10(indsza) + szafrac*g10(indsza+1)
              gc08 = (1.-szafrac)*g08(indsza) + szafrac*g08(indsza+1)
              gcor(icnt) = -1. * ((1.-pfrac)*gc10 + pfrac*gc08)
            else
              szafrac = (sza - szanode(4) / (szanode(5)) - szanode(4))
              gc10 = (1.-szafrac)*g10(4) + szafrac*g10(5)
              gc08 = (1.-szafrac)*g08(4) + szafrac*g08(5)
              gcor(icnt) = -1. * ((1.-pfrac)*gc10 + pfrac*gc08)
            endif
            iloc(icnt) = iloc(icnt) - 1
            if (iloc(icnt).lt.1) iloc(icnt) = 1
            if (iloc(icnt).gt.7) iloc(icnt) = 7
            icnt = icnt+1
          endif
        enddo
        icnt = icnt-1

c       -- flux calculation
        if (ssza(imn).le.88.) then
c         -- local noon 305, 310, 325, and 380 nm uv fluxes
c         call uvf(iloc(imn),indprs,indprf,tcl,ssza(imn),pteran,prfrac,
c         -- Use average Brewer SZA for uv fluxes
          call uvf(iloc(imn),indprs,indprf,tcl,br_sza_av,pteran,prfrac,
     1     fr_pr,totoz,ler(1),salb,res,gcor(imn),uvclear,uv)

c         -- local noon erythemal flux and daily erythemal dose
c	  CHANGE to Brewer sza?
          call ery(imn,icnt,iloc,indprs,indtau,indprf,ssza,pteran,
     1     prfrac,fr_pr,fr_tau,totoz,ler(1),salb,res,gcor,eclear,
     2     edclear,eflux,edflux)     

c         write (25,'(3f8.2,4f8.2,2f9.2,4f8.2,2f9.2)')	!  global OMI output
c    1     lat,long,100.*ler(1),uvclear,eclear,edclear,uv,eflux,edflux

        write(6,*) nsite, year, day ,lat_toms, long_toms            

          write (25,8)
     &                 nsite,           ! site number
     &                 int(year),             !4-digit year number yyyy
     &                 int(day),              !day number  doy
     &                 int(time),             !gmt seconds from midnight
     &                 isp          ,    !scan position
c    &                 sat_toms     ,    !satellite zenith angle [deg]
     &                 lat_toms,              !fov latitude
     &                 long_toms,             !fov longitude
c    &       sol_dist_corr      ,             !sol dist corr
     &                 o3bst,            !total ozone
     &                 100.*mler,        ! Mler                  
     &                 100.*salb   ,     ! TOMS assumed surface albedo  (with snow correction)
     &                 100.*ler(1)    ,  ! LER @ 360 at Pteran: ler(1)
     &                 1.  *ref_toms,    !	Official V7 reflectivity ( from L2)
     &           (uvclear(j),j=1,3)  ,   ! TOMS clear sky irradiance @ 305, 310, 324
     &                 eclear  ,   	 ! TOMS noon clear sky erythemal flux 
     &                 tcl,              ! TOMS Cloud optical thickness at 380nm
     &           (uv(j),j=1,3)  ,        ! TOMS irradiance @ 305, 310, 324
     &                 eflux   ,   	 ! TOMS noon erythemal flux 
     &                 algflg,               ! TOMS alg.flag (if >10 Snow probability >=50%)
c    &                 ct_mean,            ! mean CT for weights 1:5
c    &                 0.01*ref   ,               ! V7 R380 (<=1993)  or R360 (=>1996)
c    &                      ref_v6_pc ,           ! V6 at Pcloud R380 (<=1993)  or R360 (=>1996)
c    &                 t325 ,                ! F325 from UVRSR
c    &                 idist,            !distance from fov center to site
c    &                 pteran,           !terrain pressure at fov center
     &                 ai_toms   ,          !  Official TOMS AI (cloud fraction algorithm, 340/380 and 331/360)
     &                 res       ,          !  331 residue with 360 ler
c    &                 isoi,             !sulfur oxide index
     &                 br_sza_av,              !solar zenith angle
     &  1000.*br305,
     &  1000.*br310,
     &  1000.*br324,
     &  1000.*br_cie,
     &                 dist,                 ! dist, km
     &                 sn_br                 ! Brewer snow depth, cm




        endif	!   write only for SZA < 88 deg


      enddo

 800  continue
      close (24)
      close (25)

      stop
998   format(' site,    0       ! WOUV site number [nn]              '/
     &       ' year,    1       ! TOMS: 4-digit year number [yyyy]'/
     &  'day,           2       ! TOMS: day number,(day of year) '/
     &  'time,          3       ! TOMS: overpass GMT [seconds] from midnight'/
c    &  'satZA          4       ! TOMS: satellite zenith angle [deg] '/
     &  'isp            4       ! TOMS: scan position (18-nadir; right to left - N7; left to right -EP)'/
     &  'lat,           5       ! TOMS latitude of the center of FOV'/
     &  'lon,           6       ! TOMS longitude of the center of FOV'/
c    &  '1/r2           8       ! TOMS Solar distance correction     '/
     &  'o3bst,         7       ! TOMS: total colunm ozone'/
     &  'mler           8       ! MLER (%)                                        '/
     &  'salb    ,      9       ! TOMS: surface albedo (with snow correction) '/
     &  'ler(1)  ,      10      ! TOMS LER V6 at Pteran (1996)'/
     &  'ref (V7),      11      !	Official V7 reflectivity ( from L2)'/
     &  't305_clear     10      ! TOMS: Fclear at 305nm FWHM=0.55nm and br_av_sza'/
     &  't310_clear     11      ! TOMS: Fclear at 310nm FWHM=0.55nm and br_av_sza'/
     &  't324_clear     12      ! TOMS: Fclear at 324nm FWHM=0.55nm and br_av_sza'/
c    &  't380_clear     15      ! TOMS: Fclear at 380nm FWHM=0.55nm and br_av_sza'/
     &  'tcie_clear     13      ! TOMS: Fclear , erythemally weighted (CIE weighting)'/
     &  'tcl            14      ! TOMS: Effective cloud optical thickness            '/
     &  't305_cloud     15      ! TOMS: Fcloud at 305nm'/
     &  't310_cloud     16      ! TOMS: Fcloud at 310nm'/
     &  't324_cloud     17      ! TOMS: Fcloud at 324nm '/
c    &  't380_cloud     21      ! TOMS: Fcloud at 380nm '/
     &  'tcie_cloud     18      ! TOMS: Fcloud, CIE weighting'/
     &  'algflg,        21      ! TOMS alg.flag (if >10 Snow probability >=50%)'/
     &  'ai_toms,       19      ! TOMS AI (V7)                                 '/
     &  'res   ,        20      ! TOMS residue @ 331nm (ler @ 360nm)           '/
     &  'br_sza_av      20      ! Brewer: SZA_AV '/
     &  'br305          21      ! Brewer: Fcloud at 305nm'/
     &  'br310          22      ! Brewer: Fcloud at 310nm'/
     &  'br324          23      ! Brewer: Fcloud at 324nm'/
     &  'br_cie         24      ! Brewer: Fcloud, CIE weighting'/
     &  'dist           25      ! distance between station and the center of the TOMS FOV, km'/
     &  'snow           26      ! Brewer snow depth; -1 -no snow; 0 - data     '/
     &    )
    8 format(
     &       1x, I3,         ! site #
     &       1x, I4.4,       !year
     &       1x, i3.3,       !day
     &       1x, I5,         !gmt seconds
     &       2x, I2,         !scan position
c    &       1x, F5.2,       !satellite zenith angle
     &       1x, F6.2,       !lat_toms
     &       1x, F7.2,       !lon
c    &       1x, F6.2,       !sol_dist_corr
     &       1x, F5.1,       !total ozone
     &       1x, F6.2,      ! MLER, %                                                  
     &       1x, F6.2,      ! TOMS assumed surface albedo (with snow correction)
     &       1x, F6.2    ,  ! LER @ 360 at Pteran: ler(1)
     &       1x, F6.2,      !	Official V7 reflectivity ( from L2)
     &       4(1x, F6.2)  , ! Overpass Clear sky irradiance for weights 1:5
     &       1x, F6.2  ,    ! Overpass tcl at 360nm
     &       4(1x, F6.2)  ,    ! Overpass All sky  irradiance for weights 1:5
     &       1x, I4,        ! TOMS alg.flag (if >10 Snow probability >=50%)
c    &       2x, F4.2,      ! CT =F_cloud/F_clear mean for weights 1:5
c    &       1x, F5.3,      ! V7 R380 or R360
c    &       1x, F6.3,      ! f325 from UVRSR (total)
c    &       1x, F6.3,      ! V6 pcloud R380 or R360
c    &       1x, F5.1,      ! distance, km
c    &       1x, F5.2,      !terrain pressure at fov
     &       1x, F5.2,       !  TOMS V7 aerolsol index
     &       1x, F5.2,       !  residue 331 (360 ler)
c    &       1x, I4,         !sulfur oxide index
     &       1x, F5.2,       ! Brewer solar zenith angle
     &       4(1x, F6.2)  ,  ! Brewer fluxes
     &       1x, f5.1,         !distance
     &       1x, f5.1)      ! Brewer snow depth, cm

c-- error handling
 
 900    continue
        WRITE(6,*) 'main: could not open sodankyla_all.input file '
        STOP
 901    continue
        WRITE(6,*) 'main: error while reading data in sodankyla_all.input file '
        STOP
 
 899    continue
        write(6,*) 'rdtables: could not read TABLE.REF tables'
        stop
 
      end
