anchor.diff:          where(relative_humidity>(1-MINIMUM_RH)) relative_humidity=1-MINIMUM_RH
           where(relative_humidity<MINIMUM_RH) relative_humidity=MINIMUM_RH
       end function relative_humidity
  -    
  +
       function shortwave(day_frac, cloud_cover, 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
           real, dimension(nx) :: shortwave
           real, intent(in), dimension(nx) :: day_frac, cloud_cover, solar_elevation
  -        
  +
           real, dimension(nx) :: sin_solar_elev
           integer, intent(in) :: nx
  -        
  +
           sin_solar_elev = sin(solar_elevation)
           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
  -    
  +
       function longwave(T_air, cloud_cover,nx)
   !       compute longwave down at the surface based on air temperature and cloud fraction
   !       based on Idso and Jackson (1969)
  
positive.code:         real :: coolingrate
          integer :: nx,ny,j,k,nz
          real, allocatable, dimension(:) :: rh,T_air,solar_elevation, hydrometeors,day_frac
          
          
          !$omp parallel private(nx,ny,nz,j,rh,T_air,solar_elevation,hydrometeors,day_frac,coolingrate) &
          !$omp shared(theta,pii,qv,p,qc,qs,qr,date,lon,cloud_cover,swdown,lwdown)
          nx=size(lat,1)
          ny=size(lat,2)
          nz=size(qv,2)
          
          allocate(rh(nx))
          allocate(T_air(nx))
          allocate(solar_elevation(nx))
          allocate(hydrometeors(nx))
          allocate(day_frac(nx))
          
          ! 1.5K/day radiative cooling rate (300 = W/m^2 at 270K)
          coolingrate=1.5*(dt/86400.0) *stefan_boltzmann / 300.0 
          
          
          !$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)
              
          end do
          !$omp end do
          
          deallocate(rh,T_air,solar_elevation, hydrometeors, day_frac)
          !$omp end parallel
          
      end subroutine ra_simple
  end module module_ra_simple
  
negative.code:         real,intent(in),dimension(nx)::rh,qc
          integer, intent(in) :: nx
          real,dimension(nx)::temporary
  
          temporary=((1-rh)*qc)**0.25
          where(temporary >1) temporary=1
          where(temporary <0.0001) temporary=0.0001
  
          cloudfrac=qc-qcmin
          where(cloudfrac<0) cloudfrac=0
  
          cloudfrac=(rh**0.25) * (1-exp((-2000*(cloudfrac)) / temporary))
          where(cloudfrac<0) cloudfrac=0
          where(cloudfrac>1) cloudfrac=1
  
      end function
  
      function calc_solar_elevation(date,lon,j,nx,day_frac)
          implicit none
          real, dimension(nx) :: calc_solar_elevation
  