C ------------------------------------------------------------------------
C    GENERAZIONE AUTOMATICA DI UNA MESH RETTANGOLARE
C
C    RETTANGOLO DI LATI  F1 (ASSE X),  F2(ASSE Y), CON UN VERTICE
C    IN X_0, Y_0, SITUATO NEL PRIMO QUADRANTE
C
C  DATI DI INPUT:
C    F1, F2  : lati del rettangolo lungo l'asse x e y
C    M1      : numero suddivisioni lato F1
C    M2      : numero suddivisioni lato F2
C    X_0, Y_0: coordinate  vertice del rettangolo
C
C    ------------------------------
C   |                              |
C   |                              |
C   |                              |
C   |                              |
C   |                              |  F2
C   |                              |
C   |                              |
C   |                              |
C   |                              |
C   |                              |
C    ------------------------------
C   X_0,Y_0        F1
C
C  DATI DI OUTPUT
C
C  NN    : numero nodi della triangolazione
C  NT    : numero triangoli della griglia
C  X(NN), Y(NN): coordinate dei nodi
C
C
C
C ----------------------------------------------------------------------------
C     
      program griglia
      implicit none
      integer  M1,M2, E(3),A(3)
      integer  N1,N2,NN,NT,I,J,I7, IA1
      real*8 f1,f2
      real*8 X(1000000),Y(1000000)
      real*8 X_0, Y_0
      real*8 z1,z2,s1,s2
C  file di input: mesh.data
C  questo file di input deve essere un file che contiene riga per riga
C  i seguenti valori:
C  F1, F2
C  M1, M2
C  X_0, Y_0
      OPEN(5,FILE='mesh.data',STATUS='OLD')
C  files di output: 
C  mesh.out: in questo file vengono scritti i dati di input del problema
C            e i valori di output N e NT 
C  meshxytria.out:  file di output che contiene i valori N e NT,
C                   la successione dei nodi in senso antiorario per
C                   ciascun triangolo e le coordinate dei nodi
      OPEN(6,FILE='mesh.out',STATUS='UNKNOWN')
      OPEN(13,FILE='meshxytria.out',STATUS='UNKNOWN')
      READ(5,*)F1,F2
      WRITE(6,22)F1,F2
      READ(5,*)M1,M2
      WRITE(6,255)M1,M2
      READ(5,*) X_0, Y_0
      WRITE(6,256) X_0,  Y_0
      N1=M1+1
      N2=M2+1
      NN=N1*N2
      NT=2*M1*M2
      WRITE(6,24)NN,NT
      WRITE(13,*)NN,NT
C
C   CALCOLO DEI VERTICI DEI TRIANGOLI
C      PRIMO CICLO
      DO  I=1,M2
        E(3)=I
        DO J=1,M1
           E(1)=E(3)
           E(2)=E(1)+N2+1
           E(3)=E(2)-1
           WRITE(13,125)(E(I7),I7=1,3)
        END DO
        A(3)=I+1
        DO J=1,M1
           A(2)=A(3)
           A(1)=A(3)-1
           A(3)=A(2)+N2
           WRITE(13,125)(A(I7),I7=1,3)
        END DO  
      END DO

C
C
C   CALCOLO DELLE COORDINATE DEI NODI
C
      Z1=F1/DFLOAT(M1)
      Z2=F2/DFLOAT(M2)
      IA1=0
      DO  J=1,N1
          DO  I=1,N2
            IA1=IA1+1
            S1=DFLOAT(I)
            S2=DFLOAT(J)
            X(IA1)= X_0 + (S2-1.0D0)*Z1
            IF (DABS(X(IA1)) .LT. 1.0D-13) X(IA1)=0.0D+00
            Y(IA1)= Y_0+ F2-(S1-1.0D0)*Z2
            IF (DABS(Y(IA1)) .LT. 1.0D-13) Y(IA1)=0.0D+00
            WRITE(13,47)X(IA1),Y(IA1)
          END DO
      END DO
 22     FORMAT(' LATO F1(ASSEX)=',E12.5,3X,  'LATO F2 (ASSE Y) =',E12.5)
 23     FORMAT(2I3)
 255    FORMAT(' N.SUDD. LATO F1=',I12,3X,'N. SUDD. LATO F2=',I12)
 256   FORMAT('              X_0=', E12.5,3x,'             Y_0=', E12.5)
 24    FORMAT('           N.NODI=',I12,3X,'          N.EL. FIN.=',I12)
 47     FORMAT(2X,2(3X,1PE15.6))
 125  FORMAT(3I15)
      CLOSE(5)
      CLOSE(6)
      CLOSE(13)
      CLOSE(15)
      STOP
      END
