爱气象，爱气象家园！ 
·设为首页 开启辅助访问 充值积分新帖子新回复热门贴精华帖

气象家园

 找回密码
 立即注册

QQ登录

只需一步，快速开始

搜索
气象家园»论坛 › 编程作图 › Fortran › fortran编程问题
12下一页
返回列表 发新帖
查看: 6758|回复: 14

[求助] fortran编程问题

[复制链接]
liyf
发表于 2014-5-5 23:00:26 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看，没有帐号？立即注册

x
本帖最后由 liyf 于 2014-5-7 22:19 编辑

! read 20090416/high/height;temper;t-td;uv
! lon 32~160                                                                                  interval: 4 degree
! lat 12~80                                                                                          interval: 4 degree
! time 1~12:  2009.04.16.20:00~04.22.08:00            interval: 12 hours
! U,v,height,t: 1000,925,850,700,500,400,300,250,200,150,100  hpa



REAL,PARAMETER    :: Omega=7.292e-5,R=6371e+3,PI=3.1415926,Delta=4,DeltaT=12
INTEGER,PARAMETER :: nx=33,ny=18,nz=11,nt=12
REAL F,EE,AA,BB,sign
INTEGER WtTopLevel(NX,NY,NT)
REAl P(NZ),deltP(10),lat(ny),lon(nx),sigmadeltD
CHARACTER timefile(12)*12 ,levelfile(11)*4
real,allocatable :: temper(:,:,:,:),q(:,:,:,:),u(:,:,:,:),     &
                                        v(:,:,:,:),height(:,:,:,:),w(:,:,:,:),     &
                                        qu(:,:,:,:),qv(:,:,:,:),adq(:,:,:,:),      &
                                        adqv(:,:,:,:),theta(:,:,:,:),Wt(:,:,:,:),  &
                                        vorg(:,:,:,:),voro(:,:,:,:),div(:,:,:,:),  &
                                        deltD(:,:,:),ttd(:,:,:,:),thetaP(:,:,:,:)
allocate (temper(NX,NY,NZ,NT),q(NX,NY,nz,NT),u(NX,NY,NZ,NT),   &
                 v(NX,NY,NZ,NT),height(NX,NY,NZ,NT),w(NX,NY,NZ,NT),    &
                 qu(NX,NY,NZ,NT),qv(NX,NY,NZ,NT),adq(NX,NY,NZ,NT),     &
                 adqv(NX,NY,NZ,NT),theta(NX,NY,NZ,NT),Wt(NX,NY,NZ,NT), &
                 vorg(NX,NY,NZ,NT),voro(NX,NY,NZ,NT),div(NX,NY,NZ,NT), &
                 deltD(NX,NY,NT),ttd(NX,NY,NZ,NT),thetaP(NX,NY,NZ,NT))

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
DATA P    /1000,925,850,700,500,400,300,250,200,150,100/
DATA deltP/75,  75, 50, 200,100,100,50, 50, 50, 50/
DATA timefile/'09041620.000','09041708.000','09041720.000','09041808.000','09041820.000','09041908.000', &
              '09041920.000','09042008.000','09042020.000','09042108.000','09042120.000','09042208.000'/
DATA levelfile/'1000','925','850','700','500','400','300','250','200','150','100'/


