pro tomrad,pressure,sza,satza,phi,ref,wav,c0,c1,c2,bet, $
 rhon,prfoz,prftemp,iter,pol,i0,i1,i2,t,sb,itot,nval
;
; procedure matmlt
;
; author:
;
;    colin seftor
;
; date:
;
;    22-february-96
;
; purpose:
;
;    provide a front-end interface to the mateer radiative transfer
;    algorithm
;
; category:
;
;    radiative transfer routine
;
; calling sequence:
;
;    pro tomrad,pressure,sza,satza,phi,ref,wav,c0,c1,c2,bet, $
;     rhon,prfoz,prftemp,iter,pol,i0,i1,i2,t,sb,nval
;
; input:
;
;    pressure:      pressure
;
;    sza:           solar zenith angle(s); in degrees; max of 10
;
;    satza:         satellite zenith angle(s); in degrees; max of 6
;
;    phi:           phi angle(s); in degrees; in degrees; no maximum
;
;    ref:           reflectivity(s); in percent; max of 15
;
;    wav:           wavelength(s); in angstroms; max of 2500
;
;    c0:            absorption coefficient(s); max of 2500
;
;    c1:            1st temperature coefficient(s); max of 2500
;
;    c2:            2nd temperature coefficient(s); max of 2500
;
;    beta:          rayleigh scattering coefficient(s); max of 2500
;
;    rhon:          molecular depolarization factor(s); max of 2500
;
;    prfoz:         ozone profile: real array of 11
;
;    prftemp:       temperature profile: real array of 11
;
;    iter:          number of iterations for each wavelength; number
;                   is equal to numwav; max num of iterations is 12 
;
;    pol:           turns on (1) or off (0) molecular anisotropy
;
; output:
;
;    i0:            i0
;
;    i1:            equals z1 * q1
;
;    i2:            equals z2 * q2
;
;    t:             t
;
;    sb:            sb
;
;    itot:          radiance 
;
;    nval:          the nvalue
;
;    the above are used in the following formula:
;
;                                              t*ref
;    itot = i0 + i1*cos(phi) + i2*cos(2phi) + --------
;                                             1-sb*ref
;
;    the n value is then
;
;    n = -100. * log(itot/pi)
;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
input_file  = 'PROF'
jprint      = [0,0,0,0,0,0,0,0,0,1,1] ; jprint(11) MUST be left on
numsza      = n_elements(sza)
numsatza    = n_elements(satza)
numwav      = n_elements(wav)
numref      = n_elements(ref)
numphi      = n_elements(phi)
sb          = dblarr(numwav)
q1          = dblarr(numsza,numsatza)
q2          = dblarr(numsza,numsatza)
i0          = dblarr(numwav,numsza,numsatza)
i1          = dblarr(numwav,numsza,numsatza)
i2          = dblarr(numwav,numsza,numsatza)
z1          = dblarr(numwav,numsza,numsatza)
z2          = dblarr(numwav,numsza,numsatza)
t           = dblarr(numwav,numsza,numsatza)
ground_term = dblarr(numwav,numsza,numsatza,numref)
itot        = dblarr(numwav,numsza,numsatza,numref,numphi)
nval        = dblarr(numwav,numsza,numsatza,numref,numphi)
dum         = lonarr(3)
dumwav      = dblarr(numwav)
; 
; -- open file to store absorption, temperature, and rayleigh
; -- scattering coefficients and molecular depolarization factor
;  openw,20,'ABS_COEF.DAT'
  openw,20,'coe.dat'
  printf,20,'       WLVAC         C0          C1          C2             ' $
  +'BETA         RHO'
; -- store them
  for i = 0,numwav-1 do begin
     printf,20,10.*wav(i),c0(i),c1(i),c2(i),bet(i),rhon(i)
  endfor
;
; -- build a new environmental control file, ENV
;
;spawn,'rm ENV'
;openw,19,'ENV'
;printf,19,'lv7tab = T'
;printf,19,'lv7tabout = T'
;printf,19,'prf_type = 1'
;printf,19,'gc_type=1'
;printf,19,'lprtflx = F'
;
;close,19

  close,20
; -- open file to store input parameters to mateer code
  openw,21,input_file
; -- store parameters
  two=2.0
  printf,21,'INPUTPRF'
  printf,21,pressure,format='(f4.2)'
  printf,21,numsza,format='(i2)'
  printf,21,sza,format='(10f6.2)'
  printf,21,numsatza,format='(i2)'
  printf,21,satza,format='(10f6.2)'
  printf,21,numphi,format='(i2)'
  printf,21,phi,format='(10f6.1)'
  printf,21,numref,format='(i2)'
  printf,21,ref/100.,format='(10f6.3)'
  printf,21,10.*wav(0),10.*wav(numwav-1),format='(2f8.2)'
  printf,21,prfoz,format='(11f7.2)'
  printf,21,prftemp,format='(11f6.1)'
  printf,21,jprint,format='(11i2)'
  printf,21,two,format='(i2)'
  printf,21,2400.0,2900.0,format='(12f9.2)'
  printf,21,iter,format='(12i3)'
  printf,21,pol,format='(i3)'
  close,21
; -- call tomrad code to do calculations
   spawn,'src/tomrad > /dev/null'
; -- open output file that tomrad code produced
  openr,22,'tomnval.dat',/f77_unformatted
; -- read in parameters 
  readu,22,dum
  readu,22,dumwav
  readu,22,i0
  readu,22,z1 
  readu,22,z2
  readu,22,t
  readu,22,sb
  close,22
;
; -- calculate q1 and q2 
  q1 = (-3./8.) * ((cos(sza*!dtor) * sin(sza*!dtor)) # sin(satza*!dtor))
  q2 = (3./32.) * (sin(sza*!dtor))^2 $
   # ((sin(satza*!dtor))^2 / cos(satza*!dtor))
; -- calculate i1, i2, and ground term 
  for iwav = 0,numwav-1 do begin
     i1(iwav,*,*) = z1(iwav,*,*) * q1(*,*)
     i2(iwav,*,*) = z2(iwav,*,*) * q2(*,*)
     for iref = 0,numref-1 do begin
;       -- turn reflectivity into a fraction
        frac_ref = ref(iref)/100.
;       -- calculate ground term
        ground_term(iwav,*,*,iref) = t(iwav,*,*)*frac_ref        $
                                   / (1.-frac_ref*sb(iwav))
     endfor
  endfor
; -- calculate itotal
  for iphi = 0,numphi-1 do begin
     for iref = 0,numref-1 do begin
        itot(*,*,*,iref,iphi) = i0(*,*,*)                        $ 
                              + i1(*,*,*)*cos(phi(iphi)*!dtor)   $
                              + i2(*,*,*)*cos(2*phi(iphi)*!dtor) $
                              + ground_term(*,*,*,iref)
     endfor
  endfor
; -- calculate n value
  nval = -100.*alog10(itot/!pi)
;
return
end