!PROGRAMA PRINCIPAL
PROGRAM PROGRAMAPRINCIPAL
USE LAPACK95
IMPLICIT NONE
INTEGER,PARAMETER :: LONG=8!selected_int_kind(15)
INTEGER :: i,j,ppx,ppy,m,n,t,cont,k,cont1,total,cargasdistribuidas,fila,w
!INTEGER,DIMENSION(4) :: CC = (/1,1,1,1/) 
REAL (KIND=LONG) :: E,poisson,D,D1,espesor,alpha,a=0.,b=0.
REAL (KIND=LONG) :: pasox
REAL (KIND=LONG),ALLOCATABLE,DIMENSION(:,:) :: s,MX,MY,MT
REAL (KIND=LONG),ALLOCATABLE,DIMENSION(:) :: P,X
REAL (KIND=LONG),ALLOCATABLE,DIMENSION(:,:) :: coef
INTEGER,ALLOCATABLE,DIMENSION(:) :: ipiv
poisson=0.2 !se debe pedir en el java
E=29240177.4 !kN/m^2
espesor=0.3!revisar
pasox=0.1 !viene del mallado
alpha=1.
cont=1
k=22952.96 !kN/m^3
D=E*(espesor**3)/(12*(1-(poisson**2)))
WRITE(*,*)'D=',D,'kN/m^3'
D1=k/D
cont1=0
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!BORRAR LOS VALORES DESPUÉS!!!!!!!!!!!
CALL PASO(pasox)
WRITE(*,*)'TRAS SALIR DE LA SUBRUTINA DE PASO',pasox
CALL GENERACIONMALLADO(a,b,pasox,cargasdistribuidas)
!a=1
!b=1
!cargasdistribuidas=64 
write(*,*)'a=',a,'b=',b,'pasox=',pasox,'cargasdistribuidas',cargasdistribuidas
write(*,1005)pasox
1005 FORMAT(F10.3)
READ(*,*)
ppx=NINT((a/pasox)+1) !Debe venir implícito desde la funcion de mallado
ppy=NINT((b/pasox)+1)
WRITE(*,*)'PPX=',ppx,'ppy=',ppy
total=ppx*ppy
ALLOCATE(s(ppy,ppx))
ALLOCATE(MX(ppy,ppx))
ALLOCATE(MY(ppy,ppx))
ALLOCATE(MT(ppy,ppx))
ALLOCATE(P(total))
ALLOCATE(X(total))
ALLOCATE(ipiv(total))
ALLOCATE(coef(ppx*ppy,ppx*ppy))
s=0.0
P=0.0
!OPEN(10,FILE="/home/Placas/mimatriz.dat")
!OPEN(11,FILE="/home/Placas/mimatrizcoef.dat")
!OPEN(12,FILE="/home/Placas/mimatrizcoefdef.dat")

FILAS: DO m=1,ppy
  COLUMNAS: DO n=1,ppx
  write(*,*)'Vuelta=',cont1
  cont1=cont1+1
   if(MOD(cont1,1000)==2) then
   write(*,*)'Vuelta=',cont1
   end if
	    !Estudiamos primero los puntos centrales
	 centrales:IF (((n>1).AND.(n<(ppx-1))).AND.((m>1).AND.(m<(ppy-1)))) THEN !Todo lo que debía ser dividido por paso^4 está   
	           !Línea central
	           s(m,n)=((6+6*(alpha)**4+8*(alpha**2))/(pasox**4))+D1
	           s(m,n+1)=(-4*(1+(alpha**2)))/(pasox**4)
	           s(m,n+2)=1/(pasox**4)
	           s(m,n-1)=-4*(1+(alpha**2))/(pasox**4)
	           s(m,n-2)=1/(pasox**4)
	           !Línea por encima de la central
	           s(m-1,n)=-4*(alpha**2)*(1+(alpha**2))/(pasox**4)
	           s(m-1,n+1)=2*(alpha**2)/(pasox**4)
	           s(m-1,n-1)=2*(alpha**2)/(pasox**4)
	           !Línea por debajo de la central
	           s(m+1,n)=-4*(alpha**2)*(1+(alpha**2))/(pasox**4)
	           s(m+1,n+1)=2*(alpha**2)/(pasox**4)
	           s(m+1,n-1)=2*(alpha**2)/(pasox**4)
	           !Líneas superior e inferior
	           s(m+2,n)=alpha**4/(pasox**4)
	           s(m-2,n)=alpha**4/(pasox**4)
	           END IF centrales
!!!!!!!!!CONDICIONES DE CONTORNO!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!CONDICIÓN DE CONTORNO INFERIOR
	!FILA ANTERIOR AL BORDE INFERIOR (PUNTO IZQUIERDO)
	alinfizq:IF((m==(ppy-1)).AND.(n==2)) THEN
	      !Línea central
	      s(m,n)=((5+5*(alpha**4)+8*(alpha**2))/(pasox**4))+D1
	      s(m,n-1)=(-2*(2*(alpha**2)-poisson*(alpha**2)+1))/(pasox**4)
	      s(m,n+1)=(-4*(1+alpha**2)/(pasox**4))
	      s(m,n+2)=(1.)/(pasox**4)
	      !Línea por encima de la línea central
	      s(m-1,n-1)=(alpha**2)*(2-poisson)/(pasox**4)
	      s(m-1,n)=-4*(alpha**2)*(1+(alpha**2))/(pasox**4)
	      s(m-1,n+1)=2*(alpha**2)/(pasox**4)
	      !Línea superior
	      s(m-2,n)=alpha**4/(pasox**4) !DEBERÍA IR UN MENOS Y HABÍA UN MAS
	      !Línea por debajo de la línea central
	      s(m+1,n)=-2*(alpha**2)*(2-poisson+(alpha**2))/(pasox**4)
	      s(m+1,n+1)=(alpha**2)*(2-poisson)/(pasox**4)
	      s(m+1,n-1)=2*(alpha**2)*(1-poisson)/(pasox**4)
	END IF alinfizq
	!FILA ANTERIOR AL BORDE INFERIOR (PUNTO DERECHO)
	alinfder:IF((m==(ppy-1)).AND.(n==ppx-1)) THEN
	      !Línea central
	      s(m,n)=((5+5*(alpha**4)+8*(alpha**2))/(pasox**4))+D1
	      s(m,n+1)=-2*(2*(alpha**2)-poisson*(alpha**2)+1)/(pasox**4)
	      s(m,n-1)=-4*(1+alpha**2)/(pasox**4)
	      s(m,n-2)=1./(pasox**4)
	      !Línea por encima de la línea central
	      s(m-1,n+1)=(alpha**2)*(2-poisson)/(pasox**4)
	      s(m-1,n)=-4*(alpha**2)*(1+(alpha**2))/(pasox**4)
	      s(m-1,n-1)=2*(alpha**2)/(pasox**4)
	      !Línea superior
	      s(m-2,n)=alpha**4/(pasox**4)
	      !Línea por debajo de la línea central
	      s(m+1,n)=-2*(alpha**2)*(2-poisson+(alpha**2))/(pasox**4)
	      s(m+1,n-1)=(alpha**2)*(2-poisson)/(pasox**4)
	      s(m+1,n+1)=2*(alpha**2)*(1-poisson)/(pasox**4)
	END IF alinfder
