Badly formatted Ada code from science paper

badly_formatted.adb

with Text_IO,Ada.Numerics.Generic_Elementary_Functions;
with Ada.Numerics;
use Text_IO,Ada.Numerics;

procedure Eutecticmain is
   package Real_Io is new FLOAT_IO(Float);
   use Real_Io;
   package Integer_IO is new Text_IO.INTEGER_IO(Integer); use Integer_IO;
   use Integer_Io;
   package Elem_Fct is new Ada.Numerics.Generic_Elementary_Functions (Float);
   use Elem_Fct;

   Relax: constant := 0.1; -- Under - relaxation
   Rg: constant := 8.314510; -- R
   MaxConst : constant := 20; -- Max number of components
   Precision : constant := 0.00001; -- Accuracy
   To : constant := 273.15; -- Temperature reference
   N_Const: Integer;
   CaractH : array (1 .. MaxConst ) of Float;
   CaractT : array (1 .. MaxConst ) of Float;
   I,J : Integer;
   Jc : array (1 .. MaxConst,1 .. MaxConst ) of Float ;
   Dx,X,F : array (1 .. MaxConst ) of Float ;
   Coef,G,Residu,Tp : Float ;

   -- Jc jacobian matrix inversion then multiply by F vector
   procedure SolvJc
   is
      I,J,K,L,N:Integer;
      X: Float ;
   begin
      for I in 1 ..  N_Const loop
         begin
            X:=Jc(I,I);
            Jc(I,I):=1.0;
            for J in 1 ..  N_Const loop
               Jc(I,J):=Jc(I,J)/X;
            end loop;
            for K in 1 ..  N_Const loop
               if K/=I then
                  begin
                     X:=Jc(K,I);
                     Jc(K,I):=0.0;
                     for L in 1 ..  N_Const loop
                        Jc(K,L):=Jc(K,L)-X*Jc(I,L);
                     end loop;
                  end;
               end if;
            end loop;
         end;
      end loop;
      for I in 1 ..  N_Const loop
         begin
            Dx(I):=0.0;
            for J in 1 ..  N_Const loop
               Dx(I):=Dx(I)+Jc(J,I)*F(J);
            end loop;
         end;
      end loop;
   end;
begin
   Put(".....................................................................");New_Line;
   Put(".....................................................................");New_Line;
   Put(".....................................................................");New_Line;
   Put(" ");New_Line;
   Put(" N Components Eutectic ");New_Line;
   Put(" ");New_Line;
   Put(" L.Brunet 2002 ");New_Line;
   Put(".....................................................................");New_Line;
   Put(".....................................................................");New_Line;
   Put(".....................................................................");New_Line;
   Put("N components= ");
   Get(N_Const);
   for I in 1.. N_Const loop
      begin
         Put("Enter Hfus and Tfus for ");Put(I); New_Line;
         Put(" H(kJ/mol)=");
         Get(CaractH(I));
         Put(" T(K)=");
         Get(CaractT(I));
      end;
   end loop;
   Put("Running...");New_Line;
   -- Create jacobian matrix
   for I in 1.. N_Const loop
      begin
         for J in 1.. N_Const loop
            Jc(I,J):=0.0;
         end loop;
         Dx(I):=0.0;
         F(I):=0.0;
         X(I):=1.0/Float(N_Const);
      end;
   end loop;
   X(N_Const):=CaractT(1);
   loop
      for I in 1.. N_Const loop
         begin
            if I<N_Const then
               Jc(I,I):=1.0/X(I)
                 ;
            else
               begin
                  Coef:=1.0;
                  for J in 1.. N_Const-1 loop
                     Coef:=Coef-X(J);
                  end loop;
                  for J in 1.. N_Const-1 loop
                     Jc(J,I):=-1.0/Coef;
                  end loop;
               end;
            end if;
            Jc(N_Const,I):=-CaractH(I)/Rg/(X(N_Const)**2);
         end;
      end loop;
      -- Create vector F
      for I in 1.. N_Const loop
         begin
            if I<N_Const then
               G:=LOG(X(I))+
                 CaractH(I)/Rg*(1.0/X(N_Const)-1.0/CaractT(I));
            else
               G:=LOG(Coef)+
                 CaractH(I)/Rg*(1.0/X(N_Const)-1.0/CaractT(I));
            end if;
            F(I):=G;
         end;
      end loop;
      -- Compute DeltaX
      SolvJc;
      -- Newton-Raphson method
      for I in 1.. N_Const loop
         X(I):=X(I)-Relax*Dx(I);
      end loop;
      -- Calcul du residu
      Residu:=0.0;
      for I in 1.. N_Const-1 loop
         Residu:=Residu+abs(Dx(I));
      end loop;
      exit when Residu<Precision;
   end loop;
   -- Display results
   Put("------RESULT------");
   New_Line;
   Coef:=1.0;
   for J in 1.. N_Const-1 loop
      Coef:=Coef-X(J);
      Put("X"); Put(J);Put("=");Put(X(J));New_Line;
   end loop;
   Put("X");Put(N_Const);Put("=");Put(Coef);New_Line;
   Put("T(K)=");Put(X(N_Const));Put(" = ");Put(X(N_Const)-To,5,2);Put(" Degrees C");
end Eutecticmain;