您的位置：360文档中心› Simple算法_fortran语言编写

Simple算法_fortran语言编写

合集下载

fortran 语言编程

fortran 语言编程

fortran 语言编程Fortran 语言编程Fortran（Formula Translation）是一种面向科学计算和工程计算的高级编程语言。

它于1957年诞生于IBM，是最早被广泛采用的科学计算语言之一，目前已经发展到第四个版本（Fortran 2018）。

Fortran是一种编译型语言，它通过编写源代码并使用编译器将其转换成机器语言来执行。

本文将详细介绍Fortran语言的基础知识、语法规则和常用的编程技巧，以帮助读者了解和掌握这门强大的科学计算语言。

第一步：安装Fortran编译器要开始编写和运行Fortran程序，首先需要安装Fortran编译器。

有多种Fortran编译器可供选择，其中最常用的是GNU Fortran（gfortran）和Intel Fortran Compiler（ifort）。

可以从官方网站或其他可信的来源获得这些编译器的安装程序，并按照提示进行安装。

第二步：编写并编译Fortran程序在开始编写Fortran程序之前，需要了解Fortran的基本语法规则。

Fortran使用固定格式或自由格式，固定格式的源代码按照列格式排列，每行的前6列被保留用于行号和注释，从第7列开始是可执行代码。

自由格式没有列格式的限制，更加灵活，但在编译阶段需要指定自由格式。

下面是一个简单的Fortran程序示例，用于计算并输出两个数的和：fortranprogram additionimplicit noneinteger :: a, b, sumprint *, "Enter two numbers:"read *, a, bsum = a + bprint *, "The sum is:", sumend program addition将以上代码保存为一个以.f90为后缀名的文件（例如addition.f90），然后使用编译器将其编译成可执行程序。

Fortran科学计算编程课件

Fortran科学计算编程课件

Fortran科学计算编程课件Fortran（即Formula Translation的缩写）是一种用于科学计算和数值分析的编程语言。

由于其高效性和可靠性，Fortran在科学计算领域广泛应用。

本课件将为您介绍Fortran编程语言，并提供相关的科学计算示例和实践。

第一部分：Fortran简介Fortran是一种面向过程的编程语言，由IBM在1950年代开发，旨在为科学计算提供高性能和可移植性。

Fortran以其强大的数学计算和数组处理能力闻名，是科学家和工程师们首选的编程语言之一。

第二部分：Fortran基本语法2.1 变量和数据类型在Fortran中，变量的类型必须在声明时指定，并且不同类型的变量可以存储不同类型的数据，如整数、实数和字符等。

这一部分将详细介绍Fortran的常用数据类型和变量声明方法。

2.2 控制结构Fortran支持常见的控制结构，例如条件控制语句（if-else语句）、循环语句（do循环）和跳转语句（goto语句）。

我们将逐一介绍这些结构的语法和使用方法，并给出实例演示。

第三部分：科学计算实例为了帮助您更好地理解Fortran的应用，本课件将提供一些科学计算的实例，包括线性代数运算、梯度下降算法等。

通过这些实例，您将学会如何用Fortran编写科学计算程序，并了解其在实践中的应用。

第四部分：Fortran编程工具和资源4.1 编译器Fortran编程需要一个符合语言标准的编译器来将代码转化为可执行文件。

我们将向您介绍几个常用的Fortran编译器，并给出安装和配置的指南。

4.2 开发环境为了提高编程效率，许多IDE（集成开发环境）和编辑器提供了对Fortran的支持。

我们将为您介绍几个常用的Fortran开发环境，并给出使用方法和建议。

4.3 学习资源如果您想深入学习Fortran，我们还准备了一些在线教程和学习资源的链接，方便您进一步扩展知识和技能。

结语通过本课件的学习，您将对Fortran科学计算编程有更深入的了解，并掌握基本的语法和应用技巧。

fortran编程语言概述

fortran编程语言概述

fortran编程语言概述Fortran编程语言概述【引言】Fortran（Formula Translation，公式翻译）是一种高级程序设计语言，最初为科学和工程计算而设计。

自20世纪50年代问世以来，Fortran已成为一种广泛使用的编程语言，为数值计算和科学计算提供了强大的功能和效率。

本文将深入探讨Fortran编程语言的多个方面，包括历史背景、语法特性、应用领域和优缺点等。

【历史背景】Fortran的发展可以追溯到20世纪50年代，当时计算机主要用于数值计算和科学研究。

最早的Fortran编译器于1957年发布，它为科学家和工程师提供了一种简单而强大的编程工具。

随着计算机技术的进步和计算需求的不断增长，Fortran逐渐演化为各种版本，包括Fortran 77、Fortran 90、Fortran 95、Fortran 2003和Fortran 2018等。

【语法特性】Fortran具有一套独特的语法规则和特性，使其适用于科学和工程计算。

首先，Fortran使用固定格式的源代码布局，即每行代码的前六个字符用于行标签，接着是区分符号，然后是代码语句。

其次，Fortran允许使用自然语言类似的语法，使得代码易于阅读和理解。

此外，Fortran 提供了丰富的数学和科学函数，以及数组和矩阵操作，用于处理复杂的数值计算。

【应用领域】由于其优秀的数值计算性能和丰富的科学计算库，Fortran被广泛应用于各个领域。

首先，天文学和物理学领域利用Fortran进行天体力学模拟和物理过程建模。

其次，气象学和气候研究使用Fortran编写气候模型和天气预测程序。

此外，工程学、生物学、化学、经济学等领域也从Fortran的数值计算能力中受益匪浅。

【优点】Fortran具有以下几个明显的优点。

首先，Fortran在数值计算和科学计算方面具有出色的性能，其编译器能够针对各种硬件架构进行优化。

其次，Fortran带有一套丰富的科学计算库，提供了许多常用的数学和科学函数，方便开发人员进行复杂的计算。

fortran编程的步骤

fortran编程的步骤

fortran编程的步骤Fortran编程的步骤一、引言Fortran（Formula Translation）是一种高级程序设计语言，特别适用于科学计算和数值计算。

本文将介绍Fortran编程的步骤，帮助初学者了解如何使用Fortran进行程序开发。

二、编写程序的基本步骤1. 确定程序的目标：在开始编写Fortran程序之前，需要明确程序的目标和需求。

确定程序的输入和输出，以及所需的计算或处理步骤。

这有助于确保编写的程序能够满足预期的功能和要求。

2. 设计算法和数据结构：根据程序的目标，设计合适的算法和数据结构。

