source: branches/iLoveclim/SOURCES/courbures.f90 @ 244

!> \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