Hi All,

I am trying to figure out how to allocate the size for Array dynamically, as in I want it to be from 0 to a Number till grid independent solution is reached by a Condition.

So below is my program, I run it many time for different values of N. But then I don't want to do this, could you please help me with this


Thank You!

P.S: I have tried different ways but then since the entire program depends on N. I am not able to figure out what to do!




program HW1_P1_a
Integer, Parameter :: N=250 !Globally used variable, hence PARAMETER
Double Precision Del_X,A(N),B(N),C(N),D(N),PA,QA,RA,PB,QB,RB
Double Precision L,M,Theta_B,Theta_L, THETA(N)
Integer I,T_B,T_L,T_INF
! "Details about Arguments used: All must be Double Precision (except N)"
! N Number of mesh points
! Del_X Mesh size
! A() Lower diagonal of tridiagonal matrix
! B() Main diagonal of tridiagonal matrix
! C() Upper diagonal of tridiagonal matrix
! D() Right side vector (do not adjust for bcs)
! PA Left b.c.(K=1): Coefficient of Theta
! QA Left b.c.(K=1): Coefficient of dTheta/Del_X
! RA Left b.c.(K=1): Right hand side
! PB Right b.c.(K=N): Coefficient of Theta
! QB Right b.c.(K=N): Coefficient of dTheta/Del_X
! RB Right b.c.(K=N): Right hand side
! L Length of Fin
! T_B Temperature at Base of Fin
! T_L Temperature at the Tip of Fin
! T_inf Temperature of surrounding air
! Theta_B Temperature excess at Fin Base
! Theta_L Temperature excess at Fin Tip
! P Perimeter
! h Heat Transfer Coeffecient
! Ac Cross Sectional Area of Fin
! K Thermal Conductivity



! Inititalising values of the variables, "Given Conditions"
T_B=200
T_L=60
T_INF=20
L=0.1
P=0.09
Ac=7.E-4
h=15
K=80



! Equation for Excess Temperature

Theta_B=T_B-T_INF
Theta_L=T_L-T_INF

! Equation for Delta X and M

Del_X=L/(N-1) !Mesh size
M=SQRT((h*P)/(K*Ac))

! Fin equation derivatives are approximated using second order accurate finite difference and multiplying with (Del_X)^2

g=M**2
f=0

! Elements of the Matrix

DO I=1,N

A(i) = 1 - (Del_X/2)*f
B(i) = -(2 + (Del_X**2)*g)
C(i) = 1 + (Del_X/2)*f
D(i) = 0

ENDDO


! Dirchlet Boundary Condition at the Base

PA=1
QA=0
RA=Theta_B

! Dirclet Boundary condition at the Tip

PB=1
QB=0
RB=Theta_L



call THOMAS(N,Del_X,A,B,C,D,PA,QA,RA,PB,QB,RB,THETA)

Write(*,*) 'Values of Theta by Numerical Solutions for N:',N


DO I=1,N

write (*,*) (I*Del_X),THETA(I)

ENDDO
End HW1_P1_a








!************************************************* *********************
!
! SUBROUTINE: THOMAS
!
! DESCRIPTION: Use the Thomas algorithm to solve a system of
! linear algebraic equations where the
! coefficient matrix is tridiagonal. Boundary
! conditions can be Dirichlet, Neumann or Mixed.
!
! Left boundary condition is (X = Xmin):
!
! PA*Theta + QA*dTheta/Del_X = RA
!
! Right boundary condition is (X = Xmax):
!
! PB*Theta + QB*dTheta/Del_X = RB
!
! System of Equations:
!
! B1 C1 0 0 . . 0 = D1
! A2 B2 C2 0 . .
! 0 A3 B3 C3 . .
! . . . .
! . Ai Bi Ci = Di
! 0 . . . . An Bn = Dn
!
!
! RETURN ARGUMENTS:
! THETA() Function value
!
! VARIABLES:
!
! ************************************************** ********************
!
SUBROUTINE THOMAS(N,Del_X,A,B,C,D,PA,QA,RA,PB,QB,RB,THETA)

IMPLICIT DOUBLE PRECISION (A-H,O-Z)

PARAMETER (KDIM=10001)

DIMENSION A(*),B(*),C(*),D(*),THETA(*)
DIMENSION F(KDIM),DELTA(KDIM)

! -- Adjust coefficients for boundary conditions

B(1) = QA*B(1) + 2.D0*Del_X*PA*A(1)
C(1) = QA*(C(1) + A(1))
D(1) = QA*D(1) + 2.D0*Del_X*RA*A(1)

A(N) = QB*(A(N) + C(N))
B(N) = QB*B(N) - 2.D0*Del_X*PB*C(N)
D(N) = QB*D(N) - 2.D0*Del_X*RB*C(N)

! -- Forward elimination

F(1) = C(1)/B(1)
DELTA(1) = D(1)/B(1)
DO K=2,N
X1 = B(K) - A(K)*F(K-1)
X1 = 1.D0/X1
F(K) = C(K)*X1
DELTA(K) = X1*(D(K) - A(K)*DELTA(K-1))
END DO

! -- Back substitution

THETA(N) = DELTA(N)
DO K=N-1,1,-1
THETA(K) = DELTA(K) - F(K)*THETA(K+1)
END DO

RETURN
END