Fortran求台风中心的算法

抹茶里

lqouc 发表于 2013-7-25 10:40

                               
登录/注册后可看大图

fnl在grads里面有专用的解码工具，你需要生成一个ctl，然后再提取数据，数据要的就是海平面气压。存方的格式 ...

只要气压就可以了吗？是不是变量里面的pressfc变量？之前还准备找个大风区域，在大风区域里面找最低气压。目前还没有头疼到程序这一步，目前已经将需要的变量用grads读取成grd格式的文件了，但是感觉好像读取错误，在fortran里面貌似不能正常使用这个grd数据文件{:3_60:}，是不是要再转换成十进制的。。

抹茶里

lqouc 发表于 2013-7-25 10:44

                               
登录/注册后可看大图

既然你想写个fortran的，目前好像论坛还没有现成的，如果写出来了，欢迎分享给大家。会给你加分的。

出来了可以贴上，额，但感觉好遥远的样纸{:3_53:}

lqouc

抹茶里 发表于 2013-7-25 11:17

                               
登录/注册后可看大图

只要气压就可以了吗？是不是变量里面的pressfc变量？之前还准备找个大风区域，在大风区域里面找最低气压。 ...

不需要换文本格式，就是用表面气压presssfc就行。
程序有问题可以吧脚本贴上来。

抹茶里

lqouc 发表于 2013-7-25 13:55

                               
登录/注册后可看大图

不需要换文本格式，就是用表面气压presssfc就行。
程序有问题可以吧脚本贴上来。

目前程序不报错了，但是print出来之后是一系列的0{:3_55:}。。有个地方尚不明白，我得出来的是关于一个pressfc变量的grd文件，如果不转格式，在fortran中直接运算的话，怎么知道它里面的数据是什么样的呢？最后我要得出经纬度的话，这个经纬度在fortran里面又怎样识别呢？因为之前是对txt文档数据进行处理，里面很容易就能看出来多少行多少列，我要的因素在哪几列这样子。。

lqouc

抹茶里 发表于 2013-7-25 14:43

                               
登录/注册后可看大图

目前程序不报错了，但是print出来之后是一系列的0。。有个地方尚不明白，我得出来的是关于一个pr ...

为什么是0的话没看你的程序我也说不好，可能是你的数据类型的问题，不知道你fortran什么程度。
想知道的话你就print好了，或输出一个文本文件，不可能什么都用文本格式的，太不方便了。
至于经维度你就用格点坐标乘以格距就好了啊，这些都是一个数组函数就出来了。

抹茶里

lqouc 发表于 2013-7-25 17:48

                               
登录/注册后可看大图

为什么是0的话没看你的程序我也说不好，可能是你的数据类型的问题，不知道你fortran什么程度。
想知道的 ...

fortran刚上手，程序没有报错，（数据网格是1*1的）但是结果输出来如下图：

fortran程序如下：
program main
implicit none
integer lon,lat,t
parameter(lon=360,lat=181,t=124)
integer i,j,k
real p(lon,lat,t)
open(10,file='f:\2013722\fnl\2004080131typhoon.grd',form='binary')
open(11,file='f:\2013722\wind2004\2004080131typhoon.dat',form='formatted')
do i=1,lon
   do j=1,lat
       do k=1,t
        read(10) P(i,j,k)
     end do
        end do
end do
do i=1,lon
   do j=1,lat
       do k=1,t
       write(11,*) P(i,j,k)
           print*, p(i,j,k)
           end do
        end do  
end do
stop
end

lqouc

抹茶里 发表于 2013-7-26 09:47

                               
登录/注册后可看大图

fortran刚上手，程序没有报错，（数据网格是1*1的）但是结果输出来如下图：

fortran程序如下：

你的数据小于你的循环数，这是报错的问题，吧grads脚本贴上来。还有循环顺序也不对，时间要在最外面。你的grad怎么写的fortran就怎么读。

抹茶里

lqouc 发表于 2013-7-26 10:32

                               
登录/注册后可看大图

你的数据小于你的循环数，这是报错的问题，吧grads脚本贴上来。还有循环顺序也不对，时间要在最外面。你的 ...

这是grd文件的ctl：
dset F:\2013722\fnl\2004080131typhoon.grd
undef 1e+15
title Sample GRID Data
xdef 360 linear 0 1
ydef 181 linear -90  1
tdef 124 linear 00Z01aug2004 360mn
zdef 26 levels 1000 975 950 925 900 850 800 750 700 650 600 550 500 450 400 350 300 250 200 150 100 70 50 30 20 10
vars 1  
pressfc     0 0 surface pressure.
endvars

读数据的gs文件写得很简单：
*'reinit'
'open f:\2013722\fnl\batch11.ctl'
'set fwrite  F:\2013722\fnl\2004080131typhoon.grd'
'set gxout fwrite'

*d=01
*h=00,06,12,18
*while(d<=31)
*'sdfopen  F:\2013722\fnl\fnl_0408'%d'_'%h'_00'
*say 'F:\2013722\fnl\fnl_0408'%d'_'%h'_00'
*'set x 100 180'
*'set y 0 50'
*'set x 0 475'
*'set y 181 281'
'd pressfc'
*'d ugrdmwl'
*'d vgrdmwl'
*'d tlml'
*'d qlml'
*'d ulml'
*'d vlml'
*'d prectot'
*'d pgentot'

*d=d+1
*'close 1'
*endwhile
'disable fwrite'

lqouc

抹茶里 发表于 2013-7-26 10:58

                               
登录/注册后可看大图

这是grd文件的ctl：
dset F:\2013722\fnl\2004080131typhoon.grd
undef 1e+15

grads脚本完全不对，高度和时间是需要用循环写的。当然你根本不需要高度。
所以你的ctl也就不对。
那么fortran自然是要出错的。

抹茶里

lqouc 发表于 2013-7-26 11:10

                               
登录/注册后可看大图

grads脚本完全不对，高度和时间是需要用循环写的。当然你根本不需要高度。
所以你的ctl也就不对。
那么 ...

原数据是fnl格式的，对应的ctl部分如下：
dset F:\2013722\fnl\fnl_0408%d2_%h2_00
index F:\2013722\fnl\112.idx
options template
undef 9.999E+20
*title fnl_20090715_00_00_c.c
*  produced by grib2ctl v0.9.12.5p33e
dtype grib 3
options yrev
ydef 181 linear -90.000000 1
xdef 360 linear 0.000000 1.000000
tdef 124 linear 00Z01aug2004 6hr
zdef 26 levels
1000 975 950 925 900 850 800 750 700 650 600 550 500 450 400 350 300 250 200 150 100 70 50 30 20 10
vars 93
no4LFTXsfc  0 132,1,0  ** surface Best (4-layer) lifted index [K]
no5WAVA500mb  0 230,100,500 ** 500 mb 5-wave geopot. height anomaly [gpm]
no5WAVH500mb  0 222,100,500 ** 500 mb 5-wave geopotential height [gpm]
ABSVprs 26 41,100,0 ** (profile) Absolute vorticity

再用gs文件读取成grd的话，还要用到循环吗？第一个ctl是对应的读取后的grd文件的，感觉只有气压一个变量了