计算程序

计算程序,第1张

c fengfor五层准三维流(潜/弱/承/弱/承,对弱透水层只考虑越流不考虑d性释放),调参用

c 使潜、承、承上下节点、上下单元完全重合,且总单元数ms=3×ms1,总节点数ids=3×ids1

c 潜单元编号由1-ms1,中承单元编号由ms1+1-ms2,深承单元编号由ms2+1-ms

c 10-13行:可调数组的各个参数预先赋值,不同计算区只要改变这些参数即可用此程序,其中潜单元个数=ms1个,中承单元个数=ms2-ms1个,深承单元个数=ms-ms2个;潜节点编号由1-ids1,中承节点编号由ids1+1ids2,深承节点编号由ids2+1-ids,n:未知水头节点数;mi:开采节点数;igs:拟合点总数;igs1:潜拟合点数,igs2igs1:中承拟合点数,icqs:参数区个数:ihv:计算时段数,jbn:纯外引地表水灌溉节点个数;nn1:潜已知水头节点个数;nn:潜+中+深已知水头节点总数

parameter(ms1=108,ids1=67,n1=59,mi1=17,igs1=11,icqs1=8)

parameter(ms2=216,ids2=134,n2=126,mi2=47,igs2=25,icqs2=16)

parameter(ms=324,ids=201,n=193,mi=69,igs=35,icqs=24)

parameter(ihv=24,jbn=9,nn1=8,nn=24)

c h0:时段初水头(m);hed:时段末水头(m);in:各单元三节点编号(必须按小中大顺序排)

c fh:拟合点各时段计算水位(m);sh:拟合点各时段实测水位(m)

dimension h0(ids),hed(ids),in(ms,3),fh(igs,ihv),sh(igs,ihv)

c idyh:导水矩阵工作单元;zb:各节点x,y坐标(输入坐标为图面mm数);igdh:拟合点节点号;q:开采井各时段水量(抽为正,注为负,流出边界为正,流入为负)(m3/d)

dimension idyh(ids,9),zb(ids,2),igdh(igs),q(mi,ihv)

c d:各节点导水矩阵;ifdh:开采井节点号;s:拟合点拟合误差(m);mqh:各单元的参数区号;e:各节点释(储)水矩阵dimension ifdh(mi),d(ids,9),s(igs,ihv),mqh(ms),e(ids)

c fqc:各参数区最多4个参数之值;sss:各节点水头总变幅,由t0时刻起算,(m);h:初始流场,gh潜水层独有的灌溉入渗系数

dimension fqc(icqs,4),sss(ids),h(ids),gh(icqs1)

c yo:各节点越流矩阵工作单元;jbiao:纯引地表水灌溉的节点编号,共有jbn个节点;qxia:灌溉抽取地下水均铺在灌溉区的水层厚度;gq:灌溉抽出的地下水水层厚度;gq=qxia;w:各时段潜水各节点的降水灌溉回归水量

dimension yo(ids),jbiao(jbn),qxia(ihv),gq(ihv),w(ids1,ihv)

c zzz:潜水各节点底板标高(m);ep0:潜水各节点含水层厚(m);ep1:潜水各单元含水层厚均值(m);x:各时段降水量(m/时段)

dimension zzz(ids1),ep0(ids1),ep1(ms1),x(ihv)

c tl:各时段步长值(天);bc:三角单元几何量(bi,ci,bj,cj,bk,ck)

dimension tl(ihv),bc(3,2)

c 分别是潜(h0nn1),中(h0nn2),深(h0nn3)层已知变水头节点各时段的水头

dimension h0nn1(nn1,ihv),h0nn2(nn1,ihv),h0nn3(nn1,ihv)

c file1(igs)8:定义8字符的字符串igs个;filee8:定义8字符的文件名

character file1(igs)8,file2(igs)8,filee8

write(,23)

23 format(1x,’zz1=?zz2=?迭代因子:zz1潜水0—1,zz2承压水1—2’)

