anchor.diff:      else
         udfBc=0
       end if
  -    close(FID)
       open(FID,file='ic',action='read')
  -    read(FID,*)udfIC
  -    if(udfIc==0)then
  -      read(FID,*)pInit
  -      read(FID,*)TInit
  -      read(FID,*)uInit(1)
  -      read(FID,*)uInit(2)
  -      read(FID,*)uInit(3)
  -    else
  -      read(FID,*)tmpD
  -      iUdf(1)=int(tmpD)
  -      read(FID,*)tmpD
  -      iUdf(2)=int(tmpD)
  -      read(FID,*)tmpD
  -      iUdf(3)=int(tmpD)
  -      read(FID,*)tmpD
  -      iUdf(4)=int(tmpD)
  -      read(FID,*)tmpD
  -      iUdf(5)=int(tmpD)
  -    end if
  +      read(FID,*)udfIC
  +      if(udfIc==0)then
  +        read(FID,*)pInit
  +        read(FID,*)TInit
  +        read(FID,*)uInit(1)
  +        read(FID,*)uInit(2)
  +        read(FID,*)uInit(3)
  +      else
  +        read(FID,*)tmpD
  +        iUdf(1)=int(tmpD)
  +        read(FID,*)tmpD
  +        iUdf(2)=int(tmpD)
  +        read(FID,*)tmpD
  +        iUdf(3)=int(tmpD)
  +        read(FID,*)tmpD
  +        iUdf(4)=int(tmpD)
  +        read(FID,*)tmpD
  +        iUdf(5)=int(tmpD)
  +      end if
       close(FID)
       open(FID,file='sim',action='read')
  -    read(FID,*)tFinal
  -    read(FID,*)dt
  -    read(FID,*)tInt
  +      read(FID,*)tFinal
  +      read(FID,*)dt
  +      read(FID,*)tInt
       close(FID)
       open(FID,file='fl',action='read')
  -    read(FID,*)Rgas
  -    read(FID,*)gamm
  +      read(FID,*)Rgas
  +      read(FID,*)gamm
       close(FID)
       if(udfIc/=0.or.udfBc/=0)then
         open(FID,file='udf',action='read')
  -      call readUDFTab(FID,udf)
  +        call readUDFTab(FID,udf)
         close(FID)
       end if
  -    ! indexes of boundary conditions
  -    call mapCondTab(grid,bc,iBC)
       ! initialize algebraic solver
       call momentumEq%init(grid%nC*DIMS,momentumRHS,maa=MAXIT_MOMENTUM)
       call momentumEq%setMaxIt(MAXIT_MOMENTUM)
  
positive.code:   double precision,parameter::RTOL_MOMENTUM=1d-8 !< momentum prediction relative tolerance
    double precision,parameter::RTOL_PRESSURE=1d-8 !< pressure correction relative tolerance
    double precision,parameter::RTOL_DENSITY=1d-7 !< density equation relative tolerance
    double precision,parameter::RTOL_ENERGY=1d-8 !< energy equation relative tolerance
    
    integer,parameter::BC_WALL_TEMP=0 !< wall boundary with prescribed temperature
    integer,parameter::BC_WALL_TEMP_UDF=5 !< wall boundary with prescribed temperature by UDF
  
negative.code:         p(n)=pGst
          u(:,n)=uGst(:)
          temp(n)=TGst
          rho(n)=rho(m)!pGst/Rgas/TGst
          rhou(:,n)=-rhou(:,m)!rho(n)*uGst(:)
          rhoE(n)=rhoE(m)!rho(n)*(1d0/(gamm-1d0)*Rgas*TGst+0.5d0*dot_product(uGst,uGst))
        end if
      end do
    end subroutine
  