program MYOWNDGD

       IMPLICIT NONE
       INTEGER MXGAM, DIMX
       PARAMETER (MXGAM=60) !option du type de matrice
       PARAMETER (DIMX=10000) !taille max de la matrice sous forme vecteur
       DOUBLE PRECISION X(DIMX),GAMMA(MXGAM),DET !X matrice var-cov, type matrice, determinant
       CHARACTER *72 BASENAME ! nom fichier
       INTEGER NV,NGAM,OPT,IFRM !NGAM=nb parametres NV=nb niveau var matrice OPT=option sortie .own et .gdg

       INTEGER I,J,K,LL,LENBN
       DOUBLE PRECISION C0,C1,C2

!modif (david,2013) pour faire passer modele anedependant
	Integer :: NGAMMAANTE,col,rtest,jd,datatemps,nbcar !nb de para en tout	
	Integer, allocatable ::  temps(:) !vecteur de stockage des valeurs de temps
	Integer, allocatable :: tante(:),NGAMADM(:),oante(:),ovar(:) !ordre antedep, ordre var 
	DOUBLE PRECISION, allocatable :: GAMMAANTE(:)
	character :: tmodel*40,file*80,BASEFILE*40
	integer, allocatable :: typepara(:),typeorder(:)
	integer, allocatable :: typeKr(:),typeKrant(:)
integer,parameter::   io_p=40          ! unit number for parameter file

!modif(david,2016), pour seconde vaiante antedep multicaractere
	Integer, allocatable :: tantecross(:),oantecross(:),typecross(:),mattricross(:,:)
	Integer, allocatable :: debantecross(:),deltaantecross(:)
	integer :: valcross
	
	
!modif(david,2016), pour rajouter une correlation initiale entre les caracteres au temps 0
	Integer,allocatable :: prescorrinit(:),mattriinit(:,:),prescorrinit2(:)
	DOUBLE PRECISION , allocatable:: corrinit(:)	

!modif pour faire multi ou uni
 integer :: n,xx(40)
 character :: xc(40)*40
Integer :: oanteuni,ovaruni 

!      ASREML writes  BASENAME.own
!             runs    MYOWNGDG BASENAME
!
!             reads   BASENAME.gdg
! ------------------------------------------
!      First obtain BASENAME
! Unix & DIgital Fortran
       CALL GETARG(1,BASENAME)  !permet recup de l argument 1 ecrit apres executable pour creer les .gdg et savoir quel .own lire
! Lahey
!      CALL GETCL(BASENAME)
       
       LENBN=INDEX(BASENAME,' ')

! ------------------------------------------
!      NOW READ   BASENAME.own
      BASENAME(LENBN:)='.own'

      OPEN(21,file=BASENAME,status='old',err=920)
      READ(21,21,ERR=930,END=930) NV,NGAM,OPT !NV calcule tout seul par asreml, NGAM=chiffre apres OWN, OPT, donnee apres !F
21   FORMAT(4I5)

      IF(NGAM.GT.MXGAM) GOTO 930
 
BASEFILE='para'
if (OPT.ne.1) then
BASEFILE(5:)=str(OPT)
endif
 
!ma partie avec fonction speciale
!---------------------------------
!---------------------------------
open(io_p,file=BASEFILE)
call readline(io_p,xx,xc,n)
print*,xc(1)
if (xc(1).eq.'NBCHARACTER') then
      print*,'multiple trait model'
else if (xc(1).eq.'MODEL') then
      print*,'single trait model' 
      nbcar=1
else if (xc(1).ne.'NBCHARACTER'.and.xc(1).ne.'MODEL') then
 print*,'error para file'
 stop           
  endif 
close(io_p)

 
if(nbcar.eq.1) then
!print*,'unicaractere avant lecture'
allocate(oante(1), ovar(1))
allocate(temps(NV))
 print*,'size time vector',NV
 temps=0
call lectureparamuni
print*,'parameter file OK'
else
call lectureparam
print*,'parameter file OK' 
endif
print*,'time vector', temps(:)
if (tmodel.eq.'SAD') then
print*,'number of parameters', ngammaante
!print*,'ngam',ngam	
if (NGAMMAANTE.ne.NGAM) then
STOP 'ERROR NB PARAMETERS'
endif	

 	READ(21,22,ERR=930,END=930) (GAMMAANTE(I),I=1,NGAM)
 
 22   FORMAT(6D14.7)
      CLOSE(21)
 
!debut ecriture du .gdg

  IFRM=-1
	BASENAME(LENBN:)='.gdg'
      OPEN(21,file=BASENAME,status='unknown',err=950)
      WRITE(21,21) NV*nbcar,NGAM,OPT,IFRM
 
      WRITE(21,22) (GAMMAANTE(I),I=1,NGAM)

if(nbcar.eq.1) then
 call makeinverseSAD
 call makeinversederivSAD
else
call makeinverseSADmulti
call makeinversederivSADmulti
endif

	close(21)

else if (tmodel.eq.'AR') then
!partie de Gilmour
print*,'it is not a SAD model, use another software'     


endif


!     This bit of code designed to make sure file is properly closed before returning

      OPEN(22,FILE=BASENAME,STATUS='OLD',ERR=950)
      READ(22,21,err=950) I
      CLOSE(22)

! ------------------------------------------
       STOP
910    STOP 'MYOWNGDG has no BASENAME argument on command line'
920    STOP 'MYOWNGDG could not open BASENAME.own'
930    STOP 'MYOWNGDG failed to read BASENAME.own'
940    STOP 'MYOWNGDG has invalid/inconsistent data in BASENAME.own'
950    STOP 'MYOWNGDG could not open BASENAME.gdg'
960    STOP 'MYOWNGDG has no BASENAME argument on command line'
 
 
 
       
contains
subroutine lectureparamuni 
 integer n, x(40),ef,stat,curr_eff,i,j,k,ii
 character xc(40)*40
integer, dimension(size(temps)) :: ttemps




open(io_p,file=BASEFILE)

 call readline(io_p,x,xc,n)
   if (xc(1) /= 'MODEL') then
      write(201,*)'MODEL not found'
      stop
  endif 
      
  call readline(io_p,x,xc,n) 
  tmodel=xc(1)
  

 if (tmodel.eq.'SAD') then
 
 call readline(io_p,x,xc,n)
   if (xc(1) /= 'ANTEDEP_ORDER') then
      write(201,*)'ANTEDEP_ORDER not found'
      stop
  endif 
  
      
  call readline(io_p,x,xc,n) 
  oanteuni=x(1)
  
  
allocate(tante(oanteuni))
  
  
   call readline(io_p,x,xc,n)
   if (xc(1) /= 'MODEL_ANTEDEP') then
       write(201,*)'MODEL_ANTEDEP not found'
       stop
   endif
   
   
	call readline(io_p,x,xc,n)
	tante(1:oanteuni)=x(1:oanteuni)
	