c zz1:潜水为亚松弛系数,一般取085为好;zz2:承压水为超松弛系数,当ids<300时取12-13,当ids>300时取13-15为好

zz1=085

zz2=125

write(,24)ms1,ids1,n1,mi1,igs1,icqs1

write(,24)ms2,ids2,n2,mi2,igs2,icqs2,ihv

write(,24)ms,ids,n,mi,igs,icqs

24 format(5x,12i5)

open(1,file=’fqc’,status=’old’)

c 首先从’fqc’中读取潜水的区号m10(空读),参数1(K),参数2(μ),参数3潜-中的越流系数(k’/m’),参数4降水入渗系数(α),gh灌溉入渗系数

read(1,)(m10,fqc(i,1),fqc(i,2),fqc(i,3),fqc(i,4),gh(i),

i=1,icqs1)

c 接着从’fqc’中读取中层水的区号m10(空读),参数1(T),参数2(μ)

read(1,)(m10,fqc(i,1),fqc(i,2),i=icqs1+1,icqs2)

c 接着从’fqc’中读取深层水的区号m10(空读),参数1(T),参数2(μ),参数3:中-深的越流系数(k’/m’)

read(1,)(m10,fqc(i,1),fqc(i,2),fqc(i,3),i=icqs2+1,icqs)

close(1)

c 参数3:中-深的越流系数(k’/m’)同样用于中层水

do 54 i=1,icqs1

fqc(icqs1+i,3)=fqc(i,3)

54 continue

c qxia:灌溉抽取地下水平均铺在灌溉区的水层厚度

open(1,file=’gq’,status=’old’)

read(1,)(qxia(i),i=1,ihv)

close(1)

c gq:灌溉抽出的地下水水层厚度

open(1,file=’gq’,status=’old’)

read(1,)(gq(i),i=1,ihv)

close(1)

c jbiao:纯引地表水灌溉的节点编号,共有jbn个节点

open(1,file=’jbiao’,status=’old’)

read(1,)(jbiao(i),i=1,jbn)

close(1)

c’in’:各单元三节点编号文件(必须按小中大顺序排),ia为单元号空读

open(1,file=’in’,status=’old’)

read(1,)(ia,(in(i,j),j=1,3),i=1,ms)

close(1)

c igdh:拟合点节点编号

open(1,file=’igdh’,status=’old’)

read(1,)(igdh(i),i=1,igs)

close(1)

c mqh:各单元的参数区号,ia为单元号空读

open(1,file=’mqh’,status=’old’)

read(1,)(ia,mqh(i),i=1,ms)

close(1)

c ifdh:开采井节点号

open(1,file=’ifdh’,status=’old’)

read(1,)(ifdh(i),i=1,mi)

close(1)

c tl(i):各时段值(天/时段,每月为1个时段)

open(1,file=’tl1’,status=’old’)

read(1,)(tl(i),i=1,ihv)

close(1)

c 每个时段的大气降水量(m),x(i)/dt则换算为(m/d)

open(1,file=’x’,status=’old’)

read(1,)(x(i),i=1,ihv)

close(1)

c’zb’:各节点x,y坐标文件(此处为1:25万图面的mm数),ia为节点号空读

open(1,file=’zb’,status=’old’)

read(1,)(ia,(zb(i,j),j=1,2),i=1,ids)

close(1)

c 把各节点x,y坐标1:25万图面的mm数换算为实地的m数

do 888 i=1,ids

do 888 j=1,2

zb(i,j)=zb(i,j)25

888 continue

c 潜水开采井各时段的开采量(m3/d),空读ia开采井个数,空读ib开采井所在节点编号

open(1,file=’qq’,status=’old’)

read(1,)ia,((ib,q(i,j),i=1,mi1),j=1,ihv)

close(1)

c 中层水开采井第一时段的开采量(m3/d),空读ia开采井个数,空读ib开采井所在节点编号

open(1,file=’qz1’,status=’old’)

read(1,)ia,(ib,q(i,1),i=mi1+1,mi2)

close(1)

