fortran版diamond4转grads的源程序

mofangbao

前两天看到论坛有朋友在问这方面的问题，看了一下四类数据的格式，发现还是算好转的，自描述比较清楚。
这是Micaps第四类数据的介绍：
============================

文件头：

diamond  4  数据说明（字符串）  年  月  日  时次  时效  层次（均为整数）经度格距  纬度格距  起始经度  终止经度  起始纬度  终止纬度（均为浮点数） 纬向格点数  经向格点数（均为整数）  等值线间隔  等值线起始值  终止值  平滑系数加粗线值（均为浮点数）

注：此类数据用于画格点数据的等值线。网格可以为经纬度网格，也可以为直角坐标网格。

1.      当使用直角坐标网格数据时：1）将等值线终止值改为-1（直角坐标在兰勃托投影下）或-2（直角坐标在麦开托投影下）或-3（直角坐标在北半球投影下）。2）把网格经度间隔和纬度间隔改为格点数据第一行最后一个点的经纬度。3）把起始经度和起始纬度改为格点数据第一行第一个点的经纬度。4）把终止经度和终止纬度改为格点数据最后一行最后一个点的经纬度。

2.      第4类数据文件可以直接用于填格点值。文件头中可以指定填图方式。指定方法为：1）把加粗线值改为-1，表示画等值线同时填图，2）改为-2表示只填图，不画等值线。

数据：

数据按先纬向后经向放（直角坐标网格时为先X方向后Y方向），均为浮点数。

===========================================================
    我首先使用了Meteoinfo进行转换，但是有一点小问题，等射月楼主放出新版后大家就可以直接用MeteoInfo画图了，转换格式有时候确实有点麻烦
==================
    这个转换要注意的主要是数据的存放顺序和grads中是否一致的问题，其他的就是简单的读入输出的问题了，因为本身就是格点数据，也不需要插值什么的，程序会自动生成CTL文件。
注意点：
1、如果您输入的是纯文件名而不是带路径的文件名，则需要在CTL中修改路径才能使用该CTL；
2、程序没有对时间进行自动分析，有兴趣的朋友可以写个子程序处理一下
    为了节约大家的金钱，下面是fortran的源代码和一些简单的说明：
===========
program main
real head(19) !用于存储自描述内容
real,allocatable::gdata(:,:)  !存储数据，由头文件中的格点数决定长度
integer x,y,i,j      !x,y方向的格点数
character filepath*100,desc(3)*50 !文件路径和文件第一行内容，desc主要在生成ctl时使用
!****输出提示内容，读入文件路径******
print*,"1、请输入文件路径，若为当前路径可直接输入文件名，路径中请不要出现空格"
print*,"2、如果输入的是文件名而非路径，请在运行结束后自行修改生成的CTL文件的路径"
print*,"3、没有对生成的CTL文件的时间进行处理，需要的请自行修改"
print*,""
print*,"输入文件名或路径："
read(*,'(a)')filepath
!*********打开文件循环读取********
open(1,file=''//trim(filepath)//'',status='old',err=101)
read(1,*)desc
read(1,*)head
x=int(head(13)+0.0001)
y=int(head(14)+0.0001)
allocate(gdata(y,x))    !申请内存
read(1,*)((gdata(j,i),i=1,x),j=1,y)
close(1)
!************根据数据存放格式选择生成的数据存放顺序************
open(1,file=''//trim(filepath)//'.grd',status='replace',form='binary',err=102)
!*********************************************
!  head(7)<0 经度（沿x）自东向西；
!  head(8)<0 纬度沿y自北向南，反之亦然
!*********************************************
if(head(7)<0)then   
if(head(8)<0)then
  do j=y,1,-1
   do i=x,1,-1
    write(1)gdata(j,i)
   enddo
  enddo
else
  do j=1,y
   do i=x,1,-1
    write(1)gdata(j,i)
   enddo
  enddo
endif
else
   if(head(8)<0)then    !比较常用的一种情况（如欧洲500hpa预报）
  do j=y,1,-1
   do i=1,x
    write(1)gdata(j,i)
   enddo
  enddo
else
  do j=1,y
   do i=1,x
    write(1)gdata(j,i)
   enddo
  enddo
endif
endif
close(1)
!************生成CTL，时间固定了，需要的自己改******************
open(1,file=''//trim(filepath)//'.ctl',status='replace',err=103)
write(1,*)"dset "//trim(filepath)//".grd"
write(1,*)"undef 9999"
write(1,*)"title "//trim(desc(3))
write(1,'(1x,a,i5,a,f10.4,1x,f8.4)')"xdef ",x," linear ",head(9),head(7)
write(1,'(1x,a,i5,a,f10.4,1x,f8.4)')"ydef ",y," linear ",head(12),head(8)*(-1)
write(1,'(1x,a,i5)')"zdef 1 levels ",int(head(6)+0.0001)
write(1,*)"tdef 1 linear 01jan1950 1dy"
write(1,*)"vars 1"
write(1,*)"m2g 1 99 diamond4tograds"
write(1,*)"endvars"
close(1)
goto 100
!*****************出错处理***********************************
101 print*,"读取源数据文件：",trim(filepath),"出错！"
goto 104
102 print*,"新建GrADS文件：",trim(filepath)//".grd","出错！"
goto 104
103 print*,"新建CTL文件：",trim(filepath)//".ctl","出错！"
104 close(1)
stop
!******************转换完成，给出提示************************
100 continue
print*,"转换完成，生成的新文件为："
print*,trim(filepath)//".grd"
print*,trim(filepath)//".ctl"
print*,"如果您输入的是纯文件名，请您补全ctl中的dset路径"
pause "按任意键退出"   !此处可以改为继续循环，进行批处理，但是上面也要做适当修改
deallocate(gdata)
end
===================代码结束================
嫌复制麻烦的朋友可以回复后下载

游客，如果您要查看本帖隐藏内容请回复

下面是两张对比图（还是Micaps1出的，请包涵，各位也可以用自己的数据在MeteoInfo出了之后对比一下）：

传说中的谁

很快可以开一个micaps转grads的专题了。

喜上眉梢

发表于 昨天 23:33 |只看该作者
很快可以开一个micaps转grads的专题了。

对，赞同，开个专题，将所有转GRADS的程序汇总与此，方便使用、方便学习，但首先得感谢版主和程序开发者。

gslcr

很好，学习了

zdhai

感谢楼主的无私奉献

nyqxtlh

很好，学习了

joannahere

支持！

滇西boy

顶啊~~~~~~

LEO

好东西！

mofangbao

LEO 发表于 2011-8-16 17:10

                               
登录/注册后可看大图

好东西！

呵呵 也是看了你的提问才写的 多来看看哦