call readline(io_p,x,xc,n)
   if (xc(1) /= 'MODEL_INNOVATION_VARIANCE') then
       write(201,*)'MODEL_INNOVATION_VARIANCE not found'
       stop
   endif	
   call readline(io_p,x,xc,n)
	ovaruni=x(1)
	
	
   NGAMMAANTE=oanteuni+sum(tante(1:oanteuni))+1+ovaruni
   allocate(GAMMAANTE(oanteuni+sum(tante(1:oanteuni))+1+ovaruni),typepara(NGAM),typeorder(NGAM))
   GAMMAANTE=0.0d0
   typepara=0
   typeorder=0
do i=1,oanteuni
do j=0,tante(i)
typepara(sum(tante(1:i-1))+i+j)=i
typeorder(sum(tante(1:i-1))+i+j)=j
enddo
enddo   
do i=0,ovaruni
typeorder(sum(tante(1:oanteuni))+oanteuni+i+1)=i
enddo


call readline(io_p,x,xc,n)
   if (xc(1) /= 'DATA_FILE') then
       write(201,*)'DATA_FILE not found'
       stop
   endif	
   call readline(io_p,x,xc,n)
	file=xc(1)

call readline(io_p,x,xc,n)
   if (xc(1) /= 'COLUMN_FACTOR') then
       write(201,*)'COLUMN_FACTOR not found'
       stop
   endif	
   call readline(io_p,x,xc,n)
	col=x(1)

	
open(12,file=file)
jd=1

do
call readline(12,x,xc,n)
  if (n == -1) exit
 

 datatemps=x(col)
 
 call test
 
 if(jd.eq.NV+1) exit
enddo   
close (io_p)

! ttemps=QuickSort(temps) 
 !  temps=ttemps
 call tri(temps) 
    endif

end subroutine lectureparamuni

!###########################################################################################
!###########################################################################################
!###########################################################################################

subroutine lectureparam 
 integer n, x(40),ef,stat,curr_eff,i,j,k,ii,l,curseur,pointeurint,pointeurext,pointeurintbis
 character xc(40)*40




open(io_p,file=BASEFILE)
call readline(io_p,x,xc,n)

if (xc(1) /= 'NBCHARACTER') then
      write(201,*)'NBCHARACTER not found'
      stop
  endif 
     
  call readline(io_p,x,xc,n)
  
 
  nbcar=x(1)
print*,'NB CHARACTER :',nbcar
allocate(oante(nbcar),ovar(nbcar),oantecross(nbcar*(nbcar-1)),&
mattricross(nbcar*(nbcar-1),3),typecross(nbcar*(nbcar-1)),prescorrinit(nbcar*(nbcar-1)/2),&
mattriinit(nbcar*(nbcar-1)/2,3),prescorrinit2(nbcar*(nbcar-1)/2),&
debantecross(nbcar*(nbcar-1)),deltaantecross(nbcar*(nbcar-1)))
NV=NV/nbcar
allocate(temps(NV))
 temps=0
oante=0
oantecross=0
debantecross=0
deltaantecross=0
ovar=0

call readline(io_p,x,xc,n)


   if (xc(1) /= 'MODEL') then
      write(201,*)'MODEL not found'
      stop
  endif 
      
  call readline(io_p,x,xc,n) 
  tmodel=xc(1)
 
 if (tmodel.eq.'SAD') then
 

 call readline(io_p,x,xc,n)
   if (xc(1) /= 'ANTEDEP_ORDER') then
      write(201,*)'ANTEDEP_ORDER not found'
      stop
  endif 
  
      
  call readline(io_p,x,xc,n)
  if (n.ne.nbcar) then
	print*,'erreur nb antedep'
	stop
	endif 
  oante(1:nbcar)=x(1:nbcar)
 
 allocate(tante(sum(oante(1:nbcar*(nbcar+1)/2))))
 tante=0
  call readline(io_p,x,xc,n)
   if (xc(1) /= 'CROSS_ANTEDEP_ORDER') then
      write(201,*)'CROSS_ANTEDEP_ORDER not found'
      stop
  endif 
  
      
  call readline(io_p,x,xc,n)
  if (n.ne.nbcar*(nbcar-1)) then
	print*,'erreur nb cross_antedep'
	stop
	endif 
  oantecross(1:nbcar*(nbcar-1))=x(1:nbcar*(nbcar-1))
 
  valcross=0
  do i=1,nbcar*(nbcar-1)
  typecross(i)=valcross
  valcross=abs(valcross-1)
  enddo
 
 
 call readline(io_p,x,xc,n)
   if (xc(1) /= 'BEGIN_CROSS_ANTEDEP_ORDER') then
      write(201,*)'BEGIN_CROSS_ANTEDEP_ORDER not found'
      stop
  endif 
  
      
  call readline(io_p,x,xc,n)
  if (n.ne.nbcar*(nbcar-1)) then
	print*,'erreur nb begin_cross_antedep'
	stop
	endif 
  debantecross(1:nbcar*(nbcar-1))=x(1:nbcar*(nbcar-1))
 deltaantecross(1:nbcar*(nbcar-1))=oantecross(1:nbcar*(nbcar-1))-debantecross(1:nbcar*(nbcar-1))

!verif pas d erreur i.e. pas delta negatif
do i=1,nbcar*(nbcar-1)
  if (deltaantecross(i).lt.0) then
  print*,'error begin_cross_antedep>order_cross_antedep'
  stop
  endif
  enddo
 
 allocate(tantecross(sum(deltaantecross(1:nbcar*(nbcar-1)))))
 tantecross=0

   call readline(io_p,x,xc,n)
   if (xc(1) /= 'MODEL_ANTEDEP') then
       write(201,*)'MODEL_ANTEDEP not found'
       stop
   endif
   
   
	call readline(io_p,x,xc,n)
	if (n.ne.sum(oante(1:nbcar))) then
	print*,'erreur nb val antedep'
	stop
	endif
	tante(1:sum(oante(1:nbcar)))=x(1:sum(oante(1:nbcar)))
	
	call readline(io_p,x,xc,n)
   if (xc(1) /= 'MODEL_CROSS_ANTEDEP') then
       write(201,*)'MODEL_CROSS_ANTEDEP not found'
       stop
   endif
  
	call readline(io_p,x,xc,n)
	if (n.ne.sum(deltaantecross(1:nbcar*(nbcar-1))).and.&
	sum(deltaantecross(1:nbcar*(nbcar-1))).ne.0) then
	print*,'erreur nb val cross_antedep'
	stop
	endif
	tantecross(1:sum(deltaantecross(1:nbcar*(nbcar-1))))=x(1:sum(deltaantecross(1:nbcar*(nbcar-1))))

	
	
call readline(io_p,x,xc,n)
   if (xc(1) /= 'MODEL_INNOVATION_VARIANCE_COVARIANCE') then
       write(201,*)'MODEL_INNOVATION_VARIANCE_COVARIANCE not found'
       stop
   endif	
   call readline(io_p,x,xc,n)
   if (n.ne.nbcar) then
	print*,'erreur nb ovar'
	stop
	endif
	ovar(1:nbcar)=x(1:nbcar)
	