算法描述了解决问题的步骤和逻辑，而数据结构则定义了程序中使用的数据类型和数据组织方式。

3. 编写代码：根据算法和数据结构的设计，开始编写Fortran代码。

Fortran使用特定的语法和语句结构，需要熟悉其语法规则和常用的编程技巧。

代码的编写应遵循良好的编码风格，包括适当的缩进、注释和命名规范。

4. 编译程序：编写完Fortran代码后，需要使用Fortran编译器将源代码转换成可执行的机器代码。

编译过程将检查代码中的语法错误和逻辑错误，并生成可执行文件。

Fortran编译器通常会提供丰富的编译选项，可以根据需要进行调整。

5. 调试和测试：编译成功后，可以对程序进行调试和测试。

调试是指查找和修复程序中的错误和问题，测试是指验证程序的正确性和性能。

调试和测试是编程过程中不可或缺的环节，可以使用调试器和测试框架等工具辅助进行。

6. 优化和性能调整：在程序调试和测试完成后，可以考虑对程序进行优化和性能调整。

优化旨在改进程序的执行效率和资源利用率，可以通过改进算法、调整编译选项和使用高级优化技术来实现。

7. 文档撰写：在编程过程中，应及时记录程序的设计和实现细节。

文档可以包括程序的功能描述、算法和数据结构的说明、代码注释和使用说明等。

良好的文档可以提高代码的可读性和可维护性，并方便其他人理解和使用程序。

FORTRAN 程序设计01

FORTRAN 程序设计01
FORTRAN 程序设计
一, Fortran 简介 二,程序流程 三,循环结构 四,数据类型 五,数组 六,过程和模块 七,输入,输出,文件 八,复习,总结
Fortran是目前国际上广泛流行的一种高级语言, 适用于科学计算.Fortran是英文FORmula TRANslation的缩写,意为"公式翻译".它 是为科学,工程问题中的那些能够用数学公式 表达的问题而设计的语言,主要用于数值计算. 这种语言简单易学,因为可以像抄写数学教科 书里的公式一样书写数学公式,它比英文书写 的自然语言更接近数学语言.Fortran语言是 第一个真正推广的高级语言.
Fortran77 (I ～ N)规则; 以(I,j,k,m,l,n)字母开头的 变量为整型. 建议编程时使用声明语句,并在 程序中加入implicit non 语句, 屏蔽(I ～ N)规则;
asb exp sin cos Asin Acos Sqrt ……
asb(x) |x| 指数运算 exp(x) exp(x) 正弦函数 sin(x) sinx 余弦函数 cos(x) cosx 反正弦 Asin(x) arcsinx 反余弦 Acos(x) arccosx 开平方 Sqrt(x) x ………………
系统为每一个常量,变量分配一个存储单元,放它的 值.
6, 数据类型
整型:(数学上的整数集合) integer(n),属性列表,变量列表 例:integer (kind=4) i,t2,pop -2147483638 ～ 2147483637 例:integer (1):: k,m -128 ～127
实型(数学上的实数集合) real(n),属性列表,变量列表 单精度: real(4), real 双精度: real(8), double real(4),:: dx,dy real(8),:: d_p,y2

SIMPLE算法编程实例

SIMPLE算法编程实例

c Demonstration of the SIMPLE (overlapping grid) scheme for a flow problem in a duct cc Dec. 2007 Hangzhou Dianzi Universitydimension u(32),uf(33),a(33),p(32),pp(32),aa(32),su(32)dimension aup(33),aue(33),auw(33),ap(32),ae(32),aw(32),s(32)c u --- velocity on main cell centerc uf --- velocity on main cell facec a --- cell side areac p --- pressurec pp --- pressure correctionc aa --- cell center cross-section areac su --- source term for momentum eqationc aup --- main coefficient for momentum equationc aue --- east coefficient for momentum equationc auw --- west coefficient for momentum equationc ap --- main coefficient for pressure correction equationc ae --- east coefficient for pressure correction equationc aw --- west coefficient for pressure correction equationc s --- source term for pressure correction equationc grid-generation for the converging-diverging ductdx=0.01do i=1,33aaa=0.1/16.**2bb=0.1a(i)=aaa*float(i-17)**2+bbend dodo i=1,32aa(i)=.5*(a(i)+a(i+1))end doc parameter settingden=1000.c boundary conditionsu(1)=10.u(32)=10.viscos=0.001c boundary condition for face velocityuf(1)=10.uf(33)=10.c initial conditiondo i=2,32u(i)=10.uf(i)=10.p(i)=0.end doc do-loop starts (iteration starts)do N=1,1000c solving momentum equation on 2nd griddo i=2,32aue(i)=amax1(abs(den*u(i)*AA(i)/2.),viscos*aa(i)/dx)* -den*u(i)*AA(i)/2.auw(i)=amax1(abs(den*u(i-1)*AA(i-1)/2.),viscos*aa(i-1)/dx)* +den*u(i-1)*AA(i-1)/2.su(i)=a(i)*(p(i-1)-p(i))aup(i)=auw(i)+aue(i)end doc write(6,*)aue(10),auw(10),aup(10),su(10),' aue,auw,aup,su'c small iterationsdo i=2,32uf(i)=(aue(i)*uf(i+1)+auw(i)*uf(i-1)+su(i))/aup(i)end doc extrapolationc aue(1)=aue(2)c auw(1)=auw(2)c aue(33)=aue(32)c auw(33)=auw(32)c solving pressure correction equation on main griddo i=2,31aw(i)=den*a(i)**2/aup(i)ae(i)=den*a(i+1)**2/aup(i+1)ap(i)=aw(i)+ae(i)s(i)=den*a(i)*uf(i)-den*a(i+1)*uf(i+1)end doc set pp=0do i=1,32pp(i)=0.end doc small iterationsdo its=1,30do i=2,31pp(i)=(ae(i)*pp(i+1)+aw(i)*pp(i-1)+s(i))/ap(i)end doend doc correct velocitydo i=2,32uf(i)=uf(i)+(pp(i-1)-pp(i))*a(i)/aup(i)end doc interpolationdo i=1,32u(i)=.5*(uf(i+1)+uf(i))end doc correct pressuredo i=2,31p(i)=p(i)+0.1*pp(i)end doc extrapolationp(1)=2.*p(2)-p(3)p(32)=2.*p(31)-p(30)c monitoring the iterationswrite(6,*)'N=',N, (' u(',i,')=',u(i),i=10,20,10)write(6,*)'N=',N, (' uf(',i,')=',uf(i),i=10,20,10) write(6,*)'N=',N, (' p(',i,')=',p(i),i=10,20,10)end doc results writing into a fileopen(1,file='result.txt')do i=1,32write(1,*)'i=',i,'p(',i,')=',p(i),' u(',i,')=',u(i) * ,' uf(',i,')=',uf(i)end dostopend。

