>我首先编写了一个包含 6 个子例程的代码，我在主程序中编写了输出，我得到了正确的文本文件，但现在我还想在同一文本文件中添加来自另一个子例程的一个变量的输出这是我的代码示例
! code for the particle tracing in magnetic field
implicit none
include "mydata.inc"
External ad,bsstep
real*8 x(ndim,nv),b(4)
real*8 pwE(ndim,3),pwB(ndim,3),ep
real*8 ke,gamma,kew,epx,hmin,rr,epx0
real*8 eps(3),eta0(5)
real*8 m,d_eta
real*8 t,tw,dt,rt
integer i,j,k,nok,nbad,filen
character*80 fname1
! real,parameter :: pi
t=0.d0
tw=0.d0
m=me
x=0.d0
pwE=0.d0
pwB=0.d0
!
!=====loop================!
k=1
120 call ini_v(x,ke,pwB)
x(2,3)=eta0(k)
call dispersion(pwE,pwB,eps,ep)
write(fname1,1005) k
open(1,file=fname1)
write(1,*)'time,eta,phi,pitch angle,ke,deta'
!====normalization========!
ke=ke/(m*c**2)
gamma=ke+1
x(1,2)=cos(alpha)*dsqrt(ke**2+2*ke)
x(2,2)=sin(alpha)*dsqrt(ke**2+2*ke)
dt=0.001
epx=1.d-6 ! required accuracy
epx0=epx
hmin=0.d0
!================================
rt=0
kew=ke*(m*c**2)/1.0d6/e !ke in the units of electron volts
100 write(*,*)rt
if(rt.gt.15)then
goto 110
endif
if(rt.eq.0)then
write(1,1000)t/(2*pi),
&(x(i,3),i=1,2),alpha*180.d0/pi,kew,d_eta
else
write(1,1000)t/(2*pi),
&(x(i,3),i=1,2),alpha*180.d0/pi,kew,d_eta
endif
!=========adaptivestepsize method==============!
call odeint(x,t,t+dt,epx,dt,hmin,nok,nbad,ad,bsstep,
&pwE,pwB,gamma,ep)
x(1:2,3)=mod(x(1:2,3),2*pi) !this is the range of angles
alpha=dacos(x(1,2)/dsqrt(x(1,2)**2+x(2,2)**2))
call go(KE,x(1:2,2),gamma)
KEw=KE*(m*c**2)/1.0d6/e
t=t+dt
rt=rt+dt
goto 100
110 write(*,*) eta0(k)
if(k.le.size(eta0)) then
k=k+1
t=0
goto 120
endif
1000 format(6(E10.4,2x),4(f15.6,2x))
1005 format("output",i2.2,".txt")
stop
end
我想从哪里获取值的子程序是这个
====
================================================! subroutine to calculate all the derivatives for this code
subroutine ad(step,kin,kout,pwE,pwB,gamma,ep)
include "mydata.inc"
real*8 step
real*8 kin(ndim,nv),kout(ndim,nv)
real*8 pwE(ndim,3),pwB(ndim,3),ep
real*8 gamma
real*8 tm(3)
kout(1,1)=kin(2,2)*cos(kin(2,3)-kin(1,3))/gamma
kout(2,1)=kin(2,2)*sin(kin(2,3)-kin(1,3))/gamma
kout(3,1)=kin(1,2)/gamma
kout(1,2)=-pwE(3,1)*sin(kin(1,3))
&+(pwB(1,1)*kin(2,2)/(2.*gamma))
&*((1+ep)*sin(kin(2,3))+(1-ep)*sin(kin(2,3)-2*kin(1,3)))
tm(1)=-pwB(1,1)*kin(1,2)/2./gamma
tm(2)=pwB(1,1)*emicwave/wavek/cos(si)
tm(3)=tan(si)*pwE(3,1)
kout(2,2)=sin(kin(2,3))*(tm(1)*(1+ep)+1./2.*((1+ep)*tm(2)-tm(3)))
&+sin(kin(2,3)-2*kin(1,3))
&*(tm(1)*(1-ep)+1./2.*((1-ep)*tm(2)+tm(3)))
kout(1,3)=EMICwave-kin(1,2)*wavek*cos(si)/gamma
&-kin(2,2)*wavek*sin(si)*cos(kin(2,3)-kin(1,3))/gamma
!----------------------d(eta)/dt --------------------------!
! kin(2,3)=eta, kout(2,3)=d(eta)/dt
! time derivatives of phase(phi)
!----------------------------------------------------------!
tm(1)=-2*tm(1)
kout(2,3)=EMICwave-kin(1,2)*wavek*cos(si)/gamma
&-kin(2,2)*wavek*sin(si)*cos(kin(2,3)-kin(1,3))/gamma
&+1./gamma-cos(kin(2,3))/(2.*kin(2,2))
&*((1+ep)*tm(1)-(1+ep)*tm(2)+tm(3))
&-cos(kin(2,3)-2*kin(1,3))/(2.*kin(2,2))
&*((1-ep)*tm(1)-(1-ep)*tm(2)-tm(3))
&-pwB(3,1)/gamma*cos(kin(1,3))
return
end
首先：我会远离小于 10 的文件单元 - 你永远不会知道你的编译器是否恰好将那个特定的数字用于标准输出或类似的东西(大多数使用单元 0 表示标准错误，5 和 6 表示标准输入和标准输出，但不能保证(。
最安全的方法是在打开期间使用newunit参数：
integer :: unt
open(newunit=unt, file=filename)
write(unt, *) "Some text"
close(unt)
如果在调用子例程时文件处于打开状态，则子例程理论上应该只需要该文件单元即可将输出写入该文件(假设子例程未pure(。
如果需要，可以将文件 I/O 传递给模块，然后在主程序和子例程中use该模块。