在cygwin / gfortran下编译一个Fortran程序给出了“_WinMain @ 16”的未定义引用“ [英] Compiling a Fortran program under cygwin/gfortran gives "undefined reference to `_WinMain@16'"
本文介绍了在cygwin / gfortran下编译一个Fortran程序给出了“_WinMain @ 16”的未定义引用“的处理方法，对大家解决问题具有一定的参考价值，需要的朋友们下面随着小编来一起学习吧！
问题描述
在Cygwin终端中输入
$ gfortran -o threed_euler_fluxes_v3.exe threed_euler_fluxes_v3.f90
并且我得到编译器错误
/usr/lib/gcc/i686-pc-cygwin/4.5.3 /../../../ libcygwin.a（libcmain.o）：在函数`main'中：
/ usr / src / debug / cygwin-1.7.17-1 / winsup / cygwin / lib / libcmain.c：39：对`_WinMain @ 16'的未定义引用
collect2：ld返回1退出状态
我也试过像这样编译
$ gfortran -o threed_euler_fluxes_v3.exe threed_euler_fluxes_v3.f90 -shared
但是当我尝试运行时，我得到了一个错误说它不是一个有效的Windows程序？
下面是完整的Fortran代码。我删除了一些评论，以便将字数限制在30k以下。这里是原创。
$ $ $
$ integer $ parameter $参数： p2 = selected_real_kind（15）！双精度
！输入
real（p2），intent（in）:: primL（5），primR（5）！输入：原始变量
real（p2），intent（in）:: njk（3）！输入：面向法向量
！输出
实数（p2），意图（输出）:: num_flux（5）！输出：数值通量
！一些常数
real（p2）:: zero = 0.0_p2
real（p2）:: one = 1.0_p2
real（p2 ）:: two = 2.0_p2
real（p2）:: half = 0.5_p2
real（p2）:: fifth = 0.2_p2
！局部变量
真实（p2）:: nx，ny，nz！法线向量
real（p2）:: mx，my，mz！正交切向量
real（p2）:: lx，ly，lz！另一个正交切线向量
real（p2）:: abs_n_cross_l！ （p2）:: uL，uR，vL，vR，wL，wR！速度组件。
real（p2）:: rhoL，rhoR，pL，pR！原始变量。
real（p2）:: qnL，qnR，qmL，qmR，qlL，qlR！正常速度和切线速度
real（p2）:: aL，aR，HL，HR！声速，总焓
real（p2）:: RT，rho，u，v，w，H，a，qn，ql，qm！ Roe-平均值
real（p2）:: drho，dqn，dql，dqm，dp，LdU（5）！波浪强度
real（p2）:: ws（5），R（5,5）！波速和右特征向量
real（p2）:: dws（5）！用于熵修正的抛物线拟合宽度
real（p2）:: fL（5），fR（5），diss（5）！通量耗散项
real（p2）:: gamma = 1.4_p2！比热比
real（p2）:: temp，tempx，tempy，tempz！ Temoprary变量
！面法向量（单位向量）
nx = njk（1）
ny = njk（2）
nz = njk（3）
tempx = ny * ny + nz * nz
tempy = nz * nz + nx * nx
tempz = nx * nx + ny * ny
如果（tempx> = tempy。并且tempx> = tempz），则
lx =零
ly = -nz
lz = ny
elseif（tempy> = tempx然后
lx = -nz
ly =零
lz = nx
elseif（tempz> = tempx .and tempz> = tempz） tempy）then
lx = -ny
ly = nx
lz = zero
else
！不可能发生
写入（*，*）子程序inviscid_roe：不可能发生，请报告问题。
stop
endif
！使它成为单位矢量。
temp = sqrt（lx * lx + ly * ly + lz * lz）
lx = lx / temp
ly = ly / temp
lz = lz / temp
mx = ny * lz - nz * ly
my = nz * lx - nx * lz
mz = nx * ly - ny * lx
abs_n_cross_l = sqrt（mx ** 2 + my ** 2 + mz ** 2）
mx = mx / abs_n_cross_l
my = my / abs_n_cross_l
mz = mz / abs_n_cross_l
rhoL = primL（1）
uL = primL（2）
vL = primL（3）
wL = primL（4）
qnL = uL * nx + vL * ny + wL * nz
qlL = uL * lx + vL * ly + wL * lz
qmL = uL * mx + vL * my + wL * mz
pL = primL（5）
aL = sqrt（gamma * pL / rhoL）
HL = aL * aL /（gamma-one）+ half *（uL * uL + vL * vL + wL * wL）
！右状态
rhoR = primR（1）
uR = primR（2）
vR = primR（3）
wR = primR（4）
qnR = uR * nx + vR * ny + wR * nz
qlR = uR * lx + vR * ly + wR * lz
qmR = uR * mx + vR * my + wR * mz
pR = primR （5）
aR = sqrt（gamma * pR / rhoR）
HR = aR * aR /（gamma-one）+ half *（uR * uR + vR * vR + wR * wR）
RT = sqrt（rhoR / rhoL）
rho = RT * rhoL！Roe平均密度
u =（uL + RT * uR）/（一个RT）！Roe平均的X速度
v =（vL + RT * vR）/（一+ RT）！Roe-平均的y速度
w =（wL + RT * wR）/（一（+ RT））！Roe平均z湍流
H =（HL + RT * HR）/（一+ RT）！Roe平均总焓
a = sqrt（（gamma-one）*（ H-half *（u * u + v * v + w * w）））！Roe平均的声音速度
qn = u * nx + v * ny + w * nz！ mal velocity velocity mal mal b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b b =切线速度
！波动强度
drho = rhoR - rhoL！密度差
dp = pR - pL！压差
dqn = qnR - qnL ！正常速度差
dql = qlR - qlL！正切速度差l
dqm = qmR - qmL！正切速度差m
LdU（1）=（dp - rho * a * dqn）/（two * a * a）！左移声波强度
LdU（2）= drho - dp /（a * a）！熵波强度
LdU ）=（dp + rho * a * dqn）/（two * a * a）！右移声波强度
LdU（4）= rho * dql！剪波强度
LdU（5） = rho * dqm！剪波强度
！波速的绝对值
ws（1）= abs（qn-a）！左移声波速度
ws（2） = abs（qn）！熵波速度
ws（3）= abs（qn + a）！右移声波速度
ws（4）= abs（qn）！剪波速度
ws（5）= abs（qn）！剪波速度
！Harten的Entropy Fix JCP（1983），49，pp357-393：只适用于非线性场。
！注意：它通过在ws = 0附近进行抛物线拟合来避免波速消失。
dws（1）= five
if（ws（1）
if（ws（3） （3））ws（3）= half *（ws（3）* ws（3）/ dws（3）+ dws（3））
！Right特征向量
！左移声波
R（1,1）= 1
R（2,1）= u - a * nx
R（3,1）= v - a * ny
R（4,1）= w - a * nz
R（5,1）= H - a * qn
！ （3,2）= v
R（4,2）= w（2,2）= u
$
R（5,2）= half *（u * u + v * v + w * w）
！右移声波
R（1,3）=一
R（2,3）= u + a * nx
R（3,3）= v + a * ny $ （5）= H + a * qn
！剪切波
R（1,4）=零
R（2,4）= 1x
R（3,4）= ly
R（4,4）= lz
R（5,4）= ql
！剪切波
R（1,5）=零
R（2,5）= mx
R（3,5）= my
R（4,5）= mz
R（5,5）= qm
diss（:) = ws（1）* LdU（1）* R（：，1）+ ws（2）* LdU（2）* R（：2）+ ws（3）* LdU（3）* R（：3）&
+ ws（4）* LdU（4）* R（：4）+ ws（5）* LdU（5）* R（：，5）
fL（1）= rhoL * qnL
fL（2）= rhoL * qnL * uL + pL * nx
fL（3）= rhoL * qnL * vL + pL * ny $ b （5）= rhoL * qnL * HL
fR（1）= rhoR * qnR
（4）= rhoL * qnL * wL + pL * nz
fL fR（2）= rhoR * qnR * uR + pR * nx
fR（3）= rhoR * qnR * vR + pR * ny
fR（4）= rhoR * qnR * wR + pR * nz
fR（5）= rhoR * qnR * HR
num_flux = half *（fL + fR - diss）
子程序inviscid_roe_n（primL，primR，njk，num_flux）
隐式无
整型，参数:: p2 = selected_real_kind（15）！双精度
！输入
real（p2），intent（in）:: primL（5），primR（5）！输入：原始变量
real（p2），intent（in）:: njk（3）！输入：面向法向量
！输出
实数（p2），意图（输出）:: num_flux（5）！输出：数值通量
！一些常数
real（p2）:: zero = 0.0_p2
real（p2）:: one = 1.0_p2
real（p2 ）:: two = 2.0_p2
real（p2）:: half = 0.5_p2
real（p2）:: fifth = 0.2_p2
！局部变量
真实（p2）:: nx，ny，nz！法向量
real（p2）:: uL，uR，vL，vR，wL，wR！速度组件。
real（p2）:: rhoL，rhoR，pL，pR！原始变量。
real（p2）:: qnL，qnR！正常速度
实数（p2）:: aL，aR，HL，HR！声速，总焓
real（p2）:: RT，rho，u，v，w，H，a，qn！ Roe-平均值
real（p2）:: drho，dqn，dp，LdU（4）！波浪强度
real（p2）:: du，dv，dw！速度差异
real（p2）:: ws（4），R（5,4）！波速和右特征向量
real（p2）:: dws（4）！用于熵修正的抛物线拟合宽度
real（p2）:: fL（5），fR（5），diss（5）！通量耗散项
real（p2）:: gamma = 1.4_p2！比热比
！面法向量（单位向量）
nx = njk（1）
ny = njk（2）
nz = njk（3）
！原始变量和其他变量。
！左状态
rhoL = primL（1）
uL = primL（2）
vL = primL（3）
wL = primL（4）
qnL = uL * nx + vL * ny + wL * nz
pL = primL（5）
aL = sqrt（gamma * pL / rhoL）
HL = aL * aL /（gamma-one）+ half *（uL * uL + vL * vL + wL * wL）
！右状态
rhoR = primR（1）
uR = primR（2）
vR = primR（3）
wR = primR（4）
qnR = uR * nx + vR * ny + wR * nz
pR = primR（5）
aR = sqrt（gamma * pR / rhoR）
HR = aR * aR /（gamma-one）+ half *（uR * uR + vR * vR + wR * wR）
！首先计算Roe平均数量
！注意：有关
的信息，请参阅http://www.cfdnotes.com/cfdnotes_roe_averaged_density.html！ Roe平均密度。
RT = sqrt（rhoR / rhoL）
rho = RT * rhoL！Roe平均密度
u =（uL + RT * uR）/（one + RT）！Roe平均速度b bw =（wL + RT * wR）/（1 + RT）！Roe平均速度b bw =（wL + RT * wR）/（1 + RT）平均z速度
H =（HL + RT * HR）/（1 + RT）！Roe平均总焓
a = sqrt（（gamma-1）*（H-half *（u * u + v * v + w * w）））！Roe平均声速
qn = u * nx + v * ny + w * nz！ b $ b！波动强度
drho = rhoR - rhoL！密度差
dp = pR - pL！压差
dqn = qnR - qnL！正常速度差$ b $ （2）=（dp-rho * a * dqn）/（2 * a * a）！左移声波强度
LdU（2）= drho - dp / a）！波浪强度
LdU（ 3）=（dp + rho * a * dqn）/（2 * a * a）！右移声波强度
LdU（4）=ρ！剪切波强度（不是真的，只是一个因子）
！波的绝对值速度
ws（1）= abs（qn-a）！左移声波
ws（2）= abs（ （3）= abs（qn + a）！右移声波
ws（4）= abs（qn）！剪波
！Harten的Entropy Fix JCP（1983），49，pp357-393：仅适用于非线性领域。
！注意：它通过在ws = 0附近进行抛物线拟合来避免波速消失。
dws（1）= five
if（ws（1）
if（ws（3） （3）= dws（3））ws（3）= half *（ws（3）* ws（3）/ dws（3）+ dws（3））
R 1,1）= 1
R（2,1）= u - a * nx
R（3,1）= v - a * ny
R（4,1）= w - a * nz
R（5,1）= H - a * qn
R（1,2）=一
R（2,2） = u
R（3,2）= v
R（4,2）= w
R（5,2）= half *（u * u + v * v + w * w）
！右移声波
R（1,3）=一
R（2,3）= u + a * nx
R（3,3）= v + a * ny $ （5）= H + a * qn
！两个剪切波分量合并为一个（波强度合并）。
du = uR - uL
dv = vR - vL
dw = wR - wL
R（1,4）=零
R（2,4）= du = dqn * nx
R（3,4）= dv - dqn * ny
R（4,4）= dw - dqn * nz
R（5,4）= u *耗散项：| An |（UR-UL）= R |λ| L * dU = [ws（k）* sum_k） R（：，k）* L * dU（k）]
diss（:) = ws（1）* LdU（1）* R（：，1）+ ws（2）* LdU （2）* R（：，2）&
+ ws（3）* LdU（3）* R（：3）+ ws（4）* LdU（4）* R（：4）
！计算物理通量：fL = Fn（UL）和fR = Fn（UR）
fL（1）= rhoL * qnL
fL（2）= rhoL * qnL * uL + pL * nx
fL（3）= rhoL * qnL * vL + pL * ny
fL（4）= rhoL * qnL * wL + pL * nz
fL（5）= rhoL * qnL * HL $ （1）= rhoR * qnR
fR（2）= rhoR * qnR * uR + pR * nx
fR（3）= rhoR * qnR * vR + pR * （5）= rhoR * qnR * HR
！这是数字通量：Roe通量= 1/2 * [Fn（UL）+ Fn（UR） - | An |（UR-UL）]
num_flux = half *（fL + fR - diss）
！法线方向上的正常最大波速。
！ wsn = abs（qn）+ a
结束子程序inviscid_roe_n
子程序inviscid_rotated_rhll（primL，primR，njk，num_flux）
隐式无
整数，参数:: p2 = selected_real_kind（15）！双精度
！输入
real（p2），intent（in）:: primL（5），primR（5）！输入：原始变量
real（p2），intent（in）:: njk（3）！输入：面向法向量
！输出
实数（p2），意图（输出）:: num_flux（5）！输出：数值通量
！一些常数
real（p2）:: zero = 0.0_p2
real（p2）:: one = 1.0_p2
real（p2 ）:: two = 2.0_p2
real（p2）:: half = 0.5_p2
real（p2）:: fifth = 0.2_p2
！局部变量
真实（p2）:: nx，ny，nz！面对法向量
real（p2）:: uL，uR，vL，vR，wL，wR！速度组件。
real（p2）:: rhoL，rhoR，pL，pR！原始变量。
real（p2）:: qnL，qnR！正常速度
实数（p2）:: aL，aR，HL，HR！声速，总焓
real（p2）:: RT，rho，u，v，w，H，a，qn！ Roe-平均值
real（p2）:: drho，dqn，dp，LdU（4）！波浪强度
real（p2）:: du，dv，dw！速度成分差异
real（p2）:: eig（4）！特征值
real（p2）:: ws（4），R（5,4）！绝对波速度和右特征向量
real（p2）:: dws（4）！用于熵修正的抛物线拟合宽度
real（p2）:: fL（5），fR（5），diss（5）！通量耗散项
real（p2）:: gamma = 1.4_p2！比热比
real（p2）:: SRp，SLm！ HLL部分的波速
real（p2）:: nx1，ny1，nz1！沿着哪个向量应用HLL
real（p2）:: nx2，ny2，nz2！沿着哪个向量应用Roe
real（p2）:: alpha1，alpha2！新法线的预测
real（p2）:: abs_dq！速度差的大小
real（p2）:: temp，tempx，tempy，tempz！临时变量
！面法向量（单位向量）
nx = njk（1）
ny = njk（2）
nz = njk（3）
！原始变量和其他变量。
！左状态
rhoL = primL（1）
uL = primL（2）
vL = primL（3）
wL = primL（4）
qnL = uL * nx + vL * ny + wL * nz
pL = primL（5）
aL = sqrt（gamma * pL / rhoL）
HL = aL * aL /（gamma-one）+ half *（uL * uL + vL * vL + wL * wL）
！右状态
rhoR = primR（1）
uR = primR（2）
vR = primR（3）
wR = primR（4）
qnR = uR * nx + vR * ny + wR * nz
pR = primR（5）
aR = sqrt（gamma * pR / rhoR）
HR = aR * aR /（gamma-one）+ half *（uR * uR + vR * vR + wR * wR）
！计算物理通量：fL = Fn（UL）和fR = Fn（UR）
fL （1）= rhoL * qnL
fL（2）= rhoL * qnL * uL + pL * nx
fL（3）= rhoL * qnL * vL + pL * ny
fL ）= rhoL * qnL * wL + pL * nz
fL（5）= rhoL * qnL * HL
fR（1）= rhoR * qnR
fR（2）= rHR * qnR * uR + pR * nx
fR（3）= rhoR * qnR * vR + pR * ny
fR（4）= rhoR * qnR * wR + pR * nz
fR （5）= rhoR * qnR * HR
abs_dq = sqrt（（uR-uL）** 2 +（vR-vL）** 2 +（wR- wL）** 2）
if（abs_dq> 1.0e-12_p2）then
nx1 =（uR-uL）/ abs_dq
ny1 =（vR-vL）/ abs_dq
nz1 =（wR-wL）/ abs_dq
tempx = ny * ny + nz * nz
tempy = nz * nz + nx * nx
tempz = nx * nx + ny * ny
if（tempx> = tempy .and。 tempx> = tempz），则
nx1 =零
ny1 = -nz
nz1 = ny
elseif（tempy> = tempx .and tempy> = tempz）then
nx1 = -nz
ny1 =零
nz1 = nx
elseif（tempz> = tempx .and tempz> = tempy）then
nx1 = - ny
ny1 = nx
nz1 =零
else
！不可能发生
写（*，*）inviscid_rotated_rhll：不可能发生，请报告问题。
stop
endif
！使它成为单位矢量。
temp = sqrt（nx1 * nx1 + ny1 * ny1 + nz1 * nz1）
nx1 = nx1 / temp
ny1 = ny1 / temp
nz1 = nz1 / temp
endif
alpha1 = nx * nx1 + ny * ny1 + nz * nz1
！使alpha1总是正面的。
temp = sign（one，alpha1）
nx1 = temp * nx1
ny1 = temp * ny1
nz1 = temp * nz1
alpha1 = temp * alpha1
！n2 =垂直于n1的方向。
！注意：这个矢量有无限多的选择。
！未来可能会发现最好的选择。
！在这里，我们在文中使用公式（4.4）：
！ （nx2，ny2，nz2）=（n1xn）xn1 / |（n1xn）xn1 | （'x'是矢量产品。）
！ （tempx，tempy，tempz）= n1xn
tempx = ny1 * nz - nz1 * ny
tempy = nz1 * nx - nx1 * nz
tempz = nx1 * ny - ny1 * nx
！ （nx2，ny2，nz2）=（n1xn）xn1
nx2 = tempy * nz1 - tempz * ny1
ny2 = tempz * nx1 - tempx * nz1
nz2 = tempx * ny1 - tempy * nx1
！使n2为单位向量
temp = sqrt（nx2 * nx2 + ny2 * ny2 + nz2 * nz2）
nx2 = nx2 / temp
ny2 = ny2 / temp
nz2 = nz2 / temp
alpha2 = nx * nx2 + ny * ny2 + nz * nz2
！让alpha2总是正面的。
temp = sign（one，alpha2）
nx2 = temp * nx2
ny2 = temp * ny2
nz2 = temp * nz2
alpha2 = temp * alpha2
！------------------------------------------- -------------------------------------
！现在我们要计算Roe以n2作为正常的流量并修改
！波速（5.12）。注意：这里的Roe通量是不带切线向量计算的。
！详情请参阅我喜欢CFD，VOL.1：第57页，方程（3.6.31）。
！首先计算Roe平均数量
！注意：有关
的信息，请参阅http://www.cfdnotes.com/cfdnotes_roe_averaged_density.html！ Roe平均密度。
RT = sqrt（rhoR / rhoL）
rho = RT * rhoL！Roe平均密度。
u =（uL + RT * uR）/（one + RT）！Roe-averaged x-velocity
v =（vL + RT * vR）/（one + RT）！Roe-averaged y-velocity
w =（wL + RT * wR）/（one + RT）！Roe平均z速度
H =（HL + RT * HR）/（1 + RT）！Roe-平均总焓
a = sqrt（（gamma-one）*（H-half *（u * u + v * v + w * w）））！Roe平均声音速度
！ - -------------------------------------------------- -
！计算HLL部分的波速估计值，
！跟随Einfeldt：
！
！ B. Einfeldt，关于Godunov型气体动力学方法，
！ SIAM数值分析杂志25（2）（1988）294-318。
！
！注意：HLL实际上应用于n1，但这是
！我们需要整合HLL。参见JCP2008论文。
qn = u * nx1 + v * ny1 + w * nz1
qnL = uL * nx1 + vL * ny1 + wL * nz1
qnR = uR * nx1 + vR * ny1 + wR * nz1
SLm = min（零，qn - a，qnL - aL）！最小波速估计值
SRp = max（zero，qn + a，qnR + aR）估计
！这是唯一使用n1 =（nx1，ny1，nz1）的地方。
！下面从未使用n1 =（nx1，ny1，nz1）。
！--------------------------------------------- -------
！波动强度
qn = u * nx2 + v * ny2 + w * nz2
qnL = uL * nx2 + vL * ny2 + wL * nz2
qnR = uR * nx2 + vR * ny2 + wR * nz2
drho = rhoR - rhoL！密度差
dp = pR - pL！压力差值
dqn = qnR - qnL！正常速度差
LdU（1）=（dp-rho * a * dqn）/（two * a * a）！左移声波强度
LdU（2）= drho - dp /（a * a）！熵波强度
LdU（3）=（dp + rho * a * dqn）/（two * a * a）！右移声波强度
LdU（4）= rho！剪波强度（不是真的，只是一个因子）
！波速（特征值）
eig（1）= qn-a！左移声波速度
eig（2）= qn！熵波速度
eig（3）= qn + a！右移声波速度
eig（4）= qn！剪波速度
！波速的绝对值（Eige （1）= abs（qn-a）！左移声波速度
ws（2）= abs（qn）！熵波速度
ws （3）= abs（qn + a）！右移声波速度
ws（4）= abs（qn）！剪波速度
！Harten的Entropy Fix JCP（1983） ，49，pp357-393：只适用于非线性领域。
！注意：它通过在ws = 0附近进行抛物线拟合来避免波速消失。
dws（1）= five
if（ws（1）
if（ws（3） （3））ws（3）= half *（ws（3）* ws（3）/ dws（3）+ dws（3））
！将波速Rotated-RHLL：原始JCP2008论文中的公式（5.12）。
ws = alpha2 * ws - （alpha1 * two * SRp * SLm + alpha2 *（SRp + SLm）* eig）/（SRp-SLm）
！我们用上述修正波速计算方向为n2
！的Roe耗散项。 HLL波速的作用类似于
！熵固定或特征值限制;他们只贡献
！给出的分数α1（小于或等于1.0）。参见JCP2008论文。
！右特征向量：
！注意：两个剪切波分量合并为一个，所以切向量
！不是必需的。这就是为什么这里只有4个向量。
！左移声波
R（1,1）= 1
R（2,1）= u - a * nx2
R（3,1）= v - a * ny2
R（4,1）= w - a * nz2
R（5,1）= H - a * qn
！ （3,2）= v
R（4,2）= w（2,2）= u
$
R（5,2）= half *（u * u + v * v + w * w）
！右移声波
R（1,3）=一
R（2,3）= u + a * nx2
R（3,3）= v + a * ny2
R（4,3）= w + a * nz2
R（5,3）= H + a * qn
！两个剪切波分量合并为一个（波强度合并）。
du = uR - uL
dv = vR - vL
dw = wR - wL
R（1,4）=零
R（2,4）= （4,4）= dw - dqn * nz2
R（5,4）= u * du + v * dv + w * dw - qn * dqn
！耗散项：Roe耗散与修改的波速。
！ | w |（k）* R（：，k）* L * dU（k）]，其中n = n2，| An | dU = R | L * dU = sum_k。 （1）* LdU（1）* R（：，1）+ ws（2）* LdU（2）* R（：2）&
+ ws（3）* LdU（3）* R（：3）+ ws（4）* LdU（4）* R（：4）
！计算旋转-RHLL通量。 （它看起来像带有Roe耗散的HLL流量）。
num_flux =（SRp * fL - SLm * fR）/（SRp-SLm） - half * diss
！正常方向上的最大正常波速。
！ wsn = abs（qn）+ a
结束子程序inviscid_rotated_rhll
！-------------------------- -------------------------------------------------- ----
解决方案
您的文件不是所有！它是一个子程序的集合。你不能编译它作为程序运行，只能作为一个目标文件或一个库（试试 -c 或 -shared ）。您必须添加主程序主体才能将其编译为程序并运行它！
In the Cygwin terminal I enter
$ gfortran -o threed_euler_fluxes_v3.exe threed_euler_fluxes_v3.f90
and I get the compiler error
/usr/lib/gcc/i686-pc-cygwin/4.5.3/../../../libcygwin.a(libcmain.o): In function `main':
/usr/src/debug/cygwin-1.7.17-1/winsup/cygwin/lib/libcmain.c:39: undefined reference to `_WinMain@16'
collect2: ld returned 1 exit status
I also tried compiling like this
$ gfortran -o threed_euler_fluxes_v3.exe threed_euler_fluxes_v3.f90 -shared
but when I tried running I got an error saying it wasn't a valid windows program?
Heres the complete fortran code. I removed some comments inorder to keep word limit below 30k Heres the original.
subroutine inviscid_roe(primL, primR, njk, num_flux)
implicit none
integer , parameter :: p2 = selected_real_kind(15) ! Double precision
!Input
real(p2), intent( in) :: primL(5), primR(5) ! Input: primitive variables
real(p2), intent( in) :: njk(3) ! Input: face normal vector
!Output
real(p2), intent(out) :: num_flux(5) ! Output: numerical flux
!Some constants
real(p2) :: zero = 0.0_p2
real(p2) :: one = 1.0_p2
real(p2) :: two = 2.0_p2
real(p2) :: half = 0.5_p2
real(p2) :: fifth = 0.2_p2
!Local variables
real(p2) :: nx, ny, nz ! Normal vector
real(p2) :: mx, my, mz ! Orthogonal tangent vector
real(p2) :: lx, ly, lz ! Another orthogonal tangent vector
real(p2) :: abs_n_cross_l ! Magnitude of n x l
real(p2) :: uL, uR, vL, vR, wL, wR ! Velocity components.
real(p2) :: rhoL, rhoR, pL, pR ! Primitive variables.
real(p2) :: qnL, qnR, qmL, qmR, qlL, qlR ! Normal and tangent velocities
real(p2) :: aL, aR, HL, HR ! Speed of sound, Total enthalpy
real(p2) :: RT,rho,u,v,w,H,a,qn, ql, qm ! Roe-averages
real(p2) :: drho,dqn,dql,dqm,dp,LdU(5) ! Wave strengths
real(p2) :: ws(5), R(5,5) ! Wave speeds and right-eigenvectors
real(p2) :: dws(5) ! Width of a parabolic fit for entropy fix
real(p2) :: fL(5), fR(5), diss(5) ! Fluxes ad dissipation term
real(p2) :: gamma = 1.4_p2 ! Ratio of specific heats
real(p2) :: temp, tempx, tempy, tempz ! Temoprary variables
! Face normal vector (unit vector)
nx = njk(1)
ny = njk(2)
nz = njk(3)
tempx = ny*ny + nz*nz
tempy = nz*nz + nx*nx
tempz = nx*nx + ny*ny
if ( tempx >= tempy .and. tempx >= tempz ) then
lx = zero
ly = -nz
lz = ny
elseif ( tempy >= tempx .and. tempy >= tempz ) then
lx = -nz
ly = zero
lz = nx
elseif ( tempz >= tempx .and. tempz >= tempy ) then
lx = -ny
ly = nx
lz = zero
else
! Impossible to happen
write(*,*) "subroutine inviscid_roe: Impossible to happen. Please report the problem."
stop
endif
! Make it the unit vector.
temp = sqrt( lx*lx + ly*ly + lz*lz )
lx = lx/temp
ly = ly/temp
lz = lz/temp
mx = ny*lz - nz*ly
my = nz*lx - nx*lz
mz = nx*ly - ny*lx
abs_n_cross_l = sqrt(mx**2 + my**2 + mz**2)
mx = mx / abs_n_cross_l
my = my / abs_n_cross_l
mz = mz / abs_n_cross_l
rhoL = primL(1)
uL = primL(2)
vL = primL(3)
wL = primL(4)
qnL = uL*nx + vL*ny + wL*nz
qlL = uL*lx + vL*ly + wL*lz
qmL = uL*mx + vL*my + wL*mz
pL = primL(5)
aL = sqrt(gamma*pL/rhoL)
HL = aL*aL/(gamma-one) + half*(uL*uL+vL*vL+wL*wL)
! Right state
rhoR = primR(1)
uR = primR(2)
vR = primR(3)
wR = primR(4)
qnR = uR*nx + vR*ny + wR*nz
qlR = uR*lx + vR*ly + wR*lz
qmR = uR*mx + vR*my + wR*mz
pR = primR(5)
aR = sqrt(gamma*pR/rhoR)
HR = aR*aR/(gamma-one) + half*(uR*uR+vR*vR+wR*wR)
RT = sqrt(rhoR/rhoL)
rho = RT*rhoL !Roe-averaged density
u = (uL + RT*uR)/(one + RT) !Roe-averaged x-velocity
v = (vL + RT*vR)/(one + RT) !Roe-averaged y-velocity
w = (wL + RT*wR)/(one + RT) !Roe-averaged z-velocity
H = (HL + RT*HR)/(one + RT) !Roe-averaged total enthalpy
a = sqrt( (gamma-one)*(H-half*(u*u + v*v + w*w)) ) !Roe-averaged speed of sound
qn = u*nx + v*ny + w*nz !Roe-averaged face-normal velocity
ql = u*lx + v*ly + w*lz !Roe-averaged face-tangent velocity
qm = u*mx + v*my + w*mz !Roe-averaged face-tangent velocity
!Wave Strengths
drho = rhoR - rhoL !Density difference
dp = pR - pL !Pressure difference
dqn = qnR - qnL !Normal velocity difference
dql = qlR - qlL !Tangent velocity difference in l
dqm = qmR - qmL !Tangent velocity difference in m
LdU(1) = (dp - rho*a*dqn )/(two*a*a) !Left-moving acoustic wave strength
LdU(2) = drho - dp/(a*a) !Entropy wave strength
LdU(3) = (dp + rho*a*dqn )/(two*a*a) !Right-moving acoustic wave strength
LdU(4) = rho*dql !Shear wave strength
LdU(5) = rho*dqm !Shear wave strength
!Absolute values of the wave speeds
ws(1) = abs(qn-a) !Left-moving acoustic wave speed
ws(2) = abs(qn) !Entropy wave speed
ws(3) = abs(qn+a) !Right-moving acoustic wave speed
ws(4) = abs(qn) !Shear wave speed
ws(5) = abs(qn) !Shear wave speed
!Harten's Entropy Fix JCP(1983), 49, pp357-393: only for the nonlinear fields.
!NOTE: It avoids vanishing wave speeds by making a parabolic fit near ws = 0.
dws(1) = fifth
if ( ws(1) < dws(1) ) ws(1) = half * ( ws(1)*ws(1)/dws(1)+dws(1) )
dws(3) = fifth
if ( ws(3) < dws(3) ) ws(3) = half * ( ws(3)*ws(3)/dws(3)+dws(3) )
!Right Eigenvectors
! Left-moving acoustic wave
R(1,1) = one
R(2,1) = u - a*nx
R(3,1) = v - a*ny
R(4,1) = w - a*nz
R(5,1) = H - a*qn
! Entropy wave
R(1,2) = one
R(2,2) = u
R(3,2) = v
R(4,2) = w
R(5,2) = half*(u*u + v*v + w*w)
! Right-moving acoustic wave
R(1,3) = one
R(2,3) = u + a*nx
R(3,3) = v + a*ny
R(4,3) = w + a*nz
R(5,3) = H + a*qn
! Shear wave
R(1,4) = zero
R(2,4) = lx
R(3,4) = ly
R(4,4) = lz
R(5,4) = ql
! Shear wave
R(1,5) = zero
R(2,5) = mx
R(3,5) = my
R(4,5) = mz
R(5,5) = qm
diss(:) = ws(1)*LdU(1)*R(:,1) + ws(2)*LdU(2)*R(:,2) + ws(3)*LdU(3)*R(:,3) &
+ ws(4)*LdU(4)*R(:,4) + ws(5)*LdU(5)*R(:,5)
fL(1) = rhoL*qnL
fL(2) = rhoL*qnL * uL + pL*nx
fL(3) = rhoL*qnL * vL + pL*ny
fL(4) = rhoL*qnL * wL + pL*nz
fL(5) = rhoL*qnL * HL
fR(1) = rhoR*qnR
fR(2) = rhoR*qnR * uR + pR*nx
fR(3) = rhoR*qnR * vR + pR*ny
fR(4) = rhoR*qnR * wR + pR*nz
fR(5) = rhoR*qnR * HR
num_flux = half * (fL + fR - diss)
subroutine inviscid_roe_n(primL, primR, njk, num_flux)
implicit none
integer , parameter :: p2 = selected_real_kind(15) ! Double precision
!Input
real(p2), intent( in) :: primL(5), primR(5) ! Input: primitive variables
real(p2), intent( in) :: njk(3) ! Input: face normal vector
!Output
real(p2), intent(out) :: num_flux(5) ! Output: numerical flux
!Some constants
real(p2) :: zero = 0.0_p2
real(p2) :: one = 1.0_p2
real(p2) :: two = 2.0_p2
real(p2) :: half = 0.5_p2
real(p2) :: fifth = 0.2_p2
!Local variables
real(p2) :: nx, ny, nz ! Normal vector
real(p2) :: uL, uR, vL, vR, wL, wR ! Velocity components.
real(p2) :: rhoL, rhoR, pL, pR ! Primitive variables.
real(p2) :: qnL, qnR ! Normal velocities
real(p2) :: aL, aR, HL, HR ! Speed of sound, Total enthalpy
real(p2) :: RT,rho,u,v,w,H,a,qn ! Roe-averages
real(p2) :: drho,dqn,dp,LdU(4) ! Wave strengths
real(p2) :: du, dv, dw ! Velocity differences
real(p2) :: ws(4), R(5,4) ! Wave speeds and right-eigenvectors
real(p2) :: dws(4) ! Width of a parabolic fit for entropy fix
real(p2) :: fL(5), fR(5), diss(5) ! Fluxes ad dissipation term
real(p2) :: gamma = 1.4_p2 ! Ratio of specific heats
! Face normal vector (unit vector)
nx = njk(1)
ny = njk(2)
nz = njk(3)
!Primitive and other variables.
! Left state
rhoL = primL(1)
uL = primL(2)
vL = primL(3)
wL = primL(4)
qnL = uL*nx + vL*ny + wL*nz
pL = primL(5)
aL = sqrt(gamma*pL/rhoL)
HL = aL*aL/(gamma-one) + half*(uL*uL+vL*vL+wL*wL)
! Right state
rhoR = primR(1)
uR = primR(2)
vR = primR(3)
wR = primR(4)
qnR = uR*nx + vR*ny + wR*nz
pR = primR(5)
aR = sqrt(gamma*pR/rhoR)
HR = aR*aR/(gamma-one) + half*(uR*uR+vR*vR+wR*wR)
!First compute the Roe-averaged quantities
! NOTE: See http://www.cfdnotes.com/cfdnotes_roe_averaged_density.html for
! the Roe-averaged density.
RT = sqrt(rhoR/rhoL)
rho = RT*rhoL !Roe-averaged density
u = (uL + RT*uR)/(one + RT) !Roe-averaged x-velocity
v = (vL + RT*vR)/(one + RT) !Roe-averaged y-velocity
w = (wL + RT*wR)/(one + RT) !Roe-averaged z-velocity
H = (HL + RT*HR)/(one + RT) !Roe-averaged total enthalpy
a = sqrt( (gamma-one)*(H-half*(u*u + v*v + w*w)) ) !Roe-averaged speed of sound
qn = u*nx + v*ny + w*nz !Roe-averaged face-normal velocity
!Wave Strengths
drho = rhoR - rhoL !Density difference
dp = pR - pL !Pressure difference
dqn = qnR - qnL !Normal velocity difference
LdU(1) = (dp - rho*a*dqn )/(two*a*a) !Left-moving acoustic wave strength
LdU(2) = drho - dp/(a*a) !Entropy wave strength
LdU(3) = (dp + rho*a*dqn )/(two*a*a) !Right-moving acoustic wave strength
LdU(4) = rho !Shear wave strength (not really, just a factor)
!Absolute values of the wave Speeds
ws(1) = abs(qn-a) !Left-moving acoustic wave
ws(2) = abs(qn) !Entropy wave
ws(3) = abs(qn+a) !Right-moving acoustic wave
ws(4) = abs(qn) !Shear waves
!Harten's Entropy Fix JCP(1983), 49, pp357-393: only for the nonlinear fields.
!NOTE: It avoids vanishing wave speeds by making a parabolic fit near ws = 0.
dws(1) = fifth
if ( ws(1) < dws(1) ) ws(1) = half * ( ws(1)*ws(1)/dws(1)+dws(1) )
dws(3) = fifth
if ( ws(3) < dws(3) ) ws(3) = half * ( ws(3)*ws(3)/dws(3)+dws(3) )
R(1,1) = one
R(2,1) = u - a*nx
R(3,1) = v - a*ny
R(4,1) = w - a*nz
R(5,1) = H - a*qn
R(1,2) = one
R(2,2) = u
R(3,2) = v
R(4,2) = w
R(5,2) = half*(u*u + v*v + w*w)
! Right-moving acoustic wave
R(1,3) = one
R(2,3) = u + a*nx
R(3,3) = v + a*ny
R(4,3) = w + a*nz
R(5,3) = H + a*qn
! Two shear wave components combined into one (wave strength incorporated).
du = uR - uL
dv = vR - vL
dw = wR - wL
R(1,4) = zero
R(2,4) = du - dqn*nx
R(3,4) = dv - dqn*ny
R(4,4) = dw - dqn*nz
R(5,4) = u*du + v*dv + w*dw - qn*dqn
!Dissipation Term: |An|(UR-UL) = R|Lambda|L*dU = sum_k of [ ws(k) * R(:,k) * L*dU(k) ]
diss(:) = ws(1)*LdU(1)*R(:,1) + ws(2)*LdU(2)*R(:,2) &
+ ws(3)*LdU(3)*R(:,3) + ws(4)*LdU(4)*R(:,4)
!Compute the physical flux: fL = Fn(UL) and fR = Fn(UR)
fL(1) = rhoL*qnL
fL(2) = rhoL*qnL * uL + pL*nx
fL(3) = rhoL*qnL * vL + pL*ny
fL(4) = rhoL*qnL * wL + pL*nz
fL(5) = rhoL*qnL * HL
fR(1) = rhoR*qnR
fR(2) = rhoR*qnR * uR + pR*nx
fR(3) = rhoR*qnR * vR + pR*ny
fR(4) = rhoR*qnR * wR + pR*nz
fR(5) = rhoR*qnR * HR
! This is the numerical flux: Roe flux = 1/2 *[ Fn(UL)+Fn(UR) - |An|(UR-UL) ]
num_flux = half * (fL + fR - diss)
!Normal max wave speed in the normal direction.
! wsn = abs(qn) + a
end subroutine inviscid_roe_n
subroutine inviscid_rotated_rhll(primL, primR, njk, num_flux)
implicit none
integer , parameter :: p2 = selected_real_kind(15) ! Double precision
!Input
real(p2), intent( in) :: primL(5), primR(5) ! Input: primitive variables
real(p2), intent( in) :: njk(3) ! Input: face normal vector
!Output
real(p2), intent(out) :: num_flux(5) ! Output: numerical flux
!Some constants
real(p2) :: zero = 0.0_p2
real(p2) :: one = 1.0_p2
real(p2) :: two = 2.0_p2
real(p2) :: half = 0.5_p2
real(p2) :: fifth = 0.2_p2
!Local variables
real(p2) :: nx, ny, nz ! Face normal vector
real(p2) :: uL, uR, vL, vR, wL, wR ! Velocity components.
real(p2) :: rhoL, rhoR, pL, pR ! Primitive variables.
real(p2) :: qnL, qnR ! Normal velocity
real(p2) :: aL, aR, HL, HR ! Speed of sound, Total enthalpy
real(p2) :: RT,rho,u,v,w,H,a,qn ! Roe-averages
real(p2) :: drho,dqn,dp,LdU(4) ! Wave strengths
real(p2) :: du, dv, dw ! Velocity conponent differences
real(p2) :: eig(4) ! Eigenvalues
real(p2) :: ws(4), R(5,4) ! Absolute Wave speeds and right-eigenvectors
real(p2) :: dws(4) ! Width of a parabolic fit for entropy fix
real(p2) :: fL(5), fR(5), diss(5) ! Fluxes ad dissipation term
real(p2) :: gamma = 1.4_p2 ! Ratio of specific heats
real(p2) :: SRp,SLm ! Wave speeds for the HLL part
real(p2) :: nx1, ny1, nz1 ! Vector along which HLL is applied
real(p2) :: nx2, ny2, nz2 ! Vector along which Roe is applied
real(p2) :: alpha1, alpha2 ! Projections of the new normals
real(p2) :: abs_dq ! Magnitude of the velocity difference
real(p2) :: temp, tempx, tempy, tempz ! Temporary variables
! Face normal vector (unit vector)
nx = njk(1)
ny = njk(2)
nz = njk(3)
!Primitive and other variables.
! Left state
rhoL = primL(1)
uL = primL(2)
vL = primL(3)
wL = primL(4)
qnL = uL*nx + vL*ny + wL*nz
pL = primL(5)
aL = sqrt(gamma*pL/rhoL)
HL = aL*aL/(gamma-one) + half*(uL*uL+vL*vL+wL*wL)
! Right state
rhoR = primR(1)
uR = primR(2)
vR = primR(3)
wR = primR(4)
qnR = uR*nx + vR*ny + wR*nz
pR = primR(5)
aR = sqrt(gamma*pR/rhoR)
HR = aR*aR/(gamma-one) + half*(uR*uR+vR*vR+wR*wR)
!Compute the physical flux: fL = Fn(UL) and fR = Fn(UR)
fL(1) = rhoL*qnL
fL(2) = rhoL*qnL * uL + pL*nx
fL(3) = rhoL*qnL * vL + pL*ny
fL(4) = rhoL*qnL * wL + pL*nz
fL(5) = rhoL*qnL * HL
fR(1) = rhoR*qnR
fR(2) = rhoR*qnR * uR + pR*nx
fR(3) = rhoR*qnR * vR + pR*ny
fR(4) = rhoR*qnR * wR + pR*nz
fR(5) = rhoR*qnR * HR
abs_dq = sqrt( (uR-uL)**2 + (vR-vL)**2 + (wR-wL)**2 )
if ( abs_dq > 1.0e-12_p2) then
nx1 = (uR-uL)/abs_dq
ny1 = (vR-vL)/abs_dq
nz1 = (wR-wL)/abs_dq
tempx = ny*ny + nz*nz
tempy = nz*nz + nx*nx
tempz = nx*nx + ny*ny
if ( tempx >= tempy .and. tempx >= tempz ) then
nx1 = zero
ny1 = -nz
nz1 = ny
elseif ( tempy >= tempx .and. tempy >= tempz ) then
nx1 = -nz
ny1 = zero
nz1 = nx
elseif ( tempz >= tempx .and. tempz >= tempy ) then
nx1 = -ny
ny1 = nx
nz1 = zero
else
! Impossible to happen
write(*,*) "inviscid_rotated_rhll: Impossible to happen. Please report the problem."
stop
endif
! Make it the unit vector.
temp = sqrt( nx1*nx1 + ny1*ny1 + nz1*nz1 )
nx1 = nx1/temp
ny1 = ny1/temp
nz1 = nz1/temp
endif
alpha1 = nx*nx1 + ny*ny1 + nz*nz1
! Make alpha1 always positive.
temp = sign(one,alpha1)
nx1 = temp * nx1
ny1 = temp * ny1
nz1 = temp * nz1
alpha1 = temp * alpha1
!n2 = direction perpendicular to n1.
! Note: There are infinitely many choices for this vector.
! The best choice may be discovered in future.
! Here, we employ the formula (4.4) in the paper:
! (nx2,ny2,nz2) = (n1xn)xn1 / |(n1xn)xn1| ('x' is the vector product.)
! (tempx,tempy,tempz) = n1xn
tempx = ny1*nz - nz1*ny
tempy = nz1*nx - nx1*nz
tempz = nx1*ny - ny1*nx
! (nx2,ny2,nz2) = (n1xn)xn1
nx2 = tempy*nz1 - tempz*ny1
ny2 = tempz*nx1 - tempx*nz1
nz2 = tempx*ny1 - tempy*nx1
! Make n2 the unit vector
temp = sqrt( nx2*nx2 + ny2*ny2 + nz2*nz2 )
nx2 = nx2/temp
ny2 = ny2/temp
nz2 = nz2/temp
alpha2 = nx*nx2 + ny*ny2 + nz*nz2
! Make alpha2 always positive.
temp = sign(one,alpha2)
nx2 = temp * nx2
ny2 = temp * ny2
nz2 = temp * nz2
alpha2 = temp * alpha2
!--------------------------------------------------------------------------------
!Now we are going to compute the Roe flux with n2 as the normal with modified
!wave speeds (5.12). NOTE: the Roe flux here is computed without tangent vectors.
!See "I do like CFD, VOL.1" for details: page 57, Equation (3.6.31).
!First compute the Roe-averaged quantities
! NOTE: See http://www.cfdnotes.com/cfdnotes_roe_averaged_density.html for
! the Roe-averaged density.
RT = sqrt(rhoR/rhoL)
rho = RT*rhoL !Roe-averaged density.
u = (uL + RT*uR)/(one + RT) !Roe-averaged x-velocity
v = (vL + RT*vR)/(one + RT) !Roe-averaged y-velocity
w = (wL + RT*wR)/(one + RT) !Roe-averaged z-velocity
H = (HL + RT*HR)/(one + RT) !Roe-averaged total enthalpy
a = sqrt( (gamma-one)*(H-half*(u*u + v*v + w*w)) ) !Roe-averaged speed of sound
!----------------------------------------------------
!Compute the wave speed estimates for the HLL part,
!following Einfeldt:
!
! B. Einfeldt, On Godunov-type methods for gas dynamics,
! SIAM Journal on Numerical Analysis 25 (2) (1988) 294–318.
!
! Note: HLL is actually applied to n1, but this is
! all we need to incorporate HLL. See JCP2008 paper.
qn = u *nx1 + v *ny1 + w *nz1
qnL = uL*nx1 + vL*ny1 + wL*nz1
qnR = uR*nx1 + vR*ny1 + wR*nz1
SLm = min( zero, qn - a, qnL - aL ) !Minimum wave speed estimate
SRp = max( zero, qn + a, qnR + aR ) !Maximum wave speed estimate
! This is the only place where n1=(nx1,ny1,nz1) is used.
! n1=(nx1,ny1,nz1) is never used below.
!----------------------------------------------------
!Wave Strengths
qn = u *nx2 + v *ny2 + w *nz2
qnL = uL*nx2 + vL*ny2 + wL*nz2
qnR = uR*nx2 + vR*ny2 + wR*nz2
drho = rhoR - rhoL !Density difference
dp = pR - pL !Pressure difference
dqn = qnR - qnL !Normal velocity difference
LdU(1) = (dp - rho*a*dqn )/(two*a*a) !Left-moving acoustic wave strength
LdU(2) = drho - dp/(a*a) !Entropy wave strength
LdU(3) = (dp + rho*a*dqn )/(two*a*a) !Right-moving acoustic wave strength
LdU(4) = rho !Shear wave strength (not really, just a factor)
!Wave Speed (Eigenvalues)
eig(1) = qn-a !Left-moving acoustic wave velocity
eig(2) = qn !Entropy wave velocity
eig(3) = qn+a !Right-moving acoustic wave velocity
eig(4) = qn !Shear wave velocity
!Absolute values of the wave speeds (Eigenvalues)
ws(1) = abs(qn-a) !Left-moving acoustic wave speed
ws(2) = abs(qn) !Entropy wave speed
ws(3) = abs(qn+a) !Right-moving acoustic wave speed
ws(4) = abs(qn) !Shear wave speed
!Harten's Entropy Fix JCP(1983), 49, pp357-393: only for the nonlinear fields.
!NOTE: It avoids vanishing wave speeds by making a parabolic fit near ws = 0.
dws(1) = fifth
if ( ws(1) < dws(1) ) ws(1) = half * ( ws(1)*ws(1)/dws(1)+dws(1) )
dws(3) = fifth
if ( ws(3) < dws(3) ) ws(3) = half * ( ws(3)*ws(3)/dws(3)+dws(3) )
!Combine the wave speeds for Rotated-RHLL: Eq.(5.12) in the original JCP2008 paper.
ws = alpha2*ws - (alpha1*two*SRp*SLm + alpha2*(SRp+SLm)*eig)/(SRp-SLm)
!Below, we compute the Roe dissipation term in the direction n2
!with the above modified wave speeds. HLL wave speeds act something like
!the entropy fix or eigenvalue limiting; they contribute only by the amount
!given by the fraction, alpha1 (less than or equal to 1.0). See JCP2008 paper.
!Right Eigenvectors:
!Note: Two shear wave components are combined into one, so that tangent vectors
! are not required. And that's why there are only 4 vectors here.
! Left-moving acoustic wave
R(1,1) = one
R(2,1) = u - a*nx2
R(3,1) = v - a*ny2
R(4,1) = w - a*nz2
R(5,1) = H - a*qn
! Entropy wave
R(1,2) = one
R(2,2) = u
R(3,2) = v
R(4,2) = w
R(5,2) = half*(u*u + v*v + w*w)
! Right-moving acoustic wave
R(1,3) = one
R(2,3) = u + a*nx2
R(3,3) = v + a*ny2
R(4,3) = w + a*nz2
R(5,3) = H + a*qn
! Two shear wave components combined into one (wave strength incorporated).
du = uR - uL
dv = vR - vL
dw = wR - wL
R(1,4) = zero
R(2,4) = du - dqn*nx2
R(3,4) = dv - dqn*ny2
R(4,4) = dw - dqn*nz2
R(5,4) = u*du + v*dv + w*dw - qn*dqn
!Dissipation Term: Roe dissipation with the modified wave speeds.
! |An|dU = R|Lambda|L*dU = sum_k of [ ws(k) * R(:,k) * L*dU(k) ], where n=n2.
diss(:) = ws(1)*LdU(1)*R(:,1) + ws(2)*LdU(2)*R(:,2) &
+ ws(3)*LdU(3)*R(:,3) + ws(4)*LdU(4)*R(:,4)
!Compute the Rotated-RHLL flux. (It looks like the HLL flux with Roe dissipation.)
num_flux = (SRp*fL - SLm*fR)/(SRp-SLm) - half*diss
!Normal max wave speed in the normal direction.
! wsn = abs(qn) + a
end subroutine inviscid_rotated_rhll
!--------------------------------------------------------------------------------
解决方案
Your file is not a program at all! It is a collection of subprograms. You cannot compile it for running as a program, only as an object file or a library (try -c or -shared). You must add the main program body to be able to compile it as a program and run it!
这篇关于在cygwin / gfortran下编译一个Fortran程序给出了“_WinMain @ 16”的未定义引用“的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
查看全文
