anchor.diff:        !        points in order to be able to add other lines, it is a
         !        good idea to make NFODIM quite a bit larger
   
  -!      COMMON // R(NDDIM),ENTOT(NDDIM),T(NDDIM)
  -!     $ ,XJC(NDDIM),XJCARR(NDDIM,NFDIM),XJL(NDDIM,MAXIND)
  -!     $ ,EDDI(3,NDDIM),EDDARR(3,NDDIM,NFDIM),TAUROSS(NDDIM)
  -!     $ ,RNE(NDDIM),VELO(NDDIM),GRADI(NDDIM)
  -!     $ ,XLAMBDA(NFDIM),FWEIGHT(NFDIM),EMFLUX(NFDIM),AKEY(NFDIM)
  -!     $ ,ENLTE(NDIM)
  -!     $ ,P(NPDIM),Z(NDDIM,NPDIM),POPNUM(NDDIM,NDIM)
  -!     $ ,HTOT(NDDIM),GTOT(NDDIM),XTOT(NDDIM),ETOT(NDDIM)
  -!     $ ,POP1(NDDIM,NDIM),POP2(NDDIM,NDIM),POP3(NDDIM,NDIM)
  -!     $ ,ETA(NDDIM),OPA(NDDIM),THOMSON(NDDIM),TAUTHOM(NDDIM)
  -!     $ ,A(NPDIM),B(NPDIM,NPDIM),C(NPDIM),W(NPDIM)
  -!     $ ,BX(NPDIM,NPDIM,NDDIM),WX(NPDIM,NDDIM)
  -!     $ ,U(NDDIM,NPDIM)
  -
  -!     $ VDU(NDDIM)
  -!     $ ,ZRAY(NDADDIM),XCMF(NDADDIM),RRAY(NDADDIM),WLK(NDDIM,NPDIM)
  -!     $ ,Tion_pot(nddim,3),dil(nddim)
  -
  -!     $ ,ITNE(NDDIM),NFEDGE(NDIM)
  -!     $ ,IWARN(NDDIM),LEVELPL(NDIM)
  -
         real*8, allocatable, dimension(:) ::       R
         real*8, allocatable, dimension(:) ::       entot, enlte
         real*8, allocatable, dimension(:) ::       T
  
positive.code:       NP = read_param('NP')
        NF = read_param('NF')
  
        print*, 'fioss here: ', ND, NP, NF
  
        stop
  
        IFL = 3; open(IFL, file = 'MODFILE', STATUS='OLD'); rewind IFL
  
  
negative.code:       CALL READMOD(IFL,N,ND,TEFF,R,NP,P,Z,ENTOT,VELO,GRADI,RSTAR,VDOP,NF,
       $             XLAMBDA(1 : NF),FWEIGHT(1 : NF),AKEY(1 : NF),
       $             ABXYZ,NATOM,MODHEAD,JOBNUM,LBLANK)
  
        close(IFL)
  
        IFL = 3; open(IFL, file = 'POPNUM', STATUS='OLD')
        !***MH readpop also reads number of depth points, temperature, rne
  