anchor.diff:          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 
  -        
  -        
  +        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
  
positive.code:         real, dimension(:,:), intent(in) :: lon
          integer,intent(in) :: j, nx
          real, dimension(:), intent(out) :: day_frac
          
          integer :: i
          real, dimension(nx) :: declination,day_of_year,hour_angle
          
          do i=1,nx
              day_of_year(i) = floor(calc_day_of_year(date + lon(i,j)/360.0))
              hour_angle(i) = 2*pi* mod(date+0.5 + lon(i,j)/360.0,1.0)
          end do
          day_frac=day_of_year/365.25
          
          ! fast approximation see : http://en.wikipedia.org/wiki/Position_of_the_Sun
          declination = (-0.4091) * cos(2.0*pi/365.0*(day_of_year+10))
          calc_solar_elevation = sin_lat_m(:,j) * sin(declination) + & 
                                 cos_lat_m(:,j) * cos(declination) * cos(hour_angle)
          calc_solar_elevation = asin(calc_solar_elevation)
  
          ! if the sun is below the horizon just set elevation to 0
          where(calc_solar_elevation<0) calc_solar_elevation=0
      end function calc_solar_elevation
      
      subroutine ra_simple(theta,pii,qv,qc,qs,qr,p,swdown,lwdown,cloud_cover,lat,lon,date,options,dt)
          implicit none
          real,dimension(:,:,:), intent(inout) :: theta
  
negative.code: !>----------------------------------------------------------
  !! Very simple radiation code modeled after the description in Finch and Best(2004?)
  !! of Reiff et al. 1984 shortwave and Idso and Jackson (1969) longwave.
  !!
  !! Clearsky Shortwave radiation is calculated as a function of day of year and time of day.
  !! Cloudy Shortwave is calculated as clearsky SW * f(cloud cover) [0.25-1]
  !!
  !! Cloud cover is calculated as in Xu and Randal (1996) as f(surface_RH, qc+qs+qr)
  !!
  !! Clearsky Longwave radiation is calculated as f(Tair)
  !! Cloudy longwave is scaled up by a f(cloud cover) [1-1.2]
  !!
  !! The entry point to the code is ra_simple.
  !!
  !! <pre>
  !! Call tree graph :
  