anchor.diff:  		shortwave=So * (1+0.035*cos(day_frac * 2*pi)) * sin_solar_elev * (0.48+0.29*sin_solar_elev)
   		
   		! note it is cloud_cover**3.4 in Reiff, but this makes almost no difference and integer powers are fast
  -		shortwave=shortwave * (1 - 0.75 * cloud_cover**3) 
  +		shortwave=shortwave * (1 - (0.75 * (cloud_cover**3)) )
   		
   	end function shortwave
   	
  
positive.code: 		where(relative_humidity<MINIMUM_RH) relative_humidity=MINIMUM_RH
  	end function relative_humidity
  	
  	function shortwave(cloud_cover, day_frac, solar_elevation,nx)
  ! 		compute shortwave down at the surface based on solar elevation, fractional day of the year, and cloud fraction
  ! 		based on Reiff et al. (1984)
  		implicit none
  
negative.code: 			! apply a simple radiative cooling to the atmosphere
  			theta(2:nx-1,:,j)=theta(2:nx-1,:,j) - (((theta(2:nx-1,:,j)*pii(2:nx-1,:,j))**4) * coolingrate)
  			
  		end do
  		!$omp end do
  		
  