do i=1,nx
lon(i)=32+(i-1)*Delta
lon(i)=lon(i)*pi/180.
enddo
do i=1,ny
lat(i)=12+(i-1)*Delta
lat(i)=lat(i)*pi/180.
enddo
!--------------read temperature----------------
DO IZ=1,NZ        ! level
DO IT=1,NT        ! time
OPEN(11,file='d:/tianfen2/micaps/temper/'//levelfile(iz)//'/'//timefile(it))
    do i=1,4
    read(11,*)
        enddo
    do ij=NY,1,-1
        read(11,*)(temper(ii,ij,iz,it),ii=1,NX)           ! C
        enddo
CLOSE(11)
ENDDO
ENDDO
!--------------read u,v -------------------------
DO IZ=1,NZ        ! level
DO IT=1,NT        ! time
OPEN(11,file='d:/tianfen2/micaps/uv/'//levelfile(iz)//'/'//timefile(it))
    do i=1,3
    read(11,*)
        enddo
    do ij=36,1,-1
          if(ij>ny) then
          read(11,*)(u(ii,ij-ny,iz,it),ii=1,NX)           ! C
          else
      read(11,*)(v(ii,ij,iz,it),ii=1,NX)
          endif
        enddo
CLOSE(11)
ENDDO
ENDDo
!--------------read t-td -------------------------
DO IZ=1,NZ        ! level
DO IT=1,NT        ! time
OPEN(11,file='d:/tianfen2/micaps/t-td/'//levelfile(iz)//'/'//timefile(it))
    do i=1,4
    read(11,*)
        enddo
    do ij=NY,1,-1
        read(11,*)(ttd(ii,ij,iz,it),ii=1,NX)           ! C
        enddo
CLOSE(11)
ENDDO
ENDDO
!--------------read height -------------------------
DO IZ=1,NZ        ! level
DO IT=1,NT        ! time
OPEN(11,file='d:/tianfen2/micaps/height/'//levelfile(iz)//'/'//timefile(it))
    do i=1,4
    read(11,*)
        enddo
    do ij=NY,1,-1
        read(11,*)(height(ii,ij,iz,it),ii=1,NX)           ! C
        enddo
CLOSE(11)
ENDDO
ENDDO
!@@@@@@@@@@@@@@@@@@@@ read end @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!======== t,td => specific humidity:
do IZ=1,NZ
do IT=1,NT
  do i=1,nx
  do j=1,ny
  AA=temper(i,j,iz,it)-ttd(i,j,iz,it)
  BB=6.1078*exp(17.3*AA/(273.16+AA-35.86))
  q(i,j,iz,it)=622*BB/(p(iz)-0.378*BB)
  enddo
  enddo
enddo
enddo
!######################## geostrophic wind vorticity ###########################
do IZ=1,NZ
do IT=1,NT
  do i=2,nx-1
  do j=2,ny-1
  vorg(i,j,iz,it)=((height(i+1,j,iz,it)+height(i-1,j,iz,it)-2*height(i,j,iz,it))/((cos(lat(j))*Delta)**2)+&
  (height(i+1,j,iz,it)+height(i-1,j,iz,it)-2*height(i,j,iz,it))/Delta**2-&
  (height(i,j+1,iz,it)-height(i,j-1,iz,it))/(2*Delta*tan(lat(j))))/(9.8*Omega*R**2)
  enddo
  enddo
enddo
enddo
!######################## vorticity of observed wind     ########################
do IZ=1,NZ
do IT=1,NT
  do i=2,nx-1
  do j=2,ny-1
  voro(i,j,iz,it)=((v(i+1,j,iz,it)-v(i-1,j,iz,it))/(cos(lat(j))*Delta)-&
  (u(i,j+1,iz,it)-u(i,j-1,iz,it))/Delta+2*u(i,j,iz,it)*tan(lat(j)))/2/R
  enddo
  enddo
enddo
enddo
!######################## divergence                 ###########################
do IZ=1,NZ
do IT=1,NT
  do i=2,nx-1
  do j=2,ny-1
  div(i,j,iz,it)=((u(i+1,j,iz,it)-u(i-1,j,iz,it))/(cos(lat(j))*Delta)+&
  (v(i,j+1,iz,it)-v(i,j-1,iz,it))/Delta-2*v(i,j,iz,it)*tan(lat(j)))/2/R
  enddo
  enddo
enddo
enddo
!######################## vertical velocity          ###########################
do i=2,nx-1
do j=2,ny-1
  do IT=1,NT
  w(i,j,1,it)=0
  do IZ=2,NZ
  w(i,j,iz,it)=w(i,j,iz-1,it)+(div(i,j,iz,it)+div(i,j,iz-1,it))/2*(p(iz-1)-p(iz))
  enddo
  enddo
enddo
enddo
!--------          correct w      --------------
!- - - - - thermaldynamic method for w at the top level - - - - - - -
!- - - - - 100,150hpa:  reference to 100,150,200 hpa
do i=2,nx-1
do j=2,ny-1
  do IT=1,NT
  F=w(i,j,NZ,it)/(NZ*(NZ+1))
  do IZ=1,NZ-1
  w(i,j,iz,it)=w(i,j,iz,it)-iz*(iz-1)/F
  enddo
  w(i,j,NZ,it)=0
  enddo
enddo
enddo
!###############   vapor flux &  vapor flux divergence      ##################
do IZ=1,NZ
do IT=1,NT
  do i=2,nx-1
  do j=2,ny-1
  qu(i,j,iz,it)=u(i,j,iz,it)*q(i,j,iz,it)/9.8
  qv(i,j,iz,it)=v(i,j,iz,it)*q(i,j,iz,it)/9.8
  AA=u(i,j,iz,it)*(q(i+1,j,iz,it)-q(i-1,j,iz,it))/(cos(lat(j))*Delta)+v(i,j,iz,it)*(q(i,j+1,iz,it)-q(i,j-1,iz,it))/Delta
  adqv(i,j,iz,it)=(AA/(2*R)+q(i,j,iz,it)*div(i,j,iz,it))/9.8
  enddo
  enddo
enddo
enddo
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@


        OPEN(8,FILE='d:\tianfen2\chengxu\out2.grd',form='binary')
        do it=1,nt

        WRITE(8) (((u(i,j,iz,it),i=1,nx),j=1,ny),iz=1,nz)
        WRITE(8) (((v(i,j,iz,it),i=1,nx),j=1,ny),iz=1,nz)
        WRITE(8) (((q(i,j,iz,it),i=1,nx),j=1,ny),iz=1,nz)
        WRITE(8) (((temper(i,j,iz,it),i=1,nx),j=1,ny),iz=1,nz)
        WRITE(8) (((height(i,j,iz,it),i=1,nx),j=1,ny),iz=1,nz)

        WRITE(8) (((voro(i,j,iz,it),i=1,nx),j=1,ny),iz=1,nz)
        WRITE(8) (((vorg(i,j,iz,it),i=1,nx),j=1,ny),iz=1,nz)
        WRITE(8) (((div(i,j,iz,it),i=1,nx),j=1,ny),iz=1,nz)

        WRITE(8) (((adqv(i,j,iz,it),i=1,nx),j=1,ny),iz=1,nz)
    WRITE(8) (((qu(i,j,iz,it),i=1,nx),j=1,ny),iz=1,nz)
        WRITE(8) (((qv(i,j,iz,it),i=1,nx),j=1,ny),iz=1,nz)

        WRITE(8) (((w(i,j,iz,it),i=1,nx),j=1,ny),iz=1,nz)

        enddo
    CLOSE(8)

!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
deallocate (temper,q,u,v,height,w,qu,qv,adq,adqv,theta,Wt, &
                    vorg,voro,div,deltD,ttd,thetaP)
END
上面是我的程序，但出错了，错误如图所示，大概意思是没有文件，但文件存在，如图、。。。请教如何改正。。。谢谢

O(KI7HCEHWFND`4Z0Q4[3`8.jpg
H1AN4XRBP%X6C1C7)Q]%JXH.jpg
O4M4)9PV}9B54}C_J_)VTA8.jpg
密码修改失败请联系微信：mofangbao
回复

