anchor.diff:  		allocate(hydrometeors(nx))
   		allocate(day_frac(nx))
   		
  -		coolingrate=1.5*(dt/86400.0) *stefan_boltzmann / 300.0 !1.5K/day radiative cooling rate (300 = W/m^2 at 270K)
  +		! 1.5K/day radiative cooling rate (300 = W/m^2 at 270K)
  +		coolingrate=1.5*(dt/86400.0) *stefan_boltzmann / 300.0 
   		
   		
   		!$omp do
  
positive.code: 		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) 
  		
  	end function shortwave
  	
  
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
  		
  