anchor.diff:                  ! accumulate soil moisture over the entire column
                   domain%soil_totalmoisture%data_2d = domain%soil_water_content%data_3d(:,1,:) * DZS(1) * 1000
                   do i = 2,num_soil_layers
  -                    domain%soil_totalmoisture%data_2d = domain%soil_totalmoisture%data_2d + domain%soil_water_content%data_3d(:,i,:) * DZS(i)
  +                    domain%soil_totalmoisture%data_2d = domain%soil_totalmoisture%data_2d + domain%soil_water_content%data_3d(:,i,:) * DZS(i) * 1000
                   enddo
   
                   ! 2m Air T and Q are not well defined if Tskin is not coupled with the surface fluxes
  -                call surface_diagnostics(domain%sensible_heat%data_2d,    &
  -                                         QFX,                             &
  -                                         domain%skin_temperature%data_2d, &
  -                                         QSFC,                            &
  -                                         CHS2,                            &
  -                                         CQS2,                            &
  -                                         domain%temperature_2m%data_2d,   &
  -                                         domain%humidity_2m%data_2d,      &
  -                                         domain%surface_pressure%data_2d)
  +                call surface_diagnostics(domain%sensible_heat%data_2d,          &
  +                                         QFX,                                   &
  +                                         domain%skin_temperature%data_2d,       &
  +                                         QSFC,                                  &
  +                                         CHS2,                                  &
  +                                         CQS2,                                  &
  +                                         domain%temperature_2m%data_2d,         &
  +                                         domain%humidity_2m%data_2d,            &
  +                                         domain%surface_pressure%data_2d,       &
  +                                         VEGFRAC,                               &
  +                                         domain%veg_type,                       &
  +                                         domain%land_mask,                      &
  +                                         domain%temperature_2m_veg%data_2d,     &
  +                                         domain%temperature_2m_bare%data_2d,    &
  +                                         domain%mixing_ratio_2m_veg%data_2d,    &
  +                                         domain%mixing_ratio_2m_bare%data_2d)
  +
               endif
           endif
           if (options%physics%landsurface>0) then
  
positive.code:         where(exchange_C < MIN_EXCHANGE_C) exchange_C=MIN_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
          implicit none
          REAL, DIMENSION(ims:ime, jms:jme ), INTENT(IN)    ::  HFX, QFX, TSK, QSFC
          REAL, DIMENSION(ims:ime, jms:jme ), INTENT(INOUT) ::  Q2, T2
          REAL, DIMENSION(ims:ime, jms:jme ), INTENT(IN)    ::  PSFC, CHS2, CQS2
          integer :: i,j, nx,ny
          real :: rho
  
          !$omp parallel default(shared), private(nx,ny,i,j,rho)
          nx=size(HFX,1)
          ny=size(HFX,2)
          !$omp do
          do j=jts,jte
              do i=its,ite
                  RHO = PSFC(I,J)/(Rd * TSK(I,J))
                  if(CQS2(I,J).lt.1.E-3) then
                     Q2(I,J) = QSFC(I,J)
                  else
                     Q2(I,J) = QSFC(I,J) - QFX(I,J)/(RHO*CQS2(I,J))
                  endif
                  if(CHS2(I,J).lt.1.E-3) then
                     T2(I,J) = TSK(I,J)
                  else
                     T2(I,J) = TSK(I,J) - HFX(I,J)/(RHO*CP*CHS2(I,J))
                  endif
                  ! TH2(I,J) = T2(I,J)*(1.E5/PSFC(I,J))**ROVCP
              enddo
          enddo
  
negative.code:                          kVARS%sensible_heat, kVARS%latent_heat, kVARS%u_10m, kVARS%v_10m, kVARS%temperature_2m,        &
                           kVARS%humidity_2m, kVARS%surface_pressure, kVARS%longwave_up, kVARS%ground_heat_flux,          &
                           kVARS%soil_totalmoisture, kVARS%soil_deep_temperature, kVARS%roughness_z0, kVARS%ustar,        &
                           kVARS%snow_height, kVARS%lai, kVARS%temperature_2m_veg,                                        &
                           kVARS%veg_type, kVARS%soil_type, kVARS%land_mask])
  
               call options%advect_vars([kVARS%potential_temperature, kVARS%water_vapor])
  