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:     integer,allocatable, dimension(:,:) :: IVGTYP,ISLTYP
      integer :: ITIMESTEP, update_interval, cur_vegmonth
      
      real, parameter :: kappa=0.4
      real, parameter :: freezing_threshold=273.15
      real, parameter :: SMALL_PRESSURE=0.1 !note: 0.1Pa is very small 1e-10 wouldn't affect a single-precision float
      real, parameter :: SMALL_QV=1e-10
  
negative.code:                 endif
                  
                  ! update T2m and Q2m prior to LSM call to make sure everything is in sync
  !                 call surface_diagnostics(domain%sensible_heat, QFX, domain%skin_t, QSFC,  &
  !                                          CHS2, CQS2,domain%T2m, domain%Q2m, domain%psfc)
                                           
                  call lsm_noah(domain%dz_inter,domain%qv,domain%p_inter,domain%th*domain%pii,domain%skin_t,  &
                              domain%sensible_heat,QFX,domain%latent_heat,domain%ground_heat, &
  