(原创)二维圆形颗粒流离散元计算程序--Fortran语言

module global

!-------------------------------------------------------------------!

! 第一部分参数: 划分网格相关 !

!-------------------------------------------------------------------!

real bj_ptx(20),bj_pty(20),bj_ptnum !ptx,pty边界折点的x,y坐标值,边界点的个数

real bj_ptxmin,bj_ptxmax,bj_ptymin,bj_ptymax !边界点坐标ptx,pty的极值

real bj_areasum !由折线形成的闭合区域的总面积

integer ran_num,c_num !ran_num随机种子个数,c_num指的是定位好的单元个数

real c_ox(10000),c_oy(10000),c_R(10000),theta(10000) !单元的特征值,分别为坐标值和半径值,转角

real ran_Rmean !随机单元半径的平均值

integer jc !与单元i交叉的单元个数

integer c_cnumJC(10000),c_cNoJC(1000) !c_cnumJC记录与单元i交叉的单元个数,c_cNoJC与单元i交叉的单元号码

real c_codist,c_cRsum !c_codist:两单元的圆心距离,c_cRsum:两单元的半径相加,c_cRsum=Ri+Rmk

real c_Rmin !单元半径的最小值

integer i,j,k !全局变量

real mesh_layer,mesh_lapse(10) !圆孔筛的层数,存放各层圆孔筛孔径round-hole mesh diameter的递减系数

real,parameter::pi=3.1415926

!-------------------------------------------------------------------!

! 第二部分参数:离散元参数的确定 !

!-------------------------------------------------------------------!

integer js_step,js_stepsum !js_step专门用来控制计算时步的变量,js_stepsum循环计算的总时步数

integer seekJC_step,save_step !seekJC_step每计算多少时步检索一次单元的交叉情况,save_step:每计算多少时步保存一次

character knksfix_flag !所有单元的刚度kn,ks是否统一

character damp_flag,rebound_flag !damp_flag:单元是否考虑阻尼,rebound_flag:单元是否反弹？

character verticalbj_flag,liju_flag !是否考虑竖向边界,是否考虑力矩作用

real c_thk,c_niu,c_E,c_dent,c_fric,c_cohe !单元的厚度,泊松比,弹模,密度,摩擦系数,粘聚力

real c_mass(10000),c_weight(10000) !质量,平均质量,重力

integer kn,ks,c_kn(10000),c_ks(10000) !kn和ks为法向和切向刚度系数,如果所有单元的刚度系数都一致,则采用kn,ks,否则采用c_kn(),c_ks() integer cw_nstif,cw_sstif !单元与边界墙之间的法向和切向刚度系数

real cw_fric,cw_coh !单元与边界墙的摩擦系数,粘聚力

real c_Ip(10000) !圆盘的极惯性距。应该是转动惯量 polar monent of inertia

real ray_a,ray_b,damp_cn,damp_cs !Rayleigh常数的a,b值(适用于绝对速度),接触阻尼系数(适用于相对速度)

integer c_cnumJCsum !单元i之前所有单元(i-1)的交叉单元个数

integer c_bjNoJC(10000) !单元交叉的边界编号

integer fig_color !定义图形的颜色。

!-------------------------------------------------------------------!

! 第三部分参数:离散元计算相关 !

!-------------------------------------------------------------------!

real laplength !两单元交叉时,圆心连线方向上交叉部分的长度

real lap_nlength,lap_slength,lap_xlength,lap_ylength !单元与单元之间的法向、切向,单元与墙之间的叠合量real vx_XD,vy_XD,vn_XD,vs_XD !单元圆心的x,y,法向n,切向s的相对速度

real un_XD,us_XD !单元圆心的相对位移的法向和切向分量

real delta_Fx,delta_Fy,delta_M !在一个deltat间隔内,单元i所受的法向力和切向力在X和Y方向的合力,及力矩real cohe_F !临时的粘聚力,注意此时cohe_F的单位为(N)。

real cos_gejiao,sin_gejiao !gejiao单元i与其它单元相交弦线在单元i的圆心角的(一半),注意是一半

real c_gelength !单元i与其它单元或墙体相交弦线的长度

real cos_alpha,sin_alpha !沿圆心连线方向单位矢量的x,y分量,!!!!注意,这里的相对一律是单元i相对于其它单元real Fs_max !最大的剪切力。

real damp_Fn,damp_Fs !圆盘单元接触处法向和切向阻尼分量(接触阻尼系数*相对速度)

real deltat_yita !时步影响系数系数。deltat_yita=0.01-0.1之间。

real deltat_ok,deltat_max !step time: 时步的大小。deltat_max:参考时步的大小值。

real c_Fxi,c_Fyi,c_Mi,w_Fxi,w_Fyi,w_Mi !单元i所受单元作用力和墙壁作用力

real,parameter::ag=9.8

!nt时刻:即初时刻,各形参的声明, initial

real vx_ini(10000),vy_ini(10000),angv_ini(10000),ax_ini(10000),ay_ini(10000),anga_ini(10000)

real ux_ini(10000),uy_ini(10000),u_ini(10000)

real Fx_ini(10000),Fy_ini(10000),F_ini(10000),M_ini(10000)

!nt+sti/2时刻:即半个时步时,各形参的声明, middle

real vx_mid(10000),vy_mid(10000),angv_mid(10000),ax_mid(10000),ay_mid(10000),anga_mid(10000)

real ux_mid(10000),uy_mid(10000),u_mid(10000)

real Fx_mid(10000),Fy_mid(10000),F_mid(10000),M_mid(10000)

!nt+sti时刻:即时步末时,各形参的声明, end

real vx_end(10000),vy_end(10000),angv_end(10000),ax_end(10000),ay_end(10000),anga_end(10000)

real ux_end(10000),uy_end(10000),u_end(10000)

end module

!------------------------------------------------------------

program DEM_XP !二维离散元计算程序

use global

implicit none

call DEMcell() !在CAD中产生随机离散单元

call cell_Information()

call DEMcomputation() !离散元计算

end

!---------------------------------------------------------------------!

! (1) 主子程序: 离散元网格的划分 !

!---------------------------------------------------------------------!

subroutine DEMcell()

use global

implicit none

integer cellmode_flag !离散元的生成方式

open (10001,file='region.txt') !region是事先存放边界点的txt文件

bj_ptnum=5 !形成一个闭合矩形

write (*,*) "请顺时针方向输入各折点的x,y坐标(>0):"

do i=1,bj_ptnum,1

read (10001,*) bj_ptx(i),bj_pty(i)

write(*,"(A2,I1,A23,F6.3,A1,F6.3)") "第",i,"个边界折点的x和y坐标为:",bj_ptx(i),",",bj_pty(i) end do

close (10001)

open (10002,file='bj_box.scr',status='replace') !将区域边界输出到CAD

do i=1,bj_ptnum-1,1

write(10002,"(A4)") "LINE"

write(10002,'(E12.6,",",E12.6)') bj_ptx(i),bj_pty(i)

write(10002,'(E12.6,",",E12.6)') bj_ptx(i+1),bj_pty(i+1)

write(10002,'(A1)') !这一行很重要,其作用是换行,并且不产生任何字符。

close (10002)

bj_ptxmin=minval(bj_ptx)

bj_ptxmax=maxval(bj_ptx)

bj_ptymin=minval(bj_pty)

bj_ptymax=maxval(bj_pty)

bj_areasum=(bj_ptxmax-bj_ptxmin)*(bj_ptymax-bj_ptymin)

10003 write(*,*) "请选择离散单元的生成方式:1随机,2均匀布置,3CAD导入"

read(*,*) cellmode_flag

if (cellmode_flag==1) then

write(*,*) "请输入要产生离散单元的个数(推荐值800)?" !以便知道单元的近似半径,并初步给定网格的边长

read(*,*) ran_num

ran_Rmean=sqrt(bj_areasum/(ran_num*pi)) !孔隙率为0时,离散单元的平均半径,因为要对半径进行筛,所以半径大点没关系 call randomcell()

elseif (cellmode_flag==2) then

call uniformcell()

elseif (cellmode_flag==3) then

call CADcell()

else

write(*,*) "请注意选择单元生成方式:1-3"

goto 10003

end if

return

end subroutine

!---------------------------------------------------------------------!

! (1.1) 直接法随机网格划分子程序 !

!---------------------------------------------------------------------!

subroutine randomcell()

use global

implicit none

call rancellradius() !产生所有单元的随机半径

call celltoCAD() !将定位好的数据写到CAD格式的文件

return

end subroutine

!---------------------------------------------------------------------!

! (1.1.1) 随机单元半径的产生 !

!---------------------------------------------------------------------!

subroutine rancellradius() !蒙特卡洛方法产生随机半径

use global

implicit none

real ran_areasum !随机半径的所有单元的面积

