function [csza] = cossza(glat,glon,sbsllat,sbsllon)

%   Computes cosine of solar zenith angle from geographic coordinates of point
%   in question and of subsolar point.
%   940730 A. D. Richmond, NCAR

%    Inputs:
%    	glat = geographic latitude of point (degrees)
%    	glon = geographic longitude of point (degrees)
%    	sbsllat = geographic latitude of subsolar point (degrees)
%    	sbsllon = geographic longitude of subsolar point (degrees)
%    Output:
%    	csza = cosine of solar zenith angle

%Constants:
rtod = 5.72957795130823e1;
dtor = 1.745329251994330e-2;

%Caluculating cosine of solar zenith angle:

      ct   = sin(glat*dtor);
      st   = sqrt(1. - ct*ct);
      cs   = sin(sbsllat*dtor);
      ss   = sqrt(1. - cs*cs);
      ang  = (glon - sbsllon)*dtor;
      csza = ct*cs + st*ss*cos(ang);