!CONDICIÓN DE CONTORNO 
	!FILA ANTERIOR AL BORDE INFERIOR (PUNTO IZQUIERDO)	
	alsupizq: IF ((m==2).AND.(n==2)) THEN
	!Línea central
          s(m,n)=((5+5*(alpha**4)+8*(alpha**2))/(pasox**4))+D1
	      s(m,n-1)=-2*(2*(alpha**2)-poisson*(alpha**2)+1)/(pasox**4)
	      s(m,n+1)=-4*(1+alpha**2)/(pasox**4)
	      s(m,n+2)=1./(pasox**4)
	!Línea por encima de la línea central
	      s(m-1,n+1)=(alpha**2)*(2-poisson)/(pasox**4)
	      s(m-1,n)=-2*(alpha**2)*(2-poisson+(alpha**2))/(pasox**4)
	      s(m-1,n-1)=2*(alpha**2)*(1-poisson)/(pasox**4)
	!Línea por debajo de la línea central
		  s(m+1,n)=-4*(alpha**2)*(1+(alpha**2))/(pasox**4)
	      s(m+1,n-1)=(alpha**2)*(2-poisson)/(pasox**4)
	      s(m+1,n+1)=2*(alpha**2)/(pasox**4)
	  !Línea superior
	      s(m+2,n)=alpha**4/(pasox**4)
	END IF alsupizq
	alsupder: IF((m==(2)).AND.(n==(ppx-1))) THEN!REVISAR ECUACIONES 
	!Línea central
          s(m,n)=((5+5*(alpha**4)+8*(alpha**2))/(pasox**4))+D1
	      s(m,n+1)=-2*(2*(alpha**2)-poisson*(alpha**2)+1)/(pasox**4)
	      s(m,n-1)=-4*(1+alpha**2)/(pasox**4)
	      s(m,n-2)=1./(pasox**4)
	!Línea por encima de la línea central
	      s(m-1,n-1)=(alpha**2)*(2-poisson)/(pasox**4)
	      s(m-1,n)=-2*(alpha**2)*(2-poisson+(alpha**2))/(pasox**4)
	      s(m-1,n+1)=2*(alpha**2)*(1-poisson)/(pasox**4)
	!Línea por debajo de la línea central
		  s(m+1,n)=-4*(alpha**2)*(1+(alpha**2))/(pasox**4)
	      s(m+1,n+1)=(alpha**2)*(2-poisson)/(pasox**4)
	      s(m+1,n-1)=2*(alpha**2)/(pasox**4)
	!Línea superior
	      s(m+2,n)=alpha**4/(pasox**4)
	END IF alsupder
	!CONDICIONES DE CONTORNO EN LOS LADOS
	!LADO PARALELO AL IZQUIERDO,CORRIDA UNA POSICIÓN, n==2
	pladoizq: IF ((n==2).AND.((m>2).AND.(m<(ppy-1)))) THEN !CAMBIÉ LOS LÍMITES. Y LO VOLVÍ A DEJAR COMO ANTES
	!Línea central
	      s(m,n)=((5+6*(alpha**4)+8*(alpha**2))/(pasox**4))+D1
	      s(m,n+1)=-4*(1+(alpha**2))/(pasox**4)
	      s(m,n+2)=1./(pasox**4)
	      s(m,n-1)=-2*(2*(alpha**2)-poisson*(alpha**2)+1)/(pasox**4)
	!Línea por encima de la línea central
	      s(m-1,n)=-4*(alpha**2)*(1+alpha**2)/(pasox**4)
	      s(m-1,n+1)=2*(alpha**2)/(pasox**4)
	      s(m-1,n-1)=(alpha**2)*(2-poisson)/(pasox**4)
	!Línea por debajo de la línea central
          s(m+1,n)=-4*(alpha**2)*(1+alpha**2)/(pasox**4)
          s(m+1,n+1)=2*(alpha**2)/(pasox**4)
          s(m+1,n-1)=(alpha**2)*(2-poisson)/(pasox**4)
    !Línea superior
		  s(m-2,n)=alpha**4/(pasox**4)
	!Línea inferior
          s(m+2,n)=alpha**4/(pasox**4)
	END IF pladoizq
	
	ladoizq: IF (((m>2).AND.(m<(ppy-1))).AND.(n==1)) THEN !Revisados coeficientes Y TODAS LAS ECUACIONES
	!Línea central
	      s(m,n)=(((1+4*(alpha**2)*(1-poisson)+3*(alpha**4)*(1-(poisson**2)))*2.)/(pasox**4))+D1
	      s(m,n+1)=((-2*(2*(alpha**2)-poisson*(alpha**2)+1))*2.)/(pasox**4)
	      s(m,n+2)=((1.)*2.)/(pasox**4)
	!Línea por encima de la central
          s(m-1,n)=(((-2*(alpha**2)*(1-poisson+(alpha**2)*(1-(poisson**2))))*2.)/(pasox**4))
          s(m-1,n+1)=(((alpha**2)*(2-poisson))*2.)/(pasox**4)
    !Línea por debajo de la central
          s(m+1,n)=((-2*(alpha**2)*(1-poisson+(alpha**2)*(1-(poisson**2))))*2.)/(pasox**4)
          s(m+1,n+1)=(((alpha**2)*(2-poisson))*2.)/(pasox**4)
    !Línea superior
          s(m-2,n)= ((0.5*(alpha**4)*(1-(poisson**2)))*2.)/(pasox**4)
          s(m+2,n)= ((0.5*(alpha**4)*(1-(poisson**2)))*2.)/(pasox**4)
	END IF ladoizq
 !LADO PARALELO AL DERECHO,CORRIDA UNA POSICIÓN, n==ppx-1
	pladoder: IF (n==(ppx-1).AND.((m>2).AND.(m<(ppy-1)))) THEN !CAMBIÉ LOS LIMITES Y LOS VOLVÍ A CAMBIAR
	!Línea central
	      s(m,n)=((5+6*(alpha**4)+8*(alpha**2))/(pasox**4))+D1
	      s(m,n-1)=-4*(1+(alpha**2))/(pasox**4)
	      s(m,n-2)=1./(pasox**4)
	      s(m,n+1)=-2*(2*(alpha**2)-poisson*(alpha**2)+1)/(pasox**4)
	 !Línea por encima de la línea central
	      s(m-1,n)=-4*(alpha**2)*(1+alpha**2)/(pasox**4)
	      s(m-1,n-1)=2*(alpha**2)/(pasox**4)
	      s(m-1,n+1)=(alpha**2)*(2-poisson)/(pasox**4)
	 !Línea por debajo de la línea central
          s(m+1,n)=-4*(alpha**2)*(1+alpha**2)/(pasox**4)
          s(m+1,n-1)=2*(alpha**2)/(pasox**4)
          s(m+1,n+1)=(alpha**2)*(2-poisson)/(pasox**4)
	 !Línea superior
	      s(m-2,n)=alpha**4/(pasox**4)
	 !Línea inferior
	      s(m+2,n)=alpha**4/(pasox**4)
	END IF pladoder   
 !LADO DERECHO
   ladoder: IF (((m>2).AND.(m<(ppy-1))).AND.(n==ppx)) THEN !CAMBIÉ LOS LIMITES Y LOS VOLVI A PONER IGUAL
	!Línea central
	      s(m,n)=(((1+4*(alpha**2)*(1-poisson)+3*(alpha**4)*(1-(poisson**2)))*2.)/(pasox**4))+D1
	      s(m,n-1)=((-2*(2*(alpha**2)-poisson*(alpha**2)+1))*2.)/(pasox**4)
	      s(m,n-2)=((1.)*2.)/(pasox**4)
	!Línea por encima de la central
          s(m-1,n)=((-2*(alpha**2)*(1-poisson+(alpha**2)*(1-(poisson**2))))*2.)/(pasox**4)
          s(m-1,n-1)=(((alpha**2)*(2-poisson))*2.)/(pasox**4)
    !Línea por debajo de la central
          s(m+1,n)=((-2*(alpha**2)*(1-poisson+(alpha**2)*(1-(poisson**2))))*2.)/(pasox**4)
          s(m+1,n-1)=(((alpha**2)*(2-poisson))*2.)/(pasox**4)
    !Línea superior
          s(m-2,n)= ((0.5*(alpha**4)*(1-poisson))*2.)/(pasox**4)
          s(m+2,n)= ((0.5*(alpha**4)*(1-poisson))*2.)/(pasox**4)
	END IF ladoder
 !LADO PARALELO AL SUPERIOR,CORRIDA UNA POSICIÓN, m==2
 pladosup: IF((m==2).AND.((n>2).AND.(n<(ppx-1)))) THEN !CAMBIÉ LOS LIMITES Y LOS VOLVÍ A CAMBIAR
    !Línea central     
          s(m,n)=(5+6*(alpha**4)+8*(alpha**2))/(pasox**4)+D1
          s(m,n+1)=(-4*(alpha**2)*(1+(alpha**2)))/(pasox**4)
          s(m,n+2)=alpha**4/(pasox**4)
          s(m,n-1)=(-4*(alpha**2)*(1+(alpha**2)))/(pasox**4)
          s(m,n-2)=(alpha**4)/(pasox**4)
   !Línea por encima de la central. Coincide con el borde superior
          s(m-1,n)=(-2*(2*(alpha**2)-poisson*(alpha**2)+1))/(pasox**4)
          s(m-1,n+1)=((alpha**2)*(2-poisson))/(pasox**4)
          s(m-1,n-1)=(alpha**2)*(2-poisson)/(pasox**4)
   !Línea por debajo de la central.
          s(m+1,n)=-4*(1+(alpha**2))/(pasox**4) !Había un menos y lo puse bien,debe ser un mas
          s(m+1,n+1)= 2*(alpha**2)/(pasox**4)
          s(m+1,n-1)=2*(alpha**2)/(pasox**4)
   !Línea inferior
          s(m+2,n)=1./(pasox**4)
	END IF pladosup
ladosup:IF((m==1).AND.((n>2).AND.(n<ppx-1))) THEN !CAMBIÉ LOS LIMITES Y LOS VOLVÍ A CAMBIAR
   !Línea central
          s(m,n)=(((1+4*(alpha**2)*(1-poisson)+3*(alpha**4)*(1-(poisson**2)))*2.)/(pasox**4))+D1
          s(m,n+1)=((-2*(alpha**2)*(1-poisson+(alpha**2)*(1-(poisson**2))))*2.)/(pasox**4)
          s(m,n-1)=((-2*(alpha**2)*(1-poisson+(alpha**2)*(1-(poisson**2))))*2.)/(pasox**4)
          s(m,n+2)=((0.5*(alpha**2)*(1-(poisson**2)))*2.)/(pasox**4)
          s(m,n-2)=((0.5*(alpha**2)*(1-(poisson**2)))*2.)/(pasox**4)
   !Línea por debajo de la central
          s(m+1,n)=((-2*(2*alpha**2-(poisson*(alpha**2))+1))*2.)/(pasox**4)
          s(m+1,n+1)=(((alpha**2)*(2-poisson))*2.)/(pasox**4)
          s(m+1,n-1)=(((alpha**2)*(2-poisson))*2.)/(pasox**4)
   !Línea inferior
          s(m+2,n)=1./(pasox**4)
	END IF ladosup
 !LADO PARALELO AL INFERIOR,CORRIDA UNA POSICIÓN, m==ppy-1
 pladoinf: IF((m==ppy-1).AND.((n>2).AND.(n<(ppx-1)))) THEN !CAMBIÉ LOS LIMITES Y LOS VOLVÍ A CAMBIAR
     
    !Línea central     
          s(m,n)=((5+6*(alpha**4)+8*(alpha**2))/(pasox**4))+D1
          s(m,n+1)=(-4*(alpha**2)*(1+(alpha**2)))/(pasox**4)
          s(m,n+2)=alpha**4/(pasox**4)
          s(m,n-1)=-4*(alpha**2)*(1+(alpha**2))/(pasox**4)
          s(m,n-2)=alpha**4/(pasox**4)
   !Línea por debajo de la central. Coincide con el borde inferior
          s(m+1,n)=-2*(2*(alpha**2)-poisson*(alpha**2)+1)/(pasox**4)
          s(m+1,n+1)=(alpha**2)*(2-poisson)/(pasox**4)
          s(m+1,n-1)=(alpha**2)*(2-poisson)/(pasox**4)
   !Línea por encima de la central. !he modificado estas lineas pues aparecia un mas (m+1) y queda fuera
          s(m-1,n)=-4*(1+(alpha**2))/(pasox**4)!Había un menos y lo puse bien,debe ser un mas
          s(m-1,n+1)= 2*(alpha**2)/(pasox**4)
          s(m-1,n-1)=2*(alpha**2)/(pasox**4)
   !Línea superior !he modificado estas lineas pues aparecia un mas (m+1) y queda fuera
          s(m-2,n)=1./(pasox**4)
	END IF pladoinf