real ran_porosum !随机单元的初始孔隙率 porosity, 1-c_areasum/bj_areasum

integer ran_step !产生随机半径的步数

real ran_miu,ran_sigma !随机数服从正态分布的均值和方差

real ran_bin,ran_binsum !ran_bin为0-1之间的随机数，和为ran_binsum,

real ran_ok !ran_ok最终合乎要求的随机数

character agree_flag !子程序最后用来存储对本子程序产生随机数的满意与否"Y","N"

character mesh_flag !离散元的生成方式

11101 write(*,*) "请选择对单元半径的递减处理,默认层数10,默认递减系数0.85,是否接受,是(Y)否(N)" read(*,*) mesh_flag

if (mesh_flag=="Y") then

mesh_layer=10 !圆孔筛的层数,默认值为5,最大值为10

do i=1,mesh_layer,1

mesh_lapse(i)=0.85 !圆孔筛各层孔径的递减系数,默认值为0.85

end do

elseif (mesh_flag=="N") then

write(*,*) "请选择对单元半径的递减层数(1~10)和递减系数(0.01~0.99)"

read(*,*) mesh_layer

do i=1,mesh_layer,1

read(*,*) mesh_lapse(i)

end do

write(*,*) "请注意选择Y或N,递减层数的范围为:1-10,递减系数的范围为:0.01-0.99"

goto 11101

end if

11102 ran_step=0 !控制总的产生随机数的次数,直至产生单元的总面积小于整个区域的面积为止 write(*,"(A29,F6.4)") "孔隙率为0时离散单元的半径(m):",ran_Rmean

write(*,*) "请输入服从正态分布的单元半径的均值和方差"

read(*,*) ran_miu,ran_sigma

11103 ran_areasum=0

do i=1,ran_num,1

11104 ran_binsum=0 !ran_binsum存放68个0-1的随机数的总和

do j=1,100,1 !100个0-1的随机数

!call random_seed () ! 系统根据日期和时间随机地提供种子

call random_number (ran_bin) !调用库随机函数,产生0到1之间的随机数

ran_binsum=ran_binsum+ran_bin

end do

ran_ok=ran_miu+ran_sigma*(ran_binsum-100/2)/((100/12)**0.5)

if (ran_ok<=0.0) then

goto 11104 !最终随机数<=0,重新生成随机单元。

end if

c_R(i)=ran_ok

ran_areasum=ran_areasum+pi*c_R(i)**2

end do

ran_porosum=(bj_areasum-ran_areasum)/bj_areasum !全部单元排列(之间不交叉)后产生的孔隙率 if (ran_porosum<0) then

ran_step=ran_step+1

if (ran_step<100) then

goto 11103

else

write(*,"(A4,I3,1x,A37)") "经过",100,"步自动选取,仍然没有合适的,需调整参数!"

goto 11102

end if

write(*,"(A11,F6.3,A11,F6.3,A13,F5.2,A1)") "区域总面积:",bj_areasum,"单元总面积:",ran_areasum,"单元的孔隙率:",100*ran_porosum,"%" write(*,*) , "生成的单元是否符合要求?是(Y)或否(N)?"

read(*,*) agree_flag

if (agree_flag=="N") then

goto 11102

end if

return

end subroutine

!---------------------------------------------------------------------!

! (1.1.2) 根据单元的交叉情况重新定位单元 !

!---------------------------------------------------------------------!

subroutine cellrelocate()

use global

implicit none

integer jc_num !所有单元交叉的总次数,也就是goto 11201总步数

integer button(10) !圆孔筛的开关,以免产生累乘

button=0

jc_num=0 !记录整个单元定位过程,单元交叉的总次数

do i=1,ran_num,1

11201 call rancellcoordinate() !产生单元i的随机圆心坐标

jc=0 !如果存在单元与i单元交叉,则jc=1

do j=1,i-1,1 !只需判断1~i-1 !如果i=1,则这一步不执行

c_codist=sqrt((c_ox(i)-c_ox(j))**2+(c_oy(i)-c_oy(j))**2)

c_cRsum=c_R(i)+c_R(j)