c 中层水开采井第13时段的开采量(m3/d),空读ia开采井个数,空读ib开采井所在节点编号

open(1,file=’qz2’,status=’old’)

read(1,)ia,(ib,q(i,13),i=mi1+1,mi2)

close(1)

c 深层水开采井第一时段的开采量(m3/d),空读ia开采井个数,空读ib开采井所在节点编号

open(1,file=’qs1’,status=’old’)

read(1,)ia,(ib,q(i,1),i=mi2+1,mi)

close(1)

c 深层水开采井第13时段的开采量(m3/d),空读ia开采井个数,空读ib开采井所在节点编号

open(1,file=’qs2’,status=’old’)

read(1,)ia,(ib,q(i,13),i=mi2+1,mi)

close(1)

c 中、深层水开采井第2-12时段的开采量等于第一时段的开采量

do 887 j=2,12

do 886 i=mi1+1,mi

q(i,j)=q(i,1)

886 continue

887 continue

c 中、深层水开采井第14-ihv时段的开采量等于第13时段的开采量

do 884 j=14,ihv

do 883 i=mi1+1,mi

q(i,j)=q(i,13)

883 continue

884 continue

c sh:拟合点各时段实测水位(m),空读ia拟合点所在节点编号

open(3,file=’sh’,status=’old’)

read(3,)(ia,(sh(i,j),j=1,ihv),i=1,igs)

close(3)

c qc:潜水初始流场,空读is节点编号

open(1,file=’qc’,status=’old’)

read(1,)(is,h(i),i=1,ids1)

close(1)

c zc:中层水初始流场,空读is节点编号

open(1,file=’zc’,status=’old’)

read(1,)(is,h(i),i=ids1+1,ids2)

close(1)

c sc:深层水初始流场,空读is节点编号

open(1,file=’sc’,status=’old’)

read(1,)(is,h(i),i=ids2+1,ids)

close(1)

c h0nn:分别是潜(h0nn1),中(h0nn2),深(h0nn3)层已知变水头节点各时段的水头,空读ii节点编号

open(1,file=’h0nn’)

read(1,)(ii,(h0nn1(i,ikv),ikv=1,ihv),i=1,nn1)

read(1,)(ii,(h0nn2(i,ikv),ikv=1,ihv),i=1,nn1)

read(1,)(ii,(h0nn3(i,ikv),ikv=1,ihv),i=1,nn1)

close(1)

c zzz:潜水各节点底板标高(m),空读ii节点编号

open(1,file=’zzz’)

read(1,)(ii,zzz(i),i=1,ids1)

close(1)

c’file1’存放igs个字符串(拟合节点名称占4字符,及观测井的原编号再占4字符)

open(1,file=’file1’,status=’old’)

read(1,110)(file1(i),i=1,igs)

close(1)

c’file2’存放igs个字符串(拟合节点名称占4字符,加dat后缀再占4字符)

open(1,file=’file2’,status=’old’)

read(1,110)(file2(i),i=1,igs)

close(1)

110 format(10a8)

c 计算开始前先把全部节点的时段末刻水头hed(i),时段初刻水头h0(i)赋初值h(i),以使开始迭代时hed(i),h0(i)不为零

do 1993 i=1,ids

hed(i)=h(i)

1993 h0(i)=h(i)

c 对中、深层承压水导水矩阵工作单元idyh,释水矩阵e,越流矩阵yo,导水矩阵d,赋0

do 25 i=1+ids1,ids

do 25 j=1,9

idyh(i,j)=0

e(i)=00

yo(i)=00

25 d(i,j)=00

c sum2用来累计计算区总面积(用承压水总面积代替),先赋零

sum2=0

c 对承压水逐个单元计算几何量及导水、释水、越流矩阵

c 对承压水第ip单元的三个节点号依次赋给i,j,k及i1,j1,k1

do 80 ip=ms1+1,ms

i=in(ip,1)

j=in(ip,2)

k=in(ip,3)

i1=i

j1=j

k1=k