fortran 编程指南

fortran 编程指南

fortran 编程指南英文回答：Fortran is a programming language that was developed in the 1950s and is still widely used today, especially in scientific and engineering applications. It stands for Formula Translation and was designed to be used for numerical and scientific computing. Fortran has evolved over the years and the latest version is Fortran 2018.One of the main advantages of Fortran is its efficiency in handling numerical computations. It has built-in support for arrays and mathematical functions, making it ideal for scientific calculations. Fortran also has a strong static typing system, which helps catch errors at compile-time and improves performance.Another benefit of Fortran is its extensive library of mathematical and scientific functions. These libraries provide a wide range of pre-written code for tasks such aslinear algebra, numerical integration, and solving differential equations. This saves programmers time and effort, as they don't have to write these functions from scratch.Fortran also has a long history and a large communityof users, which means there are plenty of resourcesavailable for learning and troubleshooting. There are numerous books, tutorials, and online forums dedicated to Fortran programming. This makes it easier for beginners to get started and for experienced programmers to find help when needed.One of the drawbacks of Fortran is its syntax, whichcan be seen as outdated compared to modern programming languages. It uses a fixed-format style, where each linehas a specific structure and indentation is not significant. This can make the code look less readable and harder to maintain.Furthermore, Fortran is not as versatile as some other languages when it comes to non-numerical tasks. It lacksbuilt-in support for string manipulation and file I/O operations, which can be limiting in certain applications. However, there are ways to work around these limitations by using external libraries or interfacing with other languages.In conclusion, Fortran is a powerful language for numerical and scientific computing, with a long history and a large community of users. It offers efficiency, extensive libraries, and plenty of resources for learning and troubleshooting. While its syntax may be seen as outdated, it remains a popular choice for scientific and engineering applications.中文回答：Fortran是一种编程语言，于1950年代开发，并且至今仍广泛应用于科学和工程领域。

SIMPLE算法讲解

SIMPLE算法讲解

SIMPLE算法讲解在计算机科学中，算法是一组解决问题的明确指令步骤。

算法可以用于执行各种复杂任务，从图像处理到数据排序。

而SIMPLE（The Symbolic Instruction Manipulation Program for Libraries and Education）算法是为了简化低级机器指令的编写而产生的一种高级编程语言。

本篇文章将简要介绍SIMPLE算法及其应用。

SIMPLE算法最初由美国化学家Daniel D. McCracken于1960年代开发，旨在为非计算机专业人士提供一种易于学习和使用的编程语言。

SIMPLE算法的语法类似于英语，因此非常易于理解和阅读。

它被广泛应用于教育领域，并被用作初学者编程课程的教学工具。

```BEGININPUTa,bc:=a+bOUTPUTcEND```上述代码的功能是输入两个数值a和b，然后将它们相加，并将结果输出。

通过这个简单的例子，我们可以看到SIMPLE算法注重指令的可读性和可理解性。

除了基本的输入输出操作，SIMPLE算法还提供了条件语句和循环语句等控制结构，使得编写复杂程序也变得相对容易。

例如，我们可以使用条件语句来实现一个简单的判断奇偶数的程序：```BEGININPUT numIF num % 2 = 0 THENOUTPUT "Even"ELSEOUTPUT "Odd"ENDIFEND```上述代码通过使用IF-ELSE语句对给定的数字进行判断，并输出相应的结果。

通过使用条件语句，我们可以根据不同的条件执行不同的指令，从而实现程序的分支控制。

SIMPLE算法还提供了多种数据类型和操作符，包括整型、字符型、浮点型、数组等。

这使得编写更复杂的程序成为可能。

同时，SIMPLE算法不仅可以用于简单的算术和逻辑操作，还可以进行文件读写、排序等更高级的计算操作。

尽管SIMPLE算法已经有一段时间没有被广泛使用，但它的原则和设计理念仍然对现代编程有一定的影响。

  1. 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

CcccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccCccc This computer program was copied from the graduate student course programCccc of the University of Minnesota. Part of it was re-formulated to meet theCccc personal computer environment. Some inappropriate expressions were alsoCccc corrected. The program is used only for teaching purpose. No part of itCccc may be published. You may use it as a frame to re-develop your own codeCccc for research purpose. XJTU Instructor, 1995.11**************************************************************************** *----------------------------MAIN PROGRAM-----------------------------------**************************************************************************** LOGICAL LSTOPCOMMON/CNTL/LSTOP**************************************************************************** OPEN(08,FILE='teresul')CALL SETUP0CALL GRIDCALL SETUP1CALL START10 CALL DENSECALL BOUNDCALL OUTPUTIF(.NOT.LSTOP) GO TO 15CLOSE(08)STOP15 CALL SETUP2GO TO 10END*---------------------------------------------------------------------------SUBROUTINE DIFLOW**************************************************************************** COMMON/COEF/FLOW,DIFF,ACOF**************************************************************************** ACOF=DIFFIF(FLOW .EQ.0.0)RETURNTEMP=DIFF-ABS(FLOW)*0.1ACOF=0.IF(TEMP .LE. 0. ) RETURNTEMP=TEMP/DIFFACOF=DIFF*TEMP**5RETURNEND*--------------------------------------------------------------------------SUBROUTINE SOLVE****************************************************************************DOUBLE PRECISION TITLELOGICAL LSOLVE,LPRINT,LBLK,LSTOPCOMMON F(22,22,10),P(22,22),RHO(22,22),GAM(22,22),CON(22,22),& AIP(22,22),AIM(22,22),AJP(22,22),AJM(22,22),AP(22,22),& X(22),XU(22),XDIF(22),XCV(22),XCVS(22),& Y(22),YV(22),YDIF(22),YCV(22),YCVS(22),& YCVR(22),YCVRS(22),ARX(22),ARXJ(22),ARXJP(22),& R(22),RMN(22),SX(22),SXMN(22),XCVI(22),XCVIP(22)COMMON DU(22,22),DV(22,22),FV(22),FVP(22),& FX(22),FXM(22),FY(22),FYM(22),PT(22),QT(22)COMMON /INDX/NF,NFMAX,NP,NRHO,NGAM,L1,L2,L3,M1,M2,M3,& IST,JST,ITER,LAST,TITLE(13),RELAX(13),TIME,DT,XL,YL,&IPREF,JPREF,LSOLVE(10),LPRINT(13),LBLK(10),MODE,NTIMES(10),RHOCON**************************************************************************** ISTF=IST-1JSTF=JST-1IT1=L2+ISTIT2=L3+ISTJT1=M2+JSTJT2=M3+JST**************************************************************************** DO 999 NT=1,NTIMES(NF)DO 999 N=NF,NF*---------------------------------------------------------------------------IF(.NOT. LBLK(NF)) GO TO 10PT(ISTF)=0.QT(ISTF)=0.DO 11 I=IST,L2BL=0.BLP=0.BLM=0.BLC=0.DO 12 J=JST,M2BL=BL+AP(I,J)IF(J .NE. M2) BL=BL-AJP(I,J)IF(J .NE. JST) BL=BL-AJM(I,J)BLP=BLP+AIP(I,J)BLM=BLM+AIM(I,J)BLC=BLC+CON(I,J)+AIP(I,J)*F(I+1,J,N)+AIM(I,J)*F(I-1,J,N)& +AJP(I,J)*F(I,J+1,N)+AJM(I,J)*F(I,J-1,N)-AP(I,J)*F(I,J,N)12 CONTINUEDENOM=BL-PT(I-1)*BLMDENO=1.E15IF(ABS(DENOM/BL) .LT. 1.E-10) DENOM=1.E20*DENOQT(I)=(BLC+BLM*QT(I-1))/DENOM11 CONTINUEBL=0.DO 13 II=IST,L2I=IT1-IIBL=BL*PT(I)+QT(I)DO 13 J=JST,M213 F(I,J,N)=F(I,J,N)+BL*--------------------------------------------------------------------------- PT(JSTF)=0.QT(JSTF)=0.DO 21 J=JST,M2BL=0.BLP=0.BLM=0.BLC=0.DO 22 I=IST,L2BL=BL+AP(I,J)IF(I .NE. L2) BL=BL-AIP(I,J)IF(I .NE. IST) BL=BL-AIM(I,J)BLP=BLP+AJP(I,J)BLM=BLM+AJM(I,J)BLC=BLC+CON(I,J)+AIP(I,J)*F(I+1,J,N)+AIM(I,J)*F(I-1,J,N)& +AJP(I,J)*F(I,J+1,N)+AJM(I,J)*F(I,J-1,N)-AP(I,J)*F(I,J,N)22 CONTINUEDENOM=BL-PT(J-1)*BLMIF (ABS(DENOM/BL) .LT. 1E-10) DENOM=1.E20*DENOPT(J)=BLP/DENOMQT(J)=(BLC+BLM*QT(J-1))/DENOM21 CONTINUEBL=0.DO 23 JJ=JST,M2J=JT1-JJBL=BL*PT(J)+QT(J)DO 23 I=IST,L223 F(I,J,N)=F(I,J,N)+BL10 CONTINUE*-----------------------------------------------------------------------DO 90 J=JST,M2PT(ISTF)=0.QT(ISTF)=F(ISTF,J,N)DO 70 I=IST,L2DENOM=AP(I,J)-PT(I-1)*AIM(I,J)TEMP=CON(I,J)+AJP(I,J)*F(I,J+1,N)+AJM(I,J)*F(I,J-1,N)QT(I)=(TEMP+AIM(I,J)*QT(I-1))/DENOM70 CONTINUEDO 80 II=IST,L2I=IT1-II80 F(I,J,N)=F(I+1,J,N)*PT(I)+QT(I)90 CONTINUE*----------------------------------------------------------------------- DO 190 JJ=JST,M3J=JT2-JJPT(ISTF)=0.QT(ISTF)=F(ISTF,J,N)DO 170 I=IST,L2DENOM=AP(I,J)-PT(I-1)*AIM(I,J)PT(I)=AIP(I,J)/DENOMTEMP=CON(I,J)+AJP(I,J)*F(I,J+1,N)+AJM(I,J)*F(I,J-1,N)QT(I)=(TEMP+AIM(I,J)*QT(I-1))/DENOM170 CONTINUEDO 180 II=IST,L2I=IT1-II180 F(I,J,N)=F(I+1,J,N)*PT(I)+QT(I)190 CONTINUE*----------------------------------------------------------------------- DO 290 I=IST,L2PT(JSTF)=0.QT(JSTF)=F(I,JSTF,N)DO 270 J=JST,M2DENOM=AP(I,J)-PT(J-1)*AJM(I,J)PT(J)=AJP(I,J)/DENOMTEMP=CON(I,J)+AIP(I,J)*F(I+1,J,N)+AIM(I,J)*F(I-1,J,N)QT(J)=(TEMP+AJM(I,J)*QT(J-1))/DENOM270 CONTINUEDO 280 JJ=JST,M2J=JT1-JJ280 F(I,J,N)=F(I,J+1,N)*PT(J)+QT(J)290 CONTINUE*----------------------------------------------------------------------- DO 390 II=IST,L3I=IT2-IIPT(JSTF)=0.QT(JSTF)=F(I,JSTF,N)DO 370 J=JST,M2DENOM=AP(I,J)-PT(J-1)*AJM(I,J)TEMP=CON(I,J)+AIP(I,J)*F(I+1,J,N)+AIM(I,J)*F(I-1,J,N)QT(J)=(TEMP+AJM(I,J)*QT(J-1))/DENOM370 CONTINUEDO 380 JJ=JST,M2J=JT1-JJ380 F(I,J,N)=F(I,J+1,N)*PT(J)+QT(J)390 CONTINUE************************************************************************ 999 CONTINUEDO 400 J=2,M2DO 400 I=2,L2CON(I,J)=0.AP(I,J)=0.400 CONTINUERETURNEND************************************************************************ SUBROUTINE SETUP************************************************************************ DOUBLE PRECISION TITLELOGICAL LSOLVE,LPRINT,LBLK,LSTOPCOMMON F(22,22,10),P(22,22),RHO(22,22),GAM(22,22),CON(22,22),& AIP(22,22),AIM(22,22),AJP(22,22),AJM(22,22),AP(22,22),& X(22),XU(22),XDIF(22),XCV(22),XCVS(22),& Y(22),YV(22),YDIF(22),YCV(22),YCVS(22),&YCVR(22),YCVRS(22),ARX(22),ARXJ(22),ARXJP(22),&R(22),RMN(22),SX(22),SXMN(22),XCVI(22),XCVIP(22)COMMON DU(22,22),DV(22,22),FV(22),FVP(22),& FX(22),FXM(22),FY(22),FYM(22),PT(22),QT(22)COMMON /INDX/NF,NFMAX,NP,NRHO,NGAM,L1,L2,L3,M1,M2,M3,& IST,JST,ITER,LAST,TITLE(13),RELAX(13),TIME,DT,XL,YL,& IPREF,JPREF,LSOLVE(10),LPRINT(13),LBLK(10),MODE,NTIMES(10),RHOCONCOMMON/CNTL/LSTOPCOMMON/SORC/SMAX,SSUMCOMMON/COEF/FLOW,DIFF,ACOFDIMENSION U(22,22),V(22,22),PC(22,22)EQUIVALENCE (F(1,1,1),U(1,1)),(F(1,1,2),V(1,1)),(F(1,1,3),PC(1,1))************************************************************************1 FORMAT(//15X,'COMPUTATION IN CARTISIAN COORDINATES')2 FORMAT(//15X,'COMPUTATION FOR AXISYMMETRICAL SITUATION')3 FORMAT(//15X,' COMPUTATION IN POLAR COORDINATES ')4 FORMAT(1X,14X,40(1H*),//)*-----------------------------------------------------------------------ENTRY SETUP0NFMAX=10NP=11NRHO=12NGAM=13LSTOP=.FALSE.DO 779 I=1,10LSOLVE(I)=.FALSE.LBLK(I)=.TRUE.779 NTIMES(I)=1DO 889 I=1,13LPRINT(I)=.FALSE.889 RELAX(I)=1.MODE=1LAST=5TIME=0.ITER=0DT=1.0E+10IPREF=1JPREF=1RHOCON=1RETURN*----------------------------------------------------------------------- ENTRY SETUP1L2=L1-1L3=L2-1M2=M1-1M3=M2-1X(1)=XU(2)DO 5 I=2,L25 X(I)=0.5*(XU(I+1)+XU(I))X(L1)=XU(L1)Y(1)=YV(2)DO 10 J=2,M210 Y(J)=0.5*(YV(J+1)+YV(J))Y(M1)=YV(M1)DO 15 I=2,L115 XDIF(I)=X(I)-X(I-1)DO 18 I=2,L218 XCV(I)=XU(I+1)-XU(I)DO 20 I=3,L220 XCVS(I)=XDIF(I)XCVS(3)=XCVS(3)+XDIF(2)XCVS(L2)=XCVS(L2)+XDIF(L1)DO 22 I=3,L3XCVI(I)=0.5*XCV(I)22 XCVIP(I)=XCVI(I)XCVIP(2)=XCV(2)XCVI(L2)=XCV(L2)DO 35 J=2,M135 YDIF(J)=Y(J)-Y(J-1)DO 40 J=2,M240 YCV(J)=YV(J+1)-YV(J)DO 45 J=3,M245 YCVS(J)=YDIF(J)YCVS(3)=YCVS(3)+YDIF(2)YCVS(M2)=YCVS(M2)+YDIF(M1)IF (MODE .NE. 1) GO TO 55DO 52 J=1,M1RMN(J)=1.52 R(J)=1.GO TO 5655 DO 50 J=2,M150 R(J)=R(J-1)+YDIF(J)RMN(2)=R(1)DO 60 J=3,M260 RMN(J)=RMN(J-1)+YCV(J-1)RMN(M1)=R(M1)56 CONTINUEDO 57 J=1,M1SX(J)=1.SXMN(J)=1.IF(MODE .NE. 3) GO TO 57SX(J)=R(J)IF(J .NE. 1) SXMN(J)=RMN(J)57 CONTINUEDO 62 J=2,M2YCVR(J)=R(J)*YCV(J)ARX(J)=YCVR(J)IF (MODE .NE. 3) GO TO 62ARX(J)=YCV(J)62 CONTINUEDO 64 J=4,M364 YCVRS(J)=0.5*(R(J)+R(J-1))*YDIF(J)YCVRS(3)=0.5*(R(3)+R(1))*YCVS(3)YCVRS(M2)=0.5*(R(M1)+R(M3))*YCVS(M2) IF(MODE .NE. 2) GO TO 67DO 65 J=3,M3ARXJ(J)=0.25*(1.+RMN(J)/R(J))*ARX(J)65 ARXJP(J)=ARX(J)-ARXJ(J)GO TO 6867 DO 66 J=3,M3ARXJ(J)=0.5*ARX(J)66 ARXJP(J)=ARXJ(J)68 ARXJP(2)=ARX(2)ARXJ(M2)=ARX(M2)DO 70 J=3,M3FV(J)=ARXJP(J)/ARX(J)70 FVP(J)=1.-FV(J)DO 85 I=3,L2FX(I)=0.5*XCV(I-1)/XDIF(I)85 FXM(I)=1.-FX(I)FX(2)=0.FXM(2)=1.FX(L1)=1.FXM(L1)=0.DO 90 J=3,M2FY(J)=0.5*YCV(J-1)/YDIF(J)90 FYM(J)=1.-FY(J)FY(2)=0.FYM(2)=1.FY(M1)=1.FYM(M1)=0.*---CON,AP,U,V,RHO,PC AND P ARRAYS ARE INITIALIZED HERE---- DO 95 J=1,M1DO 95 I=1,L1PC(I,J)=0.U(I,J)=0.V(I,J)=0.CON(I,J)=0.AP(I,J)=0.RHO(I,J)=RHOCONP(I,J)=0.95 CONTINUEIF(MODE .EQ. 1) WRITE(8,1)IF(MODE .EQ. 2) WRITE(8,2)IF(MODE .EQ. 3) WRITE(8,3)WRITE(8,4)RETURN*----------------------------------------------------------------------ENTRY SETUP2*---COEFFICIENTS FOR THE U EQUATION----NF=1IF(.NOT. LSOLVE(NF)) GO TO 100IST=3JST=2CALL GAMSORREL=1.-RELAX(NF)DO 102 I=3,L2FL=XCVI(I)*V(I,2)*RHO(I,1)FLM=XCVIP(I-1)*V(I-1,2)*RHO(I-1,1)FLOW=R(1)*(FL+FLM)DIFF=R(1)*(XCVI(I)*GAM(I,1)+XCVIP(I-1)*GAM(I-1,1))/YDIF(2)CALL DIFLOW102 AJM(I,2)=ACOF+AMAX1(0.,FLOW)DO 103 J=2,M2FLOW=ARX(J)*U(2,J)*RHO(1,J)DIFF=ARX(J)*GAM(1,J)/(XCV(2)*SX(J))CALL DIFLOWAIM(3,J)=ACOF+AMAX1(0.,FLOW)DO 103 I=3,L2IF(I .EQ. L2) GO TO 104FL=U(I,J)*(FX(I)*RHO(I,J)+FXM(I)*RHO(I-1,J))FLP=U(I+1,J)*(FX(I+1)*RHO(I+1,J)+FXM(I+1)*RHO(I,J))FLOW=ARX(J)*0.5*(FL+FLP)DIFF=ARX(J)*GAM(I,J)/(XCV(I)*SX(J))GO TO 105104 FLOW=ARX(J)*U(L1,J)*RHO(L1,J)DIFF=ARX(J)*GAM(L1,J)/(XCV(L2)*SX(J))105 CALL DIFLOWAIM(I+1,J)=ACOF+AMAX1(0.,FLOW)AIP(I,J)=AIM(I+1,J)-FLOWIF (J .EQ. M2) GOTO 106FL=XCVI(I)*V(I,J+1)*(FY(J+1)*RHO(I,J+1)+FYM(J+1)*RHO(I,J))FLM=XCVIP(I-1)*V(I-1,J+1)*(FY(J+1)*RHO(I-1,J+1)+FYM(J+1)*& RHO(I-1,J))GM=GAM(I,J)*GAM(I,J+1)/(YCV(J)*GAM(I,J+1)+YCV(J+1)*GAM(I,J)+ & 1.0E-30)*XCVI(I)GMM=GAM(I-1,J)*GAM(I-1,J+1)/(YCV(J)*GAM(I-1,J+1)+YCV(J+1)*& GAM(I-1,J)+1.E-30)*XCVIP(I-1)DIFF=RMN(J+1)*2.*(GM+GMM)GO TO 107106 FL=XCVI(I)*V(I,M1)*RHO(I,M1)FLM=XCVIP(I-1)*V(I-1,M1)*RHO(I-1,M1)DIFF=R(M1)*(XCVI(I)*GAM(I,M1)+XCVIP(I-1)*GAM(I-1,M1))/YDIF(M1) 107 FLOW=RMN(J+1)*(FL+FLM)CALL DIFLOWAJM(I,J+1)=ACOF+AMAX1(0.,FLOW)AJP(I,J)=AJM(I,J+1)-FLOWVOL=YCVR(J)*XCVS(I)APT=(RHO(I,J)*XCVI(I)+RHO(I-1,J)*XCVIP(I-1))& /(XCVS(I)*DT)AP(I,J)=AP(I,J)-APTCON(I,J)=CON(I,J)+APT*U(I,J)AP(I,J)=(-AP(I,J)*VOL+AIP(I,J)+AIM(I,J)+AJP(I,J)+AJM(I,J))& /RELAX(NF)CON(I,J)=CON(I,J)*VOL+REL*AP(I,J)*U(I,J)DU(I,J)=VOL/(XDIF(I)*SX(J))CON(I,J)=CON(I,J)+DU(I,J)*(P(I-1,J)-P(I,J))DU(I,J)=DU(I,J)/AP(I,J)103 CONTINUECALL SOLVE100 CONTINUE*---COEFFICIENTS FOR THE V EQUATION----NF=2IF(.NOT. LSOLVE(NF)) GO TO 200IST=2JST=3CALL GAMSORREL=1.-RELAX(NF)DO 202 I=2,L2AREA=R(1)*XCV(I)FLOW=AREA*V(I,2)*RHO(I,1)DIFF=AREA*GAM(I,1)/YCV(2)CALL DIFLOW202 AJM(I,3)=ACOF+AMAX1(0.,FLOW)DO 203 J=3,M2FL=ARXJ(J)*U(2,J)*RHO(1,J)FLM=ARXJP(J-1)*U(2,J-1)*RHO(1,J-1)FLOW=FL+FLMDIFF=(ARXJ(J)*GAM(1,J)+ARXJP(J-1)*GAM(1,J-1))/(XDIF(2)*SXMN(J)) CALL DIFLOWAIM(2,J)=ACOF+AMAX1(0.,FLOW)DO 203 I=2,L2IF(I .EQ. L2)GO TO 204FL=ARXJ(J)*U(I+1,J)*(FX(I+1)*RHO(I+1,J)+FXM(I+1)*RHO(I,J))FLM=ARXJP(J-1)*U(I+1,J-1)*(FX(I+1)*RHO(I+1,J-1)+FXM(I+1)*& RHO(I,J-1))GM=GAM(I,J)*GAM(I+1,J)/(XCV(I)*GAM(I+1,J)+XCV(I+1)*GAM(I,J)+ & 1.E-30)*ARXJ(J)GMM=GAM(I,J-1)*GAM(I+1,J-1)/(XCV(I)*GAM(I+1,J-1)+XCV(I+1)*& GAM(I,J-1)+1.0E-30)*ARXJP(J-1)DIFF=2.*(GM+GMM)/SXMN(J)GO TO 205204 FL=ARXJ(J)*U(L1,J)*RHO(L1,J)FLM=ARXJP(J-1)*U(L1,J-1)*RHO(L1,J-1)DIFF=(ARXJ(J)*GAM(L1,J)+ARXJP(J-1)*GAM(L1,J-1))/(XDIF(L1)*SXMN(J)) 205 FLOW=FL+FLMCALL DIFLOWAIM(I+1,J)=ACOF+AMAX1(0.,FLOW)AIP(I,J)=AIM(I+1,J)-FLOWIF(J .EQ. M2) GO TO 206AREA=R(J)*XCV(I)FL=V(I,J)*(FY(J)*RHO(I,J)+FYM(J)*RHO(I,J-1))*RMN(J)FLP=V(I,J+1)*(FY(J+1)*RHO(I,J+1)+FYM(J+1)*RHO(I,J))*RMN(J+1)FLOW=(FV(J)*FL+FVP(J)*FLP)*XCV(I)DIFF=AREA*GAM(I,J)/YCV(J)GO TO 207206 AREA=R(M1)*XCV(I)FLOW=AREA*V(I,M1)*RHO(I,M1)DIFF=AREA*GAM(I,M1)/YCV(M2)207 CALL DIFLOWAJM(I,J+1)=ACOF+AMAX1(0.,FLOW)AJP(I,J)=AJM(I,J+1)-FLOWVOL=YCVRS(J)*XCV(I)SXT=SX(J)IF(J .EQ. M2) SXT=SX(M1)SXB=SX(J-1)IF(J .EQ. 3) SXB=SX(1)APT=(ARXJ(J)*RHO(I,J)*0.5*(SXT+SXMN(J))+ARXJP(J-1)*RHO(I,J-1)*& 0.5*(SXB+SXMN(J)))/(YCVRS(J)*DT)AP(I,J)=AP(I,J)-APTCON(I,J)=CON(I,J)+APT*V(I,J)AP(I,J)=(-AP(I,J)*VOL+AIP(I,J)+AIM(I,J)+AJP(I,J)+AJM(I,J))& /RELAX(NF)CON(I,J)=CON(I,J)*VOL+REL*AP(I,J)*V(I,J)DV(I,J)=VOL/YDIF(J)CON(I,J)=CON(I,J)+DV(I,J)*(P(I,J-1)-P(I,J))DV(I,J)=DV(I,J)/AP(I,J)203 CONTINUECALL SOLVE200 CONTINUE*---COEFIICIENTS FOR THE PRESSURE CORRECTION EQUATION----NF=3IF(.NOT. LSOLVE(NF)) GO TO 500IST=2JST=2CALL GAMSORSMAX=0.SSUM=0.DO 410 J=2,M2DO 410 I=2,L2VOL=YCVR(J)*XCV(I)410 CON(I,J)=CON(I,J)*VOLDO 402 I=2,L2ARHO=R(1)*XCV(I)*RHO(I,1)CON(I,2)=CON(I,2)+ARHO*V(I,2)402 AJM(I,2)=0.DO 403 J=2,M2ARHO=ARX(J)*RHO(1,J)CON(2,J)=CON(2,J)+ARHO*U(2,J)AIM(2,J)=0.DO 403 I=2,L2IF(I .EQ. L2) GO TO 404ARHO=ARX(J)*(FX(I+1)*RHO(I+1,J)+FXM(I+1)*RHO(I,J))FLOW=ARHO*U(I+1,J)CON(I,J)=CON(I,J)-FLOWCON(I+1,J)=CON(I+1,J)+FLOWAIP(I,J)=ARHO*DU(I+1,J)AIM(I+1,J)=AIP(I,J)GO TO 405404 ARHO=ARX(J)*RHO(L1,J)CON(I,J)=CON(I,J)-ARHO*U(L1,J)AIP(I,J)=0.405 IF(J .EQ. M2) GO TO 406ARHO=RMN(J+1)*XCV(I)*(FY(J+1)*RHO(I,J+1)+FYM(J+1)*RHO(I,J)) FLOW=ARHO*V(I,J+1)CON(I,J)=CON(I,J)-FLOWCON(I,J+1)=CON(I,J+1)+FLOWAJP(I,J)=ARHO*DV(I,J+1)AJM(I,J+1)=AJP(I,J)GO TO 407406 ARHO=RMN(M1)*XCV(I)*RHO(I,M1)CON(I,J)=CON(I,J)-ARHO*V(I,M1)AJP(I,J)=0.407 AP(I,J)=AIP(I,J)+AIM(I,J)+AJP(I,J)+AJM(I,J)PC(I,J)=0.SMAX=AMAX1(SMAX,ABS(CON(I,J)))SSUM=SSUM+CON(I,J)403 CONTINUECALL SOLVE*---COMEE HERE TO CORRECT THE PRESSURE AND VELOCITIES DO 501 J=2,M2DO 501 I=2,L2P(I,J)=P(I,J)+PC(I,J)*RELAX(NP)IF(I .NE. 2) U(I,J)=U(I,J)+DU(I,J)*(PC(I-1,J)-PC(I,J))IF(J .NE. 2) V(I,J)=V(I,J)+DV(I,J)*(PC(I,J-1)-PC(I,J))501 CONTINUE500 CONTINUE*---COEFFICIENTS FOR OTHER EQUATIONS----IST=2JST=2DO 600 N=4,NFMAXNF=NIF(.NOT. LSOLVE(NF)) GO TO 600CALL GAMSORREL=1.-RELAX(NF)DO 602 I=2,L2AREA=R(1)*XCV(I)FLOW=AREA*V(I,2)*RHO(I,1)DIFF=AREA*GAM(I,1)/YDIF(2)CALL DIFLOW602 AJM(I,2)=ACOF+AMAX1(0.,FLOW)DO 603 J=2,M2FLOW=ARX(J)*U(2,J)*RHO(1,J)DIFF=ARX(J)*GAM(1,J)/(XDIF(2)*SX(J))CALL DIFLOWAIM(2,J)=ACOF+AMAX1(0.,FLOW)DO 603 I=2,L2IF(I .EQ. L2) GO TO 604FLOW=ARX(J)*U(I+1,J)*(FX(I+1)*RHO(I+1,J)+FXM(I+1)*RHO(I,J))DIFF=ARX(J)*2.*GAM(I,J)*GAM(I+1,J)/((XCV(I)*GAM(I+1,J)+& XCV(I+1)*GAM(I,J)+1.0E-30)*SX(J))GO TO 605604 FLOW=ARX(J)*U(L1,J)*RHO(L1,J)DIFF=ARX(J)*GAM(L1,J)/(XDIF(L1)*SX(J))605 CALL DIFLOWAIM(I+1,J)=ACOF+AMAX1(0.,FLOW)AIP(I,J)=AIM(I+1,J)-FLOWAREA=RMN(J+1)*XCV(I)IF(J .EQ. M2) GO TO 606FLOW=AREA*V(I,J+1)*(FY(J+1)*RHO(I,J+1)+FYM(J+1)*RHO(I,J))DIFF=AREA*2.*GAM(I,J)*GAM(I,J+1)/(YCV(J)*GAM(I,J+1)+& YCV(J+1)*GAM(I,J)+1.0E-30)GO TO 607606 FLOW=AREA*V(I,M1)*RHO(I,M1)DIFF=AREA*GAM(I,M1)/YDIF(M1)607 CALL DIFLOWAJM(I,J+1)=ACOF+AMAX1(0.,FLOW)AJP(I,J)=AJM(I,J+1)-FLOWVOL=YCVR(J)*XCV(I)APT=RHO(I,J)/DTAP(I,J)=AP(I,J)-APTCON(I,J)=CON(I,J)+APT*F(I,J,NF)AP(I,J)=(-AP(I,J)*VOL+AIP(I,J)+AIM(I,J)+AJP(I,J)+AJM(I,J))& /RELAX(NF)CON(I,J)=CON(I,J)*VOL+REL*AP(I,J)*F(I,J,NF)603 CONTINUECALL SOLVE600 CONTINUETIME=TIME+DTITER=ITER+1IF(ITER .GE. LAST) LSTOP=.TRUE.RETURNEND*-----------------------------------------------------------------------SUBROUTINE SUPPLY************************************************************************ DOUBLE PRECISION TITLELOGICAL LSOLVE,LPRINT,LBLK,LSTOPCOMMON F(22,22,10),P(22,22),RHO(22,22),GAM(22,22),CON(22,22),& AIP(22,22),AIM(22,22),AJP(22,22),AJM(22,22),AP(22,22),& X(22),XU(22),XDIF(22),XCV(22),XCVS(22),& Y(22),YV(22),YDIF(22),YCV(22),YCVS(22),&YCVR(22),YCVRS(22),ARX(22),ARXJ(22),ARXJP(22),&R(22),RMN(22),SX(22),SXMN(22),XCVI(22),XCVIP(22)COMMON DU(22,22),DV(22,22),FV(22),FVP(22),&FX(22),FXM(22),FY(22),FYM(22),PT(22),QT(22)COMMON/INDX/NF,NFMAX,NP,NRHO,NGAM,L1,L2,L3,M1,M2,M3,& IST,JST,ITER,LAST,TITLE(13),RELAX(13),TIME,DT,XL,YL,& IPREF,JPREF,LSOLVE(10),LPRINT(13),LBLK(10),MODE,NTIMES(10),RHOCONDIMENSION U(22,22),V(22,22),PC(22,22)EQUIVALENCE (F(1,1,1),U(1,1)),(F(1,1,2),V(1,1)),(F(1,1,3),PC(1,1))************************************************************************10 FORMAT(1X,26(1H*),3X,A10,3X,26(1H*))20 FORMAT(1X,4H I =,I6,6I9)30 FORMAT(1X,1HJ)40 FORMAT(1X,I2,3X,1P7E9.2)50 FORMAT(1X,1H )51 FORMAT(1X,' I =',2X,7(I4,5X))52 FORMAT(1X,' X =',1P7E9.2)53 FORMAT(1X,'TH =',1P7E9.2)54 FORMAT(1X,'J =',2X,7(I4,5X))55 FORMAT(1X,'Y =',1P7E9.2)************************************************************************ ENTRY UGRIDXU(2)=0.DX=XL/FLOAT(L1-2)DO 1 I=3,L11 XU(I)=XU(I-1)+DXYV(2)=0.DY=YL/FLOAT(M1-2)DO 2 J=3,M12 YV(J)=YV(J-1)+DYRETURN************************************************************************ ENTRY PRINTIF(.NOT. LPRINT(3)) GO TO 80*---CALCULATE THE STREAM FUNTION----------------------------------------F(2,2,3)=0.DO 82 I=2,L1IF(I .NE. 2) F(I,2,3)=F(I-1,2,3)-RHO(I-1,1)*V(I-1,2)& *R(1)*XCV(I-1)DO 82 J=3,M1RHOM=FX(I)*RHO(I,J-1)+FXM(I)*RHO(I-1,J-1)82 F(I,J,3)=F(I,J-1,3)+RHOM*U(I,J-1)*ARX(J-1)80 CONTINUE*IF( .NOT. LPRINT(NP)) GO TO 90**---CONSTRUCT BOUNDARY PRESSURES BY EXTRAPOLATIONDO 91 J=2,M2P(1,J)=(P(2,J)*XCVS(3)-P(3,J)*XDIF(2))/XDIF(3)91 P(L1,J)=(P(L2,J)*XCVS(L2)-P(L3,J)*XDIF(L1))/XDIF(L2)DO 92 I=2,L2P(I,1)=(P(I,2)*YCVS(3)-P(I,3)*YDIF(2))/YDIF(3)92 P(I,M1)=(P(I,M2)*YCVS(M2)-P(I,M3)*YDIF(M1))/YDIF(M2)P(1,1)=P(2,1)+P(1,2)-P(2,2)P(L1,1)=P(L2,1)+P(L1,2)-P(L2,2)P(1,M1)=P(2,M1)+P(1,M2)-P(2,M2)P(L1,M1)=P(L2,M1)+P(L1,M2)-P(L2,M2)PREF=P(IPREF,JPREF)DO 93 J=1,M1DO 93 I=1,L193 P(I,J)=P(I,J)-PREF90 CONTINUE*IF(TIME.GT.0.5*DT) GOTO 320WRITE (8,50)IEND=0301 IF(IEND .EQ. L1) GO TO 310IBEG=IEND+1IEND=IEND+7IEND=MIN0(IEND,L1)WRITE (8,50)WRITE(8,51) (I,I=IBEG,IEND)IF(MODE .EQ. 3) GO TO 302WRITE(8,52) (X(I),I=IBEG,IEND)GO TO 303302 WRITE (8,53) (X(I),I=IBEG,IEND)303 GO TO 301310 JEND=0WRITE(8,50)311 IF(JEND .EQ. M1) GO TO 320JBEG=JEND+1JEND=JEND+7JEND=MIN0(JEND,M1)WRITE(8,50)WRITE(8,54) (J,J=JBEG,JEND)WRITE(8,55) (Y(J),J=JBEG,JEND)GO TO 311320 CONTINUE*DO 999 N=1,NGAMNF=NIF(.NOT. LPRINT(NF)) GO TO 999WRITE(8,50)WRITE(8,10) TITLE(NF)IFST=1JFST=1IF(NF .EQ. 1 .OR. NF .EQ. 3) IFST=2IF(NF .EQ. 2 .OR. NF .EQ. 3) JFST=2IBEG=IFST-7110 CONTINUEIBEG=IBEG+7IEND=IBEG+6IEND=MIN0(IEND,L1)WRITE(8,50)WRITE(8,20) (I,I=IBEG,IEND)WRITE(8,30)JFL=JFST+M1DO 115 JJ=JFST,M1J=JFL-JJWRITE(8,40) J,(F(I,J,NF),I=IBEG,IEND) 115 CONTINUEIF(IEND .LT. L1) GO TO 110999 CONTINUERETURNEND。

相关文档
最新文档