!LADO INFERIOR
ladoinf:IF((m==ppy).AND.((n>2).AND.(n<ppx-1))) THEN !CAMBIÉ LOS LIMITES Y LOS VOLVÍ A CAMBIAR
   !Línea central
          s(m,n)=(((1+4*(alpha**2)*(1-poisson)+3*(alpha**4)*(1-(poisson**2)))*2.)/(pasox**4))+D1
          s(m,n+1)=((-2*(alpha**2)*(1-poisson+(alpha**2)*(1-(poisson**2))))*2.)/(pasox**4)
          s(m,n-1)=((-2*(alpha**2)*(1-poisson+(alpha**2)*(1-(poisson**2))))*2.)/(pasox**4)
          s(m,n+2)=((0.5*(alpha**2)*(1-(poisson**2)))*2.)/(pasox**4)
          s(m,n-2)=((0.5*(alpha**2)*(1-(poisson**2)))*2.)/(pasox**4)
   !Línea por encima de la central
          s(m-1,n)=((-2*(2*alpha**2-(poisson*(alpha**2))+1))*2.)/(pasox**4)
          s(m-1,n+1)=(((alpha**2)*(2-poisson))*2.)/(pasox**4)
          s(m-1,n-1)=(((alpha**2)*(2-poisson))*2.)/(pasox**4)
   !Línea inferior
          s(m-2,n)=(1.*2.)/(pasox**4)
	END IF ladoinf
!!PUNTOS PUNTOS CERCANOS A LAS ESQUINAS
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!PUNTO CERCANO A LA ESQUINA INFERIOR IZQUIERDA (EJE VERTICAL)

!cesqizqinfver
	IF ((m==(ppy-1)).AND.(n==1)) THEN 

   !Borde vertical correspondiente al borde izquierdo
          s(m,n)=(((1+4*(alpha**2)*(1-poisson)+2.5*(alpha**4)*(1-(poisson**2)))*2.)/(pasox**4))+D1
          s(m-1,n)=((-2*(alpha**2)*(1-poisson+(alpha**2)*(1-(poisson**2))))*2.)/(pasox**4)
          s(m-2,n)=((0.5*(alpha**4)*(1-(poisson**2)))*2.)/(pasox**4)
          s(m+1,n)=((-2*(alpha**2)*(1-poisson+0.5*(alpha**2)*(1-(poisson**2))))*2.)/(pasox**4)
   !Línea vertical paralela al borde izquierdo
          s(m,n+1)=((-2*(2*(alpha**2)-poisson*(alpha**2)+1))*2.)/(pasox**4) 
          s(m-1,n+1)=(((alpha**2)*(2.-poisson))*2.)/(pasox**4)
          s(m+1,n+1)=(((alpha**2)*(2.-poisson))*2.)/(pasox**4)
   !Línea más exterior
          s(m,n+2)=(1.*2.)/(pasox**4)

	END IF 
!!PUNTO CERCANO A LA ESQUINA SUPERIOR IZQUIERDA (EJE VERTICAL)
!cesqizqsupver
	IF ((m==2).AND.(n==1)) THEN !CAMBIÉ LOS LÍMITES.arregle el cuadrado en poisson
   !Borde vertical correspondiente al borde izquierdo
          s(m,n)=(((1+4*(alpha**2)*(1-poisson)+2.5*(alpha**4)*(1-(poisson**2)))*2.)/(pasox**4))+D1
          s(m-1,n)=(((-2*(alpha**2)*(1-poisson+0.5*(alpha**2)*(1-(poisson**2))))*2.)/(pasox**4))
          s(m+2,n)=(((0.5*(alpha**4)*(1-(poisson**2)))*2.)/(pasox**4))
          s(m+1,n)=((-2*(alpha**2)*(1-poisson+(alpha**2)*(1-(poisson**2))))*2.)/(pasox**4)
   !Línea vertical paralela al borde izquierdo
          s(m,n+1)=((-2*(2*(alpha**2)-poisson*(alpha**2)+1))*2.)/(pasox**4) 
          s(m-1,n+1)=(((alpha**2)*(2.-poisson))*2.)/(pasox**4)
          s(m+1,n+1)=(((alpha**2)*(2.-poisson))*2.)/(pasox**4)
   !Línea más exterior
          s(m,n+2)=(1.*2.)/(pasox**4)
	END IF 
!!PUNTO CERCANO A LA ESQUINA SUPERIOR DERECHA (EJE VERTICAL)
!cesqdersupver:
	IF ((n==(ppx)).AND.(m==2)) THEN !CAMBIÉ LOS LIMITES.arreglé el cuadrado en poisson
   !Borde vertical correspondiente al borde izquierdo
          s(m,n)=(((1+4*(alpha**2)*(1-poisson)+2.5*(alpha**4)*(1-(poisson**2)))*2.)/(pasox**4))+D1
          s(m-1,n)=((-2*(alpha**2)*(1-poisson+0.5*(alpha**2)*(1-(poisson**2))))*2.)/(pasox**4)
          s(m+2,n)=((0.5*(alpha**4)*(1-(poisson**2)))*2.)/(pasox**4)
          s(m+1,n)=((-2*(alpha**2)*(1-poisson+(alpha**2)*(1-(poisson**2))))*2.)/(pasox**4)
   !Línea vertical paralela al borde izquierdo
          s(m,n-1)=((-2*(2*(alpha**2)-poisson*(alpha**2)+1))*2.)/(pasox**4) 
          s(m-1,n-1)=(((alpha**2)*(2.-poisson))*2.)/(pasox**4)
          s(m+1,n-1)=(((alpha**2)*(2.-poisson))*2.)/(pasox**4)
   !Línea más exterior
          s(m,n-2)=(1.*2.)/(pasox**4)
	END IF 
!!PUNTO CERCANO A LA ESQUINA INFERIOR DERECHA (EJE VERTICAL)
!cesqderinfver:
	IF ((m==(ppy-1)).AND.(n==ppx)) THEN !CAMBIÉ LOS LIMITES.arreglé el cuadrado en poisson
   !Borde vertical correspondiente al borde izquierdo
          s(m,n)=(((1+4*(alpha**2)*(1-poisson)+2.5*(alpha**4)*(1-(poisson**2)))*2.)/(pasox**4))+D1
          s(m-1,n)=(((-2*(alpha**2)*(1-poisson+(alpha**2)*(1-(poisson**2))))*2.)/(pasox**4))
          s(m-2,n)=((0.5*(alpha**4)*(1-(poisson**2)))*2.)/(pasox**4)
          s(m+1,n)=((-2*(alpha**2)*(1-poisson+0.5*(alpha**2)*(1-(poisson**2))))*2.)/(pasox**4)
   !Línea vertical paralela al borde izquierdo
          s(m,n-1)=((-2*(2*(alpha**2)-poisson*(alpha**2)+1))*2.)/(pasox**4) 
          s(m-1,n-1)=(((alpha**2)*(2.-poisson))*2.)/(pasox**4)
          s(m+1,n-1)=(((alpha**2)*(2.-poisson))*2.)/(pasox**4)
   !Línea más exterior
          s(m,n-2)=(1.*2.)/(pasox**4)
	END IF 
!!PUNTO CERCANO A LA ESQUINA SUPERIOR IZQUIEDA (EJE HORIZONTAL)
!cesqizqsuphor:
	IF ((m==1).AND.(n==2)) THEN !CAMBIÉ LOS LÍMITES arreglé el cuadrado en s(m+1,n+-1)
   !Borde horizontal correspondiente al borde superior
          s(m,n)=(((1+4*(alpha**2)*(1-poisson)+2.5*(alpha**4)*(1-(poisson**2)))*2.)/(pasox**4))+D1
          s(m,n-1)=((-2*(alpha**2)*(1-poisson+0.5*(alpha**2)*(1-(poisson**2))))*2.)/(pasox**4)
          s(m,n+2)=((0.5*(alpha**4)*(1-(poisson**2)))*2.)/(pasox**4)
          s(m,n+1)=((-2*(alpha**2)*(1-poisson+(alpha**2)*(1-(poisson**2))))*2.)/(pasox**4)
   !Línea vertical paralela al borde izquierdo
          s(m+1,n)=((-2*(2*(alpha**2)-poisson*(alpha**2)+1))*2.)/(pasox**4) 
          s(m+1,n+1)=(((alpha**2)*(2.-poisson))*2.)/(pasox**4)
          s(m+1,n-1)=(((alpha**2)*(2.-poisson))*2.)/(pasox**4)
   !Línea más exterior
          s(m+2,n)=(1.*2.)/(pasox**4)
	END IF 