c idyh(i1,9)存放与i1点同单元的所有节点号,最多9个,可以用不完,即i1点的idyh可以有几个idyh(i1,i2)=0;当计算第ip单元时,i1点的idyh由i1占第1个(i2=1)位置,j1,k1只能占i2=2,3,…,9 中的位置;且先占者排前,193行:使j1,k1分两轮到idyh中找位置;当j1找时,195行发现i1的idyh中已有j1,则跳到j1对i1的96循环体头,换为k1对i1循环;当196行没发现已有j1,且196行判断此位不空时返回96循环体头找下一个位置,当碰到第1个空位时,由j1占据,然后跳到j1对i1循环体头,换为k1对i1循环

do 90 iv=1,3

idyh(i1,1)=i1

do 94 iu=j1,k1,k1-j1

do 96 i2=2,9

if(idyh(i1,i2)eqiu)goto 94

if(idyh(i1,i2)ne0)goto 96

if(idyh(i1,i2)eq0)idyh(i1,i2)=iu

goto 94

96 continue

94 continue

c 把ip单元的下一个节点j提为i1,k提为j1,i降为k1,然后返回90循环头处理ip单元的下一个节点,循环3次则ip单元中i,j,k 3个点的idyh全部找完

iuu=i1

i1=j1

j1=k1

k1=iuu

90 continue

c 第ip单元中i,j,k三节点的x坐标赋给xi,xj,xk,y坐标赋给yi,yj,yk

xi=zb(i,1)

xj=zb(j,1)

xk=zb(k,1)

yi=zb(i,2)

yj=zb(j,2)

yk=zb(k,2)

c 第ip单元三角形面积ss

ss=abs((xj-xi)(yk-yi)-(xk-xi)(yj-yi))05

c 累加各单元面积(这里的sum2累加之后是中、深层两层面积之和)

sum2=sum2+ss

c 第ip单元的bi~bc(1,1),ci~bc(1,2),bj~b(2,1),cj~b(2,2),bk~b(3,1),ck~b(3,2)

bc(1,1)=yj-yk

bc(1,2)=xj-xk

bc(2,1)=yk-yi

bc(2,2)=xk-xi

bc(3,1)=yi-yj

bc(3,2)=xi-xj

c 第ip单元所在参数区号赋给jv,参数1(T)赋给txy,参数2(μ)赋给ts,参数3(k’/m’)赋给rmk

jv=mqh(ip)

txy=fqc(jv,1)

ts=fqc(jv,2)

rmk=fqc(jv,3)

c 第ip单元三节点i,j,k的释水矩阵元素c(ii,ii)及c(ii,jj)……式中没包括1/(2Δt)时,随着ip=ms1+1,ms的单元循环而对有关单元求其和

e(i)=-tsss/30+e(i)

e(j)=-tsss/30+e(j)

e(k)=-tsss/30+e(k)

c越流矩阵元素,第ip单元1m水头差时的越流量(m3/d/m)均分给三节点i,j,k,随着ip=ms1+1,ms的单元循环而对有关单元求其和

yo(i)=yo(i)+rmkss/30

yo(j)=yo(j)+rmkss/30

yo(k)=yo(k)+rmkss/30

c 计算ip单元各节点的导水矩阵元素d(i,j),(i=1,2,3)(j=1,2,3)

do 100 iu=1,3

do 104 iuu=1,3

c 当iu=1时,即ip单元的i节点,分别计算iuu=1,2,3,即ip单元的i,j,k节点对i节点的ai,ai即公式的没求和部分

i=in(ip,iu)

j=in(ip,iuu)

ai=txy(bc(iu,1)bc(iuu,1)+bc(iu,2)bc(iuu,2))/ss/40

c 在ip单元的三个节点中,排除第3点,只让第1,2点(i,j点)的ai加入i点的导水矩阵元素d(i,j)中,第243行的j可分别轮到i,j,k三点,但第245句的1~9个中,仅有i,j点在i点的idyh中,此句排除了第3点加入i点的导水矩阵元素d(i,j)中随着ip=ms1+1,ms的循环,到251句时承压水导水矩阵已完全形成

