exphfeq.f90

一维热传导方程的求解

!************************************************************************
PROGRAM EXPHFEQ
      
      IMPLICIT NONE

!     I           RUNTIME INDEX
!     T           TIME
!     jmax       NR. OF GRID-POINTS
!     TMAX        ABSOLUTE TIME = 1 A
!     LAMBDA      HEAT CONDUCTION W M^-1 K^-1
!     DELTAT      TIME STEP  
!     Q           HEAT FLOW (W/M^2)       
!     ROHC        VOLUMETRIC HEAT CAPACITY (J/(M^3K)) = 
!                             1650 KG/M^3 * 837,2 J/(KG K)      
!     TEMP        TEMPERATURE IN THE PRESENT       
!     TEMPZ       TEMPERATURE IN THE FUTURE       
!     DELTAX      GRID DISTANCE (M)       
! **********************************************************************
      INTEGER           I,T,JMAX,JMAXS,TMAX
      INTEGER           AUSG, FAKTOR
      INTEGER           DELTAT
      PARAMETER(        JMAXS = 200)
      DOUBLE PRECISION  LAMBDA,Q,ROHC,RZ,INITEMP,BC
      DOUBLE PRECISION  TEMP(JMAXS),TEMPZ(JMAXS)
      DOUBLE PRECISION  DELTAX;
      DOUBLE PRECISION  REST,CFL
! **********************************************************************
!     READ PARAMETER
      OPEN(49,FILE='parameter.dat')
1     READ(49,*) LAMBDA
2     READ(49,*) ROHC
3     READ(49,*) Q
4     READ(49,*) RZ
5     READ(49,*) TMAX
6     READ(49,*) DELTAT
7     READ(49,*) INITEMP
8     READ(49,*) DELTAX
9     READ(49,*) jmax
10    READ(49,*) AUSG
11    READ(49,*) BC

! **********************************************************************
!     CHECK FOR CFL-STABILITY
      CFL=(LAMBDA/ROHC)*(DELTAT/DELTAX**2)
      WRITE(*,*) 'CFL = ', CFL
      IF(CFL.GT.0.5D0)THEN
      WRITE(*,*) 'ERROR WITH COURANT-FRIEDRICH-LEVY-CRITERIA'
      STOP
      ENDIF
! **********************************************************************
      OPEN(50,FILE='EXPL.OUT') 
!     INITIALISING OF THE TEMPERATURMATRIX DO T = 283,15 K = 10癈 
!     AVERAGE TEMPERTURE IN THE GROUND 
      DO I=1, jmax
      TEMP(I)=INITEMP
      ENDDO
! **********************************************************************
!     NUMBER OF OUTPUTS TO FILE = 20 
      AUSG=20; 
      FAKTOR = TMAX/AUSG
!     START CALCULATION 
      DO	T=0, TMAX, DELTAT
			REST = MOD(T,FAKTOR)
            DO I=1, JMAX
!     DIRICHLET BORDER CONDITION FIXED TEMPERATURE AT HIGHEST/FIRST GRID POINT 
            TEMP(1) = 283.15+BC
            TEMP(2) = 283.15+BC
            TEMP(3) = 283.15+BC
!     NOW THE MAIN CALCULATION 
            TEMPZ(I)=TEMP(I)+((LAMBDA*DELTAT)/(ROHC*DELTAX*DELTAX))*(TEMP(I-1)-2.D0*TEMP(I)+TEMP(I+1))+RZ
!     NEUMANN BORDER CONDITION: BORDER FLOW AT LOWEST/LAST GRID POINT 
!     SETTING LOCAL DERIVATION (ORTS ABLEITUNG) TO A CONCRETE VALUE = 0 
            TEMPZ(JMAX)=TEMPZ(JMAX-1)
!     THATS IT - THE NEUMANN BORDER CONDITION       
!     WRITE OUTPUT FILE                   
            IF(REST.EQ.0.0)THEN
                  WRITE(50,'(I8,1X,F20.8)') I, TEMP(I)
            ENDIF
!     END OF GRID POINTS LOOP       
            ENDDO
!     DELIVERY OF THE PROVISIONAL RESULT TO THE TIME LOOP 
            DO	I=1, JMAX
				TEMP(I)=TEMPZ(I)
            ENDDO
            IF (REST.EQ.0.0)THEN 
				WRITE(*,*) 'WRITE OUTPUT DO T (DAY) = ',(T/(60*60*24)),' T (SEC): ',T
!     SEE WHAT HAPPENS             
            ENDIF
!     END OF TIME LOOP       
      ENDDO      
END