PROGRAM EXPHFEQ IMPLICIT NONE C I RUNTIME INDEX C T TIME C jmax NR. OF GRID-POINTS C TMAX ABSOLUTE TIME = 1 A C LAMBDA HEAT CONDUCTION W M^-1 K^-1 C DELTAT TIME STEP C Q HEAT FLOW (W/M^2) C ROHC VOLUMETRIC HEAT CAPACITY (J/(M^3K)) = C 1650 KG/M^3 * 837,2 J/(KG K) C TEMP TEMPERATURE IN THE PRESENT C TEMPZ TEMPERATURE IN THE FUTURE C DELTAX GRID DISTANCE (M) C ********************************************************************** 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 C ********************************************************************** C 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 C ********************************************************************** C 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 C ********************************************************************** OPEN(50,FILE='EXPL.OUT') C INITIALISING OF THE TEMPERATURMATRIX DO T = 283,15 K = 10°C C AVERAGE TEMPERTURE IN THE GROUND DO I=1, jmax TEMP(I)=INITEMP ENDDO C ********************************************************************** C NUMBER OF OUTPUTS TO FILE = 20 AUSG=20; FAKTOR = TMAX/AUSG C START CALCULATION DO T=0, TMAX, DELTAT REST = MOD(T,FAKTOR) DO I=1, JMAX C 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 C NOW THE MAIN CALCULATION TEMPZ(I)=TEMP(I)+((LAMBDA*DELTAT)/(ROHC*DELTAX*DELTAX))* & (TEMP(I-1)-2.D0*TEMP(I)+TEMP(I+1))+RZ C NEUMANN BORDER CONDITION: BORDER FLOW AT LOWEST/LAST GRID POINT C SETTING LOCAL DERIVATION (ORTS ABLEITUNG) TO A CONCRETE VALUE = 0 TEMPZ(JMAX)=TEMPZ(JMAX-1) C THATS IT - THE NEUMANN BORDER CONDITION C WRITE OUTPUT FILE IF(REST.EQ.0.0)THEN WRITE(50,'(I8,1X,F20.8)') I, TEMP(I) ENDIF C END OF GRID POINTS LOOP ENDDO C 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 C SEE WHAT HAPPENS ENDIF C END OF TIME LOOP ENDDO END