do 106 k=1,9

if(jeqidyh(i,k))d(i,k)=ai+d(i,k)

106 continue

104 continue

100 continue

80 continue

c 至此,中、深层几何量计算完,以下开始时段循环;时段循环中把潜水也加进来

c 在时段循环中,潜水的导水矩阵,释水矩阵,越流矩阵与含水层厚度有关,所以在时段循环中完成,另外,承压水与时段有关的部分也在时段循环中完成

do 280 ikv=1,ihv

write(,)’ikv=’,ikv

c 该时段步长赋给dt

dt=tl(ikv)

c 以下开始计算潜水末刻流场

c 潜水已知水头节点的已知水头赋给本时段初、末刻

do 32 i=1,nn1

hed(n1+i)=h0nn1(i,ikv)

h0(n1+i)=h0nn1(i,ikv)

32 continue

c 先计算潜水几何量,因为潜水的含水层厚度M 随时段而变,所以放在时段循环内

do 33 i=1,ids1

h0(i)=hed(i)

c 本时段潜水各节点的含水层厚度ep0(i)用时段初刻水位=hed(i)减去底板标高zzz(i)求得

ep0(i)=hed(i)-zzz(i)

c本时段潜水各节点的灌溉回归水量加降水回归水量w(i,ikv)先赋零,以备求和之用

w(i,ikv)=00

33 continue

c 该时段潜水各单元的含水层厚度用三节点之均值

do 299 i=1,ms1

299 ep1(i)=(ep0(in(i,1))+ep0(in(i,2))+ep0(in(i,3)))/3

do 2050 i=1,ids1

c 同前,对潜水导水矩阵工作单元idyh,释水矩阵e,越流矩阵yo,导水矩阵d,先赋0

do 2050 j=1,9

idyh(i,j)=0

e(i)=00

yo(i)=00

2050 d(i,j)=00

do 8015 ip=1,ms1

i=in(ip,1)

j=in(ip,2)

k=in(ip,3)

i1=i

j1=j

k1=k

do 9015 iv=1,3

idyh(i1,1)=i1

do 9415 iu=j1,k1,k1-j1

do 9615 i2=2,9

if(idyh(i1,i2)eqiu)goto 9415

if(idyh(i1,i2)ne0)goto 9615

if(idyh(i1,i2)eq0)idyh(i1,i2)=iu

goto 9415

9615 continue

9415 continue

iuu=i1

i1=j1

j1=k1

k1=iuu

9015 continue

xi=zb(i,1)

xj=zb(j,1)

xk=zb(k,1)

yi=zb(i,2)

yj=zb(j,2)

yk=zb(k,2)

ss=abs((xj-xi)(yk-yi)-(xk-xi)(yj-yi))05

jvv=mqh(ip)

dyk=fqc(jvv,4)

dyk1=gh(jvv)

do 776 ig=1,3

ngh=in(ip,ig)

do 777 jb=1,jbn

if(ngheqjbiao(jb))goto 778

777 continue

w(ngh,ikv)=w(ngh,ikv)-qxia(ikv)ss/3

778 w(ngh,ikv)=w(ngh,ikv)+x(ikv)/dtdykss/3+gq(ikv)dyk1ss/3

776 continue

bc(1,1)=yj-yk

bc(1,2)=xj-xk

bc(2,1)=yk-yi

bc(2,2)=xk-xi

bc(3,1)=yi-yj

bc(3,2)=xi-xj

jv=mqh(ip)

txy=fqc(jv,1)ep1(ip)

ts=fqc(jv,2)

rmk=fqc(jv,3)

e(i)=-tsss/30+e(i)

e(j)=-tsss/30+e(j)

e(k)=-tsss/30+e(k)

yo(i)=yo(i)+rmkss/30

yo(j)=yo(j)+rmkss/30

yo(k)=yo(k)+rmkss/30

do 1001 iu=1,3

do 1041 iuu=1,3

