Tek-Tips is the largest IT community on the Internet today!

Members share and learn making Tek-Tips Forums the best source of peer-reviewed technical information on the Internet!

Register Log in
  • Congratulations strongm on being selected by the Tek-Tips community for having the most helpful posts in the forums last week. Way to Go!

Creating Subroutines/Subfunctions

Status
Not open for further replies.
takafoo

takafoo

Programmer
Apr 15, 2011
12
AU
I am playing around with the below code and would like to 'clean' it up a little by making a subroutine or function but I am having trouble when I try. Can anyone help me with what part I can/should break down and how I would go about doing it?

Code: 
PROGRAM A1004

   IMPLICIT NONE

! Declare Variables

   COMPLEX(KIND=8) :: mat_i = (0.0d0,1.0d0)
   COMPLEX(KIND=8) :: temp, r, t
   COMPLEX(KIND=8), DIMENSION(:,:), ALLOCATABLE :: mat
   COMPLEX(KIND=8), DIMENSION(:), ALLOCATABLE :: vec

   REAL(KIND=8) :: n_x, n_0, lambda, L
   REAL(KIND=8) :: k, beta, pi, exp, step, x_i, x_0
   REAL(KIND=8) :: R0, T0, start_time, end_time
   REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: x

   INTEGER :: m, numEq, i, j
   INTEGER :: allocate_status=0, IO_status
   INTEGER :: NRHS, LD_mat, LD_vec, N, INFO
   INTEGER, DIMENSION(:), ALLOCATABLE :: IPIV

   CHARACTER(LEN=20) :: filename, output

! Request File Information

100 CONTINUE

   PRINT *, "Please specify parameter filename: "
   READ (*,*) filename
   WRITE(*,*)

   PRINT *, "Please specify output filename: "
   READ (*,*) output
   WRITE(*,*)

! Load Parameters File

  OPEN(1, file=filename, status='old', IOSTAT=IO_status)
  IF (IO_status /= 0) THEN
    WRITE(*,*) "----- ERROR -----"
    WRITE(*,*) "Cannot open file ", TRIM(filename), " . File cannot be located."
    WRITE(*,*) "Please check file and try again."
    WRITE(*,*)
    GOTO 100
  ENDIF

    READ(1,*) lambda
    READ(1,*) n_x
    READ(1,*) n_0
    READ(1,*) L
    READ(1,*) m
  CLOSE(1)

! Output File Contents

   WRITE(*,*) 'Wavelength = ' , lambda
   WRITE(*,*) 'Refractive Index Function = ' , n_x
   WRITE(*,*) 'Refractive Index = ' , n_0
   WRITE(*,*) 'Length = ' , L
   WRITE(*,*) 'Grid Point Parameter = ',  m
   WRITE(*,*)

! Set Variables

   pi = 3.141592653589793239
   k = 2.0d0 * pi / lambda
   beta = n_0 * k
   step = L / m
   numEq = m + 3
   x_0 = 0.0d0

! Error Checking

   IF (lambda <= 0) THEN
     WRITE(*,*) "ERROR! lambda <= 0"
     GOTO 200
   ENDIF

   IF (n_x <= 0) THEN
     WRITE(*,*) "ERROR! Refractive Index Function <= 0"
     GOTO 200
   ENDIF

   IF (n_0 <= 0) THEN
     WRITE(*,*) "ERROR! Refractive Index <= 0"
     GOTO 200
   ENDIF

  IF (L <= 0) THEN
     WRITE(*,*) "ERROR! l <= 0"
     GOTO 200
   ENDIF

  IF (m <= 2) THEN
     WRITE(*,*) "ERROR! Grid Point Parameter < 2"
     GOTO 200
   ENDIF

! Bulk Coding

   WRITE(*,*) "====================================="
   WRITE(*,*)
   WRITE(*,*) "EXECUTING PROGRAM"
   WRITE(*,*)

   CALL cpu_time(start_time)

   ALLOCATE(mat(numEq,numEq),vec(numEq),IPIV(numEq), STAT=allocate_status)

   IF (allocate_status /= 0) then
      WRITE(*,*) "Not enough memory for the allocation"
      GOTO 200
   ENDIF

   vec = (0.0d0, 0.0d0)
   vec(1) = CMPLX(0.0d0,step*beta)
   vec(numEq - 1) = (1.0d0, 0.0d0)

   mat = (0.0d0,0.0d0)
   mat(1,1) = (-1.0d0,0.0d0)
   mat(1,2) = (1.0d0,0.0d0)
   mat(1,numEq-1) = CMPLX(0.0d0, step*beta)
   mat(numEq-2,numEq-3) = (-1.0d0,0d0)
   mat(numEq-2,numEq-2) = (1.0d0,0d0)

   temp = exp(CMPLX(0.0d0,beta*l))
   mat(numEq-2,numEq) = CMPLX(-mat_i*step*beta*temp)
   mat(numEq-1,1) = (1.0d0,0.0d0)
   mat(numEq-1,numEq-1) = (-1.0d0,0.0d0)
   mat(numEq,numEq-2) = (1.0d0,0.0d0)
   mat(numEq,numEq) = CMPLX(-temp)

   DO i=2, numEq - 3
      mat(i,i) = -2.0d0 + step**2 * n_x**2 * k**2
      mat(i,i-1) = (1.0d0,0.0d0)
      mat(i,i+1) = (1.0d0,0.0d0)
   END DO

