anchor.diff:      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]
  
positive.code:         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)
      end function longwave
  
      function cloudfrac(rh,qc,nx)
  
negative.code:         cos_lat_m=cos(domain%lat/360.0 * 2*pi)
          sin_lat_m=sin(domain%lat/360.0 * 2*pi)
  
          nrad_layers = 5
  
      end subroutine ra_simple_init
  
  !   day_of_year=calc_day_of_year(date)
  !   time_of_day=calc_time_of_day(date)
  
  
      function relative_humidity(t,qv,p,nx)
          implicit none
          real,               dimension(nx) :: relative_humidity
          real,   intent(in), dimension(nx) :: t
          real,   intent(in), dimension(nx)  :: qv, p
          integer,intent(in) :: nx
          real,               dimension(nx) :: mr, e, es
  
          ! convert sensible humidity to mixing ratio
          mr = qv / (1-qv)
          ! convert mixing ratio to vapor pressure
          e = mr * p / (0.62197+mr)
          ! convert temperature to saturated vapor pressure
          es = 611.2 * exp(17.67 * (t - 273.15) / (t - 29.65))
          ! finally return relative humidity
          relative_humidity = e / es
          ! because it is an approximation things could go awry and rh outside or reasonable bounds could break something else.
          ! alternatively air could be supersaturated (esp. on boundary cells) but cloud fraction calculations will break.
          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)
  