i=in(ip,iu)

j=in(ip,iuu)

ai=txy(bc(iu,1)bc(iuu,1)+bc(iu,2)bc(iuu,2))/ss/40

do 1061 k=1,9

if(jeqidyh(i,k))d(i,k)=ai+d(i,k)

1061 continue

1041 continue

1001 continue

8015 continue

c 至此,该时段潜水几何量计算完

c 时段末水头迭代计数器ikv2,最大误差记录amax

ikv2=0

9881 amax=00

ikv2=ikv2+1

write(,34)ikv,ikv2

34 format(3x,i4,20x,’ikv2=’,i4)

c 开采井点录用计数器iq

iq=1

c 对潜水n1个未知水头节点逐点计算该时段末刻水头

do 4002 i=1,n1

c i节点常数项res:减该时段i点降水灌溉回归量,减i点该时段来自承压水的越流量(越流量计算采用时段初承压水水头,时段末潜水水位),减i点该时段储存量的增加量

res=-w(i,ikv)-(hed(i+ids1)-hed(i))yo(i)-e(i)(hed(i)-h0(i))/dt

if(ieqifdh(iq))then

c 如果i节点恰是开采节点,则i节点常数项res中再加上该时段i点的开采量,然后开采井点录用计数器iq加1,开采井点只能被1个节点录用

res=res+q(iq,ikv)

iq=iq+1

endif

c 把与i节点共单元的k节点的导水矩阵元素(k=1,2,…,最多可9,当k=1时即i点自己)依次加进i节点常数项res中,k节点的导水矩阵元素采用上一轮迭代出的hed(k)计算

do 3002 j=1,9

k=idyh(i,j)

if(keq0)goto 3002

res=res+d(i,j)hed(k)

3002 continue

c 把常数项除以:(i节点对i节点自身导水矩阵元素d(i,1)的负数,加,i节点储水阵元素e(i)/dt)

res=res/(-d(i,1)+e(i)/dt)

c 把i节点的res乘以亚松弛系数zz1,加给上一轮迭代出的hed(i)中,作为这一轮迭代出的hed(i):

hed(i)=reszz1+hed(i)

c 把1-n1个节点中最大的res挑出,并把其点号记到imax1中:

if(abs(res)leamax)goto 4002

amax=abs(res)

imax1=i

c 循环返回,继续下一个节点:

4002 continue

988 amax=00

232 continue

write(,)’amax1=’,amax,imax1

if(amaxgt9999999e-03)goto 9881

c 至此,潜水本时段末刻流场计算完,以下开始中层水本时段末刻流场的计算:

do 232 i=1,nn1

hed(n2+i)=h0nn2(i,ikv)

h0(n2+i)=h0nn2(i,ikv)

iqq=mi1+1

do 400 i=ids1+1,n2

c 其中的越流流入量采用(潜水本时段末水位hed(i-ids1)减承压水本时段末水位hed(i))之差hedd计算

hedd=hed(i-ids1)-hed(i)

heddd=hed(i)-hed(i+ids1)

res=-heddyo(i)+hedddyo(i+ids1)-e(i)(hed(i)-h0(i))/dt

50 if(ieqifdh(iqq))then

res=res+q(iqq,ikv)

iqq=iqq+1

endif

do 300 j=1,9

k=idyh(i,j)

if(keq0)goto 300

res=res+d(i,j)hed(k)

300 continue

res=res/(-d(i,1)+e(i)/dt)

hed(i)=reszz2+hed(i)

if(abs(res)leamax)goto 400

amax=abs(res)

imax2=i

400 continue

write(,)’amax2=’,amax,imax2

if(amaxgt9999999e-03)goto 988

c 至此,中层水全部节点末刻水头迭代完一轮,节点误差最大者如果 >001m,则988 句进行下一轮迭代,如果≤001m,则迭代结束,中层水本时段末刻流场计算完,以下开始深层水本时段末刻流场的计算:

989 amax=0

iqqq=mi2+1

do 332 i=1,nn1

