anchor.diff:      double precision::uAvg(DIMS),HAvg,cAvg
       
       call grid%up()
  -    if(.not.(allocated(JacP)))then
  +    if(.not.allocated(JacP))then
         allocate(JacP(10,10,grid%nP))
       end if
  -    if(.not.(allocated(JacC)))then
  +    if(.not.allocated(JacC))then
         allocate(JacC(5,size(grid%neib,1)*5+5,grid%nC))
       end if
       JacP(:,:,:)=0d0
  
positive.code:     double precision,allocatable::gradF(:,:,:)
      
      call grid%up()
      if(.not.(allocated(adv)))then
        allocate(adv(size(s,1),grid%nC))
      end if
      adv(:,:)=0d0
  
negative.code:     double precision,intent(in)::f(:,:,:) !< fluxes
      double precision,allocatable,intent(inout)::adv(:,:) !< advection output
      
      if(.not.allocated(adv))then
        allocate(adv(size(s,2),grid%nC))
      end if
      adv(:,:)=0d0
  