call readline(io_p,x,xc,n)
   if (xc(1) /= 'CORRINIT') then
       write(201,*)'CORRINIT not found'
       stop
   endif	
   call readline(io_p,x,xc,n)
  if (n.ne.nbcar*(nbcar-1)/2) then
	print*,'erreur nb corrinit'
	stop
	endif
   prescorrinit=x(1:nbcar*(nbcar-1)/2)
	
do l=1,nbcar*(nbcar-1)/2
if (prescorrinit(l).eq.1) then
prescorrinit2(l)=1
	if (debantecross((l-1)*2+1).eq.0.and.deltaantecross((l-1)*2+1).ne.0.or. &
	debantecross((l-1)*2+2).eq.0.and.deltaantecross((l-1)*2+2).ne.0) then
	print*,'error correlation between error term ==> surparametrized model, increase begin_cross_antedep_order'
	stop
	endif

else
prescorrinit2(l)=0
endif
enddo	
	
   NGAMMAANTE=sum(tante(1:sum(oante(1:nbcar)))+1)+sum(tantecross(1:sum(deltaantecross(1:nbcar*(nbcar-1))))+1) &
   +sum(ovar(1:nbcar)+1)+sum(prescorrinit2(1:(nbcar*(nbcar-1)/2)))
  
   allocate(GAMMAANTE(NGAMMAANTE),typepara(NGAMMAANTE),typeorder(NGAMMAANTE),&
   typeKr(NGAMMAANTE),typeKrant(NGAMMAANTE))
   GAMMAANTE=0.0d0
   typepara=0
   typeorder=0
   typeKr=0
   typeKrant=0
   pointeurint=0
   pointeurext=0
   pointeurintbis=0
   mattricross=0
!print*,'nbcar',nbcar,'oante',oante(:)
curseur=0
!print*,'boucle pour identif'
do l=1,nbcar*nbcar
if (l.le.nbcar) then
	do i=1,oante(l)
	do j=0,tante(sum(oante(1:(l-1)))+i)
	
	curseur=curseur+1
	
	typepara(curseur)=i
	typeorder(curseur)=j
	typeKr(curseur)=l
	typeKrant(curseur)=l
	enddo
	enddo

elseif (l.gt.nbcar) then
	if (typecross(l-nbcar).eq.0) then
		if (pointeurintbis.eq.pointeurext) then
		pointeurintbis=pointeurintbis+1
		pointeurext=0
		endif
  		 if (pointeurext.lt.pointeurintbis) then
   		 pointeurext=pointeurext+1
		 endif
	pointeurint=pointeurintbis+1
	endif
	if (typecross(l-nbcar).eq.0) then
	mattricross(l-nbcar,1)=(l-nbcar)
	mattricross(l-nbcar,2)=pointeurint
	mattricross(l-nbcar,3)=pointeurext
	else if (typecross(l-nbcar).eq.1) then
	mattricross(l-nbcar,1)=(l-nbcar)
	mattricross(l-nbcar,2)=pointeurext
	mattricross(l-nbcar,3)=pointeurint
	endif
	
	do i=1,deltaantecross(l-nbcar)
	do j=0,tantecross(sum(deltaantecross(1:(l-1-nbcar)))+i)
	
	curseur=curseur+1
	typepara(curseur)=i !+debantecross(l-1-nbcar)-1
	typeorder(curseur)=j
	
	if (typecross(l-nbcar).eq.0) then
	typeKr(curseur)=pointeurint
	typeKrant(curseur)=pointeurext
	else if (typecross(l-nbcar).eq.1) then
	typeKr(curseur)=pointeurext
	typeKrant(curseur)=pointeurint
	endif
	enddo
	enddo
endif

enddo   
!print*,'finboucle identif'
!print*,'typeKr',typeKr(:)
!print*,'typeKrant',typeKrant(:)
!print*,'typepara',typepara(:)
!print*,'typeorder',typeorder(:)
do l=1,nbcar
do i=0,ovar(l)
curseur=curseur+1

!print*,'l',l,'ovar(l)',ovar(l),'i',i,'curseur',curseur
!print*,'localisation', sum(tante(:)+1)+sum(ovar(1:l-1)+1)+i+1
typeorder(curseur)=i
typeKr(curseur)=l
typeKrant(curseur)=l
enddo
!print*,'curseur l',l,curseur
enddo



pointeurint=0
   pointeurext=0
   pointeurintbis=0

!mise a jour de corrinit

do l=1,nbcar*(nbcar-1)/2
if (pointeurintbis.eq.pointeurext) then
		pointeurintbis=pointeurintbis+1
		pointeurext=0
		endif
  		 if (pointeurext.lt.pointeurintbis) then
   		 pointeurext=pointeurext+1
		 endif
	pointeurint=pointeurintbis+1
mattriinit(l,1)=l
mattriinit(l,2)=pointeurint
mattriinit(l,3)=pointeurext
if (prescorrinit2(l).ne.0) then
curseur=curseur+1
typepara(curseur)=100
typeorder(curseur)=0
typeKr(curseur)=pointeurext
typeKrant(curseur)=pointeurint

endif
enddo

!print*,'typecross',typecross(:)
!print*,'typepara',typepara(:)
!print*,'typerorder',typeorder(:)
!print*,'debantecross',debantecross(:)
!print*,'finantecross',oantecross(:)
!print*,'nbtempscross',deltaantecross(:)
!print*,'typeKr',typeKr(:)
!print*,'typeKrant',typeKrant(:)
!print*,'mattricross',mattricross(1,:)
!print*,'mattricross',mattricross(2,:)

call readline(io_p,x,xc,n)
   if (xc(1) /= 'DATA_FILE') then
       write(201,*)'DATA_FILE not found'
       stop
   endif	
   call readline(io_p,x,xc,n)
	file=xc(1)
!print*,'fichier',file
call readline(io_p,x,xc,n)
   if (xc(1) /= 'COLUMN_FACTOR') then
       write(201,*)'COLUMN_FACTOR not found'
       stop
   endif	
   call readline(io_p,x,xc,n)
	col=x(1)	
	
open(12,file=file)

jd=1
do
call readline(12,x,xc,n)
!print*,'lecture temps data',x(:)
  if (n == -1) exit
 
 datatemps=x(col)
 !print*,'datatemps',datatemps
 !print*,'datatemps avant test',datatemps
if (NV.gt.1) then
 !print*,'test pour temps'
 call test
 else
! print*,'gestion temps OK'
 temps(1)=datatemps
 jd=jd+1
endif 
 !print*,'datatemps apres test',datatemps
 !print*,'temps avant test',temps(:)
 if(jd.eq.NV+1) exit
enddo   
close (io_p)

 !ttemps=QuickSort(temps) 
 !  temps=ttemps
 call tri(temps) 
    endif
!attention il manque un test pour savoir si c est la bonne colonne
!print*,'NV pour stop temps', NV

