anchor.diff:          allocate(Z0(ime,jme))
           Z0=0.01 !?
           allocate(QSFC(ime,jme))
  -        QSFC=0
  +        QSFC=domain%qv(:,1,:)
           allocate(Ri(ime,jme))
           Ri=0
           allocate(z_atm(ime,jme))
  
positive.code:         type(domain_type), intent(inout) :: domain
          real, intent(in) :: dt
          integer :: nx,ny
          nx=ime
          ny=jme
          ! 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. 
          where(domain%qv<SMALL_QV) domain%qv=SMALL_QV
          
      end subroutine apply_fluxes
  
negative.code:         do j=1,ny
              do i=1,nx
                  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))
  