I am developing a Fortran code (.f90) which shall calculate some matrix in some time interval (dt1=0.001) and these matrices have to be integrated in some time steps (dt=0.1). Though I am experienced in FORTRAN 77, new to Fortran 90. I am unable to make dimension of matrix real (I think that is the problem, I may be wrong!). Below is the part of a long program. I am trying (from last 1 month) different methods but no success, output are NaN. I am using gfortran in Ubuntu 16.04
PROGRAM MODEL
IMPLICIT NONE
REAL::DT,DT1,DK2
INTEGER::I,J,IL,IM,E(100000),jm,DK1
REAL::XN2(1),XO2(1)
REAL,DIMENSION(100000)::AN,BN,AN1,BN1
REAL,DIMENSION(21)::XA,XB
REAL,DIMENSION(21,21)::WA,WB,WAB,WBA
REAL,DIMENSION(21)::SUMA,SUMAB,SUMB,SUMBA
XN2=1E-7
XO2=1E-8
OPEN(1,FILE='TEST1.dat')
OPEN(2,FILE='TEST2.dat')
OPEN(3,FILE='TEST3.dat')
OPEN(4,FILE='TEST4.dat')
DO I=1,21
read(1,*)(WA(I,J),J=1,21)
read(2,*)(WB(I,J),J=1,21)
read(3,*)(WAB(I,J),J=1,21)
read(4,*)(WBA(I,J),J=1,21)
!enddo
IM=1
IL=1
AN(IM)=0.0
BN(IM)=0.0
!AN1(IL)=0.0
!BN1(IL)=0.0
JM=1
E(:)=0.0
DT=0.1
DT1=0.01
DK1=1.0
DK2=1.0
DO WHILE(DK1<=10)
!DO I=1,21
DO WHILE(DK2<=100)
do j=1,21
CALL XAA(WA,WAB,SUMA,SUMAB,XN2,XO2,AN(IM),BN(IM),XA)
CALL XBB(WB,WBA,SUMB,SUMBA,XN2,XO2,AN(IM),BN(IM),XB)
enddo
AN(IM+1)=AN(IM)+XA(I)*DT1
BN(IM+1)=BN(IM)+Xb(I)*DT1
AN1(IL)=AN(IM+1)
BN1(IL)=BN(IM+1)
AN(IM)=AN1(IL)
BN(IM)=BN1(IL)
IM=IM+1
IL=IL+1
WRITE(1112,203)I,an(Im),bn(Im)
DK2=DK2+DT1
ENDDO
203 FORMAT(' ',i2,1000000E13.5)
E(JM)=DK1
AN(DK1)=AN(IM)
BN(DK1)=BN(IM)
WRITE(1111,203)I,AN(DK1),BN(DK1)
DK1=DK1+DT
JM=JM+1
ENDDO
ENDDO
END PROGRAM MODEL
!!SUBROUTINES XAA
SUBROUTINE XAA(WA,WAB,SUMA,SUMAB,XN2,XO2,AN,BN,XA)
INTEGER::I,J
REAL,DIMENSION(100000)::AN,BN
REAL,intent(out)::XA(21)
REAL::XN2(1),XO2(1)
REAL,DIMENSION(21,21)::WA,WAB
REAL,DIMENSION(21)::SUMA,SUMAB
SUMA(:)=0.0
SUMAB(:)=0.0
DO I=1,21
DO J=1,21
SUMA(I)=SUMA(I)+WA(I,J)*XN2(1)
SUMAB(I)=SUMAB(I)+WAB(I,J)*XO2(1)
ENDDO
XA(I)=1e-5+SUMA(I)*AN(1)+SUMAB(I)*BN(1)
ENDDO
RETURN
END SUBROUTINE XAA
!!SUBROUTINES XBB
SUBROUTINE XBB(WB,WBA,SUMB,SUMBA,XN2,XO2,AN,BN,XB)
INTEGER::I,J
REAL::XN2(1),XO2(1)
REAL,DIMENSION(100000)::AN,BN
REAL,intent(out)::XB(21)
REAL,DIMENSION(21,21)::WB,WBA
REAL,DIMENSION(21)::SUMB,SUMBA
SUMB(:)=0.0
SUMBA(:)=0.0
DO I=1,21
DO J=1,21
SUMB(I)=SUMB(I)+WB(I,J)*XN2(1)
SUMBA(I)=SUMBA(I)+WBA(I,J)*XO2(1)
ENDDO
XB(I)=1e-5+SUMB(I)*AN(1)+SUMBA(I)*BN(1)
ENDDO
RETURN
END SUBROUTINE XBB
Further to previous comments in the chat and changes to the code:
For such a uniformly discretized problem, it is very easy (and best) to use integers for all the manipulations of array indices, and recover real values of time by multiplying numbers of steps and time discretization intervals.
In this case, I think it is easiest if you use two nested loops: an external loop over "long" integrated time steps, and an internal one where you perform the integration over "short" time steps. At the beginning of each long time step you copy the value of your vectors from the previous long time step, and for each short time step in the inner loop you add what you need to add to them.
In your case, you can use array syntax to avoid so many row and column indices, both in the main program and in the subroutines.
Please do read all the coefficients of your input matrices at the beginning of your program, before starting to use their values.
If you are going to integrate real values, you are likely going to add or subtract real numbers with different orders of magnitude. In order to do this accurately, you should ensure a high enough precision in your reals, e.g., using a KIND obtained via the Fortran 95 SELECTED_REAL_KIND function.
Provide explicit interfaces for your subroutines, e.g., by making them internal. Then you can pass arrays as implicit-shape, which simplifies the declarations.
If an intent can be defined for each of the arguments of your subprograms, please provide them. If the compiler has access to intents (via explicit interfaces) it can do more checks and more optimisations.
The usual Fortran 77 practice of passing workspaces to subroutines is no longer needed, and for small workspaces it should not affect performance in any noticeable way.
It is good practice to define parameters only once, rather than having code filled with 21s.
The resulting code is:
PROGRAM MODEL
IMPLICIT NONE
INTEGER, PARAMETER :: DP = SELECTED_REAL_KIND(15,300)
INTEGER, PARAMETER :: N = 21 ! Dimension of arrays.
INTEGER, PARAMETER :: Nlong = 1000 ! Number of integrated ("long")
! time steps.
INTEGER, PARAMETER :: Nshort = 100 ! Number of short time steps
! within each long time step.
REAL(kind=DP), PARAMETER :: DT1 = 0.001_DP ! Small time step
REAL(kind=DP), PARAMETER :: DT = Nshort * DT1 ! Long time step
! REAL(kind=DP), PARAMETER :: T = DT * Nlong ! Total time
REAL(kind=DP), DIMENSION(N,N)::WA,WB,WAB,WBA ! Input arrays.
REAL(kind=DP) :: XN2(1), XO2(1) ! Other coefficients.
REAL(kind=DP), DIMENSION(N,0:Nlong) :: AN, BN ! Output vectors.
REAL(kind=DP), DIMENSION(N) :: XA, XB ! Vectors used in integrations.
INTEGER :: I, J ! Indices for rows and columns of arrays.
INTEGER :: IM ! Index for long time steps.
INTEGER :: IL ! Index for short time steps.
REAL(kind=DP) :: time_now
! Set coefficient values.
XN2=1E-7_DP
XO2=1E-8_DP
! Open input files,
OPEN(1,FILE='TEST1.dat')
OPEN(2,FILE='TEST2.dat')
OPEN(3,FILE='TEST3.dat')
OPEN(4,FILE='TEST4.dat')
! read input arrays
DO I=1,N
read(1,*)(WA(I,J),J=1,N)
read(2,*)(WB(I,J),J=1,N)
read(3,*)(WAB(I,J),J=1,N)
read(4,*)(WBA(I,J),J=1,N)
END DO
! and close input files.
DO I = 1, 4
CLOSE(I)
END DO
! Initialize output vectors for time=0.
AN(:,0) = 0.0_DP
BN(:,0) = 0.0_DP
! Loop over long time steps.
DO IM = 1, Nlong
! Initialize vectors from final values at previous time step.
AN(:,IM) = BN(:,IM-1)
BN(:,IM) = BN(:,IM-1)
! Loop over short time steps.
DO IL = 1, Nshort
time_now = (IM-1) * DT + IL * DT1
CALL XAA(WA,WAB,XN2,XO2,AN(:,IM),BN(:,IM),XA)
CALL XBB(WB,WBA,XN2,XO2,AN(:,IM),BN(:,IM),XB)
! Update vectors.
AN(:,IM) = AN(:,IM) + XA * DT1
BN(:,IM) = BN(:,IM) + XB * DT1
! Print their values.
WRITE(1112,'(F12.5,2(1X,i7),42(1X,E13.5))') &
time_now, IM, IL, AN(:,IM), BN(:,IM)
END DO
! Print integrated values.
WRITE(1111,'(F12.5,1X,i7,42(1X,E13.5))') &
time_now, IM, AN(:,IM), BN(:,IM)
END DO
CONTAINS
!!SUBROUTINES XAA
SUBROUTINE XAA(WA,WAB,XN2,XO2,AN,BN,XA)
! Arguments
REAL(kind=DP), DIMENSION(:,:) :: WA,WAB
REAL(kind=DP), DIMENSION(:), INTENT(IN) :: AN,BN
REAL(kind=DP), INTENT(OUT) :: XA(:)
REAL(kind=DP), INTENT(IN) :: XN2(1),XO2(1)
! Local variables.
REAL(kind=DP), DIMENSION(size(XA)) :: SUMA,SUMAB
INTEGER::I
DO I=1,N
SUMA(I)=SUM(WA(I,:))*XN2(1)
SUMAB(I)=SUM(WAB(I,:))*XO2(1)
END DO
XA(:)=1e-5_DP+SUMA(:)*AN(:)+SUMAB(:)*BN(:)
END SUBROUTINE XAA
!!SUBROUTINES XBB
SUBROUTINE XBB(WB,WBA,XN2,XO2,AN,BN,XB)
! Arguments
REAL(kind=DP), DIMENSION(:,:), INTENT(IN) :: WB,WBA
REAL(kind=DP), INTENT(IN) ::XN2(1),XO2(1)
REAL(kind=DP), DIMENSION(:), INTENT(IN) ::AN,BN
REAL(kind=DP), INTENT(OUT) :: XB(:)
! Local variables.
INTEGER::I
REAL(kind=DP), DIMENSION(size(XB)) :: SUMB,SUMBA
DO I=1,N
SUMB(I)=SUM(WB(I,:))*XN2(1)
SUMBA(I)=SUM(WBA(I,:))*XO2(1)
ENDDO
XB(:)=1e-5_DP+SUMB(:)*AN(:)+SUMBA(:)*BN(:)
END SUBROUTINE XBB
END PROGRAM MODEL