anchor.diff:              
               ! exchange coefficients
               windspd=sqrt(domain%u10**2+domain%v10**2)
  -            call calc_exchange_coefficient(windspd,domain%skin_t,domain%T2m,CHS)
  +            if (exchange_term==1) then
  +                call calc_exchange_coefficient(windspd,domain%skin_t,domain%T,CHS)
  +            elseif (exchange_term==2) then
  +                call calc_mahrt_holtslag_exchange_coefficient(windspd,domain%skin_t,domain%T,domain%znt,CHS)
  +            endif
  +!             print*, CHS(128,103)
               CHS2=CHS
               CQS2=CHS
               
  
positive.code:         !from : http://www.srh.noaa.gov/images/epz/wxcalc/mixingRatio.pdf
          sat_mr=0.6219907*e_s/(p-e_s) !(kg/kg)
      end function sat_mr
      
      subroutine calc_exchange_coefficient(wind,tskin,airt,exchange_C)
          implicit none
          real, dimension(:,:),intent(inout) :: wind,tskin,airt
          real,dimension(:,:),intent(inout) :: exchange_C
          
          ! Richardson number 
          where(wind==0) wind=1e-10
          
          Ri = gravity/airt * (airt-tskin)*z_atm/(wind**2)
          
          where(Ri<0)  exchange_C=lnz_atm_term * (1.0-(15.0*Ri)/(1.0+(base_exchange_term * sqrt((-1.0)*Ri))))
          where(Ri>=0) exchange_C=lnz_atm_term * 1.0/((1.0+15.0*Ri)*sqrt(1.0+5.0*Ri))
          
          where(exchange_C > MAX_EXCHANGE_C) exchange_C=MAX_EXCHANGE_C
      end subroutine calc_exchange_coefficient
      
      subroutine surface_diagnostics(HFX, QFX, TSK, QSFC, CHS2, CQS2,T2, Q2, PSFC)
          ! taken almost directly / exactly from WRF's module_sf_sfcdiags.F
  
negative.code:                 
                  ! now that znt has been updated, we need to recalculate terms
                  lnz_atm_term = log((z_atm+domain%znt)/domain%znt)
                  if (exchange_term==1) then
                      base_exchange_term=(75*karman**2 * sqrt((z_atm+domain%znt)/domain%znt)) / (lnz_atm_term**2)
                      lnz_atm_term=(karman/lnz_atm_term)**2
                  endif
                                  
                  ! note this is more or less just diagnostic and could be removed
                  domain%lwup=stefan_boltzmann*EMISS*domain%skin_t**4
                  RAINBL=domain%rain
                  
  !                 i=128
  !                 j=103
  !                 print*,"            ---------------------------"
  !                 print*, "   soil_t ", "       skin_t ", "         T2m"
  !                 print*, domain%soil_t(i,1,j), domain%skin_t(i,j), domain%T2m(i,j)
  !                 print*, "   Tair ", "        sensible ", "        CHS"
  !                 print*, domain%pii(i,1,j)*domain%th(i,1,j), domain%sensible_heat(i,j), CHS(i,j)
  !                 print*, "    SWD  ", "          LWD ", "         LWU "
  !                 print*, domain%swdown(i,j), domain%lwdown(i,j), domain%lwup(i,j)
  !                 print*, "ln z term,        base exchange term        wind"
  !                 print*, lnz_atm_term(i,j), base_exchange_term(i,j), windspd(i,j)
  
  !                 print*, "WARNING!! enforcing surface sensible heat flux greater than 0!!"
  !                 where(domain%sensible_heat<0) domain%sensible_heat=0
  !                 do j=1,ny
  