hed(n+i)=h0nn3(i,ikv)

h0(n+i)=h0nn3(i,ikv)

332 continue

do 401 i=ids2+1,n

hedd=hed(i-ids1)-hed(i)

res=-heddyo(i)-e(i)(hed(i)-h0(i))/dt

5011 if(ieqifdh(iqqq))then

res=res+q(iqqq,ikv)

iqqq=iqqq+1

endif

do 301 j=1,9

k=idyh(i,j)

if(keq0)goto 301

res=res+d(i,j)hed(k)

301 continue

res=res/(-d(i,1)+e(i)/dt)

hed(i)=reszz2+hed(i)

if(abs(res)leamax)goto 401

amax=abs(res)

imax3=i

401 continue

write(,)’amax3=’,amax,imax3

if(amaxgt9999999e-03)goto 989

c至此,深层水全部节点本时段末刻水头迭代完一轮,节点误差最大者如果>001m,则返回989句进行下一轮迭代,如果≤001m,则迭代结束,深层水本时段末刻流场计算完

c以下开始把本时段末刻水位潜,中,深层水全部节点本时段末刻水头hed赋给下时段初刻水位h0,三层拟合点本时段末刻水位hed存入fh

189 do 281 i=1,ids

h0(i)=hed(i)

281 continue

do 190 j=1,igs

mp=igdh(j)

fh(j,ikv)=hed(mp)

190 continue

280 continue

c 至此,三层全包括的时段循环完

c 屏幕输出计算区总面积sum2/2,(m2);sum2是中、深层两层的累加和:

write(,)’sum=’,sum2/2,’m2’

do 1926 i=1,ids

1926 sss(i)=hed(i)-h(i)

c’ok1’输出潜水的最终流场:节点号i、坐标zb、最终流场hed、该节点0-ihv时间(最终流场比初始流场)的水位总变幅:

open(1,file=’ok1’)

write(1,1927)(i,zb(i,1),zb(i,2),hed(i),sss(i),i=1,ids1)

close(1)

c’ok2’输出中层水的最终流场:节点号i、坐标zb、最终流场hed、该节点0-ihv时间(最终流场比初始流场)的水位总变幅:

open(1,file=’ok2’)

write(1,1927)(i,zb(i,1),zb(i,2),hed(i),sss(i),i=ids1+1,ids2)

close(1)

c’ok3’输出深层水的最终流场:节点号i、坐标zb、最终流场hed、该节点0-ihv时间(最终流场比初始流场)的水位总变幅:

open(1,file=’ok3’)

write(1,1927)(i,zb(i,1),zb(i,2),hed(i),sss(i),i=ids2+1,ids)

close(1)

1927 format(i4,2f70,2f72)

c 计算出拟合点各时段误差(计算水位-实测水位)s(i,iuv)

do 8 iuv=1,ihv

do 2003 i=1,igs

s(i,iuv)=fh(i,iuv)-sh(i,iuv)

2003 continue

8 continue

c 输出拟合点各时段误差

open(1,file=’gandat’)

do 3511 i=1,igs,5

write(1,2993)file1(i),file1(i+1),file1(i+2),file1(i+3),file1(i+4)

2993 format(5(3x,a8,3x))

do 3778 iv=1,ihv

f0=fh(i,iv)

f1=fh(i+1,iv)

f2=fh(i+2,iv)

f3=fh(i+3,iv)

f4=fh(i+4,iv)

s0=s(i,iv)

s1=s(i+1,iv)

s2=s(i+2,iv)

s3=s(i+3,iv)

s4=s(i+4,iv)

write(1,3556)iv,f0,s0,f1,s1,f2,s2,f3,s3,f4,s4

3778 continue

3556 format(i3,10f72)

3511 continue

close(1)

c 输出拟合节点计算的历时水位fh(i,iuv)、实测的历时水位sh(i,iuv)、拟合误差s(i,iuv),拟合节点名称file1(i)do 2013 i=1,igs

tt=0

filee=file2(i)

