anchor.diff:          endif
           bc%dsst_dt   = bc%dsst_dt   / nsteps
   
  -        
  +        ! parallel version didn't seem to work
   !         !$omp parallel firstprivate(ny,nsteps) &
   !         !$omp private(j) &
   !         !$omp shared(bc)
  
positive.code:                 domain%latent_heat(:,j)    = domain%latent_heat(:,j)  + bc%dlh_dt(:,j)
              endif
              
              ! these only get updated if we are using the fluxes derived from the forcing model
              if (options%physics%boundarylayer==kPBL_BASIC) then
                  domain%pbl_height(:,j) = domain%pbl_height(:,j)+ bc%dpblh_dt(:,j)
              endif
  
negative.code:         call diagnostic_update(domain,options)
      end subroutine forcing_update
      
      !>------------------------------------------------------------
      !! Update model diagnostic fields
      !!
      !! Calculates most model diagnostic fields such as Psfc, 10m height winds and ustar
      !!
      !! @param domain    Model domain data structure to be updated
      !! @param options   Model options (not used at present)
      !!
      !!------------------------------------------------------------
      subroutine diagnostic_update(domain,options)
          implicit none
          type(domain_type), intent(inout) :: domain
  