!!PUNTO CERCANO A LA ESQUINA SUPERIOR DERECHA (EJE HORIZONTAL)
!cesqdersuphor:
	IF ((m==1).AND.(n==(ppx-1))) THEN !CAMBIÉ LOS LÍMITES ARREGLÉ EL CUADRADO
   !Borde horizontal correspondiente al borde superior
          s(m,n)=(((1+4*(alpha**2)*(1-poisson)+2.5*(alpha**4)*(1-(poisson**2)))*2.)/(pasox**4))+D1
          s(m,n+1)=((-2*(alpha**2)*(1-poisson+0.5*(alpha**2)*(1-(poisson**2))))*2.)/(pasox**4)
          s(m,n-2)=((0.5*(alpha**4)*(1-(poisson**2)))*2.)/(pasox**4)
          s(m,n-1)=((-2*(alpha**2)*(1-poisson+(alpha**2)*(1-(poisson**2))))*2.)/(pasox**4)
   !Línea vertical paralela al borde izquierdo
          s(m+1,n)=((-2*(2*(alpha**2)-poisson*(alpha**2)+1))*2.)/(pasox**4) 
          s(m+1,n+1)=(((alpha**2)*(2.-poisson))*2.)/(pasox**4)
          s(m+1,n-1)=(((alpha**2)*(2.-poisson))*2.)/(pasox**4)
   !Línea más exterior
          s(m+2,n)=(1.*2.)/(pasox**4)
	END IF 
!!PUNTO CERCANO A LA ESQUINA INFERIOR IZQUIERDA (EJE HORIZONTAL)
!cesqizqinfhor:
	IF ((m==(ppy)).AND.(n==2)) THEN !CAMBIÉ LOS LIMITES
   !Borde correspondiente al borde inferior
          s(m,n)=(((1+4*(alpha**2)*(1-poisson)+2.5*(alpha**4)*(1-(poisson**2)))*2.)/(pasox**4))+D1
          s(m,n+1)=((-2*(alpha**2)*(1-poisson+(alpha**2)*(1-(poisson**2))))*2.)/(pasox**4)
          s(m,n+2)=((0.5*(alpha**4)*(1-(poisson**2)))*2.)/(pasox**4)
          s(m,n-1)=((-2*(alpha**2)*(1-poisson+0.5*(alpha**2)*(1-(poisson**2))))*2.)/(pasox**4)
   !Línea  paralela al borde inferior
          s(m-1,n)=((-2*(2*(alpha**2)-poisson*(alpha**2)+1))*2.)/(pasox**4) 
          s(m-1,n-1)=(((alpha**2)*(2.-poisson))*2.)/(pasox**4)
          s(m-1,n+1)=(((alpha**2)*(2.-poisson))*2.)/(pasox**4)
   !Línea más exterior
          s(m-2,n)=(1.*2.)/(pasox**4)
	END IF 
!!PUNTO CERCANO A LA ESQUINA INFERIOR DERECHA (EJE HORIZONTAL)
!cesqderinfhor:
	IF ((m==(ppy)).AND.(n==(ppx-1))) THEN !CAMBIÉ LOS LIMITES Y LOS VOLVÍ A CAMBIAR
   !Borde correspondiente al borde inferior
          s(m,n)=(((1+4*(alpha**2)*(1-poisson)+2.5*(alpha**4)*(1-(poisson**2)))*2.)/(pasox**4))+D1
          s(m,n-1)=((-2*(alpha**2)*(1-poisson+(alpha**2)*(1-(poisson**2))))*2.)/(pasox**4)
          s(m,n-2)=((0.5*(alpha**4)*(1-(poisson**2)))*2.)/(pasox**4)
          s(m,n+1)=((-2*(alpha**2)*(1-poisson+0.5*(alpha**2)*(1-(poisson**2))))*2.)/(pasox**4)
   !Línea  paralela al borde inferior
          s(m-1,n)=((-2*(2*(alpha**2)-poisson*(alpha**2)+1))*2.)/(pasox**4) 
          s(m-1,n-1)=(((alpha**2)*(2.-poisson))*2.)/(pasox**4)
          s(m-1,n+1)=(((alpha**2)*(2.-poisson))*2.)/(pasox**4)
   !Línea más exterior
          s(m-2,n)=(1.*2.)/(pasox**4)
	END IF 
!PUNTOS ESQUINA
!ESQUINA INFERIOR IZQUIERDA
esqinfizq:IF((m==ppy).AND.(n==1)) THEN
   !Línea central
          s(m,n)=(((2*(alpha**2)*(1-poisson)+0.5*(1+(alpha**4))*(1-(poisson**2)))*4.)/(pasox**4))+D1
          s(m,n+1)=((-2*((alpha**2)*(1-poisson)+0.5*(1-(poisson**2))))*4.)/(pasox**4)
          s(m,n+2)=((0.5*(1-(poisson**2)))*4.)/(pasox**4)
   !Línea por encima de la central
          s(m-1,n)=((-2*(alpha**2)*(1-poisson+0.5*(alpha**2)*(1-(poisson**2))))*4.)/(pasox**4) 
          s(m-1,n+1)=((2*(alpha**2)*(1-poisson))*4.)/(pasox**4)
   !Línea superior
          s(m-2,n)=((0.5*(alpha**4)*(1-(poisson**2)))*4.)/(pasox**4)      
END IF esqinfizq
esqinfder: IF((m==ppy).AND.(n==ppx)) THEN
  !Línea central
          s(m,n)=(((2*(alpha**2)*(1-poisson)+0.5*(1+(alpha**4))*(1-(poisson**2)))*4.)/(pasox**4))+D1
          s(m,n-1)=((-2*((alpha**2)*(1-poisson)+0.5*(1-(poisson**2))))*4.)/(pasox**4)
          s(m,n-2)=((0.5*(1-(poisson**2)))*4.)/(pasox**4)
   !Línea por encima de la central
          s(m-1,n)=((-2*(alpha**2)*(1-poisson+0.5*(alpha**2)*(1-(poisson**2))))*4.)/(pasox**4) 
          s(m-1,n-1)=((2*(alpha**2)*(1-poisson))*4.)/(pasox**4)
   !Línea superior
          s(m-2,n)=((0.5*(alpha**4)*(1-(poisson**2)))*4.)/(pasox**4)   
END IF esqinfder
esqsupizq: IF((m==1).AND.(n==1)) THEN
   !Línea central
          s(m,n)=(((2*(alpha**2)*(1-poisson)+0.5*(1+(alpha**4))*(1-(poisson**2)))*4.)/(pasox**4))+D1
          s(m,n+1)=((-2*((alpha**2)*(1-poisson)+0.5*(1-(poisson**2))))*4.)/(pasox**4)
          s(m,n+2)=((0.5*(1-(poisson**2)))*4.)/(pasox**4)
   !Línea por debajo de la central
          s(m+1,n)=((-2*(alpha**2)*(1-poisson+0.5*(alpha**2)*(1-(poisson**2))))*4.)/(pasox**4) 
          s(m+1,n+1)=((2*(alpha**2)*(1-poisson))*4.)/(pasox**4)
   !Línea inferior
          s(m+2,n)=((0.5*(alpha**4)*(1-(poisson**2)))*4.)/(pasox**4)  
END IF esqsupizq 
esqsupder: IF ((m==1).AND.(n==ppx)) THEN
  !Línea central
          s(m,n)=(((2*(alpha**2)*(1-poisson)+0.5*(1+(alpha**4))*(1-(poisson**2)))*4.)/(pasox**4))+D1
          s(m,n-1)=(((-2*((alpha**2)*(1-poisson)+0.5*(1-(poisson**2))))*4.)/(pasox**4))
          s(m,n-2)=((0.5*(1-(poisson**2)))*4.)/(pasox**4)
   !Línea por debajo de la central
          s(m+1,n)=((-2*(alpha**2)*(1-poisson+0.5*(alpha**2)*(1-(poisson**2))))*4.)/(pasox**4)
          s(m+1,n-1)=((2*(alpha**2)*(1-poisson))*4.)/(pasox**4)
   !Línea inferior
          s(m+2,n)=((0.5*(alpha**4)*(1-(poisson**2)))*4.)/(pasox**4) 
END IF esqsupder
!HACEMOS EL VOLCADO A UN VECTOR Y LO ESCRIBIMOS EN UN FICHERO
          t=1
          DO j=1,ppy
		    DO i=1,ppx
		   
		      P(t)=s(j,i)

		      t=t+1
!		     
		    END DO
          END DO

          DO t=1,total
             coef(cont1,t)=P(t)
          END DO
          s=0.0
          w=w+1
 END DO COLUMNAS
