cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c	Programme qui cree un maillage sol. Plusieurs couches peuvent
c	etre superposees. Le maillage est logarithmique pres des 
c	interfaces et regulier au milieu de chaque couche.
c	On cree un fichier qui contient le nombre total de noeuds, le 
c	nombre de noueds par couches et les dz entre chaque noeud.
c
c	Reference:	These A. Antonino (1992)
c
c	Auteur	: I.Braud, LTHE 
c
c	Routines: bibliotheque de routines subr_anto.f
c
c	Les significations de chaque variable peuvent etre trouvees
c	dans subr_anto.f ou les arguments des diverses routines sont
c	decrits.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

c	Parametres:	nnoeum= nombre maxi de noeuds dans le sol
c			ncoum= nombre maxi de couches de sol
c				de la vegetation
	
	parameter(nnoeum=200,ncoum=5)

	character*80 filename,filename1

c	Caracteristiques du sol
	dimension epaiscouch(ncoum),epaish(ncoum),epaisb(ncoum),
     &		dzhmin(ncoum),z(nnoeum),
     &		nnmil(ncoum),nntot(ncoum),dzbmin(ncoum),dz(nnoeum)

	data iti,ito /5,6/
	data idon,ilist,ires /10,11,12/

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c	Lecture des donnees
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

	open(idon,file='maillage.par',access='sequential',status=
     &		'unknown',form='formatted')
c	call ouvrir(iti,ito,idon,0,0,'r',filename,'Donnees        ')

c	Lecture des donnees relatives aux caracteristiques  du maillage

c	Nombre de couches
	read(idon,*) ncouch
	n=1

	do 100 i=1,ncouch

c	Lecture des caracteristiques du maillage pour chaque couche

	read(idon,*) epaiscouch(i),epaish(i),epaisb(i),dzhmin(i),
     &			dzbmin(i),nnmil(i)

c	Calcul du maillage pour la couche. n sera le nombre total de noeuds
c	dz contiendra les epaisseurs entre chaque noeud.

	call maillagecouche(nnmil(i),dzhmin(i),dzbmin(i),epaiscouch(i),
     &			epaish(i),epaisb(i),nntot(i),dz(n))
	n=n+nntot(i)
100	continue
	nntot(ncouch)=nntot(ncouch)+1

c	Test si n>nnoeum

	nn=nnoeum
	if(n.gt.nn) then
		write(iti,1000)
1000		format('nombre de noeuds trop grand.Arret.')
		stop
	end if

c	Ecriture du fichier du maillage

	call ouvrir(iti,ito,ilist,0,0,'w',filename1,'Maillage      ')
	write(ilist,*) ncouch
	write(ilist,*) (nntot(i),i=1,ncouch)
	write(ilist,*) n
	write(ilist,*) (dz(i),i=1,n-1)
	close(ilist)

c	Calcul de la somme des dz et des profondeurs

	som=0.
	z(1)=0. 
	do 200 i=2,n
	som=som+dz(i-1)
	z(i)=z(i-1)+dz(i-1)
200	continue

	write(ito,*) som
	write(ito,*) (z(i),i=1,n)

	stop 900
	end

	subroutine maillagecouche(nnmil, dzhmin,dzbmin,epaiscouch,
     &				epaish, epaisb,nntot,dz)

	dimension dz(1)
	
	gamh=1.25
	gamb=0.8

c	Calcul du nombre de noeuds sur les interfaces haut et bas

	xnnh=alog(1.-epaish*(1.-gamh)/dzhmin)/alog(gamh)
	nnh=int(xnnh)
	xnnh=real(nnh)
	xnnb=alog(1.-epaisb*(1.-gamh)/dzbmin)/alog(gamh)
	nnb=int(xnnb)
	xnnb=real(nnb)

c	On recalcule l'epaisseur du premier compartiment des interfaces 
c	haut et bas. Elles ne valent pas tout a fait dzhmin et dzbmax car 
c	on arrondit nnh et nnb aux entiers les plus proches.

	dzh1=epaish*(1.-gamh)/(1.-gamh**xnnh)
	dzb1=epaisb*(1.-gamh)/(1.-gamh**xnnb)

c	Calcul des epaisseurs de l'interface haut.

	do 100 i=1, nnh
	xi=real(i)
	dz(i)=dzh1*gamh**(xi-1.)
100	continue

c	Calcul des epaisseurs intermediaires (maillage regulier)

	epaismil=epaiscouch-epaish-epaisb
	nfin=nnh+nnmil
	xnnmil=real(nnmil)
	do 110 i=nnh+1, nfin
	dz(i)=epaismil/xnnmil
110	continue

c	Calcul des epaisseurs de l'interface inferieure

	ndeb=nfin+1
	nfin=nfin+nnb
	xi=0
	do 120 i=ndeb,nfin
	xi=xi+1.
	dz(i)=dzb1*gamh**(xnnb-xi)
120	continue

	nntot=nnh+nnmil+nnb
	return
	end