c234567
      program princip

c calcul de la solution de l equation du ptentiel
c par elements finis quadratiques Q9
c conditions aux limites :
c de Dirichlet sur Gamma_1  u = g
c de Neumann sur Gamma_n    du/dn = 0

c reservation memoire

c234567
      implicit real*8(c)
      implicit real*8(a-b,d-h,o-z)
      parameter(nmem=1700000)
c#include "param"
      parameter(ntanak=80,nhmai=260,nvmai=60,along=13.d0,
     & h1=1.d0)

      dimension atab(nmem)
      dimension itab(nmem)

c calcul des dimensions des tableaux

c nombre de noeuds
      nbsom = (nhmai+1)*(nvmai+1)
c nombre d'éléments
      nbel = nhmai/2*nvmai/2      

      write(*,*)'nombre d elements',nbel
      write(*,*) 'nombre de noeuds',nbsom

c nombre de ddl (1 par noeud)

      ndl=nbsom

c calcul des pointeurs entiers
      write(*,*) 'calcul des pointeurs'

c connectivite elts->ddl (9*nbel)
      iptie=1
c tableau des noeuds imposes (nbsom)
      ipticlim=iptie+9*nbel
c hauteurs des colonnes (ndl)
      iptkhj=ipticlim+nbsom
c tableau des debuts de colonnes (ndl+1)
      iptkld=iptkhj+ndl

c calcul des pointeurs reels et complexes
c (un complexe = 2 reels)

c coordonnees des noeuds (2*nbsom)
      iptx=1
      ipty=iptx+nbsom
c valeurs des noeuds imposes (ndl)
      iptclim=ipty+2*nbsom
c diagonale de K (ndl)
      iptvkgd=iptclim+nbsom
c triangle sup. de K (lvkgs : inconnu !)
      iptvkgs=iptvkgd+ndl

c lecture du maillage

      write(*,*) 'lecture du maillage'

      call lecm(nbel,nbsom,nbclim,itab(iptie),
     $  atab(iptx),atab(ipty),
     $  itab(ipticlim),atab(iptclim),xom)

c initialisations de :
c hauteurs de colonnes
c et tables des debuts de colonnes

      write(*,*) 'initialisations'


      call init(nbel,nbsom,ndl,itab(iptie)
     &,itab(iptkhj),
     & itab(iptkld),lvkgs)

c test de memoire
      write(*,*) 'longueur de vkgs ',lvkgs
      write(*,*) 'reste',nmem-lvkgs-iptvkgs,' de libre'
      if (lvkgs+iptvkgs.gt.nmem) then
      write(*,*) 'pas assez de memoire'
      write(*,*) 'modifier nmem'
      stop
      endif

c vecteur second membre (ndl)
c234567
      iptforce=iptvkgs+lvkgs

c assemblage de la matrice de rigidite
c assemblage du vecteur second membre
c imposition des c-l par penalisation

      write(*,*) 'assemblage'

      call assemb(nbclim,nbel,ndl,itab(iptie),itab(iptkld)
     &, atab(iptvkgd),atab(iptvkgs)
     &, lvkgs,atab(iptforce),atab(iptx),atab(ipty)
     &,itab(ipticlim),atab(iptclim),xom)

c pointeur de la solution (ndl)

      iptsolu=iptforce+ndl

c factorisation de la matrice de rigidite
      write(*,*) 'factorisation'

      call sol(atab(iptvkgs),atab(iptvkgd)
     &,atab(iptvkgs),atab(iptforce),
     &itab(iptkld),atab(iptsolu),ndl,6,1,0,0,energ,ier)

c calcul de la solution
      write(*,*) 'resolution'

      call sol(atab(iptvkgs),atab(iptvkgd)
     &,atab(iptvkgs),atab(iptforce),
     & itab(iptkld),atab(iptsolu),ndl,6,0,1,0,energ,ier)

c sauvegarde de la solution et calcul de l'erreur

c pointeur de la vrai sol (ndl)

      iptvraisol=iptsolu+ndl
    
      write(*,*) 'calcul du gradient'

      call sauv(nbel,nbsom,itab(iptie),
     & atab(iptx),atab(ipty),atab(iptsolu)
     &,atab(iptvraisol))

      lastptr=iptvraisol+ndl
      write(*,*) 'reels : mem total',nmem
     &,'  utilise',lastptr
      lastpti=iptkld+ndl+1
      write(*,*) 'entiers : mem total',nmem
     &,'  utilise',lastpti

      end

     
c234567
      subroutine assemb(nbclim,nel,ndl,ie,kld,
     & vkgd,vkgs,lvkgs,force,xx,yy,
     & iclim,clim,xom)

c assemblage du second membre et 
c de la matrice de rigidite

      implicit real*8(c)
      implicit real*8(a-b,d-h,o-z)

      parameter(npg=9,ndle=9)

      integer ie(ndle,nel),kld(ndl+1)
      dimension force(ndl),vkgs(lvkgs),vkgd(ndl)
     &,xx(ndl),yy(ndl)
      dimension xpg(npg),xref(ndle),yref(ndle)
     & ,vkelem(ndle,ndle),ypg(npg),hpg(npg),
     & dmat(3,3),xp(2),vni(ndle,3,npg),ajac(2,2),
     & ajinv(2,2),bmat(3,ndle),dbmat(3,ndle),btdb(ndle,ndle)
     &,iclim(ndl),clim(ndl)
      
