PROGRAM decay
!**********************************************
!* Fortran code dealing with the decay problem.
!*
!* Author: J. Kaempf, 2008
!* Change: Daimu, 2023
!**********************************************
! Text after pronounciation marks are treated as comments.
! Every FORTRAN program starts with the "program" statement.
! The program name, "decay" here, should be different from
! the file name. Lower case or upper case doesn't matter in
! FORTRAN.
!
! Now comes the declaration section.
! Every symbol used later on needs to be declared.
INTEGER :: n ! time level
INTEGER :: ntot ! total number of iterations
INTEGER :: nout ! output every nout's iteration step
INTEGER :: mode ! use either explicit or implicit scheme
REAL :: C_explicit ! concentration at time level n
REAL :: CN_implicit ! "N" for new, concentration at next time level n+1
REAL :: C_implicit ! concentration at time level n
REAL :: CN_explicit ! "N" for new, concentration at next time level n+1
REAL :: CTRUE ! exact analytical solution for comparison
REAL :: CZERO ! Initial value of C
REAL :: dt ! time step
REAL :: kappa ! decay constant
REAL :: time ! time counter
REAL :: fac_explicit ! factor used in the scheme
REAL :: fac_implicit ! factor used in the scheme
! Now comes the initialization of variables and parameters.
CZERO = 100.0 ! initial concentration is 100
C_explicit = CZERO ! initial value is used for exact solution
C_implicit = CZERO
kappa = 0.0001 ! value of decay constant; you can also write this as kappa = 1.0e-4
dt = 3600. ! value of time step; dt = 60 seconds here
! Here is an example of an IF-statement that writes a message to the command screen,
! if the stability condition is violated.
!mode = 2 ! mode = 1 (explicit scheme); mode = 2 (implicit scheme)
fac_explicit = 1.0-dt*kappa
fac_implicit = 1.0/(1.+dt*kappa)
IF(mode == 1)THEN
if(fac<=0.0)write(6,*)'STABITY CRITERION ALERT: REDUCE TIME STEP'
END IF
! Total simulation time is 24 hours. How many interation steps is this?
ntot = 15.0*3600./dt ! Remember that an hour has 60 minutes of 60 seconds each
! Data output at every hour of the simulation.
nout = 1.0*3600./dt
! Open file for output.
! The first (unit) number is a referene for output (see below).
! The "file" statement specifies the desired file name.
! The statement "formatted" means ascii output.
! The status "unknown" implies new creation of a file if this does not exist,
! otherwise an existing file will be overwritten.
! Be careful not to overwrite files that you might need in the future.
open(10, file = 'daimu_decay.txt', form = 'formatted', status = 'unknown')
! Write initial time and concentration to this file.
WRITE(10,*)0,100.0,100.0,100.0
! Now comes the start of the iteration loop.
! It means that the statements between the DO-loop start and
! the corresponding "END DO" statement are repeated for n = 1 to
! n = ntot at steps of 1. In this case, statements are repeated
! a total number of 864000 times. Try to do this on a piece of paper...
!****** Start of iteration *******
DO n = 1,ntot
!*********************************
CN_explicit = C_explicit*fac_explicit ! prediction for next time step
CN_implicit = C_implicit*fac_implicit ! prediction for next time step
time = REAL(n)*dt ! time counter
CTRUE = CZERO*exp(-kappa*time) ! exact analytical solution
C_explicit = CN_explicit ! updating for upcoming time step
C_implicit = CN_implicit ! updating for upcoming time step
! Data output if the counter i is a constant integer multiple of nout.
! For instance, nout = 30 produces outputs at i = 30, 60. 90, ....
! Note that two equal signs (==) have to be used in the IF statement.
IF(mod(i,nout) == 0)THEN
! Output to unit 10. The associated file name is given above. The "*" creates an automated format.
! Here we produce three output columns: the first is the time in hours, the second our prediction, and the
! third the exact solution. IF-statements with more than 1 line have to use a "THEN" and an "END IF".
WRITE(10,*)time/3600.0, C_explicit, C_implicit, CTRUE
! We can also write a note a message to the screen here.
write(6,*)"Data output at time = ",time/3600.0," hours"
END IF
!****** End of iteration *******
END DO
!*******************************
STOP
END PROGRAM decay
import numpy as np
import matplotlib.pyplot as plt
plt.figure('Decay Problem')
data=np.loadtxt('E:\公众号\Ocean Modelling for Beginners\第二章\daimu_decay.txt').T;
res_explicit = data[1,:]
res_implicit = data[2,:]
analytical = data[3,:]
time = np.arange(0, 16, 1)
plt.plot(time, analytical, color='black', linestyle='dashed', label='analytical')
plt.plot(time, res_explicit, label='explicit')
plt.plot(time, res_implicit, label='implicit')
plt.legend()
plt.title('Decay Problem')
plt.xlabel('Time (hours)')
plt.ylabel('Concentration (100%)')
plt.ylim(0, 100)
plt.xlim(0, 15)
plt.xticks(time)
plt.grid(True)
plt.show()
联系客服