%
%	In the CERES SYN Ed. 4 flux calculations, the layer cloud fractions (defined by the postion of cloud top) 
% are overlapped (under certain conditions) such that the cloud base used in the flux calculations is not 
% necessarily the cloud base associated with each layer.  For example, if a high level cloud is present with a 
% Cloud Top at 200 hPA and a Cloud base at 300 hPA and a CF = 0.5 (exposed to space) the flux calculations 
% will include profile with a high cloud at 200 hPA covering have the domain, BUT cloud base may not be at 
% 300 hPA, it may be lower.  Specicifially if the overlap of an upper layer with a lower layers occurs more often
% than single-layer conditions for the upper cloud, the cloud base will be set to the base of the lowest layer with 
% the largest overlap, as long as the upper layer OD is > 6.
%
%	The inputs are arrays of length 4, where:
%
%		SYN_CF = the four SYN cloud fractions (CF) for each of the 4 SYN layers. Layers with no cloud should have CF==0;
%	     SYN_CB = the cloud base (CB) associated with each SYN layer
%		SYN_OD = the optical depth (OD) associated with each SYN layer.  
% 
%	where SYN_CF(1) = high cloud fracitonj and SYN_CF(4) each low cloud fraction.
%
function [Overlap_CB,Overlap_index]=SYN_cloud_bases_after_overlap(SYN_CF,SYN_CB,SYN_OD)

	% initialize outputs 

		Overlap_CB    = NaN*ones(1,4);		% Cloud Base values used in flux caculations 
		Overlap_index = NaN*ones(1,4);		% index which indicates which overlap possibility is most common (if any overlap is used).

	% check input array lengths

		if(length(SYN_CF)~=4 | length(SYN_CB)~=4 | length(SYN_OD)~=4)
			disp('SYN input arrays should be arrays of length 4')
			return
		end


	% check total input CF is less than = 1

		if(sum(SYN_CF)>1)
			disp('sum of SYN Layer Cloud Fractions must be 1 or less')
			return
		end


	% calculate the "True" fraction (T) for each layer given the SYN fractions are the fraction exposed to space, S (assuming radome overla)

		T(1) = SYN_CF(1);
		T(2) = SYN_CF(2) / (1 - SYN_CF(1));
		T(3) = SYN_CF(3) / (1 - SYN_CF(1) - SYN_CF(2));
		T(4) = SYN_CF(4) / (1 - SYN_CF(1) - SYN_CF(2) - SYN_CF(3));


	% based on the True cloud fraction (T) for each of the 16 possible cloud overlaps (c) (since clouds can exist in up to 4 layers) and associated cloud base (cb) 

		% overlap with a high layer cloud on top

		c(1) = T(1)*T(2)*T(3)*T(4);				
		c(2) = T(1)*T(2)*T(3)*(1-T(4));
		c(3) = T(1)*T(2)*(1-T(3))*T(4);
		c(4) = T(1)*T(2)*(1-T(3))*(1-T(4));
		c(5) = T(1)*(1-T(2))*T(3)*T(4);
		c(6) = T(1)*(1-T(2))*T(3)*(1-T(4));
		c(7) = T(1)*(1-T(2))*(1-T(3))*T(4);
		c(8) = T(1)*(1-T(2))*(1-T(3))*(1-T(4));
		
		cb(1) = SYN_CB(4);
		cb(2) = SYN_CB(3);
		cb(3) = SYN_CB(4);
		cb(4) = SYN_CB(2);
		cb(5) = SYN_CB(4);
		cb(6) = SYN_CB(3);
		cb(7) = SYN_CB(4);
		cb(8) = SYN_CB(1);

			
		% overlap with a high-mid cloud layer on top
		
		c(9)  = (1-T(1))*T(2)*T(3)*T(4);
		c(10) = (1-T(1))*T(2)*T(3)*(1-T(4));
		c(11) = (1-T(1))*T(2)*(1-T(3))*T(4);
		c(12) = (1-T(1))*T(2)*(1-T(3))*(1-T(4));
		
		cb(9)  = SYN_CB(4);
		cb(10) = SYN_CB(3);
		cb(11) = SYN_CB(4);
		cb(12) = SYN_CB(2);

		% overlap with a low-mid cloud layer on top
		
		c(13)  = (1-T(1))*(1-T(2))*T(3)*T(4);
		c(14)  = (1-T(1))*(1-T(2))*T(3)*(1-T(4));

		cb(13) = SYN_CB(4);
		cb(14) = SYN_CB(3);
			

		% only low level cloud

		c(15) = (1-T(1))*(1-T(2))*(1-T(3))*T(4);

		cb(15) = SYN_CB(4);
		

	% check the above yields the correct total CF

		if( abs(sum(SYN_CF) - sum(c)) >1E-6)

			disp('Coding Error or Invalid SYN_CF input')
			return
		end

	% Now select final four categories.   Final categories used in the flux calculations are required to include 1 category with high, high-mid, mid-low, and low top, respectively and 

		% high cloud

		if(SYN_CF(1)>0)

			if(SYN_OD(1)<6)

					% ignore overlap
					Overlap_CB(1) = SYN_CB(1);
			else
					% find overlap type with largest cloud fraction
					max_cf = max(c(1:8));
					ind = find(c(1:8)==max_cf);

					% select the first type (if two or more are tied)
					Overlap_index(1) = ind(1);

					% use the cloud base for the this overlay type
					Overlap_CB(1) = cb(Overlap_index(1));
			end

		end


		% high-mid cloud

		if(SYN_CF(2)>0)

			if(SYN_OD(2)<6)

					% ignore overlap
					Overlap_CB(2) = SYN_CB(2);
			else
					% find overlap type with largest cloud fraction
					max_cf = max(c(9:12));
					ind = find(c(9:12)==max_cf);

					% select the first type (if two or more are tied)
					Overlap_index(2) = 8 + ind(1);

					% use the cloud base for the this overlay type
					Overlap_CB(2) = cb(Overlap_index(2));
			end

		end


		% mid-low cloud

		if(SYN_CF(3)>0)

			if(SYN_OD(3)<6)

					% ignore overlap
					Overlap_CB(3) = SYN_CB(3);
			else
					% find overlap type with largest cloud fraction
					max_cf = max(c(13:14));
					ind = find(c(13:14)==max_cf);

					% select the first type (if two or more are tied)
					Overlap_index(3) = 12 + ind(1);

					% use the cloud base for the this overlay type
					Overlap_CB(3) = cb(Overlap_index(3));
			end
		end

		% low-only cloud

		if(SYN_CF(4)>0 )

			Overlap_index(4) = 15;
			Overlap_CB(4) = cb(Overlap_index(4));

		end
return

		