END DO FILAS
open(26,FILE="/home/felipe/Escritorio/coef.dat")
do j=1,total
do i=1,total
write(26,*)coef(j,i)
end do
end do
close(26)
P=0.0
X=0.0
write(*,*)'total=',total,'A=',a,'b=',b,'Cargasdistribuidas',cargasdistribuidas
!!Preparamos el vector de cargas
CALL SUBVECTORCARGAS(pasox,a,b,cargasdistribuidas,P,total)
open(26,FILE="/home/felipe/Escritorio/ld.dat")
DO i=1,total
write(26,*)P(i)
end do
write(*,*)'hemos escrito en todos los ficheros'
close(26)
!!Empezamos la resolución
!CALL GAUSSSEIDEL(coef,P,X,total)
CALL getrf(coef,ipiv) 
CALL getrs(coef,ipiv,P,'N') 
open(26,FILE="/home/Placas/soluciones.dat")
DO i=1,total
write(26,*)P(i)
end do
write(*,*)'hemos escrito en todos los ficheros'
close(26)
WRITE(*,*)'CALCULAMOS LOS FLECTORES'
!!Calculamos los momentos
MX=0.0
MY=0.0
MT=0.0
s=0.0
CALL MOMENTOS(P,s,MX,MY,MT,ppx,ppy,poisson,D,pasox)
!CALL MOMENTOS(X,s,MX,MY,MT,ppx,ppy,poisson,D,pasox)
!CLOSE(10)
!CLOSE(11)
!CLOSE(12)
!Limpiamos la memoria para poder reutilizarla
!DEALLOCATE(coef)
!DEALLOCATE(P)
!DEALLOCATE(X)
!DEALLOCATE(s)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
END PROGRAM PROGRAMAPRINCIPAL
!!!!!!!!!!!!!!!!!!!!!!!!!!!FIN DEL PROGRAMA PRINCIPAL!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!FUNCIONES Y SUBRUTINAS!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE PASO(pasodef)
IMPLICIT NONE
INTEGER,PARAMETER :: DOBLE=8
INTEGER:: numerocargas,i,mcdb1,mcdb2,mcdb,mcda1,mcda2,mcda,mcdpilares,DIST,h1,h2
INTEGER:: Euclides
INTEGER:: DISTX,DISTX1,DISTX2,cont,MCDY,h,hant,contportico,cambiopaso,porticos
INTEGER :: hante,he,contporticoe
REAL(KIND=DOBLE),INTENT (INOUT) :: pasodef
REAL(KIND=DOBLE),ALLOCATABLE :: DATOS(:,:)
REAL(KIND=DOBLE) :: ymin,ymax
!REAL,ALLOCATABLE :: DISTX(:) !Creo un vector para las distancias en x
!REAL,ALLOCATABLE :: DISTY(:) !Creo un vector para las distancias en x
OPEN(100,FILE="/home/Placas/cargas.dat")
OPEN(101,FILE="/home/Placas/MisCargasOrdenadas.dat")
READ(100,*) numerocargas
WRITE(*,*) 'NUMERO DE CARGAS=',numerocargas
ALLOCATE(DATOS(numerocargas,5))
!ALLOCATE(DISTX(numerocargas-1))
!ALLOCATE(DISTY(numerocargas-1))
cont=0
contportico=0
contporticoe=0
hante=0
he=0
porticos=0
DISTX=0
DISTX1=0
DISTX2=0
MCDY=0
hant=0
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!VOLCAMOS LOS DATOS DEL FICHERO A LA MATRIZ DATOS
DO i=1,numerocargas
    READ(100,*) DATOS(i,1),DATOS(i,2),DATOS(i,3),DATOS(i,4),DATOS(i,5)
    WRITE(*,12) DATOS(i,:)
    12 FORMAT(5(F10.3))
END DO
!ORDENAMOS EN EL EJE Y
CALL ORDENACIONY(numerocargas,DATOS)
!ORDENAMOS EN EL EJE X
CALL ORDENACIONX(numerocargas,DATOS)
DO i=1,numerocargas
    WRITE(101,*) DATOS(i,:) !ESCRIBIMOS EN EL FICHERO MISCARGASORDENADAS.DAT PARA COMPROBAR
END DO

!YA TENEMOS LA MATRIZ ORDENADA DE MENOR A MAYOR EN X Y EN Y
DO i=1,numerocargas
    IF (i==1) THEN 
       mcdb1=NINT(1000*DATOS(i,3)) !BASE DE PILARES
       mcdb2=NINT(1000*DATOS(i+1,3)) !BASE DE PILARES
       mcda1=NINT(1000*DATOS(i,4)) !ALTO DE PILARES 
       IF (contportico==1) THEN
       			hant=Euclides(h1,h2)
	       		write(*,*)'hant=',hant
	       ELSE
	       		h=Euclides(h1,h2)
	       		write(*,*)'h=',h
	       		hant=Euclides(hant,h)
	       		write(*,*)'hant=',hant
	       END IF
       mcda2=NINT(1000*DATOS(i+1,4)) !ALTO DE PILARES
       mcdb=Euclides(mcdb1,mcdb2) !MCD DE LA BASE DE LOS PILARES 
       mcda=Euclides(mcda1,mcda2) !MCD DE LA BASE DE LOS PILARES
     
    ELSE
       !EN CUALQUIER OTRO CASO, COGER EL ANTERIOR MCD Y EL CORRESPONDIENTE A I
        mcdb1=NINT(1000*DATOS(i,3)) !BASE DE PILARES
        mcdb=Euclides(mcdb1,mcdb)
        mcda1=NINT(1000*DATOS(i,4)) !ALTO DE PILARES
        mcda=Euclides(mcda1,mcda) !MCD DE LA BASE DE LOS PILARES
     
    END IF
    
 AnalisisDatosEjeX: IF ((i<(numerocargas)).AND.(numerocargas>2)) THEN
 CentrosEnLaMismaLineaEjeX: IF (DATOS(i,1)<DATOS(i+1,1)) THEN
  IF (cont==0) THEN !Se halla el mcd de las dos primeras distancias
      DISTX1=NINT(ABS(ABS(DATOS(i,1)-DATOS(i+1,1))-(DATOS(i,3)/2.)-(DATOS(i+1,3)/2.))*1000)
      DISTX2=NINT(ABS(ABS(DATOS(i+1,1)-DATOS(i+2,1))-(DATOS(i+1,3)/2.)-(DATOS(i+2,3)/2.))*1000)
      DISTX=Euclides(DISTX1,DISTX2)
            !WRITE(*,*)'En el if distx=',DISTX	
      cont=1
 ELSE
      IF (i>2) THEN !Empieza con el mcd de los dos primeros y continua con el tercer dato
       DISTX1=NINT(ABS(ABS(DATOS(i,1)-DATOS(i+1,1))-(DATOS(i,3)/2.)-(DATOS(i+1,3)/2.))*1000)
       DISTX=Euclides(DISTX1,DISTX)
           !write(*,*)'En el else distx=',DISTX
      END IF
  END IF
 END IF CentrosEnLaMismaLineaEjeX
 END IF AnalisisDatosEjeX
   
      !SÓLO ENTRA SI ES UNA BANDA INTERIOR
      ymin=DATOS(1,2)
      ymax=DATOS(numerocargas,2)
      IF ((DATOS(i,2)>ymin).AND.(DATOS(i,2)<ymax)) THEN  
      	       contportico=contportico+1
      	       porticos=1
	       CALL DISTANCIASEJEY(DATOS,i,numerocargas,h1,h2,porticos)
	      
	       IF (contportico==1) THEN
	       		hant=Euclides(h1,h2)
	       		write(*,*)'hant=',hant
	       ELSE
	       		h=Euclides(h1,h2)
	       		write(*,*)'h=',h
	       		hant=Euclides(hant,h)
	       		write(*,*)'hant=',hant
	       END IF
      ELSE										!Estoy aquí
                write(*,*)'porticos=',porticos
                write(*,*)'ymin=',ymin,'DATOS(i,2)=',DATOS(i,2)
      		
      	        IF((porticos==0).AND.(DATOS(i,2)==ymin)) THEN
      	        	WRITE(*,*)'ESTOY DENTRO?,tengo una i=',i
      	              	CALL DISTANCIASEJEY(DATOS,i,numerocargas,h1,h2,porticos)
      	              	contporticoe=contporticoe+1
      	              	 IF (contporticoe==1) THEN
				hante=Euclides(h1,h2)
				write(*,*)'hante=',hante
			ELSE
				he=Euclides(h1,h2)
				write(*,*)'he=',he
				hante=Euclides(hante,he)
				write(*,*)'hante=',hante
			END IF	
		END IF
      END IF
END DO
!Una vez salimos ya tenemos el mcd de los pilares 
mcdpilares=Euclides(mcda,mcdb)
write(*,*)'mcdpilares=',mcdpilares
WRITE(*,*)'HANT=',hant,'hante=',hante
hant=Euclides(hante,hant)
write(*,*)'hant=',hant
write(*,*)'DISTX=',DISTX
DIST=Euclides(DISTX,hant)
write(*,*)'DIST=',DIST
pasodef=Euclides(mcdpilares,DIST)
WRITE(*,*)'EL PASO ADECUADO ES ',pasodef/1000.,',metros.'
pasodef=pasodef/1000.
WRITE(*,*)'¿Desea cambiar el paso?.Pulse 1 si desea cambiarlo y cero si no quiere'
READ(*,*)cambiopaso
IF (cambiopaso==1) THEN 
	WRITE(*,*)'Introduzca el paso'
	READ(*,*) pasodef
	WRITE(*,*) 'El nuevo paso es ',pasodef,'metros'
END IF
CLOSE(101)
END SUBROUTINE PASO

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!							FUNCIONES												 !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
FUNCTION Euclides(a,b)
IMPLICIT NONE
INTEGER :: a,b !Variables que entran
INTEGER :: Euclides !Variable que sale
INTEGER :: c
IF (a<b) THEN 
c=b
b=a
a=c
END IF
DO WHILE(b>0)
 c=mod(a,b)
 a=b
 b=c
END DO
Euclides=a
END FUNCTION Euclides
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!							SUBRUTINAS      										 !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

