只适用于二维情况下的矩形，超级简单的程序，带的参数也少，是自己写的第一个程序~哈哈

初温为0，边界条件自己设定

        PROGRAM HEATTRANSFER_2D
        IMPLICIT NONE

        CHARACTER :: FILENAME*30,STEP*6
        INTEGER :: XMAX,YMAX,ZMAX,N,NMAX,STEPD
        INTEGER :: I,J,K,S(4)
        REAL :: DX,DY,DT,DT0,Kr,T0(1:4),A,B,T
        REAL,ALLOCATABLE :: X(:,:,:),Y(:,:,:),Z(:,:,:),TEMP(:,:,:)

        WRITE(*,*) "PLEASE INPUT FILENAME:"
        READ(*,*) FILENAME
        !----------------------------------------------------读取网格文件--------------------------------------------------------------------------
        OPEN(10,FILE=TRIM(FILENAME))       !Read Mesh File
        READ(10,*) XMAX,YMAX               !Get Mesh Size          
        ZMAX=1                             
        ALLOCATE(X(XMAX,YMAX,ZMAX),Y(XMAX,YMAX,ZMAX),Z(XMAX,YMAX,ZMAX))   !生成网格坐标存储空间
        ALLOCATE(TEMP(XMAX,YMAX,ZMAX))                                    !生成温度存储空间
        READ(10,*) (((X(I,J,K),I=1,XMAX),J=1,YMAX),K=1,ZMAX),(((Y(I,J,K),I=1,XMAX),J=1,YMAX),K=1,ZMAX),(((Z(I,J,K),I=1,XMAX),J=1,YMAX),K=1,ZMAX)
        CLOSE(10)
        !-----------------------------------------------------确定步长时间------------------------------------------------------------------------
        DX=(MAXVAL(X)-MINVAL(X))/XMAX
        DY=(MAXVAL(Y)-MINVAL(Y))/YMAX
        
        WRITE(*,*) 'Please input the rate of heat transfer k'
        READ(*,*) Kr
        DT0=(1/2.0/K)/(1/DX**2+1/DY**2)/2.0
		WRITE(*,*) 'Please input the step time, suggest it smaller than ',DT0
		READ(*,*) DT
		
		!------------------------------------------------------确定边界\初值条件-----------------------------------------------------------------------
		TEMP(1:XMAX,1:YMAX,1:ZMAX)=0
		!------------------确定边界条件类型------------------
		DO I=1,4
				WRITE(*,*) "Select Boundary Condition Type For"
				SELECT CASE(I)
				    CASE (1) 
				        WRITE(*,*) "Left side"
				    CASE (2)
				        WRITE(*,*) "Bottom"
				    CASE (3)
				        WRITE(*,*) "Right side"
				    CASE (4) 
				        WRITE(*,*) "Top"
				END SELECT
				WRITE(*,*) "1\Temperature"
				WRITE(*,*) "2\Rate of Temperatur Spread"
				READ(*,*) S(I)
		END DO
		!-----------------确定边界条件数值-------------------
		WRITE(*,*) "Please input the number for boundaries of left, bottom, right and up,use comma to seperate./The number should be wrote in 7 chracter and 3 float number(including comma)"
		READ(*,"(F7.2,','F7.3,',',F7.3,','F7.3)") T0(1:4)
		!----------------------------------------------------------------输入迭代次数\显示代数------------------------------------------------------------
		WRITE(*,*) 'Please input the number of steps'
		READ(*,*) NMAX
		WRITE(*,*) "Please input the number between two steps to show"
		READ(*,*) STEPD
		N=0
		T=0