!print*,'le temps',temps(:)
end subroutine lectureparam

character(len=20) function str(k)
!   "Convert an integer to string."
    integer, intent(in) :: k
    write (str, *) k
    str = adjustl(str)
end function str

subroutine test
integer :: k

rtest=0
do k=1,jd

if (temps(k).ne.datatemps) then
rtest=rtest+1
endif
enddo

if (rtest.eq.jd) then
temps(jd)=datatemps

jd=jd+1
endif

end subroutine test

subroutine tri(t)
implicit none
integer, dimension (:), intent(inout) :: t
integer :: temp
integer :: i,j
do i=1,size(t)-1
	do j=i+1,size(t)
	if (t(i).gt.t(j)) then
	temp=t(i)
	t(i)=t(j)
	t(j)=temp
	endif
	enddo
enddo
end subroutine tri


recursive function quicksort(InList) result(OutList)
    
    INTEGER,DIMENSION(:) :: InList
    INTEGER,DIMENSION(size(InList)) :: OutList
    INTEGER,DIMENSION(size(InList)) :: SupList, OrderedSupList, InfList, OrderedInfList
    INTEGER :: pivot
    INTEGER :: i,j, InfListSize, SupListSize


if(size(InList,1) < 2) then
       OutList(1) = Inlist(1)
    else
       ! Le pivot sera le premier élément de la liste
       pivot = InList(1)
 
       ! On trie la liste
       InfListSize = 0
       SupListSize = 0
       do i = 2, size(InList,1)
          if(InList(i) < Pivot) then
             InfListSize = InfListSize + 1
             InfList(InfListSize) = InList(i)
          elseif(InList(i) >= Pivot) then
             SupListSize = SupListSize + 1
             SupList(SupListSize) = InList(i)
          end if
       enddo
 
       ! On recompose la liste
       if(InfListSize < 1) then
          OrderedSupList = QuickSort(SupList(1:SupListSize))
          OutList = (/ Pivot, (OrderedSupList(j), j=1,SupListSize) /)
       elseif(SupListSize < 1) then
          OrderedInfList = QuickSort(InfList(1:InfListSize))
          OutList = (/ (OrderedInfList(j), j=1,InfListSize), Pivot /)
       else
          OrderedInfList = QuickSort(InfList(1:InfListSize))
          OrderedSupList = QuickSort(SupList(1:SupListSize))
 
          OutList = (/ (OrderedInfList(j), j=1,InfListSize), Pivot, (OrderedSupList(j), j=1,SupListSize) /)
       endif
    end if

end function quicksort


function localisation(car1,car2) result(loc)
integer :: car1,car2,loc
integer :: pointeurint, pointeurext,pointeurintbis,pointexttemp,pointinttemp,zz
	
	if (car1.eq.car2) then
	loc=-1
	else if (car1.ne.car2) then	
	do zz=1,nbcar*(nbcar-1)
	pointeurint=mattricross(zz,2)
	pointeurext=mattricross(zz,3)	
		if (pointeurint.eq.car1.and.pointeurext.eq.car2) then
		loc=zz
		exit
	endif
	enddo
	endif
	

end function localisation

function locinit(car1,car2) result(loc)
integer :: car1,car2,loc
integer :: pointeurint, pointeurext,pointeurintbis,pointexttemp,pointinttemp,zz
	
	if (car1.eq.car2) then
	loc=-1
	else if (car1.ne.car2) then	
	do zz=1,nbcar*(nbcar-1)/2
	pointeurint=mattriinit(zz,2)
	pointeurext=mattriinit(zz,3)	
		if (pointeurint.eq.car1.and.pointeurext.eq.car2) then
		loc=zz
		exit
	endif
	enddo
	endif
	

end function locinit

subroutine makeinverseSAD
!subroutine qui calcul G-1=L'D-1L, selon formule G-1(ij)=sum(l=1 to NV)(sum(m=1 to NV)(L(mi)*L(lj)*D-1(ml)))
integer :: i,j,k,l,m
real :: val,theta,thetap,var
do i=1,NV
do j=max(i-oanteuni,1),i
val=0.0d0

	do l=1,NV
	theta=0.0d0
	if (l.eq.j) then
	theta=1.0d0
	elseif (l.gt.j.and.l.le.j+oanteuni) then
	do k=0,tante(l-j)
	theta=theta-gammaante(sum(tante(1:l-j-1))+l-j+k)*temps(l)**(k)
	enddo
	endif
	thetap=0.0d0
	if (l.eq.i) then
	thetap=1.0d0
	elseif (l.gt.i.and.l.le.i+oanteuni) then
	do k=0,tante(l-i)
	thetap=thetap-gammaante(sum(tante(1:l-i-1))+l-i+k)*temps(l)**(k)
	enddo
	endif
		var=1.0d0
		do k=0,ovaruni
		var=var*exp(gammaante(sum(tante(1:oanteuni))+oanteuni+k+1)*temps(l)**(k))
		enddo
		!print*,'variance au temps',temps(l),':',var
		!print*,'para utilise',gammaante(sum(tante(1:oante))+oante+1),gammaante(sum(tante(1:oante))+oante+2)
	val=val+(theta*thetap/var)
	enddo

		if (val.ne.0.0d0) then
		WRITE(21,'(I5,I5,G18.10)') i,j,val
		!print*,'i',i,'j',j,'val',val
		endif
enddo
enddo


end subroutine makeinverseSAD



subroutine makeinversederivSAD
!subroutine qui calcul G-1=L'D-1L, selon formule G-1(ij)=sum(l=1 to NV)(sum(m=1 to NV)(L(mi)*L(lj)*D-1(ml)))
integer :: i,j,k,l,m,ii
real :: val,theta,thetap,var

