anchor.diff:                  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 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, &
  
positive.code:             
              ! exchange coefficients
              windspd=sqrt(domain%u10**2+domain%v10**2)
              call calc_exchange_coefficient(windspd,domain%skin_t,domain%T2m,CHS)
              CHS2=CHS
              CQS2=CHS
              
  
negative.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
          real, dimension(:,:,:),intent(inout) :: airt
          real,dimension(:,:),intent(inout) :: exchange_C
  
          ! Richardson number
          where(wind==0) wind=1e-10
  
          Ri = gravity/airt(:,1,:) * (airt(:,1,:)-tskin)*z_atm/(wind**2)
  
  !         print*,"--------------------------------------------------"
  !         print*, "Surface Richardson number"
  !         print*, Ri(128,103), airt(128,1,103), tskin(128,103), z_atm(128,103), wind(128,103)
          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))
  !         print*, exchange_C(128,103), lnz_atm_term(128,103), base_exchange_term(128,103)
  
          where(exchange_C > MAX_EXCHANGE_C) exchange_C=MAX_EXCHANGE_C
      end subroutine calc_exchange_coefficient
  
  
  ! eqn A11 in Appendix A.2 of Chen et al 1997 (see below for reference)
      subroutine F2_formula(F2, z_atm, zo, Ri)
          real, dimension(:,:), intent(inout) :: F2, z_atm, zo, Ri
          
          ! for the stable case from Mahrt (1987)
          where(Ri>=0) F2=exp(-Ri)
          ! for the unstable case from Holtslag and Beljaars (1989)
          where(Ri<0)  F2=(1-(15*Ri)/(1+((70.5*karman**2 * sqrt(-Ri * z_atm/zo))/(lnz_atm_term**2))) )
          
      end subroutine F2_formula
  !From Appendix A.2 in Chen et al 1997 
  ! Impact of Atmospheric Surface-layer Parameterizations in the new Land-surface Scheme of the Ncep Mesoscale ETA Model
  ! Boundary-Layer Meteorology 85:391-421
      subroutine calc_mahrt_holtslag_exchange_coefficient(wind,tskin,airt,znt, exchange_C)
          implicit none
          real, dimension(:,:),intent(inout) :: wind,tskin
          real, dimension(:,:,:),intent(inout) :: airt
          real,dimension(:,:),intent(inout) :: exchange_C, znt
      
          ! Richardson number 
          where(wind==0) wind=1e-10
          Ri = gravity/airt(:,1,:) * (airt(:,1,:)-tskin)*z_atm/(wind**2)
          
          call F2_formula(base_exchange_term, z_atm,znt,Ri)
          
          exchange_C = karman**2 * base_exchange_term / lnz_atm_term**2
          
          where(exchange_C > MAX_EXCHANGE_C) exchange_C=MAX_EXCHANGE_C
      end subroutine calc_mahrt_holtslag_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
  