100     N=N+1
        T=T+DT

		SELECT CASE(S(1))
		    CASE(1) 
		        TEMP(1,1:YMAX,1:ZMAX)=T0(1)
		    CASE(2) 
		        TEMP(1,1:YMAX,1:ZMAX)=(0-T0(1)*2*DX-TEMP(3,1:YMAX,1:ZMAX)+4*TEMP(2,1:YMAX,1:ZMAX))/3
        END SELECT
        SELECT CASE(S(2))
		    CASE(1) 
		        TEMP(1:XMAX,1,1:ZMAX)=T0(2)
		    CASE(2) 
		        TEMP(1:XMAX,1,1:ZMAX)=(0-T0(2)*2*DY-TEMP(1:XMAX,3,1:ZMAX)+4*TEMP(1:XMAX,2,1:ZMAX))/3
        END SELECT
        SELECT CASE(S(3))
		    CASE(1) 
		        TEMP(XMAX,1:YMAX,1:ZMAX)=T0(3)
		    CASE(2) 
		        TEMP(XMAX,1:YMAX,1:ZMAX)=(0-T0(3)*2*DX-TEMP(XMAX-2,1:YMAX,1:ZMAX)+4*TEMP(XMAX-1,1:YMAX,1:ZMAX))/3
        END SELECT
        SELECT CASE(S(4))
		    CASE(1) 
		        TEMP(1:XMAX,YMAX,1:ZMAX)=T0(4)
		    CASE(2) 
		        TEMP(1:XMAX,YMAX,1:ZMAX)=(0-T0(1)*2*DY-TEMP(1:XMAX,YMAX-2,1:ZMAX)+4*TEMP(1:XMAX,YMAX-3,1:ZMAX))/3
        END SELECT

        !-----------------------------------------------------二阶中心差分格式显式推进计算---------------------------------------------------------------------------
       	A=Kr*DT/DX**2
		B=Kr*DT/DY**2
		K=1
		XDO:DO I=2,XMAX-1
			YDO:DO J=2,YMAX-1
				TEMP(I,J,K)=A*(TEMP(I+1,J,K)+TEMP(I-1,J,K))+B*(TEMP(I,J+1,K)+TEMP(I,J-1,K))+(1-2*(A+B))*TEMP(I,J,K)
			END DO YDO
		END DO XDO
        !------------------------------------------------------收敛之判断-----------------------------------------------------------------------------------
        IF (MOD(N,STEPD) == 0) THEN
			WRITE(STEP,"(I6)") N
			OPEN (1,FILE='OUTPUT'//STEP//'_'//TRIM(FILENAME))
			WRITE(1,*) 'VARIABLES=','"X","Y","T"'
			WRITE(1,*) 'ZONE,I=',XMAX,',J=',YMAX,',F=POINT'
			DO I=1,XMAX
				DO J=1,YMAX
					WRITE(1,*) X(I,J,K),Y(I,J,K),TEMP(I,J,K)
				END DO
			END DO
			CLOSE(1)
		ENDIF
		
        WRITE(6,*) 'TIME=',T,',TEMPERATURE=',TEMP(XMAX/2,YMAX/2,K)
		IF(N < NMAX)  GOTO 100
		IF(N == NMAX) GOTO 101

101		WRITE(*,*) "THE COMPUTATION IS OVER!"
        
        WRITE(STEP,'(I6)') NMAX
        OPEN (2,FILE='OUTPUT'//STEP//'_'//TRIM(FILENAME))
        WRITE(2,*) 'VARIABLES=','"X","Y","T"'
        WRITE(2,*) 'ZONE,I=',XMAX,',J=',YMAX,',F=POINT'
        DO I=1,XMAX
            DO J=1,YMAX
            	WRITE(2,*) X(I,J,K),Y(I,J,K),TEMP(I,J,K)
            END DO
        END DO
        CLOSE(2)
        

        STOP
        END PROGRAM

PLot3D是一种在CFD领域很常见的格式。几乎所有的网格生成软件都支持这种格式的输出，因此搞清楚这种文件的格式的读入和输出对CFD研究人员的工作带来很大的帮助。这里给出的都是读入ASCII格式文件的代码，输出方式类似。

PLot3D的网格文件后缀名为"xyz"，可能是二进制（bin）格式或者文本文件格式(ASCII)。

1.单块网格（Single Grid）

 读入3维单块网格

2.多块网格（Multiple Grids）

读入3维多块网格（似乎GRIDGEN读入网格时需采用这种格式的文件，你可以将其改成输出程序，所输出的网格就可以用GRIDGEN读入了）