使用道具 举报

lqouc
发表于 2014-5-6 08:01:57 | 显示全部楼层
之前老五说的隐藏后缀的问题你检查了么？
可以在文件名那里加一个trim试试。
密码修改失败请联系微信：mofangbao
回复 支持 反对

使用道具 举报

三块表
发表于 2014-5-6 08:19:59 | 显示全部楼层
应该是后缀名隐藏了，000并不是你的后缀名
密码修改失败请联系微信：mofangbao
回复 支持 反对

使用道具 举报

nunu18
发表于 2014-5-6 10:52:20 | 显示全部楼层
没找到11所指的文件，你检查一下路径，后缀名被隐藏了吧
密码修改失败请联系微信：mofangbao
回复 支持 反对

使用道具 举报

一个不小心
发表于 2014-5-6 12:34:26 | 显示全部楼层
好高端的样子。。。。
密码修改失败请联系微信：mofangbao
回复 支持 反对

使用道具 举报

liyf
 楼主| 发表于 2014-5-7 22:15:46 | 显示全部楼层
lqouc 发表于 2014-5-6 08:01
之前老五说的隐藏后缀的问题你检查了么？
可以在文件名那里加一个trim试试。

检查了，后缀名没隐藏
密码修改失败请联系微信：mofangbao
回复 支持 反对

使用道具 举报

liyf
 楼主| 发表于 2014-5-7 22:23:21 | 显示全部楼层
nunu18 发表于 2014-5-6 10:52
没找到11所指的文件，你检查一下路径，后缀名被隐藏了吧

没有隐藏后缀名。。。当寻找文件d:\tianfen2\micaps\temper\1000没有出错，但求925时就出错了，表示不理解。。。
密码修改失败请联系微信：mofangbao
回复 支持 反对

使用道具 举报

liyf
 楼主| 发表于 2014-5-7 22:26:53 | 显示全部楼层
lqouc 发表于 2014-5-6 08:01
之前老五说的隐藏后缀的问题你检查了么？
可以在文件名那里加一个trim试试。

没有隐藏后缀名。。。当寻找文件d:\tianfen2\micaps\temper\1000没有出错，但求925时就出错了，表示不理解。。。
密码修改失败请联系微信：mofangbao
回复 支持 反对

使用道具 举报

liyf
 楼主| 发表于 2014-5-7 22:27:56 | 显示全部楼层
一个不小心 发表于 2014-5-6 12:34
好高端的样子。。。。

只不过是编码多了些。。。分开看就简单了
密码修改失败请联系微信：mofangbao
回复 支持 反对

使用道具 举报

lqouc
发表于 2014-5-7 22:29:10 | 显示全部楼层
liyf 发表于 2014-5-7 22:15
检查了，后缀名没隐藏

我说的trim你试过了？你的字符声明是四位的，但是实际路径只有三位。
密码修改失败请联系微信：mofangbao
回复 支持 反对

使用道具 举报

下一页 »
12下一页
返回列表 发新帖
高级模式
B Color Image Link Quote Code Smilies
您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

小黑屋|帮助中心|清除痕迹

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表，不代表气象家园立场，禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表