[110] | 1 | ;+ |
---|
[232] | 2 | ; |
---|
[136] | 3 | ; @file_comments |
---|
[125] | 4 | ; compute the weight and address needed to interpolate data from |
---|
| 5 | ; an "irregular 2D grid" (defined as a grid made of quadrilateral cells) |
---|
| 6 | ; to any grid using the bilinear method |
---|
| 7 | ; |
---|
[226] | 8 | ; @categories |
---|
[157] | 9 | ; Interpolation |
---|
[110] | 10 | ; |
---|
[202] | 11 | ; @param olonin {in}{required}{type=2d array} |
---|
[136] | 12 | ; longitude of the input data |
---|
[110] | 13 | ; |
---|
[202] | 14 | ; @param olat {in}{required}{type=2d array} |
---|
[136] | 15 | ; latitude of the input data |
---|
| 16 | ; |
---|
[202] | 17 | ; @param omsk {in}{required}{type=2d array or -1} |
---|
[136] | 18 | ; land/sea mask of the input data |
---|
[202] | 19 | ; put -1 if input data are not masked |
---|
[136] | 20 | ; |
---|
[202] | 21 | ; @param alonin {in}{required}{type=2d array} |
---|
[136] | 22 | ; longitude of the output data |
---|
| 23 | ; |
---|
[202] | 24 | ; @param alat {in}{required}{type=2d array} |
---|
[136] | 25 | ; latitude of the output data |
---|
| 26 | ; |
---|
[202] | 27 | ; @param amsk {in}{required}{type=2d array or -1} |
---|
[136] | 28 | ; land/sea mask of the output data |
---|
[202] | 29 | ; put -1 if output data are not masked |
---|
[136] | 30 | ; |
---|
[202] | 31 | ; @param weig {out}{type=2d array} |
---|
| 32 | ; (see ADDR) |
---|
| 33 | ; |
---|
| 34 | ; @param addr {out}{type=2d array} |
---|
[118] | 35 | ; 2D arrays, weig and addr are the weight and addresses used to |
---|
| 36 | ; perform the interpolation: |
---|
| 37 | ; dataout = total(weig*datain[addr], 1) |
---|
| 38 | ; dataout = reform(dataout, jpia, jpja, /over) |
---|
[110] | 39 | ; |
---|
| 40 | ; @restrictions |
---|
[125] | 41 | ; - the input grid must be an "irregular 2D grid", defined as a grid made |
---|
[271] | 42 | ; of quadrilateral cells |
---|
[110] | 43 | ; - We supposed the data are located on a sphere, with a periodicity along |
---|
[125] | 44 | ; the longitude |
---|
[110] | 45 | ; - to perform the bilinear interpolation within quadrilateral cells, we |
---|
[125] | 46 | ; first morph the cell into a square cell and then compute the bilinear |
---|
| 47 | ; interpolation. |
---|
[295] | 48 | ; - if some corners of the cell are land points, their weights are set to 0 |
---|
[125] | 49 | ; and the weight is redistributed on the remaining "water" corners |
---|
[110] | 50 | ; - points located out of the southern and northern boundaries or in cells |
---|
[295] | 51 | ; containing only land points are set the same value as their closest neighbors |
---|
[125] | 52 | ; |
---|
[110] | 53 | ; @history |
---|
[118] | 54 | ; |
---|
[477] | 55 | ; - pinsardf 20120813T123457Z curie51.c-curie.tgcc.ccc.cea.fr (Linux) |
---|
| 56 | ; |
---|
| 57 | ; * less crypted control print |
---|
| 58 | ; |
---|
| 59 | ; - June 2006: Sebastien Masson (smasson\@lodyc.jussieu.fr) |
---|
| 60 | ; |
---|
[226] | 61 | ; @version |
---|
| 62 | ; $Id$ |
---|
[118] | 63 | ; |
---|
[110] | 64 | ;- |
---|
[495] | 65 | PRO compute_fromirr_bilinear_weigaddr, olonin, olat, omsk, alonin, alat, amsk $ |
---|
[509] | 66 | , weig, addr, lonref1 = lonref1, onplane = onplane, notuse_neighbor = notuse_neighbor |
---|
[110] | 67 | ; |
---|
[125] | 68 | compile_opt idl2, strictarrsubs |
---|
[110] | 69 | ; |
---|
[477] | 70 | s = size(SCOPE_TRACEBACK(/STRUCTURE),/DIMENSION) |
---|
| 71 | routine = (SCOPE_TRACEBACK(/STRUCTURE))[s-1].Routine |
---|
| 72 | ; |
---|
[110] | 73 | jpia = (size(alonin, /dimensions))[0] |
---|
| 74 | jpja = (size(alonin, /dimensions))[1] |
---|
| 75 | ; |
---|
| 76 | jpio = (size(olonin, /dimensions))[0] |
---|
[226] | 77 | jpjo = (size(olonin, /dimensions))[1] |
---|
[509] | 78 | ; |
---|
[202] | 79 | ; mask check |
---|
[257] | 80 | IF n_elements(omsk) EQ 1 AND omsk[0] EQ -1 THEN BEGIN |
---|
| 81 | omsk = replicate(1b, jpio, jpjo) |
---|
[271] | 82 | ; if this is ORCA2 grid... |
---|
| 83 | IF (jpio EQ 180 OR jpio EQ 182) AND (jpjo EQ 149 OR jpjo EQ 148 OR jpjo EQ 147 ) THEN BEGIN |
---|
| 84 | ; we look for ill defined cells. |
---|
| 85 | IF jpio EQ 182 THEN lontmp = olonin[1:180, *] ELSE lontmp = olonin |
---|
| 86 | ; longitudinal size of the cells... |
---|
| 87 | a = (lontmp-shift(lontmp, 1, 0)) |
---|
| 88 | d1 = 0.5 * max( [[[a]], [[360+a]]] MOD 360, dimension = 3) |
---|
| 89 | a = (lontmp-shift(lontmp, 1, 1)) |
---|
| 90 | d2 = 0.5 * max( [[[a]], [[360+a]]] MOD 360, dimension = 3) |
---|
| 91 | a = (shift(lontmp, 0, 1)-shift(lontmp, 1, 0)) |
---|
| 92 | d3 = 0.5 * max( [[[a]], [[360+a]]] MOD 360, dimension = 3) |
---|
| 93 | a = (shift(lontmp, 0, 1)-shift(lontmp, 1, 1)) |
---|
| 94 | d4 = 0.5 * max( [[[a]], [[360+a]]] MOD 360, dimension = 3) |
---|
| 95 | md = [[[d1]], [[d2]], [[d3]], [[d4]]] |
---|
| 96 | bad = max(md, dimension = 3) GE 1.5 |
---|
| 97 | bad = (bad + shift(bad, -1, -1) + shift(bad, 0, -1) + shift(bad, -1, 0)) < 1 |
---|
| 98 | IF jpio EQ 182 THEN BEGIN |
---|
| 99 | omsk[1:180, 80:120] = 1b - bad[*, 80:120] |
---|
| 100 | omsk[0, *] = omsk[180, *] |
---|
| 101 | omsk[181, *] = omsk[1, *] |
---|
| 102 | ENDIF ELSE omsk[*, 80:120] = 1b - bad[*, 80:120] |
---|
| 103 | ENDIF |
---|
[495] | 104 | ENDIF ELSE BEGIN |
---|
[271] | 105 | IF n_elements(omsk) NE jpio*jpjo THEN BEGIN |
---|
| 106 | ras = report('input grid mask do not have the good size') |
---|
| 107 | stop |
---|
| 108 | ENDIF |
---|
[495] | 109 | ENDELSE |
---|
| 110 | IF n_elements(amsk) EQ 1 AND amsk[0] EQ -1 THEN amsk = replicate(1b, jpia, jpja) ELSE BEGIN |
---|
[271] | 111 | IF n_elements(amsk) NE jpia*jpja THEN BEGIN |
---|
| 112 | ras = report('output grid mask do not have the good size') |
---|
| 113 | stop |
---|
| 114 | ENDIF |
---|
[495] | 115 | ENDELSE |
---|
[110] | 116 | ; |
---|
[509] | 117 | IF n_elements(lonref1) EQ 0 THEN lonref1 = 0. |
---|
| 118 | lonref2 = lonref1 + 360. |
---|
[110] | 119 | ; longitude, between 0 and 360 |
---|
[509] | 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. |
---|
[110] | 133 | ; |
---|
[297] | 134 | ; we work only on the output grid water points |
---|
[110] | 135 | awater = where(amsk EQ 1) |
---|
[125] | 136 | nawater = n_elements(awater) |
---|
[110] | 137 | ; |
---|
| 138 | ; define all cells that have corners located at olon, olat |
---|
| 139 | ; we define the cell with the address of all corners |
---|
| 140 | ; |
---|
| 141 | ; 3 2 |
---|
| 142 | ; +------+ |
---|
| 143 | ; | | |
---|
| 144 | ; | | |
---|
| 145 | ; | | |
---|
| 146 | ; +------+ |
---|
| 147 | ; 0 1 |
---|
| 148 | ; |
---|
| 149 | alladdr = lindgen(jpio, jpjo-1) |
---|
[303] | 150 | alladdrm1 = shift(alladdr, -1) |
---|
| 151 | IF ((jpio EQ 92 ) AND (jpjo EQ 76 OR jpjo EQ 75 OR jpjo EQ 74 )) OR $ |
---|
| 152 | ((jpio EQ 182 ) AND (jpjo EQ 149 OR jpjo EQ 148 OR jpjo EQ 147 )) OR $ |
---|
| 153 | ((jpio EQ 722 ) AND (jpjo EQ 522 OR jpjo EQ 521 OR jpjo EQ 520 )) OR $ |
---|
| 154 | ((jpio EQ 1442) AND (jpjo EQ 1021 OR jpjo EQ 1020 OR jpjo EQ 1019)) THEN alladdrm1[jpio-1, *] = alladdr[2, *] |
---|
[495] | 155 | |
---|
[110] | 156 | celladdr = lonarr(4, jpio*(jpjo-1)) |
---|
| 157 | celladdr[0, *] = alladdr |
---|
[303] | 158 | celladdr[1, *] = alladdrm1 |
---|
| 159 | celladdr[2, *] = temporary(alladdrm1) + jpio |
---|
| 160 | celladdr[3, *] = temporary(alladdr) + jpio |
---|
[110] | 161 | ; |
---|
| 162 | ; list the cells that have at least 1 ocean point as corner |
---|
| 163 | good = where(total(omsk[celladdr], 1) GT 0) |
---|
[125] | 164 | ; keep only those cells |
---|
[110] | 165 | celladdr = celladdr[*, temporary(good)] |
---|
| 166 | ; |
---|
| 167 | xcell = olon[celladdr] |
---|
| 168 | minxcell = min(xcell, dimension = 1, max = maxxcell) |
---|
| 169 | ycell = olat[celladdr] |
---|
| 170 | minycell = min(ycell, dimension = 1, max = maxycell) |
---|
| 171 | ; foraddr: address of the ocean water cell associated to each atmosphere water point |
---|
| 172 | foraddr = lonarr(nawater) |
---|
| 173 | ; forweight: x/y position of the atmosphere water point in the ocean water cell |
---|
| 174 | forweight = dblarr(nawater, 2) |
---|
| 175 | ; |
---|
| 176 | ; Loop on all the water point of the atmosphere |
---|
| 177 | ; We look for which ocean water cell contains the atmosphere water point |
---|
| 178 | ; |
---|
[303] | 179 | delta = max([(360./jpio), (180./jpjo)])* 3. |
---|
[509] | 180 | print, routine, ' total loop length = ' , nawater |
---|
[110] | 181 | FOR n = 0L, nawater-1 DO BEGIN |
---|
| 182 | ; control print |
---|
[477] | 183 | IF (n MOD 5000) EQ 0 THEN print, routine, ' n = ' , n |
---|
[121] | 184 | ; longitude and latitude of the atmosphere water point |
---|
[110] | 185 | xx = alon[awater[n]] |
---|
| 186 | yy = alat[awater[n]] |
---|
| 187 | ; 1) we reduce the number of ocean cells for which we will try to know if |
---|
[125] | 188 | ; xx,yy is inside. |
---|
[110] | 189 | CASE 1 OF |
---|
[295] | 190 | ; if we are near the north pole |
---|
[110] | 191 | yy GE (90-delta):BEGIN |
---|
| 192 | lat1 = 90-2*delta |
---|
| 193 | good = where(maxycell GE lat1) |
---|
| 194 | onsphere = 1 |
---|
| 195 | END |
---|
| 196 | ; if we are near the longitude periodicity area |
---|
[509] | 197 | xx LE (lonref1+delta) OR xx GE (lonref2-delta):BEGIN |
---|
[110] | 198 | lat1 = yy-delta |
---|
| 199 | lat2 = yy+delta |
---|
[509] | 200 | good = where((minxcell LE (lonref1+2*delta) OR maxxcell GE (lonref2-2*delta)) AND maxycell GE lat1 AND minycell LE lat2) |
---|
[110] | 201 | onsphere = 1 |
---|
| 202 | END |
---|
| 203 | ; other cases |
---|
| 204 | ELSE:BEGIN |
---|
| 205 | lon1 = xx-delta |
---|
| 206 | lon2 = xx+delta |
---|
| 207 | lat1 = yy-delta |
---|
| 208 | lat2 = yy+delta |
---|
| 209 | good = where(maxxcell GE lon1 AND minxcell LE lon2 AND maxycell GE lat1 AND minycell le lat2) |
---|
| 210 | ; ORCA cases : orca grid is irregular only northward of 40N |
---|
| 211 | CASE 1 OF |
---|
[303] | 212 | (jpio EQ 90 OR jpio EQ 92 ) AND (jpjo EQ 76 OR jpjo EQ 75 OR jpjo EQ 74 ):onsphere = yy GT 40 |
---|
| 213 | (jpio EQ 180 OR jpio EQ 182 ) AND (jpjo EQ 149 OR jpjo EQ 148 OR jpjo EQ 147 ):onsphere = yy GT 40 |
---|
| 214 | (jpio EQ 720 OR jpio EQ 722 ) AND (jpjo EQ 522 OR jpjo EQ 521 OR jpjo EQ 520 ):onsphere = yy GT 40 |
---|
[271] | 215 | (jpio EQ 1440 OR jpio EQ 1442) AND (jpjo EQ 1021 OR jpjo EQ 1020 OR jpjo EQ 1019):onsphere = yy GT 40 |
---|
[509] | 216 | ELSE:onsphere = 1 - keyword_set(onplane) |
---|
[110] | 217 | ENDCASE |
---|
| 218 | END |
---|
| 219 | ENDCASE |
---|
| 220 | ; we found a short list of possible ocean water cells containing the atmosphere water point |
---|
| 221 | IF good[0] NE -1 THEN BEGIN |
---|
| 222 | ; in which cell is located the atmosphere water point? |
---|
| 223 | ; Warning! inquad use clockwise quadrilateral definition |
---|
| 224 | ind = inquad(xx, yy $ |
---|
| 225 | , xcell[0, good], ycell[0, good] $ |
---|
| 226 | , xcell[3, good], ycell[3, good] $ |
---|
| 227 | , xcell[2, good], ycell[2, good] $ |
---|
| 228 | , xcell[1, good], ycell[1, good] $ |
---|
[303] | 229 | , onsphere = onsphere, newcoord = newcoord, /noprint, delta = delta, /double) |
---|
[110] | 230 | ; keep only the first cell (if the atmospheric point was located in several ocean cells) |
---|
| 231 | ind = ind[0] |
---|
| 232 | ; we found one ocean water cell containing the atmosphere water point |
---|
| 233 | IF ind NE -1 THEN BEGIN |
---|
| 234 | ind = good[ind] |
---|
| 235 | ; now, we morph the quadrilateral ocean cell into the reference square (0 -> 1) |
---|
[226] | 236 | ; in addition we get the coordinates of the atmospheric point in this new morphed square |
---|
[110] | 237 | IF onsphere THEN BEGIN |
---|
| 238 | ; Warning! quadrilateral2square use anticlockwise quadrilateral definition |
---|
| 239 | xy = quadrilateral2square(newcoord[0, 0], newcoord[1, 0] $ |
---|
| 240 | , newcoord[0, 3], newcoord[1, 3] $ |
---|
| 241 | , newcoord[0, 2], newcoord[1, 2] $ |
---|
| 242 | , newcoord[0, 1], newcoord[1, 1] $ |
---|
[282] | 243 | , newcoord[0, 4], newcoord[1, 4], /double) |
---|
[110] | 244 | ENDIF ELSE BEGIN |
---|
| 245 | xy = quadrilateral2square(xcell[0, ind], ycell[0, ind] $ |
---|
| 246 | , xcell[1, ind], ycell[1, ind] $ |
---|
| 247 | , xcell[2, ind], ycell[2, ind] $ |
---|
[282] | 248 | , xcell[3, ind], ycell[3, ind], xx, yy, /double) |
---|
[125] | 249 | ENDELSE |
---|
[110] | 250 | ; take care of rounding errors... |
---|
[282] | 251 | zero = where(abs(xy) LT 1e-6) |
---|
| 252 | IF zero[0] NE -1 THEN xy[zero] = 0.d |
---|
| 253 | one = where(abs(1.d - xy) LT 1e-6) |
---|
| 254 | IF one[0] NE -1 THEN xy[one] = 1.d |
---|
[202] | 255 | ; some checks... |
---|
| 256 | tmpmsk = omsk[celladdr[*, ind]] |
---|
| 257 | CASE 1 OF |
---|
| 258 | xy[0] LT 0 OR xy[0] GT 1 : stop |
---|
| 259 | xy[1] LT 0 OR xy[1] GT 1 : stop |
---|
| 260 | xy[0] EQ 0 AND xy[1] EQ 0 AND tmpmsk[0] EQ 0 : foraddr[n] = -1 |
---|
| 261 | xy[0] EQ 1 AND xy[1] EQ 0 AND tmpmsk[1] EQ 0 : foraddr[n] = -1 |
---|
| 262 | xy[0] EQ 1 AND xy[1] EQ 1 AND tmpmsk[2] EQ 0 : foraddr[n] = -1 |
---|
| 263 | xy[0] EQ 0 AND xy[1] EQ 1 AND tmpmsk[3] EQ 0 : foraddr[n] = -1 |
---|
| 264 | xy[0] EQ 0 AND (tmpmsk[0]+tmpmsk[3]) EQ 0 : foraddr[n] = -1 |
---|
| 265 | xy[0] EQ 1 AND (tmpmsk[1]+tmpmsk[2]) EQ 0 : foraddr[n] = -1 |
---|
| 266 | xy[1] EQ 0 AND (tmpmsk[0]+tmpmsk[1]) EQ 0 : foraddr[n] = -1 |
---|
| 267 | xy[1] EQ 1 AND (tmpmsk[2]+tmpmsk[3]) EQ 0 : foraddr[n] = -1 |
---|
[226] | 268 | ELSE: BEGIN |
---|
[202] | 269 | ; we keep its address |
---|
[271] | 270 | foraddr[n] = ind |
---|
[110] | 271 | ; keep the new coordinates |
---|
[202] | 272 | forweight[n, 0] = xy[0] |
---|
| 273 | forweight[n, 1] = xy[1] |
---|
| 274 | END |
---|
| 275 | ENDCASE |
---|
[110] | 276 | ENDIF ELSE foraddr[n] = -1 |
---|
[509] | 277 | ENDIF ELSE foraddr[n] = -1 |
---|
[110] | 278 | ENDFOR |
---|
[509] | 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 | |
---|
[110] | 287 | ; do we have some water atmospheric points that are not located in an water oceanic cell? |
---|
| 288 | bad = where(foraddr EQ -1) |
---|
| 289 | IF bad[0] NE -1 THEN BEGIN |
---|
| 290 | ; yes? |
---|
| 291 | ; we look for neighbor water atmospheric point located in water oceanic cell |
---|
| 292 | badaddr = awater[bad] |
---|
| 293 | good = where(foraddr NE -1) |
---|
| 294 | ; list the atmospheric points located in water oceanic cell |
---|
| 295 | goodaddr = awater[good] |
---|
| 296 | ; there longitude and latitude |
---|
| 297 | goodlon = alon[goodaddr] |
---|
| 298 | goodlat = alat[goodaddr] |
---|
| 299 | ; for all the bad points, look for a neighbor |
---|
| 300 | neig = lonarr(n_elements(bad)) |
---|
[303] | 301 | onsphere = 1 |
---|
[509] | 302 | print, routine, ' total loop length = ', n_elements(bad) |
---|
[303] | 303 | FOR i = 0L, n_elements(bad)-1L DO BEGIN |
---|
[477] | 304 | IF (i MOD 500) EQ 0 THEN print, routine, ' i = ', i |
---|
[303] | 305 | xtmp = alon[badaddr[i]] |
---|
| 306 | ytmp = alat[badaddr[i]] |
---|
| 307 | CASE 1 OF |
---|
| 308 | ; if we are near the north pole |
---|
| 309 | ytmp GE (90-delta):keep = where(goodlat GE 90-3*delta, cnt) |
---|
| 310 | ; if we are near the longitude periodicity area |
---|
[509] | 311 | xtmp LE (lonref1+delta) OR xtmp GE (lonref2-delta):keep = where((goodlon LE (lonref1+3*delta) OR goodlon GE (lonref2-3*delta)) $ |
---|
[303] | 312 | AND goodlat GE (ytmp-3*delta) AND goodlat LE (ytmp+3*delta), cnt) |
---|
| 313 | ; other cases |
---|
| 314 | ELSE:BEGIN |
---|
| 315 | keep = where(goodlon GE (xtmp-3*delta) AND goodlon LE (xtmp+3*delta) $ |
---|
| 316 | AND goodlat GE (ytmp-3*delta) AND goodlat le (ytmp+3*delta), cnt) |
---|
| 317 | ; ORCA cases : orca grid is irregular only northward of 40N |
---|
| 318 | CASE 1 OF |
---|
| 319 | (jpio EQ 90 OR jpio EQ 92 ) AND (jpjo EQ 76 OR jpjo EQ 75 OR jpjo EQ 74 ):onsphere = yy GT 40 |
---|
| 320 | (jpio EQ 180 OR jpio EQ 182 ) AND (jpjo EQ 149 OR jpjo EQ 148 OR jpjo EQ 147 ):onsphere = yy GT 40 |
---|
| 321 | (jpio EQ 720 OR jpio EQ 722 ) AND (jpjo EQ 522 OR jpjo EQ 521 OR jpjo EQ 520 ):onsphere = yy GT 40 |
---|
| 322 | (jpio EQ 1440 OR jpio EQ 1442) AND (jpjo EQ 1021 OR jpjo EQ 1020 OR jpjo EQ 1019):onsphere = yy GT 40 |
---|
| 323 | ELSE: |
---|
| 324 | ENDCASE |
---|
| 325 | END |
---|
| 326 | ENDCASE |
---|
[495] | 327 | IF cnt NE 0 THEN BEGIN |
---|
[303] | 328 | neig[i] = keep[(neighbor(xtmp, ytmp, goodlon[keep], goodlat[keep], sphere = onsphere))[0]] |
---|
| 329 | ENDIF ELSE neig[i] = (neighbor(alon[badaddr[i]], alat[badaddr[i]], goodlon, goodlat, sphere = onsphere))[0] |
---|
[110] | 330 | ENDFOR |
---|
| 331 | ; get the address regarding foraddr |
---|
| 332 | neig = good[neig] |
---|
| 333 | ; associate each bad point with its neighbor (get its address and weight) |
---|
| 334 | foraddr[bad] = foraddr[neig] |
---|
| 335 | forweight[bad, *] = forweight[neig, *] |
---|
[303] | 336 | ENDIF |
---|
[110] | 337 | ; transform the address of the ocean cell into the address of its 4 corners |
---|
| 338 | newaaddr = celladdr[*, temporary(foraddr)] |
---|
| 339 | ; now we compute the weight to give at each corner |
---|
| 340 | newaweig = dblarr(4, nawater) |
---|
| 341 | a = reform(forweight[*, 0], 1, nawater) |
---|
| 342 | b = reform(forweight[*, 1], 1, nawater) |
---|
[271] | 343 | forweight = -1 ; free memory |
---|
[309] | 344 | newaweig = [(1.d - a)*(1.d - b), (1.d - b)*a, a*b, (1.d - a)*b] |
---|
[271] | 345 | a = -1 & b = -1 ; free memory |
---|
[110] | 346 | ; mask the weight to suppress the corner located on land |
---|
[307] | 347 | newaweig = newaweig*(omsk[newaaddr]) |
---|
[110] | 348 | ; for cell with some land corner, |
---|
[282] | 349 | ; we have to redistribute the weight on the remaining water corners |
---|
[110] | 350 | ; weights normalization |
---|
[282] | 351 | totalweig = total(newaweig, 1, /double) |
---|
[309] | 352 | ;; IF min(totalweig, max = ma) LE 0.d then stop |
---|
| 353 | ;; IF ma- 1.d GT 1.e-6 then stop |
---|
[282] | 354 | newaweig = newaweig/(replicate(1.d, 4)#totalweig) |
---|
[110] | 355 | ; weights |
---|
| 356 | weig = dblarr(4, jpia*jpja) |
---|
| 357 | weig[*, awater] = temporary(newaweig) |
---|
| 358 | ; address |
---|
[307] | 359 | addr = lonarr(4, jpia*jpja) |
---|
[110] | 360 | addr[*, awater] = temporary(newaaddr) |
---|
| 361 | ; |
---|
| 362 | RETURN |
---|
| 363 | END |
---|