fortran+grads
fortran指令大全
附录C SCILAB 部分函数指令表(c)LIAMA. All rights reserved.(注解:本指令表只收集了部分常用指令, 有关全部指令请参照文档文件)+ 加- 减* 矩阵乘数组乘*.1. 通用指令^ 矩阵乘方数组乘方^.\ 反斜杠或左除help 在线帮助/ 斜杠或右除apropos 文档中关键词搜寻或.\ 数组除/.ans 缺省变量名以及最新表达式的运算结果== 等号~= 不等号clear 从内存中清除变量和函数< 小于exit 关闭SCILAB> 大于quit 退出SCILAB<= 小于或等于save 把内存变量存入磁盘>= 大于或等于exec 运行脚本文件&,and 逻辑与mode 文件运行中的显示格式|,or 逻辑或getversion 显示SCILAB 版本~,not 逻辑非ieee 浮点运算溢出显示模式选择: 冒号who 列出工作内存中的变量名( ) 园括号edit 文件编辑器[ ] 方括号type 变量类型{ } 花括号what 列出SCILAB 基本命令小数点.format 设置数据输出格式, 逗号chdir 改变当前工作目录; 分号getenv 给出环境值// 注释号mkdir 创建目录= 赋值符号pwd 显示当前工作目录' 引号evstr 执行表达式' 复数转置号转置号'.ans 最新表达式的运算结果2.运算符和特殊算符%eps 浮点误差容限, =2-52≈2.22×10-16%i 虚数单位= √(-1)%inf 正无穷大%pi 圆周率,π=3.1415926535897....3. 编程语言结构abort 中止计算或循环break 终止最内循环case 同select 一起使用continue 将控制转交给外层的for或while循环else 同if一起使用elseif 同if一起使用end 结束for,while,if 语句for 按规定次数重复执行语句if 条件执行语句otherwise 可同switch 一起使用pause 暂停模式return 返回select 多个条件分支then 同if一起使用while 不确定次数重复执行语句eval 特定值计算feval 函数特定值计算或多变量计算function 函数文件头global 定义全局变量isglobal 检测变量是否为全局变量error 显示错误信息lasterror 显示最近的错误信息sprintf 按格式把数字转换为串warning 显示警告信息4.基本数学函数acos 反余弦acosh 反双曲余弦acot 反余切acoth 反双曲余切acsc 反余割acsch 反双曲余割asin 反正弦asinh 反双曲正弦atan 反正切atanh 反双曲正切cos 余弦cosh 双曲余弦cotg 余切coth 双曲余切sin 正弦sinh 双曲正弦tan 正切tanh 双曲正切exp 指数log 自然对数log10 常用对数log2 以2为底的对数sqrt 平方根abs 绝对值conj 复数共轭imag 复数虚部real 复数实部ceil 向上(正无穷大方向)取整fix 向零方向取整floor 向下(负无穷大方向)取整round 四舍五入取整sign 符号函数gsort 降次排序erf 误差函数erfc 补误差函数gamma gamma 函数interp 插值函数interpln 线性插值函数intsplin 样条插值函数smooth 样条平滑函数spline 样条函数quarewave 方波函数sign 符号函数double 将整数转换为双精度浮点数5.基本矩阵函数和操作eye 单位阵zeros 全零矩阵ones 全1 矩阵rand 均匀分布随机阵genmarkov 生成随机Markov 矩阵linspace 线性等分向量logspace 对数等分向量logm 矩阵对数运算cumprod 矩阵元素累计乘cumsum 矩阵元素累计和toeplitz Toeplitz 矩阵disp 显示矩阵和文字内容length 确定向量的长度size 确定矩阵的维数diag 创建对角阵或抽取对角向量find 找出非零元素1的下标matrix 矩阵变维rot90 矩阵逆时针旋转90度sub2ind 据全下标换算出单下标tril 抽取下三角阵triu 抽取上三角阵conj 共轭矩阵companion 伴随矩阵det 行列式的值norm 矩阵或向量范数nnz 矩阵中非零元素个数null 清空向量或矩阵中的某个元素orth 正交基rank 矩阵秩trace 矩阵迹cond 矩阵条件数rcond 逆矩阵条件数inv 矩阵的逆lu LU分解或高斯消元法pinv 伪逆qr QR分解givens Givens 变换linsolve 求解线性方程lyap Lyapunov 方程hess Hessenberg 矩阵poly 特征多项式schur Schur 分解expm 矩阵指数expm1 矩阵指数的Pade逼近expm2 用泰勒级数求矩阵指数expm3 通过特征值和特征向量求矩阵指数funm 计算一般矩阵函数logm 矩阵对数sqrtm 矩阵平方根6. 特性值与奇异值spec 矩阵特征值gspec 矩阵束特征值bdiag 块矩阵, 广义特征向量eigenmarkov 正则化Markov 特征向量pbig 特征空间投影svd 奇异值分解sva 奇异值分解近似7. 矩阵元素运算cumprod 元素累计积cumsum 元素累计和hist 统计频数直方图max 最大值mean 平均值median 中值min 最小值prod 元素积sort 由大到小排序std 标准差sum 元素和trapz 梯形数值积分corr 求相关系数或方差8. 稀疏矩阵运算sparse 稀疏矩阵(只存储非零元素)adj2sp 邻接矩阵转换为稀疏矩阵full 稀疏矩阵转换为全矩阵mtlb_sparse 将SCILAB 稀疏矩阵转换为MA TLAB稀疏矩阵格式sp2adj 稀疏矩阵转换为邻接矩阵speye 稀疏矩阵方式单位阵sprand 稀疏矩阵方式随机矩阵spzeros 稀疏矩阵方式全零阵lufact 稀疏矩阵LU分解lusolve 稀疏矩阵方程求解spchol 稀疏矩阵Cholesky分解9. 输入输出函数diary 生成屏幕文本记录disp 变量显示file 文件管理input 用户键盘输入load 读已存的变量mclose 关闭文件mget 读二进制文件mgetl 按行读ASCII码文件mgetstr 读字符串中单个字mopen 打开文件mput 写二进制文件mfscanf 读ASCII 码文件print 将变量记录为文件read 读矩阵变量save 存变量为二进制文件strartup 启动文件write 按格式存文件xgetfile 对话方式获取文件路径x_dialog 建立Xwindow参数输入对话框Tk_Getvar 得到Tk文件变量Tk_EvalFile 执行Tk 文件10. 函数与函数库操作deff 在线定义函数edit 函数编辑器function 打开函数定义functions SCILAB 函数或对象genlib 在给定目录下建立所有文件的函数库get_function_path 读函数库的文件存储目录路径getd 读函数库中的全部文件getf 在文件中定义一个函数lib 函数库定义macro SCILAB函数或对象macrovar 输入变量个数newfun 输出变量个数11. 字符串操作code2str 将SCILAB数码转换为字符串convstr 字母大小转换emptystr 清空字符串grep 搜寻相同字符串part 字符提取str2code 将字符串转换为SCILAB数码string 字符串转换strings SCILAB 对象, 字符串strcat 连接字符strindex 字符串的字符位置搜寻strsubst 字符串中的字符替换12. 日期与时间date 日期getdate 读日期与时间timer CPU时间计时13. 二维图形函数plot2d 直角坐标下线性刻度曲线champ 2 维向量场champ1 由颜色箭头表示的2维向量场contour2d 等高线图errbar 曲线上增加误差范围框线条grayplot 应用颜色表示的表面xgrid 画坐标网格线histplot 统计频数直方图Matplot 散点图阵列14. 三维图形函数plot3d 三维表面plot3d1 用颜色或灰度表示的三维表面param3d 三维中单曲线param3d1 三维中多曲线contour 三维表面上的等高线图hist3d 三维表示的统计频数直方图geom3d 三维向二维上的投影15. 线条类图形xpoly 单线条或单多边形xpolys 多线条或多各多边形xrpoly 正多边形xsegs 非连接线段xfpoly 单个多边形内填充xfpolys 多个多边形内填充xrect 矩形xfrect 单个矩形内填充xrects 多个矩形内填充xarc 单个弧线段或弧园xarcs 多个弧线段或弧园xfarc 单个弧线段或弧园填充xfarcs 多个弧线段或弧园填充xarrows 多箭头16. 图形注释, 变换xstring 图形中字符xstringb 框内字符xtitle 图形标题xaxis 轴名标注plotframe 图形加框并画坐标网格线isoview 等尺寸比例显示(原图形窗口不改变)square 等尺寸比例显示(原图形窗口改变)xsetech 设置小窗口xchange 转换实数为图形象素坐标值subplot 设置多个子窗口17. 图形颜色及图形文字colormap 应用颜色图getcolor 交互式选择颜色图addcolor 增加新色于颜色图graycolormap 线性灰度图hotcolormap 热色(红到黄色)颜色图xset 图形显示方式设定xget 读当前图形显示方式设定getsymbol 交互式选择符号和尺寸18. 图形文件及图形文字xsave 将图形存储为文件xload 从磁盘中读出图形文件xbasimp 将图形按PS文件打印或存储为文件xs2fig 将图形生成Xfig 格式文件xbasc 取消图形窗及其相关内容xclear 清空图形窗driver 选择图形驱动器xinit 图形驱动器初始化xend 关闭图形xbasr 图形刷新replot 更改显示范围后的图形刷新xdel 关闭图形xname 改变当前图形窗名称19. 控制分析用图形bode 伯德图坐标gainplot 幅值图坐标(伯德图中的幅值图) nyquist 奈奎斯特图m_circle M-圆图chart 尼库拉斯图black Black-图evans 根轨迹图sgrid s 平面图plzr 零-极点图zgrid z 平面图20. 图形应用中的其它指令graphics 图形库指令表xclick 等待鼠标在图形上的点击输入locate 由鼠标点击读入图形中的多点位置坐标xgetmouse 由鼠标点击读入图形中的当前点位置坐标21. 系统与控制abcd 状态空间矩阵cont_mat 可控矩阵csim 线性系统时域响应dsimul 状态空间的离散时域响应feedback 反馈操作符flts 时域响应(离散、采样系统〕frep2tf 基于传递函数的频域响应freq 频域响应g_margin 幅值裕量imrep2ss 基于状态空间的脉冲响应lin 线性化操作lqe Kalman 滤波器lqg LQG补偿器lqr LQ补偿器ltitr 基于状态空间的离散时域响应obscont 基于观测器的控制器observer 观测器obsv_mat 观测矩阵p_margin 相位裕量phasemag 相位与幅值计算ppol 极点配置repfreq 频域响应ricc Riccati 方程rtitr 基于传递函数的离散时域响应sm2ss 系统矩阵到状态空间变换ss2ss 反馈连接的状态空间到状态空间变换ss2tf 状态空间到传递函数变换stabil 稳定性计算tf2ss 传递函数到状态空间变换time_id SISO系统最小方差辨识22. 鲁棒控制augment 被控对象增广操作bstap Hankel 矩阵近似ccontrg H∞控制器dhnorm 离散H∞范数h2norm H2 范数h_cl 闭环矩阵h_inf H∞控制器h_norm H∞范数hankelsv Hankel 矩阵奇异值leqr H∞控制器的LQ增益linf 无穷范数riccati Riccati 矩阵sensi 敏感函数23. 动态系统arma ARMA模型arma2p 基于AR模型中获得多项式矩阵armac ARMAX 辨识arsimul ARMAX系统仿真noisegen 噪声信号发生器odedi 常微分方程仿真检测prbs_a 伪随机二进制序列发生器reglin 线性拟合24. 系统与控制实例artest Arnold 动态系统bifish 鱼群人口发展的离散时域模型boucle 具有观测器的动态系统相位图chaintest 生物链模型gpech 渔业模型fusee 登陆火箭问题lotest Lorennz 吸引子mine 采矿问题obscontl可控可观系统portr3d 三维相位图portrait 二维相位图recur 双线性回归方程systems 动态系统tangent 动态系统的线性化tadinit 动态系统的交互初始化25. 非线性工具(优化与仿真〕bvode 边界值问题的常微分方程dasrt 隐式微分方程过零解dassl 代数微分方程datafit 基于测量数据的参数辨识derivative 导数计算fsolve 非线性函数过零解impl 线性微分方程int2d 二维定积分int3d 三维定积分intg 不定积分leastsq 非线性最小二乘法linpro 线性规划lmisolver 线性不等矩阵ode 常微分方程ode_discrete 离散常微分方程ode_root 常微分方程根解odedc 连续/离散常微分方程optim 非线性优化quapro 线性二次型规划semidef 半正定规划26. 多项式计算coeff 多项式系数coffg 多项式矩阵逆degree 多项式阶数denom 分母项derivat 有理矩阵求导determ 矩阵行列式值factors 因式分解hermit Hermit 型horner 多项式计算invr 有理矩阵逆lcm 最小公倍数ldiv 多项式矩阵长除numer 分子项pdiv 多项式矩阵除pol2des 多项式矩阵到表达式变换pol2str 多项式到字符串变换polfact 最小因式residu 余量roots 多项式根simp 多项式化简systmat 系统矩阵27. 信号处理%asn 椭圆积分%k Jacobi完全椭圆积分%sn Jacobi 椭圆函数analpf 模拟量低通滤波器buttmag Butterworth 滤波器响应cepstrum 倒谱计算cheb1mag Chebyshev 一型响应cheb2mag Chebyshev 二型响应chepol Chebyshev 多项式convol 卷积corr 相关, 协方差cspect 谱估计(应用相关法)dft 离散富立叶变换fft 快速富立叶变换filter 滤波器建模fsfirlin FIR滤波器设计hank 协方差矩阵到Hankel矩阵变换hilb Hilbert 变换iir IIR数字滤波器intdec 信号采样率更改kalm Kalman 滤波器更新mese 最大熵谱估计mfft 多维快速富立叶变换mrfit 频率响应拟合phc Markov 过程srkf Kalman 滤波器平方根sskf 稳态Kalman 滤波器system 观测更新wfir 线性相位FIR滤波器weiener Weiener(维纳)滤波器window 对称窗函数yulewalk 最小二乘滤波器zpbutt Buthererworth 模拟滤波器zpch1 Chebyshev 模拟滤波器28. 音频信号analyze 音频信号频域图auread 读*.au 音频文件auwrite 写*.au 音频文件lin2mu 将线性信号转换为µ率码信号loadwave 取*.wav 音频文件mapsound 音频信号图示mu2lin 将µ率码信号转换为线性信号playsnd 音频信号播放savewave 存*.wav 音频文件wavread 读*.wav 音频文件wavwrite 写*.wav 音频文件29. 语言与数据转换工具ascii 字符串的ASCII码excel2sci 读ASCII 格式的Excel 文件fun2string 将SCILAB 函数生成ASCII 码mfile2sci 将MA TLAB 的M 格式文件转换为SCI格式文件mtlb_load 取MA TLAB第4版本文件中变量matlb_save 按MA TLAB 第 4 版本文件格式存变量pol2tex 将多项式转换为TeX格式sci2for 将SCILAB 函数转换为FORTRAN格式文件texprint 按TeX 格式输出SCILAB 对象translatepaths 将子目录下的所有MA TLAB 文件转换为SCI文件格式一个公式写成Fortran语言代码program baiduinteger::I,J,Nreal*8::Cr,Treal*8,dimension(:),allocatable ::P,XN=3!变量X的个数Cr=5.0d0!常量Cr,自己设定T=4.0d0!常量T,自己设定allocate(P(N),X(N))! =======读入变量X的值do I=1,Nwrite(*,*)"请输入第",I," 个变量的值:"read(*,*)X(I)enddo! =======读入变量X的值do I=1,NP(I)=(-4.2d0/Cr**2*X(I)+2.9/Cr)*Twrite(*,*)“第”,I," 个变量X对应结果:",P(I)enddoend。
fortran 教学大纲
fortran 教学大纲Fortran 教学大纲Fortran(Formula Translation)是一种面向科学和工程计算的编程语言。
它在计算机科学的历史中扮演了重要角色,被广泛应用于科学计算、数值分析和大规模计算等领域。
本文将探讨 Fortran 教学的大纲,以帮助初学者系统地学习和掌握这门编程语言。
一、引言在本节中,我们将介绍 Fortran 的起源和发展,以及它在科学计算领域的重要性。
我们将讨论 Fortran 的特点,如其面向数值计算和高性能计算的优势,以及它对于科学家和工程师的实际应用。
二、基本语法和数据类型在这一部分,我们将介绍 Fortran 的基本语法规则和常用数据类型。
我们将讨论变量的声明和赋值,运算符的使用,以及控制流语句如条件语句和循环语句的编写方法。
此外,我们还将介绍 Fortran 中的基本数据类型,如整数、实数和字符类型,并讨论它们的使用场景和注意事项。
三、数组和矩阵运算Fortran 是一种强大的数组和矩阵运算语言。
在这一部分,我们将学习如何声明和操作一维和多维数组,以及如何进行矩阵运算。
我们将介绍 Fortran 中的数组索引和切片操作,以及常用的矩阵运算函数。
此外,我们还将讨论数组和矩阵的内存布局和性能优化技巧。
四、函数和子程序函数和子程序是 Fortran 中的重要概念,它们可以帮助我们组织和重用代码。
在这一部分,我们将学习如何声明和调用函数,以及如何编写和调用子程序。
我们将介绍函数的返回值和参数传递方式,以及子程序的参数传递和变量作用域。
此外,我们还将讨论递归函数和模块化编程的技巧。
五、文件操作和输入输出在科学计算中,数据的读取和保存是非常重要的。
在这一部分,我们将学习如何使用 Fortran 进行文件操作和输入输出。
我们将介绍如何打开和关闭文件,以及如何读取和写入数据。
此外,我们还将讨论格式化输入输出和二进制文件的处理方式,以及异常处理和错误处理的方法。
IDL、GRAD、NCL绘图命令对应关系一览表
IDL、NCL、GRADS、MATLAB绘图命令对应关系一览表1.grads的数据文件与ncl的什么文件对应?grads只支持按照一定顺序存储的二进制数据文件,后缀名以.grd或者.dat或者.grb结束。
这种数据必须以时间为最外层,然后每个变量按照向量形式存储,每个变量由外向内的存储顺序是高度(或等压面)==》纬度==》经度。
这必须注意,否则画图容易出现一堆一堆的乱线条,这就说明你的数据没有按照grads的要求存储。
再看看ncl,可以说ncl支持绝大多数各种数据的读写,包括netcdf,hdf,以及二进制数据甚至ASCII码(如果说是十进制数据或许你会更熟悉),前两种数据一般都有头文件,不能用C 语言或者Fortran读取,都需要插件才可以读取,ncl可以直接读取,matlab中也可以直接读取NETCDF格式(.nc)的数据。
grads中可以读取按照说明存取的NC数据,这种数据必须又正确的时间说明,也就是说时间必须是真实的,有些模式模拟出来的数据grads的sdfopen 命令是打不开的,因为一般模式都是nonleap run,都是平年,没有设定闰年,造成了时间说明不真实,grads就会报错。
那么grads'如何使用NC数据呢?所以建议使用ncl转换数据,将NC数据,hdf数据或者十进制数据转换成grd数据,供grads使用。
这样说明是在是太空洞了,那么下面我举个例子吧。
eg1) 使用grads将netcdf数据转换成grd数据'reinit'var.1='air';var.2='hgt';var.3='uwnd'var.4='vwnd';var.5='omega';var.6='shum'j=6while (j<=6)'set fwrite/disk3/users/Rao_Jian/ERA-Interim-daily/entropy/'var.j'.daily.1979-2010.grd' 'set gxout fwrite'i=1979while (i<=2010)'sdfopen/disk3/users/lbq/ERA-Interim-daily/pressure/'var.j'.interim.'i'.nc'tt=1if(i=1980|i=1984|i=1988|i=1992|i=1996|i=2000|i=2004|i=2008)while (tt<=366)'set t 'ttzz=1while (zz<=37)'set x 1 240''set y 1 121''set z 'zz'd 'var.j''zz=zz+1endwhilett=tt+1endwhileelsewhile (tt<=365)'set t 'ttzz=1zz=1while (zz<=37)'set x 1 240''set y 1 121''set z 'zz'd 'var.j''zz=zz+1endwhilett=tt+1endwhileendifi=i+1name='/disk3/users/lbq/ERA-Interim-daily/pressure/'var.j'.interim.'i'.nc''close 1'endwhile'disable fwrite'j=j+1endwhileeg2)使用ncl将数据输出成二进制数据,grads可以使用(截取部分)patho = "/disk3/users/Rao_Jian/Hadley/"system("rm "+patho+"ev_ts.grd")system("rm "+patho+"ev_ts2.grd")system("rm "+patho+"ev_ts.txt")system("rm"+patho+"ev_ts2.txt")do nt=0,719fbindirwrite(patho+"ev_ts.grd",ev_ts(:,nt));写成二进制数据,供grads使用end dofbindirwrite(patho+"ev_ts2.grd",ev_ts(time|:,evn|:))asciiwrite(patho+"ev_ts.txt",ev_ts(time|:,evn|:));写成十进制数据,可以贴到EXCEL中使用asciiwrite(patho+"ev_ts2.txt",ev_ts);print(ev_ts(0,:))此外ncl中还有其他的读写函数,如fbinrecread,fbinrecwrite,fbinwrite,fbinreadncl读取netcdf3/4、hdf4、grib-1/2也是通过文件操作的eg3)ncl支持直接读取格式文件的读法fi = addfile(pathi+"HadISST_sst.nc","r")sst0 =fi->sst(960:1679,:,:) ; load 50 year data duringfrom 1950 to 2009注意:这类文件的后缀名一般为.nc /.hdf/ .grb/.hdfeos/.ccm2.grads中的描述文件与ncl中的什么对应描述文件(.ctl)是一个纯文本文件,我们有数据grads还是不能出图,需要一个描述文件来指定他的存储数据个数,维度(时间长度、层数和经纬度信息)。
气象预报预警系统开发中VB与GRADS、FORTRAN混合编程技术
fortran矩阵运算 -回复
fortran矩阵运算-回复Fortran矩阵运算是一种用于处理线性代数方程的编程语言。
它为数值计算和科学计算提供了强大的工具,特别适用于处理大规模的矩阵运算。
在这篇文章中,我们将逐步回答一些关于Fortran矩阵运算的常见问题和操作。
第一步,我们将介绍Fortran中的矩阵表示。
在Fortran中,矩阵可以用二维数组来表示,其中行和列的索引分别从1开始。
例如,一个3×3的矩阵可以用以下方式在Fortran中表示:fortranreal, dimension(3, 3) :: matrix在这个声明中,“real”表示矩阵元素的数据类型为实数类型,dimension(3, 3)表示矩阵的行数和列数分别为3.第二步,我们将讨论矩阵的初始化和赋值。
在Fortran中,我们可以通过循环将数值赋给矩阵的各个元素。
例如,以下代码将给矩阵赋予一些随机数值:fortrando i = 1, 3do j = 1, 3call random_number(matrix(i, j))end doend do在这个例子中,random_number()函数用于生成随机数,并将其赋值给矩阵的每个元素。
使用循环结构可以方便地对矩阵进行赋值和初始化。
第三步,我们将介绍矩阵的基本运算。
Fortran提供了一系列的内置函数和运算符,能够执行矩阵的加法、减法、乘法以及转置等操作。
例如,以下代码演示了如何计算两个矩阵的乘法:fortranreal, dimension(3, 3) :: matrix1, matrix2, result...result = matmul(matrix1, matrix2)在这个例子中,matmul()函数用于计算两个矩阵的乘积,并将结果存储在result矩阵中。
Fortran还提供了elemental运算符,它允许对矩阵的每个元素进行逐个操作,例如逐个相加或相乘。
第四步,我们将介绍矩阵的求解和行列式计算。
fortran编程的步骤
fortran编程的步骤Fortran编程的步骤一、引言Fortran(Formula Translation)是一种高级程序设计语言,特别适用于科学计算和数值计算。
本文将介绍Fortran编程的步骤,帮助初学者了解如何使用Fortran进行程序开发。
二、编写程序的基本步骤1. 确定程序的目标:在开始编写Fortran程序之前,需要明确程序的目标和需求。
确定程序的输入和输出,以及所需的计算或处理步骤。
这有助于确保编写的程序能够满足预期的功能和要求。
2. 设计算法和数据结构:根据程序的目标,设计合适的算法和数据结构。
算法描述了解决问题的步骤和逻辑,而数据结构则定义了程序中使用的数据类型和数据组织方式。
3. 编写代码:根据算法和数据结构的设计,开始编写Fortran代码。
Fortran使用特定的语法和语句结构,需要熟悉其语法规则和常用的编程技巧。
代码的编写应遵循良好的编码风格,包括适当的缩进、注释和命名规范。
4. 编译程序:编写完Fortran代码后,需要使用Fortran编译器将源代码转换成可执行的机器代码。
编译过程将检查代码中的语法错误和逻辑错误,并生成可执行文件。
Fortran编译器通常会提供丰富的编译选项,可以根据需要进行调整。
5. 调试和测试:编译成功后,可以对程序进行调试和测试。
调试是指查找和修复程序中的错误和问题,测试是指验证程序的正确性和性能。
调试和测试是编程过程中不可或缺的环节,可以使用调试器和测试框架等工具辅助进行。
6. 优化和性能调整:在程序调试和测试完成后,可以考虑对程序进行优化和性能调整。
优化旨在改进程序的执行效率和资源利用率,可以通过改进算法、调整编译选项和使用高级优化技术来实现。
7. 文档撰写:在编程过程中,应及时记录程序的设计和实现细节。
文档可以包括程序的功能描述、算法和数据结构的说明、代码注释和使用说明等。
良好的文档可以提高代码的可读性和可维护性,并方便其他人理解和使用程序。
fortran+grads
fortran+gradsWRFCould not chdir to home directory /home/walilai: Stale NFS file handle/usr/X11R6/bin/xauth: error in locking authority file/home/walilai/.XauthorityBugINF ! 1/0=無限大(可用程式寫寫看)NaN ! 1/0=無限大read(unit,*) ! 注意unit數字使用限制Segmentation fault!矩陣沒定義到的空間(I or j…)或從第6行開始寫(須從7)doi=i+1 ! 寫在這會出現Segmentation faultstation(i)=' ' ! 可能是if沒被定義到I,stationread(10,21,end=31) st,lon,lat,mm,dd,hh,rainif(mm.eq.7.and.dd.eq.13.and.hh.eq.24)theni=i+1 !OKstation(i)=' ' !OKdo j=1,idoi=i+1c station(i)=' ' !若沒c掉,會出現Segmentation fault read(11,110,end=97) st,yy,mm,dd,hh,c,raindo_ud: off end of record:k=0 !a.out OK!!c do j=1,a !c處出現do_ud: off end of record:c k=j-1 do i=1,awrite(11,rec=3*k+1)station(i)write(11,rec=3*k+2)st_lon(i)write(11,rec=3*k+3)st_lat(i)k=k+1enddoopen: null file name!沒有讀到檔案,或是原本要輸入end才可結束程式的地方案到其他按鍵list in: end of file!檔案讀不完全,可能status=unknownopen: 'new' file exists:status!因檔案存在,不可用newdo_ud: end of file!字串字數沒寫對apparent state: unit 10 named f060200_112.5_35.datfmt: read unexpected character !找以下檔案會發現錯誤apparent state: unit 11 named c-cwb200807.txtlast format: (A6,I4,I2,I2,I2,A47,F5.1)!所讀數字出現文字T(雨跡)!format格式讀錯!可能計算出現nan而失敗invalid integer: read unexpected characterapparent state: unit 87 named autop2010-07.txtlast format: list iolately reading sequential formatted external IOAborted! read(87,100,end=920)st,yy,mm,dd,hh,c,rain將格式由*改為100就可fmt: end of fileapparent state: unit 92 named remove_2.txtlast format: (A6,1X,F6.1,1X,F5.1,2X,F7.1,2X,F7.1)lately reading sequential formatted external IOAborted!沒有close(??)或是無線迴圈的end=900,或是開檔時只有11個,卻讀到do i=1,12,或是1D矩陣寫成3D矩陣,或是已開啟資料寫兩次read()list in: end of fileapparent state: unit 31 named fort.31last format: list iolately reading direct formatted external IOAborted!txt_name(31)='xxx'之後,後面出現read(31,*) txt_name(31) /tmp/ccI6mztx.o(.text+0x139d): In function `MAIN__': : undefined reference to `filt25_'副程式最後要加include ‘filt25.f’do_ud: off end of record !recl那nx和ny寫錯apparent state: unit 85 named grads_5.datstnmap -i xxx.ctlName of binary data set:/work4/tony1321/rain/200807/0710_daily.datNumber of times in the data set: 24Number of surface variables: 7Number of level dependent variables: 0Starting scan of station data binary file. Binary data file open: /work4/tony1321/rain/200807/0710_daily.datProcessing time = 1Invalid station hdr found in station binary filePossible causes: Invalid level count in hdrDescriptor file mismatchFile not station dataInvalid relative timelevs = -1027080205 flag = -998639206 time = 1.4013e-45 Processing time = 24 !okTime = 24 has stn count = 415Max reports per time: 415 reports at t = 1Max data elements in largest report: 5Version 2 Station map file created: rain.mapstnmap: WARNING!! This stnmap file can only be accessed by GrADS Version 1.9b4stnmap: WARNING!! However, GrADS Version 1.9b4 can read both versionsstnmap: COMMENT -- use the -1 command line option to create a map for older GrADS versionsgradsshaded會把contour蓋過去,因此要先寫shaded,再寫contourq file 1.2……:不需開ctl即可看到各種變數若ctl觀測為一段時間,沒再gs定義set t XX,他會自動跑下去Constant field. Value = 2.60788e-09'set font 0-5 調整字體'set cstyle 0-7 調整線的類型'set t 'tt 要空格(空格 = 等於)'set digsize 數字'設定風標長度'set mpdraw off' 將grads內建地圖關掉Grads裡都用小寫較不會有問題('d Tem' ok但'd 'Tem不行)會出現畫出來的圖亂讀,譬如讀溫度卻跑出壓力,那可能式氣壓層數寫錯,如有25層卻只寫到21層Define error: Define not yet valid for station dataDefault file is a station data file! station data必放在第一個Low Level I/O Error: Read error on data fileData file name = grads.datError reading 81 bytes at location 9882 !此處排序錯誤Data Request Error: Error for variable 'inttem'Error ocurred at column 1Low Level I/O Error: Read error on data fileData file name = 20091012.dat !只有一層但ctl卻寫21層Error reading 141 bytes at location 417501Data Request Error: Error for variable 'dbz'Error ocurred at column 1DISPLAY error: Invalid expressionExpression = dbzCannot plot color bar: No shading informationDISPLAY error: Invalid expression Expression = inttemData Request Warning: Request beyond file limitsCannot plot data - all undefined values !set dfile順序要放對Open Error: Can't open description file !開到不存在的ctl檔SET Error: No files open yetDISPLAY error: no file open yetSyntax Error: Invalid Operand !ctl檔名稱寫錯,沒讀到'ccslp' not a variable or function nameError ocurred at column 1DISPLAY error: Invalid expressionExpression = ccslpLow Level I/O Error: Read error on data fileData file name = 20091012_RCWF.datError reading 141 bytes at location 79524Data Request Error: Error for variable 'dbz'Error ocurred at column 1DISPLAY error: Invalid expressionExpression = dbzCannot plot color bar: No shading information! 只有一層z=1,卻寫到set z 5而讀錯Open Error: Can't open description fileSET Error: No files open yetSET Error: No files open yet !ctl檔名寫錯SET Error: No files open yetSET Error: No files open yetDISPLAY error: no file open yetCannot plot color bar: No shading informationOpen Error: Can't open description fileSyntax Error: Invalid Operand'u' not a variable or function nameError ocurred at column 1DISPLAY error: Invalid expressionExpression = u!gs裡ctl超過一個時,需每個都要正確存在,否則會出現上面這樣,也就是全部ctl的路徑要全寫和高度(層)有關的變數,用'uprs'代表比較好在畫ec資料時須注意1.Data Request Warning: Request beyond file limits'set t 'tt位置順序寫錯有兩個ctl以上若時間定義不同,地形(t=1)2.如果gradsnc 裡面 t while迴圈順序寫錯,他會亂讀ex:tt=1while(tt<=124)第一個ctl#tt=tt+1 <---這裡必須關掉迴圈,否則會亂讀#endwhile#tt=1#while(tt<=12)第二個ctltt=tt+1 endwhile <---這迴圈只能寫在這'd rain''set gxout stat''d rain'#print resultline=sublin(result,8)word=subwrd(line,5)rmax=substr(word,1,4)try try 看'draw title case1A 9-21 rainfall ('data')''printim case1a_9_21_rainfall.gif x1600 y1200 white'Set clab %gK 在數字後寫上單位KWarning from LOG: Data has 2 values <= zeroThese were set to the undefined value !同時定義到lev和zOpen Error: Inconsistent time countCount in station map file = 18Count in descriptor file = 1 !須做stnmap –I XXX.ctl The data file was not opened. !開啟的ctl時間沒有連續Cat完dat檔時網格點可用4*nx*ny*時間來驗證若已經有dfile時,則不可再對變數做類似var.2 or var.3SET error: Missing or invalid arguments for GXOUT option Data Request Error: File number out of range ! 開錯ctl檔Variable = ter.3Error ocurred at column 1DISPLAY error: Invalid expressionExpression = ter.3(t=1)Test2.fimplicit noneinteger ii,jj,i,j,a,bparameter (ii=2,jj=2,a=2,b=2)real x(ii,jj),y(ii,jj),z(ii,jj),w(ii,jj)open(10,file='test2.txt')do i=1,iiread(10,*)(x(i,j),j=1,jj)enddodo i=1,iic do j=1,jjc write(*,*)x(i,j) !只能寫成一排write(*,*)(x(i,j),j=1,jj) !可以寫成矩陣c enddoenddoc write(*,*)x(1,2)close(10)end programa.out11. 12. !x(1,2)=1222. 23. 991.fprogram qqimplicit noneinteger a,b,c(9,9)do a=1,9do b=1,9c(a,b)=a*benddoenddodo b=1,9write(*,*)(a,'×',b,'=',c(a,b),a=1,9) enddo! c(a,b) 寫成一列stopend./a.out !製作九九乘法1X 1= 1 2X 1= 2 3X 1= 3 4X 1= 4 5X 1= 5 6X 1= 6 7X 1= 7 8X 1= 8 9X 1= 91X 9= 9 2X 9= 18 3X 9= 27 4X 9= 36 5X 9= 45 6X 9= 54 7X 9= 63 8X 9= 72 9X 9= 81992.fProgram matrix tttttimplicit noneinteger a,b,c(9,9)open(6,file='992.txt',status='unknown',form='formatted')! 6=給予這個檔案的編號do b=1,9do a=1,9c(a,b)=a*benddoenddo! ------------------------------do b=1,9write(6,*)(a,'x',b,'=',c(a,b),a=1,9)enddo! 此處6是從6檔案寫出來stopendvi 992.txt1x 1= 1 2x 1= 2 3x 1= 3 4x 1= 4 5x 1= 5 6x 1= 6 7x 1= 7 8x 1= 8 9x 1= 91x 9= 9 2x 9= 18 3x 9= 27 4x 9= 36 5x 9= 45 6x 9= 54 7x 9= 63 8x 9= 72 9x 9= 81993.fProgram win or go homeimplicit noneinteger a,b,cwrite(*,*)"a="read(*,*)awrite(*,*)"b="read(*,*)bc=a*b write(*,*) a,'X',b,'=',cif(c<60 .or. a<10) then write(*,*) "go home" elsewrite(*,*) "win"endifstopend./a.outa=50b=3250X 32= 1600WinTest1.fProgram summer or playoff implicit none integer a,b,c,dreal :: e,fwrite(*,*)"San Antonio " write(*,*)"win game" read(*,*) ab=82-awrite(*,*)'lost game'write(*,*)be=a/bwrite(*,*)'winnig',ewrite(*,*)"Final Standings"write(*,*)a,' -',bwrite(*,*)"Phoenix "write(*,*)"win game"read(*,*)cd=82-cwrite(*,*)'lost game'write(*,*)df=c/dwrite(*,*)'winnig',fwrite(*,*)"Final Standings"write(*,*)c,' -',dif(e>f) thenwrite(*,*)'San Antonio in the playoff game' else if (e .eq. f) thenwrite(*,*)'The same record'elsewrite(*,*)'Phoneix in the playoff game'endifstopend ./a.outSan Antoniowin game50lost game32winnig 1.Final Standings50 - 32Phoenixwin game43lost game39winnig 1.Final Standings43 - 39The same recordFile.f 快速開啟許多檔案方法implicit noneinteger i,j,kcharacter file_name*10open(10,file='file_name.txt',status='old')do k=1,3read(10,*)file_nameopen(k,file=file_name,status='unknown') if(k.eq.1) then write(k,*) 'one'else if(k.eq.2) thenwrite(k,*) 'two'else if(k.eq.3) thenwrite(k,*) 'third'endifenddoclose(k)close(10)end programHw3.fimplicit noneinteger jreal k,dtcomplex*16 ya(21)complex*16 yf(21)complex*16 yb(21)complex*16 ireal rya(21),ryf(21),ryb(21)real iya(21),iyf(21),iyb(21)i=cmplx(0.,1.)open(10,file='hw3.dat',form='unformatted',access='direct' ,r ecl=4*21) !21代表畫圖時x軸從時間0到20 do j=1,21 ya(j)=exp(i*0.5*(j-1.))enddodo j=1,21rya(j)=dreal(ya(j)) !取複數實部 enddodo j=1,21iya(j)=dimag(ya(j))enddoyf(1)=1. !ICyf(2)=1+0.5*i !ICdo j=2,20yf(j+1)=yf(j-1)+i*yf(j)enddodo j=1,21ryf(j)=dreal(yf(j))enddodo j=1,21iyf(j)=dimag(yf(j))enddoyb(1)=1.yb(2)=1+0.5*ido j=2,21yb(j+1)=yb(j)+0.25*i*(3.*yb(j)-yb(j-1))enddodo j=1,21ryb(j)=dreal(yb(j))enddodo j=1,21iyb(j)=dimag(yb(j))enddowrite(10,rec=1) (rya(j),j=1,21) !不能用do迴圈寫write(10,rec=2) (ryf(j),j=1,21) !用do寫跑出來的圖 write(10,rec=3) (ryb(j),j=1,21) !會怪怪的write(10,rec=4) (iya(j),j=1,21)write(10,rec=5) (iyf(j),j=1,21)write(10,rec=6) (iyb(j),j=1,21)close(10)end programhw3.ctldset ^hw3.datundef -999.9title hw3xdef 21 linear 0 1 # 0:x軸的起始點,間隔為1 ydef 1 linear 0 1tdef 1 linear 00z01nov2010 12hrzdef 1 linear 1 1vars 6rya 0 99 qryf 0 99 wryb 0 99 eiya 0 99 riyf 0 99 tiyb 0 99 yendvarshw3.gs'reinit''open hw3.ctl''set button bcol1' #'set xlab auto''set grads off'#'set xaxis 0 20' #會造成座標軸重複覆蓋#'set yaxis -2 2 1' #'d ryb'#'d ryf'#'d rya'#'draw title real''d iyb''d iyf''d iya''draw title imagine''draw xlab time''draw ylab pi(t)'1000.gs'reinit''open 198706_part_all.ctl''set grads off' # 把圖下方的grads去除'set lon 90 150' # 取經度90~150'set lat 5 45' # 取緯度5~45'set mpdset hires''set lev 500' # 500hPa高度'set t 1 30' # 所取時間為1~30tt=1while(tt<=30) # t在30以下都可畫出'c' # 防止在每個t時間所畫出的圖被一直蓋上'set t 'tt # t=上面的tt'set gxout stream''set ccolor 2''d uprs*5;vprs*5''set gxout contour''d zprs/9.8'pull aaqtt=tt+1 # 按下任意鍵後再畫出下一張圖endwhile08+09.gs#average (08+09)/2 height'reinit''open 200806_part_all.ctl' # 自動等於dfile 1 'open 200906_part_all.ctl' # 自動等於dfile 2#'set mproj nps''set grads off''set grid off''set lon 90 150''set lat 5 45''set mpdset hires'zz=1while(zz<=10) 'clear' # 清除上一張圖'set z 'zz'set gxout contour''set dfile 1''define avea=ave(zprs,t=1,t=60)' # avea是自行取的任意變數'set dfile 2''define aveb=ave(zprs,t=1,t=60)''define avez=(avea+aveb)/2''d avez/9.8''draw xlab longitude''draw ylab latitude''draw title (08+09)/2'pull abczz=zz+1endwhile19870623.gs'reinit''open 198706_part_all.ctl'#'enable print test'#============================= 'set grads off''set lon 115 127''set lat 20 28''set lev 850''set t 45''set mpdset hires'#============================= #'set gxout shaded'#============================= #'set ccolor 4''set cstyle 7''d tprs'#============================= 'set gxout contour''d zprs/9.8'#============================= #'set ccolor 6''set gxout barb''d uprs*2;vprs*2''set string 1'#'set strsiz 2 2''draw title 850hPa 1987-06-23-00:00'1.gs'reinit''open 198706_part_all.ctl''set grads off''set grid off''set lon 90 150''set lat 5 45''set mpdset hires' zz=1while(zz<=10)'c''set z 'zz'set ccolor 5''set cstyle 7''define avet=ave(tprs,t=1,t=60)' 'd avet'#'set gxout contour'#'define avez=ave(zprs,t=1,t=60)' #'d avez/9.8'*'set gxout barb'*'set ccolor 1'*'define aveu=ave(uprs,t=1,t=60)' *'define avev=ave(vprs,t=1,t=60)' *'d aveu*5;avev*5'pull qqqif(z=1)'draw title 123';endifif(z=2)'draw title 123';endifzz=zz+1endwhile2_years_average.gs'reinit''open 198706_part_all.ctl''open 200806_part_all.ctl''set grads off''set lon 105 130''set lat 10 40''set mpdset hires''set dfile 1''set dfile 2'zz=1while(zz<=9)'c''set gxout contour''define avez=ave(zprs.1,t=1,t=60)' #兩年六月份平均'define avea=ave(zprs.2,t=1,t=60)''set z 'zz'd (avea+avez)/9.8'pull abczz=zz+1endwhileaverage.gs'reinit''open 198706_part_all.ctl''set grads off' 'set grid off''set lon 90 150''set lat 5 40''set mpdset hires'zz=1while(zz<=9)'c''set grads off''set gxout barb''define aveu=ave(uprs,t=1,t=60)' 'define avev=ave(vprs,t=1,t=60)' 'd aveu*4;avev*4''set gxout contour''set z 'zz'define avez=ave(zprs,t=1,t=60)' 'd avez/9.8'*'set cstyle 2'*'define avet=ave(tprs,t=1,t=60)' *'d avet'*no sensepull abczz=zz+1endwhilewhile z=1,9if(z=1)'draw title 1000hPa';endif if(z=2)'draw title 850hPa';endif endwhilediscrete.gs'reinit''open 198706_part_all.ctl''set lon 100 140''set lat 15 40''set mpdset hires''set t 1'zz=1while(zz<=9)'c''set grads off'*'set gxout barb'*'d uprs*4;vprs*4''set gxout contour''set z 'zz'd zprs/9.8''set cstyle 2' 'set ccolor 2''d tprs'pull abczz=zz+1endwhilewhile(zz<=9)if(z=1);fmon='1000hpa';endifif(z=2);fmon='850hpa';endifendwhilef1987.gs'reinit''open 198706_part_all.ctl''set fwrite f1987.dat' 開起.dat 'set gxout fwrite' 將資料寫到.dat 'set lon 120''set lat 25''set t 1 60''set z 1 14'tt=1while(tt<=60)'c''set t 'ttzz=1while(zz<=14)'c''set z 'zz'd zprs' # 重力位高度,需除9.8'd tprs''d uprs''d vprs''d rh'zz=zz+1endwhilett=tt+1endwhileimplicit noneinteger k,t,zrealzprs(14,60),tprs(14,60),uprs(14,60),vprs(14,60),rh(14,60) open(10,file='f1987.dat',form='unformatted',status='old' ,access='direct',recl=4)open(11,file='f1987.txt')k=0do t=1,60 #和f1987.gs檔時間順序一樣do z=1,14read(10,rec=5*k+1) zprs(z,t) #注意rec寫法,不是從10開始read(10,rec=5*k+2) tprs(z,t)read(10,rec=5*k+3) uprs(z,t)read(10,rec=5*k+4) vprs(z,t)read(10,rec=5*k+5) rh(z,t)k=k+1enddo !zenddo !tdo t=1,60do z=1,14write(11,'(f10.2)') zprs(z,t)/9.8 write(11,'(f10.2)') tprs(z,t)write(11,'(f10.2)') uprs(z,t)write(11,'(f10.2)') vprs(z,t)write(11,'(f10.2)') rh(z,t)enddo !zenddo !tclose(10)close(11)end program畫台灣地形 ter2.gs'reinit''open ter2.ctl'#============================='set grads off''set lon 120.5 122''set lat 21.5 24''set mpdset hires''set xlint 0.5' #x軸數字間隔'set ylint 0.5' #y軸數字間隔'set strsiz 5 8'#'set ccolor 3'#============================='set gxout shaded'#'set cint 1013'#'set cmax 1''set cmin 1' #第一層以下無顏色'set cstyle 4''d ter'#'set gxout contour'#'d ter''set display color white''printim ter2.gif gif'2009-10.fimplicit noneinteger i,j,k,a,b,c,d !a.d=總資料數parameter(a=285936,b=384,c=930,d=30) !b.c=測站數character (20) station_auto(a),station_cwb(c)character (11) x1,x2,x3,x4,x5,x6,x7,x8,x9,x10character (11) y1,y2,y3,y4,y5,y6,y7,y8,y9,y10character (11) z1,z2,z3,z4,z5,z6,z7,z8,z9,z10character (11) w1,w2,w3,w4,w5,w6,w7,w8,w9,w10 character (11) time,ti(a),time_cwb,ti_cwb(c)real total,total_cwb,mr_cwb(c)real p,t,xx,xxx,rain,xxxxreal pp(a),tt(a),yy(a),yyy(a),rr(a),yyyy(a),dr(a)real rain_cwbreal rr_cwb(c)open(10,file='autop2009-10.txt')open(12,file='d-cwb200910.txt')open(11,file='2009-10.txt')do i=1,aread(10,*) station_auto(i),time,p,t,xx,xxx,rain,xxxxti(i)=timerr(i)=rainenddodo i=1,c !read:讀一整列(.txt橫的資料)read(12,*)station_cwb(i),time_cwb,x3,x4,x5,x6,x7,x8,x9,x10, !土法煉鋼$ y1,y2,y3,y4,y5,y6,y7,y8,y9,y10,$ z1,z2,z3,z4,z5,z6,z7,z8,z9,z10,$ w1,w2,w3,w4,w5,w6,w7,w8,w9,w10,rain_cwbti_cwb(i)=time_cwbrr_cwb(i)=rain_cwbenddoc=================================================auto station do k=1,btotal=0.do i=1.+(k-1.)*744.,744.*kif(rr(i).eq.-9996.)thenrr(i)=0.else if(rr(i).eq.-9991.)thenrr(i)=0.else if(rr(i).eq.-9997.)thenrr(i)=0.endiftotal=total+rr(i)dr(i)=totalenddoenddodo i=744,a,744 !字串需用自由格式*c write(11,'(45x,f7.1)')dr(i) !45x:空45格write(11,*)station_auto(i),ti(i),i,dr(i)enddoc====================================== ============cwb station do k=1,dtotal_cwb=0.do i=1.+(k-1.)*31.,31.*kif(rr_cwb(i).eq.-9998.)thenrr_cwb(i)=0.endiftotal_cwb=total_cwb+rr_cwb(i)mr_cwb(i)=total_cwbenddoenddodo i=31,c,31c write(11,'(45x,f7.1)') station_cwb(i),mr_cwb(i)write(11,*) station_cwb(i),ti_cwb(i),i,mr_cwb(i)enddoc====================================== =================== close(10)close(11)close(12)end programa.outC0A520 2009103124 744 199.C0A530 2009103124 1488 813.C0A540 2009103124 2232 782.C0A550 2009103124 2976 1541.C0A560 2009103124 3720 860.5C0A570 2009103124 4464 746.C0A580 2009103124 5208 512.5C0A590 2009103124 5952 806.5C0A860 2009103124 6696 1052.5C0A870 2009103124 7440 405.5C0A880 2009103124 8184 561.C0A890 2009103124 8928 898.C0A920 2009103124 9672 443.5C0A930 2009103124 10416 502.1011rain.fimplicit noneinteger i,j,k,a,b,ccharacter station(550)*6,st*6 !數字超過測站即可integer mm,dd,hh,st_dd(10000),st_hh(10000),st_mm(10000) real rain,lon,lat,st_lon(10000),st_lat(10000550)real st_rain(10000)open(11,file='1011rain.txt',status='old')open(10,file='sister.txt',status='old')c====================================== ======================= c stc COA520i=0doc i=i+1 !這個位置會出Segmentation faultread(10,*,end=31)st,lon,lat,mm,dd,hh,rainif(mm.eq.10.and.dd.eq.11.and.hh.eq.24)then !必寫,要不然會亂讀時間i=i+1 !這個位置OKstation(i)=' 'c if(mm.eq.10.and.dd.eq.11.and.hh.eq.24)thendo j=1,iif(st.eq.station(j).and.lon.eq.st_lon(j)/doc/0 914886775.html,t.eq.st_lat(j $).and.mm.eq.st_mm(j))then !測站、經緯度都重複,需用矩陣 i=i-1goto 32endifenddostation(i)=stc print*, station(i)st_lon(i)=lonst_lat(i)=latst_mm(i)=mmst_dd(i)=ddst_hh(i)=hhst_rain(i)=rain write(*,21)station(i),st_lon(i),st_lat(i),$ st_mm(i),st_dd(i),st_hh(i) !都用矩陣寫才可c write(*,21)station(i),st_lon(i),st_lat(i),c $ st_mm(i),st_dd(i),st_hh(i)32 continuecprint*,station(i),st_lon(i),st_lat(i),st_mm(i),st_dd(i),st_hh( i)enddocprint*,station(i),st_lon(i),st_lat(i),st_mm(i),st_dd(i),st_hh(i)31 continuek=iprint*,kc====================================== ======================= i=0doi=i+1read(10,21,end=41)station(i),st_lon(i),st_lat(i),mm,st_dd(i) $ ,st_hh(i),rainif(st_dd(i).eq.10.and.st_hh(i).eq.24)thenc print*,st_hh(i)cwrite(11,21)station(i),st_lon(i),st_lat(i),mm,dd,hh,rainendifenddo41 continue21 format(a6,1x,f7.3,1x,f6.3,1x,i2,1x,i2,1x,i2,1x,f7.2) close(11) ! 1011rain.txtclose(10) ! ister.txtend programC1V450 120.376 22.540 10 11 24 0.00C1V460 120.884 23.438 10 11 24 6.00C1X010 120.252 23.332 10 11 24 0.00在地圖上標記號在grads下1. reinit2. open ter2.ctl3. set lon和lat4. d ter5. q w2xy 經度緯度6. X = ? Y = ?2008_2009.gs'reinit''open 200806_part_all.ctl''open 200906_part_all.ctl''set mpdset mres''set grads off''set lon 100 140' 'set lat 7 37''set lev 1000''set xlint 5''set xlint 5''set gxout contour''set dfile 1'#'define tprs.1=avg(tars,t=1,t=60)' #可不用設這個'set t 10 ' #有設這個時'define ctone=tprs.1-273' #變數用全英文比較好'set gxout contour''set dfile 2''set t 10 ''define cttwo=tprs.2-273' 'd (ctone+cttwo)/2'。
VB与FORTRAN、GrADS混合编程开发绘制降水分布图软件
VB与FORTRAN、GrADS混合编程开发绘制降水分布图软件张富龙;刘爽;兰明胜【摘要】通过研究VB、FORTRAN和GrADS三者之间的相互调用方法,利用VB 编程指令代码开发操作界面,FORTRAN进行数据处理转换,GrADS绘制图形。
在三者有效结合的编程技巧下,实现了绘制乡镇加密自动站降水分布图软件的开发。
【期刊名称】《气象灾害防御》【年(卷),期】2015(022)001【总页数】3页(P35-37)【关键词】GrADS 降水分布图乡镇加密自动站混合编程【作者】张富龙;刘爽;兰明胜【作者单位】松原市气象台,松原138000【正文语种】中文【中图分类】P409绘制降水分布图软件有很多,但是以往绘制往往是每次手动输入指令,调试数据,比较麻烦,而且浪费时间,即便有专门的降水分布图绘制软件,也基于省级地区以上的边界为底图,不适合市、县局的业务应用。
现根据实际业务需要,为了更方便地做好降水服务工作,特利用VB、FORTRAN和GrADS三者的混合编程,开发一款以市、县边界为底图的,能够自动处理降水数据、可操作性强的降水分布图软件,用于日常业务工作,使日常工作流程化、规范化和具有可操作性,实现雨情分布任意时段一键出图的功能,方便数据查询和服务材料的制作。
FORTRAN语言是世界上广泛流行的、最适用于数值计算的一种计算机语言,具有强大的数值计算功能与数学分析能力,长期以来在气象领域做出了重大的贡献[1]。
但其在可视化程序设计方面比较欠缺。
GrADS(Grid Analysis and Display System,格点分析和显示系统)是当今气象界广泛使用的一种数据处理和显示软件系统,其提供了一个全32位交互操作的气象格点数据与站点数据的分析与显示环境,再加上该软件具有操作简单、功能强大、显示快速、出图类型多样化、图形美观等特点,使其一直以来备受气象同行的青睐[2]。
但其在数据处理方面比较差,所以FORTRAN和GrADS经常一起在气象中使用互补不足[3-4]。
优秀的fortran程序编程规范
Programming Guidelines for PARAMESH Software Development(NOTE: This document is heavily based upon theIntroductionThis document describes the programming guidelines to be used by software developers wishing to contribute software to the PARAMESH, parallel, adaptive mesh refinement software. We welcome people to contribute software and/or bug fixes to the PARAMESH AMR software. Software to be added to PARAMESH can come in 2 forms:∙Improvements to the basic PARAMESH kernal software found in the mpi_source, source and hearders directories.∙Software the addes additional functionality to PARAMESH. This type of software should be added as separate entities within the utilities directory.Complete applications should not be added as part of PARAMESH. PARAMESH is only meant to be a tool which supports parallel adaptive mesh applications and any software which supports this goal will be considered for acceptance into PARAMESH. For instance, a solver for the poisson equation that works with PARAMESH would be acceptable, but an application that solves the equation of gas dynamics would not.The PARAMESH software is slowly being evolved to be consistent with this document. Any new software which is contributed should follow these guidlines. If not, it will be rejected. This document deals mainly with Fortran 90, since most new PARAMESH software will probably be written in that language. [Throughout this document, the term "Fortran" should be understood to mean Fortran 90.] Since we expect C and C++ also to be used, a separate document dealing with them will be developed in thefuture. In the meantime, this document can serve as a general guideline for developing code to be used with PARAMESH in those programming languages.The guidlines in this document should be adhered to by ANY software which will be released as part of the PARAMESH package ofsource code. This includes software 'utilities' (stored in the paramesh/utilities directory) which add functionality to PARAMESH for different algorithms. It also should be applied to any new code developed and added to the main source code for PARAMESH in the paramesh/source, paramesh/mpi_source, or paramesh/headers directories.The guidelines are intended to enhance the following aspects of the final product, listed in decreasing order of importance:∙Maintainability- refers to how easy it is to understand the purpose of each element of the program, and to modify and extend the program.∙Portability- refers to how easily the program can be ported to new computational platforms.∙Efficiency- refers to the amount of computer resources (CPU time, memory, disk storage, etc.) required to run the program.Program Development and DesignItems in this section are fairly general and fundamental in nature. They impact all three of the items listed above - maintainability, portability, and efficiency.LanguageTry to use ANSI standard Fortran 90 exclusively. If you must, you can use C or C++, but it must work with PARAMESH and be callable from a Fortran 90 program.Organization∙Write modular code.∙In general, put each subprogram in a separate file, using the subprogram name as the file name, with a .F90 extension.∙Within each routine, use interface blocks to explicitly specify the interface to your contributed routines.∙Group related files in a single directory.∙Names of files and directories should reflect their purpose. Common Blocks∙Don't use common blocks, use Modules instead. Period !Data Types∙Use Implicit none in each program unit, and explicitly declare all variables and parameters. Common variables and parameters should be declared in the relevant include file.∙Don't use *'ed forms, like Real*8. Declare variables using Real ::∙Don't compare arithmetic expressions of different types; convert the type explicitly.Dynamic Memory∙Assign memory for arrays dynamically, using automatic arrays, allocatable arrays, and/or array pointers. Explicitly deallocate memory used by allocatable arrays and array pointers when they're no longer needed.Coding Style (see mpi_source/mpi_amr_guardcell.F90 for a complete example)Items in this section are fairly specific, and primarily impact the readability, and thus the maintainability, of the final product. It is recognized that rules for "good coding style" are somewhat subjective. Program Units∙Begin main programs with a Program statement.∙Don't use multiple entries or alternate returns.∙Use the intent attribute in the type declaration statement for all variables passed into our out of subroutines and functions. Make sure to include these in the interface block that you create for the surbroutine.∙Match the arguments in the calling (sub)program to those of the called subprogram in both number and type.∙Use the following order for statements within each subprogram: o Standard header section. This should be comments in the format used with the robodoc code documenation software (Seethe PARAMESH source code for examples).o Use moduleso Parameter definitionso Type declarations for subprogram argumentso Type declarations for local variableso Executable code∙Functions should not have side effects. (I.e., don't change the arguments or any common variables inside the function.) ∙Use generic names for library functions, rather thanprecision-specific ones.∙Name external functions in an External statement.Statement Form∙Use free-form formatting, but for readability:o Keep line lengths below 80 characters.o Start each line in column 7 or higher.o Reserve columns 1-5 for statement labels.o Don't use the optional continuation character (i.e., &) at the start of continuing lines. Avoid splitting keywords andcharacter strings between lines.Note that with free-form formatting, an & must be the last character (except for comments) in a line that is to be continued.∙Use a ! in column 1 for non-blank comment lines.∙Split long lines before or after an operator, preferably a + or -.∙Don't write more than one statement per line.Statement Labels∙Minimize the use of statement labels, where appropriate.∙Don't use unreferenced labels.∙Use labels in ascending order.Upper/Lower Case∙Use upper case for parameters, upper case for subrotines and functions from libraries outside of PARMAESH (such as MPI), lower case with an initial capital letter for Fortran keywords, and lower case for everything else except comments and character strings.∙Write comments as normal text, with normal capitalization rules. Spacing∙Use spacing to enhance readability.∙Indent contents of code blocks (i.e., do loops, block if, etc.).Suggested amount is three spaces.∙Don't use tabs.∙Use spacing in equations to clarify precedence of operators. I.e., normally put one space on either side of =, +, and - operators(except in subscripts), but none around *, /, or ** operators. For example, this:y1 = (-b + Sqrt(b**2 - 4.*a*c))/(2.*a)is easier to read than this:y1=(-b+Sqrt(b**2-4.*a*c))/(2.*a)or this:y1 = ( - b + Sqrt ( b ** 2 - 4. * a * c ) ) / ( 2. * a ) ∙Use spacing to reveal patterns in continuation lines and in separate but logically related statements. For example, this:dum1 = Sqrt((fr (i,j) - fr ( 1, j))**2 + &(fth(i,j) - fth( 1, j))**2)dum2 = Sqrt((fr (i,j) - fr (n1, j))**2 + &(fth(i,j) - fth(n1, j))**2)dum3 = Sqrt((fr (i,j) - fr ( i, 1))**2 + &(fth(i,j) - fth( i, 1))**2)dum4 = Sqrt((fr (i,j) - fr ( i,n2))**2 + &(fth(i,j) - fth( i,n2))**2)is easier to read than this:dum1 = Sqrt((fr(i,j) - fr(1,j))**2 + (fth(i,j) - &fth(1,j))**2)dum2 = Sqrt((fr(i,j) - fr(n1,j))**2 + (fth(i,j) - &fth(n1,j))**2)dum3 = Sqrt((fr(i,j) - fr(i,1))**2 + (fth(i,j) - &fth(i,1))**2)dum4 = Sqrt((fr(i,j) - fr(i,n2))**2 + (fth(i,j) - &fth(i,n2))**2)Variable Names∙Use names that are descriptive of the entity being represented, and/or are consistent with the standard notation in the field.∙In general, follow standard Fortran convention for the variable type. I.e., integers start with i, j, k, l, m, or n, all others are real.∙Don't use keyword, subprogram, or module names for variables.∙Don't give a local variable the same name as any common variable. Arrays∙Dimension arrays in the type declaration statement, not in a separate Dimension statement.∙When passing character variables into a subprogram, use the assumed-length form in the type declaration statement inside the subprogram. I.e.,Subroutine sub (c)Character*(*) c∙Don't exceed the bounds of the array dimensions.Control Statements∙Short do loops may be written using simple Do and End do statements, without labels.∙Long do loops and if blocks (more than a page or so), should mark the end of the construct in some way that "connects" it with the start. One convenient and readable method is to use an in-linecomment on the ending statement that repeats the beginningstatement. E.g.,If ( bccode == 13 ) then[Lines and lines of code]Do i = 1,nzonesIf ( zondim(1,i) > 0 ) then[More lines and lines of code]End if ! If ( zondim(1,i) > 0 ) thenEnd do ! Do i = 1,nzonesEnd if ! If ( bccode == 13 ) then∙Minimize the use of Go to statements, especially where they can be replaced by short'ish if blocks, but don't create convoluted code just to avoid using them. Don't be afraid to use a Go to where it makes sense. An example might be a long (more than a page)conditional section of code. In this case a well-commented Go to block, which ends with an easily-noticed statement label, may be more readable than an indented if block without an ending statement label. Also consider making a long conditional section a separate subprogram.Calls to other Libraries (such as MPI).<>Capitalize the entire subroutine name when making the call to the libary routine, e.g.Call MPI_BARRIER(MPI_COMM_WORLD,ierr)Comments∙Use comments liberally to describe what's being done. Where code may be confusing, use longer comments to describe why something's being done the way that it is.∙Make each comment meaningful; don't simply re-iterate what's already obvious from the coding itself. As an obvious example, this:!-----Fill the guardcells of all PARAMESH blocksCall amr_guardcellis more meaningful than this:!-----Call amr_guardcellCall amr_guardcell∙Use a consistent method to help the reader distinguish comments from code, such as the "-----" leaders in the examples above.∙Start the text of comments at the same indentation level as the code being described.∙Use a standard header section at the beginning of each subprogram defining its purpose.∙Use in-line comments, with ! as the delimiter, where appropriate for short explanations or clarifications. Start in-line comments far enough to the right (e.g., three spaces or more from the end of the statement) to help distinguish comments from code. Whereappropriate, align them vertically with nearby in-line comments.∙Define each common block variable using an in-line comment on its type statement in the include file. Each common variable will thus have a separate type statement.∙Define key local variables using in-line comments on the type statements in the subprogram.Obsolete/Forbidden FeaturesThe following Fortran features are either formally declared as obsolete, or widely considered to be poor programming practice, and should not be used:∙Arithmetic if statements∙Do loops with non-integer indices∙Shared do loop termination statements∙Pause statements∙Assigned and computed Go to statements∙Hollerith edit descriptors and Hollerith character strings∙Equivalence statements∙Alternate return statements。
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年代开发,并且至今仍广泛应用于科学和工程领域。
Grads站点资料的多层次绘图
Grads站点资料的多层次绘图看了很多关于如何使用grads和fortran转换站点资料的帖子,大多数都集中在地面或者单层的,对于高空的多个层次很少有提到,而且大多数有些模糊,不好理解,我把自己做的一些整理传上来,供大家参考参考。
欢迎大家指出理解和程序的错误。
一,首先是数据的说明选用micaps高空观测资料850hpa,700hpa,600 hpa,500hpa,200hpa,100 hpa 六个层次作为例子。
为了方便,时间层次就选择1个时间。
micaps高空观测资料的说明如下:diamond 2 10年07月28日20时200百帕高空观测10 07 28 20 200 419因此定义mt=419,nz=4数据每行的排列分别是:站号,经度,纬度,高度……..具体详见micaps 的使用手册,这里不再赘述。
根据这些定义变量。
二,fortran读取站点资料注:1.根据micaps的高空资料定义相应的变量,由于850-100这6个层次的站点是相同的,站号,经纬度等只需定义成一维的,如lat(mt),而lev(nz)表示有6个层次,也为一维;高度场等又有mt个站的资料又有nz个层次,因此定义为2维的。
2.由于micaps资料分别在不同的文件内,文件名重复,我修改为100_2820,200_2820,500_2820…….然后使用trim函数顺序读取资料。
3.为了更好的说明地面和高空资料的读取,将地形高度high(mt)设置为地面的量,nlev定义为7,即高空六层,地面一层4.以输出地面一个变量,高空三个变量为例。
先进行站点的循环,再进行高空层次的循环。
5.前面我们定义了lev(nz),现在给lev(nz)赋值,使用select case.,分别赋值为850,700,600,500,200,100。
在使用是发现如果不进行赋值,后续工作也无法进行,会出现错误。
程序如下:program stationparameter (mt=419,nz=6) ! mt 总站数,nz总层次character stid(mt)*8integer flag,nlevreal lat(mt),lon(mt),tim,high(mt),number(mt),lev(nz) realheight(mt,nz),temp(mt,nz),ttd(mt,nz),ud(mt,nz),us(mt,nz),u(mt,nz),v(m t,nz),td(mt,nz)!顺序读取资料character stationfile(6)*8data stationfile/'850','700','600','500','200','100'/do k=1,nzopen(k,file=''//trim(stationfile(k))//'_2820.000')read(k,*)read(k,*)!跳过micaps开头两行的说明文字do i=1,mtread(k,*)stid(i),lon(i),lat(i),high(i),number(i),height(i,k),temp(i,k ),ttd(i,k),ud(i,k),us(i,k)end doend do !结束k的循环tim=0.0flag=1nlev=7 !高空+地面open(40,file='20.grd',form='binary')do i=1,mt !站点循环write(40)stid(i),lat(i),lon(i),tim,nlev,flag,high(i)do k=1,nz !高空循环!对lev(nz)进行赋值select case(k)case(1)lev(k)=850.0case(2)lev(k)=700.0case(3)lev(k)=600.0case(4)lev(k)=500.0case(5)lev(k)=200.0case defaultlev(k)=100.0end select!在lev(k)后将要输出的所有高空变量一次性输出write(40)lev(k),height(i,k),temp(i,k),ttd(i,k)end doend donlev=0write(40)stid(mt-1),lat(mt-1),lon(mt-1),tim,nlev,flag close(40)end三,ctl及map文件的生成注:看到相关帖子说在站点资料的ctl文件中也要设置zdef,我试了下,设置和不设置都是可以的Ctl文件如下:dset f:\20.grddtype stationstnmap f:\20.mapundef 9999.0title 20zdef 6 levels 850 700 600 500 200 100tdef 1 linear 28jul2010 1DYvars 4high 0 99 Terrain Heighth 6 99 heightt 6 99 temperaturettd 6 99 ttdendvarsCtl文件编写完成后,就可以在grads中生成map文件了,相关命令如下:! stnmap -i f:/20.ctl注意,斜杠一定是“/”,方向不能反第三步结束后,站点资料的读取就结束了,可以通过grads打开Ctl文件看到相关站点的值了。
fortran实习4grades
实习四数据转换及数据描述文件建立实习目的:熟悉GrADS数据格式,掌握数据描述文件的建立。
实习内容:现有ASCII码(十进制数据格式)数据资料文件u850.dat和v850.dat,其空间范围:60-150°E,0-40°N;层次:u、v为850Pa;时段:1982.1-1985.12;分辨率:2.5°*2.5°;数据排放顺序满足GrADS要求。
请编写Fortran程序将这2个文件转换成1个二进制(binary)文件,并给出相应数据描述文件(即CTL文件)。
实习要求:根据实习内容和资料,撰写Fortran程序,并运行检查。
根据数据文件描述信息,打开文本编辑器,编写CTL文件。
打开GrADS,参考GrADS常用命令,编写gs文件绘制水平风场图,并存图。
一.程序Fortran:program shixi3implicit noneinteger,parameter::x=37,y=17,tim=48real v(x,y,tim),u(x,y,tim)integer i,j,kopen(1,file='c:\u850.dat')open(2,file='c:\v850.dat')do k=1,timread(1,100) !为跳过年份和时间行read(1,200)((u(i,j,k),i=1,x),j=1,y) !读入u方向各经纬度风速read(2,100)read(2,200)((v(i,j,k),i=1,x),j=1,y)enddo100 format(2I7)200 format(37F6.2)open (3,file='c:\wind850.grd',form='binary')do k=1,timwrite(3) ((u(i,j,k),i=1,x),j=1,y) !写入文件write(3) ((v(i,j,k),i=1,x),j=1,y)enddoclose(1)close(2)close(3)end program shixi3Ctl文件:dset C:/wind850.grdundef -9.99E+33title horizontal wind fild xdef 37 linear 60.0 2.5 ydef 17 linear 0.0 2.5zdef 1 levels 850tdef 48 linear JAN1982 1mo vars 2u 1 99 u wind(m/s)v 1 99 v wind(m/s) endvars二.结果截图。
Fortran语言快速学习指南
Fortran语言快速学习指南Fortran是一种高级程序设计语言,专门用于科学和工程计算。
它是最早的计算机编程语言之一,在科学计算领域一直具有广泛的应用。
本篇文章将提供一份Fortran语言的快速学习指南,帮助读者快速掌握该语言的基本概念和用法。
一、Fortran简介Fortran是一种编译型语言,最初于1957年由IBM公司开发。
它的名字源自“Formula Translation”的缩写,最早用于数值计算和科学研究。
Fortran的设计目标是简化数值计算的编程过程,提高程序的效率和可读性。
二、Fortran的基本语法1. 变量和数据类型在Fortran中,变量是用于存储数据的容器,每个变量都必须先声明后使用。
Fortran提供了多种数据类型,包括整型、实型、字符型等,用于存储不同类型的数据。
2. 数组Fortran中的数组是一种特殊的变量类型,可以同时存储多个值。
数组的声明方式为:TYPE, DIMENSION(下标范围) :: 数组名。
例如,声明一个包含10个整数的数组可以写作:INTEGER, DIMENSION(10) :: array。
3. 条件语句和循环Fortran提供了条件语句(IF语句)和循环结构(DO循环)来控制程序的执行流程。
条件语句根据某个条件的真假来选择不同的执行路径,而循环结构可以重复执行一段代码块多次。
4. 子程序和函数Fortran支持子程序和函数的定义和调用。
子程序用于执行一段特定的代码块,而函数则返回一个值给调用者。
子程序和函数的定义方式为:SUBROUTINE 子程序名和 FUNCTION 函数名。
三、Fortran的特性和优势1. 高效性:Fortran在科学计算领域具有很高的效率和性能优势。
它能够直接利用计算机硬件的特性,进行高速、大规模的数值计算。
2. 数值计算支持:Fortran提供了丰富的数值计算库和函数,方便进行复杂的科学计算。
这些库包括矩阵运算、线性代数、微分方程求解等。
fortran常用算法程序集
Fortran常用算法程序集简介Fortran是一种面向科学和工程计算的编程语言,通常用于数值计算和数据分析。
它有着强大的数学计算能力和高性能,被广泛应用于科学计算、工程仿真、天气预报等领域。
本文将介绍一些常用的Fortran算法程序集,包括数值积分、矩阵运算、排序算法等。
数值积分数值积分是求解定积分的一种方法,用于计算曲线下面积、求解微分方程等。
Fortran提供了一些常用的数值积分算法,如梯形法则、辛普森法则等。
梯形法则梯形法则是数值积分中最简单的算法之一,基本思想是将曲线下面积近似为一系列梯形的和。
下面是使用Fortran编写的梯形法则算法示例:! 梯形法则real function trapezoidal_rule(f, a, b, n)real, external:: freal:: a, b, n, hreal:: x, suminteger:: ih = (b - a) / nsum= f(a) + f(b)do i =1, n-1x = a + i * hsum=sum+2* f(x)end dotrapezoidal_rule = (h/2) *sumend function trapezoidal_rule辛普森法则辛普森法则是一种更精确的数值积分算法,基于多项式插值的思想。
它将曲线分成若干小段,每段近似为一个二次函数,然后对每个二次函数进行积分。
下面是使用Fortran编写的辛普森法则算法示例:! 辛普森法则real function simpsons_rule(f, a, b, n)real, external:: freal:: a, b, n, hreal:: x, sum1, sum2integer:: ih = (b - a) / nsum1 = f(a) + f(b)sum2 =0do i =1, n-1, 2x = a + i * hsum2 = sum2 +4* f(x)end dodo i =2, n-2, 2x = a + i * hsum2 = sum2 +2* f(x)end dosimpsons_rule = (h/3) * (sum1 + sum2)end function simpsons_rule矩阵运算矩阵运算是科学计算中常用的一个重要环节,Fortran提供了丰富的矩阵运算库,包括矩阵乘法、矩阵转置、矩阵求逆等。
fortran语言的基本概念
fortran语言的基本概念Fortran语言的基本概念概述Fortran(Formula Translation)是一种用于科学计算和数值分析的高级编程语言。
它是最早的编程语言之一,由IBM公司在20世纪50年代开发。
Fortran语言以其高效、可靠和快速的特点,长期被广泛应用于科学和工程计算领域。
基本结构Fortran程序以程序单元作为基本的执行单位,程序单元分为主程序和子程序两种类型。
主程序•主程序是Fortran程序的入口点,每个Fortran程序必须包含一个主程序。
•主程序由PROGRAM关键字定义,后接程序名称。
•主程序内包含了程序的执行逻辑和控制流程,通过调用子程序完成具体的计算任务。
•子程序是主程序的辅助部分,用于定义重复使用的计算任务或功能模块。
•子程序由SUBROUTINE或FUNCTION关键字定义。
•SUBROUTINE用于定义过程(Procedure),只执行特定的计算任务,不返回结果。
•FUNCTION也是定义过程,但它还可以返回一个值作为结果。
变量与数据类型Fortran语言使用变量来存储和操作数据,变量需要先进行声明。
Fortran提供了多种数据类型,用于存储不同类型的数据。
声明变量•变量声明使用DIMENSION或INTEGER/REAL/COMPLEX/LOGICAL等关键字。
•变量声明语句通常出现在主程序的开头或子程序的参数列表中。
数据类型•INTEGER用于表示整数类型的数据。
•REAL用于表示浮点数类型的数据。
•COMPLEX用于表示复数类型的数据。
•LOGICAL用于表示逻辑类型的数据,可取值为TRUE或FALSE。
Fortran语言提供了多种控制流程语句,用于在程序中实现条件判断和循环操作。
条件判断•IF-THEN语句用于根据条件执行不同的代码块。
•IF-THEN-ELSE语句可根据条件选择执行不同的代码块。
循环•DO语句用于实现循环操作,可指定循环的起始值、结束值和步进值。
fortran 教程
fortran 教程Fortran是一种古老而强大的编程语言,最初在1957年开发。
它被广泛用于科学和工程计算,特别是对大型和复杂的计算任务。
Fortran之所以如此受欢迎,是因为它在数学计算领域表现出色。
它拥有丰富的数学函数和运算符,并且支持高精度计算。
此外,Fortran还具有强大的数组处理能力,可以轻松处理大规模数据。
Fortran的语法相对简单,易于学习和使用。
它使用英语类似的语法,语句以换行符结束。
Fortran中的语句通常以关键字开始,例如"PROGRAM","SUBROUTINE"和"DO"等。
Fortran具有自己的变量类型,包括整数(INTEGER)、实数(REAL)和字符(CHARACTER)等。
变量必须在使用之前先声明,并且可以指定其大小和精度。
Fortran还支持过程式编程,包括子程序和函数的定义。
子程序可以接受输入参数,并返回结果。
这种模块化的编程方法可以提高代码的可读性和可维护性。
Fortran程序通常由一个主程序(PROGRAM)和若干个子程序(SUBROUTINE)组成。
主程序是程序的入口点,而子程序则可以被主程序或其他子程序调用。
Fortran还提供了许多控制结构,包括条件语句(IF-THEN-ELSE)和循环语句(DO)等。
这些结构可以帮助程序在不同的情况下做出不同的决策和重复执行特定的代码块。
在写Fortran程序时,编码风格非常重要。
良好的编码风格可以使程序更易于阅读和理解,减少错误的发生。
在Fortran中,常用的编码风格包括正确缩进、适当的变量命名和注释的使用等。
总结起来,Fortran是一种强大而易于学习的编程语言,特别适用于数学计算和科学工程领域。
通过掌握Fortran的基本语法和编码风格,您将能够编写高效且可靠的程序。
fortran常用算法程序集
fortran常用算法程序集Fortran是一种高级编程语言,广泛应用于科学计算和数值分析领域。
它的强大之处在于它提供了丰富的算法库,使程序开发人员能够快速实现各种常见算法。
本文将介绍一些Fortran常用的算法程序集,帮助读者更好地理解和应用这些算法。
一、线性代数算法线性代数是科学计算和数值分析的基础,Fortran提供了许多用于求解线性方程组、矩阵分解和矩阵运算的算法。
其中一些常用的算法包括:1. 高斯消元法高斯消元法是一种求解线性方程组的方法,可以将线性方程组转化为上三角或下三角矩阵,并进一步求解。
Fortran提供了多种高斯消元法的实现,如LU分解法和托伯利兹矩阵法。
2. 特征值与特征向量计算特征值与特征向量计算是矩阵分解的一种重要问题。
Fortran提供了多种算法来计算矩阵的特征值与特征向量,如幂法、反幂法、QR算法等。
3. 矩阵乘法和矩阵求逆矩阵乘法和矩阵求逆是线性代数中常见的操作。
Fortran提供了多种高效的算法来实现矩阵乘法和矩阵求逆,如Strassen算法、LU分解法等。
二、数值计算算法数值计算算法广泛应用于科学计算、数值模拟和数据分析等领域。
Fortran提供了多种数值计算算法的实现,如数值积分、函数逼近、插值算法等。
以下是一些常用的数值计算算法:1. 数值积分数值积分可以用来对函数进行近似求积,求解曲线下面积或计算定积分。
Fortran提供了多种数值积分方法的实现,如梯形法则、辛普森法则和龙贝格方法等。
2. 函数逼近函数逼近是将多项式或其他数学函数与给定函数进行拟合,用于简化函数计算或数据分析。
Fortran提供了多种函数逼近的方法,如最小二乘逼近、最大误差逼近等。
3. 插值算法插值算法用于根据已知的离散数据点估计未知点的值。
Fortran提供了多种插值算法的实现,如拉格朗日插值法、牛顿插值法和样条插值法等。
三、优化算法优化算法用于求解最优化问题,如寻找函数最大值或最小值的点。
三对角矩阵fortran语言
三对角矩阵fortran语言什么是三对角矩阵,如何在Fortran语言中实现它,以及它在计算机科学和工程中的应用。
一、什么是三对角矩阵?三对角矩阵是一种具有特殊结构的矩阵。
它的主对角线上的元素是非零的,而其上方和下方的对角线上元素均为零。
换句话说,只有主对角线和两个相邻的对角线上的元素不为零,其余元素均为零。
三对角矩阵可以表示为以下形式:[A B 0 0 …0][C A B 0 …0][0 C A B …0][………………][0 0 …C A B][0 0 …0 C A]其中,A、B和C代表矩阵中的非零元素。
三对角矩阵通常是方阵。
二、如何在Fortran语言中实现三对角矩阵?在Fortran语言中,我们可以使用数组来表示矩阵,并利用循环结构来计算和操作矩阵的元素。
为了实现三对角矩阵,我们可以利用Fortran中的二维数组来存储矩阵,同时使用循环结构来对主对角线和相邻对角线上的元素进行赋值。
下面是一个简单的示例代码:fortranprogram tridiagonal_matriximplicit noneinteger, parameter :: n = 5 ! 矩阵的维度integer :: i, jreal :: A(n,n) ! 存储三对角矩阵的数组! 初始化矩阵do i = 1, ndo j = 1, nA(i, j) = 0.0 ! 将所有元素初始化为零end doend do! 赋值主对角线和相邻对角线上的元素do i = 1, nA(i, i) = 2.0 ! 主对角线上的元素if (i < n) thenA(i, i+1) = 1.0 ! 上方相邻对角线上的元素A(i+1, i) = 1.0 ! 下方相邻对角线上的元素end ifend do! 打印矩阵write(*, '(5F5.2)') (A(i, j), i = 1, n, j = 1, n)end program tridiagonal_matrix通过上述示例代码,我们可以创建一个5x5的三对角矩阵,并将其打印出来。
FORTRAN数值方法及其在物理学中应用7
FORTRAN数值方法及其在物理学中应用7 FORTRAN数值方法及其在物理学中应用7Fortran(Formula Translation)是一种古老但仍然被广泛应用的编程语言,特别适用于数值计算和科学计算。
在物理学中,Fortran经常被用于实施各种数值方法,帮助解决复杂的物理问题。
本文将介绍Fortran数值方法的一些常见应用。
1. 数值积分方法:Fortran在物理学中广泛用于数值积分方法的实现。
例如,将连续函数转化为离散值的数值积分可以通过Simpson法则或梯形法则来实现。
这些方法需要通过将积分区间划分为若干小区间,然后在每个小区间上进行数值积分。
Fortran提供了丰富的数值计算库,例如BLAS和LAPACK,用于实现高效的数值积分方法。
2. 常微分方程求解:物理学中有许多问题可以建模为常微分方程(ODEs)。
Fortran数值方法可以用于求解这些ODEs,并获得系统的解析解。
例如,Euler法,龙格-库塔法和阻尼牛顿法(Damped Newton's method)等都是常见的Fortran数值求解方法。
数值ODE求解器在物理学中广泛应用于模拟和预测各种物理系统的行为。
3. 矩阵运算:物理学中的很多问题可以表示为矩阵运算。
Fortran提供了高效的矩阵操作库,例如BLAS和LAPACK,可以用于求解矩阵方程,计算特征值和特征向量等。
这些方法在量子力学、光学和电磁场模拟等领域中非常重要。
4. 最小二乘拟合:物理学中常常需要从实验数据中拟合出最佳的理论模型。
最小二乘法是一种常用的拟合方法,可以用于找到最适合实验数据的参数。
Fortran中的最小二乘法库可以用于最小二乘求解,并提供了各种拟合指标,例如平方和误差等。
5. 数值优化方法:物理学中的许多问题可以视为找到使一些目标函数最小化或最大化的最优化问题。
Fortran数值优化方法可以用于解决这些问题。
例如,梯度下降法、粒子群优化算法和遗传算法等都是常见的Fortran数值优化方法,用于寻找最优解或全局最优解。
《2024年基于MapReduce编程模型的Fortran代码重构》范文
《基于MapReduce编程模型的Fortran代码重构》篇一一、引言随着计算技术的发展,高性能计算逐渐成为科学研究与工程应用中不可或缺的环节。
Fortran作为一种传统的科学计算语言,在高性能计算领域有着广泛的应用。
然而,传统的Fortran代码在处理大规模数据时,往往存在效率低下的问题。
为了解决这一问题,本文提出了一种基于MapReduce编程模型对Fortran代码进行重构的方法,以提高代码的并行处理能力和计算效率。
二、MapReduce编程模型MapReduce是一种编程模型,用于大规模数据集的并行处理。
其核心思想是将复杂的计算过程分解为两个步骤:Map阶段和Reduce阶段。
在Map阶段,程序将输入数据分解为多个键值对,并对每个键值对执行相同的操作。
在Reduce阶段,程序将Map 阶段输出的中间结果进行合并和归约,以得到最终结果。
三、Fortran代码重构策略1. 代码拆分:将原有的Fortran代码拆分为多个函数或子程序,每个函数或子程序负责处理一部分数据或执行特定的计算任务。
2. 并行化:利用OpenMP、MPI等并行计算技术,将拆分后的代码进行并行化处理,以提高计算效率。
3. Map阶段实现:在Map阶段,将输入数据按照一定的规则划分为多个键值对,并调用相应的函数或子程序进行处理。
4. Reduce阶段实现:在Reduce阶段,将Map阶段输出的中间结果进行合并和归约,以得到最终结果。
四、Fortran代码重构实例以一个典型的科学计算问题为例,如大规模矩阵乘法运算。
原Fortran代码可能采用串行方式实现,计算效率低下。
通过基于MapReduce编程模型的Fortran代码重构,我们可以将矩阵乘法运算拆分为多个小规模的矩阵乘法任务,并利用并行计算技术进行加速处理。
具体实现步骤如下:1. 代码拆分:将原始矩阵拆分为多个子矩阵,每个子矩阵负责一部分计算任务。
2. 并行化:利用OpenMP或MPI等并行计算技术,对每个子矩阵的乘法运算进行并行化处理。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
WRFCould not chdir to home directory /home/walilai: Stale NFS file handle/usr/X11R6/bin/xauth: error in locking authority file/home/walilai/.XauthorityBugINF ! 1/0=無限大(可用程式寫寫看)NaN ! 1/0=無限大read(unit,*) ! 注意unit數字使用限制Segmentation fault!矩陣沒定義到的空間(I or j…)或從第6行開始寫(須從7)doi=i+1 ! 寫在這會出現Segmentation faultstation(i)=' ' ! 可能是if沒被定義到I,stationread(10,21,end=31) st,lon,lat,mm,dd,hh,rainif(mm.eq.7.and.dd.eq.13.and.hh.eq.24)theni=i+1 !OKstation(i)=' ' !OKdo j=1,idoi=i+1c station(i)=' ' !若沒c掉,會出現Segmentation fault read(11,110,end=97) st,yy,mm,dd,hh,c,raindo_ud: off end of record:k=0 !a.out OK!!c do j=1,a !c處出現do_ud: off end of record:c k=j-1 do i=1,awrite(11,rec=3*k+1)station(i)write(11,rec=3*k+2)st_lon(i)write(11,rec=3*k+3)st_lat(i)k=k+1enddoopen: null file name!沒有讀到檔案,或是原本要輸入end才可結束程式的地方案到其他按鍵list in: end of file!檔案讀不完全,可能status=unknownopen: 'new' file exists:status!因檔案存在,不可用newdo_ud: end of file!字串字數沒寫對apparent state: unit 10 named f060200_112.5_35.datfmt: read unexpected character !找以下檔案會發現錯誤apparent state: unit 11 named c-cwb200807.txtlast format: (A6,I4,I2,I2,I2,A47,F5.1)!所讀數字出現文字T(雨跡)!format格式讀錯!可能計算出現nan而失敗invalid integer: read unexpected characterapparent state: unit 87 named autop2010-07.txtlast format: list iolately reading sequential formatted external IOAborted! read(87,100,end=920)st,yy,mm,dd,hh,c,rain將格式由*改為100就可fmt: end of fileapparent state: unit 92 named remove_2.txtlast format: (A6,1X,F6.1,1X,F5.1,2X,F7.1,2X,F7.1)lately reading sequential formatted external IOAborted!沒有close(??)或是無線迴圈的end=900,或是開檔時只有11個,卻讀到do i=1,12,或是1D矩陣寫成3D矩陣,或是已開啟資料寫兩次read()list in: end of fileapparent state: unit 31 named fort.31last format: list iolately reading direct formatted external IOAborted!txt_name(31)='xxx'之後,後面出現read(31,*) txt_name(31) /tmp/ccI6mztx.o(.text+0x139d): In function `MAIN__':: undefined reference to `filt25_'副程式最後要加 include ‘filt25.f’do_ud: off end of record !recl那nx和ny寫錯apparent state: unit 85 named grads_5.datstnmap -i xxx.ctlName of binary data set:/work4/tony1321/rain/200807/0710_daily.datNumber of times in the data set: 24Number of surface variables: 7Number of level dependent variables: 0Starting scan of station data binary file. Binary data file open:/work4/tony1321/rain/200807/0710_daily.datProcessing time = 1Invalid station hdr found in station binary filePossible causes: Invalid level count in hdrDescriptor file mismatchFile not station dataInvalid relative timelevs = -1027080205 flag = -998639206 time = 1.4013e-45Processing time = 24 !okTime = 24 has stn count = 415Max reports per time: 415 reports at t = 1Max data elements in largest report: 5Version 2 Station map file created: rain.mapstnmap: WARNING!! This stnmap file can only be accessed by GrADS Version 1.9b4stnmap: WARNING!! However, GrADS Version 1.9b4 can read both versionsstnmap: COMMENT -- use the -1 command line option to create a map for older GrADS versionsgradsshaded會把contour蓋過去,因此要先寫shaded,再寫contourq file 1.2……:不需開ctl即可看到各種變數若ctl觀測為一段時間,沒再gs定義set t XX,他會自動跑下去Constant field. Value = 2.60788e-09'set font 0-5 調整字體'set cstyle 0-7 調整線的類型'set t 'tt 要空格(空格 = 等於)'set digsize 數字'設定風標長度'set mpdraw off' 將grads內建地圖關掉Grads裡都用小寫較不會有問題('d Tem' ok但'd 'Tem不行)會出現畫出來的圖亂讀,譬如讀溫度卻跑出壓力,那可能式氣壓層數寫錯,如有25層卻只寫到21層Define error: Define not yet valid for station dataDefault file is a station data file! station data必放在第一個Low Level I/O Error: Read error on data fileData file name = grads.datError reading 81 bytes at location 9882 !此處排序錯誤Data Request Error: Error for variable 'inttem'Error ocurred at column 1Low Level I/O Error: Read error on data fileData file name = 20091012.dat !只有一層但ctl卻寫21層Error reading 141 bytes at location 417501Data Request Error: Error for variable 'dbz'Error ocurred at column 1DISPLAY error: Invalid expressionExpression = dbzCannot plot color bar: No shading informationDISPLAY error: Invalid expression Expression = inttemData Request Warning: Request beyond file limitsCannot plot data - all undefined values !set dfile順序要放對Open Error: Can't open description file !開到不存在的ctl檔SET Error: No files open yetDISPLAY error: no file open yetSyntax Error: Invalid Operand !ctl檔名稱寫錯,沒讀到'ccslp' not a variable or function nameError ocurred at column 1DISPLAY error: Invalid expressionExpression = ccslpLow Level I/O Error: Read error on data fileData file name = 20091012_RCWF.datError reading 141 bytes at location 79524Data Request Error: Error for variable 'dbz'Error ocurred at column 1DISPLAY error: Invalid expressionExpression = dbzCannot plot color bar: No shading information! 只有一層z=1,卻寫到set z 5而讀錯Open Error: Can't open description fileSET Error: No files open yetSET Error: No files open yet !ctl檔名寫錯SET Error: No files open yetSET Error: No files open yetDISPLAY error: no file open yetCannot plot color bar: No shading informationOpen Error: Can't open description fileSyntax Error: Invalid Operand'u' not a variable or function nameError ocurred at column 1DISPLAY error: Invalid expressionExpression = u!gs裡ctl超過一個時,需每個都要正確存在,否則會出現上面這樣,也就是全部ctl的路徑要全寫和高度(層)有關的變數,用'uprs'代表比較好在畫ec資料時須注意1.Data Request Warning: Request beyond file limits'set t 'tt位置順序寫錯有兩個ctl以上若時間定義不同,地形(t=1)2.如果gradsnc 裡面 t while迴圈順序寫錯,他會亂讀ex:tt=1while(tt<=124)第一個ctl#tt=tt+1 <---這裡必須關掉迴圈,否則會亂讀#endwhile#tt=1#while(tt<=12)第二個ctltt=tt+1 endwhile <---這迴圈只能寫在這'd rain''set gxout stat''d rain'#print resultline=sublin(result,8)word=subwrd(line,5)rmax=substr(word,1,4)try try 看'draw title case1A 9-21 rainfall ('data')''printim case1a_9_21_rainfall.gif x1600 y1200 white'Set clab %gK 在數字後寫上單位KWarning from LOG: Data has 2 values <= zeroThese were set to the undefined value !同時定義到lev和zOpen Error: Inconsistent time countCount in station map file = 18Count in descriptor file = 1 !須做stnmap –I XXX.ctl The data file was not opened. !開啟的ctl時間沒有連續Cat完dat檔時網格點可用4*nx*ny*時間來驗證若已經有dfile時,則不可再對變數做類似var.2 or var.3SET error: Missing or invalid arguments for GXOUT optionData Request Error: File number out of range ! 開錯ctl檔 Variable = ter.3Error ocurred at column 1DISPLAY error: Invalid expressionExpression = ter.3(t=1)Test2.fimplicit noneinteger ii,jj,i,j,a,bparameter (ii=2,jj=2,a=2,b=2)real x(ii,jj),y(ii,jj),z(ii,jj),w(ii,jj)open(10,file='test2.txt')do i=1,iiread(10,*)(x(i,j),j=1,jj)enddodo i=1,iic do j=1,jjc write(*,*)x(i,j) !只能寫成一排write(*,*)(x(i,j),j=1,jj) !可以寫成矩陣c enddoenddoc write(*,*)x(1,2)close(10)end programa.out11. 12. !x(1,2)=1222. 23. 991.fprogram qqimplicit noneinteger a,b,c(9,9)do a=1,9do b=1,9c(a,b)=a*benddoenddodo b=1,9write(*,*)(a,'×',b,'=',c(a,b),a=1,9)enddo! c(a,b) 寫成一列stopend./a.out !製作九九乘法1X 1= 1 2X 1= 2 3X 1= 3 4X 1= 4 5X 1= 5 6X 1= 6 7X 1= 7 8X 1= 8 9X 1= 91X 9= 9 2X 9= 18 3X 9= 27 4X 9= 36 5X 9= 45 6X 9= 54 7X 9= 63 8X 9= 72 9X 9= 81992.fProgram matrix tttttimplicit noneinteger a,b,c(9,9)open(6,file='992.txt',status='unknown',form='formatted')! 6=給予這個檔案的編號do b=1,9do a=1,9c(a,b)=a*benddoenddo! ------------------------------do b=1,9write(6,*)(a,'x',b,'=',c(a,b),a=1,9)enddo! 此處6是從6檔案寫出來stopendvi 992.txt1x 1= 1 2x 1= 2 3x 1= 3 4x 1= 4 5x 1= 5 6x 1= 6 7x 1= 7 8x 1= 8 9x 1= 91x 9= 9 2x 9= 18 3x 9= 27 4x 9= 36 5x 9= 45 6x 9= 54 7x 9= 63 8x 9= 72 9x 9= 81993.fProgram win or go homeimplicit noneinteger a,b,cwrite(*,*)"a="read(*,*)awrite(*,*)"b="read(*,*)bc=a*b write(*,*) a,'X',b,'=',cif(c<60 .or. a<10) then write(*,*) "go home"elsewrite(*,*) "win"endifstopend./a.outa=50b=3250X 32= 1600WinTest1.fProgram summer or playoff implicit noneinteger a,b,c,dreal :: e,fwrite(*,*)"San Antonio " write(*,*)"win game"read(*,*) ab=82-awrite(*,*)'lost game'write(*,*)be=a/bwrite(*,*)'winnig',ewrite(*,*)"Final Standings"write(*,*)a,' -',bwrite(*,*)"Phoenix "write(*,*)"win game"read(*,*)cd=82-cwrite(*,*)'lost game'write(*,*)df=c/dwrite(*,*)'winnig',fwrite(*,*)"Final Standings"write(*,*)c,' -',dif(e>f) thenwrite(*,*)'San Antonio in the playoff game' else if (e .eq. f) thenwrite(*,*)'The same record'elsewrite(*,*)'Phoneix in the playoff game'endifstopend ./a.outSan Antoniowin game50lost game32winnig 1.Final Standings50 - 32Phoenixwin game43lost game39winnig 1.Final Standings43 - 39The same recordFile.f 快速開啟許多檔案方法implicit noneinteger i,j,kcharacter file_name*10open(10,file='file_name.txt',status='old')do k=1,3read(10,*)file_nameopen(k,file=file_name,status='unknown') if(k.eq.1) thenwrite(k,*) 'one'else if(k.eq.2) thenwrite(k,*) 'two'else if(k.eq.3) thenwrite(k,*) 'third'endifenddoclose(k)close(10)end programHw3.fimplicit noneinteger jreal k,dtcomplex*16 ya(21)complex*16 yf(21)complex*16 yb(21)complex*16 ireal rya(21),ryf(21),ryb(21)real iya(21),iyf(21),iyb(21)i=cmplx(0.,1.)open(10,file='hw3.dat',form='unformatted',access='direct' ,recl=4*21) !21代表畫圖時x軸從時間0到20 do j=1,21ya(j)=exp(i*0.5*(j-1.))enddodo j=1,21rya(j)=dreal(ya(j)) !取複數實部 enddodo j=1,21iya(j)=dimag(ya(j))enddoyf(1)=1. !ICyf(2)=1+0.5*i !ICdo j=2,20yf(j+1)=yf(j-1)+i*yf(j)enddodo j=1,21ryf(j)=dreal(yf(j))enddodo j=1,21iyf(j)=dimag(yf(j))enddoyb(1)=1.yb(2)=1+0.5*ido j=2,21yb(j+1)=yb(j)+0.25*i*(3.*yb(j)-yb(j-1))enddodo j=1,21ryb(j)=dreal(yb(j))enddodo j=1,21iyb(j)=dimag(yb(j))enddowrite(10,rec=1) (rya(j),j=1,21) !不能用do迴圈寫 write(10,rec=2) (ryf(j),j=1,21) !用do寫跑出來的圖 write(10,rec=3) (ryb(j),j=1,21) !會怪怪的write(10,rec=4) (iya(j),j=1,21)write(10,rec=5) (iyf(j),j=1,21)write(10,rec=6) (iyb(j),j=1,21)close(10)end programhw3.ctldset ^hw3.datundef -999.9title hw3xdef 21 linear 0 1 # 0:x軸的起始點,間隔為1 ydef 1 linear 0 1tdef 1 linear 00z01nov2010 12hrzdef 1 linear 1 1vars 6rya 0 99 qryf 0 99 wryb 0 99 eiya 0 99 riyf 0 99 tiyb 0 99 yendvarshw3.gs'reinit''open hw3.ctl''set button bcol1' #'set xlab auto''set grads off'#'set xaxis 0 20' #會造成座標軸重複覆蓋#'set yaxis -2 2 1' #'d ryb'#'d ryf'#'d rya'#'draw title real''d iyb''d iyf''d iya''draw title imagine''draw xlab time''draw ylab pi(t)'1000.gs'reinit''open 198706_part_all.ctl''set grads off' # 把圖下方的grads去除'set lon 90 150' # 取經度90~150'set lat 5 45' # 取緯度5~45'set mpdset hires''set lev 500' # 500hPa高度'set t 1 30' # 所取時間為1~30tt=1while(tt<=30) # t在30以下都可畫出'c' # 防止在每個t時間所畫出的圖被一直蓋上'set t 'tt # t=上面的tt'set gxout stream''set ccolor 2''d uprs*5;vprs*5''set gxout contour''d zprs/9.8'pull aaqtt=tt+1 # 按下任意鍵後再畫出下一張圖endwhile08+09.gs#average (08+09)/2 height'reinit''open 200806_part_all.ctl' # 自動等於dfile 1 'open 200906_part_all.ctl' # 自動等於dfile 2#'set mproj nps''set grads off''set grid off''set lon 90 150''set lat 5 45''set mpdset hires'zz=1while(zz<=10) 'clear' # 清除上一張圖'set z 'zz'set gxout contour''set dfile 1''define avea=ave(zprs,t=1,t=60)' # avea是自行取的任意變數'set dfile 2''define aveb=ave(zprs,t=1,t=60)''define avez=(avea+aveb)/2''d avez/9.8''draw xlab longitude''draw ylab latitude''draw title (08+09)/2'pull abczz=zz+1endwhile19870623.gs'reinit''open 198706_part_all.ctl'#'enable print test'#============================='set grads off''set lon 115 127''set lat 20 28''set lev 850''set t 45''set mpdset hires'#=============================#'set gxout shaded'#=============================#'set ccolor 4''set cstyle 7''d tprs'#============================='set gxout contour''d zprs/9.8'#=============================#'set ccolor 6''set gxout barb''d uprs*2;vprs*2''set string 1'#'set strsiz 2 2''draw title 850hPa 1987-06-23-00:00'1.gs'reinit''open 198706_part_all.ctl''set grads off''set grid off''set lon 90 150''set lat 5 45''set mpdset hires' zz=1while(zz<=10)'c''set z 'zz'set ccolor 5''set cstyle 7''define avet=ave(tprs,t=1,t=60)' 'd avet'#'set gxout contour'#'define avez=ave(zprs,t=1,t=60)' #'d avez/9.8'*'set gxout barb'*'set ccolor 1'*'define aveu=ave(uprs,t=1,t=60)' *'define avev=ave(vprs,t=1,t=60)' *'d aveu*5;avev*5'pull qqqif(z=1)'draw title 123';endifif(z=2)'draw title 123';endifzz=zz+1endwhile2_years_average.gs'reinit''open 198706_part_all.ctl''open 200806_part_all.ctl''set grads off''set lon 105 130''set lat 10 40''set mpdset hires''set dfile 1''set dfile 2'zz=1while(zz<=9)'c''set gxout contour''define avez=ave(zprs.1,t=1,t=60)' #兩年六月份平均'define avea=ave(zprs.2,t=1,t=60)''set z 'zz'd (avea+avez)/9.8'pull abczz=zz+1endwhileaverage.gs'reinit''open 198706_part_all.ctl''set grads off' 'set grid off''set lon 90 150''set lat 5 40''set mpdset hires'zz=1while(zz<=9)'c''set grads off''set gxout barb''define aveu=ave(uprs,t=1,t=60)' 'define avev=ave(vprs,t=1,t=60)' 'd aveu*4;avev*4''set gxout contour''set z 'zz'define avez=ave(zprs,t=1,t=60)' 'd avez/9.8'*'set cstyle 2'*'define avet=ave(tprs,t=1,t=60)' *'d avet'*no sensepull abczz=zz+1endwhilewhile z=1,9if(z=1)'draw title 1000hPa';endif if(z=2)'draw title 850hPa';endif endwhilediscrete.gs'reinit''open 198706_part_all.ctl''set lon 100 140''set lat 15 40''set mpdset hires''set t 1'zz=1while(zz<=9)'c''set grads off'*'set gxout barb'*'d uprs*4;vprs*4''set gxout contour''set z 'zz'd zprs/9.8''set cstyle 2' 'set ccolor 2''d tprs'pull abczz=zz+1endwhilewhile(zz<=9)if(z=1);fmon='1000hpa';endifif(z=2);fmon='850hpa';endifendwhilef1987.gs'reinit''open 198706_part_all.ctl''set fwrite f1987.dat' 開起.dat 'set gxout fwrite' 將資料寫到.dat'set lon 120''set lat 25''set t 1 60''set z 1 14'tt=1while(tt<=60)'c''set t 'ttzz=1while(zz<=14)'c''set z 'zz'd zprs' # 重力位高度,需除9.8'd tprs''d uprs''d vprs''d rh'zz=zz+1endwhilett=tt+1endwhileimplicit noneinteger k,t,zrealzprs(14,60),tprs(14,60),uprs(14,60),vprs(14,60),rh(14,60) open(10,file='f1987.dat',form='unformatted',status='old',access='direct',recl=4)open(11,file='f1987.txt')k=0do t=1,60 #和f1987.gs檔時間順序一樣do z=1,14read(10,rec=5*k+1) zprs(z,t) #注意rec寫法,不是從10開始 read(10,rec=5*k+2) tprs(z,t)read(10,rec=5*k+3) uprs(z,t)read(10,rec=5*k+4) vprs(z,t)read(10,rec=5*k+5) rh(z,t)k=k+1enddo !zenddo !tdo t=1,60do z=1,14write(11,'(f10.2)') zprs(z,t)/9.8 write(11,'(f10.2)') tprs(z,t)write(11,'(f10.2)') uprs(z,t)write(11,'(f10.2)') vprs(z,t)write(11,'(f10.2)') rh(z,t)enddo !zenddo !tclose(10)close(11)end program畫台灣地形 ter2.gs'reinit''open ter2.ctl'#============================='set grads off''set lon 120.5 122''set lat 21.5 24''set mpdset hires''set xlint 0.5' #x軸數字間隔'set ylint 0.5' #y軸數字間隔'set strsiz 5 8'#'set ccolor 3'#============================='set gxout shaded'#'set cint 1013'#'set cmax 1''set cmin 1' #第一層以下無顏色'set cstyle 4''d ter'#'set gxout contour'#'d ter''set display color white''printim ter2.gif gif'2009-10.fimplicit noneinteger i,j,k,a,b,c,d !a.d=總資料數 parameter (a=285936,b=384,c=930,d=30) !b.c=測站數 character (20) station_auto(a),station_cwb(c)character (11) x1,x2,x3,x4,x5,x6,x7,x8,x9,x10character (11) y1,y2,y3,y4,y5,y6,y7,y8,y9,y10character (11) z1,z2,z3,z4,z5,z6,z7,z8,z9,z10character (11) w1,w2,w3,w4,w5,w6,w7,w8,w9,w10 character (11) time,ti(a),time_cwb,ti_cwb(c)real total,total_cwb,mr_cwb(c)real p,t,xx,xxx,rain,xxxxreal pp(a),tt(a),yy(a),yyy(a),rr(a),yyyy(a),dr(a)real rain_cwbreal rr_cwb(c)open(10,file='autop2009-10.txt')open(12,file='d-cwb200910.txt')open(11,file='2009-10.txt')do i=1,aread(10,*) station_auto(i),time,p,t,xx,xxx,rain,xxxxti(i)=timerr(i)=rainenddodo i=1,c !read:讀一整列(.txt橫的資料)read(12,*)station_cwb(i),time_cwb,x3,x4,x5,x6,x7,x8,x9,x10, !土法煉鋼$ y1,y2,y3,y4,y5,y6,y7,y8,y9,y10,$ z1,z2,z3,z4,z5,z6,z7,z8,z9,z10,$ w1,w2,w3,w4,w5,w6,w7,w8,w9,w10,rain_cwbti_cwb(i)=time_cwbrr_cwb(i)=rain_cwbenddoc=================================================auto station do k=1,btotal=0.do i=1.+(k-1.)*744.,744.*kif(rr(i).eq.-9996.)thenrr(i)=0.else if(rr(i).eq.-9991.)thenrr(i)=0.else if(rr(i).eq.-9997.)thenrr(i)=0.endiftotal=total+rr(i)dr(i)=totalenddoenddodo i=744,a,744 !字串需用自由格式*c write(11,'(45x,f7.1)')dr(i) !45x:空45格write(11,*)station_auto(i),ti(i),i,dr(i)enddoc==================================================cwb station do k=1,dtotal_cwb=0.do i=1.+(k-1.)*31.,31.*kif(rr_cwb(i).eq.-9998.)thenrr_cwb(i)=0.endiftotal_cwb=total_cwb+rr_cwb(i)mr_cwb(i)=total_cwbenddoenddodo i=31,c,31c write(11,'(45x,f7.1)') station_cwb(i),mr_cwb(i)write(11,*) station_cwb(i),ti_cwb(i),i,mr_cwb(i)enddoc========================================================= close(10)close(11)close(12)end programa.outC0A520 2009103124 744 199.C0A530 2009103124 1488 813.C0A540 2009103124 2232 782.C0A550 2009103124 2976 1541.C0A560 2009103124 3720 860.5C0A570 2009103124 4464 746.C0A580 2009103124 5208 512.5C0A590 2009103124 5952 806.5C0A860 2009103124 6696 1052.5C0A870 2009103124 7440 405.5C0A880 2009103124 8184 561.C0A890 2009103124 8928 898.C0A920 2009103124 9672 443.5C0A930 2009103124 10416 502.1011rain.fimplicit noneinteger i,j,k,a,b,ccharacter station(550)*6,st*6 !數字超過測站即可integer mm,dd,hh,st_dd(10000),st_hh(10000),st_mm(10000) real rain,lon,lat,st_lon(10000),st_lat(10000550)real st_rain(10000)open(11,file='1011rain.txt',status='old')open(10,file='sister.txt',status='old')c============================================================= c stc COA520i=0doc i=i+1 !這個位置會出Segmentation faultread(10,*,end=31)st,lon,lat,mm,dd,hh,rainif(mm.eq.10.and.dd.eq.11.and.hh.eq.24)then !必寫,要不然會亂讀時間i=i+1 !這個位置OKstation(i)=' 'c if(mm.eq.10.and.dd.eq.11.and.hh.eq.24)thendo j=1,iif(st.eq.station(j).and.lon.eq.st_lon(j)t.eq.st_lat(j $).and.mm.eq.st_mm(j))then !測站、經緯度都重複,需用矩陣 i=i-1goto 32endifenddostation(i)=stc print*, station(i)st_lon(i)=lonst_lat(i)=latst_mm(i)=mmst_dd(i)=ddst_hh(i)=hhst_rain(i)=rain write(*,21)station(i),st_lon(i),st_lat(i),$ st_mm(i),st_dd(i),st_hh(i) !都用矩陣寫才可endifc write(*,21)station(i),st_lon(i),st_lat(i),c $ st_mm(i),st_dd(i),st_hh(i)32 continuecprint*,station(i),st_lon(i),st_lat(i),st_mm(i),st_dd(i),st_hh( i)enddocprint*,station(i),st_lon(i),st_lat(i),st_mm(i),st_dd(i),st_hh(i)31 continuek=iprint*,kc============================================================= i=0doi=i+1read(10,21,end=41)station(i),st_lon(i),st_lat(i),mm,st_dd(i) $ ,st_hh(i),rainif(st_dd(i).eq.10.and.st_hh(i).eq.24)thenc print*,st_hh(i)cwrite(11,21)station(i),st_lon(i),st_lat(i),mm,dd,hh,rainendifenddo41 continue21 format(a6,1x,f7.3,1x,f6.3,1x,i2,1x,i2,1x,i2,1x,f7.2) close(11) ! 1011rain.txtclose(10) ! ister.txtend programa.outC1V450 120.376 22.540 10 11 24 0.00C1V460 120.884 23.438 10 11 24 6.00C1X010 120.252 23.332 10 11 24 0.00在地圖上標記號在grads下1. reinit2. open ter2.ctl3. set lon和lat4. d ter5. q w2xy 經度緯度6. X = ? Y = ?2008_2009.gs'reinit''open 200806_part_all.ctl''open 200906_part_all.ctl''set mpdset mres''set grads off''set lon 100 140' 'set lat 7 37''set lev 1000''set xlint 5''set xlint 5''set gxout contour''set dfile 1'#'define tprs.1=avg(tars,t=1,t=60)' #可不用設這個'set t 10 ' #有設這個時'define ctone=tprs.1-273' #變數用全英文比較好'set gxout contour''set dfile 2''set t 10 ''define cttwo=tprs.2-273''d (ctone+cttwo)/2'。