if (c_codist

jc=1 ! 返回单元i交叉的结果，找到一个交叉单元即返回,所以jc=0或1

end if

if (jc==1) goto 11202 !找到一个交叉单元,则单元i完成检索

end do

do while (c_cnumJC(i)/=0) !即c_cnumJC(i)=1，存在与单元i交叉的单元

jc_num=jc_num+1

select case(jc_num)

case(1:20000) !单元交叉次数为20000(单元i交叉时,则重新产生中心坐标,如此循环,因此交叉次数会很多) goto 11201

case(20001:40000)

button(1)=button(1)+1

if (button(1)==1) then

write(*,*) "进入第1个孔筛,半径递减系数为:",mesh_lapse(1)

do j=i,ran_num,1

c_R(j)=c_R(j)*mesh_lapse(1)

end do

end if

goto 11201

case(40001:60000)

button(2)=button(2)+1

if (button(2)==1) then

write(*,*) "进入第2个孔筛,半径递减系数为:",mesh_lapse(2)

do j=i,ran_num,1

c_R(j)=c_R(j)*mesh_lapse(2)

end do

end if

goto 11201

case(60001:80000)

button(3)=button(3)+1

if (button(3)==1) then

write(*,*) "进入第3个孔筛,半径递减系数为:",mesh_lapse(3)

do j=i,ran_num,1

c_R(j)=c_R(j)*mesh_lapse(3)

end do

end if

goto 11201

button(4)=button(4)+1

if (button(4)==1) then

write(*,*) "进入第4个孔筛,半径递减系数为:",mesh_lapse(4) do j=i,ran_num,1

c_R(j)=c_R(j)*mesh_lapse(4)

end do

end if

goto 11201

case(100001:110000)

button(5)=button(5)+1

if (button(5)==1) then

write(*,*) "进入第5个孔筛,半径递减系数为:",mesh_lapse(5) do j=i,ran_num,1

c_R(j)=c_R(j)*mesh_lapse(5)

end do

end if

goto 11201

case(110001:120000)

button(6)=button(6)+1

if (button(6)==1) then

write(*,*) "进入第6个孔筛,半径递减系数为:",mesh_lapse(6) do j=i,ran_num,1

c_R(j)=c_R(j)*mesh_lapse(6)

end do

end if

goto 11201

case(120001:130000)

button(7)=button(7)+1

if (button(7)==1) then

write(*,*) "进入第7个孔筛,半径递减系数为:",mesh_lapse(7) do j=i,ran_num,1

c_R(j)=c_R(j)*mesh_lapse(7)

end if

goto 11201

case(130001:140000)

button(8)=button(8)+1

if (button(8)==1) then

write(*,*) "进入第8个孔筛,半径递减系数为:",mesh_lapse(8) do j=i,ran_num,1

c_R(j)=c_R(j)*mesh_lapse(8)

end do

end if

goto 11201

case(140001:150000)

button(9)=button(9)+1

if (button(9)==1) then

write(*,*) "进入第9个孔筛,半径递减系数为:",mesh_lapse(9) do j=i,ran_num,1

c_R(j)=c_R(j)*mesh_lapse(9)

end do

end if

goto 11201

case(150001:160000)

button(10)=button(10)+1

if (button(10)==1) then

write(*,*) "进入第10个孔筛,半径递减系数为:",mesh_lapse(10) do j=i,ran_num,1

c_R(j)=c_R(j)*mesh_lapse(10)

end do

end if

goto 11201

case default

goto 11203

end select

end do

11203 if (i

c_num=i

else !此时循环进行完毕,但返回的是i=ran_num+1

c_num=ran_num

end if

return

end subroutine

!---------------------------------------------------------------------!

! (1.1.3) 随机单元中心坐标的产生 !

!---------------------------------------------------------------------!

subroutine rancellcoordinate()

use global

implicit none

real xy_x,xy_y !x和y方向0-1之间随机数

real xy_a,xy_b,xy_c,xy_d !网格中心离4个边界的距离

11301 call random_number (xy_x) !调用库随机函数,产生随机的x,y坐标

call random_number (xy_y)

c_ox(i)=xy_x*(bj_ptxmax-bj_ptxmin)+bj_ptxmin !注意如果区域不是矩形,此处应特殊处理 c_oy(i)=xy_y*(bj_ptymax-bj_ptymin)+bj_ptymin

xy_a=c_ox(i)-(c_R(i)+bj_ptxmin) !左x边，用数组来控制单元不能超越边界

xy_b=c_oy(i)-(c_R(i)+bj_ptymin) !下y边

xy_c=bj_ptxmax-(c_ox(i)+c_R(i)) !右x边

xy_d=bj_ptymax-(c_oy(i)+c_R(i)) !上y边

if (xy_a<0.0 .or. xy_b<0.0 .or. xy_c<0.0 .or. xy_d<0.0 ) then

goto 11301

end if

return

!---------------------------------------------------------------------!

! (1.1.4) 将单元的信息写成 CAD 格式 !

!---------------------------------------------------------------------!

subroutine celltoCAD()

use global

implicit none

real c_areasum !最终定位后单元的面积

real c_porosum !单元全部定位后的孔隙率

c_areasum=0

do i=1,c_num,1

c_areasum=c_areasum+pi*c_R(i)**2 !所有单元的面积。

end do

c_porosum=1-c_areasum/bj_areasum

jc=0 !定位完后,找出交叉单元的个数

do i=1,c_num,1

do j=i+1,c_num,1 !只需要将单元i与单元i+1~c_num进行判断即可,该步实际执行到i=c_num-1 c_codist=sqrt((c_ox(i)-c_ox(j))**2+(c_oy(i)-c_oy(j))**2)

c_cRsum=c_R(i)+c_R(j)

if (c_codist

jc=jc+1

end if

end do

end do

write(*,"(A15,I3,A10)") "单元定位后,更有",jc,"个单元相交"

if (jc==0) then

write(*,*) "单元定位成功! 随机单元的信息保存在cell_division.scr中,可以调用查看"

write(*,"(A24,I5,A21,F5.2,A1)") "实际生成了单元总个数:",c_num,"定位后单元的孔隙率为:",100*c_porosum,"%" else

write(*,"(A17,I4,A12)") "单元定位失败,存在",jc,"个单元相交"

end if

call cellsort() !对单元进行排序,这样画图时方便查看

do i=1,c_num,1

write(11401,"(A6)") 'CIRCLE' !生成单元的SCR文件。

write(11401,'(E12.6,",",E12.6)') c_ox(i),c_oy(i)

write(11401,"(E12.6)") c_R(i)

write(11401,"(A5)") 'ENTER' !输入文字

write(11401,"(A4)") 'TEXT' !输入文字

write(11401,"(A1)") 'J' !对正

write(11401,"(A1)") 'C' !中心

write(11401,'(E12.6,",",E12.6)') c_ox(i),c_oy(i) !输入文字的位置 write(11401,"(E12.6)") (bj_ptymax-bj_ptymin)/200 !默认字体的高度 write(11401,"(A1)") '0' !旋转角度

if (i<10) then

write(11401,"(I1)") i

else if (i<99) then

write(11401,"(I2)") i

else if (i<999) then

write(11401,"(I3)") i

else if (i<9999) then

write(11401,"(I4)") i

else if (i<99999) then

write(11401,"(I5)") i

end if

write(11401,"(A1)") !产生一空行，目的是结束CAD中的text命令 end do

close(11401)

return

end subroutine

!---------------------------------------------------------------------! ! (1.1.4.1) 将单元的按 Y 坐标进行排序 ! !---------------------------------------------------------------------! subroutine cellsort() !对单元按Y由大到小进行排序,这样画图时方便查看

implicit none

real y_max,temp !随机单元y坐标的最大值

do i=1,c_num,1

y_max=c_oy(i)

do j=i+1,c_num,1

if (y_max

temp=c_oy(j)

c_oy(j)=c_oy(i)

c_oy(i)=temp

temp=c_ox(j)

c_ox(j)=c_ox(i)

c_ox(i)=temp

temp=c_R(j)

c_R(j)=c_R(i)

c_R(i)=temp

y_max=c_oy(i)

else if(y_max==c_oy(j)) then !假如y坐标相同,则以x坐标由大到小排序 if(c_ox(i)

temp=c_oy(j)

c_oy(j)=c_oy(i)

c_oy(i)=temp

temp=c_ox(j)

c_ox(j)=c_ox(i)

c_ox(i)=temp

temp=c_R(j)

c_R(j)=c_R(i)

c_R(i)=temp

y_max=c_oy(i)

end if

end if

end do

return

end subroutine

!---------------------------------------------------------------------!

! (1.2) 均匀离散单元的生成 !

!---------------------------------------------------------------------!

subroutine uniformcell() !根据单元的层数确定离散单元

use global

implicit none

real c_layer !单元在竖向上的层数

character bj_shape !边界形状,矩形或其它

integer lay_cnum(100) !每层的单元个数

!确定单元的层数。

write(*,*) "请确定要生成多少单元层数(参考值20)"

read (*,*) c_layer

ran_Rmean=(bj_ptymax-bj_ptymin)/(2*c_layer) !确定半径的大小。

!注意：要解决任意折线区域内的单元生成。包括负斜率的直线。

!write(*,*) "区域是矩形区域吗？Y or N"

!read(*,*) bj_shape

bj_shape="Y"

c_num=0 !生成的单元总数

if (bj_shape=="Y") then

do i=1,c_layer,1

lay_cnum(i)=(bj_ptxmax-bj_ptxmin)/(2*ran_Rmean)

!开始生成第i行的单元。从右到左生成。

do j=1,lay_cnum(i),1

c_num=c_num+1

c_ox(c_num)=bj_ptxmax-((j-1)+0.5)*2*ran_Rmean ! x(j)=pointx(1)+(2*j-1)*Ra c_oy(c_num)=bj_ptymin+((i-1)+0.5)*2*ran_Rmean !即各层y的值分别为Ra,3Ra,5Ra,......

c_R(c_num)=ran_Rmean*0.999 !单元的半径全部相同

end do

call celltoCAD()

return

end subroutine

!---------------------------------------------------------------------! ! (1.3) 直接由文件 CADcell.txt 导入单元 ! !---------------------------------------------------------------------! subroutine CADcell() !由CAD图形导入产生单元

use global

implicit none

write(*,*) "从已知文件CADcell.txt中提取单元"

write(*,*) "总共有多少个单元："

read (*,*) c_num

open (14001,file='CADcell.txt')

do i=1,c_num,1

read (14001,*) c_ox(i),c_oy(i),c_R(i)

end do

close(14001)

call celltoCAD()

return

end subroutine

!---------------------------------------------------------------------! ! (2) 主子程序: 离散元参数的确定 ! !---------------------------------------------------------------------! subroutine cell_Information()

use global

implicit none

real c_massmin !单元质量mass最小值

character(len=79)::restraint_file

integer,parameter::file_id=20001

logical alive

integer restraint_num !需要设置外部约束的单元数

character readcellinfo_flag !直接读入数据标志

character agree_flag !单元的物理力学参数是否符合要求

open (20002,file='indata.txt')

write(*,*) "单元需施加外部约束吗?是(Y)或否(N):"

read(*,*) restraint_flag

if (restraint_flag=="Y") then !输入外部约束文件

20003 write(*,*) "输入施加外部约束的文件名(包括后缀名):" !读入外部约束文件

read (*,"(A79)") restraint_file

k=1 !对应restraint_file中数据行数

inquire(file=restraint_file, exist=alive)

if (alive) then

write(*,*) "输入需要施加外约束的单元个数:"

read(*,*) restraint_num

open(unit=file_id,file=restraint_file,access="sequential",status="old") !"old" Indicates an existing file do while (k

read(file_id,*) k,vx_ini(k),vy_ini(k)

k=k+1

end do

read(*,*)

else

write(*,*) TRIM(restraint_file)," doesn't exist." !返回字符restraint_file,但后面的空格全部取消 goto 20003

end if

end if

close(20002)

if(readcellinfo_flag=="N") then

write(*,*) "请输入单元的厚度c_thk(m,参考值0.01)"

read(*,*) c_thk

write(*,*) "请输入单元的密度c_dent(kg/m^3,参考值2150)"

read(*,*) c_dent

write(*,*) "输入单元的泊松比值c_niu(参考值0.25):"

read(*,*) c_niu

write(*,*) "请输入单元的弹性模量(Pa,参考值1000000000):"

read(*,*) c_E

write(*,*) "请输入单元之间的粘聚力c_cohe(Pa,参考值450000):"

read(*,*) c_cohe

write(*,*) "请输入单元之间摩擦系数(参考值0.51):" !tan(pi/6)=0.36

read(*,*) c_fric

write(*,*) "请输入单元与边界墙之间的粘聚力cw_coh(Pa,参考值160000):"

read(*,*) cw_coh

write(*,*) "请输入单元与边界墙之间的摩擦系数cw_fric(参考值0.45):" !tan(pi/6)

read(*,*) cw_fric

write(*,*) "请输入单元接触点的法向和切向刚度系数kn(N/m,参考值20000),ks(N/m,参考值10000):" read(*,*) kn,ks

write(*,*) "请输入单元与边界墙(刚性地面)之间的法向刚度系数cw_nstif(N/m,参考值150000)"

read(*,*) cw_nstif

write(*,*) "请输入单元与边界墙(刚性地面)之间的切向刚度系数cw_sstif(N/m,参考值100000):"

read(*,*) cw_sstif

write(*,*) "输入Rayleigh阻尼系数Ray_a(参考值0.3),Ray_b(参考值0.6):"

read(*,*) ray_a,ray_b

write(*,*) "单元触地是否反弹,是(Y),否(N):"

read(*,*) rebound_flag

write(*,*) "单元之间的刚度系数是一固定值?是(Y)或否(N)"

read(*,*) knksfix_flag

write(*,*) "是否考虑阻尼影响？是(Y)否(N)"

read(*,*) damp_flag

read(*,*) liju_flag

write(*,*) "是否考虑两竖向边界? 是(Y)否(N)"

read(*,*) verticalbj_flag

write(*,*) "确定时步系数的值(0.01):" !确定计算的时步。!先求出时步的推荐值。

read(*,*) deltat_yita

write(*,*) "请输入计算的总步数(推荐值100000):" !确定计算的步数。

read(*,*) js_stepsum

write(*,*) "请输入每计算多少步存储一次结果数据(推荐值100):"

read(*,*) save_step

write(*,*) "每经过多少时步查找单元的相交情况(推荐值1):"

read(*,*) seekJC_step

write(*,*) "设置CAD图形显示的格式,1:红色,2:黄色,3:绿色,4:青色,5:蓝色,6:品红,7:白色" read(*,*) fig_color

elseif(readcellinfo_flag=="Y") then

open(20005,file='cell_info.txt') !从cell_info.txt中读单元信息

read(20005,*) c_thk

read(20005,*) c_dent

read(20005,*) c_niu

read(20005,*) c_E

read(20005,*) c_cohe

read(20005,*) c_fric

read(20005,*) cw_coh

read(20005,*) cw_fric

read(20005,*) kn

read(20005,*) ks

read(20005,*) cw_nstif

read(20005,*) cw_sstif

read(20005,*) Ray_a

read(20005,*) Ray_b

read(20005,*) rebound_flag

read(20005,*) knksfix_flag

read(20005,*) damp_flag

read(20005,*) verticalbj_flag

read(20005,*) deltat_yita

read(20005,*) js_stepsum

read(20005,*) save_step

read(20005,*) seekJC_step

read(20005,*) fig_color

else

write(*,*) "请正确选择是(Y)或否(N)"

goto 20004

end if

close(20005)

write(*,*)"单元的厚度(m)",c_thk

write(*,*)"单元的密度(kg/m^3)",c_dent

write(*,*)"单元的泊松比值",c_niu

write(*,*)"单元的弹性模量(pa)",c_E

write(*,*)"单元之间的粘聚力(N)",c_cohe

write(*,*)"单元之间摩擦系数",c_fric

write(*,*)"单元与边界墙之间的粘聚力(N)",cw_coh

write(*,*)"单元与边界墙之间的摩擦系数",cw_fric

write(*,*)"单元接触点的法向和切向刚度系数",kn

write(*,*)"单元接触点的法向和切向刚度系数",ks

write(*,*)"单元与边界墙(刚性地面)之间的法向刚度系数",cw_nstif

write(*,*)"单元与边界墙(刚性地面)之间的切向刚度系数",cw_sstif

write(*,*)"Rayleigh阻尼常数值",Ray_a

write(*,*)"Rayleigh阻尼常数值",Ray_b

write(*,*)"单元触地是否反弹,是(Y),否(N):", rebound_flag

write(*,*)"单元之间的刚度系数是一固定值?是(Y)或否(N)", knksfix_flag write(*,*)"是否考虑阻尼影响?是(Y)否(N)", damp_flag

write(*,*)"是否考虑单元的力矩作用? 是(Y)否(N)", liju_flag

write(*,*)"是否考虑两竖向边界? 是(Y)否(N)",verticalbj_flag

write(*,*)"确定时步系数的值:",deltat_yita

FORTRAN程序设计复习题及答案

FORTRAN程序设计复习题 一、选择题 B (1)下列各FORTRAN表达式中合法的是 A) S+T*2P >= B) .NOT. (A*B+C) C) A2+B2/(C+D) <= D) (A+B).NOT．A*B．GT．(.NOT.只跟一个表达式) C (2)数学式(3/5)ex+y的FORTRAN表达式是 A) 3*EXP(X+Y)/5 B) 3*E* *(X+Y)/ C) (3/5)*EXP(X+Y)D) EXP(X+Y) D (3)下列FORTRAN77表达式中不合法的是 A) A．GT．B．EQV．C．GT．D B) A．AND．B．AND．C．AND．D C) .NOT．(X．LE．D) A．LT．B．LT．C．LT．D D(4)下列叙述中不正确的是 A) FORTRAN子程序可以单独编译 B) 对一个FORTRAN源程序进行编译和连接无误后可生成可执行文件 C) 即使编译和连接都正确无误，FORTRAN程序运行时仍可能出错 D) FORTRAN连接的主要任务是把函数库中的函数翻译成机器指令(正确描述：主要任务为连接目标文件) B (5)在下列FORTRAN77运算符中，优先级最高的是 A) .AND. B) .NOT. C) .OR. D) .EQ. B (6)FORTRAN表达式"6/5+9/2**3/2"的值为 A) 33 B) 1 C) 5 D) 3 A (7)下列FORTRAN77表达式中，合法的是： A) .AND.. B) 10.0 C) D) 提示：A）相当于 .AND.（.NOT.（）） D (8)关于编译一个FORTRAN源程序文件，下列说法中错误的是 A) 允许编译只有一个主程序而没有子程序的源文件 B) 允许编译有多个子程序的源文件 C) 允许编译只有一个子程序而没有主程序的源文件 D) 允许编译有多个主程序的源文件 C (9)在FORTRAN77源程序中，续行标志符必须放在 A) 第1列 B) 第1-6列C) 第6列D) 第5列 D (10)下列关于"SUBROUTIN E MAP(X，Y)"语句行的叙述中，不正确的是 A) 这是子程序的第一个语句 B) 字符串"MAP"是子程序名 C) 变量X是子程序的形参D) 子程序执行后，MAP将返回整型数据 提示：子程序无返回值，自定义函数才有） A (11)FORTRAN表达式"2/4+"的值是 A) B) 1 C) D) 0 提示：2/4默认等于整型，=》 D (12)FORTRAN表达式"MOD，"的值是 A) B)0.0 C) D) A (13下列FORTRAN运算符中，优先级最低的是 A）逻辑运算符.AND. B）算术运算符*

(完整)Fortran经典编程语言笔记(你值得拥有)

FORTRAN笔记 2014.10.07 目录 第七讲_FORTRAN的基本知识.ppt (2) FORTRAN语言程序设计初步 (2) FORTRAN源程序的书写格式（以77为例） (2) 变量 (2) 变量类型 (2) 算术运算符和运算优先级 (3) 赋值语句 (3) 参数语句（PARAMETER语句） (3) END语句 (3) PAUSE语句 (3) 逻辑运算和选择结构 (4) 关系表达式 (4) FORTRAN中数组的定义及使用 (4) 其他 (5) 1. fortran语言定义CHARACTER*6 TTL(14,3),CNAM(400)是什么意思? (5) 2. fortran里character*10 是什么意思 (5) 3. Fortran中kind是什么函数？ (5)

第七讲_FORTRAN的基本知识.ppt FORTRAN语言程序设计初步 FORTRAN是Formula Translation的缩写，意为“公式翻译”，它是为科学、工程问题或企事业管理中的那些能够用 数学公式表达的问题而设计的，其数值计算的功能较强。 常用的是FORTRAN77和FORTRAN90两种标准。 1、一个程序由若干个程序单位组成。主程序和每一个子程序分别是一个独立的程序单位。 2、每一个程序单位都是以“END”结束的。 3、一个程序单位包括若干行。 1）语句行。由一个FORTRAN语句组成。 2）非语句行，即注释行。 4、FORTRAN程序中的语句可以没有标号，也可以有标号，根据需要而定。标号的作用是标志一个语句以便被其 他语句引用。 5、一个程序单位中各类语句的位置是有一定规定的。 6、FORTRAN源程序必须按一定的格式书写。 FORTRAN源程序的书写格式（以77为例） 每一行有80列，分别如下： 1、第1-5列为标号区。一行中第一列为“C”或“*”，该行即被认为是注释行。 2、第6列为“续行标志区”，如果在一行的第6列上写一个非空格和非零的字符，则该行作为其上一行的续行。 3、第7-72列为语句区。 4、第73-80列，注释区。 变量 变量名：一个变量需要用一个名字（变量名）来识别。在同一个程序单位中不能用同一个变量名代表两个不同的变 量。 FORTRAN的变量名按以下规则选定： 1）第一个字符必须是字母，即变量名必须以字母开头； 2）在一个字母后面可以跟1-5为数字或字母。 如果选定的变量名超过6个字符，则只有前面6个字符有效。 注：在变量名中大写与小写字母是等价的。 变量类型 整型变量Integer、实型变量Real、双精度变量Double Precision、复型变量Complex、逻辑型变量Logical和字符型变量Character。 1、隐含约定（I-N规则） FORTRAN规定：在程序中的变量名，凡以字母I，J，K，L，M，N六个字母开头的，即认为该变量为整型变量。 在程序中，凡是变量名以字母I,J,K,L,M,N,i,j,k,l,m,n开头的变量被默认为整型变量，以其他字母开头的变量被 默认为实型变量。 2、用类型说明语句确定变量类型 1）INTEGER语句（整型说明语句） 2）REAL语句（实型说明语句） 3）DOUBLE PRECISION语句（双精度说明语句） 4）COMPLEX语句（复型说明语句） 5）LOGICAL语句（逻辑型说明语句）

MATLAB 与C C + + 、FORTRAN语言混合编程

MATLAB 与C/ C + + 、FORTRAN语言混合编程摘要:对MATLAB 与C/ C + + 和FORTRAN 语言进行混合编程的常用方法进行了介绍,分析了其实现方式和各自的利弊,并用实例对MEX 文件实现方式进行了较详细的论述. 关键词: MATLAB ; C/ C + + ; FORTRAN ; 混合编程 中图分类号: TP313 文献标识码: A 文章编号:16722948X(2004) 0620547205 1 混合编程的意义及其实现方式 1. 1 混合编程的意义 MATLAB 语言具有功能强大、开发效率高等诸 多优点, 已在工程实际中得到广泛应用, 但是与 FORTRAN、C/ C + + 等语言相比,其执行效率较低, 特别是当程序中含有大量循环语句(例如迭代计算) 时,MATLAB 就有些力不从心, 速度很慢, 而运用 FORTRAN 等擅长数值计算语言进行编程,其运行效 率高. 一方面,如果放弃MATLAB 强大功能和数量 众多的应用工具箱,无疑是资源的极大浪费. 另一方 面,针对工程实际,人们用FORTRAN、C/ C + + 语言 已编写了大量实用程序,如果将其重新改写成M 文 件移植到MATLAB 环境中,不仅要花费大量的时间 和精力,而且有时还降低了其运行效率. 可否将二者 优势互补呢? 混合编程就是其有效的解决途径. 1. 2 混合编程的实现 正是考虑到上面这些原由,MATLAB 系统提供 了其应用程序接口(Application Program Interface) 来 解决这些问题. API 主要包括3 部分:MEX 文件——— 外部程序调用接口,用来在MATLAB 环境下调用 FORTRAN、C/ C + + 语言编写的程序;MAT 文件应 用程序———数据输入输出接口,用于MATLAB 系统 与外部环境之间互传数据; 计算引擎函数库——— MATLAB 处于后台作为一个计算引擎,与其它应用 程序建立客户机/ 服务器关系,在其它应用程序中调 用[1 ,2 ] . 1. 2. 1 MEX 文件 MEX 文件是按照一定格式,用FORTRAN 或C/ C + + 语言编写的源程序,在MATLAB 下借助相应 的编译器,生成的动态链接函数的统称. 在Windows 操作系统下,是用MATLAB 附带的批处理mex. bat 来编译生成文件后缀名为. dll (Dynamic Link Li2 brary) 动态链接库文件,该文件可在MATLAB 环境 下,像命令函数一样直接运行和调用,使用起来极为 方便. 采取MEX 文件方式,是重复利用已有FOR2 TRAN、C/ C + + 程序,让MATLAB 和FORTRAN、

fortran语言语法

FORTRAN是世界上最早出现的高级编程语言，是工程界最常用的编程语言，它在科学计算中（如航空航天、地质勘探、天气预报和建筑工程等领域）发挥着极其重要的作用。经过40多年的发展，伴随着FORTRAN语言多次版本的更新及相应开发系统的出现，其功能不断完善，最新版本的开发系统几乎具备了VC、VB的所有特点，如图形界面编程、数据库等。目前，工科院校开设的计算机编程语言课首选仍然是FORTRAN :< 说实话，从科技发展的趋势来说这不是好事。您可以设想一下，如果需要用鹅毛笔抄写大量的古籍是什么感受！ 强烈建议阅读《发掘C#特性赋予科学计算项目以威力》 1 FORTRAN77四则运算符 + - * / ** (其中**表示乘方) 在表达式中按优先级次序由低到高为： +或-→*或/→**→函数→() 2 FORTRAN77变量类型 隐含约定：I-N规则 凡是以字母I，J，K，L，M，N六个字母开头的，即认为是整型变量，其它为实型变量。 用类型说明语句确定变量类型：可以改变I-N规则

用IMPLICIT语句将某一字母开头的全部变量指定为所需类型 如IMPLICIT REAL (I,J) 三种定义的优先级别由低到高顺序为：I-N规则→IMPLICIT语句→类型说明语句，因此，在程序中IMPLICIT语句应放在类型说明语句之前。 数组的说明与使用 使用I-N规则时用DIMENSION说明数组，也可在定义变量类型同时说明数组，说明格式为：数组名(下标下界,下标上界），也可省略下标下界，此时默认为1，例：DIMENSION IA(0:9),ND(80:99),W(3,2),NUM(-1:0),A(0:2,0:1,0:3)

C语言转换为fortran语言

C/C++采用的是缺省调用约定是STDCALL约定.在C程序中,可以在函数原型的声明中使用_stdcall关键字来指明过程采用STDCALL调用约定。 Fortran过程采用的缺省标识符是，全部大写的过程名加上“_”前缀和“@n”后缀。在C程序中保留标识符的大小写。编译程序会给采用STDCALL约定的过程标识符加上“_”前缀和“@n”后缀。 Fortran过程缺省的参数传递方式是引用方式是。对于下面这个Fortarn过程：SUBROUTINE ForSub(ivar,rvar) INTEGER ivar REAL rvar WRITE(*,*) ivar,rvar END 在C语言程序中应给出过程的函数原型及调用方式为： void main() { extern void__stdcall FORSUB(int*I,float*f); int iCV AR=1; float rCV AR=2.0; FORSUB(&iCV AR,&rCV AR); } 在C++中调用Fortan的过程，在声明函数原型时需要用extern“C”语句，以避免C++编译程序对标识符的修饰；并且C++也可以通过引用方式传递参数。对于上面的Fortran过程，C++程序应给出的函数原型及调用方法是： void main() { extern “C”{void__stdcall FORSUB(int*I,float*f);} int iCV AR=1; float rCV AR=2.0; FORSUB(&iCV AR,&rCV AR); } 另外，也可以在Fortran中用!MS$ATTRIBUTES编译伪指令来改变Fortran子过程的调用约定，以便于被其他语言的程序调用。在下面的例子中，过程ForSub具有C语言的调用约定。 SUBROUTINE ForSub(ivar,rvar) !MS$ATTRIBUTES C：：ForSub INTEGER ivar REAL rvar WRITE(*,*) ivar,rvar END 这样，这个过程使用的是C调用约定，并且参数传递方式也变为传值方式，过程的标识符变为全部小写且有_前缀而无后缀的方式。在C语言源程序中的函数原型及调用方法为：void main() { extern void FORSUB(int ivar,float rvar); int iVar=1;

FORTRAN 95 语法速查

FORTRAN 95 语法速查 ----------白云、李学哲、陈国新、贾波编著《FORTRAN95程序设计》读书笔记 目录：一、应用程序的创建与运行/FORTRAN 95所用的字符/ 变量类型及其声明，常量声明/表达式与运算符 二、输入与输出：表控、有格式 三、选择语句与结构：IF语句、CASE结构 四、DO循环结构 五、数组：数组的声明，数组的引用，数组的算术运算，数组的输入、输出，给数组赋初值， 动态数组，WHERE、FORALL语句 六、子程序：语句函数，内部子程序，调用子程序时的虚实结合：形参为数组、非定界数组、子 程序名、星号，递归子程序，外部子程序，纯子程序，逐元子程序 七、派生数据类型与结构体 八、指针与动态链表 九、文件：存取方式，基本操作语句，各类文件的读写操作 十、接口、模块 十一、公用区、存储关联、数据块子程序 十二、绘图：坐标系、设置图形颜色、创建图形程序/ 常用过程：设置线型、绘一像素点、设置当前位置、绘直线、绘弧线、绘矩形、绘多边形、绘制扇形（圆、椭圆）/ 文字信息的显示 附/录：标准函数与标准子例行程序 一、基础部份 1-1 FORTRAN 95 应用程序的创建与运行 创建或运行FORTRAN 95程序必须在Microsoft Developer Studio平台上进行。尽管程序文本及相关文件的编辑可以在任一文本编辑器上进行，然后再拷到Studio的文档窗口中。但最好还是一开始就进入Studio环境。创建FORTRAN 95 程序的步骤大致如下： 1）启动Microsoft Developer Studio 可以通过不同方式运行dfdev.exe程序以启动Microsoft Developer Studio [开始] \ Compaq Visual Fortran 6 \ Developer Studio \ dfdev.exe：或 ……\CVF66 \https://www.wendangku.net/doc/ed759878.html,\MSDEV98\dfdev.exe Microsoft Developer Studio的界面如下图所示： 文档窗口 工作空间窗口 输出窗口

C与C++与FORTRAN混合编程

C/C++/FORTRAN 混合编程 混合编程在软件编程中是经常遇到的问题，尤其是C/C++/FORTRAN的混合编程，本文主要说明以上三种语言混合编程中经常遇到的问题，同时，也说明了不同平台下混合编程应注意的问题。 混合语言编程要注意的问题主要体现在：函数调用和数据结构的存储。 1 Windows平台 函数：由于Fortran编程语言没有大小写之分，Windows平台下的混合语言编程要注意的主要是大小写的问题。考虑到编译器的差异，可以用下面的方式进行跨平台编程的函数声明。（ C/C++编译器使用Microsoft Visual C++ 6.0, Fortran编译器使用 Digital Visual Fortran 6.0）。 假设一个C的函数为 void cFunction(); 那么，只需要在它的头文件里面进行如下定义即可： #ifdef __cplusplus extern “C” void { #endif extern void __stdcall CFunction(); ＃define cFunction CFUNCTION #ifdef __cplusplus } #endif 这样，在Fortran或者C++的程序里面就可以直接调用了。 假设是一个Fortran函数SUBROUTINE FFUNCTION(); 那么，在C++头文件里进行如下的定义就可以了： #ifdef __cplusplus extern “C” void { #endif extern void __stdcall ffunction(); ＃define ffunction FFUNCTION #ifdef __cplusplus } #endif 这样，就可以在C++的程序里面直接调用。由于C编译器里面，没有定义__cplusplus这个环境变量，因此，C文件里面，也可以直接使用这个头文件。如果是一个C++函数，如： void cPlusplusFunction（）；和c函数一样，进行下面的定义即可： #ifdef __cplusplus extern “C” void { #endif extern void __stdcall cPlusplusFunction (); ＃define cPlusplusFunction CPLUSPLUSFUNCTION #ifdef __cplusplus }

计算机程序设计语言(FORTRAN语言)

计算机程序设计语言(FORTRAN语言) (总分：36.00，做题时间：90分钟) 一、 (总题数：36，分数：36.00) 1.编译程序能将高级语言编写的源程序转换成( )。 A．解释程序 B．汇编程序 C．映象程序 D．目标程序 （分数：1.00） A. B. C. D. √ 解析： 2.一个完整的FORTRAN源程序( )。 A．至少包括一个主程序 B．至少包括一个主程序和一个子程序 C．由一个主程序与一个以上的子程序组成 D．由一个主程序与一个子程序组成 （分数：1.00） A. √ B. C. D. 解析： 3.语句函数定义语句在程序内合法的位置是( )。 A．在程序块开头语句之后，END语句之前 B．在程序块开头语句之后，可执行语句之前 C．在说明语句之后，END语句之前 D．在说明语句之后，可执行语句之前 （分数：1.00） A. B. C. D. √ 解析： 4.下列关于函数子程序虚实参数的错误说法是( )。 A．可以没有形参数 B．虚实结合的数组长度可以不同 C．实参表与虚参表类型可以不同 D．函数名可以作为虚参

（分数：1.00） A. B. C. √ D. 解析： 5.下列叙述中正确的是( )。 A．FORTRAN程序块中，无名公用语句只能有一个B．FORTRAN子程序中，至少应有一个RETURN语句C．FORTRAN程序块中，最后一行必须是END语句D．FORTRAN程序块中，必须有变量说明语句 （分数：1.00） A. B. C. √ D. 解析： 6.运行下面的程序时得不到所需的结果，其主要原因是( )。INTEGER X(11) DATA X/9，8，7，6，5，4，3，2，1，0，-1/ DO 10 1=1，X(1) ，-1 WRITE(*，*)1.0/SQRT(25.0-REAL(X(1) )* * 2) 10 CONTINUE END A．没有给X(11) 赋初值 B．发生除以零的情况 C．发生负数开平方的情况 D．循环参数设置错误 （分数：1.00） A. B. C. D. √ 解析： 7.下列数据中，不符合FORTRAN常量表示法的是( )。 A．-25.6 B．2.758D3 C．'FOOT"=' D．TRUE （分数：1.00） A. B. C. D. √ 解析：

fortran语言内部函数

附录 FORTRAN 90标准函数 符号约定： ●I代表整型;R代表实型;C代表复型;CH代表字符型;S代表字符串;L代表逻辑型;A代表数组;P代表指针;T代表派生类型;AT为任意类型。 ●s:P表示s类型为P类型(任意kind值)。s:P(k)表示s类型为P类型(kind值=k)。 ●[…]表示可选参数。 ●*表示常用函数。 表1 数值和类型转换函数 函数名说明 ABS(x)*求x的绝对值∣x∣。x:I、R,结果类型同x; x:C,结果:R AIMAG(x)求x的实部。x:C,结果:R AINT(x[,kind])*对x取整,并转换为实数(kind)。x:R, kind:I,结果:R(kind) AMAX0(x1,x2,x3,…)*求x1,x2,x3,…中最大值。x I:I,结果:R AMIN0(x1,x2,x3,…)*求x1,x2,x3,…中最小值。x I:I,结果:R ANINT(x[,kind])*对x四舍五入取整,并转换为实数(kind)。x:R, kind:I,结果:R(kind) CEILING(x)*求大于等于x的最小整数。x:R,结果:I CMPLX(x[,y][,kind]))将参数转换为x、(x,0.0)或(x,y)。x:I、R、C, y:I、R,kind:I,结果:C(kind) CONJG(x)求x的共轭复数。x:C,结果:C DBLE(x)*将x转换为双精度实数。x:I、R、C,结果:R(8) DCMPLX(x[,y])将参数转换为x、(x,0.0)或(x,y)。x:I、R、C, y:I、R,结果:C(8) DFLOAT(x)将x转换为双精度实数。x:I,结果:R(8) DIM(x,y)*求x-y和0中最大值,即MAX(x-y,0)。x:I、R, y的类型同x,结果类型同x DPROD(x,y)求x和y的乘积,并转换为双精度实数。x:R, y:R,结果:R(8)

fortran语法手册

1F O R T R A N77四则运算符+ - */ ** (其中**表示乘方) 在表达式中按优先级次序由低到高为：+或-→*或/→**→函数→() 2 FORTRAN77变量类型 隐含约定：I-N规则 凡是以字母I，J，K，L，M，N六个字母开头的，即认为是整型变量，其它为实型变量。 如IMPLICIT REAL (I,J) 三种定义的优先级别由低到高顺序为：I-N规则→IMPLICIT语句→类型说明语句，因此，在程序中IMPLICIT语句应放在类型说明语句之前。 数组的说明与使用 使用I-N规则时用DIMENSION说明数组，也可在定义变量类型同时说明数组，说明格式为：数组名(下标下界,下标上界），也可省略下标下界，此时默认为1，例：DIMENSION IA(0:9),ND(80:99),W(3,2),NUM(-1:0),A(0:2,0:1,0:3) REAL IA(10),ND(80:99)使用隐含DO循环进行数组输入输出操作：例如 WRITE(*,10) ('I=',I,'A=',A(I),I=1,10,2) 10FORMAT(1X,5(A2,I2,1X,A2,I4)) 使用DATA语句给数组赋初值 变量表中可出现变量名,数组名,数组元素名,隐含DO循环，但不许出现任何形式的表达式：例如 DATA A,B,C/,, DATA A/,B/,C/ DATA A,B,C/3*CHARACTER*6 CHN(10) DATA CHN/10*''/INTEGER NUM(1000) DATA (NUM(I),I=1,500)/500*0/,(NUM(I),I=501,1000)/500*1/ 3 FORTRAN77程序书写规则

[转载]Fortran 77, C, C++ 和 Fortran 90 的比较

[转载]Fortran 77, C, C++ 和 Fortran 90 的比较 收藏 发信人: quasar (飞贼克斯), 信区: Fortran 标题: Fortran 77, C, C++ 和 Fortran 90 的比较(转载) 发信站: 南京大学小百合站 (Tue Jun 1 10:59:14 2004) 瀚海星云 -- 文章阅读 [讨论区: MathTools] 发信人: HuiCai (老灰菜), 信区: SciComp 标题: Fortran 77, C, C++ 和 Fortran 90 的比较(转载) 发信站: 瀚海星云 (2002年12月19日10:40:38 星期四), 站内信件 【以下文字转载自 Fortran 讨论区】 【原文由 HuiCai 所发表】 Fortran 77, C, C++ 和 Fortran 90 的比较 https://www.wendangku.net/doc/ed759878.html,/develop/article/16/16085.shtm 三十年来, 从 Fortran 77 开始, Fortran 成为了计算科学的主要语言.在这段时间里, Fortran 的数值能力变得非常稳定而且优于其它计算机语言; 最大的改变来自于不断增长的各种可靠的数值过程库的种类. Fortran 联合(union), 它的使用技巧, 扩充的数值库为计算科学赋予了良好的基础. 可是在过去十几年中, 动态数据结构(特别是动态数组)的重要性不窜上升, UNIX 工作站, 复杂的交互式可视化工具, 以及更近的并行体系结构--Fortran 77 都没有实现--刺激了其它语言作为计算语言的使用, 最明显的一个例子是C. 最近C++ 也已经引起人们的兴趣, Fortran 通过发展到 Fortran 90来弥补它在现代科学计算 方面的不足. 这部分的一个通常的工作是比较四种语言对科学计算的适应性的, 这四种语言是两个C 的代表(C, C++) 和两个Fortran的代表(Fortran 77, Fortran 90). 下面的表格总结了这种比较, 后面的内容试图合 理地解释这种等级排序, 从最好(1)到最差(4).. 功能 ------------ F77 - C - C++ - F90 数值健壮性 ---- 2 ---- 4 --- 3 ----- 1 数据并行性 ---- 3 ---- 3 --- 3 ----- 1 数据抽象 ------- 4 ---- 3 --- 2 ----- 1 面向对象编程 - 4 ---- 3 --- 1 ----- 2 函数型编程 ---- 4 ---- 3 --- 2 ----- 1 平均等级 ------ 3.4 - 3.2 - 2.2 -- 1.2

fortran语法手册

1 FORTRAN77四则运算符 + - * / ** (其中**表示乘方) 在表达式中按优先级次序由低到高为： +或-→*或/→**→函数→() 2 FORTRAN77变量类型 2.1 隐含约定：I-N规则 凡是以字母I，J，K，L，M，N六个字母开头的，即认为是整型变量，其它为实型变量。 如 IMPLICIT REAL (I,J) 三种定义的优先级别由低到高顺序为：I-N规则→IMPLICIT语句→类型说明语句，因此，在程序中IMPLICIT语句应放在类型说明语句之前。 2.4 数组的说明与使用 使用I-N规则时用DIMENSION说明数组，也可在定义变量类型同时说明数组，说明格式为：数组名(下标下界,下标上界），也可省略下标下界，此时默认为1，例： DIMENSION IA(0:9),ND(80:99),W(3,2),NUM(-1:0),A(0:2,0:1,0:3) REAL IA(10),ND(80:99)使用隐含DO循环进行数组输入输出操作：例如 WRITE(*,10) ('I=',I,'A=',A(I),I=1,10,2) 10FORMAT(1X,5(A2,I2,1X,A2,I4)) 2.5 使用DATA语句给数组赋初值 变量表中可出现变量名,数组名,数组元素名,隐含DO循环，但不许出现任何形式的表达式：例如 DATA A,B,C/-1.0,-1.0,-1.0/ DATA A/-1.0/,B/-1.0/,C/-1.0/ DATA A,B,C/3*-1.0/CHARACTER*6 CHN(10)

DATA CHN/10*' '/INTEGER NUM(1000) DATA (NUM(I),I=1,500)/500*0/,(NUM(I),I=501,1000)/500*1/ 3 FORTRAN77程序书写规则 程序中的变量名，不分大小写； 变量名称是以字母开头再加上1到5位字母或数字构成，即变更名字串中只有前6位有效； 一行只能写一个语句； 程序的第一个语句固定为PROGRAM 程序名称字符串 某行的第１个字符至第５个字符位为标号区，只能书写语句标号或空着或注释内容； 某行的第１个字符为Ｃ或＊号时，则表示该行为注释行，其后面的内容为注释内容； 某行的第６个字符位为非空格和非０字符时，则该行为上一行的续行，一个语句最多可有19个续行； 某行的第７至７２字符位为语句区，语句区内可以任加空格以求美观； 某行的第７３至８０字符位为注释区，８０字符位以后不能有内容。 4 FORTRAN77关系运算符 .GT. 大于 .GE. 天于或等于 .LT. 小于 .LE. 小于或等于 .EQ. 等于 .NE. 不等于 .AND. 逻辑与 .OR. 逻辑或 .NOT. 逻辑非 .EQV. 逻辑等 .NEQV. 逻辑不等 运算符优先级由高到低顺序为:()→**→*或/→+或-→.GT.或.GE.或.LT.或.LE.或.EQ.或.NE.→.NOT.→.AND.→.OR.→.EQV.或.NEQV 5 FORTRAN77语句

Simple算法_fortran语言编写

Cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc Cccc This computer program was copied from the graduate student course program Cccc of the University of Minnesota. Part of it was re-formulated to meet the Cccc personal computer environment. Some inappropriate expressions were also Cccc corrected. The program is used only for teaching purpose. No part of it Cccc may be published. You may use it as a frame to re-develop your own code Cccc for research purpose. XJTU Instructor, 1995.11 **************************************************************************** *----------------------------MAIN PROGRAM----------------------------------- **************************************************************************** LOGICAL LSTOP COMMON/CNTL/LSTOP **************************************************************************** OPEN(08,FILE='teresul') CALL SETUP0 CALL GRID CALL SETUP1 CALL START 10 CALL DENSE CALL BOUND CALL OUTPUT IF(.NOT.LSTOP) GO TO 15 CLOSE(08) STOP 15 CALL SETUP2 GO TO 10 END *--------------------------------------------------------------------------- SUBROUTINE DIFLOW **************************************************************************** COMMON/COEF/FLOW,DIFF,ACOF **************************************************************************** ACOF=DIFF IF(FLOW .EQ.0.0)RETURN TEMP=DIFF-ABS(FLOW)*0.1 ACOF=0. IF(TEMP .LE. 0. ) RETURN TEMP=TEMP/DIFF ACOF=DIFF*TEMP**5 RETURN END *-------------------------------------------------------------------------- SUBROUTINE SOLVE ****************************************************************************

C与FORTRAN的差别

https://www.wendangku.net/doc/ed759878.html,/c.htm 2. C语言和Fortran语言的差异 由于两者产生的背景不同，它们是存在差异的，在比较了几组源代码之后，主要有以下体会： C 最大的优点在于灵活，不但可以藉由 struct 来定义新的数据结构，同时 C 的pointer 更可以让我们自由而且有效率地处理大数据。而在 UNIX 系统中，由于整个操作系统绝大部分就是 C 写出来的，故我们也有方便的 C 函数库，直接使用系统资源与享受系统带来的服务，以做到一些低阶、快速的动作。而FORTRAN从一开始就用于科学计算，它与C的差异主要表现为： * 复数运算的速度 * 程序参数与字串 * 内存的动态管理 * 多维阵列的处理 * 函数调用与参数传递 2.1. 复数运算的速度 在进行复数运算的时候，C++ 可以定义复数的 class，还可以重新定义所有的四则运算式，复杂的算式也可以做到由一个表达式来解决。但它的重新定义复数四则运算是用函数来做的，使用函数来调用其速度很慢，除非采用 inline function 的方式，但会遇到以

下的问题：要先将这个算式拆解，分别算过后再重组结果，故表面上程序代码很简洁，但实际上是 compiler做了很多工作，还是要付出相当的计算时间代价的。 至于 Fortran，最大的优点在于复数 (complex number) 的运算，复数是 Fortran 的基本数据类型之一，这正是 C 所缺乏的 (C 基本上只有实型与整型类型而已)。虽然C 也可以由 struct 的定义，达到复数四则运算的目的，但却很可能牺牲了程序效能，或者是程序写起来相当繁杂降低可读性。因此，在大量而且要求高速的复数运算场合， Fortran 实际上比 C 还要适合。 然而既然复数已是 Fortran 基本数据类型之一，则 Fortran compiler在设计上可以做到对复数特别的 optimization，例如如果遇到较短的复数运算式，它可以用“心算”直接得出 real_part 与 imag_part 的 expression，像这样: real(a) =……；imag(a) = ……. 如此只需两步就得到结果。直到遇到太长太复杂的算式，才去做拆解的动作。 这样使用 C 来做复数运算可能需要绕圈圈，而且绕出来的圈圈可能还不小。不过如果程序中需要复合的数据结构，如一个自定义的数据结构中既有浮点数、整数、还有字符串时， Fortran 只有举白旗投降了。当然， Fortran 如果要做还是可以做，只是不太方便，而且可能也需要绕圈圈。但如果使用 Fortran 90 则不成问题了，因为 Fortran 90 也有类似 C 的 struct 结构以定义复合的数据类型。 2.2. 程序参数与字串 C 程序可以有参数串列， Fortran 则没有。例如，当程序执行时，必须输入 a, b, c 三个参数，在 C 可以这样写: int main(int argc, char **argv) { int a, b, c;

Fortran语言

#LOCAL REAL Uds,Uqs,angles,vs,w,o,ang,ta,tb,t0,Ud1,Uq1,E #LOCAL REAL tc,t,mt,pi #LOCAL INTEGER sign,M,nl,n t=TIME tc=$TC nl=DNINT(tc/delt) pi=3.1415926 E=$Udc w=2*pi*50 if(t.lt.delt) then sign=0 M=0 Uds=0 Uqs=0 angles=0 vs=0 o=0 ang=0 ta=0 tb=0 t0=0 Ud1=0 Uq1=0 n=1 end if if(mod(n,nl).lt.1E-7) then !初始化变量 M=M+1 Uds=$Ud Uqs=$Uq angles=$angle !变换到alphabeta坐标直角坐标

Ud1=cos(angles)*Uds+sin(angles)*Uqs Uq1=sin(angles)*Uds-cos(angles)*Uqs ! Ud1=$Ud ! Uq1=$Uq !生成矢量的模和角度 vs=sqrt(Ud1**2+Uq1**2) if(Ud1.gt.0.and.Uq1.ge.0) then ang=atan(Uq1/Ud1) else ang=atan(Uq1/Ud1)+2*pi end if if(Ud1.lt.0) then ang=atan(Uq1/Ud1)+pi end if if(Ud1.eq.0) then if(Uq1.ge.0) then ang=pi/2 else ang=3*pi/2 end if end if !判断矢量的所属,计算角度 if(ang.ge.0.and.ang.lt.pi/3) then sign=1 o=ang ! o=w*t end if if(ang.ge.pi/3.and.ang.lt.2*pi/3) then sign=2 o=ang-pi/3 ! o=w*t-pi/3 end if if(ang.ge.2*pi/3.and.ang.lt.3*pi/3) then sign=3 o=ang-2*pi/3 ! o=w*t-2*pi/3

fortran初步学习资料

第一章: Fortran语言程序设计初步 Fortran语言的发展概况 本节介绍Fortran的起源与发展历史，讲述Fortran由产生到形成标准FortranIV、Fortran77，并进一步形成新标准Fortran90/95的发展历程。 1.1.1 Fortran的历史 a)a)FortranI FortranIV Fortran是目前国际上广泛流行的一种高级语言，适用于科学计算。Fortran 是英文FORmula TRANslatio n的缩写，意为“公式翻译”。它是为科学、工程问题中的那些能够用数学公式表达的问题而设计的语言，主要用于数值计算。这种语言简单易学，因为可以像抄写数学教科书里的公式一样书写数学公式，它比英文书写的自然语言更接近数学语言。Fortran语言是第一个真正推广的高级语言。至今它已有四十多年历史，但仍历久不衰，始终是数值计算领域所使用的主要语言。Fortran语言问世以来，根据需要几经发展，先后推出形成了很多版本。 第一代Fortran语言是在1954年提出来的，称为FortranI。它于1956年在IBM 704计算机上得以实现。在此之前编写计算机程序是极为繁琐的，程序员需要详细了解为之编写代码的计算机的指令、寄存器和中央处理器(CPU)等方面的知识。源程序本身是用数学符号(八进制码)编写的，后来采用了助记符，即所谓机器码或汇编码，这些编码由汇编程序转换为指令字。在50年代书写和调试一个程序要很长时间，因为用这种方式编写程序显然是很不方便的，尽管它能使CPU高效地工作。正是这些原因，促使由John Backus率领的IBM公司的一个小组研究开发最早的高级程序设计语言Fortran。其目的是开发一种容易理解、简单易学又能几乎像汇编一样高效运行的语言，他们取得了极大的成功。Fortran 语言作为第一种高级语言不仅是一次创新，也是一次革命。它使程序员摆脱了使用汇编语言的冗长乏味的负担，而且它使得不再只是计算机专家才能编写计算机程序，任何一名科学家或工程技术人员，只要稍加努力学习和使用Fortran，就能按自己的意图编写出用于科学计算的程序。 经过不断发展，FortranI形成了很多不同版本，其中最为流行的是1958年出现的FortranII，它对FortranI进行了很多扩充(如引进了子程序)，FortranII 在很多机器上得以实现。其后出现的FortranIII未在任何计算机上实现。1962年出现的FortranIV对原来的Fortran作了一些改变，使得FortranII源程序在FortranIV编译程序下不能全部直接使用，导致了语言不兼容的问题。这样就形成了当时同时使用FortranII和FortranIV两种程序设计语言的局面。

fortran与c语言接口参数传递混合编程

Sun Studio 12：Fortran 编程指南 ?Previous: 第10 章并行化 第11 章C-Fortran 接口 本章论述Fortran 与C 的互操作性方面的问题，内容仅适用于Sun Studio Fortran 95 和C 编译器的特定情况。 11.9 Fortran 2003 与C 的互操作性简要说明了Fortran 2003 标准第15 部分中提到的C 绑定功能。（此标准可以从国际Fortran 标准Web 站点获得）。Fortran 95 编译器实现了标准中所述的这些功能。 如不特别注明，32 位x86 处理器视为与32 位SPARC 处理器等同。对于64 位 x86 处理器和64 位SPARC 处理器也是如此，只是x86 系统未定义REAL*16 和COMPLEX*32 数据类型，这些数据类型只能用于SPARC。 11.1 兼容性问题 大多数 C-Fortran 接口必须在以下这些方面全部保持一致： ?函数和子例程的定义及调用 ?数据类型的兼容性 ?参数传递（按引用或按值） ?参数的顺序 ?过程名（大写、小写或带有结尾下划线(_)） ?向链接程序传递正确的库引用 某些C-Fortran 接口还必须符合： ?数组索引及顺序 ?文件描述符和stdio ?文件权限 11.1.1 函数还是子例程？ 函数一词在 C 和Fortran 中有不同的含义。根据具体情况做出选择很重要：?在C 中，所有的子程序都是函数；但void函数不会返回值。 ?在Fortran 中，函数会传递一个返回值，但子例程一般不传递返回值。

当Fortran 例程调用C 函数时： ?如果被调用的C 函数返回一个值，则将其作为函数从Fortran 中调用。 ?如果被调用的C 函数不返回值，则将其作为子例程调用。 当C 函数调用Fortran 子程序时： ?如果被调用的Fortran 子程序是一个函数，则将其作为一个返回兼容数据类型的函数从 C 中调用。 ?如果被调用的Fortran 子程序是一个子例程，则将其作为一个返回int（与Fortran INTEGER*4兼容）或void值的函数从C 中调用。如果Fortran 子例程使用交替返回，则会返回一个值，这种情况下它是RETURN语句中的表达式的值。如果RETURN语句中没有出现表达式，但在SUBROUTINE语句中声明了交替返回，则会返回零。 11.1.2 数据类型的兼容性 表11–2总结了Fortran 95（与 C 比较）数据类型的数据大小和缺省对齐。该表假设未应用影响对齐或提升缺省数据大小的编译选项。请注意以下事项： ? C 数据类型int、long int和long在32 位环境下是等同的（4 字节）。 但是，在64 位环境下long和指针为8 字节。这称为LP64 数据模型。 ?在64 位SPARC 环境下，当用任意-m64选项进行编译时， REAL*16和COMPLEX*32与16 字节边界对齐。 ?标有4/8 的对齐表示缺省情况下与8 字节边界对齐，但在COMMON 块中与 4 字节边界对齐。COMMON 中的最大缺省对齐为4 字节。当用-m64选项进 行编译时，4/8/16 表示与16 字节边界对齐。 ?REAL(KIND=16)、REAL*16、COMPLEX(KIND=16)、COMPLEX*32只能用于SPARC 平台。 ?数组和结构的元素及字段必须兼容。 ?不能按值传递数组、字符串或结构。 ?可以在调用点使用%VAL(arg)，按值将参数从Fortran 95 例程传递到C 例程。假如Fortran 例程具有一个显式接口块，该接口块用VALUE属性声明了伪参数，则可以按值将参数从 C 传递到Fortran 95。 ?数值序列类型的组件的对齐方式与通用块的对齐方式相同，也会受到-aligncommon选项的影响。数值序列类型是这样一种序列类型：其中所 有组件的类型为缺省整数、缺省实数、双精度实数、缺省复数或缺省逻辑，而不是指针。 ?在大多数情况下，非数值序列类型的数据类型组件以自然对齐的方式对齐，但QUAD 变量除外。对于四精度变量，32 位SPARC 平台和64 位SPARC 平台之间的对齐方式不同。 ?在所有平台上，用BIND(C) 属性定义的VAX 结构和数据类型的组件始终与Ｃ结构具有相同的对齐方式。