我在这个 fortran90 代码中面临一个难以调试的问题——保留或注释掉
write(*,*)操作系统：Ubuntu 16.04
gfortran --version
GNU Fortran (Ubuntu 5.4.0-6ubuntu1~16.04.12) 5.4.0 20160609
版权所有 (C) 2015 Free Software Foundation, Inc.
FORTRAN90 代号：gillsp_notworking.f90
编译：gfortran -ffree-line-length-none gillsp_notworking.f90
这是两个有问题的行（标记为！***** error write ****）：
W1：在第 416 行写语句
W2：在第 417 行写语句
代码：
program code
implicit none
real*8, parameter :: time=3600.d0, t0=0.d0
real*8, parameter :: er=1.d-6
real*8, parameter :: vton=600.d0
real*8, parameter :: dl=1.d0/370.d0
real*8,parameter :: kp1be = 11.6d0/vton, km1be = 1.4d0
real*8,parameter :: kp2be = 3.4d0/vton, km2be = 0.2d0
real*8,parameter :: kp3be = 2.9d0/vton, km3be = 5.4d0
real*8,parameter :: kp1pe = 0.d0/vton, km1pe = 0.d0
real*8,parameter :: kp2pe = 0.d0/vton, km2pe = 0.d0
real*8,parameter :: kp3pe = 0.d0/vton, km3pe = 0.d0
real*8,parameter :: kpdom = 16.d-5/vton
real*8,parameter :: kpcf = 13.d0/vton, kmcf = 0.7d0
real*8,parameter :: kmbcf = 4.4, kmpcf = 0.d0
real*8,parameter :: ktd = 0.003d0
real*8,parameter :: ktdpi = 0.3d0
real*8,parameter :: ksever = 0.01d0
real*8,parameter :: kcp0=12.8d0/vton, kcm0=2.d-4
real*8,parameter :: kcp1=0.3d0/vton, kcm1=6.d-3
real*8,parameter :: alpha = 20.d0
real*8,parameter :: rho = 0.3d0
real*8,parameter :: rhocf = 0.8d0
real*8,parameter :: rhocp = 0.d0
real*8,parameter :: V = 1.d0
integer, parameter :: Nmin=100,Npos=nint(60.d0/dl)
integer, parameter :: N=nint(rho*vton*V)
integer, parameter :: Ncf=nint(rhocf*vton*V)
integer, parameter :: Ncp=nint(rhocp*vton*V/1000.d0)
integer, parameter :: nstep=300
integer, parameter :: ndmax=3000
integer, parameter :: minsevL=10
integer, parameter :: plotmovie=1
real*8, parameter :: delt=time/nstep
real*8 :: t,r1,r2,pl,ph
real*8 :: tau
real*8 :: beta1,beta2,beta3,beta4,beta5,beta6,beta7,beta8,beta_sev
real*8 :: beta_cp,beta_cm
real*8 :: beta_sum
real*8 :: alpha0
real*8 :: rsev1,rsev2
real*8 :: pifactr
integer :: m, m1(-Nmin:Npos), m10(-Nmin:Npos)
integer :: mcf, domL(1:ndmax),domR(1:ndmax), ndom,cofcount,domLen
integer :: cap
integer :: domL0(1:ndmax),domR0(1:ndmax),drcount,domrem(1:ndmax)
integer :: pesever(1:10),besever(1:10),besvcount,pesvcount, sevindex, dumind
integer :: i,j,k,nk
integer :: l1(-Nmin:Npos)
integer :: l1bind, l1pind, l1bind0, l1pind0
integer :: gatp1, gadppi1, gadp1, glcf1, gzero
integer :: l1sum
integer :: k1,k2, dgcount, dscount, nmerge, domflag
integer :: nfile
integer :: dumdl
integer :: Bstate
integer :: nPicof
character*80::fname, Rid
write(*,*) 'Array size=', -Nmin, Npos
!------------------------------------------------------
t=0
nk=0
nfile=0
!----------- initialisation -----------
l1=0
m1=0 !! monomer states: 1 ATP, 2 ADP.Pi, 3 ADP, 0 empty
l1bind=0 !! index of incoming monomer in barbed end
l1pind=-1 !! index of incoming monomer in pointed end
m=N !! actin monomers
cap=Ncp !! free capping
Bstate = 1 !! Free barbed end
pifactr = 1.d0 !! factor increase in Pi release rate
ndom=0
cofcount=0
domL=0
domR=0
nPicof=0 !! number of ADP.Pi and cofilin interface
mcf=Ncf-sum(domR-domL)-ndom
write(200,*)t,sum(l1),m, l1pind, l1bind, cap, mcf, cofcount
!********************************************************************
!********************************************************************
!------------------- time loop ------------------------
tloop: do while (t.le.time)
!----------------- counting G-atp, G-adp, G-cofilin in filament --------------
gzero=count(m1<er)
gatp1=count(m1<1+er) - gzero
gadppi1=count(m1<2+er) - (gatp1+gzero)
gadp1=count(m1<3+er) - (gatp1+gzero+gadppi1)
glcf1=count(m1>3+er)
!----------------- current length ---------------------
l1sum = sum(l1)
!------------------ propensity sum ---------------------
!---- store monomer states, domain info etc -----
m10 = m1
l1bind0=l1bind
l1pind0=l1pind
domL0=domL
domR0=domR
!------------------- capping reaction ------------------
if (m10(l1bind-1).eq.4) then
if(Bstate.eq.1)then
beta_cp = kcp1*cap/V
beta_cm = 0.d0
elseif (Bstate.eq.2) then
beta_cp = 0.d0
beta_cm = kcm1
endif
else
if(Bstate.eq.1)then
beta_cp = kcp0*cap/V
beta_cm = 0.d0
elseif (Bstate.eq.2) then
beta_cp = 0.d0
beta_cm = kcm0
endif
endif
!-------------------- barbed end -----------------------
if(Bstate.eq.1)then ! Barbed-end state if
if(m10(l1bind-1).eq.1)then
beta1 = km1be
beta2 = kp1be*(m/V)
elseif (m10(l1bind-1).eq.2) then
beta1 = km2be
beta2 = kp2be*(m/V)
elseif (m10(l1bind-1).eq.3) then
beta1 = km3be
beta2 = kp3be*(m/V)
elseif (m10(l1bind-1).eq.4) then
beta1 = kmbcf
beta2 = 0.d0
endif
elseif (Bstate.eq.2) then
beta1 = 0.d0
beta2 = 0.d0
endif
!----- the first monomer incorporation------!
if(t.lt.30)then
if(m10(l1bind-1).eq.0)then
beta2 = kp1be*(m/V)
endif
endif
!-------------------- pointed end -----------------------
if(m10(l1pind+1).eq.1) then
beta3 = km1pe
beta4=kp1pe*(m/V)
elseif (m10(l1pind+1).eq.2) then
beta3 = km2pe
beta4=kp2pe*(m/V)
elseif (m10(l1pind+1).eq.3) then
beta3 = km3pe
beta4 = kp3pe*(m/V)
elseif (m10(l1pind+1).eq.4) then
beta3 = kmpcf
beta4 = 0.d0
endif
!------------------- bulk reactions ---------------------
!------------------ cofilin reactions ------------------
if(ndom.gt.0)then
nmerge=0
do k1=1,ndom
do k2=1,ndom
if(domL(k1)-1.eq.domR(k2)) nmerge=nmerge+1 ! single free site between two cof domains
enddo
enddo
endif
!---------------- Pi-cofilin interface ------------------
nPicof=0
do k1=l1pind0,l1bind0-1
if(abs(m10(k1)-m10(k1+1)).eq.2.and.m10(k1)+m10(k1+1).eq.6)then
nPicof=nPicof+1
endif
enddo
!------------------- bulk reactions ---------------------
beta5 = ktdpi*gatp1 + ktd*(gadppi1-nPicof) + alpha*ktd*nPicof !---- total rate of hydrolysis
beta6 = kpdom*(mcf/V)*max(0.d0,(gadp1-2.d0*ndom-2.d0*nmerge)) !---- domain creation rate
!----- count domain growth probability -----
if(ndom.gt.0)then
dgcount=0
dscount=0
do k=1,ndom
dumdl=domR(k)-domL(k)-1 ! current domain length
if(m10(domL(k)).eq.3) dgcount=dgcount+1
if(m10(domR(k)).eq.3) dgcount=dgcount+1
if(m10(domL(k)).le.3.and.dumdl.gt.2) dscount=dscount+1
if(m10(domR(k)).le.3.and.dumdl.gt.2) dscount=dscount+1
enddo
else
dgcount=0
dscount=0
endif
beta7 = kpcf*(mcf/V)*dgcount !---- domain growth rate
beta8 = kmcf*dscount !---- domain shrinkage rate
!----- severing probability -----
if(ndom.gt.0)then
besvcount=0
pesvcount=0
besever=-1
pesever=-1
do k=1,ndom
domLen = domR(k)-domL(k)-1
if (domLen.ge.minsevL) then !---- domain length condition
if(m10(domL(k)).lt.4.and.m10(domL(k)).gt.0) then
pesvcount=pesvcount+1
pesever(pesvcount)=domL(k)
endif
if(m10(domR(k)).lt.4.and.m10(domR(k)).gt.0) then
besvcount=besvcount+1
besever(besvcount)=domR(k)
endif
else !------ domain length condition
pesvcount=0
besvcount=0
endif !------ domain length condition
enddo
else
pesvcount=0
besvcount=0
endif
beta_sev = ksever*(besvcount+pesvcount) !---- severing rate
beta_sum = beta1+beta2+beta3+beta4+beta5+beta6+beta7+beta8+beta_sev+beta_cp+beta_cm
alpha0 = beta_sum
!----------------- next reaction time: tau --------------
r1 = rand()
do while(r1.lt.1.d-8)
r1 = rand()
enddo
tau = (1.d0/alpha0)*log(1.d0/r1)
!----------------------- all reactions ------------------------
r2 = rand()
!---------------- growth and decay ------------------
!---------------- filament 1 ----------------
!---------------- barbed end ----------------
pl=0.d0
ph=beta1/alpha0
if (r2.ge.pl.and.r2.lt.ph.and.l1sum.gt.er) then ! barbed end disassembly
if (m1(l1bind-1).eq.4) then
! mcf = mcf+1 Pool conc. remains same
do k=1,ndom
if(domR(k).eq.l1bind) domR(k)=domR(k)-1
enddo
endif
l1bind = l1bind - 1
l1(l1bind) = 0
m1(l1bind) = 0
! m = m+1 Pool conc. remains same
Rid="BE-"
endif
pl=beta1/alpha0
ph=(beta1+beta2)/alpha0
if (r2.ge.pl.and.r2.lt.ph) then ! barbed end assembly
m1(l1bind) = 1
l1(l1bind) = 1
l1bind = l1bind + 1
! m = m-1 Pool conc. remains same
Rid="BE+"
endif
!---------------- pointed end ---------------
pl=(beta1+beta2)/alpha0
ph=(beta1+beta2+beta3)/alpha0
if (r2.ge.pl.and.r2.lt.ph.and.l1sum.gt.er) then ! pointed end disassembly
if (m1(l1pind+1).eq.4) then
! mcf = mcf+1 Pool conc. remains same
do k=1,ndom
if(domL(k).eq.l1pind) domL(k)=domL(k)+1
enddo
endif
l1pind = l1pind + 1
l1(l1pind) = 0
m1(l1pind) = 0
! m = m+1 Pool conc. remains same
Rid="PE-"
endif
pl=(beta1+beta2+beta3)/alpha0
ph=(beta1+beta2+beta3+beta4)/alpha0
if (r2.ge.pl.and.r2.lt.ph) then ! pointed end assembly
m1(l1pind) = 1
l1(l1pind) = 1
l1pind = l1pind - 1
! m = m-1 Pool conc. remains same
Rid="PE+"
endif
!---------- severing of filament -------------
pl=(beta1+beta2+beta3+beta4)/alpha0
ph=(beta1+beta2+beta3+beta4+beta_sev)/alpha0
if (r2.ge.pl.and.r2.lt.ph) then ! filament severing if
rsev1=rand()
rsev2=rand()
if(rsev1.le.0.2d0.and.besvcount.gt.0)then ! find severing dir if loop
dumind = 1 + nint(rsev2*(besvcount-1)) ! use besever
sevindex = besever(dumind) ! assumption: barbed end is not cofilin-bound (not good with capping)
else
if (pesvcount.gt.0) then
dumind = 1 + nint(rsev2*(pesvcount-1)) ! use pesever, if pesevcount=0, use besever
sevindex = pesever(dumind)
else
dumind = 1 + nint(rsev2*(besvcount-1)) ! use besever
sevindex = besever(dumind)
endif
endif ! find severing dir if loop
l1bind=sevindex ! keep olp pointed end
Bstate=1 ! set new barbed end free
! m=m+(l1bind0-sevindex) Pool conc. remains same
m1(sevindex:l1bind0)=0
l1(sevindex:l1bind0)=0
drcount=0
domrem=0
do j=sevindex,l1bind0 ! j-loop
! if(m10(j).eq.4) then Pool conc. remains same
! mcf=mcf+1
! endif
do k=1,ndom
if(domL(k).eq.j)then
drcount=drcount+1
domrem(drcount)=j
domL(k)=0
domR(k)=0
endif
enddo
enddo ! j-loop
Rid="SEVER"
write(*,*)'sever', m10(sevindex-1) ,sevindex, l1bind0, besvcount, pesvcount !***** error write ****
write(*,*)'sever', m10(sevindex-1) !***** error write ****
endif ! filament severing if
!---------------- capping ------------------
pl=ph
ph=(beta1+beta2+beta3+beta4+beta_sev+beta_cp)/alpha0
if (r2.ge.pl.and.r2.lt.ph) then
! cap=cap-1 Pool conc. remains same
Bstate = 2
Rid="CAPPING +"
endif
!---------------- de-capping ------------------
pl=ph
ph=(beta1+beta2+beta3+beta4+beta_sev+beta_cp+beta_cm)/alpha0
if (r2.ge.pl.and.r2.lt.ph) then
! cap=cap+1 Pool conc. remains same
Bstate = 1
Rid="CAPPING -"
endif
!---------- hydrolysis in filament -----------
pl=(beta1+beta2+beta3+beta4+beta_sev+beta_cp+beta_cm)/alpha0
do j=l1pind0,l1bind0 ! j monomer loop
if(m10(j).eq.1)then ! Hydrolysis: G-ATP -> G-ADP.Pi
ph=pl+(ktdpi/alpha0)
if (r2.ge.pl.and.r2.lt.ph) then
m1(j) = 2
Rid="HYD atp"
endif
pl=ph
elseif(m10(j).eq.2)then ! Hydrolysis: G-ADP.Pi -> G-ADP
if(m10(j-1).eq.4.or.m10(j+1).eq.4)then ! enhanced Pi release at domain boundaries
pifactr=alpha
else
pifactr=1.d0
endif
ph=pl+(pifactr*ktd/alpha0)
if (r2.ge.pl.and.r2.lt.ph) then
m1(j) = 3
Rid="Pi Release"
endif
pl=ph
elseif(m10(j).eq.3)then ! single Cofilin binding G-ADP -> G-ADP-Cofilin
domflag=0
if(ndom.gt.0)then ! ndom-if
do k=1,ndom
if(j.eq.domL(k)) then
domflag=1
endif
if(j.eq.domR(k)) then
domflag=1
endif
enddo
endif ! ndom-if
if (domflag.lt.1) then
ph=pl+((kpdom*mcf/V)/alpha0) ! Domian creation
if (r2.ge.pl.and.r2.lt.ph) then
m1(j) = 4
! mcf=mcf-1 Pool conc. remains same
ndom=ndom+1
domL(ndom)=j-1
domR(ndom)=j+1
Rid="DOM"
write(*,*) 'domain created', ndom
endif
pl=ph
endif
endif ! hydrolysis if
enddo ! j monomer loop
! DOMAINS ONLY GROW IF ADP ACTIN IS PRESENT
if(ndom.gt.0)then ! ndom if ! Cofilin domain growth
do k=1,ndom
dumdl=domR(k)-domL(k)-1 ! current domain length
!--------------------- domain growth --------------------
if (m10(domL(k)).eq.3) then ! Left boundary ADP-if
ph=pl+((kpcf*mcf/V)/alpha0)
if (r2.ge.pl.and.r2.lt.ph) then
m1(domL(k))=4
! mcf=mcf-1 Pool conc. remains same
domL(k)=domL(k)-1
Rid="DOML+"
endif
pl=ph
endif
if (m10(domR(k)).eq.3) then ! Right boundary ADP-if
ph=pl+((kpcf*mcf/V)/alpha0)
if (r2.ge.pl.and.r2.lt.ph) then
m1(domR(k))=4
! mcf=mcf-1 Pool conc. remains same
domR(k)=domR(k)+1
Rid="DOMR+"
endif
pl=ph
endif
! --------------------- domain shrinkage -------------------
if (m10(domL(k)).le.3.and.dumdl.gt.2) then ! Left boundary shrink-if
ph=pl+(kmcf/alpha0)
if (r2.ge.pl.and.r2.lt.ph) then
m1(domL(k)+1)=3
! mcf=mcf+1 Pool conc. remains same
domL(k)=domL(k)+1
Rid="DOML-"
endif
pl=ph
endif
if (m10(domR(k)).le.3.and.dumdl.gt.2) then ! Right boundary shrink-if
ph=pl+(kmcf/alpha0)
if (r2.ge.pl.and.r2.lt.ph) then
m1(domR(k)-1)=3
! mcf=mcf+1 Pool conc. remains same
domR(k)=domR(k)-1
Rid="DOMR-"
endif
pl=ph
endif
enddo ! k ndom loop
endif ! ndom if
!----------------- count bound cofilin --------------------
cofcount=0
do j=l1pind,l1bind
if (m1(j).eq.4) then
cofcount=cofcount+1
endif
enddo
!------------------- ERROR CHECKS ----------------------
if (tau.lt.0.d0) exit tloop
if (m1(l1bind).ne.0.or.m1(l1pind).ne.0) exit tloop
if (m.ne.N) exit tloop
if (mcf.ne.Ncf) exit tloop
if (cap.ne.Ncp) exit tloop
! if (sum(l1).eq.0) then
! write(*,*) 'FULLY DEPOLYMERIZED'
! t = time
! endif
!----------------- FULL DECORATION ----------------
! if (cofcount.eq.sum(l1)) write(*,*) 'fully decorated filament'
!-------------------------------------------------------
t = t + tau
enddo tloop !***** time loop *****
!------------------ time loop ended -----------------------
write(*,*) 'domains created', ndom
do k=1,ndom
write(*,*) k, 'size', domR(k)-domL(k)-1, 'L-R', domL(k), domR(k)
enddo
do k=l1pind,l1bind
write(101,*) k, m1(k)
enddo
if(t.ge.time) write(*,*) 'CODE ENDED WITHOUT ERROR'
write(*,*) 'CODE STOPPED ','t=',t,'id=',Rid
if (tau.lt.0.d0) write(*,*) 'ERROR: tau negative', tau, alpha0, beta1+beta2+beta3+beta4+beta5, beta6, beta7
if (m1(l1bind).ne.0.or.m1(l1pind).ne.0) write(*,*) 'ERROR: l1bind/l1pind not zero', m1(l1bind), m1(l1pind)
if (m.ne.N) write(*,*) 'ERROR: constant conc. fails', m, N
if (mcf.ne.Ncf) write(*,*) 'ERROR: constant conc. fails', mcf, Ncf
if (cap.ne.Ncp) write(*,*) 'ERROR: constant conc. fails', cap, Ncp
write(*,*) 'pe',l1pind,'be',l1bind, 'domains', ndom
stop
end program code
W1和W2被注释掉时的O/P：
Array size= -100 22200
domain created 1
domain created 2
domain created 3
domain created 4
domain created 5
domain created 6
domain created 7
domains created 7
1 size -1 L-R 0 0
2 size -17 L-R 50 34
3 size -1 L-R 0 0
4 size -1 L-R 0 0
5 size -1 L-R 0 0
6 size -1 L-R 0 0
7 size -1 L-R 0 0
CODE ENDED WITHOUT ERROR
CODE STOPPED t= 3600.2563350269443 id=BE-
pe -1 be 0 domains 7
Note: The following floating-point exceptions are signalling: IEEE_DENORMAL
W1 和 W2 未注释掉时的 O/P：
Array size= -100 22200
domain created 1
domain created 2
domain created 3
domain created 4
domain created 5
domain created 6
sever 4 138 620 2 2
sever 4
domain created 7
sever 2 106 390 1 1
sever 2
domain created 8
domain created 9
domain created 10
domain created 11
sever 3 237 584 1 1
sever 3
domain created 12
domain created 13
sever 2 227 362 2 2
sever 2
sever 3 198 303 1 1
sever 3
domain created 14
sever 2 346 622 1 1
sever 2
domain created 15
domain created 16
sever 3 320 783 2 1
sever 3
domain created 17
domain created 18
domain created 19
sever 3 405 862 1 1
sever 3
domain created 20
sever 4 607 893 1 1
sever 4
domain created 21
sever 3 512 867 2 2
sever 3
domain created 22
domain created 23
domain created 24
sever 3 576 1023 3 2
sever 3
domain created 25
domain created 26
domain created 27
sever 4 584 887 2 2
sever 4
domain created 28
domain created 29
domain created 30
domain created 31
domain created 32
sever 3 1003 1124 4 4
sever 3
domain created 33
sever 2 596 1128 2 2
sever 2
domains created 33
1 size -10 L-R 34 25
2 size -23 L-R 47 25
3 size -45 L-R 79 35
4 size 24 L-R 0 25
5 size -1 L-R 0 0
6 size 0 L-R 0 1
7 size -1 L-R 0 0
8 size -1 L-R 0 0
9 size 44 L-R 95 140
10 size 132 L-R 139 272
11 size -1 L-R 0 0
12 size -1 L-R 0 0
13 size -1 L-R 0 0
14 size -1 L-R 0 0
15 size 197 L-R 271 469
16 size -1 L-R 0 0
17 size -1 L-R 0 0
18 size -1 L-R 0 0
19 size -1 L-R 0 0
20 size -1 L-R 0 0
21 size -1 L-R 0 0
22 size 87 L-R 468 556
23 size -1 L-R 0 0
24 size -1 L-R 0 0
25 size 0 L-R 569 570
26 size -1 L-R 0 0
27 size 5 L-R 555 561
28 size -1 L-R 0 0
29 size -1 L-R 0 0
30 size -1 L-R 0 0
31 size -1 L-R 0 0
32 size -1 L-R 0 0
33 size -1 L-R 0 0
CODE ENDED WITHOUT ERROR
CODE STOPPED t= 3600.0229776615342 id=Pi Release
pe -1 be 860 domains 33
Note: The following floating-point exceptions are signalling: IEEE_DENORMAL
没有一个结果是非物理的，并且可以在 mac 上使用不同的编译器重现不同的结果。我也用
gfortran -Wall -Wno-tabs -g -fcheck=all -fbacktrace gillsp_notworking.f90