缺测值处理的一个小程序

呼啦呼啦

   我使用了NCEP的大气温度逐月资料，缺测值为-9.96921e+36，我编了一个小程序来计算每年的一个平均值，之前一直用的这个程序，今天突然发现这个缺测值的处理上（就是红色字体的地方），我是不是编的有问题啊？麻烦各位帮我看一下


program main
integer,parameter::  nx=105,ny=17,nt=384,mo=12,my=32
integer i,j,t,kk,iy,im
real tmsu(nx,ny,nt),tmsu2(mo,my,nx,ny),tmsu_ave(my,nx,ny),tmsu_sum(my,nx,ny)
open(1,file='D:\by\evalution_MSU_AMSU_BT\evalution\MSU_AMSU\AMSU\long_term\total_ocean\air_lev700_jp.grd',form='binary')
do t=1,nt
read(1) ((tmsu(i,j,t),i=1,nx),j=1,ny)
enddo
close(1)

do iy=1,my
do im=1,mo
kk=(iy-1)*12+im
do i=1,nx
do j=1,ny
tmsu2(im,iy,i,j)=tmsu(i,j,kk)
end do
end do
end do
end do

do iy=1,my
do i=1,nx
do j=1,ny
tmsu_sum(iy,i,j)=0.0
tmsu_ave(iy,i,j)=0.0
do im=1,mo
if(tmsu2(im,iy,i,j).ne.-9.96921e+36) then
tmsu_sum(iy,i,j)=tmsu_sum(iy,i,j)+tmsu2(im,iy,i,j)
else
tmsu_sum(iy,i,j)=-9.96921e+36
endif
enddo
if(tmsu_sum(iy,i,j).ne.-9.96921e+36)then
tmsu_ave(iy,i,j)=tmsu_sum(iy,i,j)/real(mo)
else
tmsu_ave(iy,i,j)=-9.96921e+36
endif
enddo
enddo
enddo

open(10,file='D:\by\evalution_MSU_AMSU_BT\evalution\MSU_AMSU\AMSU\long_term\total_ocean\air_lev700_aveyear.grd',form='binary')
do iy=1,my
write(10) ((tmsu_ave(iy,i,j),i=1,nx),j=1,ny)
enddo
close(10)
end

lqouc

唉，最近fortran版真是没人气啊，难得有个发帖的，不过楼主你这提问的水平真的值得提高哈，有空一定要看看“提问的智慧”，地址如下http://bbs.06climate.com/forum.php?mod=viewthread&tid=4571
现在说一下你提出的问题，简单来说就是想在平均的时候考虑缺测值是吧。但是你没说清你到底要怎样处理缺测值，是只要有一次缺测平均结果就为缺测？还是缺测值不参与计算。这样算法的设计是不同的，虽然都很简单。还有你说让大家给你看，虽然标出了问题代码，但是好歹描述一下你的问题，给点错误的结果出来啊。你来一句“我是不是编的有问题啊？”，这是问谁呢？
好了不啰嗦了，好在程序短小，我看完了，简单说下程序，缺测部分不管你想用哪种处理方法，这么写肯定是不对的了，我猜测你的思路应该是只要有一次缺测就赋缺测值。改动如下：
do iy=1,my
do i=1,nx
do j=1,ny
tmsu_ave=0.0; cc=0 ! 那个sum没用了，这里先声明一个整形计数器
do im=1,mo
if(tmsu2(im,iy,i,j).ne.-9.96921e+36) then
tmsu_ave(iy,i,j)=tmsu_ave(iy,i,j)+tmsu2(im,iy,i,j)
else
cc=cc+1
endif
enddo
if(cc.eq.0)then
tmsu_ave(iy,i,j)=tmsu_ave(iy,i,j)/real(mo)
else
tmsu_ave(iy,i,j)=-9.96921e+36
endif
enddo
enddo
enddo
你可以试一下，然后友情赠送一下你前面的问题，数组的读取（红色前面的部分）简直太麻烦了，给你精简下：
do iy=1,my
do im=1,mo
read(1) ((tmsu2(im,iy,i,j),i=1,nx),j=1,ny)   ! nt和tmsu(i,j,t) 这两个变量就没用了，省内存啊
enddo
enddo

嗯就这样吧

呼啦呼啦

lqouc 发表于 2015-3-24 15:56
唉，最近fortran版真是没人气啊，难得有个发帖的，不过楼主你这提问的水平真的值得提高哈，有空一定要看看 ...

谢谢版主！！以后改进帖子的提问方式。  还想问一下，像下面这样子写法来是不是就不用计数器了？
do im=1,mo
      do iy=1,my
      do i=1,nx
      do j=1,ny
      If(tmsu2(im,iy,i,j).ne.-9.99e+33.and.df(im,i,j).ne.-9.99e+33)then
        tmsu2(im,iy,i,j)=tmsu2(im,iy,i,j)/df(im,i,j)
      else
         tmsu2(im,iy,i,j)=-9.99e+33
      endif
      end do
