disort _ function(input.list = disort.input, 
		  layer.tau = input.list$layer.tau,   
		  ssa = input.list$ssa, 
		  moments = input.list$moments, 
		  mu0 = input.list$mu0,   
		  observation.tau = input.list$observation.tau, 
		  observation.mu = input.list$observation.mu, 
		  observation.phi = input.list$observation.phi, 
		  surface.albedo = 0) {

# A vertically homogenous medium has 
# length(layer.tau) = length(ssa) = ncol(moments) = 1
# Otherwise the 3 lengths must agree. 
  nlayers _ length(layer.tau) 
  if (nlayers == 1) {
    moments _ moments/moments[1]
  } else {
    ssa _ rep(ssa, length = nlayers)
    if(!is.matrix(moments)) 
       moments _ matrix(rep(moments, nlayers), ncol = nlayers)
    moments _ moments/rep(moments[1, ], rep(nrow(moments), ncol(moments)))
  }

# Disort is quite picky about its input. Here we test for situations which 
#  will make the subroutine fail. 
  if(any(ssa  < 0) | any(ssa  > 1)) 
    stop("Single scattering albedo out of bounds.")
  if(any(moments < 0) | any(moments > 1)) 
    stop("Legendre moments out of bounds.")
  if(any(observation.mu < -1) | any(observation.mu > 1))
    stop( "Mus out of bounds.")
  if(length(observation.mu) > 1 & any(diff(observation.mu) <= 0)) 
    stop("Observation polar angles must be ascending.")
  if(any(observation.tau > sum(layer.tau)) | any(observation.tau < 0)) 
    stop("Observation level requested outside bounds of medium")
  if(length(observation.tau) > 1 & any(diff(observation.tau) <= 0)) 
    stop("Observation levels observation.mu must be in increasing order.")
  if(any(observation.phi < 0) | any(observation.phi > 360)) 
    stop("Observation azimuthal angles out of bounds")
 
# The Fortran code should be loaded when the library is attached. 
  if(!is.loaded(symbol.For("disort_driver"))) 
    stop("Fortran object code hasn't been loaded.")
  results _ .Fortran("disort_driver", 
	     nlayers   = as.integer(nlayers),
	     layer.tau = as.single (layer.tau), 
	     ssa       = as.single (ssa), 
	     moments   = as.single (array(moments,   
				 dim = c(length(moments)/nlayers, nlayers))),
             nmoments  = as.integer(length(moments)/nlayers - 1), 
	     observation.tau      = as.single (observation.tau),    
             nobservation.tau     = as.integer(length(observation.tau)),
	     mu0       = as.single (mu0 ),
	     nmu0      = as.integer(length(mu0 )),
	     observation.mu       = as.single (observation.mu ),    
             nobservation.mu      = as.integer(length(observation.mu )),
	     observation.phi      = as.single (observation.phi),    
             nphi      = as.integer(length(observation.phi)),
	     surface.albedo       = as.single (surface.albedo), 
	     direct    = single    (length(observation.tau)), 
	     diff.dn   = single    (length(observation.tau)), 
	     diff.up   = single    (length(observation.tau)), 
	     dfdt      = single    (length(observation.tau)), 
	     uavg      = single    (length(observation.tau)), 
	     intensity = as.single (array(0, 
			    dim = c(length(observation.mu ), 
				    length(observation.tau), 
				    length(observation.phi), 
				    length(mu0)))),
             azim.av.I = as.single (array(0, 
			    dim = c(length(observation.mu ), 
				    length(observation.tau),
				    length(mu0)))))

# Errors in disort get passed back through the xerror mechanism and reported.
  xerror.summary()

# Intensity and azimuthally averaged intensity get passed back as big vectors;
# here they get reshaped into arrays with appropriate dimension names.
  results$intensity _ array(results$intensity, 
                            dim = c(length(observation.mu ), 
				    length(observation.tau), 
				    length(observation.phi), length(mu0)), 
			    dimnames = list(paste("mu" , observation.mu ), 
					    paste("tau", observation.tau), 
					    paste("phi", observation.phi),
					    paste("mu0", mu0 )) )
  results$azim.av.I _ array(results$azim.av.I, 
			    dim = c(length(observation.mu ), 
				    length(observation.tau), 
				    length(mu0)), 
			    dimnames = list(paste("mu" , observation.mu ), 
					    paste("tau", observation.tau),
					    paste("mu0", mu0 )) )
  return(list(layer.tau = results$layer.tau, 
	      ssa       = results$ssa, 
	      moments   = results$moments, 
	      mu0       = results$mu0, 
	      observation.tau = results$observation.tau, 
	      observation.mu  = results$observation.mu , 
	      observation.phi = results$observation.phi, 
	      surface.albedo  = results$surface.aledbo, 
	      direct    = results$direct, 
	      diff.dn   = results$diff.dn, 
	      diff.up   = results$diff.up, 
	      dfdt      = results$dfdt, 
	      uavg      = results$uavg, 
	      intensity = results$intensity, 
	      azim.average.I  = results$azim.av.I))
}