!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Polar Night Cloud Detection:
!
!	call clouds(icloud,t22,t27,t28,t31,t32,t33,isol,iclear)
!
!		icloud(5) integer	(0-no cloud, 1-cloud for ech of 5 tests
!		txx	  real		brightness temperatures for channel xx
!               isol	  integer*2	solar zenith angle to separate day (sol<90) and night (s>90)
!		iclear	  integer	(0-fail, 1-pass, i.e., clear sky, ignore icloud settings)
!
!
! Reference: Ackerman et al., 2010: Discriminating clear-sky from cloud with MODIS algoritm
!			theoretical basis document (MOD35), CIMSS, Univ. Wisconsin, 117pp.
!	     Frey et al, 2008: Cloud detection with MODIS. Part 1: Improvement in the MODIS
!			cloud mask for collection 5, J. Atmosp., Ocean. Tech., 25, 1057-1072.
!	     Liu et al., 2004: Nighttime polar cloud detection with MODIS, Rem. Sens. Env., 92,
!			181-194.
!
! Jan.-6, 2012 Andreas Muenchow, University of Delaware
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	subroutine clouds(icloud,t22,t27,t28,t31,t32,t33,isol,iclear)
	integer icloud(5),iclear(3)
	integer*2 isol
	real t22,t27,t28,t31,t32,t33,sun_angle
!
	sun_angle = isol*0.01
	temp = t31+273.15
!
! Test-1: 6.7 um (cloud is iflag=1)
!
		 call cloud1(icloud(1),t27)
!
! Test-2: 11-12 um (cloud is iflag=1)
!
		 call cloud2(icloud(2),t31,t32)
!
! Test-4: 11-3.9 um (cloud is iflag=1)
!
		 call cloud4(icloud(4),t22,t31,temp,sun_angle)
!
		 iclear(1) = 0
		 iclear(2) = 0
		 iclear(3) = 0
!
! Polar Night Only Tests:
!
		 if (sun_angle .gt. 90.) then
!
! Test-3: 3.9-12 um (cloud is iflag=1)
!
		    call cloud3(icloud(3),t22,t32,temp)
!
! Test-5: 7.2-11 um (cloud is iflag=1)
!
			call cloud5(icloud(5),t28,t31,temp)
			call clearsky(iclear,t27,t28,t31,t33)
		 end if
	return
	end
!
! Test-1: 6.7 um (cloud is iflag=1)
!
	subroutine cloud1(iflag,t27)
	integer iflag
	real t27
!
	iflag = 0
!
! Ackerman et al. (2010) give -58, -53, and -48 K for threshold
!
	threshold = -48.5
	if (t27 .le. threshold) then
	   iflag = 1
	end if
	return
	end
!
! Test-2: 11-12 um (cloud is iflag=1)
!
	subroutine cloud2(iflag,t31,t32)
	integer iflag
	real t31,t32
!
	iflag = 0
	threshold = 2 
	if (t31-t32 .ge. threshold) then
	   iflag = 1
	end if
	return
	end
!
! Test-3: 3.9-12 um (cloud is iflag=1)
!
	subroutine cloud3(iflag,t22,t32,temp)
	integer iflag
	real t22,t32,temp
!
	iflag = 0
	temp1 = 235.
	temp2 = 265.
	threshold1 = 4.5 
	threshold2 = 2.5 
	if (temp .le. temp1) then
	    threshold = threshold1
	else if (temp .ge. temp2) then
	    threshold = threshold2
	else
	    slope = (threshold1-threshold2)/(temp1-temp2)
	    threshold = threshold1 + (temp-temp1)*slope
	end if
!
	if (t22-t32 .ge. threshold) then
	   iflag = 1
	end if
	return
	end
!
! Test-4: 11-3.9 um (cloud is iflag=1)
!
	subroutine cloud4(iflag,t22,t31,temp,sun_angle)
	integer iflag
	real t22,t31,temp,sun_angle
!
	iflag = 0
! Night
	if (sun_angle .gt. 90.) then
	   temp1 = 235.
	   temp2 = 265.
!
! Frey et al. (2008) have -0.1, -0.2, -0.3 K for threshold1
! Liu et al. (2004) have -0.9 K for threshold1
! Ackerman et al (2010) have -0.8 K for threshold1
!
	   threshold1 = -0.6
!
! Frey et al. (2008) have +1.1, +1.0, +0.0 for threshold2
! Liu et al. (2004) have +0.5 K for threshold2
! Ackerman et al (2010) have +0.6 K for threshold2
!
	   threshold2 = +0.8
! Day
	else
	   temp1 = 230
	   temp2 = 245
!
! Frey et al. (2008) have -17.5, -14.5, -11.5 for threshold1
!
	   threshold1 = -11.5
!
! Frey et al. (2008) have -10.0, -10.0, -4.0 for threshold2
!
	   threshold2 = -4.0
	end if
!
	if (temp .le. temp1) then
	    threshold = threshold1
	else if (temp .ge. temp2) then
	    threshold = threshold2
	else
	    slope = (threshold1-threshold2)/(temp1-temp2)
	    threshold = threshold1 + (temp-temp1)*slope
	end if
!
	if (t31-t22 .ge. threshold .and. sun_angle.gt.90) then
	   iflag = 1
	end if
!
	if (t31-t22 .le. threshold .and. sun_angle.lt.90) then
	   iflag = 1
	end if
	return
	end
!
! Test-5: 7.2-11 um (cloud is iflag=1)
!
	subroutine cloud5(iflag,t28,t31,temp)
	integer iflag
	real t28,t31,temp
!
	iflag = 0
	   temp1 = 220.
	   temp2 = 245.
	   temp3 = 255.
	   temp4 = 265.
	   threshold1 = +2.0
	   threshold2 = -4.5
	   threshold3 = -17.5
	   threshold4 = -21.0
!
	if (temp .le. temp1) then
	    threshold = threshold1
	else if (temp .ge. temp4) then
	    threshold = threshold4
	else if (temp .ge. temp1 .and. temp.le. temp2) then
	    slope = (threshold2-threshold1)/(temp2-temp1)
	    threshold = threshold1 + (temp-temp1)*slope
	else if (temp .ge. temp2 .and. temp.le. temp3) then
	    slope = (threshold3-threshold2)/(temp3-temp2)
	    threshold = threshold2 + (temp-temp2)*slope
	else if (temp .ge. temp3 .and. temp.le. temp4) then
	    slope = (threshold4-threshold3)/(temp4-temp3)
	    threshold = threshold3 + (temp-temp3)*slope
	end if
!
	if (t28-t31 .le. threshold) then
	   iflag = 1
	end if
	return
	end
!
! Clear-Sky Restoring Tests due to Frey et al. (2008)
!				for Night Snow, Polar Night Snow
!
	subroutine clearsky(iclear,t27,t28,t31,t33)
	integer iclear(3)
!
! 6.7-11 um:
	if (t27-t31 .gt. 10.) then
	   iclear(1) = 1
	end if
! 13.3-11 um:
	if (t33-t31 .gt. 3.) then
	   iclear(2) = 1
	end if
! 7.2-11 um:
	if (t28-t31 .gt. 5.) then
	   iclear(3) = 1
	end if
	return
	end