NCL舍入小数的一个bug

peter19192009

1. ;******************************************************
   
2. ; check and modify if a 1-d array is succession
   
3. ;******************************************************
   
4. undef("check_succession")
   
5. function check_succession(x[*]:integer, n[1]:integer)
   
6. local number_x, index, number_0, i
   
7. begin
   
8. 
   
9.     number_x = dimsizes(x)
   
10.     index    = ind(x .eq. 0)
    
11.     number_0 = dimsizes(index)
    
12. 
    
13.     do i = 0,number_0-2
    
14.         if((index(i+1) - index(i)) .le. n) then
    
15.             x(index(i):index(i+1)) = 0
    
16.         end if
    
17.     end do
    
18. 
    
19. ; check 1 and end value
    
20.    
    
21.     if((x(0) .eq. 1) .and. (index(0) .lt. n)) then
    
22.         x(0:index(0)) = 0
    
23.     end if
    
24. 
    
25.     if((x(number_x-1) .eq. 1) .and. ((number_x-1) - index(number_0-1)) .lt. n ) then
    
26.         x(index(number_0-1):(number_x-1)) = 0
    
27.     end if
    
28. 
    
29.     return(x)
    
30. end
    
31. ;**********************************************************
    
32. ; longitude index B1D for a single day
    
33. ;**********************************************************
    
34. undef("TM_OneDay")
    
35. function TM_OneDay(x[*][*][*]:numeric, mid[1]:numeric, \
    
36.                    ext[1]:numeric, dn[1]:numeric, dm[1]:numeric, \
    
37.                    ds[1]:numeric, Emin[1]:numeric, Wmin[1]:numeric)
    
38. local Zn, Zm, Zs, dims, ntime, nlon, xblock, GHGS, GHGN, i, j
    
39. begin
    
40. 
    
41.     Zn = (mid + 1.5*ext) + dn
    
42.     Zm = (mid + 0.5*ext) + dm
    
43.     Zs = (mid - 0.5*ext) + ds
    
44.    
    
45.     dims  = dimsizes(x)
    
46.     ntime = dims(0)
    
47.     nlon  = dims(2)
    
48.    
    
49.     xblock = new((/ntime, nlon/), integer)
    
50.    
    
51.     do i = 0, ntime-1
    
52.         do j = 0, nlon-1
    
53.            GHGN = (x(i,{Zn},j) - x(i,{Zm},j)) / (Zn - Zm)
    
54.            GHGS = (x(i,{Zm},j) - x(i,{Zs},j)) / (Zm - Zs)
    
55.             if (GHGS .gt. Emin .and. GHGN .lt. Wmin) then
    
56.                 xblock(i,j) = 1
    
57.             else
    
58.                 xblock(i,j) = 0
    
59.             end if
    
60.         end do
    
61.     end do
    
62. 
    
63.     return(xblock)
    
64. end
    
65. ;*********************************************************
    
66. ;*************************************************************
    
67. load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
    
68. load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
    
69. load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
    
70. begin
    
71. 
    
72. ; parameters----
    
73. 
    
74.     Minlon = 285
    
75.     Maxlon = 345
    
76. 
    
77.     B_days = 4
    
78.     B_lons = 10
    
79. 
    
80.     mid  = 60
    
81.     ext  = 15
    
82.     dn   = 0
    
83.     dm   = 0
    
84.     ds   = 0
    
85.     Emin = 0
    
86.     Wmin = -5
    
87. 
    
88. ; read data------
    
89. 
    
90.     f      = addfile("./Z500_daily_deg1_1979_2013_JJA.nc", "r")
    
91.     geo    = short2flt(f->z(:,:,{Minlon:Maxlon}))
    
92.     geo    = geo / 9.8
    
93.     dims   = dimsizes(geo)
    
94.     Zg     = onedtond(ndtooned(geo), (/35,92,dims(1),dims(2)/))
    
95. 
    
96.     Zg!0    = "year"
    
97.     Zg&year = ispan(1979, 2013, 1)
    
98.     Zg!2    = "lat"
    
99.     Zg&lat  = geo&latitude
    
100.     Zg!3    = "lon"
     
101.     Zg&lon  = geo&longitude
     
102.    
     
103. ; calculate oneday block
     
104. 
     
105.     block = TM_OneDay(Zg({2012},:,:,:), mid, ext, dn, dm, ds, Emin, Wmin)
     
106.     write_matrix(block, "61i1", False)
     
107. 
     
108. ; check succession
     
109.    
     
110.     ndim  = dimsizes(block)
     
111.     Ntime = ndim(0)
     
112.     Nlons = ndim(1)
     
113. 
     
114.     do i = 0, Nlons-1
     
115.         block(:,i) = check_succession(block(:,i), B_days)
     
116.     end do
     
117.     write_matrix(block, "61i1", False)
     
118. 
     
119.     do i = 0, Ntime-1
     
120.         block(i,:) = check_succession(block(i,:), B_lons)
     
121.     end do
     
122.     write_matrix(block, "61i1", False)
     
123. 
     
124. ;----------------------------
     
125. end

复制代码

这里读入数据时候如果纬度反向geo    = short2flt(f->z(:,::-1,{Minlon:Maxlon}))运行结果会有很大差别，因为ext变量如果取奇数，存在舍入问题。

longlivehj

纬度反向和舍入有什么联系？bug 的话，用尽量简单的代码说明。在弄清你说的bug 前，可能已经被你的代码搞晕了。