[Fortran] 纯文本查看 复制代码
USE MSFLIB
PARAMETER IR=400,JR=400
INTEGERIS(0:IR+1,0:JR+1),TMAX,ISN(1:8),NSTATE,T,NR,IX0,IY0,IX,JY
!! 根据过去状态IS，定义现在状态为IS1；为了突出边界，特别定义ISN1
INTEGER IS1(0:IR+1,0:JR+1),ISN1(1:8)
WRITE(*,*)"PLEASE INPUT THE TIMESTEP "
READ(*,*)TMAX
ISEED=RTC()
IRC=IR/2 !=IR*ran(iseed)+1
JRC=JR/2
! 定义基体体积能为10,晶粒体积能为1
IS=8
IS(IRC,JRC)=1
!! 在循环前定义现在状态数组IS1的初始值
IS1=IS
OPEN(1,FILE="E:\LUKE.DAT")
DO T=1,TMAX
!! 每次循环前重新定义过去状态数组IS
IS=IS1
! 边界条件
IS(0,0:JR+1)=IS(IR,0:JR+1)
IS(IR+1,0:JR+1)=IS(1,0:JR+1)
IS(0:IR+1,0)=IS(0:IR+1,JR)
IS(0:IR+1,JR+1)=IS(0:IR+1,1)
!! 整个基体上，各个点上的状态都要根据规则改变，而非随机取点改变
DO IX=1,IR
DO JY=1,JR
ISN=(/IS(IX-1,JY-1),IS(IX-1,JY),IS(IX-1,JY+1),IS(IX,JY-1)
& ,IS(IX,JY+1),IS(IX+1,JY-1),IS(IX+1,JY),IS(IX+1,JY+1)/)
E0=COUNT(ISN.NE.IS(IX,JY))
! 如果不是晶粒边界，则跳出，重新循环
IF(E0.EQ.0)CYCLE
! 随机寻找一个相邻点
NR=8*RAN(ISEED)+1
NSTATE=ISN(NR)
! 判断与相邻点的能量差，并决定是否改变状态
E=COUNT(ISN.NE.NSTATE)
RD=RAN(ISEED)
IG=NSTATE-IS(IX,JY)
DE=E-E0+IG+2.5*RD-1.25
!! 用现在状态数组IS1记录边界状态的改变
IF(DE.LT.0.0)IS1(IX,JY)=NSTATE
ENDDO
END DO
!! 每循环20次在显示屏幕上刷新状态（颜色）
DO IX=1,IR
DO JY=1,JR
! IF(MOD(T,20).EQ.0)THEN
ISN1=(/IS1(IX-1,JY-1),IS1(IX-1,JY),IS1(IX-1,JY+1),IS1(IX,JY-1)
& ,IS1(IX,JY+1),IS1(IX+1,JY-1),IS1(IX+1,JY),IS1(IX+1,JY+1)/)
ISRE=SETCOLOR(IS(IX,JY))
! 如果是边界，定义特别颜色15,加以区分
IF(COUNT(ISN1.NE.IS1(IX,JY)).GT.0) ISRE=SETCOLOR(15)
ISRE=SETPIXEL(IX,JY)
! END IF
ENDDO
ENDDO
!!
! 记录循环次数和对应的晶粒面积
WRITE(1,*)T, SQRT(1.0*COUNT(IS.EQ.1))
ENDDO
CLOSE(1)
END