end do
end do
end do

lqouc

呼啦呼啦 发表于 2015-3-25 22:16
谢谢版主！！以后改进帖子的提问方式。  还想问一下，像下面这样子写法来是不是就不用计数器了？
do im ...

额，没看出来你这个df哪里来的，你这个是用来完全替代之前红字部分？没看到累加过程啊。

呼啦呼啦

lqouc 发表于 2015-3-25 23:14
额，没看出来你这个df哪里来的，你这个是用来完全替代之前红字部分？没看到累加过程啊。

上面这一小段这是我编的一个类似程序里的一段，是用来求标准化距平的，完整的程序如下，里面也有很多对缺测的处理，我主要是想不考虑缺测，缺测不参与运算，还麻烦版主再帮我看一下，下面这个程序在处理缺测方面是否正确
    program main
      parameter(mo=12,my=32,nx=105,ny=16,mt=384)
      real tmsu2(mo,my,nx,ny),avf(mo,nx,ny),df(mo,nx,ny)
      real tmsu(nx,ny,mt),tmsu3(nx,ny,mt)
      integer cc

      open(1,file='D:\by\evalution_MSU_AMSU_BT\evalution\MSU_AMSU\AMSU\
     &long_term\total_ocean\lev700_TMT_1981_2012.grd',form='binary')
      do it=1,mt
      read(1) ((tmsu(i,j,it),i=1,nx),j=1,ny)
        end do
   
      do iy=1,my
      do im=1,mo
       kk=(iy-1)*12+im
       do i=1,nx
       do j=1,ny
      tmsu2(im,iy,i,j)=tmsu(i,j,kk)
        end do
        end do
        end do
        end do
     
      do im=1,mo
      do i=1,nx
      do j=1,ny
      avf(im,i,j)=0.0
        cc=0
      do iy=1,my
      if(tmsu2(im,iy,i,j).ne.-9999.000) then
      avf(im,i,j)=avf(im,i,j)+tmsu2(im,iy,i,j)
        else
      cc=cc+1
        endif
        end do
        if(cc.eq.0) then
        avf(im,i,j)=avf(im,i,j)/real(my)
        else
      avf(im,i,j)=-9999.000
        endif
        end do
        end do
        end do

      do im=1,mo
      do iy=1,my
      do i=1,nx
      do j=1,ny
      if(tmsu2(im,iy,i,j).ne.-9999.000.and.avf(im,i,j).ne.-9999.000)then
      tmsu2(im,iy,i,j)=tmsu2(im,iy,i,j)-avf(im,i,j)
      else
      tmsu2(im,iy,i,j)=-9999.000
      endif
        end do
        end do
        end do
        end do

       do  im=1,mo
       do  i=1,nx
       do  j=1,ny
       df(im,i,j)=0.0
        cc=0
       do iy=1,my
      if(tmsu2(im,iy,i,j).ne.-9999.000)THEN
        df(im,i,j)=df(im,i,j)+tmsu2(im,iy,i,j)*tmsu2(im,iy,i,j)
        else
          cc=cc+1
      endif
      end do
        if(cc.eq.0)then
           df(im,i,j)=SQRT(df(im,i,j)/real(my))
        else
         df(im,i,j)=-9999.000
        endif
      end do
        end do
        end do


      do im=1,mo
      do iy=1,my
      do i=1,nx
      do j=1,ny
      If(tmsu2(im,iy,i,j).ne.-9999.000.and.df(im,i,j).ne.-9999.000)then
        tmsu2(im,iy,i,j)=tmsu2(im,iy,i,j)/df(im,i,j)
      else
         tmsu2(im,iy,i,j)=-9999.000
      endif
      end do
        end do
        end do
        end do

      do  iy=1,my
      do  im=1,mo
      kk=(iy-1)*12+im
      do i=1,nx
      do j=1,ny
      tmsu3(i,j,kk)=tmsu2(im,iy,i,j)
      end do
        end do
        end do
        end do


      open(2,file='D:\by\evalution_MSU_AMSU_BT\evalution\MSU_AMSU\AMSU\
     &long_term\total_ocean\TMT_zhun_lev700.grd',form='binary')
       do it=1,mt
       write(2)((tmsu3(i,j,it),i=1,nx),j=1,ny)
       end do
      END

bree115

为什么我试了一下后，觉得赋初值 tmsu_ave=0.0的位置不太对？？

这只个名字

lqouc 发表于 2015-3-24 15:56
唉，最近fortran版真是没人气啊，难得有个发帖的，不过楼主你这提问的水平真的值得提高哈，有空一定要看看 ...

那请问如果想要缺测值不参与计算该怎么办呢？