anchor.diff:                  domain%lwup=stefan_boltzmann*EMISS*domain%skin_t**4
                   RAINBL=domain%rain
                   
  -                print*, "WARNING!! enforcing surface sensible heat flux greater than 0!!"
  -                where(domain%sensible_heat<0) domain%sensible_heat=0
  -                do j=1,ny
  -                    do i=1,nx
  -                        if (domain%landmask(i,j)==kLC_LAND) then
  -                            if (domain%sensible_heat(i,j)<0) domain%sensible_heat(i,j)=0
  -                        endif
  -                    end do
  -                end do
  +!                 print*, "WARNING!! enforcing surface sensible heat flux greater than 0!!"
  +!                 where(domain%sensible_heat<0) domain%sensible_heat=0
  +!                 do j=1,ny
  +!                     do i=1,nx
  +!                         if (domain%landmask(i,j)==kLC_LAND) then
  +!                             if (domain%sensible_heat(i,j)<0) domain%sensible_heat(i,j)=0
  +!                         endif
  +!                     end do
  +!                 end do
                   
               endif
               call surface_diagnostics(domain%sensible_heat, QFX, domain%skin_t, QSFC,  &
  
positive.code:             endif
              cur_vegmonth=domain%current_month
              
              call LSM_NOAH_INIT(VEGFRAC,SNOW,SNOWC,SNOWH,domain%canopy_water,domain%soil_t,    &
                                  domain%soil_vwc, SFCRUNOFF,UDRUNOFF,ACSNOW,  &
                                  ACSNOM,domain%veg_type,domain%soil_type,     &
  
negative.code:         type(domain_type), intent(inout) :: domain
          real, intent(in) :: dt
          integer :: nx,ny
          nx = size(domain%sensible_heat,1)
          ny = size(domain%sensible_heat,2)
          
          ! convert sensible heat flux to a temperature delta term
          ! (J/(s*m^2) * s / (J/(kg*K)) => kg*K/m^2) ... /((kg/m^3) * m) => K
          dTemp=(domain%sensible_heat(2:nx-1,2:ny-1)*dt/cp)  &
               / (domain%rho(2:nx-1,1,2:ny-1)*domain%dz_inter(2:nx-1,1,2:ny-1))
          ! add temperature delta and convert back to potential temperature
          domain%th(2:nx-1,1,2:ny-1)=domain%th(2:nx-1,1,2:ny-1) + (dTemp / domain%pii(2:nx-1,1,2:ny-1))
          
          ! convert latent heat flux to a mixing ratio tendancy term
          ! (J/(s*m^2) * s / (J/kg) => kg/m^2) ... / (kg/m^3 * m) => kg/kg
          lhdQV=(domain%latent_heat(2:nx-1,2:ny-1) / LH_vaporization * dt) &
               / (domain%rho(2:nx-1,1,2:ny-1) * domain%dz_inter(2:nx-1,1,2:ny-1))
          ! add water vapor in kg/kg
          domain%qv(2:nx-1,1,2:ny-1)=domain%qv(2:nx-1,1,2:ny-1) + lhdQV
          
          ! enforce some minimum water vapor content... just in case
          where(domain%qv<SMALL_QV) domain%qv=SMALL_QV
          
      end subroutine apply_fluxes
  