（已解决）NCL计算边界水汽

MacTaylor

我已经将纬向水汽通量做成了一个nc文件，时间44，纬度73，经度144，现在想做通过每个格点西边界的水汽，设想每个格点西边界的水汽定义是：假设格点A坐标为（1，2，3），那么通过格点B1（1，1，2）、B2（1，2，2）、B3（1，3，2）的水汽通量的平均值再乘以B1至B3的距离，我写的NCL脚本如下：

begin

  f = addfile("mf_ERA40.nc","r")
   system("/bin/rm "+"test.nc")
   fo      = addfile("test.nc","c")          ; create output file

  r = 6371000.00
  pi = 3.1415926
  d = r*abs((5.0-0.0)*pi/180.0)

  x = f->UQ
  xx =  month_to_season(x,"JJA")
  printVarSummary(xx)
  yy = xx(:,1:71,1:142)
   do i=1,71
     do j=1,142
      yy(:,i-1,j-1)=(xx(:,i-1,j-1)+xx(:,i,j-1)+xx(:,i+1,j-1))/3*d
     end do
   end do
  delete(xx)  
  delete(x)
  printVarSummary(yy)

  filedimdef(fo,"time",-1,True)
  fo->WQ   = yy

  end

计算出来了test.nc，但是现在有两个问题请教：
1、这样算出来的只能是从经度0开始的，那么我用grads画图，想设定set lon -40 150，其中-40至0是没有值的，这个在计算中应该怎么解决？
2、在计算中我定义的yy = xx(:,1:71,1:142)，那么数组yy的经纬度等于是从（0，-90）开始，但是根据设想yy的起始应该是（2.5，-87.5）开始，我这样定义是不是等于把得到的数组yy值的真实经纬度都向前移了一个格点？那么应该怎么定义呢？

求高手指点，谢谢各位！

mofangbao

楼主自己解决了吧？不妨说一下方法吧。
最近大家可能放假回家了，所以问题有时候会比较难解决

MacTaylor

mofangbao 发表于 2012-7-26 08:43

                               
登录/注册后可看大图

楼主自己解决了吧？不妨说一下方法吧。
最近大家可能放假回家了，所以问题有时候会比较难解决

恩，自己解决了，脚本如下：
begin

  f = addfile("mf_ERA40.nc","r")
   system("/bin/rm "+"WQ.nc")
   fo      = addfile("WQ.nc","c")          ; create output file

;  lat = f->lat
;  lon = f->lon
  r = 6371000.00
  pi = 3.1415926
  d = r*abs((5.0-0.0)*pi/180.0)

  x = f->UQ
  xx =  month_to_season(x,"JJA")
  printVarSummary(xx)
  yy = xx
  yy(:,0,:) = xx@_FillValue
  yy(:,72,:) = xx@_FillValue
      do i = 1,71,1
        do k = 0,43,1
          yy(k,i,0)=(xx(k,i-1,143)+xx(k,i,143)+xx(k,i+1,143))/3*d
        end do
      end do
   do j =1,143,1
       do i = 1,71,1
         do k = 0,43,1
           yy(k,i,j)=(xx(k,i-1,j-1)+xx(k,i,j-1)+xx(k,i+1,j-1))/3*d
         end do
       end do
   end do
  delete(xx)  
  delete(x)
  printVarSummary(yy)

  filedimdef(fo,"time",-1,True)
  fo->WQ   = yy

  end

同时请教清风，我之前将这段写成了这样：
yy = xx
  yy(:,0,:) = xx@_FillValue
  yy(:,72,:) = xx@_FillValue
   do j = 1,143,1
    if (lon(j).eq.0.0) then
      do i = 1,71,1
        do k = 0,43,1
          yy(k,i,j-1)=(xx(k,i-1,143)+xx(k,i,143)+xx(k,i+1,143))/3*d
        end do
      end do
    else
       do i = 1,71,1
         do k = 0,43,1
           yy(k,i,j)=(xx(k,i-1,j-1)+xx(k,i,j-1)+xx(k,i+1,j-1))/3*d
         end do
       end do
    end if
   end do
得到的nc文件用grads画图和现在的基本一致，但是在0经度的地方有个缺口，像拼接的，这是为什么呢？

timevigor

  d = r*abs((5.0-0.0)*pi/180.0)请问楼主5.0-0.0是怎么算的？

zongqing123

{:eb502:}{:eb502:}{:eb502:}

zjchris

是算东西南北边界的水汽收支吗

你好麦克风

timevigor 发表于 2015-10-30 19:25
  d = r*abs((5.0-0.0)*pi/180.0)请问楼主5.0-0.0是怎么算的？

应该是数据的分辨率