【已解决】缺测值过多

赤炼修

program main
parameter(nx=866,nz=10,undef5=-9.9900000e+08)
real t(nx,nz),tgrad(nx,nz),dd(nz)
! character(3)::deepth(nz)=(/'000','005','010',050','075','100','150','200','250','300'/)
open(1,file='I:\1\deepth.txt')
read(1,*)(dd(k),k=1,nz)
close(1)
open(2,file='I:\1\tsprlatavek1-10.grd',form='binary')
read(2)((t(i,k),i=1,nx),k=1,nz)
close(2)

do i=1,nx
   do k=1,nz-1
    tgrad(i,k)=(t(i,k)-t(i,k+1))/(dd(k+1)-dd(k))
   enddo
enddo
do i=1,nx
  do k=1,nz-1
    if(t(i,k)==undef5.or.t(i,k+1)==undef5)then
          tgrad(i,k)=undef5
    endif
enddo;enddo

open(3,file='tgradlatvertical.grd',form='binary')
write(3)((tgrad(i,k),i=1,nx),k=1,nz)
close(3)

end

出图如下，图1中纵坐标是海洋深度，横坐标是经度，等值线为铅直温度梯度。


深度设置浅一点的图2，0~100米深度比300m深度缺测稍微少点


图3深度0~10m




图4是原始资料30°N，set t 1 10后显示的0~300m深度上的海温情况。

river

我觉得程序的逻辑是不是有些问题啊，改一下试试
do i=1,nx    
do k=1,nz-1  
if(t(i,k)==undef5.or.t(i,k+1)==undef5)then          
tgrad(i,k)=undef5    
else

tgrad(i,k)=(t(i,k)-t(i,k+1))/(dd(k+1)-dd(k))    
endif  

enddo
enddo
 
 

赤炼修

river 发表于 2013-10-4 19:49

                               
登录/注册后可看大图

我觉得程序的逻辑是不是有些问题啊，改一下试试
do i=1,nx    
do k=1,nz-1  

试过了图没有变化。另外读取原文件 tspr17avek1-10.grd发现有不少缺测，感觉这资料水平分辨率高但垂直层上较多缺测？准备试试对各个层次先进行区域平均变为水平上2.5*2.5,再来算铅直梯度试试。你觉得可行吗？还有，我计算铅直梯度的方法有问题吗？谢谢~~~

river

赤炼修 发表于 2013-10-5 09:33

                               
登录/注册后可看大图

试过了图没有变化。另外读取原文件 tspr17avek1-10.grd发现有不少缺测，感觉这资料水平分辨率高但垂直层上 ...

区域平均是可以的，但是缺测的地方还是要想办法处理一下，缺测值不能参加区域平均的计算，要不然结果有问题

homezhaow

囧 …… 问题还没解决吗？
建立数组 rlon（N）rlat（M）放每个格点的经纬，T（N，M，10），三维数组放温度, 以2.5乘2.5马赛克为例
do k=1,10
do i=1,N
do j=1,M
if （abs(T(int(rlon(i)/2.5)+1，int(rlon(i)/2.5)+1，k) )<10000）then
   rneed(int(rlon(i)/2.5)+1，int(rlon(i)/2.5)+1，k)=need(int(rlon(i)/2.5)+1，int(rlon(i)/2.5)+1，k)+T(int(rlon(i)/2.5)+1，int(rlon(i)/2.5)+1，k) )
   nsigma(int(rlon(i)/2.5)+1，int(rlon(i)/2.5)+1，k)=nsigma(int(rlon(i)/2.5)+1，int(rlon(i)/2.5)+1，k)+1
end;end;end

do k=1,10
do i=1,NN !马赛克后的lon格子数
do j=1,MM !马赛克后的lat格子数
   rneed(i,j,k)=rneed(i,j,k)/nsigma(i,j,k)
end;end;end

以上

吐槽下，我以为龚解决了呢 ……
潜水太多，冒个泡，带走一个积分 ……

赤炼修

homezhaow 发表于 2013-10-7 12:40

                               
登录/注册后可看大图

囧 …… 问题还没解决吗？
建立数组 rlon（N）rlat（M）放每个格点的经纬，T（N，M，10），三维数组放温度 ...

赞！多谢大神！