- Timestamp:
- 12/19/19 10:27:05 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/SRC/Interpolation/compute_fromirr_bilinear_weigaddr.pro
r495 r509 64 64 ;- 65 65 PRO compute_fromirr_bilinear_weigaddr, olonin, olat, omsk, alonin, alat, amsk $ 66 , weig, addr 66 , weig, addr, lonref1 = lonref1, onplane = onplane, notuse_neighbor = notuse_neighbor 67 67 ; 68 68 compile_opt idl2, strictarrsubs … … 76 76 jpio = (size(olonin, /dimensions))[0] 77 77 jpjo = (size(olonin, /dimensions))[1] 78 ; 78 79 ; mask check 79 80 IF n_elements(omsk) EQ 1 AND omsk[0] EQ -1 THEN BEGIN … … 114 115 ENDELSE 115 116 ; 117 IF n_elements(lonref1) EQ 0 THEN lonref1 = 0. 118 lonref2 = lonref1 + 360. 116 119 ; longitude, between 0 and 360 117 alon = alonin MOD 360 118 out = where(alon LT 0) 119 IF out[0] NE -1 THEN alon[out] = alon[out]+360 120 olon = olonin MOD 360 121 out = where(olon LT 0) 122 IF out[0] NE -1 THEN olon[out] = olon[out]+360 120 alon = (alonin + 720. ) MOD 360. 121 ; longitude, between lonref1 and lonref2 122 out = where(alon LT lonref1) 123 IF out[0] NE -1 THEN alon[out] = alon[out] + 360. 124 out = where(alon GT lonref2) 125 IF out[0] NE -1 THEN alon[out] = alon[out] - 360. 126 ; longitude, between 0 and 360 127 olon = ( olonin + 720. ) MOD 360. 128 ; longitude, between lonref1 and lonref2 129 out = where(olon LT lonref1) 130 IF out[0] NE -1 THEN olon[out] = olon[out] + 360. 131 out = where(olon GT lonref2) 132 IF out[0] NE -1 THEN olon[out] = olon[out] - 360. 123 133 ; 124 134 ; we work only on the output grid water points … … 168 178 ; 169 179 delta = max([(360./jpio), (180./jpjo)])* 3. 180 print, routine, ' total loop length = ' , nawater 170 181 FOR n = 0L, nawater-1 DO BEGIN 171 182 ; control print … … 184 195 END 185 196 ; if we are near the longitude periodicity area 186 xx LE delta OR xx GE (360-delta):BEGIN197 xx LE (lonref1+delta) OR xx GE (lonref2-delta):BEGIN 187 198 lat1 = yy-delta 188 199 lat2 = yy+delta 189 good = where((minxcell LE 2*delta OR maxxcell GE (360-2*delta)) AND maxycell GE lat1 AND minycell LE lat2)200 good = where((minxcell LE (lonref1+2*delta) OR maxxcell GE (lonref2-2*delta)) AND maxycell GE lat1 AND minycell LE lat2) 190 201 onsphere = 1 191 202 END … … 203 214 (jpio EQ 720 OR jpio EQ 722 ) AND (jpjo EQ 522 OR jpjo EQ 521 OR jpjo EQ 520 ):onsphere = yy GT 40 204 215 (jpio EQ 1440 OR jpio EQ 1442) AND (jpjo EQ 1021 OR jpjo EQ 1020 OR jpjo EQ 1019):onsphere = yy GT 40 205 ELSE:onsphere = 1 216 ELSE:onsphere = 1 - keyword_set(onplane) 206 217 ENDCASE 207 218 END … … 264 275 ENDCASE 265 276 ENDIF ELSE foraddr[n] = -1 266 ENDIF ELSE foraddr[n] = -1277 ENDIF ELSE foraddr[n] = -1 267 278 ENDFOR 279 280 IF keyword_set(notuse_neighbor) THEN BEGIN ; don't do anything for bad points -> they will have no addr/weig like land points 281 ok = where(foraddr NE -1, nawater) ; keep only the points that are above a water oceanic cell. 282 foraddr = foraddr[ok] 283 forweight = forweight[ok, *] 284 awater = awater[temporary(ok)] 285 ENDIF 286 268 287 ; do we have some water atmospheric points that are not located in an water oceanic cell? 269 288 bad = where(foraddr EQ -1) … … 281 300 neig = lonarr(n_elements(bad)) 282 301 onsphere = 1 302 print, routine, ' total loop length = ', n_elements(bad) 283 303 FOR i = 0L, n_elements(bad)-1L DO BEGIN 284 304 IF (i MOD 500) EQ 0 THEN print, routine, ' i = ', i … … 289 309 ytmp GE (90-delta):keep = where(goodlat GE 90-3*delta, cnt) 290 310 ; if we are near the longitude periodicity area 291 xtmp LE delta OR xtmp GE (360-delta):keep = where((goodlon LE 3*delta OR goodlon GE (360-3*delta)) $311 xtmp LE (lonref1+delta) OR xtmp GE (lonref2-delta):keep = where((goodlon LE (lonref1+3*delta) OR goodlon GE (lonref2-3*delta)) $ 292 312 AND goodlat GE (ytmp-3*delta) AND goodlat LE (ytmp+3*delta), cnt) 293 313 ; other cases
Note: See TracChangeset
for help on using the changeset viewer.