do ii=1,NGAM
 if (typepara(ii).gt.0) then

   do i=1,NV
   do j=max(i-oanteuni,1),i
   val=0.0d0
	!dL'/dpara*G-1L, theta qui change
	do l=1,NV
	
	theta=0.0d0
	if (l-j.eq.typepara(ii)) then
	theta=-temps(l)**typeorder(ii)
	endif
	thetap=0.0d0
	if (l.eq.i) then
	thetap=1.0d0
	elseif (l.gt.i.and.l.le.i+oanteuni) then
	do k=0,tante(l-i)
	thetap=thetap-gammaante(sum(tante(1:l-i-1))+l-i+k)*temps(l)**(k)
	enddo
	endif
		var=1.0d0
		do k=0,ovaruni
		var=var*exp(gammaante(sum(tante(1:oanteuni))+oanteuni+k+1)*temps(l)**(k))
		enddo
	val=val+(theta*thetap/var)
	enddo



	!L'*G-1dL/dpara, thetap qui change
	do l=1,NV
	
	theta=0.0d0
	if (l.eq.j) then
	theta=1.0d0
	elseif (l.gt.j.and.l.le.j+oanteuni) then
	do k=0,tante(l-j)
	theta=theta-gammaante(sum(tante(1:l-j-1))+l-j+k)*temps(l)**(k)
	enddo
	endif
	
	thetap=0.0d0
	if(l-i.eq.typepara(ii)) then
	thetap=-temps(l)**typeorder(ii)
	endif
	
	
	
		var=1.0d0
		do k=0,ovaruni
		var=var*exp(gammaante(sum(tante(1:oanteuni))+oanteuni+k+1)*temps(l)**(k))
		enddo
		
	val=val+(theta*thetap/var)
	enddo
	
   !if (val.ne.0.0d0) then
   WRITE(21,'(I5,I5,G18.10)') i,j,val
   !print*,'i',i,'j',j,'val',val
   !endif
   enddo
   enddo


 elseif (typepara(ii).eq.0) then !que L'dG-1/dparaL
 do i=1,NV
   do j=max(i-oanteuni,1),i
   val=0.0d0
 
	do l=1,NV
	theta=0.0d0
	if (l.eq.j) then
	theta=1.0d0
	elseif (l.gt.j.and.l.le.j+oanteuni) then
	do k=0,tante(l-j)
	theta=theta-gammaante(sum(tante(1:l-j-1))+l-j+k)*temps(l)**(k)
	enddo
	endif
	thetap=0.0d0
	if (l.eq.i) then
	thetap=1.0d0
	elseif (l.gt.i.and.l.le.i+oanteuni) then
	do k=0,tante(l-i)
	thetap=thetap-gammaante(sum(tante(1:l-i-1))+l-i+k)*temps(l)**(k)
	enddo
	endif
		var=1.0d0
		do k=0,ovaruni
		var=var*exp(-gammaante(sum(tante(1:oanteuni))+oanteuni+k+1)*temps(l)**(k))
		enddo
		if (typeorder(ii).eq.0) then
		var=var*(-1)
		else
		var=var*(-temps(l)**(typeorder(ii)))
		endif
	
	val=val+(theta*thetap*var)
	enddo

		!if (val.ne.0.0d0) then
		WRITE(21,'(I5,I5,G18.10)') i,j,val
		!print*,'i',i,'j',j,'val',val
		!endif
 enddo
 enddo
 endif
enddo
end subroutine makeinversederivSAD


subroutine makeinverseSADmulti
!subroutine qui calcul G-1=L'D-1L, selon formule G-1(ij)=sum(l=1 to NV)(sum(m=1 to NV)(L(mi)*L(lj)*D-1(ml)))
integer :: i,j,k,l,m,tt,ii,ll,cari,ccar,carj,carl,jj,dialj,diali,pos
real :: val,theta,thetap,var

!print*,'gammanante',gammaante(:)
!print*,'mattriinit',mattriinit(1,:)
!print*,'prescorinit',prescorrinit(1)
!print*,'ovar',ovar(:)
!print*,'temps',temps(:)
!print*,'deltaantecross',deltaantecross(:)
!print*,'oantecros',oantecross(:)
!print*,'tantecross',tantecross(:)
do i=1,nbcar*NV
cari=INT((i-1)/NV)+1
ii=i-(cari-1)*NV
do j=1,i
carj=INT((j-1)/NV)+1
jj=j-(carj-1)*NV
val=0.0d0

	
	
	do l=1,nbcar*NV
	ccar=INT((l-1)/NV)+1
	ll=l-(ccar-1)*NV
	!indicateur parametre
	dialj=ccar*(ccar-3)/2+carj+1
	diali=ccar*(ccar-3)/2+cari+1

	theta=0.0d0
	
	  if (l.eq.j) then
	  theta=1.0d0
	  elseif (l.gt.j.and.l.le.j+oante(ccar).and.carj.eq.ccar) then
		do k=0,tante(sum(oante(1:(ccar-1)))+l-j)
		theta=theta-gammaante(sum(tante(1:(sum(oante(1:ccar-1))+l-j-1))+1)+k+1)*temps(ll)**(k)
		enddo
	 
	  elseif (l.gt.j.and.carj.ne.ccar.and.ll.ge.jj+debantecross(localisation(ccar,carj)).and.&
	  ll.le.jj+oantecross(localisation(ccar,carj))-1 &
	  .and.oantecross(localisation(ccar,carj)).ne.0) then
	    	pos=sum(deltaantecross(1:(localisation(ccar,carj))-1))+ll-jj-debantecross(localisation(ccar,carj))
		do k=0,tantecross(pos+1)
		theta=theta-gammaante(sum(tante(1:(sum(oante(1:nbcar))))+1)+ &
		sum(tantecross(1:pos)+1)+k+1)*temps(ll)**(k)
				enddo				
	  ! print*,'cas 2','l',l,'j',j,'theta',theta,'localidation',localisation(ccar,carj) 
	elseif (l.lt.j.and.carj.ne.ccar.and.ll.ge.jj+debantecross(localisation(ccar,carj)).and.&
	ll.le.jj+oantecross(localisation(ccar,carj))-1 &
	.and.oantecross(localisation(ccar,carj)).ne.0) then
	    	pos=sum(deltaantecross(1:(localisation(ccar,carj))-1))+ll-jj-debantecross(localisation(ccar,carj))
		
		do k=0,tantecross(pos+1)
		theta=theta-gammaante(sum(tante(1:(sum(oante(1:nbcar))))+1)+ &
		sum(tantecross(1:pos)+1)+k+1)*temps(ll)**(k)
		
		enddo
		! print*,'cas 3','l',l,'j',j,'theta',theta,'localisation',localisation(ccar,carj) 
	 elseif (l.gt.j.and.carj.ne.ccar.and.ll.eq.1.and.jj.eq.1.and.prescorrinit(locinit(ccar,carj)).eq.1) then
	 	theta=theta-gammaante(ngammaante-sum(prescorrinit2(:))+sum(prescorrinit2(1:locinit(ccar,carj))))
		 
	   endif
	  

	  
	thetap=0.0d0
	  if (l.eq.i) then
	  thetap=1.0d0
	  elseif (l.gt.i.and.l.le.i+oante(ccar).and.cari.eq.ccar) then
		
		do k=0,tante(sum(oante(1:(ccar-1)))+l-i)
		thetap=thetap-gammaante(sum(tante(1:(sum(oante(1:ccar-1))+l-i-1))+1)+k+1)*temps(ll)**(k)
		enddo
	
	  elseif (l.gt.i.and.cari.ne.ccar.and.ll.ge.ii+debantecross(localisation(ccar,cari)).and.&
	  ll.le.ii+oantecross(localisation(ccar,cari))-1 &
	  .and.oantecross(localisation(ccar,cari)).ne.0) then
	    	
		pos=sum(deltaantecross(1:(localisation(ccar,cari))-1))+ll-ii-debantecross(localisation(ccar,cari))
		do k=0,tantecross(pos+1)
		thetap=thetap-gammaante(sum(tante(1:(sum(oante(1:nbcar))))+1)+ &
		sum(tantecross(1:pos)+1)+k+1)*temps(ll)**(k)
		
		enddo
	
	  elseif (l.lt.i.and.cari.ne.ccar.and.ll.ge.ii+debantecross(localisation(ccar,cari)).and.&
	  ll.le.ii+oantecross(localisation(ccar,cari))-1 &
	  .and.oantecross(localisation(ccar,cari)).ne.0) then
	    	pos=sum(deltaantecross(1:(localisation(ccar,cari))-1))+ll-ii-debantecross(localisation(ccar,cari))
		do k=0,tantecross(pos+1)
		thetap=thetap-gammaante(sum(tante(1:(sum(oante(1:nbcar))))+1)+ &
		sum(tantecross(1:pos)+1)+k+1)*temps(ll)**(k)
		enddo
	  elseif (l.gt.i.and.cari.ne.ccar.and.ll.eq.1.and.ii.eq.1.and.prescorrinit(locinit(ccar,cari)).eq.1) then
	 	thetap=thetap-gammaante(ngammaante-sum(prescorrinit2(:))+sum(prescorrinit2(1:locinit(ccar,cari))))	
	 
	  endif
	
	
	var=1.0d0
		do k=0,ovar(ccar)
		var=var*exp(gammaante(sum(tante(1:(sum(oante(1:nbcar))))+1)+ &
		sum(tantecross(1:sum(deltaantecross(1:nbcar*(nbcar-1))))+1)+ &
		sum(ovar(1:ccar-1)+1)+k+1)*temps(ll)**(k))
		enddo
		
	val=val+(theta*thetap/var)
		
	enddo
	
	
		if (val.ne.0.0d0) then
		WRITE(21,'(I5,I5,G18.10)') i,j,val
		!print*,'i',i,'j',j,'val',val
		endif
 
