anchor.diff:                          endif
                       endif
                   endif
  +                if (Q2(i,j) < SMALL_QV) Q2(i,j) = SMALL_QV
   
                   ! TH2(I,J) = T2(I,J)*(1.E5/PSFC(I,J))**ROVCP
               enddo
  
positive.code:         !$omp end parallel
      end subroutine surface_diagnostics
  
      subroutine apply_fluxes(domain,dt)
          ! add sensible and latent heat fluxes to the first atm level
          implicit none
          type(domain_t), intent(inout) :: domain
          real, intent(in) :: dt
  
          associate(density       => domain%density%data_3d,       &
                    sensible_heat => domain%sensible_heat%data_2d, &
                    latent_heat   => domain%latent_heat%data_2d,   &
                    dz            => domain%dz_interface%data_3d,  &
                    pii           => domain%exner%data_3d,         &
                    th            => domain%potential_temperature%data_3d, &
                    qv            => domain%water_vapor%data_3d    &
              )
  
          ! 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=(sensible_heat(its:ite,jts:jte) * dt/cp)  &
               / (density(its:ite,kts,jts:jte) * dz(its:ite,kts,jts:jte))
          ! add temperature delta and convert back to potential temperature
          th(its:ite,kts,jts:jte) = th(its:ite,kts,jts:jte) + (dTemp / pii(its:ite,kts,jts:jte))
  
          ! 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=(latent_heat(its:ite,jts:jte) / LH_vaporization * dt) &
               / (density(its:ite,kts,jts:jte) * dz(its:ite,kts,jts:jte))
          ! add water vapor in kg/kg
          qv(its:ite,kts,jts:jte) = qv(its:ite,kts,jts:jte) + lhdQV
  
          ! enforce some minimum water vapor content... just in case
          where(qv(its:ite,kts,jts:jte) < SMALL_QV) qv(its:ite,kts,jts:jte) = SMALL_QV
  
          end associate
  
  
negative.code:                              domain%cosine_zenith_angle%data_2d,       &
                               domain%latitude%data_2d,                  &
                               domain%longitude%data_2d,                 &
                               domain%dz_interface%data_3d/2,              &
                               lsm_dt,                                   &
                               DZS,                                      &
                               num_soil_layers,                          &
  