我用Fortran计算EC细网格数据的TD，但写法有误？请高手指教！

deepknight

我之前用程序1写的，但计算TD数值不对；
后来改成了程序2，计算出的TD才靠谱。
不知道这两种写法有什么区别？或者说前面的一种写法错误在哪里？我是Fortran菜鸟，还请高手指教~

程序1：
        PROGRAM TD_CALCULATION
        PARAMETER (LEV=7,NX=361,NY=281)
        DIMENSION T(NX,NY,LEV),RH(NX,NY,LEV),TD(NX,NY,LEV)

      OPEN(100,FILE='t.txt')
        OPEN(200,FILE='rh.txt')
        OPEN(300,FILE='td.txt')

        DO k=1,LEV
      READ(100,*) ((t(i,j,k),j=1,361),i=1,281)
        READ(200,*) ((rh(i,j,k),j=1,361),i=1,281)
        ENDDO
        
        WRITE(*,*)RH(1,1,1)
        
        CALL CAL_RH_TD(T,RH,TD)
        DO k=1,7
      WRITE(300,'(361f12.4)') ((td(i,j,k),j=1,361),i=1,281)
        ENDDO

        END

         SUBROUTINE CAL_RH_TD(T,RH,TD)
C     THIS SUBROUTINE CALCULATE DEW POINT BY USING TEMPERATURE AND RELATIVE HUMIDITY
C     INPUT   T ========> TEMPERATURE (CELSIUS DEGREE)
C             RH =======> RELATIVE HUMIDITY  (E/ES)
C             LEV ======> LEVEL
C     OUTPUT  TD =======> DEW POINT  (CELSIUS DEGREE)
      DIMENSION T(281,361,7),RH(281,361,7),TD(281,361,7)
      DO K=1,7
        DO I=1,281
        DO J=1,361
            RH(i,j,k)=RH(i,j,k)/100.
            IF(RH(i,j,k).GT.1.0)RH(i,j,k)=1.0
            IF(RH(i,j,k).LT.0.05)RH(i,j,k)=0.05
          RMA=ALOG(RH(i,j,k))*461.*273.15/4.18/597./1000.+T(i,j,k)/
     1        (273.15+T(i,j,k))
          TD(i,j,k)=273.15*RMA/(1-RMA)
            
      ENDDO
        ENDDO
        ENDDO

      RETURN
      END



程序2：
        PROGRAM TD_CALCULATION
        parameter(ki=281,kj=361,kc=7)
        integer i,j,c
        DIMENSION T(ki,kj,kc),RH(ki,kj,kc),TD(ki,kj,kc)

      OPEN(100,FILE='t.txt')
        OPEN(200,FILE='rh.txt')
        OPEN(300,FILE='td.txt')

        do  c=1,kc
         READ(100,*)((t(i,j,c),j=1,kj),i=1,ki)
      READ(200,*)((rh(i,j,c),j=1,kj),i=1,ki)
         
      enddo

          CALL CAL_RH_TD(T,RH,TD)

        do  c=1,kc
           WRITE(300,'(361f12.4)')((td(i,j,c),j=1,kj),i=1,ki)
            ENDDO

        WRITE(*,'(3f12.4)')        T(100,200,2),RH(100,200,2),TD(100,200,2)

        END

         SUBROUTINE CAL_RH_TD(T,RH,TD)
C     THIS SUBROUTINE CALCULATE DEW POINT BY USING TEMPERATURE AND RELATIVE HUMIDITY
C     INPUT   T ========> TEMPERATURE (CELSIUS DEGREE)
C             RH =======> RELATIVE HUMIDITY  (E/ES)
C             LEV ======> LEVEL
C     OUTPUT  TD =======> DEW POINT  (CELSIUS DEGREE)
        parameter(ki=281,kj=361,kc=7)
      DIMENSION T(ki,kj,kc),RH(ki,kj,kc),TD(ki,kj,kc)
      DO c=1,7
        DO i=1,281
        DO j=1,361
            RH(i,j,c)=RH(i,j,c)/100.
            IF(RH(i,j,c).GT.1.0)RH(i,j,c)=1.0
            IF(RH(i,j,c).LT.0.05)RH(i,j,c)=0.05
          RMA=ALOG(RH(i,j,c))*461.*273.15/4.18/597./1000.+T(i,j,c)/
     1        (273.15+T(i,j,c))
          TD(i,j,c)=273.15*RMA/(1-RMA)  
      ENDDO
        ENDDO
        ENDDO
        
      RETURN
      END

kongfeng0824

感谢楼主分享！！！

kongfeng0824

感谢楼主分享！！！

Cheen

第一个程序在调用子程序的时候数据传递出错了，主程序数组定义：
PARAMETER (LEV=7,NX=361,NY=281)
DIMENSION T(NX,NY,LEV),RH(NX,NY,LEV),TD(NX,NY,LEV) ！T(361,281,7)
子程序数组定义DIMENSION T(281,361,7)……
数组传递的时候 X Y 乱了

deepknight

Cheen 发表于 2013-4-22 21:23

                               
登录/注册后可看大图

第一个程序在调用子程序的时候数据传递出错了，主程序数组定义：
PARAMETER (LEV=7,NX=361,NY=281)
DIMEN ...

谢谢Cheen帮我发现问题！
原来是我太粗心了，犯了这么低级的错误。。。
我把x,y方向数组的范围写错了，看来还是写成ki,ki这样比较简单些。

clare

请教lz，ec细网格数据是怎么得到并处理为*.txt各要素场的？

deepknight

要先用grads读nc数据输出为二进制的dat文件，然后再用fortran读取dat文件再转成txt文件~

346883050

请问楼主，cmacast下发的ec细网格二级压缩资料*.bin怎么处理成.txt文件呢？

deepknight

346883050 发表于 2015-4-30 10:28
请问楼主，cmacast下发的ec细网格二级压缩资料*.bin怎么处理成.txt文件呢？

很抱歉，这个cmacast下发的文件处理我没研究过哈~我只做了.NC格式文件的处理...

deepknight

clare 发表于 2015-1-3 11:56
请教lz，ec细网格数据是怎么得到并处理为*.txt各要素场的？

要先用grads读nc数据输出为二进制的dat文件，然后再用fortran读取dat文件再转成txt文件~