enddo
enddo

end subroutine makeinverseSADmulti











subroutine makeinversederivSADmulti
!subroutine qui calcul G-1=L'D-1L, selon formule G-1(ij)=sum(l=1 to NV)(sum(m=1 to NV)(L(mi)*L(lj)*D-1(ml)))
integer :: i,j,k,l,m,ii,nn,ll,cari,carj,ccar,jj,dialj,diali,pos
real :: val,theta,thetap,var


do nn=1,NGAM

 if (typepara(nn).gt.0) then

   do i=1,NV*nbcar
   cari=INT((i-1)/NV)+1
   ii=i-(cari-1)*NV
	do j=1,i
	carj=INT((j-1)/NV)+1
	jj=j-(carj-1)*NV
   
   
   val=0.0d0
	
	do l=1,NV*nbcar
	ccar=INT((l-1)/NV)+1
	ll=l-(ccar-1)*NV
	!indicateur parametre
	dialj=ccar*(ccar-3)/2+carj+1
	diali=ccar*(ccar-3)/2+cari+1
	
	theta=0.0d0
	!dL'/dpara*G-1L, theta qui change
	if (typepara(nn).ne.100.and.typeKr(nn).eq.ccar.and.typeKrant(nn).eq.carj) then
		!CAS 1: antedepintraKr
		if ((l.gt.j.and.l.le.j+oante(ccar).and.carj.eq.ccar.and. &
		     (ll-jj).eq.typepara(nn))&
		     .or. &
		!CAS 2: cross antedep sens 1==>2     
		     (l.gt.j.and.carj.ne.ccar.and.ll-jj.eq.typepara(nn)+debantecross(localisation(ccar,carj))-1.and.&
	 	     oantecross(localisation(ccar,carj)).ne.0) &
		     .or. &
		!CAS 3: cross antedep sens 2==>1
		     (l.lt.j.and.carj.ne.ccar.and.ll-jj.eq.typepara(nn)+debantecross(localisation(ccar,carj))-1.and.&
			oantecross(localisation(ccar,carj)).ne.0) &
			) then
		theta=-temps(ll)**typeorder(nn)
		if (l.gt.j.and.carj.ne.ccar.and.ll-jj.eq.typepara(nn)+debantecross(localisation(ccar,carj))-1.and.&
	 	     oantecross(localisation(ccar,carj)).ne.0) then
		   !  print*,'cas2 derive','nn',nn,'j',j,'l',l,'theta',theta
		elseif  (l.lt.j.and.carj.ne.ccar.and.ll-jj.eq.typepara(nn)+debantecross(localisation(ccar,carj))-1.and.&
			oantecross(localisation(ccar,carj)).ne.0) then     
		!print*,'cas3, derive','nn',nn,'j',j,'l',l,'theta',theta
		endif
		endif
	elseif (typepara(nn).eq.100) then	
		!CAS4: corrinit
		if (l.gt.j.and.carj.ne.ccar.and.ll.eq.1.and.jj.eq.1.and.prescorrinit(locinit(ccar,carj)).eq.1) then
		theta=-1
		
		endif
		
	endif	
	
	
	
	thetap=0.0d0
	  if (l.eq.i) then
	  thetap=1.0d0
	  elseif (l.gt.i.and.l.le.i+oante(ccar).and.cari.eq.ccar) then
		
		do k=0,tante(sum(oante(1:(ccar-1)))+l-i)
		thetap=thetap-gammaante(sum(tante(1:(sum(oante(1:ccar-1))+l-i-1))+1)+k+1)*temps(ll)**(k)
		enddo
	
	  elseif (l.gt.i.and.cari.ne.ccar.and.ll.ge.ii+debantecross(localisation(ccar,cari)).and.&
	  ll.le.ii+oantecross(localisation(ccar,cari))-1 &
	  .and.oantecross(localisation(ccar,cari)).ne.0) then
	    	
		pos=sum(deltaantecross(1:(localisation(ccar,cari))-1))+ll-ii-debantecross(localisation(ccar,cari))
		do k=0,tantecross(pos+1)
		thetap=thetap-gammaante(sum(tante(1:(sum(oante(1:nbcar))))+1)+ &
		sum(tantecross(1:pos)+1)+k+1)*temps(ll)**(k)
		
		enddo
	
	  elseif (l.lt.i.and.cari.ne.ccar.and.ll.ge.ii+debantecross(localisation(ccar,cari)).and.&
	  ll.le.ii+oantecross(localisation(ccar,cari))-1 &
	  .and.oantecross(localisation(ccar,cari)).ne.0) then
	    	pos=sum(deltaantecross(1:(localisation(ccar,cari))-1))+ll-ii-debantecross(localisation(ccar,cari))
		do k=0,tantecross(pos+1)
		thetap=thetap-gammaante(sum(tante(1:(sum(oante(1:nbcar))))+1)+ &
		sum(tantecross(1:pos)+1)+k+1)*temps(ll)**(k)
		enddo
	  elseif (l.gt.i.and.cari.ne.ccar.and.ll.eq.1.and.ii.eq.1.and.prescorrinit(locinit(ccar,cari)).eq.1) then
	 	thetap=thetap-gammaante(ngammaante-sum(prescorrinit2(:))+sum(prescorrinit2(1:locinit(ccar,cari))))	
	 
	  endif
	
	
	var=1.0d0
		do k=0,ovar(ccar)
		var=var*exp(gammaante(sum(tante(1:(sum(oante(1:nbcar))))+1)+ &
		sum(tantecross(1:sum(deltaantecross(1:nbcar*(nbcar-1))))+1)+ &
		sum(ovar(1:ccar-1)+1)+k+1)*temps(ll)**(k))
		enddo
	!if (i.eq.1.and.j.eq.1.and.nn.eq.1) then
	!print*,'nn',nn,'l',l,'i',i,'j',j,'theta deriv',theta,'thetap normal',thetap,'var',var
	!endif	
	val=val+(theta*thetap/var)
	
	
	
	enddo



	!L'*G-1dL/dpara, thetap qui change
	do l=1,NV*nbcar
	ccar=INT((l-1)/NV)+1
	ll=l-(ccar-1)*NV
	!indicateur parametre
	dialj=ccar*(ccar-3)/2+carj+1
	diali=ccar*(ccar-3)/2+cari+1
	
	theta=0.0d0
	
	  if (l.eq.j) then
	  theta=1.0d0
	  elseif (l.gt.j.and.l.le.j+oante(ccar).and.carj.eq.ccar) then
		do k=0,tante(sum(oante(1:(ccar-1)))+l-j)
		theta=theta-gammaante(sum(tante(1:(sum(oante(1:ccar-1))+l-j-1))+1)+k+1)*temps(ll)**(k)
		enddo
	 
	  elseif (l.gt.j.and.carj.ne.ccar.and.ll.ge.jj+debantecross(localisation(ccar,carj)).and.&
	  ll.le.jj+oantecross(localisation(ccar,carj))-1 &
	  .and.oantecross(localisation(ccar,carj)).ne.0) then
	    	pos=sum(deltaantecross(1:(localisation(ccar,carj))-1))+ll-jj-debantecross(localisation(ccar,carj))
		do k=0,tantecross(pos+1)
		theta=theta-gammaante(sum(tante(1:(sum(oante(1:nbcar))))+1)+ &
		sum(tantecross(1:pos)+1)+k+1)*temps(ll)**(k)
				enddo				
	  
	elseif (l.lt.j.and.carj.ne.ccar.and.ll.ge.jj+debantecross(localisation(ccar,carj)).and.&
	ll.le.jj+oantecross(localisation(ccar,carj))-1 &
	.and.oantecross(localisation(ccar,carj)).ne.0) then
	    	pos=sum(deltaantecross(1:(localisation(ccar,carj))-1))+ll-jj-debantecross(localisation(ccar,carj))
		
		do k=0,tantecross(pos+1)
		theta=theta-gammaante(sum(tante(1:(sum(oante(1:nbcar))))+1)+ &
		sum(tantecross(1:pos)+1)+k+1)*temps(ll)**(k)
		
		enddo
		  
	 elseif (l.gt.j.and.carj.ne.ccar.and.ll.eq.1.and.jj.eq.1.and.prescorrinit(locinit(ccar,carj)).eq.1) then
	 	theta=theta-gammaante(ngammaante-sum(prescorrinit2(:))+sum(prescorrinit2(1:locinit(ccar,carj))))
		 
	   endif
	
	thetap=0.0d0
	!dL'/dpara*G-1L, theta qui change
	if (typepara(nn).ne.100.and.typeKr(nn).eq.ccar.and.typeKrant(nn).eq.cari) then
		!CAS 1: antedepintraKr
		if ((l.gt.i.and.l.le.i+oante(ccar).and.cari.eq.ccar.and. &
		     (ll-ii).eq.typepara(nn))&
		     .or. &
		!CAS 2: cross antedep sens 1==>2     
		     (l.gt.i.and.cari.ne.ccar.and.ll-ii.eq.typepara(nn)+debantecross(localisation(ccar,cari))-1.and.&
	 	     oantecross(localisation(ccar,cari)).ne.0) &
		     .or. &
		!CAS 3: cross antedep sens 2==>1
		     (l.lt.i.and.cari.ne.ccar.and.ll-ii.eq.typepara(nn)+debantecross(localisation(ccar,cari))-1.and.&
			oantecross(localisation(ccar,cari)).ne.0) &
			) then
		thetap=-temps(ll)**typeorder(nn)
		
		
		endif
	elseif (typepara(nn).eq.100) then	
		!CAS4: corrinit
		if (l.gt.i.and.cari.ne.ccar.and.ll.eq.1.and.ii.eq.1.and.prescorrinit(locinit(ccar,cari)).eq.1) then
		thetap=-1
		endif
		
	endif	
	
	
	var=1.0d0
		do k=0,ovar(ccar)
		var=var*exp(gammaante(sum(tante(1:(sum(oante(1:nbcar))))+1)+ &
		sum(tantecross(1:sum(deltaantecross(1:nbcar*(nbcar-1))))+1)+ &
		sum(ovar(1:ccar-1)+1)+k+1)*temps(ll)**(k))
		enddo
		
	val=val+(theta*thetap/var)
	!if (i.eq.1.and.j.eq.1.and.nn.eq.1) then
	!print*,'nn',nn,'l',l,'i',i,'j',j,'theta ',theta,'thetap deriv',thetap,'var',var
	!endif
	enddo
	
	
	
  ! if (val.ne.0.0d0) then
   WRITE(21,'(I5,I5,G18.10)') i,j,val
   !print*,'i',i,'j',j,'val',val
  ! endif
   enddo
   enddo


 	elseif (typepara(nn).eq.0) then !que L'dG-1/dparaL
 	do i=1,nbcar*NV
	cari=INT((i-1)/NV)+1
	ii=i-(cari-1)*NV
	do j=1,i
	carj=INT((j-1)/NV)+1
	jj=j-(carj-1)*NV
	val=0.0d0

	
	
	do l=1,nbcar*NV
	ccar=INT((l-1)/NV)+1
	ll=l-(ccar-1)*NV
	!indicateur parametre
	dialj=ccar*(ccar-3)/2+carj+1
	diali=ccar*(ccar-3)/2+cari+1
	
	theta=0.0d0
	
	  if (l.eq.j) then
	  theta=1.0d0
	  elseif (l.gt.j.and.l.le.j+oante(ccar).and.carj.eq.ccar) then
		do k=0,tante(sum(oante(1:(ccar-1)))+l-j)
		theta=theta-gammaante(sum(tante(1:(sum(oante(1:ccar-1))+l-j-1))+1)+k+1)*temps(ll)**(k)
		enddo
	 
	  elseif (l.gt.j.and.carj.ne.ccar.and.ll.ge.jj+debantecross(localisation(ccar,carj)).and.&
	  ll.le.jj+oantecross(localisation(ccar,carj))-1 &
	  .and.oantecross(localisation(ccar,carj)).ne.0) then
	    	pos=sum(deltaantecross(1:(localisation(ccar,carj))-1))+ll-jj-debantecross(localisation(ccar,carj))
		do k=0,tantecross(pos+1)
		theta=theta-gammaante(sum(tante(1:(sum(oante(1:nbcar))))+1)+ &
		sum(tantecross(1:pos)+1)+k+1)*temps(ll)**(k)
				enddo				
	  
	elseif (l.lt.j.and.carj.ne.ccar.and.ll.ge.jj+debantecross(localisation(ccar,carj)).and.&
	ll.le.jj+oantecross(localisation(ccar,carj))-1 &
	.and.oantecross(localisation(ccar,carj)).ne.0) then
	    	pos=sum(deltaantecross(1:(localisation(ccar,carj))-1))+ll-jj-debantecross(localisation(ccar,carj))
		
		do k=0,tantecross(pos+1)
		theta=theta-gammaante(sum(tante(1:(sum(oante(1:nbcar))))+1)+ &
		sum(tantecross(1:pos)+1)+k+1)*temps(ll)**(k)
		
		enddo
		  
	 elseif (l.gt.j.and.carj.ne.ccar.and.ll.eq.1.and.jj.eq.1.and.prescorrinit(locinit(ccar,carj)).eq.1) then
	 	theta=theta-gammaante(ngammaante-sum(prescorrinit2(:))+sum(prescorrinit2(1:locinit(ccar,carj))))
		 
	   endif
	  

	  
	thetap=0.0d0
	  if (l.eq.i) then
	  thetap=1.0d0
	  elseif (l.gt.i.and.l.le.i+oante(ccar).and.cari.eq.ccar) then
		
		do k=0,tante(sum(oante(1:(ccar-1)))+l-i)
		thetap=thetap-gammaante(sum(tante(1:(sum(oante(1:ccar-1))+l-i-1))+1)+k+1)*temps(ll)**(k)
		enddo
	
	  elseif (l.gt.i.and.cari.ne.ccar.and.ll.ge.ii+debantecross(localisation(ccar,cari)).and.&
	  ll.le.ii+oantecross(localisation(ccar,cari))-1 &
	  .and.oantecross(localisation(ccar,cari)).ne.0) then
	    	
		pos=sum(deltaantecross(1:(localisation(ccar,cari))-1))+ll-ii-debantecross(localisation(ccar,cari))
		do k=0,tantecross(pos+1)
		thetap=thetap-gammaante(sum(tante(1:(sum(oante(1:nbcar))))+1)+ &
		sum(tantecross(1:pos)+1)+k+1)*temps(ll)**(k)
		
		enddo
	
	  elseif (l.lt.i.and.cari.ne.ccar.and.ll.ge.ii+debantecross(localisation(ccar,cari)).and.&
	  ll.le.ii+oantecross(localisation(ccar,cari))-1 &
	  .and.oantecross(localisation(ccar,cari)).ne.0) then
	    	pos=sum(deltaantecross(1:(localisation(ccar,cari))-1))+ll-ii-debantecross(localisation(ccar,cari))
		do k=0,tantecross(pos+1)
		thetap=thetap-gammaante(sum(tante(1:(sum(oante(1:nbcar))))+1)+ &
		sum(tantecross(1:pos)+1)+k+1)*temps(ll)**(k)
		enddo
	  elseif (l.gt.i.and.cari.ne.ccar.and.ll.eq.1.and.ii.eq.1.and.prescorrinit(locinit(ccar,cari)).eq.1) then
	 	thetap=thetap-gammaante(ngammaante-sum(prescorrinit2(:))+sum(prescorrinit2(1:locinit(ccar,cari))))	
	 
	  endif

	var=0.0d0	
		if (typeKr(nn).eq.ccar.and.typeKrant(nn).eq.ccar) then
		var=1.0d0
		do k=0,ovar(ccar)
		var=var*exp(-gammaante(sum(tante(1:(sum(oante(1:nbcar))))+1)+ &
		sum(tantecross(1:sum(deltaantecross(1:nbcar*(nbcar-1))))+1)+ &
		sum(ovar(1:ccar-1)+1)+k+1)*temps(ll)**(k))
		enddo
		
		if (typeorder(nn).eq.0) then
		var=var*(-1)
		else
		var=var*(-temps(ll)**(typeorder(nn)))
		endif
	
		endif
	!if (nn.eq.13.and.i.eq.j) print*,'nn',nn,typeKr(nn),typeKrant(nn),'l',l,'var',var,'theta',theta,'thetap',thetap	
	val=val+(theta*thetap*var)
	!if (nn.eq.13.and.i.eq.j) print*,'i',i,'j',j,'val',val
	enddo
	
	
		!if (val.ne.0.0d0) then
		WRITE(21,'(I5,I5,G18.10)') i,j,val
		!print*,'i',i,'j',j,'val',val
		!endif