c element de reference

      xref(1)=0
      xref(2)=.5d0
      xref(3)=1.d0
      xref(4)=1.d0
      xref(5)=1.d0
      xref(6)=.5d0
      xref(7)=0
      xref(8)=0
      xref(9)=.5d0

      yref(1)=0
      yref(2)=0
      yref(3)=0
      yref(4)=.5d0
      yref(5)=1.d0
      yref(6)=1.d0
      yref(7)=1.d0
      yref(8)=.5d0
      yref(9)=.5d0


c stockage des points de gauss et des poids

      xpi=4.d0*datan(1.d0)

      ab=dsqrt(3.d0/5.d0)

      xpg(1)=-ab
      xpg(2)=0.d0
      xpg(3)=ab
      xpg(4)=-ab
      xpg(5)=0.d0
      xpg(6)=ab
      xpg(7)=-ab
      xpg(8)=0.d0
      xpg(9)=ab

      ypg(1)=-ab
      ypg(2)=-ab
      ypg(3)=-ab
      ypg(4)=0
      ypg(5)=0
      ypg(6)=0
      ypg(7)=ab
      ypg(8)=ab
      ypg(9)=ab

      do i=1,9
         xpg(i)=(xpg(i)+1.d0)*.5d0
         ypg(i)=(ypg(i)+1.d0)*.5d0
      end do
         

      ap=5.d0/9.d0
      bp=8.d0/9.d0
      cp=5.d0/9.d0

      ap=.5d0*ap
      bp=.5d0*bp
      cp=.5d0*cp

      hpg(1)=ap*ap
      hpg(2)=bp*ap
      hpg(3)=cp*ap
      hpg(4)=ap*bp
      hpg(5)=bp*bp
      hpg(6)=cp*bp
      hpg(7)=ap*cp
      hpg(8)=bp*cp
      hpg(9)=cp*cp

c calcul des fonctions de base
c et des derivees aux 
c points de gauss
c234567
      un=1.d0
      de=2.d0
      qa=4.d0
  
      do ipg=1,npg

      x=xpg(ipg)
      y=ypg(ipg)

c fonction au point de gauss

      vni(1,1,ipg)= 0.4D1*(x-0.5D0)*(x-1.D0)*(y-0.5D0)*(y-1.D0)
      vni(2,1,ipg)= -0.8D1*x*(x-1.D0)*(y-0.5D0)*(y-1.D0)
      vni(3,1,ipg)= 0.4D1*x*(x-0.5D0)*(y-0.5D0)*(y-1.D0)
      vni(4,1,ipg)= -0.8D1*x*(x-0.5D0)*y*(y-1.D0)
      vni(5,1,ipg)= 0.4D1*x*(x-0.5D0)*y*(y-0.5D0)
      vni(6,1,ipg)= -0.8D1*x*(x-1.D0)*y*(y-0.5D0)
      vni(7,1,ipg)= 0.4D1*(x-0.5D0)*(x-1.D0)*y*(y-0.5D0)
      vni(8,1,ipg)= -0.8D1*(x-0.5D0)*(x-1.D0)*y*(y-1.D0)
      vni(9,1,ipg)= 0.16D2*x*(x-1.D0)*y*(y-1.D0)
     
