我写了个 将等sigma面数据插值到等压面 程序感觉有些问题 ，求助各位大神

云大小子

!注意的是sigma层的数据读取k（1-23），从模式顶开始往下存取数据的
      program sigmap
      implicit none
      integer:: ix,jy,np,it,i,j,n,kx,m,t,irec
      parameter(ix=147,jy=79,kx=23,m=1,it=28)
      parameter(np=17)
      real    plev(np),sig(kx)
      data    plev/ 100000.0 ,92500.0 ,85000.0 ,70000.0 ,60000.0
     $                ,50000.0 ,40000.0 ,30000.0 , 25000.0, 20000.0 ,15000.0 ,
     $         10000.0 ,7000.0 ,5000.0 ,3000.0, 2000.0, 1000.0/
      
        
        data    sig/ 2.5000000E-02 , 7.5000003E-02 ,0.1250000 , 0.1750000
     $        , 0.2250000   
     $  ,0.2750000 , 0.3250000 , 0.3750000 , 0.4250000 , 0.4750000   
     $  ,0.5250000 , 0.5750000 , 0.6250000 , 0.6750000 , 0.7250000   
     $  ,0.7750000 , 0.8250000 , 0.8700000 , 0.9100000 , 0.9450000   
     $   ,0.9700000 , 0.9850000 , 0.9950000/      
           integer mrec,k,y,s                           
          real       fin(ix,jy)
        real       fout(ix,jy)                           
          real       bc(ix,jy,kx,m,it)
          real       bcp(ix,jy,np,m,it)
          real       pstar(147,79)
          irec=0
         mrec=0
          open(55,file='F:\riems2.01pom-gorcart_online\pbc\radlw2.dat',
     $          form='binary')
           open(50,file='out.dat',form='binary')
     
           open (60,file='pstartxt',form='formatted')
            
                 do j=1,79
               do i=1,147
                     
                          read(60,*) pstar(i,j)
               
                   enddo
             enddo
           close(60)
        
        do t=1,28
        do n=1,1
        do k=1,23
        do j=1,79
        do i=1,147            
            read(55)bc(i,j,k,n,t)
      
            
        enddo
        enddo
        enddo
        enddo
        enddo   
           print*,bc(10,20,1,1,22)        
        print*,bc(10,20,2,1,22)
                   call intlin(bcp,bc,pstar,sig,ix,jy,kx,plev,np)
         print*,bcp(10,20,1,1,22)
           print*,bcp(10,20,2,1,22)
           
           do t=1,it
                do n=1,m
                do k=1,np
            do j=1,79
              do i=1,147
    !              fout(i,j)=bcp(i,j,k)
   
                         write(50) bcp(i,j,k,n,t)
                    
                enddo
            enddo
                enddo
                enddo
            enddo
          end        
                  SUBROUTINE INTLIN(FP,F,PSTAR,SIG,IM,JM,KM,P,KP)
      implicit none
      INTEGER IM,JM,KM,KP,it,m
      parameter (it=28,m=1)
        REAL    FP(IM,JM,KP,m,it),F(IM,JM,KM,m,it)
      REAL    PSTAR(IM,JM)
      REAL    SIG(KM),P(KP)
      REAL    PTOP,RGAS,GRAV,BLTOP,TLAPSE
      COMMON /CONST/ PTOP,RGAS,GRAV,BLTOP,TLAPSE
      INTEGER I,J,K,N,nn,t
      INTEGER K1,K1P
      REAL    SIGP,WP,W1
      
        print*,kp !17
        print*,km !23
C
C  INTLIN IS FOR VERTICAL INTERPOLATION OF U, V, AND RELATIVE HUMIDITY.
C        THE INTERPOLATION IS LINEAR IN P.  WHERE EXTRAPOLATION IS
C        NECESSARY, FIELDS ARE CONSIDERED TO HAVE 0 VERTICAL DERIVATIVE.
c        sig(1)-->sig(23),0-->1
c     
        do t=1,it
        do nn=1,m
        DO J=1,JM
        DO I=1,IM
          DO N=1,KP
            SIGP = (P(N)-PTOP) / (PSTAR(I,J)-PTOP)
            K1=0
            DO K=1,KM
              IF (SIGP.GT.SIG(K)) K1=K
            ENDDO
            IF(SIGP.LE.SIG(1)) THEN
              FP(I,J,N,nn,t) = F(I,J,1,nn,t)
            ELSE IF((SIGP.GT.SIG(1)).AND.(SIGP.LT.SIG(KM))) THEN
              K1P = K1 + 1
              WP  = (SIGP-SIG(K1))/(SIG(K1P)-SIG(K1))
              W1  = 1.-WP
              FP(I,J,N,nn,t)  = W1*F(I,J,K1,nn,t)+WP*F(I,J,K1P,nn,t)
            ELSE IF(SIGP.GE.SIG(KM)) THEN
              FP(I,J,N,nn,t)  = F(I,J,Km,nn,t)
             print*,f(i,j,km,nn,t)
               print*,n
                        print*,fp(i,j,n,nn,t)
         
                  ENDIF
          ENDDO
        ENDDO
      ENDDO
      enddo
        enddo
        RETURN
      END
                        
这个程序 有问题 差出的结果不合适 谁有慧眼 能看出问题么 问题应该出在SUBROUTINE INTLIN(FP,F,PSTAR,SIG,IM,JM,KM,P,KP)

易小凯

具体的问题是什么呢？

言深深

这个真心没法看，自己再琢磨琢磨吧，可能你的插值子程序不对

云大小子

易小凯 发表于 2012-10-29 20:54

                               
登录/注册后可看大图

具体的问题是什么呢？

插值出来的结果 28时刻没有数据

云大小子

易小凯 发表于 2012-10-29 20:54

                               
登录/注册后可看大图

具体的问题是什么呢？

然后对比同一时刻 同一层的 数据画出的图形感觉 差别比较远

云大小子

云大小子 发表于 2012-10-30 12:10

                               
登录/注册后可看大图

然后对比同一时刻 同一层的 数据画出的图形感觉 差别比较远

上面的程序没有 任何问题了，现在我确认一下，以前是我的ctl文件没写对的效果

pawpaw丁

你好，最近在用MM5模式 把MM5输出转换为sedris格式，需要将sigma层转为高度 查资料看到这些公式：
Z=-R*A/2/G*(ln(P0/P00))**2-R*TS0/G*(ln(P0/P00))
P0=PS0*σ+PTOP
PS0=P0(SURFACE)-PTOP
P0(SURFACE)计算是按照什么呢？
希望对能提供建议

半个瞎子

楼主你好，请问你这个程序里 PTOP取的什么值呢？是10hPa么？ 还有格点上的PSTAR不应该是一样的么？“ COMMON /CONST/ PTOP,RGAS,GRAV,BLTOP,TLAPSE”这句什么意思呢？ 请楼主能够帮忙解答一下，谢谢了。