enddo
enddo
endif

enddo
!PRINT*,'typeKr',typeKr(:)
!print*,'typeKrant',typeKrant(:)
!print*,'val',gammaante(:)
end subroutine makeinversederivSADmulti

subroutine readline(unit,x,xc,n,nohash)
! Reads a line from unit and decomposes it into numeric fields in x and
! alphanumeric fields in xc; the number of fields is stored in n.
! If optional variable nohash is absent, 
!  all characters after character # are ignored.

  integer unit,x(40),n,stat,i
  integer,optional::nohash
  character xc(40)*(*), a*200
 
  do
     read(unit,'(a)',iostat=stat)a
     if (stat /= 0) then
         n=-1
         return
     endif
     a=adjustl(a)               ! ignore leading spaces
     
     if (present(nohash)) then 
        exit
       else	!ignore characters after #
        i=scan(a,'#')
        if (i == 0) exit
        if (i>1) then
            a=a(1:i-1)
            exit
        endif
     endif  	
  enddo
  
  call nums2(a,n,x,xc)
  return
end subroutine


subroutine nums2(a,n,x,xc)
! separates array a into items delimited by blanks. character elements are
! put into optional character vector xc, decoded numeric values 
! into optional real vector x, and n contains the number of items. The 
! dimension of x and xc can be lower than n.
! A modification of nums() from f77 to f90
! Now accepts real numbers
! 2/23/2000

character (*)::a
character (*),optional::xc(:)
integer,optional::x(:)
integer::n,curr,first,last,lena,stat,i

curr=1;lena=len(a);n=0

do 
  ! search for first nonspace
  first=0
  do i=curr,lena
     if (a(i:i) /= ' ') then
        first=i
	exit
     endif
  enddo
  if (first == 0) exit
  
  
  ! search for first space
  curr=first+1
  last=0
  do i=curr,lena
      if (a(i:i) == ' ') then
        last=i
        exit
      endif
  enddo
  
  if (last == 0) last=lena
  
  n=n+1
  if (present(xc)) then
     if (size(xc) >= n) xc(n)=a(first:last)
  endif   
  if (present(x)) then
     if (size(x) >= n) then
        read(a(first:last),'(i12)',iostat=stat)x(n)
        if (stat /=0) x(n)=0
     endif
  endif   
	
  curr=last+1
enddo
end subroutine


end program MYOWNDGD