SUBROUTINE ORDENACIONY(NELEM,MATRIZ)
!-----------------------------------------------------
!ORDENACION POR BURBUJA ("bubble sort") de un arreglo
!unidimensional, de menor a mayor.
!NELEM = Número de filas de la matriz
!MATRIZ = Matriz multidimensional a ordenar
!-----------------------------------------------------
IMPLICIT NONE
INTEGER,PARAMETER :: DOBLE=8
INTEGER,INTENT(IN) :: NELEM
INTEGER :: K
REAL(KIND=DOBLE),DIMENSION(NELEM,5),INTENT(INOUT) :: MATRIZ
REAL(KIND=DOBLE),DIMENSION (5) :: AUX
!-----------------------------------------------------
INTEGER I,J

!-----------------------------------------------------
IF (NELEM.LT.2) RETURN
DO I=1,NELEM-1
	DO J=1,NELEM-I
		IF (MATRIZ(J,2).GT.MATRIZ(J+1,2)) THEN
		  DO K=1,5
			AUX(K)= MATRIZ(J,K) !SE VUELCA TODA LA FILA AL VECTOR AUXILIAR
			MATRIZ(J,K)= MATRIZ(J+1,K)
			MATRIZ(J+1,K) = AUX(K)
		  END DO
		END IF
	END DO
END DO
RETURN
END SUBROUTINE ORDENACIONY

SUBROUTINE ORDENACIONX(NELEM,MATRIZ)
!-----------------------------------------------------
!ORDENACION POR BURBUJA ("bubble sort") de un arreglo
!unidimensional, de menor a mayor.
!NELEM = Número de filas de la matriz
!MATRIZ = Matriz multidimensional a ordenar
!-----------------------------------------------------
IMPLICIT NONE
INTEGER,PARAMETER:: DOBLE=8
INTEGER,INTENT(IN) :: NELEM
INTEGER :: K
REAL(KIND=DOBLE),DIMENSION(NELEM,5),INTENT(INOUT) :: MATRIZ
REAL(KIND=DOBLE),DIMENSION (5) :: AUX
!-----------------------------------------------------
INTEGER I,J

!-----------------------------------------------------
IF (NELEM.LT.2) RETURN
DO I=1,NELEM-1
DO J=1,NELEM-I
IF (MATRIZ(J,2)==MATRIZ(J+1,2)) THEN !ORDENAMOS EN EL EJE X
IF (MATRIZ(J,1).GT.MATRIZ(J+1,1)) THEN
DO K=1,5
AUX(K)= MATRIZ(J,K) !SE VUELCA TODA LA FILA AL VECTOR AUXILIAR
MATRIZ(J,K)= MATRIZ(J+1,K)
MATRIZ(J+1,K) = AUX(K)
END DO
END IF
END IF 
END DO
END DO
RETURN
END SUBROUTINE ORDENACIONX
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE DISTANCIASEJEY (CARGAS,J,NELEM,h1,h2,caso)
 IMPLICIT NONE
 INTEGER,PARAMETER :: DOBLE=8
 INTEGER,INTENT(IN) :: J,NELEM,caso
 REAL(KIND=DOBLE),DIMENSION(NELEM,5),INTENT(IN) :: CARGAS
 REAL :: rant,rmin
 INTEGER,INTENT(OUT) :: H1,H2
 INTEGER :: K,BC,YC,YSUP,YINF,BSUP,BINF,CONT=0,M,YMIN,YMAX
 !INTEGER :: Euclides
 IF (caso==1) THEN 
	 K=J
	 WRITE(*,*)'I=',J
	 WRITE(*,*)'CARGAS(J,2)=',CARGAS(J,2)
	  WRITE(*,*)'CARGAS(J,1)=',CARGAS(J,1)
	  WRITE(*,*)'CARGAS(K,2)=',CARGAS(K,2)
	  WRITE(*,*)'CARGAS(K,1)=',CARGAS(K,1)
	  WRITE(*,*)'CARGAS(J,2) ENTERA=',NINT((CARGAS(J,2))*1000)
	  WRITE(*,*)'CARGAS(K,2) ENTERA=',NINT((CARGAS(K,2))*1000)
	 CONT=0
	 DO WHILE(NINT((CARGAS(J,2))*1000)<=NINT((CARGAS(K,2))*1000))
	       K=K-1
	END DO
	write(*,*)'--------------------'
	write(*,*)'k=',K
	write(*,*)'--------------------'
	! DO WHILE (NINT((CARGAS(J,1))*1000).NE.(NINT((CARGAS(K,1))*1000)))
	!      K=K-1
	! END DO
	rant=ABS(CARGAS(J,1)-CARGAS(K,1))
	!rmin=ABS(CARGAS(J,1)-CARGAS(K,1))
	 !WRITE(*,*)'K=',K
	 WRITE(*,*)'ANTES DE ENTRAR EN EL DO-WHILE'
	 WRITE(*,*)'RANT=',rant
	 M=K
	 YMIN=NINT((CARGAS(K,2))*1000)
	 write(*,*)'YMIN=',YMIN
	DO WHILE(YMIN.EQ.NINT((CARGAS(K,2))*1000)) !comparo como entero
	!WRITE(*,*)'DENTRO'
	!WRITE(*,*)'cargas(K,1)=',CARGAS(K,1)
	!WRITE(*,*)'cargas(J,1)=',CARGAS(J,1)

	rmin=ABS(CARGAS(J,1)-CARGAS(K,1))
	IF(rant>=rmin) THEN
	WRITE(*,*)'DENTRO DEL IF M=',K
	      M=K
	rant=rmin
	END IF 
	WRITE(*,*)'rmin=',rmin
	K=K-1
	END DO
	!WRITE(*,*)'CONT=',CONT
	! write(*,*)'valores definitivos de cargas a la salida para inferior'
	! WRITE(*,*)'CARGAS(K,2)=',CARGAS(K,2)
	! WRITE(*,*)'CARGAS(K,1)=',CARGAS(K,1)
	! WRITE(*,*)'K=',K
	! WRITE(*,*)'CARGAS(M,2)=',CARGAS(M,2)
	! WRITE(*,*)'CARGAS(M,1)=',CARGAS(M,1)
	! WRITE(*,*)'M=',M
	 YINF=NINT((CARGAS(M,2))*1000)
	 BINF=NINT((CARGAS(M,4))*1000)
	 K=J
	 WRITE(*,*)'ANTES DEL PRIMER BUCLE CON K=',K
	 DO WHILE(NINT((CARGAS(J,2))*1000).EQ. NINT((CARGAS(K,2))*1000))
	       K=K+1       
	 END DO
	!  DO WHILE (NINT((CARGAS(J,1))*1000).NE.(NINT((CARGAS(K,1))*1000)))
	!      K=K-1
	! END DO

	YMAX=NINT((CARGAS(K,2))*1000)
	M=K
	rant=ABS(CARGAS(J,1)-CARGAS(K,1))
	!rmin=ABS(CARGAS(J,1)-CARGAS(K,1))
	 !WRITE(*,*)'K=',K
	 WRITE(*,*)'ANTES DE ENTRAR EN EL DO-WHILE segundo'
	 WRITE(*,*)'RANT=',rant
	DO WHILE(YMAX.EQ.NINT((CARGAS(K,2))*1000)) !comparo como entero
	!WRITE(*,*)'DENTRO'
	!WRITE(*,*)'cargas(K,1)=',CARGAS(K,1)
	!WRITE(*,*)'cargas(J,1)=',CARGAS(J,1)

	rmin=ABS(CARGAS(J,1)-CARGAS(K,1))
	IF(rant>=rmin) THEN
	WRITE(*,*)'DENTRO DEL IF M=',K
	WRITE(*,*)'CARGAS(K,2)',CARGAS(K,2)
	      M=K
	      rant=rmin
	END IF 
	WRITE(*,*)'rmin=',rmin
	K=K+1
	END DO
	 WRITE(*,*)'CARGAS(M,2)=',CARGAS(M,2)
	 WRITE(*,*)'CARGAS(M,1)=',CARGAS(M,1)
	 
	 WRITE(*,*)'K=',K
	 YSUP=NINT((CARGAS(M,2))*1000)
	 BSUP=NINT((CARGAS(M,4))*1000)
	 YC=NINT((CARGAS(J,2))*1000)
	 BC=NINT((CARGAS(J,4))*1000)
	 H1=ABS(ABS(YSUP-YC)-(BSUP/2)-(BC/2))
	 H2=ABS(ABS(YINF-YC)-(BINF/2)-(BC/2))
	 WRITE(*,*)'YSUP=',YSUP,'BSUP=',BSUP,'YC=',YC,'BC=',BC,'YINF=',YINF,'BINF=',BINF
	 WRITE(*,*)'H1=',H1
	 WRITE(*,*)'H2=',H2
