!> \file courbures.f90 !! Implementation du subroutine courbure !< !> SUBROUTINE: courbure !! \author ... !! \date ... !! @note Subroutine qui permet de trouver la surface de type !! @note A_2D(x,y)=ax2+by2+cxy+dx+ey+f !! @note passant au plus pres (moindre carre) des 9 points d'une sous grille !! @note Ensuite, calcule les diverses courbures, la pente et son azimuth !! @note Le pas de la grille est 1 !! @note Fait appel a sgelsy de Lapack qui minimise || A * X -B || !> subroutine courbure(nx,ny,dy,A_2D,A_new,pente_A,teta_A,crx,cry,crxy) implicit none integer mnode,ncoef,npts ! entree routine integer,intent(in) :: nx,ny !< dimensions des matrices real,intent(in) :: dy !< maille physique real,dimension(nx,ny),intent(in) :: A_2D !< la surface dont on cherche la courbure ! sortie routine real,dimension(nx,ny),intent(out):: A_new !< altitude finale du point real,dimension(nx,ny),intent(out):: pente_A !< pente_A (amplitude) real,dimension(nx,ny),intent(out):: teta_A !< direction de la pente real,dimension(nx,ny),intent(out):: crx !< courbure dans la direction de la pente real,dimension(nx,ny),intent(out):: cry !< courbure dans la direction transversale real,dimension(nx,ny),intent(out):: crxy !< courbure croisee ! variables locales parameter(mnode=9) !< nombre de points de la sous-grille parameter (ncoef=6) !< nombre de coefficients de A_2D parameter (npts=60000) !< nombre de points a traiter (nx*ny) real :: dx_ls !< pas de la grille dans le moindre carre integer :: i,j,k,pp,iz,jz,iz1,jz1 !< compteurs real :: x,y !< coordonnees dans la sous-grille real :: aaa !< pour le calcul de l'azimut ! specifiques pour appel a sgelsy real :: Blin(mnode,npts) !< en input c'est B(mnoeud,npts) !< en output c'est X(ncoef,npts) real :: A(mnode,ncoef) !< matrice A real :: rcond !< pour appel a sgelsy ~ max(A)/min(A) integer :: rank !< pour sgelsy (output) -> rang de R11 integer :: lwork !< taille du tableau de travail parameter(lwork = 200000) !< voir sgelsy real,dimension(lwork) :: work !< tableau de travail integer,dimension(ncoef):: jpvt !< vecteur entier integer :: info !< retour des info integer :: mm,nn,npp,lw mm=mnode ! pour ne pas passer des nn=ncoef ! constantes aux routines npp=npts lw=lwork dx_ls=1 ! pas utilise dans le moindre carre ! remplissage de la matrice A (la meme pour tous les points) k=0 do j=-1,1 ! 7 8 9 do i=-1,1 ! balaye les i puis les j 4 5 6 k=k+1 ! k varie de 1 a 9 1 2 3 x=i*dx_ls y=j*dx_ls A(k,1)=x*x A(k,2)=y*y A(k,3)=x*y A(k,4)=x A(k,5)=y A(k,6)=1 end do end do rcond=1.e-30 ! boucle sur les points pp=0 ij_loop: do j=1,ny pp_loop: do i=1,nx pp=pp+1 ! pp varie de 1 a nx*ny ! Construction de la sous-grille centree sur i,j k=0 do jz=j-1,j+1 do iz=i-1,i+1 iz1=max(1,iz) ! pour éviter les débordements du tableau A_2D iz1=min(nx,iz1) jz1=max(1,jz) jz1=min(ny,jz1) k=k+1 ! k varie de 1 a 9 Blin(k,pp)=A_2D(iz1,jz1) end do end do end do pp_loop end do ij_loop !____________________________________________________________________ ! ! appel a la subroutine sgelsy !____________________________________________________________________ call sgelsy(mm,nn,npp,a,mm,blin,mm,jpvt,rcond,rank,work,lw,info) ! Maintenant les coefficients sont dans Blin(1:6,pp) ! Blin(1,pp)= a coeff de x2 ! Blin(2,pp)= b coeff de y2 ! Blin(3,pp)= c coeff de xy ! Blin(4,pp)= d coeff de x ! Blin(5,pp)= e coeff de y ! Blin(6,pp)= f terme constant pp=0 grille_loop: do j=1,ny do i=1,nx pp=pp+1 A_new(i,j)=Blin(6,pp) ! f pente_A(i,j)=Blin(4,pp)*Blin(4,pp)+Blin(5,pp)*Blin(5,pp) ! d2+e2 pente_sup0: if ((pente_A(i,j).gt.1.e-5) & .and.(i.ne.1).and.(i.ne.nx).and.(j.ne.1).and.(j.ne.ny)) then ! pente pente_A(i,j)=(pente_A(i,j)**0.5) ! attention ce calcul etait avant fait ! apres le calcul de courbure. ! Par comparaison avec Frederique, c'etait un bug. ! courbx crx(i,j)=blin(1,pp)*blin(4,pp)*blin(4,pp) & ! ad2 + blin(2,pp)*blin(5,pp)*blin(5,pp) & ! be2 + blin(3,pp)*blin(4,pp)*blin(5,pp) ! cde crx(i,j)=2.*crx(i,j)/pente_A(i,j)/dy/dy ! dy2 pour ! passage en dim reelle ! courby cry(i,j)=blin(1,pp)*blin(5,pp)*blin(5,pp) & ! ae2 + blin(2,pp)*blin(4,pp)*blin(4,pp) & ! bd2 - blin(3,pp)*blin(4,pp)*blin(5,pp) ! -cde cry(i,j)=2.*cry(i,j)/pente_A(i,j)/dy/dy ! dy2 pour ! passage en dim reelle ! courbxy crxy(i,j)=2.*blin(1,pp)*blin(4,pp)*blin(5,pp) & ! 2ade - 2.*blin(2,pp)*blin(4,pp)*blin(5,pp) & ! -2bed - blin(3,pp)*blin(4,pp)*blin(4,pp) & ! -cd2 - blin(3,pp)*blin(5,pp)*blin(5,pp) ! -ce2 crxy(i,j)=crxy(i,j)/pente_A(i,j)/dy/dy ! dy2 pour ! passage en dim reelle ! direction de la pente aaa=blin(5,pp)/pente_A(i,j) ! pente y / pente if ((aaa.gt.-1.).and.(aaa.lt.1.)) then teta_A(i,j)=asin(aaa)*180./3.14159 else if (aaa.ge.1.) then teta_A(i,j)=90. else if (aaa.le.-1.) then teta_A(i,j)=-90. endif if (blin(4,pp).lt.0.) then teta_A(i,j)=180.-teta_A(i,j) endif else ! pente nulle crx(i,j)=0. cry(i,j)=0. crxy(i,j)=0. pente_A(i,j)=0. teta_A(i,j)=0. endif pente_sup0 ! Fin du test sur la pente ! pente pente_A(i,j)=pente_A(i,j)/dy ! pour passage en dimension reelle end do end do grille_loop ! fin de la boucle sur les points de grille return end subroutine courbure