open(1,file=filee)

do 2018 iuv=1,ihv-1

tt=tt+tl(iuv)

write(1,2020)tt,fh(i,iuv),sh(i,iuv),s(i,iuv)

2018 continue

write(1,2030)tt,fh(i,iuv),sh(i,iuv),s(i,iuv),file1(i)

close(1)

2013 continue

2020 format(1x,f103,3f153,2x)

2030 format(1x,f103,3f153,2x,’"’,a8,’"’)

stop

end

问下。。你这个程序是单纯单线程程序,还是能够进行并行处理的程序?一般VC++简单的程序是只能调用一颗逻辑cpu进行运算的。。。intel i5 的处理器物理双核、逻辑四核。只有进行并行编写的c++程序才能调用全部四个逻辑cpu,即在任务管理器Cpu占用率100%、单线程程序只能调用一颗逻辑cpu,即25%,

程序员不可以按提交代码量计算工作量。代码的数量不能作为程序员工作量的主要指标,有时候写出精简的代码,比写冗长的代码更难,用代码的多少来评估程序员的业务能力是否达标,是不客观的。衡量程序员的工作量是一件有趣的事情,不清楚做的事情以及程序的现状,是没有办法评估衡量的。

当然有可能

运行大型程序,最好在每一步加上输出信息,便于知道程序的运行状态

==============

最简单的办法,每运行一个循环,输出当前的步数

for ii = 1:N

/原有的程序代码/

disp(['calculating ',num2str(ii)])

end

我以前写过一篇关于计算矩阵的秩的小论文,里面是我的一些看法,我从中摘录了一部分,附在下面,看看对你有没有什么帮助。我的看法也是通过将矩阵化成最简形来求解,以下是这么选择的原因。其实这个问题可以讨论讨论的,当时我对自己的算法也不算很满意,所以有什么问题尽管提。

本程序是为求解矩阵的秩而进行编译的。要说明其功能,首先要明白什么是矩阵的秩。设在矩阵A中有一个不等于0的r结子式D,且所有r+1阶子式(如果存在的话)全等于0,那么D称为矩阵A的最高阶非零子式,数r称为矩阵A的秩,记作R(A)。零矩阵的秩为1。根据定义推断,计算矩阵A的秩,可以转化为计算矩阵A的最大非零子式。但是,实际应用这条原理来解决此问题并不容易。因为,应用计算机计算矩阵A所对应的行列式|A|的值非常麻烦。一个m×n的矩阵,其k阶子式多达m!/[k!•(m-k)!]•n!/[k!•(n-k)!]个,这大大增加了程序的计算量。同时,由于不同阶的子式的值的算法不易通用,故也增加了程序员的编程负担,最重要的是,程序的通用性较低,不易应用于相似题目的求解。故,本程序算法并未采用这种思路。那么,本题又应当如何求解呢?实际生活中,我们一般的求解方法是应用初等变换求解。应用初等变换,将要求的矩阵A变换成行最简形或列最简形然后再进行判断,这才是我们求解矩阵的秩的常规做法。那么,编写程序求解矩阵的秩当然也可以遵循这种做法。相对于前面所讨论的原理来说,应用这种原理进行算法设计,可以减少不少的时间,同时计算机求解的速度也能大大提高。而且,再本算法的基础上稍加改进,即可适应任何阶次的矩阵的秩的求解。

以上就是关于计算程序全部的内容,包括:计算程序、为什么编写计算量大的程序时候,任务管理器上显示程序的CPU占用率保持在25%、程序员可以按提交代码量计算工作量吗等相关内容解答,如果想了解更多相关内容,可以关注我们,你们的支持是我们更新的动力!

欢迎分享,转载请注明来源:内存溢出

原文地址: http://outofmemory.cn/zz/9799768.html

(0)
打赏 微信扫一扫 微信扫一扫 支付宝扫一扫 支付宝扫一扫
上一篇 2023-05-02
下一篇 2023-05-02

发表评论

登录后才能评论

评论列表(0条)

保存