! Call LAPACK Routine

   N = numEq
   NRHS = 1
   LD_mat = numEq
   LD_vec = numEq

   CALL ZGESV(N, NRHS, mat, LD_mat, IPIV, vec, LD_vec, INFO)

   IF (INFO == 0) THEN
      r = vec(numEq -1)
      t = vec(numEq)
   END IF

   R0 = ABS(r)**2
   T0 = ABS(t)**2

   CALL cpu_time(end_time)

! Outputting Solution to screen

   WRITE(*,*) "====================================="
   WRITE(*,*)

   WRITE(*,*) "Reflectance = ", R0
   WRITE(*,*) "Trasmittance = ", T0
   WRITE(*,*) "Show that R and T satisfy the energy conservation relation: "
   WRITE(*,*) "R + T = ", R0 + T0
   WRITE(*,*)
   WRITE(*,*) "Running time : ", end_time - start_time

! Outputting Solution to file

   OPEN(2,file=output,status='unknown', iostat=IO_status)
   IF (IO_status /= 0) THEN
      WRITE(*,*) "Error opening the file ", TRIM(output)
      GOTO 200
   ELSE

   WRITE(2,*) "Using the below: "
   WRITE(2,*) 'Wavelength = ' , lambda
   WRITE(2,*) 'Refractive Index Function = ' , n_x
   WRITE(2,*) 'Refractive Index = ' , n_0
   WRITE(2,*) 'Length = ' , L
   WRITE(2,*) 'Grid Point Parameter = ',  m
   WRITE(2,*)

   WRITE(2,*) "Solution: "
   WRITE(2,*) "Reflectance = ", R0
   WRITE(2,*) "Trasmittance = ", T0
   WRITE(2,*)
   WRITE(2,*) "Show that R and T satisfy the energy conservation relation: "
   WRITE(2,*) "R + T = ", R0 + T0
   WRITE(2,*)
   WRITE(2,*) "Running time : ", end_time - start_time

   ENDIF
   CLOSE(2)

   DEALLOCATE(mat, vec, IPIV)

200 CONTINUE

END PROGRAM A1004
 
Sort by date Sort by votes
The code is rather clean, short and simple. Why do you want to split it into subroutines ?

Just few remarks :
- the variable pi has not the precision you expect : pi = 3.141592653589793239d0 is much better (don't forget the d0 to indicate double precision)
- why do you duplicate arguments of ZGESV ?
CALL ZGESV(numEq, 1, mat, numEq, IPIV, vec, numEQ, INFO) seems OK

Now, if you really want to split a code into subroutines, you must first isolate main tasks. For instance :
- reading data
- checking data
- building the equations
- solving the problem

After that, you have to pass to each task (a subroutine by task) the information needed by that task (input and output). The main program will just call the subroutines in the right order. This is not complicated.

It is also usual to put the subroutines into one or several modules.


François Jacq
 
Upvote 0 Downvote
Thanks for your points. I thought it would be cleaner if I pulled the below into a subroutine?

Code: 
   vec = (0.0d0, 0.0d0)
   vec(1) = CMPLX(0.0d0,step*beta)
   vec(numEq - 1) = (1.0d0, 0.0d0)

   mat = (0.0d0,0.0d0)
   mat(1,1) = (-1.0d0,0.0d0)
   mat(1,2) = (1.0d0,0.0d0)
   mat(1,numEq-1) = CMPLX(0.0d0, step*beta)
   mat(numEq-2,numEq-3) = (-1.0d0,0d0)
   mat(numEq-2,numEq-2) = (1.0d0,0d0)

   temp = exp(CMPLX(0.0d0,beta*l))
   mat(numEq-2,numEq) = CMPLX(-mat_i*step*beta*temp)
   mat(numEq-1,1) = (1.0d0,0.0d0)
   mat(numEq-1,numEq-1) = (-1.0d0,0.0d0)
   mat(numEq,numEq-2) = (1.0d0,0.0d0)
   mat(numEq,numEq) = CMPLX(-temp)

   DO i=2, numEq - 3
      mat(i,i) = -2.0d0 + step**2 * n_x**2 * k**2
      mat(i,i-1) = (1.0d0,0.0d0)
      mat(i,i+1) = (1.0d0,0.0d0)
   END DO
 
Upvote 0 Downvote
Status
Not open for further replies.

Similar threads

tesaoghost
  • Locked
  • Question
problem with functions and subroutines that returns arrays 1
Replies
13
Views
151
tesaoghost
tesaoghost
mopen
  • Locked
  • Question
subroutines within the same module
Replies
1
Views
38
salgerman
salgerman
billgray1234
  • Locked
  • Question
Fortran 90, PUBLIC subroutines, and &quot;local&quot; variables within subroutines 3
Replies
10
Views
284
FJacq
FJacq
billgray1234
  • Locked
  • Question
recursive subroutines and allocatable arrays 3
Replies
6
Views
73
gummibaer
gummibaer
billgray1234
  • Locked
  • Question
fortran 90/95, recursive subroutines, and optional arguments
Replies
1
Views
79
xwb
xwb

Part and Inventory Search

Sponsor

Back
Top