anchor.diff:      end if
       
       ! solution and residual scales
  -    pScale=max(maxval(p)-minval(p),maxval(0.5d0*rho*norm2(u,1)**2),1d0)
  -    uScale=max(maxval(norm2(u,1)),sqrt(2d0*pScale/minval(rho)))
  -    tempScale=max(maxval(temp)-minval(temp),pScale/minval(rho)/Rgas)
  +    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(rhoScale*RTOL_MOMENTUM)
  +    call momentumEq%setTol(rhouScale*RTOL_MOMENTUM)
       call pressureEq%setTol(pScale*RTOL_PRESSURE)
       call densityEq%setTol(rhoScale*RTOL_DENSITY)
       call energyEq%setTol(rhoEScale*RTOL_ENERGY)
  
positive.code:     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
      close(FID)
      open(FID,file='sim',action='read')
      read(FID,*)tFinal
      read(FID,*)dt
      read(FID,*)tInt
      close(FID)
      open(FID,file='fl',action='read')
      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)
        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)
  
negative.code:     double precision,allocatable::YInit(:)
      integer::udfIc,udfBc,iUdf(5)
      
      ! read input files
      open(FID,file='grid.msh',action='read')
        call readGMSH(FID,grid)
      close(FID)
      call grid%up()
      open(FID,file='bc',action='read')
        call readCondTab(FID,bc)
      close(FID)
      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.&
  