anchor.diff:      double precision,allocatable::YInit(:)
       integer::udfIc,udfBc,iUdf(5)
       
  -    ! read inputs
  +    ! read input files
       open(FID,file='grid.msh',action='read')
  -    call readGMSH(FID,grid)
  +      call readGMSH(FID,grid)
  +    close(FID)
       call grid%up()
  +    open(FID,file='bc',action='read')
  +      call readCondTab(FID,bc)
       close(FID)
  -    open(FID,file='bc',action='read')
  -    call readCondTab(FID,bc)
  +    call mapCondTab(grid,bc,iBC)
       if(any(bc%t(:)==BC_WALL_TEMP_UDF).or.&
       &  any(bc%t(:)==BC_WALL_FLUX_UDF).or.&
       &  any(bc%t(:)==BC_IN_STATIC_UDF).or.&
  
positive.code:         p(n)=pGst
          u(:,n)=uGst(:)
          temp(n)=TGst
          rho(n)=pGst/Rgas/TGst
          rhou(:,n)=rho(n)*uGst(:)
          rhoE(n)=rho(n)*(1d0/(gamm-1d0)*Rgas*TGst+0.5d0*dot_product(uGst,uGst))
        end if
      end do
    end subroutine
  
negative.code:     end if
      
      ! solution and residual scales
      pScale=max(maxval(p),maxval(0.5d0*rho*norm2(u,1)**2))
      uScale=max(maxval(norm2(u,1)),sqrt(2d0*(maxval(p)-minval(p))/minval(rho)),0.01d0)
      tempScale=maxval(temp)
      rhoScale=maxval(rho)
      rhouScale=rhoScale*uScale
      rhoEScale=max(rhoScale*Rgas*tempScale,0.5d0*rhoScale*uScale**2)
      
      ! set tolerance and maximum number of iterations, and etc.
      call momentumEq%setTol(rhouScale*RTOL_MOMENTUM)
      call pressureEq%setTol(pScale*RTOL_PRESSURE)
      call densityEq%setTol(rhoScale*RTOL_DENSITY)
      call energyEq%setTol(rhoEScale*RTOL_ENERGY)
  