Fortran 代码 (.f95) 在 Windows g95 编译器中编译良好，但在 Ubuntu gfortran 中编译错误
Fortran code (.f95) compiles fine in Windows g95 compiler but incorrectly in Ubuntu gfortran
我正在尝试编译一个 .f95 fortran 脚本，以便它可以在 Ubuntu 上 运行。该脚本可在此处获得 -> Link to zip file containing .f95 script
当我切换到 Windows 并使用 g95 编译器编译时，它编译并且 运行 没问题。通过 wine 在 Ubuntu 中生成的 .exe 文件也 运行 没问题。
但是，如果我尝试编译生成 Ubuntu 文件，它无法正常工作。我没有收到编译错误，但如果我 运行 生成的文件，要么程序陷入无限循环，要么输出全错。我很难看出哪里出了问题，因为我没有写原始代码，对Fortran的理解也很不稳定，但这似乎与数字计算错误有关，导致非常large/small/inappropriately负输出（抱歉这么含糊）。
我是 运行ning 16.04 xenial ubuntu 和 gfortran 5.4.0。
任何 help/thoughts 感谢这让我感到困惑！谢谢
下面的代码供快速参考：
! Seed dispersal model of of Duman et al. (2015)
! Instructions and expample are found in:
! https://nicholas.duke.edu/people/faculty/katul/research.html
! in 'Library of Functions and Utilities'
! Author: Tomer Duman
! Version: Version 2
! Date: October 22, 2015
! References: Duman, T., Trakhtenbrot, A., Poggi, D.,
! Cassiani, M., Katul, G., Dissipation
! Intermittency Increases Long-Distance
! Dispersal of Heavy Particles in the Canopy
! Sublayer,
! accepcted to Boundary-Layer Meteorology
program LSmodel
implicit none
real :: sec,ran,gasdev ! random generator variables
real :: x,y,z,u,v,w,ut,vt,wt,t,dt ! simulation variables
real :: wg ! seed parametes
real :: Um,sigma_u,sigma_v,sigma_w,uw ! wind statistics variables
real :: dvaru_dz,dvarv_dz,dvarw_dz,duw_dz ! wind statistics variables
real :: dissip_m,TL ! vector over the range of ustars
real :: zs,zg,zmax ! release height & boundaries
real :: Ainv,C0inv ! inverse parameters
real :: C0,A,b,au,av,aw,dt_on_TL ! LS model parameters
real :: dz_max,dt_max ! time step limit
real :: CT,beta ! Crossing Trajectories correction
real :: C_chi,chi,TKE,T_chi,omega ! DI parameters
real :: a_ln,b_ln,sigma_chi,dissip_s ! DI parameters
real :: rhop,rho,r,g,gt,Re,AIP,Cd,nu ! IP parameters
real :: up,vp,wp,upt,vpt,wpt,vr,dt_ip,alpha ! IP parameters
integer :: seed ! random generator variables
integer :: pnum ! simulation parameters
integer :: i,j,jj,n,ii ! counting parameters
integer :: n_ip,IP=1 ! IP parameters
character(len=80) :: filename
real, allocatable,dimension(:) :: z_vec,Um_vec,sigma_u_vec,sigma_v_vec,sigma_w_vec,uw_vec
real, allocatable,dimension(:) :: dvaru_dz_vec,dvarv_dz_vec,dvarw_dz_vec,duw_dz_vec,dissip_m_vec
! setting the random generator seed
seed=7654321
sec=0.0
seed=seed+2*int(secnds(sec))
! input
open (21,file='input_parameters.txt')
read (21, *), x,C0,wg,zs,zg,beta,dt_on_TL,y,sigma_chi,C_chi,r,rhop,alpha,rho,nu
close(21)
pnum = int(x) ! number of released seeds
n = int(y) ! size of the input flow stats
wg = -1.0*wg ! seed terminal velocity [m/s]
!zs ! seed release height [m]
!C0 ! universal constant
!beta ! crossing trajectories parameter
!zg ! ground height [m]
!sigma_chi ! the standard deviation of chi (dissipation intermittency)
!C_chi ! constant for T_chi calc
!r ! particle radium [m] - set to 0 for no IP
!rhop ! particle dry density [kg/m^3]
!alpha ! drag parameter (Cd = 24/Re_p*(1+alpha*Re_p))
!rho ! fluid density [kg/m^3]
!nu ! fluid viscosity [m^2/s]
C0inv = 1.0/C0
g = 9.81
if (r==0.0) then
IP = 0
end if
! limiting parameters to prevent too big jumps in a time-step
dz_max = 0.1
dt_max = 0.1
open(unit=12,file='res.dat', form='formatted')
open(unit=13,file='res_traj.dat', form='formatted')
! allocate
allocate(z_vec(n))
allocate(Um_vec(n))
allocate(sigma_u_vec(n))
allocate(sigma_v_vec(n))
allocate(sigma_w_vec(n))
allocate(uw_vec(n))
allocate(dvaru_dz_vec(n))
allocate(dvarv_dz_vec(n))
allocate(dvarw_dz_vec(n))
allocate(duw_dz_vec(n))
allocate(dissip_m_vec(n))
! load normalized stats
open (22,file='input_flow.txt',form='formatted')
read (22,*) z_vec
read (22,*) Um_vec
read (22,*) sigma_u_vec
read (22,*) sigma_v_vec
read (22,*) sigma_w_vec
read (22,*) uw_vec
read (22,*) dvaru_dz_vec
read (22,*) dvarv_dz_vec
read (22,*) dvarw_dz_vec
read (22,*) duw_dz_vec
read (22,*) dissip_m_vec
close(22)
zmax = z_vec(n)
do i=1,pnum
t=0.0 ! initiate time and location
x=0.0
y=0.0
z=zs
do j=2,n ! interpolate
if ((z>=z_vec(j-1)).and.(z<=z_vec(j))) then
sigma_u=((sigma_u_vec(j-1)-sigma_u_vec(j))/(z_vec(j-1)-z_vec(j)))*(z-z_vec(j-1))+sigma_u_vec(j-1)
sigma_v=((sigma_v_vec(j-1)-sigma_v_vec(j))/(z_vec(j-1)-z_vec(j)))*(z-z_vec(j-1))+sigma_v_vec(j-1)
sigma_w=((sigma_w_vec(j-1)-sigma_w_vec(j))/(z_vec(j-1)-z_vec(j)))*(z-z_vec(j-1))+sigma_w_vec(j-1)
end if
end do
! velocity initiation
u=gasdev(seed)*sigma_u
v=gasdev(seed)*sigma_v
w=gasdev(seed)*sigma_w
chi=sigma_chi*gasdev(seed)-0.5*sigma_chi*sigma_chi
up=0.0 ! initiating particle velocity from rest
vp=0.0
wp=0.0
do ! time loop
do j=2,n ! interpolate
if ((z>=z_vec(j-1)).and.(z<=z_vec(j))) then
Um=((Um_vec(j-1)-Um_vec(j))/(z_vec(j-1)-z_vec(j)))*(z-z_vec(j-1))+Um_vec(j-1)
uw=((uw_vec(j-1)-uw_vec(j))/(z_vec(j-1)-z_vec(j)))*(z-z_vec(j-1))+uw_vec(j-1)
duw_dz=((duw_dz_vec(j-1)-duw_dz_vec(j))/(z_vec(j-1)-z_vec(j)))*(z-z_vec(j-1))+duw_dz_vec(j-1)
sigma_u=((sigma_u_vec(j-1)-sigma_u_vec(j))/(z_vec(j-1)-z_vec(j)))*(z-z_vec(j-1))+sigma_u_vec(j-1)
sigma_v=((sigma_v_vec(j-1)-sigma_v_vec(j))/(z_vec(j-1)-z_vec(j)))*(z-z_vec(j-1))+sigma_v_vec(j-1)
sigma_w=((sigma_w_vec(j-1)-sigma_w_vec(j))/(z_vec(j-1)-z_vec(j)))*(z-z_vec(j-1))+sigma_w_vec(j-1)
dvaru_dz=((dvaru_dz_vec(j-1)-dvaru_dz_vec(j))/(z_vec(j-1)-z_vec(j)))*(z-z_vec(j-1))+dvaru_dz_vec(j-1)
dvarv_dz=((dvarv_dz_vec(j-1)-dvarv_dz_vec(j))/(z_vec(j-1)-z_vec(j)))*(z-z_vec(j-1))+dvarv_dz_vec(j-1)
dvarw_dz=((dvarw_dz_vec(j-1)-dvarw_dz_vec(j))/(z_vec(j-1)-z_vec(j)))*(z-z_vec(j-1))+dvarw_dz_vec(j-1)
dissip_m=((dissip_m_vec(j-1)-dissip_m_vec(j))/(z_vec(j-1)-z_vec(j)))*(z-z_vec(j-1))+dissip_m_vec(j-1)
end if
end do
CT = sqrt(1.0+beta*beta*wg*wg/sigma_w/sigma_w) ! crossing trajectories correction
TL = 2.0*sigma_w*sigma_w*C0inv/dissip_m/CT ! added CT effect
TKE = 0.5*(sigma_u*sigma_u+sigma_v*sigma_v+sigma_w*sigma_w)
! -------- Adding dissipation intermittency model --------
omega=dissip_m*CT/TKE ! added CT effect
T_chi=1.0/omega/C_chi
dt=min(dt_on_TL*TL,dt_on_TL*T_chi,dt_max)
a_ln = -(chi + 0.5*sigma_chi*sigma_chi)/T_chi
b_ln = sigma_chi*sqrt(2.0/T_chi)
chi = chi+a_ln*dt+b_ln*sqrt(dt)*gasdev(seed)
dissip_s = dissip_m*exp(chi)
! --------------------------------------------------------
A = 2.0*((sigma_u*sigma_u)*(sigma_w*sigma_w)- uw*uw)
Ainv = 1.0/A
b = sqrt(C0*dissip_s*CT) ! added CT effect
au = (b*b)*(uw*w - u*sigma_w*sigma_w)*Ainv + 0.5*duw_dz &
+ Ainv*(sigma_w*sigma_w*dvaru_dz*u*w - uw*dvaru_dz*w*w &
-uw*duw_dz*u*w + sigma_u*sigma_u*duw_dz*w*w)
av = (-(b*b)*v + dvarv_dz*v*w)/2.0/sigma_v/sigma_v
aw = (b*b)*(uw*u - w*sigma_u*sigma_u)*Ainv + 0.5*dvarw_dz &
+ Ainv*(sigma_w*sigma_w*duw_dz*u*w - uw*duw_dz*w*w &
-uw*dvarw_dz*u*w + sigma_u*sigma_u*dvarw_dz*w*w)
ut = u + au*dt + b*sqrt(dt)*gasdev(seed)
vt = v + av*dt + b*sqrt(dt)*gasdev(seed)
wt = w + aw*dt + b*sqrt(dt)*gasdev(seed)
u = ut
v = vt
w = wt
! -------- Adding IP model --------
if (IP==1) then
dt_ip = dt*0.01
n_ip = 100
upt = up
vpt = vp
wpt = wp
100 do ii=1,n_ip
vr = sqrt((u+Um-upt)*(u+Um-upt)+(v-vpt)*(v-vpt)+(w-wpt)*(w-wpt))
if (vr>1000.0) then
dt_ip = dt_ip*0.5
n_ip = n_ip*2
upt = up
vpt = vp
wpt = wp
goto 100
end if
Re = 2.0*r*vr/nu
Cd = 24.0*(1.0+alpha*Re)/Re
AIP = 3.0*rho*Cd/8.0/rhop/r
gt = g*(rhop - rho)/rhop
upt = upt + AIP*vr*(u+Um-upt)*dt_ip
vpt = vpt + AIP*vr*(v-vpt)*dt_ip
wpt = wpt + (AIP*vr*(w-wpt)-gt)*dt_ip
end do
up = upt
vp = vpt
wp = wpt
end if
! ----------------------------------
if (IP==0) then
up = Um+u
vp = v
wp = w+wg
end if
x = x + up*dt
y = y + vp*dt
z = z + wp*dt
if (i<50) then ! saving trajectories of 50 seeds
write(13,*) i,t,x,y,z
end if
if (z>zmax) then
exit
end if
if (z<zg) then
dt = (z-zg)/(w+wg)
z = z - (w+wg)*dt ! ensure that z = zg at landing
x = x - (u+Um)*dt
y = y - v*dt
write(12,*) i,x,y
exit
end if
dt_max = dz_max/abs(w+wg)
t = t+dt
end do
if (mod(i,100)==0) then
print *, 'wg = ',abs(wg),' zr = ',zs,' pp ',pnum-i
end if
end do
end program LSmodel
!***********************************************************************
! This function generates Gaussian Random Deviates from uniform deviates.
! The function is from Press et al. (1992 p. 280).
function gasdev(idum)
implicit none
integer :: idum, iset
real :: gasdev,fac, gset, rsq, v1, v2, ran
save :: iset, gset
data iset/0/
if (iset.eq.0) then
1 v1=2.*ran(idum)-1.
v2=2.*ran(idum)-1.
rsq=v1**2+v2**2
if (rsq.ge.1. .or. rsq .eq. 0) go to 1
fac =sqrt(-2.*log(rsq)/rsq)
gset=v1*fac
gasdev=v2*fac
iset=1
else
gasdev=gset
iset=0
end if
return
end function gasdev
!***********************************************************************
!uniform random generator between 0 and 1
function ran(idum)
implicit none
integer, parameter :: K4B=selected_int_kind(9)
integer(K4B), intent(inout) :: idum
real :: ran
integer(K4B), parameter :: IA=16807,IM=2147483647,IQ=127773,IR=2836
real, save :: am
integer(K4B), save :: ix=-1,iy=-1,k
if (idum <= 0 .or. iy < 0) then
am=nearest(1.0,-1.0)/IM
iy=ior(ieor(888889999,abs(idum)),1)
ix=ieor(777755555,abs(idum))
idum=abs(idum)+1
end if
ix=ieor(ix,ishft(ix,13))
ix=ieor(ix,ishft(ix,-17))
ix=ieor(ix,ishft(ix,5))
k=iy/IQ
iy=IA*(iy-k*IQ)-IR*k
if (iy < 0) iy=iy+IM
ran=am*ior(iand(IM,ieor(ix,iy)),1)
end function ran
我在Ubuntu中用来编译的命令是
gfortran LSmodel.f95 -o LSmodel.o
没有编译错误，编译正常，但是在 运行之后的程序上，问题开始了。
我在下面的程序 (res.dat) 中包含了 运行 的典型预期输出：
1 21.8583908 8.47351170
2 1.44100714 -8.78142548
3 1154.74109 -265.975677
4 8.41901588 2.71606803
5 84.5189209 -20.4699802
6 86.3176270 -18.4299908
7 133.826874 43.4905090
8 4.37516022 -2.50738835
9 1.31284332 -2.65105081
10 1.96412086 2.85013437
11 4.34823132 -3.83539009
12 40.1227837 -6.60268879
13 3.88699961 2.63749719
14 7.08872795 1.51467562
15 4.72280264 2.63384581
16 0.667112768 1.37209761
17 2.09094667 1.23296225
18 4.72298622 -1.43766475
19 1.04012501 -3.13314247
20 1.91324210 0.957163811
21 1.99065340 0.611227572
22 -2.09086251 -1.41756165
23 -1.46836996 -5.55722380
24 2.41403580 2.18929257E-02
25 3.96990728 -4.91323137
26 1.54687607 -0.527718127
27 8.24332428 -1.48289037
28 4.81600523 -8.87443924
29 2.39538932 0.323360980
30 192.294815 -36.7134209
31 24.6190643 21.7993126
32 -0.124062911 3.44432545
33 16.6237335 -8.54020786
34 50.0964355 -3.29175758
35 5.23409462 2.14592004
36 6.62141275 1.47515869
37 10.7572327 0.307090789
38 63.5973434 -47.7124138
39 74.9621201 2.11509633
40 4.46293068 -1.64074826
41 11.7773390 10.0654907
42 8.26941204 6.84578228
43 0.917451978 2.69560647
44 -2.21521306 15.0752039
45 8.18219483E-02 -2.06250334
46 0.279425710 -3.10328817
47 4.37736464 -1.37771702
48 -2.85058951 -1.79835165
49 5.08391476 2.68537569
50 -4.27199030 -0.642364025
用 gfortran -Wall 编译你的程序得到
test.f90:297:4:
function ran(idum)
1
Warning: 'ran' declared at (1) is also the name of an intrinsic.
It can only be called via an explicit interface or if declared
EXTERNAL. [-Wintrinsic-shadow]
这意味着 user-defined 例程 ran() 与 gfortran 中的内在 ran() 同名。所以我们需要将其声明为外部例程（告诉编译器这是一个 user-defined 例程）：
function gasdev(idum)
implicit none
integer :: idum, iset
real :: gasdev,fac, gset, rsq, v1, v2, ran
save :: iset, gset
data iset/0/
external ran !<--- here
...
有必要在所有使用 user-defined ran 的例程中包含 external ran。 （在这个程序中，只有 gasdev() 例程使用它。）为了避免这种干扰，通常最好使用比 ran、rand 等更具体的名称（例如，rand_uniform() 或 ranranran() 可能不错）。如果该例程包含在模块中以避免此类问题，那将是非常好的，因此请在网上搜索如何使用模块以获取更多详细信息（如果需要...）
我正在尝试编译一个 .f95 fortran 脚本，以便它可以在 Ubuntu 上 运行。该脚本可在此处获得 -> Link to zip file containing .f95 script
当我切换到 Windows 并使用 g95 编译器编译时，它编译并且 运行 没问题。通过 wine 在 Ubuntu 中生成的 .exe 文件也 运行 没问题。
但是，如果我尝试编译生成 Ubuntu 文件，它无法正常工作。我没有收到编译错误，但如果我 运行 生成的文件，要么程序陷入无限循环，要么输出全错。我很难看出哪里出了问题，因为我没有写原始代码，对Fortran的理解也很不稳定，但这似乎与数字计算错误有关，导致非常large/small/inappropriately负输出（抱歉这么含糊）。
我是 运行ning 16.04 xenial ubuntu 和 gfortran 5.4.0。
任何 help/thoughts 感谢这让我感到困惑！谢谢
下面的代码供快速参考：
! Seed dispersal model of of Duman et al. (2015)
! Instructions and expample are found in:
! https://nicholas.duke.edu/people/faculty/katul/research.html
! in 'Library of Functions and Utilities'
! Author: Tomer Duman
! Version: Version 2
! Date: October 22, 2015
! References: Duman, T., Trakhtenbrot, A., Poggi, D.,
! Cassiani, M., Katul, G., Dissipation
! Intermittency Increases Long-Distance
! Dispersal of Heavy Particles in the Canopy
! Sublayer,
! accepcted to Boundary-Layer Meteorology
program LSmodel
implicit none
real :: sec,ran,gasdev ! random generator variables
real :: x,y,z,u,v,w,ut,vt,wt,t,dt ! simulation variables
real :: wg ! seed parametes
real :: Um,sigma_u,sigma_v,sigma_w,uw ! wind statistics variables
real :: dvaru_dz,dvarv_dz,dvarw_dz,duw_dz ! wind statistics variables
real :: dissip_m,TL ! vector over the range of ustars
real :: zs,zg,zmax ! release height & boundaries
real :: Ainv,C0inv ! inverse parameters
real :: C0,A,b,au,av,aw,dt_on_TL ! LS model parameters
real :: dz_max,dt_max ! time step limit
real :: CT,beta ! Crossing Trajectories correction
real :: C_chi,chi,TKE,T_chi,omega ! DI parameters
real :: a_ln,b_ln,sigma_chi,dissip_s ! DI parameters
real :: rhop,rho,r,g,gt,Re,AIP,Cd,nu ! IP parameters
real :: up,vp,wp,upt,vpt,wpt,vr,dt_ip,alpha ! IP parameters
integer :: seed ! random generator variables
integer :: pnum ! simulation parameters
integer :: i,j,jj,n,ii ! counting parameters
integer :: n_ip,IP=1 ! IP parameters
character(len=80) :: filename
real, allocatable,dimension(:) :: z_vec,Um_vec,sigma_u_vec,sigma_v_vec,sigma_w_vec,uw_vec
real, allocatable,dimension(:) :: dvaru_dz_vec,dvarv_dz_vec,dvarw_dz_vec,duw_dz_vec,dissip_m_vec
! setting the random generator seed
seed=7654321
sec=0.0
seed=seed+2*int(secnds(sec))
! input
open (21,file='input_parameters.txt')
read (21, *), x,C0,wg,zs,zg,beta,dt_on_TL,y,sigma_chi,C_chi,r,rhop,alpha,rho,nu
close(21)
pnum = int(x) ! number of released seeds
n = int(y) ! size of the input flow stats
wg = -1.0*wg ! seed terminal velocity [m/s]
!zs ! seed release height [m]
!C0 ! universal constant
!beta ! crossing trajectories parameter
!zg ! ground height [m]
!sigma_chi ! the standard deviation of chi (dissipation intermittency)
!C_chi ! constant for T_chi calc
!r ! particle radium [m] - set to 0 for no IP
!rhop ! particle dry density [kg/m^3]
!alpha ! drag parameter (Cd = 24/Re_p*(1+alpha*Re_p))
!rho ! fluid density [kg/m^3]
!nu ! fluid viscosity [m^2/s]
C0inv = 1.0/C0
g = 9.81
if (r==0.0) then
IP = 0
end if
! limiting parameters to prevent too big jumps in a time-step
dz_max = 0.1
dt_max = 0.1
open(unit=12,file='res.dat', form='formatted')
open(unit=13,file='res_traj.dat', form='formatted')
! allocate
allocate(z_vec(n))
allocate(Um_vec(n))
allocate(sigma_u_vec(n))
allocate(sigma_v_vec(n))
allocate(sigma_w_vec(n))
allocate(uw_vec(n))
allocate(dvaru_dz_vec(n))
allocate(dvarv_dz_vec(n))
allocate(dvarw_dz_vec(n))
allocate(duw_dz_vec(n))
allocate(dissip_m_vec(n))
! load normalized stats
open (22,file='input_flow.txt',form='formatted')
read (22,*) z_vec
read (22,*) Um_vec
read (22,*) sigma_u_vec
read (22,*) sigma_v_vec
read (22,*) sigma_w_vec
read (22,*) uw_vec
read (22,*) dvaru_dz_vec
read (22,*) dvarv_dz_vec
read (22,*) dvarw_dz_vec
read (22,*) duw_dz_vec
read (22,*) dissip_m_vec
close(22)
zmax = z_vec(n)
do i=1,pnum
t=0.0 ! initiate time and location
x=0.0
y=0.0
z=zs
do j=2,n ! interpolate
if ((z>=z_vec(j-1)).and.(z<=z_vec(j))) then
sigma_u=((sigma_u_vec(j-1)-sigma_u_vec(j))/(z_vec(j-1)-z_vec(j)))*(z-z_vec(j-1))+sigma_u_vec(j-1)
sigma_v=((sigma_v_vec(j-1)-sigma_v_vec(j))/(z_vec(j-1)-z_vec(j)))*(z-z_vec(j-1))+sigma_v_vec(j-1)
sigma_w=((sigma_w_vec(j-1)-sigma_w_vec(j))/(z_vec(j-1)-z_vec(j)))*(z-z_vec(j-1))+sigma_w_vec(j-1)
end if
end do
! velocity initiation
u=gasdev(seed)*sigma_u
v=gasdev(seed)*sigma_v
w=gasdev(seed)*sigma_w
chi=sigma_chi*gasdev(seed)-0.5*sigma_chi*sigma_chi
up=0.0 ! initiating particle velocity from rest
vp=0.0
wp=0.0
do ! time loop
do j=2,n ! interpolate
if ((z>=z_vec(j-1)).and.(z<=z_vec(j))) then
Um=((Um_vec(j-1)-Um_vec(j))/(z_vec(j-1)-z_vec(j)))*(z-z_vec(j-1))+Um_vec(j-1)
uw=((uw_vec(j-1)-uw_vec(j))/(z_vec(j-1)-z_vec(j)))*(z-z_vec(j-1))+uw_vec(j-1)
duw_dz=((duw_dz_vec(j-1)-duw_dz_vec(j))/(z_vec(j-1)-z_vec(j)))*(z-z_vec(j-1))+duw_dz_vec(j-1)
sigma_u=((sigma_u_vec(j-1)-sigma_u_vec(j))/(z_vec(j-1)-z_vec(j)))*(z-z_vec(j-1))+sigma_u_vec(j-1)
sigma_v=((sigma_v_vec(j-1)-sigma_v_vec(j))/(z_vec(j-1)-z_vec(j)))*(z-z_vec(j-1))+sigma_v_vec(j-1)
sigma_w=((sigma_w_vec(j-1)-sigma_w_vec(j))/(z_vec(j-1)-z_vec(j)))*(z-z_vec(j-1))+sigma_w_vec(j-1)
dvaru_dz=((dvaru_dz_vec(j-1)-dvaru_dz_vec(j))/(z_vec(j-1)-z_vec(j)))*(z-z_vec(j-1))+dvaru_dz_vec(j-1)
dvarv_dz=((dvarv_dz_vec(j-1)-dvarv_dz_vec(j))/(z_vec(j-1)-z_vec(j)))*(z-z_vec(j-1))+dvarv_dz_vec(j-1)
dvarw_dz=((dvarw_dz_vec(j-1)-dvarw_dz_vec(j))/(z_vec(j-1)-z_vec(j)))*(z-z_vec(j-1))+dvarw_dz_vec(j-1)
dissip_m=((dissip_m_vec(j-1)-dissip_m_vec(j))/(z_vec(j-1)-z_vec(j)))*(z-z_vec(j-1))+dissip_m_vec(j-1)
end if
end do
CT = sqrt(1.0+beta*beta*wg*wg/sigma_w/sigma_w) ! crossing trajectories correction
TL = 2.0*sigma_w*sigma_w*C0inv/dissip_m/CT ! added CT effect
TKE = 0.5*(sigma_u*sigma_u+sigma_v*sigma_v+sigma_w*sigma_w)
! -------- Adding dissipation intermittency model --------
omega=dissip_m*CT/TKE ! added CT effect
T_chi=1.0/omega/C_chi
dt=min(dt_on_TL*TL,dt_on_TL*T_chi,dt_max)
a_ln = -(chi + 0.5*sigma_chi*sigma_chi)/T_chi
b_ln = sigma_chi*sqrt(2.0/T_chi)
chi = chi+a_ln*dt+b_ln*sqrt(dt)*gasdev(seed)
dissip_s = dissip_m*exp(chi)
! --------------------------------------------------------
A = 2.0*((sigma_u*sigma_u)*(sigma_w*sigma_w)- uw*uw)
Ainv = 1.0/A
b = sqrt(C0*dissip_s*CT) ! added CT effect
au = (b*b)*(uw*w - u*sigma_w*sigma_w)*Ainv + 0.5*duw_dz &
+ Ainv*(sigma_w*sigma_w*dvaru_dz*u*w - uw*dvaru_dz*w*w &
-uw*duw_dz*u*w + sigma_u*sigma_u*duw_dz*w*w)
av = (-(b*b)*v + dvarv_dz*v*w)/2.0/sigma_v/sigma_v
aw = (b*b)*(uw*u - w*sigma_u*sigma_u)*Ainv + 0.5*dvarw_dz &
+ Ainv*(sigma_w*sigma_w*duw_dz*u*w - uw*duw_dz*w*w &
-uw*dvarw_dz*u*w + sigma_u*sigma_u*dvarw_dz*w*w)
ut = u + au*dt + b*sqrt(dt)*gasdev(seed)
vt = v + av*dt + b*sqrt(dt)*gasdev(seed)
wt = w + aw*dt + b*sqrt(dt)*gasdev(seed)
u = ut
v = vt
w = wt
! -------- Adding IP model --------
if (IP==1) then
dt_ip = dt*0.01
n_ip = 100
upt = up
vpt = vp
wpt = wp
100 do ii=1,n_ip
vr = sqrt((u+Um-upt)*(u+Um-upt)+(v-vpt)*(v-vpt)+(w-wpt)*(w-wpt))
if (vr>1000.0) then
dt_ip = dt_ip*0.5
n_ip = n_ip*2
upt = up
vpt = vp
wpt = wp
goto 100
end if
Re = 2.0*r*vr/nu
Cd = 24.0*(1.0+alpha*Re)/Re
AIP = 3.0*rho*Cd/8.0/rhop/r
gt = g*(rhop - rho)/rhop
upt = upt + AIP*vr*(u+Um-upt)*dt_ip
vpt = vpt + AIP*vr*(v-vpt)*dt_ip
wpt = wpt + (AIP*vr*(w-wpt)-gt)*dt_ip
end do
up = upt
vp = vpt
wp = wpt
end if
! ----------------------------------
if (IP==0) then
up = Um+u
vp = v
wp = w+wg
end if
x = x + up*dt
y = y + vp*dt
z = z + wp*dt
if (i<50) then ! saving trajectories of 50 seeds
write(13,*) i,t,x,y,z
end if
if (z>zmax) then
exit
end if
if (z<zg) then
dt = (z-zg)/(w+wg)
z = z - (w+wg)*dt ! ensure that z = zg at landing
x = x - (u+Um)*dt
y = y - v*dt
write(12,*) i,x,y
exit
end if
dt_max = dz_max/abs(w+wg)
t = t+dt
end do
if (mod(i,100)==0) then
print *, 'wg = ',abs(wg),' zr = ',zs,' pp ',pnum-i
end if
end do
end program LSmodel
!***********************************************************************
! This function generates Gaussian Random Deviates from uniform deviates.
! The function is from Press et al. (1992 p. 280).
function gasdev(idum)
implicit none
integer :: idum, iset
real :: gasdev,fac, gset, rsq, v1, v2, ran
save :: iset, gset
data iset/0/
if (iset.eq.0) then
1 v1=2.*ran(idum)-1.
v2=2.*ran(idum)-1.
rsq=v1**2+v2**2
if (rsq.ge.1. .or. rsq .eq. 0) go to 1
fac =sqrt(-2.*log(rsq)/rsq)
gset=v1*fac
gasdev=v2*fac
iset=1
else
gasdev=gset
iset=0
end if
return
end function gasdev
!***********************************************************************
!uniform random generator between 0 and 1
function ran(idum)
implicit none
integer, parameter :: K4B=selected_int_kind(9)
integer(K4B), intent(inout) :: idum
real :: ran
integer(K4B), parameter :: IA=16807,IM=2147483647,IQ=127773,IR=2836
real, save :: am
integer(K4B), save :: ix=-1,iy=-1,k
if (idum <= 0 .or. iy < 0) then
am=nearest(1.0,-1.0)/IM
iy=ior(ieor(888889999,abs(idum)),1)
ix=ieor(777755555,abs(idum))
idum=abs(idum)+1
end if
ix=ieor(ix,ishft(ix,13))
ix=ieor(ix,ishft(ix,-17))
ix=ieor(ix,ishft(ix,5))
k=iy/IQ
iy=IA*(iy-k*IQ)-IR*k
if (iy < 0) iy=iy+IM
ran=am*ior(iand(IM,ieor(ix,iy)),1)
end function ran
我在Ubuntu中用来编译的命令是
gfortran LSmodel.f95 -o LSmodel.o
没有编译错误，编译正常，但是在 运行之后的程序上，问题开始了。
我在下面的程序 (res.dat) 中包含了 运行 的典型预期输出：
1 21.8583908 8.47351170
2 1.44100714 -8.78142548
3 1154.74109 -265.975677
4 8.41901588 2.71606803
5 84.5189209 -20.4699802
6 86.3176270 -18.4299908
7 133.826874 43.4905090
8 4.37516022 -2.50738835
9 1.31284332 -2.65105081
10 1.96412086 2.85013437
11 4.34823132 -3.83539009
12 40.1227837 -6.60268879
13 3.88699961 2.63749719
14 7.08872795 1.51467562
15 4.72280264 2.63384581
16 0.667112768 1.37209761
17 2.09094667 1.23296225
18 4.72298622 -1.43766475
19 1.04012501 -3.13314247
20 1.91324210 0.957163811
21 1.99065340 0.611227572
22 -2.09086251 -1.41756165
23 -1.46836996 -5.55722380
24 2.41403580 2.18929257E-02
25 3.96990728 -4.91323137
26 1.54687607 -0.527718127
27 8.24332428 -1.48289037
28 4.81600523 -8.87443924
29 2.39538932 0.323360980
30 192.294815 -36.7134209
31 24.6190643 21.7993126
32 -0.124062911 3.44432545
33 16.6237335 -8.54020786
34 50.0964355 -3.29175758
35 5.23409462 2.14592004
36 6.62141275 1.47515869
37 10.7572327 0.307090789
38 63.5973434 -47.7124138
39 74.9621201 2.11509633
40 4.46293068 -1.64074826
41 11.7773390 10.0654907
42 8.26941204 6.84578228
43 0.917451978 2.69560647
44 -2.21521306 15.0752039
45 8.18219483E-02 -2.06250334
46 0.279425710 -3.10328817
47 4.37736464 -1.37771702
48 -2.85058951 -1.79835165
49 5.08391476 2.68537569
50 -4.27199030 -0.642364025
用 gfortran -Wall 编译你的程序得到
test.f90:297:4:
function ran(idum)
1
Warning: 'ran' declared at (1) is also the name of an intrinsic.
It can only be called via an explicit interface or if declared
EXTERNAL. [-Wintrinsic-shadow]
这意味着 user-defined 例程 ran() 与 gfortran 中的内在 ran() 同名。所以我们需要将其声明为外部例程（告诉编译器这是一个 user-defined 例程）：
function gasdev(idum)
implicit none
integer :: idum, iset
real :: gasdev,fac, gset, rsq, v1, v2, ran
save :: iset, gset
data iset/0/
external ran !<--- here
...
有必要在所有使用 user-defined ran 的例程中包含 external ran。 （在这个程序中，只有 gasdev() 例程使用它。）为了避免这种干扰，通常最好使用比 ran、rand 等更具体的名称（例如，rand_uniform() 或 ranranran() 可能不错）。如果该例程包含在模块中以避免此类问题，那将是非常好的，因此请在网上搜索如何使用模块以获取更多详细信息（如果需要...）