anchor.diff:          effective_emissivity = 1 - 0.261 * exp((-7.77e-4) * (273.16-T_air)**2)
           longwave = effective_emissivity * stefan_boltzmann * T_air**4
   
  -        longwave = longwave * (1+0.2*cloud_cover)
  +        longwave = min(longwave * (1+0.2*cloud_cover), 600.0)
       end function longwave
   
       function cloudfrac(rh,qc,nx)
  
positive.code:         !$omp do
          do j=2,ny-1
  
              T_air=theta(:,1,j)*pii(:,1,j)
              rh=relative_humidity(T_air,qv,p,j,nx)
  
              hydrometeors=qc(:,1,j)+qs(:,1,j)+qr(:,1,j)
              do k=2,nz
                  hydrometeors=hydrometeors+qc(:,k,j)+qs(:,k,j)+qr(:,k,j)
              end do
              where(hydrometeors<0) hydrometeors = 0
  
              solar_elevation=calc_solar_elevation(date,lon,j,nx,day_frac)
  
              cloud_cover(:,j) = cloudfrac(rh,hydrometeors,nx)
              swdown(:,j) = shortwave(day_frac,cloud_cover(:,j),solar_elevation,nx)
              lwdown(:,j) = longwave(T_air,cloud_cover(:,j),nx)
              ! 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)
  
  
negative.code:     use time
      implicit none
      real, allocatable, dimension(:,:) :: cos_lat_m,sin_lat_m
      integer :: nrad_layers
      real, parameter :: So=1367.0 ! Solar "constant" W/m^2
      real, parameter :: qcmin=1e-5 ! arbitrarily selected minimum "cloud" water content to affect radiation
      real, parameter :: MINIMUM_RH=1e-10 ! bounds on relative humidity = [min_rh to 1-min_rh]
  