ELSE ! EN EL CASO DE QUE NO EXISTAN PÓRTICOS INTERMEDIOS, NOS INTRODUCIMOS AQUÍ

	 K=J
	 WRITE(*,*)'ANTES DEL PRIMER BUCLE CON K=',K
	 DO WHILE(NINT((CARGAS(J,2))*1000).EQ. NINT((CARGAS(K,2))*1000))
	       K=K+1       
	 END DO
	!  DO WHILE (NINT((CARGAS(J,1))*1000).NE.(NINT((CARGAS(K,1))*1000)))
	!      K=K-1
	! END DO

	YMAX=NINT((CARGAS(K,2))*1000)
	M=K
	rant=ABS(CARGAS(J,1)-CARGAS(K,1))
	!rmin=ABS(CARGAS(J,1)-CARGAS(K,1))
	 !WRITE(*,*)'K=',K
	 WRITE(*,*)'ANTES DE ENTRAR EN EL DO-WHILE segundo'
	 WRITE(*,*)'RANT=',rant
	DO WHILE(YMAX.EQ.NINT((CARGAS(K,2))*1000)) !comparo como entero
	!WRITE(*,*)'DENTRO'
	!WRITE(*,*)'cargas(K,1)=',CARGAS(K,1)
	!WRITE(*,*)'cargas(J,1)=',CARGAS(J,1)

	rmin=ABS(CARGAS(J,1)-CARGAS(K,1))
	IF(rant>=rmin) THEN
	WRITE(*,*)'DENTRO DEL IF M=',K
	WRITE(*,*)'CARGAS(K,2)',CARGAS(K,2)
	      M=K
	      rant=rmin
	END IF 
	WRITE(*,*)'rmin=',rmin
	K=K+1
	END DO
	 WRITE(*,*)'CARGAS(M,2)=',CARGAS(M,2)
	 WRITE(*,*)'CARGAS(M,1)=',CARGAS(M,1)
	 
	 WRITE(*,*)'K=',K
	 YSUP=NINT((CARGAS(M,2))*1000)
	 BSUP=NINT((CARGAS(M,4))*1000)

	 YC=NINT((CARGAS(J,2))*1000)
	 BC=NINT((CARGAS(J,4))*1000)
	 H1=ABS(ABS(YSUP-YC)-(BSUP/2)-(BC/2))
	 H2=ABS(ABS(YSUP-YC)-(BSUP/2)-(BC/2))
	 !MCDY=Euclides(H1,H2)
	 WRITE(*,*)'yc=',YC,'BC=',BC,'YINF=',YINF,'BINF=',BINF
	 WRITE(*,*)'H1=',H1
	 WRITE(*,*)'H2=',H2
END IF
 RETURN
END SUBROUTINE DISTANCIASEJEY

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE GENERACIONMALLADO (a,b,paso,contpunto)
IMPLICIT NONE
INTEGER,PARAMETER :: DOBLE=8
INTEGER  :: i,j,Na,Nb,cargas,k,m,ppilarx,ppilary,n,np,k1,k2,r,moduloa,CONT
INTEGER,INTENT (OUT) :: contpunto
REAL (KIND=DOBLE),INTENT (INOUT) :: a,b
REAL (KIND=DOBLE),INTENT (IN) :: paso !Debe venir calculada desde la subrutina calculo del paso
REAL (KIND=DOBLE) :: xp,yp,Fp
REAL (KIND=DOBLE),DIMENSION(4) :: x,y
REAL (KIND=DOBLE),DIMENSION(:),ALLOCATABLE :: F,xc,yc,ap,bp  
REAL (KIND=DOBLE),DIMENSION(:),ALLOCATABLE:: mallax,mallay
REAL (KIND=DOBLE),DIMENSION(:),ALLOCATABLE:: L

OPEN(10,file="/home/Placas/dimensioneslosa.dat")!el punto inicial debe ser el de menor x y menor y 
OPEN(11,file="/home/Placas/malla.dat",status="REPLACE")
DO i=1,4
READ(10,*) x(i),y(i) !GUARDO EN UN PAR DE VECTORES LAS ESQUINAS
END DO
CLOSE(10)
a=x(4)-x(3)
b=y(3)-y(1)
!NUMERO DE PUNTOS EN LA DIRECCION a y b
Na=NINT(a/paso)+1
Nb=NINT(b/paso)+1
ALLOCATE(mallax(Na))
ALLOCATE(mallay(Nb))
ALLOCATE(L((Na)*(Nb)))
!WRITE(*,*)'El ancho es: ',a, 'El alto es: ',b, 'puntos en la direccion a:',Na,'puntos en la b:',Nb !controlar que funciona el calculo
!READ(*,*)
L=0.0
contpunto=0
yp=0.0
 
DO j=1,Nb!+1 !he cambiado el final, le añadí un uno
	xp=0.0
	DO i=1,Na!He cambiado el final, le añadí un uno
	WRITE(11,103)xp,yp
	103 FORMAT(F10.4,F10.4)
	mallax(i)=xp
	mallay(j)=yp
	xp=xp+paso
	END DO
yp=yp+paso
END DO
!WRITE(*,*) mallax
!WRITE(*,*) '-----------------------------------'
!WRITE(*,*) mallay
CLOSE(11)
!CALL system('gnuplot mallarepre.gnu')
OPEN(12,FILE="/home/Placas/cargas.dat",STATUS="OLD")
OPEN(13,FILE="/home/Placas/cargasdistribuidas.dat",STATUS="REPLACE")
READ(12,*)cargas 
!!Leemos cuantas cargas vamos a repartir en el pilar o área
ALLOCATE(xc(cargas),yc(cargas),ap(cargas),bp(cargas),F(cargas)) 

DO i=1,(cargas)
	    READ(12,*) xc(i),yc(i),ap(i),bp(i),F(i)
	    WRITE(*,*) xc(i),yc(i),ap(i),bp(i),F(i)
	!debemos buscar el punto de la malla que esté entre los puntos del pilar

	    k=1
	    m=1
	DO WHILE (xc(i)>mallax(m)) !Buscamos la cruceta del pilar o área
	  DO WHILE (yc(i)>mallay(k)) !Buscamos la cruceta del pilar o área
	     k=k+1
	  END DO 
	m=m+1
	END DO
	WRITE(*,*)'xc=',xc(i),'yc=',yc(i)
	!Aquí tendré la posición en el vector de k y la posición en el vector de la m
	!con lo cual tengo la posicion aproximada en mi malla
	!WRITE(*,*) 'EL PUNTO ES,X=',mallax(m),'Y=',mallay(k)!comprobación visual
	ppilarx=NINT(ap(i)/paso)+1!cantidad de puntos en la direccion x del pilar	
	ppilary=NINT(bp(i)/paso)+1!cantidad de puntos en la direccion y del pilar
	np=(ppilarx*ppilary) !Número total de puntos en el pilar
	Fp=F(i)/np !Fuerza por punto de estudio del pilar o del área correspondiente
	!EL PROBLEMA ESTÁ CUANDO EL PILAR ES DE UN TAMAÑO IMPAR la primera cifra, P.EJ 0.3
	!WRITE(*,*)'ncarga=',i,'ap=',ap(i),'bp=',bp(i),'ppilarx=',ppilarx,'ppilary=',ppilary,'np=',np,'F=',F(i),'FP=',Fp
	moduloa=NINT(ap(i)*10)
	WRITE(*,*)'MODULOA=',mod(moduloa,2)
	IF(mod(ppilarx,2)==0) THEN !esto va bien si el pilar es cuadrado, hay que meterse para el caso de que no lo sea
	k1=m-((ppilarx)/2)!puntos extremos esquina inferior izq del pilar considerado
	k2=k-((ppilary)/2)
	!    if (NINT(ap(i)*10)==4) then
	!        write(*,*) 'k1=',k1,'k2=',k2,'ppilarx',ppilarx
	!    end if
	!WRITE(*,*) 'ESTOY EN EL IF'
	ELSE
	k1=m-(ppilarx/2)!puntos extremos esquina inferior izq del pilar considerado
	k2=k-(ppilary/2)
	!WRITE(*,*) 'ESTOY EN EL ELSE'
	END IF
	!write(*,*) 'm=',m,'k=',k,'k1=',k1,'k2=',k2,'ppilarx=',ppilarx,'ppilary',ppilary

	CONT=0
	IF (mod(ppilarx,2)==0) THEN !REVISAR TRUCO QUE ME SAQUÉ DE LA MANGA

		DO r=k2,k2+ppilary-1
			DO n=k1,k1+ppilarx-1
			    CONT=CONT+1
			    contpunto=contpunto+1
			    WRITE(13,100) mallax(n),mallay(r),Fp
			    100 FORMAT(F10.2,F10.2,F10.2) 
			    !WRITE(*,*)'R=',r,'N=',n,'HEMOS ESTUDIADO TANTOS PUNTOS:',CONT 
			END DO
		END DO
	ELSE
	CONT=0
		DO r=k2,(k2+ppilary-1)
			DO n=k1,(k1+ppilarx-1)
			    CONT=CONT+1
			    contpunto=contpunto+1
			    !write(*,*)'n=',n,'r=',r
			    WRITE(13,101) mallax(n),mallay(r),Fp
			    101 FORMAT(F10.2,F10.2,F10.2) 
			END DO
		END DO
	END IF

END DO

CLOSE(12)
CLOSE(13)

call system('gnuplot cargasdistribuidas.gnu')

END SUBROUTINE GENERACIONMALLADO
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!                            aqui
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE SUBVECTORCARGAS (paso,a,b,cargas,L,total)
IMPLICIT NONE
INTEGER,PARAMETER :: DOBLE=8
INTEGER :: Na,Nb,i,k,MICONTADOR=0
INTEGER,INTENT(IN) :: cargas,total
REAL (KIND=DOBLE),INTENT(IN) :: paso,a,b
REAL (KIND=DOBLE) :: xp,yp,xf,yf,F
REAL (KIND=DOBLE),DIMENSION(total),INTENT(INOUT) :: L
OPEN(15,FILE="/home/Placas/cargasdistribuidas.dat")
OPEN(16,FILE="/home/Placas/vectorcargas.dat")
!paso=0.1 !vendrá dado por la generación
!a=30.!30.
!b=20.!20.
!Na=NINT((a/paso))+2 !he cambiado esto
!Nb=NINT((b/paso))+2
write(*,*)'los datos que se pasan a subvectorcargas son:'
WRITE(*,*)'PASO=',paso,'a=',a,'b=',b,'cargas=',cargas,'total=',total
write(*,*)'Na=',Na,'Nb=',Nb,'total=',total