c derivee /xi de la fonction au point de gauss

      vni(1,2,ipg)= 0.4D1*(x-1.D0)*(y-0.5D0)*(y-1.D0)+0.4D1*(x-0.5D0)*(y
     #-0.5D0)*(y-1.D0)
      vni(2,2,ipg)= -0.8D1*(x-1.D0)*(y-0.5D0)*(y-1.D0)-0.8D1*x*(y-0.5D0)
     #*(y-1.D0)
      vni(3,2,ipg)= 0.4D1*(x-0.5D0)*(y-0.5D0)*(y-1.D0)+0.4D1*x*(y-0.5D0)
     #*(y-1.D0)
      vni(4,2,ipg)= -0.8D1*(x-0.5D0)*y*(y-1.D0)-0.8D1*x*y*(y-1.D0)
      vni(5,2,ipg)= 0.4D1*(x-0.5D0)*y*(y-0.5D0)+0.4D1*x*y*(y-0.5D0)
      vni(6,2,ipg)= -0.8D1*(x-1.D0)*y*(y-0.5D0)-0.8D1*x*y*(y-0.5D0)
      vni(7,2,ipg)= 0.4D1*(x-1.D0)*y*(y-0.5D0)+0.4D1*(x-0.5D0)
     # *y*(y-0.5D0)
      vni(8,2,ipg)= -0.8D1*(x-1.D0)*y*(y-1.D0)-0.8D1*(x-0.5D0)
     #*y*(y-1.D0)
      vni(9,2,ipg)= 0.16D2*(x-1.D0)*y*(y-1.D0)+0.16D2*x*y*(y-1.D0)

c derivee /eta de la fonction au point de gauss


      vni(1,3,ipg)= 0.4D1*(x-0.5D0)*(x-1.D0)*(y-1.D0)+0.4D1*(x-0.5D0)*(x
     #-1.D0)*(y-0.5D0)
      vni(2,3,ipg)=-0.8D1*x*(x-1.D0)*(y-1.D0)-0.8D1*x*(x-1.D0)*(y-0.5D0)
      vni(3,3,ipg)=0.4D1*x*(x-0.5D0)*(y-1.D0)+0.4D1*x*(x-0.5D0)
     #*(y-0.5D0)
      vni(4,3,ipg)= -0.8D1*x*(x-0.5D0)*(y-1.D0)-0.8D1*x*(x-0.5D0)*y
      vni(5,3,ipg)= 0.4D1*x*(x-0.5D0)*(y-0.5D0)+0.4D1*x*(x-0.5D0)*y
      vni(6,3,ipg)= -0.8D1*x*(x-1.D0)*(y-0.5D0)-0.8D1*x*(x-1.D0)*y
      vni(7,3,ipg)= 0.4D1*(x-0.5D0)*(x-1.D0)*(y-0.5D0)+0.4D1*(x-0.5D0)*(
     #x-1.D0)*y
      vni(8,3,ipg)= -0.8D1*(x-0.5D0)*(x-1.D0)*(y-1.D0)-0.8D1*(x-0.5D0)*(
     #x-1.D0)*y
      vni(9,3,ipg)= 0.16D2*x*(x-1.D0)*(y-1.D0)+0.16D2*x*(x-1.D0)*y

      enddo

c calcul de la matrice D

      do ii=1,3
      do jj=1,3
      dmat(ii,jj)=0.d0
      enddo
      enddo

      dmat(1,1)=xom**2
      dmat(2,2)=-1.d0
      dmat(3,3)=-1.d0
     
c boucle sur les elements
c assemblage de la matrice de rigidite

      do ine=1,nel

      do ii=1,ndle
      do jj=1,ndle
      vkelem(ii,jj)=0.d0
      enddo
      enddo

c calcul de la matrice de rigidite elementaire

      do ipg=1,npg

c calcul de la matrice jacobienne

      do ii=1,2
      do jj=1,2
      ajac(ii,jj)=0.d0
      enddo
      enddo

      do k=1,ndle

      xp(1)=xx(ie(k,ine))
      xp(2)=yy(ie(k,ine))
c$$$      xp(1)=yref(k)*yref(k)
c$$$      xp(2)=-xref(k)

      do ii=1,2
      do jj=1,2

      ajac(ii,jj)=ajac(ii,jj)+xp(ii)*vni(k,jj+1,ipg)

      enddo
      enddo

      enddo

c$$$      write(*,*) ajac(1,1),' ',ajac(1,2)
c$$$      write(*,*) ajac(2,1),' ',ajac(2,2)
c$$$      stop
c calcul du jacobien

      det=ajac(1,1)*ajac(2,2)-ajac(1,2)*ajac(2,1)

      if (det.le.0.d0) then
      write(*,*) 'det negatif sur l element',ine,det
      endif

c calcul de la jacobienne inverse

      ajinv(1,1)=ajac(2,2)/det
      ajinv(2,2)=ajac(1,1)/det
      ajinv(1,2)=-ajac(1,2)/det
      ajinv(2,1)=-ajac(2,1)/det

c calcul de la matrice B
      
      do jj=1,ndle
      bmat(1,jj)=vni(jj,1,ipg)
      bmat(2,jj)=0.d0
      bmat(3,jj)=0.d0
      enddo

      do j=1,ndle
      do ii=2,3
      do kk=1,2
      bmat(ii,j)=bmat(ii,j)+ajinv(kk,ii-1)*vni(j,kk+1,ipg)
      enddo
      enddo
      enddo

c calcul de D*B

      do ii=1,3
      do jj=1,ndle
      dbmat(ii,jj)=0.d0
      do kk=1,3
      dbmat(ii,jj)=dbmat(ii,jj)+dmat(ii,kk)*bmat(kk,jj)
      enddo
      enddo
      enddo

c calcul de Bt*D*B

      do ii=1,ndle
      do jj=1,ndle
      btdb(ii,jj)=0.d0
      do kk=1,3
      btdb(ii,jj)=btdb(ii,jj)+bmat(kk,ii)*dbmat(kk,jj)
      enddo
      enddo
      enddo

      do ii=1,ndle
      do jj=1,ndle
      vkelem(ii,jj)=vkelem(ii,jj)+btdb(ii,jj)*det*hpg(ipg)
      enddo
      enddo
     
c fin de la boucle sur les points de gauss

      enddo

c distribution dans la matrice ligne de ciel

      do ii=1,ndle
      do jj=1,ndle

      i=ie(ii,ine)
      j=ie(jj,ine)

      if (i.eq.j) then
      vkgd(i)=vkgd(i)+vkelem(ii,ii)
      endif

      if (i.lt.j) then
      l=kld(j+1)-j+i
      vkgs(l)=vkgs(l)+vkelem(ii,jj)
      endif

      enddo
      enddo

c fin de la boucle sur les elements

      enddo

c boucle sur les elements
c assemblage du second membre

      do ine=1,nel
      do ii=1,ndle
    
      i=ie(ii,ine)
  
      do ipg=1,npg
      force(i)=0.d0
      enddo
   
      enddo
      enddo

c imposition des cond. aux limites
c par penalisation

      do ii=1,nbclim
      i=iclim(ii)
      vkgd(i)=1.d20
      force(i)=clim(ii)*1.d20
      enddo

      return
      end

c234567
      subroutine  init(nel,nbsom,ndl,ie
     &,khj, kld,lvkgs)

c initialisations de :
c hauteurs de colonnes
c debuts de colonnes

c234567
      implicit complex*16(c)
      implicit real*8(a-b,d-h,o-z)

      integer ie(9,nel),khj(ndl),kld(ndl+1)

c calcul des hauteurs de colonnes
      do i=1,ndl
      khj(i)=0
      enddo
 
      do ine=1,nel

      do ii=1,9
      do jj=1,9

      i=ie(ii,ine)
      j=ie(jj,ine)

      if(i.lt.j) then
       if(j-i.gt.khj(j))then
       khj(j)=j-i
       endif
      endif

      enddo
      enddo

      enddo

      do i=1,ndl
c      write(*,*) 'hauteur',i,' =',khj(i)
      enddo

c calcul de kld

      kld(1)=1
      do i=2,ndl+1
      kld(i)=kld(i-1)+khj(i-1)
      enddo

c calcul dela longueur de vkgs

      lvkgs=kld(ndl+1)-1

      do i=1,ndl+1
c      write(*,*) 'debut col.',i,' =',kld(i)
      enddo
     
      return
      end
      subroutine lecm(nbel,nbsom,nbclim,ie,xx,yy,iclim,clim,xom)

c lecture du maillage et de la connectivite

      implicit real*8(c)
      implicit real*8(a-b,d-h,o-z)

c dimension : fichier tanaka, nbre de points horizontaux et verticaux
c#include "param"
      parameter(ntanak=80,nhmai=260,nvmai=60,along=13.d0,
     & h1=1.d0)

      dimension xx(nbsom),yy(nbsom),clim(nbsom)
      integer kinterm(16)
      dimension ie(9,nbel),iclim(nbsom)

      dimension xtanak(ntanak),ytanak(ntanak),faitanak(ntanak)
      dimension xmaiup(0:nhmai),ymaiup(0:nhmai),faiup(0:nhmai)

c ouverture du fichier des valeurs données par
c la méthode de tanaka

      open(1,file='tanaka.dat',status='old')
      do i=1,ntanak
         read(1,*) j,xtanak(i),ytanak(i),faitanak(i)
      end do
      close(1)

c nombre de noeuds
      nbsom = (nhmai+1)*(nvmai+1)
c nombre d'éléments
      nbel = nhmai/2*nvmai/2

c      write(*,*) 'nbre de sommets=',nbsom
c      write(*,*) 'nbre d elements=',nbel

c longueur du domaine de calcul (adim)
c      along = 13.d0
c hauteur d'eau à l'infini (adim)
c      h1=1.d0
c pas du maillage en x
      dx = along/nhmai

      do i = 0,nhmai
         xmaiup(i)=i*dx
c recherche de l'intervalle ou se trouve le point
         j=0
 555     j=j+1

         if (xmaiup(i).ge.xtanak(j)) then
            goto 555
         endif
         
c interpolation linéaire entre j-1 et j
         t=(xmaiup(i)-xtanak(j-1))/(xtanak(j)-xtanak(j-1))
         ymaiup(i)=t*ytanak(j)+(1-t)*ytanak(j-1)
         faiup(i)=t*faitanak(j)+(1-t)*faitanak(j-1)

      end do

c$$$      do i = 0,nhmai
c$$$         write(*,*) i,xmaiup(i),ymaiup(i),faiup(i)
c$$$      enddo
c$$$      stop
c calcul  des coor des sommets

      do i =0,nhmai
         do j =0,nvmai
            dy=(ymaiup(i)+h1)/nvmai
            ii=j+(nvmai+1)*i+1
            xx(ii)=xmaiup(i)
            yy(ii)=-h1+j*dy
         end do
      end do


c$$$      open(2,file='g1',status='unknown')
c$$$      do i=1,nbsom
c$$$      write(2,*) xx(i),yy(i)
c$$$      enddo
c$$$      close(2)

c calcul de la connectivite

      do i=0,nhmai-1,2
         do j=0,nvmai-1,2
          
            iel = j/2+(nvmai/2)*i/2+1
            ii1=(j+0)+(nvmai+1)*(i+0)+1
            ii2=(j+0)+(nvmai+1)*(i+1)+1
            ii3=(j+0)+(nvmai+1)*(i+2)+1
            ii4=(j+1)+(nvmai+1)*(i+2)+1
            ii5=(j+2)+(nvmai+1)*(i+2)+1
            ii6=(j+2)+(nvmai+1)*(i+1)+1
            ii7=(j+2)+(nvmai+1)*(i+0)+1
            ii8=(j+1)+(nvmai+1)*(i+0)+1
            ii9=(j+1)+(nvmai+1)*(i+1)+1

            ie(1,iel)=ii1
            ie(2,iel)=ii2
            ie(3,iel)=ii3
            ie(4,iel)=ii4
            ie(5,iel)=ii5
            ie(6,iel)=ii6
            ie(7,iel)=ii7
            ie(8,iel)=ii8
            ie(9,iel)=ii9
            
         end do
      end do
c$$$      open(3,file='g2',status='unknown')
c$$$      do i=1,nbel
c$$$      write(3,*) xx(ie(1,i)),yy(ie(1,i))
c$$$      write(3,*) xx(ie(3,i)),yy(ie(3,i))
c$$$      write(3,*) xx(ie(5,i)),yy(ie(5,i))
c$$$      write(3,*) xx(ie(7,i)),yy(ie(7,i))
c$$$      enddo
c$$$      close(3)
c$$$      stop

c$$$      do i=1,nbel
c$$$         write(*,*) i,'-->',ie(1,i),ie(2,i),ie(3,i),ie(4,i)
c$$$     &        ,ie(5,i),ie(6,i),ie(7,i),ie(8,i),ie(9,i)
c$$$      enddo
c$$$      stop


c calcul des conditions aux limites

c bord gauche

      nbclim=0
      
      do j=0,nvmai-1
         nbclim=nbclim+1
         i=0
         ii=j+(nvmai+1)*i+1
         clim(nbclim)=0
c         clim(nbclim)=1
         iclim(nbclim)=ii
      end do

c partie sup

      do i=0,nhmai
         nbclim=nbclim+1
         j=nvmai
         ii=j+(nvmai+1)*i+1
         clim(nbclim)=faiup(i)
c         clim(nbclim)=1
         iclim(nbclim)=ii
      end do

c affichage

      xom=0
      write(*,*) 'nbre de c_l=',nbclim

c$$$      do i=1,nbclim
c$$$         xa=xx(iclim(i))
c$$$         ya=yy(iclim(i))
c$$$         write(*,*) i,' no sommet=',iclim(i),'  val=',clim(i)
c$$$      enddo

      return
      end

      subroutine sauv(nbel,nbsom,ie,xx,yy,sol,vraisol)

c lecture du maillage et de la connectivite

      implicit real*8(c)
      implicit real*8(a-b,d-h,o-z)

c#include "param"
      parameter(ntanak=80,nhmai=260,nvmai=60,along=13.d0,
     & h1=1d0)

      dimension xx(nbsom),yy(nbsom)
      dimension ie(9,nbel),sol(nbsom),
     & vraisol(nbsom),vni(9,3),ajac(2,2),
     & ajinv(2,2),vnivrai(9,3)
      dimension xref(9),yref(9)

c ouverture des fichiers

c$$$      open(1,file='MOSMEF.SOL',status='unknown')
c$$$
c$$$      write(1,*)
c$$$      write(1,*)
c$$$      write(1,'(6i5)') 1,1,1,0,0,nbsom
c$$$      write(1,*)' CAS  NO :    0   PAS NUMERO    :    0  ',
c$$$     &  'ITER. NO :    0             0.00000E+00'
c$$$
c$$$c ecriture de la solution
c$$$
c$$$      xpi=4.d0*datan(1.d0)
c$$$      xk=dsqrt(2.d0)/2.d0
c$$$      yk=dsqrt(2.d0)/2.d0
c$$$
c$$$      do is=1,nbsom
c$$$
c$$$c      write(*,*) is,xx(is),yy(is)
c$$$      xa=xx(is)
c$$$      ya=yy(is)
c$$$c      val=sol(is)-dsin(2.d0*xpi*(xk*xa+yk*ya))
c$$$      val=sol(is)
c$$$
c$$$      write(1,'(2i5,4e12.5)') is,4,xx(is),yy(is),0.,val
c$$$      
c$$$      enddo
c$$$
c$$$      close(1)

c calcul du gradient du potentiel aux points du fichier
c 'pcal.dat'

c longueur du domaine de calcul (adim)
c      along = 13.d0
c hauteur d'eau à l'infini (adim)
c      h1=1.d0
c pas du maillage en x
      pasx = along/nhmai

c distance et vitesse caracteristiques
      dchar=0.310d0
      vchar=dsqrt(9.81d0*dchar)

      open(17,file='pcal.dat',status='old')
      open(18,file='plotuv.dat',status='unknown')

c nombre de points où calculer le gradient
      read(17,*) nuser
c      write(18,*) nuser

      do iu=1,nuser

         read(17,*) if,xuser,yuser

c adimensionnement
         xuser=xuser/dchar
         yuser=yuser/dchar

c         write(*,*) xuser,yuser

c         recadrage en abscisse
         xsauv=xuser
         ysauv=yuser

         if (xuser.lt.0) then
            xuser = -xuser
c            write(*,*) xuser,iu
c            stop
         endif

         if (xuser.gt.along) then
            xuser = along-1.d-8
         endif

         if (yuser.lt.-1) then
            write(*,*) 'erreur pt n°',iu,':   y=',yuser*dchar,dchar
            stop
         endif



c recherche de l'élément dans lequel se trouve le point

         i=int(xuser/2.d0/pasx)
         kuser = 0
         do j=0,nvmai/2-1
            k= j+(nvmai/2)*i+1
            if (k.gt.nbel) stop
            x1=xx(ie(1,k))
            y1=yy(ie(1,k))
            x2=xx(ie(3,k))
            y2=yy(ie(3,k))
            x3=xx(ie(5,k))
            y3=yy(ie(5,k))
            x4=xx(ie(7,k))
            y4=yy(ie(7,k))

            t=(xuser-x1)/(x2-x1)
            yb= t*y2+(1-t)*y1
            yh= t*y3+(1-t)*y4
c            write(*,*) k,x1,y1,x2,y2,x3,y3,x4,y4

            if (yuser.le.yh.and.yuser.ge.yb) then
               kuser = k
c               write(*,*) k
            endif
         end do
         if (kuser.eq.0) then
          write(*,*) 'pt n°',iu,' dans l''air ',
     & 'x=',xsauv*dchar,' y=',ysauv*dchar
         endif

         if (kuser.ne.0) then
          write(*,*) 'pt n°',iu,' dans l''eau ',
     & 'x=',xsauv*dchar,' y=',ysauv*dchar
            
            hmax=yh
c     write(*,*) hmax
c     recherche des coordonnees dans l elt de ref
c     par newton
            
            x = .5d0
            y = .5d0
            
c     image par tau de xchap,ychap
c     fonction de base
            
 111        continue
c     write(*,*) x,y,det,kuser
            
            call fbase(vni,x,y,ajac,ajinv,det,vnivrai,ie,xx,yy,kuser)
            
            
            taux=0
            tauy=0
            
            do in=1,9
               taux=taux+vni(in,1)*xx(ie(in,kuser))
               tauy=tauy+vni(in,1)*yy(ie(in,kuser))
            end do
            
            dx = ajinv(1,1)*(xuser-taux)+ajinv(1,2)*(yuser-tauy)
            dy = ajinv(2,1)*(xuser-taux)+ajinv(2,2)*(yuser-tauy)
            
            x=x+dx
            y=y+dy
            
            
            if (dsqrt(dx*dx+dy*dy).ge.1.d-10) goto 111

c    test
            if (x.gt.1.01d0.or.y.gt.1.01d0
     $           .or.x.lt.-0.01d0.or.y.lt.-0.01d0) then
               write(*,*) x,y
               stop
            endif
            
c     calcul de phi
            
            valphi=0;
            gx=0;
            gy=0;
            do in=1,9
               valphi=valphi+vni(in,1)*sol(ie(in,kuser))         
               gx=gx+vnivrai(in,2)*sol(ie(in,kuser))         
               gy=gy+vnivrai(in,3)*sol(ie(in,kuser))         
            end do
            
         else
            
            valphi=0
            gx=0
            gy=0
            
         endif
         
         if (xsauv.lt.0) then
            valphi=-valphi
            gy=-gy
         endif
         
c     calcul de la pression par Bernouilli
         xpinf=1.d5
         rho=1000
         gpes=9.81
c         vref=dsqrt(gpes*h1)
         vref=1
         
         hinf=yy(ie(7,nbel))
         
         xpress = xpinf +rho*gpes*(hmax-yuser)
c     $        -rho*0.5d0*(gx*gx+gy*gy)*vref**2
         
         xpress = xpress
         if (kuser.eq.0) then
            xpress = 1d5
            rho =1
         endif
         
         write(18,*) iu,xsauv*dchar,ysauv*dchar
     $        ,rho,gx*vchar,gy*vchar,xpress
         
         
         
      end do
      
      close(17)
      close(18)

c tracés tecplot pour vérifier
c$$$      open(19,file='soliton.dat',status='unknown')
c$$$
c$$$c element de reference
c$$$
c$$$      xref(1)=0
c$$$      xref(2)=.5d0
c$$$      xref(3)=1.d0
c$$$      xref(4)=1.d0
c$$$      xref(5)=1.d0
c$$$      xref(6)=.5d0
c$$$      xref(7)=0
c$$$      xref(8)=0
c$$$      xref(9)=.5d0
c$$$
c$$$      yref(1)=0
c$$$      yref(2)=0
c$$$      yref(3)=0
c$$$      yref(4)=.5d0
c$$$      yref(5)=1.d0
c$$$      yref(6)=1.d0
c$$$      yref(7)=1.d0
c$$$      yref(8)=.5d0
c$$$      yref(9)=.5d0
c$$$
c$$$      do iel=1,nbel
c$$$         do ii=1,7,2
c$$$            
c$$$            call fbase(vni,xref(ii),yref(ii),ajac,
c$$$     $           ajinv,det,vnivrai,ie,xx,yy,iel)
c$$$            gx=0
c$$$            gy=0
c$$$            do in=1,9
c$$$               gx=gx+vnivrai(in,2)*sol(ie(in,iel))
c$$$               gy=gy+vnivrai(in,3)*sol(ie(in,iel))
c$$$            end do
c$$$            
c$$$            
c$$$            write(19,*) xx(ie(ii,iel)),yy(ie(ii,iel)),
c$$$     $           sol(ie(ii,iel)),gx,gy
c$$$         enddo 
c$$$      enddo
c$$$
c$$$      do iel=1,nbel
c$$$         write(19,*) (iel-1)*4+1,(iel-1)*4+2,(iel-1)*4+3,(iel-1)*4+4
c$$$      enddo
c$$$
c$$$      close(19)
      return
      end

      subroutine sol  (vkgs,vkgd,vkgi,vfg,kld,vu,neq,mp,ifac,isol,
     $  nsym,energ,ier)

c   resolution d'un systeme lineaire symetrique ou non. la matrice est
c   stockee par ligne de ciel,en memoire dans les tables vkgs,vkgd,vkgi
c
c       entrees
c          vkgs,vkgd,vkgi    matrice du systeme : parties superieure,
c                            diagonale, inferieure (double precision)
c          vfg               second membre
c          kld               pointeurs vers les hauts de colonne
c          vu                vecteur solution (qui peut etre vfg)
c          neq               nombre d'equations
c          mp                unite logique d'impression
c          ifac              si ifac.eq.1 triangularisation de
c                            la matrice
c          isol              si isol.eq.1 calcul de la solution a
c                            partir de la matrice triangularisee
c          nsym              indice de probleme non symetrique
c       sorties
c          vkgs,vkgd,vkgi    matrice triangularisee (si ifac.eq.1)
c          vfg               solution (si isol.eq.1)
c          energ             energie du systeme (si nsym.eq.0)
c          ier               mis a 1 si pivot nul rencontre
c
c=========================== debut des declarations ====================
      implicit real*8 (a-h,o,q-z)
      implicit integer (p)
      dimension vkgs(1),vkgd(1),vkgi(1),vfg(1),kld(1),vu(1)
      data vzero/0.d0/ 
c=========================== debut du code executable ==================
c
c-------  traitement
c
      ik=1
      if(vkgd(1).eq.vzero) goto 800
      energ=vzero
      ier=0
      if (isol.eq.1) then
        do i = 1, neq
          vu(i) = vfg(i)
        end do
      end if
c
c-------  pour chaque colonne ik a modifier
c
      jhk=1
      do 100 ik=2,neq
c
c-------  pointeur du haut de la colonne suivante ik+1
c
      jhk1=kld(ik+1)
c
c-------  hauteur de la colonne ik (hors termes superieur et diagonal)
c
      lhk=jhk1-jhk
      lhk1=lhk-1
c
c-------  ligne du premier terme a modifier dans la colonne ik
c
      imin=ik-lhk1
      imin1=imin-1
c
c-------  ligne du dernier terme a modifier dans la colonne ik
c
      imax=ik-1
      if(lhk1.lt.0) goto 100
      if(ifac.ne.1) goto 90
      if(nsym.eq.1) vkgi(jhk)=vkgi(jhk)/vkgd(imin1)
      if(lhk1.eq.0) goto 40
c
c-------  modifier les termes non diagonaux de la colonne ik
c
      jck=jhk+1
      jhj=kld(imin)
c
c-------  pour chaque terme place en jck, correspondant a la colonne ij
c
      do 30 ij=imin,imax
      jhj1=kld(ij+1)
c
c-------  nombre de termes modificatifs du terme place en jck
c
      ic=min0(jck-jhk,jhj1-jhj)
      if(ic.le.0.and.nsym.eq.0) goto 20
      c1=vzero
      if(ic.le.0) goto 17
      j1=jhj1-ic
      j2=jck-ic
      if(nsym.eq.1) goto 15
      vkgs(jck)=vkgs(jck)-scal(vkgs(j1),vkgs(j2),ic)
      goto 20
15    vkgs(jck)=vkgs(jck)-scal(vkgi(j1),vkgs(j2),ic)
      c1=scal(vkgs(j1),vkgi(j2),ic)
17    vkgi(jck)=(vkgi(jck)-c1)/vkgd(ij)
20    jck=jck+1
30    jhj=jhj1
c
c-------  modifier le terme diagonal
c
40    jck=jhk
      cdiag=vzero
      do 70 ij=imin1,imax
      c1=vkgs(jck)
      if(nsym.eq.1) goto 50
      c2=c1/vkgd(ij)
      vkgs(jck)=c2
      goto 60
50    c2=vkgi(jck)
60    cdiag=cdiag+c1*c2
70    jck=jck+1
      vkgd(ik)=vkgd(ik)-cdiag
      if (vkgd(ik).eq.0.) goto 800
c
c-------  resolution du systeme triangulaire inferieur
c
90    if(isol.ne.1) goto 100
      if(nsym.ne.1) vu(ik)=vfg(ik)-scal(vkgs(jhk),vu(imin1),lhk)
      if(nsym.eq.1) vu(ik)=vfg(ik)-scal(vkgi(jhk),vu(imin1),lhk)
100   jhk=jhk1
      if(isol.ne.1) goto 9999
c
c-------  resolution du systeme diagonal :
c
      if(nsym.eq.1) goto 120
      do 110 ik=1,neq
      c1=vkgd(ik)
      if (c1.eq.vzero) goto 800
      c2=vu(ik)/c1
      vu(ik)=c2
110   energ=energ+c1*c2*c2
c
c-------  resolution du systeme triangulaire superieur
c
120   ik=neq+1
      jhk1=kld(ik)
130   ik=ik-1
      if(nsym.eq.1) vu(ik)=vu(ik)/vkgd(ik)
      if(ik.eq.1) goto 9999
      c1=vu(ik)
      jhk=kld(ik)
      jbk=jhk1-1
      if(jhk.gt.jbk)goto 150
      ij=ik-jbk+jhk-1
      do 140 jck=jhk,jbk
      vu(ij)=vu(ij)-vkgs(jck)*c1
140   ij=ij+1
150   jhk1=jhk
      goto 130
c
c-------  erreurs
c
800   if (mp.ne.0) write(mp,8000) ik
8000  format(' * sol pivot nul equation',i5)
      ier=1
      goto 9999
c
c-------  fin
c
9999  continue
      return
c===========================   fin du module sol    ==================
      end

      function scal(x,y,n)
c=======================================================================
c calcul du produit scalaire
c=======================================================================
      implicit real*8 (a-h,o-z)
      dimension x(1),y(1)
      data zero/0.d0/
c-----------------------------------------------------------------------
      scal=zero
      do i=1,n
      scal=scal+x(i)*y(i)
      enddo
      return
      end

      subroutine fbase(vni,x,y,ajac,ajinv,det,vnivrai,ie,xx,yy,k)

      implicit real*8(c)
      implicit real*8(a-b,d-h,o-z)

      dimension vni(9,3),ajac(2,2),ajinv(2,2)
      dimension vnivrai(9,3),ie(9,*),xx(*),yy(*)

      vni(1,1)= 0.4D1*(x-0.5D0)*(x-1.D0)*(y-0.5D0)*(y-1.D0)
      vni(2,1)= -0.8D1*x*(x-1.D0)*(y-0.5D0)*(y-1.D0)
      vni(3,1)= 0.4D1*x*(x-0.5D0)*(y-0.5D0)*(y-1.D0)
      vni(4,1)= -0.8D1*x*(x-0.5D0)*y*(y-1.D0)
      vni(5,1)= 0.4D1*x*(x-0.5D0)*y*(y-0.5D0)
      vni(6,1)= -0.8D1*x*(x-1.D0)*y*(y-0.5D0)
      vni(7,1)= 0.4D1*(x-0.5D0)*(x-1.D0)*y*(y-0.5D0)
      vni(8,1)= -0.8D1*(x-0.5D0)*(x-1.D0)*y*(y-1.D0)
      vni(9,1)= 0.16D2*x*(x-1.D0)*y*(y-1.D0)
     
c derivee /xi de la fonction au point de gauss

      vni(1,2)= 0.4D1*(x-1.D0)*(y-0.5D0)*(y-1.D0)+0.4D1*(x-0.5D0)*(y
     #-0.5D0)*(y-1.D0)
      vni(2,2)= -0.8D1*(x-1.D0)*(y-0.5D0)*(y-1.D0)-0.8D1*x*(y-0.5D0)
     #*(y-1.D0)
      vni(3,2)= 0.4D1*(x-0.5D0)*(y-0.5D0)*(y-1.D0)+0.4D1*x*(y-0.5D0)
     #*(y-1.D0)
      vni(4,2)= -0.8D1*(x-0.5D0)*y*(y-1.D0)-0.8D1*x*y*(y-1.D0)
      vni(5,2)= 0.4D1*(x-0.5D0)*y*(y-0.5D0)+0.4D1*x*y*(y-0.5D0)
      vni(6,2)= -0.8D1*(x-1.D0)*y*(y-0.5D0)-0.8D1*x*y*(y-0.5D0)
      vni(7,2)= 0.4D1*(x-1.D0)*y*(y-0.5D0)+0.4D1*(x-0.5D0)
     # *y*(y-0.5D0)
      vni(8,2)= -0.8D1*(x-1.D0)*y*(y-1.D0)-0.8D1*(x-0.5D0)
     #*y*(y-1.D0)
      vni(9,2)= 0.16D2*(x-1.D0)*y*(y-1.D0)+0.16D2*x*y*(y-1.D0)

c derivee /eta de la fonction au point de gauss


      vni(1,3)= 0.4D1*(x-0.5D0)*(x-1.D0)*(y-1.D0)+0.4D1*(x-0.5D0)*(x
     #-1.D0)*(y-0.5D0)
      vni(2,3)=-0.8D1*x*(x-1.D0)*(y-1.D0)-0.8D1*x*(x-1.D0)*(y-0.5D0)
      vni(3,3)=0.4D1*x*(x-0.5D0)*(y-1.D0)+0.4D1*x*(x-0.5D0)
     #*(y-0.5D0)
      vni(4,3)= -0.8D1*x*(x-0.5D0)*(y-1.D0)-0.8D1*x*(x-0.5D0)*y
      vni(5,3)= 0.4D1*x*(x-0.5D0)*(y-0.5D0)+0.4D1*x*(x-0.5D0)*y
      vni(6,3)= -0.8D1*x*(x-1.D0)*(y-0.5D0)-0.8D1*x*(x-1.D0)*y
      vni(7,3)= 0.4D1*(x-0.5D0)*(x-1.D0)*(y-0.5D0)+0.4D1*(x-0.5D0)*(
     #x-1.D0)*y
      vni(8,3)= -0.8D1*(x-0.5D0)*(x-1.D0)*(y-1.D0)-0.8D1*(x-0.5D0)*(
     #x-1.D0)*y
      vni(9,3)= 0.16D2*x*(x-1.D0)*(y-1.D0)+0.16D2*x*(x-1.D0)*y

c calcul de la jacobienne

      ajac(1,1)=0
      ajac(1,2)=0
      ajac(2,1)=0
      ajac(2,2)=0

      do in=1,9
         ajac(1,1)=ajac(1,1)+vni(in,2)*xx(ie(in,k))
         ajac(1,2)=ajac(1,2)+vni(in,3)*xx(ie(in,k))
         ajac(2,1)=ajac(2,1)+vni(in,2)*yy(ie(in,k))
         ajac(2,2)=ajac(2,2)+vni(in,3)*yy(ie(in,k))
         
      end do

      det=ajac(1,1)*ajac(2,2)-ajac(1,2)*ajac(2,1)

      if (det.le.0.d0) then
      write(*,*) 'det negatif',det
      stop
      endif

c calcul de la jacobienne inverse

      ajinv(1,1)=ajac(2,2)/det
      ajinv(2,2)=ajac(1,1)/det
      ajinv(1,2)=-ajac(1,2)/det
      ajinv(2,1)=-ajac(2,1)/det

c calcul du gradient dans les bonnes variables

      do ii=1,9
         vnivrai(ii,1)=vni(ii,1)
         vnivrai(ii,2)= ajinv(1,1) * vni(ii,2) + ajinv(2,1) * vni(ii,3)
         vnivrai(ii,3)= ajinv(1,2) * vni(ii,2) + ajinv(2,2) * vni(ii,3)
      enddo

      return
      end