!ALLOCATE(L(Na*Nb))
!cargas=98 !Debemos hacer la contabilización en la parte de generación de cargas distribuidas
!genero la malla
!write(*,*)'el paso es:',paso
!read(*,*)

DO i=1,cargas
	  WRITE(*,*) 'CARGA=',i
	  READ(15,*) xf,yf,F
	  WRITE(*,*) xf,yf,F
	!  write(*,*)(NINT(xp*100))
	!  write (*,*) (NINT(xf*100))
	xp=0.0
	yp=0.0
	k=1
	!WRITE(*,*)'TENEMOS UNA XF=',xf,'TENEMOS UNA YF=',yf

	DO WHILE (((NINT(xp*100)).NE.(NINT(xf*100))).OR.((NINT(yp*100)).NE.(NINT(yf*100))))
	!Recorremos en busca del punto de aplicación de la carga
	    IF (xp>a) THEN
		xp=0.0
		yp=yp+paso
		!k=k-1 !he cambiado esto
	    ELSE
		xp=xp+paso
		k=k+1
	    END IF
	    
	 
	 END DO
	  WRITE(*,*)'K=',k,'F=',F,'xp=',xp,'yp=',yp,'xf=',xf,'yf=',yf
	  WRITE(*,*)'---------------------------------------------------------------------------'
	  L(k)=L(k)+F
END DO
DO i=1,(total)
   WRITE(16,*) L(i)
   write(*,*) L(i)
   IF (L(i)>0.0) THEN
         MICONTADOR=MICONTADOR+1
   END IF
END DO
WRITE(*,*)'MICONTADOR=',MICONTADOR
CLOSE(15)
CLOSE(16)
END SUBROUTINE SUBVECTORCARGAS
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE GAUSSSEIDEL (A,B,X,N)
IMPLICIT NONE
INTEGER,PARAMETER :: DOBLE=8,LONG=8
INTEGER ,INTENT(IN) :: N
REAL (KIND=DOBLE),DIMENSION(N,N),INTENT(INOUT) :: A
REAL (KIND=DOBLE),DIMENSION(N),INTENT(INOUT) ::  B
REAL (KIND=DOBLE),DIMENSION(N),INTENT(INOUT) :: X 
INTEGER :: MAXITER,CENTINELA,ITER
INTEGER  :: I,J
REAL (KIND=DOBLE) :: TOL,LAMBDA,AII,XANT,SUMA,ERROR
REAL :: inicio, final !CÁLCULO DE TIEMPO DE COMPUTACIÓN
WRITE(*,*)'COMIENZA EL CALCULO MEDIANTE GAUSS SEIDEL'
MAXITER=10000000
LAMBDA=1.0 !PARA VALOR UNIDAD SE TRADUCE EN GAUSS SEIDEL SIN TÉCNICAS DE RELAJACIÓN
X=0.0
TOL=0.0001!Tolerancia en porcentaje
!WRITE(*,*)'EL TERMINO A(1,1) ES=)',A(1,1)
call cpu_time(inicio)
DO I=1,N
     !WRITE(*,*)'HOLA ESTOY EN GAUSS SEIDEL,I=',I
      
     AII=A(I,I)
     !WRITE(*,*)'HOLA ESTOY EN GAUSS SEIDEL,AII=',AII
     DO J=1,N
       A(I,J)=A(I,J)/AII
     END DO
     B(I)=B(I)/AII
END DO
DO I=1,N
SUMA=B(I)
     DO J=1,N
         IF (I.NE.J) THEN
             SUMA=SUMA-A(I,J)*X(J)
         END IF
     END DO
     X(I)=SUMA
END DO
ITER=1
CENTINELA=1!VARIABLE QUE USAMOS PARA SALIR DEL BUCLE EN CASO DE QUE SE CUMPLA LA TOLERANCIA
DO WHILE ((ITER<MAXITER).AND.(CENTINELA==1))!MIRAR SI SE PUEDE HACER ASI

 DO I=1,N
     XANT=X(I)
     SUMA=B(I)
     DO J=1,N
         IF (I.NE.J) THEN
             SUMA=SUMA-A(I,J)*X(J)
         END IF
     END DO
     X(I)=LAMBDA*SUMA+(1.-LAMBDA)*XANT
     IF ((CENTINELA==1).AND.(X(I).NE.0)) THEN
         ERROR=ABS((X(I)-XANT)/X(I))*100.
         IF (ERROR<TOL) THEN 
             CENTINELA=0 
         END IF
     END IF
END DO
ITER=ITER+1
END DO
WRITE (*,*) 'ITERACION FINAL=',ITER
call cpu_time(final)
print '("Tiempo de calculo del sistema es  ",f16.3," segundos.")',final-inicio
IF (CENTINELA==1) THEN
WRITE(*,*)'EL ALGORITMO NO CONVERGE'
ELSE
!ABRIMOS EL ARCHIVO Y POSTERIORMENTE GUARDAMOS LA SOLUCIÓN
OPEN(19,FILE="/home/Placas/solucioneslosa.dat")
DO I=1,N
     WRITE(19,*)X(I)
END DO
END IF
END SUBROUTINE  GAUSSSEIDEL
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE MOMENTOS(SOL,MSOL,MOMX,MOMY,MOMT,PPX,PPY,POISSON,D,PASOX)
IMPLICIT NONE
INTEGER,PARAMETER :: DOBLE=8
INTEGER :: I,J,K
INTEGER,INTENT(IN) :: PPX,PPY
REAL(KIND=DOBLE),INTENT(IN) :: POISSON,D,PASOX
REAL(KIND=DOBLE), INTENT (INOUT),DIMENSION(PPX*PPY) :: SOL
REAL(KIND=DOBLE), INTENT (INOUT),DIMENSION(PPX,PPY) :: MSOL
REAL(KIND=DOBLE), INTENT (INOUT),DIMENSION(PPX,PPY) :: MOMX,MOMY,MOMT
!Abrimos los ficheros donde guardaremos los flectores en cada caso
OPEN(120,FILE="/home/Placas/MomentosEjeX.dat")
OPEN(121,FILE="/home/Placas/MomentosEjeY.dat")
OPEN(122,FILE="/home/Placas/MomentosTotales.dat")
!Comienza el volcado del vector soluciones a una matriz, que represente la losa
K=1
DO I=1,PPY
	DO J=1,PPX
		MSOL(I,J)=SOL(K)
		K=K+1
	END DO
END DO
!Finaliza el volcado. Ahora estudiamos para cada punto, cuales son sus momentos
!Inicialmente sólo lo haremos para borde libre, donde las condiciones de contorno
!son nulas: Mx=0,My=0;Mxy
WRITE(*,*)'W(2)=',MSOL(1,2),'W(24)=',MSOL(3,2),'W(12)=',MSOL(2,1),'W(14)=',MSOL(2,3)
I=2
J=2
WRITE(*,*)'MOMENTO X DEL PUNTO 13=',(2*(1-POISSON)*MSOL(I,J)-1.*MSOL(I,J-1)-1.*MSOL(I,J+1)-POISSON*MSOL(I-1,J)-POISSON*MSOL(I+1,J))*(D/(PASOX**2))
WRITE(*,*)'MOMENTO Y DEL PUNTO 13=',(2*(1-POISSON)*MSOL(I,J)-1.*MSOL(I-1,J)-1.*MSOL(I+1,J)-POISSON*MSOL(I,J-1)-POISSON*MSOL(I,J+1))*(D/(PASOX**2))
READ(*,*)
DO I=2,(PPY-1)
	DO J=2,(PPX-1)
		!Aquí debemos estudiar la molécula para los momentos
		MOMX(I,J)=(2*(1-POISSON)*MSOL(I,J)-1.*MSOL(I,J-1)-1.*MSOL(I,J+1)-POISSON*MSOL(I-1,J)-POISSON*MSOL(I+1,J))*(D/(PASOX**2))
		MOMY(I,J)=(2*(1-POISSON)*MSOL(I,J)-1.*MSOL(I-1,J)-1.*MSOL(I+1,J)-POISSON*MSOL(I,J-1)-POISSON*MSOL(I,J+1))*(D/(PASOX**2))
		MOMT(I,J)=(MOMX(I,J)+MOMY(I,J))/(1-POISSON)
		
	END DO
END DO

DO I=1,PPY
	DO J=1,PPX
	!Escribimos en los ficheros
		WRITE(120,*) MOMX(I,J)
		WRITE(*,*)'MX=',MOMX(I,J),'MY=',MOMY(I,J),'MT=',MOMT(I,J)
		WRITE(121,*) MOMY(I,J)
		WRITE(122,*) MOMT(I,J)
	END DO
END DO
CLOSE(120)
CLOSE(121)
CLOSE(122)
END SUBROUTINE MOMENTOS
!ifort -openmp -O3 -o PROGRAMA30102010 -I$MKLINCLUDE -L$MKLPATH PROGRAMA30102010.f90 -lmkl_lapack95 -Wl,--start-group $MKLPATH/libmkl_intel.a $MKLPATH/libmkl_intel_thread.a $MKLPATH/libmkl_core.a -Wl,--end-group -liomp5 -lpthread
