| 1 | ;+ |
|---|
| 2 | ; NAME: compute_fromreg_imoms3_weigaddr |
|---|
| 3 | ; |
|---|
| 4 | ; PURPOSE: compute the weight and address neede to interpolate data from a |
|---|
| 5 | ; "regular grid" to any grid using the imoms3 method |
|---|
| 6 | ; |
|---|
| 7 | ; CATEGORY:interpolation |
|---|
| 8 | ; |
|---|
| 9 | ; CALLING SEQUENCE: |
|---|
| 10 | ; compute_fromreg_imoms3_weigaddr, alon, alat, olon, olat, weig, addr |
|---|
| 11 | ; |
|---|
| 12 | ; INPUTS: |
|---|
| 13 | ; lonin and latin: longitude/latitude of the input data |
|---|
| 14 | ; lonout and latout: longitude/latitude of the output data |
|---|
| 15 | ; |
|---|
| 16 | ; KEYWORD PARAMETERS: |
|---|
| 17 | ; |
|---|
| 18 | ; /NONORTHERNLINE and /NOSOUTHERNLINE: activate if you don't whant to take into |
|---|
| 19 | ; account the northen/southern line of the input data when perfoming the |
|---|
| 20 | ; interpolation. |
|---|
| 21 | ; |
|---|
| 22 | ; OUTPUTS: |
|---|
| 23 | ; weig, addr: 2D arrays, weig and addr are the weight and addresses used to |
|---|
| 24 | ; perform the interpolation: |
|---|
| 25 | ; dataout = total(weig*datain[addr], 1) |
|---|
| 26 | ; dataout = reform(dataout, jpio, jpjo, /over) |
|---|
| 27 | ; |
|---|
| 28 | ; COMMON BLOCKS: none |
|---|
| 29 | ; |
|---|
| 30 | ; SIDE EFFECTS: ? |
|---|
| 31 | ; |
|---|
| 32 | ; RESTRICTIONS: |
|---|
| 33 | ; - the input grid must be a "regular grid", defined as a grid for which |
|---|
| 34 | ; lontitudes and latitudes are regularly spaced. |
|---|
| 35 | ; - We supposed the data are located on a sphere, with a periodicity along |
|---|
| 36 | ; the longitude. |
|---|
| 37 | ; - points located between the first/last 2 lines are interpolated |
|---|
| 38 | ; using a imoms3 interpolation along the longitudinal direction and linear |
|---|
| 39 | ; interpolation along the latitudinal direction |
|---|
| 40 | ; - points located out of the southern and northern boundaries are interpolated |
|---|
| 41 | ; using a imoms3 interpolation only along the longitudinal direction. |
|---|
| 42 | ; |
|---|
| 43 | ; EXAMPLE: |
|---|
| 44 | ; |
|---|
| 45 | ; MODIFICATION HISTORY: |
|---|
| 46 | ; November 2005: Sebastien Masson (smasson@lodyc.jussieu.fr) |
|---|
| 47 | ; |
|---|
| 48 | ;- |
|---|
| 49 | ; |
|---|
| 50 | ;---------------------------------------------------------- |
|---|
| 51 | ;---------------------------------------------------------- |
|---|
| 52 | ; |
|---|
| 53 | PRO compute_fromreg_imoms3_weigaddr, alonin, alatin, olonin, olat, weig, addr $ |
|---|
| 54 | , NONORTHERNLINE = nonorthernline, NOSOUTHERNLINE = nosouthernline, FORCE = force |
|---|
| 55 | ; |
|---|
| 56 | compile_opt strictarr, strictarrsubs |
|---|
| 57 | ; |
|---|
| 58 | alon = alonin |
|---|
| 59 | alat = alatin |
|---|
| 60 | olon = olonin |
|---|
| 61 | ; |
|---|
| 62 | jpia = n_elements(alon) |
|---|
| 63 | jpja = n_elements(alat) |
|---|
| 64 | ; |
|---|
| 65 | jpio = (size(olon, /dimensions))[0] |
|---|
| 66 | jpjo = (size(olon, /dimensions))[1] |
|---|
| 67 | ; |
|---|
| 68 | ; alon |
|---|
| 69 | minalon = min(alon, max = maxalon) |
|---|
| 70 | IF maxalon-minalon GE 360. THEN stop |
|---|
| 71 | ; alon must be monotonically increasing |
|---|
| 72 | IF array_equal(sort(alon), lindgen(jpia)) NE 1 THEN BEGIN |
|---|
| 73 | shiftx = -(where(alon EQ min(alon)))[0] |
|---|
| 74 | alon = shift(alon, shiftx) |
|---|
| 75 | IF array_equal(sort(alon), lindgen(jpia)) NE 1 THEN stop |
|---|
| 76 | ENDIF ELSE shiftx = 0 |
|---|
| 77 | ; alon must be regularly spaced |
|---|
| 78 | step = alon-shift(alon, 1) |
|---|
| 79 | step[0] = step[0] + 360. |
|---|
| 80 | IF total((step-step[0]) GE 1.e-6) NE 0 THEN BEGIN |
|---|
| 81 | print, 'input longitude must be regularly spaced, we stop...' |
|---|
| 82 | stop |
|---|
| 83 | ENDIF |
|---|
| 84 | ; we extend the longitude range of alon (-> easy interpolation even |
|---|
| 85 | ; near minalon et maxalon) |
|---|
| 86 | toadd = 10*jpia/360+1 |
|---|
| 87 | alon = [alon[jpia-toadd:jpia-1]-360., alon[*], alon[0:toadd-1]+360.] |
|---|
| 88 | jpia = jpia+2*toadd |
|---|
| 89 | ; alat |
|---|
| 90 | revy = alat[0] GT alat[1] |
|---|
| 91 | IF revy THEN alat = reverse(alat) |
|---|
| 92 | ; alat must be monotonically increasing |
|---|
| 93 | IF array_equal(sort(alat), lindgen(jpja)) NE 1 THEN stop |
|---|
| 94 | ; alat must be regularly spaced |
|---|
| 95 | step = alat-shift(alat, 1) |
|---|
| 96 | step = step[1:jpja - 1L] |
|---|
| 97 | IF total((step-step[0]) GE 1.e-6) NE 0 AND NOT keyword_set(force) THEN BEGIN |
|---|
| 98 | print, 'input latitude must be regularly spaced, we stop...' |
|---|
| 99 | stop |
|---|
| 100 | ENDIF |
|---|
| 101 | ; |
|---|
| 102 | if keyword_set(nonorthernline) then BEGIN |
|---|
| 103 | jpja = jpja - 1L |
|---|
| 104 | alat = alat[0: jpja-1L] |
|---|
| 105 | ENDIF |
|---|
| 106 | if keyword_set(nosouthernline) then BEGIN |
|---|
| 107 | alat = alat[1: jpja-1L] |
|---|
| 108 | jpja = jpja - 1L |
|---|
| 109 | ENDIF |
|---|
| 110 | ; olon between minalon et minalon+360 |
|---|
| 111 | out = where(olon LT minalon) |
|---|
| 112 | WHILE out[0] NE -1 DO BEGIN |
|---|
| 113 | olon[out] = olon[out]+360. |
|---|
| 114 | out = where(olon LT minalon) |
|---|
| 115 | ENDWHILE |
|---|
| 116 | out = where(olon GE minalon+360.) |
|---|
| 117 | WHILE out[0] NE -1 DO BEGIN |
|---|
| 118 | olon[out] = olon[out]- 360. |
|---|
| 119 | out = where(olon GE minalon+360.) |
|---|
| 120 | ENDWHILE |
|---|
| 121 | ; make sure that all values of olon are located within values of alon |
|---|
| 122 | IF min(olon, max = ma) LT minalon THEN stop |
|---|
| 123 | IF ma GE minalon+360. THEN stop |
|---|
| 124 | ; |
|---|
| 125 | xaddr = lonarr(16, jpio*jpjo) |
|---|
| 126 | yaddr = lonarr(16, jpio*jpjo) |
|---|
| 127 | weig = fltarr(16, jpio*jpjo) |
|---|
| 128 | ; |
|---|
| 129 | indexlon = value_locate(alon, olon) |
|---|
| 130 | IF total(alon[indexlon] GT olon) NE 0 THEN stop |
|---|
| 131 | IF total(alon[indexlon + 1L] LE olon) NE 0 THEN stop |
|---|
| 132 | IF (where(indexlon LE 1L ))[0] NE -1 THEN stop |
|---|
| 133 | IF (where(indexlon GE jpia-3L))[0] NE -1 THEN stop |
|---|
| 134 | indexlat = value_locate(alat, olat) |
|---|
| 135 | ; |
|---|
| 136 | ; for the ocean points located below the atm line |
|---|
| 137 | ; jpja-2 and above the line 1 |
|---|
| 138 | ; for those points we can always find 16 neighbors |
|---|
| 139 | ; imoms interpolation along longitude and latitude |
|---|
| 140 | ; |
|---|
| 141 | short = where(indexlat LT jpja-2L AND indexlat GE 1L) |
|---|
| 142 | ilon = indexlon[short] |
|---|
| 143 | ilat = indexlat[short] |
|---|
| 144 | ; |
|---|
| 145 | delta = alat[ilat+1L]-alat[ilat] |
|---|
| 146 | IF max(abs(delta-delta[0])) GE 1.e-6 THEN stop |
|---|
| 147 | delta = delta[0] |
|---|
| 148 | ; |
|---|
| 149 | d0 = (alat[ilat-1]-olat[short])/delta |
|---|
| 150 | IF min(d0, max = ma) LE -2 THEN stop |
|---|
| 151 | IF ma GT -1 THEN stop |
|---|
| 152 | wy0 = imoms3(temporary(d0)) |
|---|
| 153 | d1 = (alat[ilat ]-olat[short])/delta |
|---|
| 154 | IF min(d1, max = ma) LE -1 THEN stop |
|---|
| 155 | IF ma GT 0 THEN stop |
|---|
| 156 | wy1 = imoms3(temporary(d1)) |
|---|
| 157 | d2 = (alat[ilat+1]-olat[short])/delta |
|---|
| 158 | IF min(d2, max = ma) LE 0 THEN stop |
|---|
| 159 | IF ma GT 1 THEN stop |
|---|
| 160 | wy2 = imoms3(temporary(d2)) |
|---|
| 161 | d3 = (alat[ilat+2]-olat[short])/delta |
|---|
| 162 | IF min(d3, max = ma) LE 1 THEN stop |
|---|
| 163 | IF ma GT 2 THEN stop |
|---|
| 164 | wy3 = imoms3(temporary(d3)) |
|---|
| 165 | ; |
|---|
| 166 | mi = min(wy0+wy1+wy2+wy3, max = ma) |
|---|
| 167 | IF abs(mi-1) GE 1.e-6 THEN stop |
|---|
| 168 | IF abs(ma-1) GE 1.e-6 THEN stop |
|---|
| 169 | ; |
|---|
| 170 | delta = alon[ilon]-alon[ilon-1] |
|---|
| 171 | IF max(abs(delta-delta[0])) GE 1.e-6 THEN stop |
|---|
| 172 | delta = delta[0] |
|---|
| 173 | ; |
|---|
| 174 | d0 = (alon[ilon-1]-olon[short])/delta |
|---|
| 175 | IF min(d0, max = ma) LE -2 THEN stop |
|---|
| 176 | IF ma GT -1 THEN stop |
|---|
| 177 | wx0 = imoms3(temporary(d0)) |
|---|
| 178 | d1 = (alon[ilon ]-olon[short])/delta |
|---|
| 179 | IF min(d1, max = ma) LE -1 THEN stop |
|---|
| 180 | IF ma GT 0 THEN stop |
|---|
| 181 | wx1 = imoms3(temporary(d1)) |
|---|
| 182 | d2 = (alon[ilon+1]-olon[short])/delta |
|---|
| 183 | IF min(d2, max = ma) LE 0 THEN stop |
|---|
| 184 | IF ma GT 1 THEN stop |
|---|
| 185 | wx2 = imoms3(temporary(d2)) |
|---|
| 186 | d3 = (alon[ilon+2]-olon[short])/delta |
|---|
| 187 | IF min(d3, max = ma) LE 1 THEN stop |
|---|
| 188 | IF ma GT 2 THEN stop |
|---|
| 189 | wx3 = imoms3(temporary(d3)) |
|---|
| 190 | ; |
|---|
| 191 | mi = min(wx0+wx1+wx2+wx3, max = ma) |
|---|
| 192 | IF abs(mi-1) GE 1.e-6 THEN stop |
|---|
| 193 | IF abs(ma-1) GE 1.e-6 THEN stop |
|---|
| 194 | ; |
|---|
| 195 | ; line 0 |
|---|
| 196 | xaddr[0, short] = ilon - 1L |
|---|
| 197 | xaddr[1, short] = ilon |
|---|
| 198 | xaddr[2, short] = ilon + 1L |
|---|
| 199 | xaddr[3, short] = ilon + 2L |
|---|
| 200 | yaddr[0, short] = ilat - 1L |
|---|
| 201 | yaddr[1, short] = yaddr[0, short] |
|---|
| 202 | yaddr[2, short] = yaddr[0, short] |
|---|
| 203 | yaddr[3, short] = yaddr[0, short] |
|---|
| 204 | weig[0, short] = wx0 * wy0 |
|---|
| 205 | weig[1, short] = wx1 * wy0 |
|---|
| 206 | weig[2, short] = wx2 * wy0 |
|---|
| 207 | weig[3, short] = wx3 * wy0 |
|---|
| 208 | ; line 1 |
|---|
| 209 | xaddr[4, short] = ilon - 1L |
|---|
| 210 | xaddr[5, short] = ilon |
|---|
| 211 | xaddr[6, short] = ilon + 1L |
|---|
| 212 | xaddr[7, short] = ilon + 2L |
|---|
| 213 | yaddr[4, short] = ilat |
|---|
| 214 | yaddr[5, short] = ilat |
|---|
| 215 | yaddr[6, short] = ilat |
|---|
| 216 | yaddr[7, short] = ilat |
|---|
| 217 | weig[4, short] = wx0 * wy1 |
|---|
| 218 | weig[5, short] = wx1 * wy1 |
|---|
| 219 | weig[6, short] = wx2 * wy1 |
|---|
| 220 | weig[7, short] = wx3 * wy1 |
|---|
| 221 | ; line 2 |
|---|
| 222 | xaddr[8, short] = ilon - 1L |
|---|
| 223 | xaddr[9, short] = ilon |
|---|
| 224 | xaddr[10, short] = ilon + 1L |
|---|
| 225 | xaddr[11, short] = ilon + 2L |
|---|
| 226 | yaddr[8 , short] = ilat + 1L |
|---|
| 227 | yaddr[9 , short] = yaddr[8 , short] |
|---|
| 228 | yaddr[10, short] = yaddr[8 , short] |
|---|
| 229 | yaddr[11, short] = yaddr[8 , short] |
|---|
| 230 | weig[8 , short] = wx0 * wy2 |
|---|
| 231 | weig[9 , short] = wx1 * wy2 |
|---|
| 232 | weig[10, short] = wx2 * wy2 |
|---|
| 233 | weig[11, short] = wx3 * wy2 |
|---|
| 234 | ; line 3 |
|---|
| 235 | xaddr[12, short] = ilon - 1L |
|---|
| 236 | xaddr[13, short] = ilon |
|---|
| 237 | xaddr[14, short] = ilon + 1L |
|---|
| 238 | xaddr[15, short] = ilon + 2L |
|---|
| 239 | yaddr[12, short] = ilat + 2L |
|---|
| 240 | yaddr[13, short] = yaddr[12, short] |
|---|
| 241 | yaddr[14, short] = yaddr[12, short] |
|---|
| 242 | yaddr[15, short] = yaddr[12, short] |
|---|
| 243 | weig[12, short] = wx0 * wy3 |
|---|
| 244 | weig[13, short] = wx1 * wy3 |
|---|
| 245 | weig[14, short] = wx2 * wy3 |
|---|
| 246 | weig[15, short] = wx3 * wy3 |
|---|
| 247 | ; |
|---|
| 248 | mi = min(total(weig[*, short], 1), max = ma) |
|---|
| 249 | IF abs(mi-1) GE 1.e-6 THEN stop |
|---|
| 250 | IF abs(ma-1) GE 1.e-6 THEN stop |
|---|
| 251 | ; |
|---|
| 252 | ; for the ocean points located between the atm lines |
|---|
| 253 | ; jpja-2 and jpja-1 or between the atm lines 0 and 1 |
|---|
| 254 | ; linear interpolation between line 1 and line 2 |
|---|
| 255 | ; |
|---|
| 256 | short = where(indexlat EQ jpja-2 OR indexlat EQ 0) |
|---|
| 257 | IF short[0] NE -1 THEN BEGIN |
|---|
| 258 | ilon = indexlon[short] |
|---|
| 259 | ilat = indexlat[short] |
|---|
| 260 | ; |
|---|
| 261 | delta = alat[ilat+1]-alat[ilat] |
|---|
| 262 | IF max(abs(delta-delta[0])) GE 1.e-6 THEN stop |
|---|
| 263 | delta = delta[0] |
|---|
| 264 | ; |
|---|
| 265 | d1 = (alat[ilat ]-olat[short])/delta |
|---|
| 266 | IF min(d1, max = ma) LE -1 THEN stop |
|---|
| 267 | IF ma GT 0 THEN stop |
|---|
| 268 | wy1 = 1.+ temporary(d1) |
|---|
| 269 | d2 = (alat[ilat+1]-olat[short])/delta |
|---|
| 270 | IF min(d2, max = ma) LE 0 THEN stop |
|---|
| 271 | IF ma GT 1 THEN stop |
|---|
| 272 | wy2 = 1.- temporary(d2) |
|---|
| 273 | ; |
|---|
| 274 | mi = min(wy1+wy2, max = ma) |
|---|
| 275 | IF abs(mi-1) GE 1.e-6 THEN stop |
|---|
| 276 | IF abs(ma-1) GE 1.e-6 THEN stop |
|---|
| 277 | ; but imoms3 along the longitude |
|---|
| 278 | delta = alon[ilon]-alon[ilon-1] |
|---|
| 279 | IF max(abs(delta-delta[0])) GE 1.e-6 THEN stop |
|---|
| 280 | delta = delta[0] |
|---|
| 281 | ; |
|---|
| 282 | d0 = (alon[ilon-1]-olon[short])/delta |
|---|
| 283 | IF min(d0, max = ma) LE -2 THEN stop |
|---|
| 284 | IF ma GT -1 THEN stop |
|---|
| 285 | wx0 = imoms3(temporary(d0)) |
|---|
| 286 | d1 = (alon[ilon ]-olon[short])/delta |
|---|
| 287 | IF min(d1, max = ma) LE -1 THEN stop |
|---|
| 288 | IF ma GT 0 THEN stop |
|---|
| 289 | wx1 = imoms3(temporary(d1)) |
|---|
| 290 | d2 = (alon[ilon+1]-olon[short])/delta |
|---|
| 291 | IF min(d2, max = ma) LE 0 THEN stop |
|---|
| 292 | IF ma GT 1 THEN stop |
|---|
| 293 | wx2 = imoms3(temporary(d2)) |
|---|
| 294 | d3 = (alon[ilon+2]-olon[short])/delta |
|---|
| 295 | IF min(d3, max = ma) LE 1 THEN stop |
|---|
| 296 | IF ma GT 2 THEN stop |
|---|
| 297 | wx3 = imoms3(temporary(d3)) |
|---|
| 298 | ; |
|---|
| 299 | mi = min(wx0+wx1+wx2+wx3, max = ma) |
|---|
| 300 | IF abs(mi-1) GE 1.e-6 THEN stop |
|---|
| 301 | IF abs(ma-1) GE 1.e-6 THEN stop |
|---|
| 302 | ; line 1 |
|---|
| 303 | xaddr[0, short] = ilon - 1L |
|---|
| 304 | xaddr[1, short] = ilon |
|---|
| 305 | xaddr[2, short] = ilon + 1L |
|---|
| 306 | xaddr[3, short] = ilon + 2L |
|---|
| 307 | yaddr[0, short] = ilat |
|---|
| 308 | yaddr[1, short] = ilat |
|---|
| 309 | yaddr[2, short] = ilat |
|---|
| 310 | yaddr[3, short] = ilat |
|---|
| 311 | weig[0, short] = wx0 * wy1 |
|---|
| 312 | weig[1, short] = wx1 * wy1 |
|---|
| 313 | weig[2, short] = wx2 * wy1 |
|---|
| 314 | weig[3, short] = wx3 * wy1 |
|---|
| 315 | ; line 2 |
|---|
| 316 | xaddr[4, short] = ilon - 1L |
|---|
| 317 | xaddr[5, short] = ilon |
|---|
| 318 | xaddr[6, short] = ilon + 1L |
|---|
| 319 | xaddr[7, short] = ilon + 2L |
|---|
| 320 | yaddr[4, short] = ilat + 1L |
|---|
| 321 | yaddr[5, short] = yaddr[4, short] |
|---|
| 322 | yaddr[6, short] = yaddr[4, short] |
|---|
| 323 | yaddr[7, short] = yaddr[4, short] |
|---|
| 324 | weig[4, short] = wx0 * wy2 |
|---|
| 325 | weig[5, short] = wx1 * wy2 |
|---|
| 326 | weig[6, short] = wx2 * wy2 |
|---|
| 327 | weig[7, short] = wx3 * wy2 |
|---|
| 328 | ; |
|---|
| 329 | mi = min(total(weig[*, short], 1), max = ma) |
|---|
| 330 | IF abs(mi-1) GE 1.e-6 THEN stop |
|---|
| 331 | IF abs(ma-1) GE 1.e-6 THEN stop |
|---|
| 332 | ; |
|---|
| 333 | ENDIF |
|---|
| 334 | ; |
|---|
| 335 | ; for the ocean points located below the line 0 |
|---|
| 336 | ; Interpolation only along the longitude |
|---|
| 337 | ; |
|---|
| 338 | short = where(indexlat EQ -1) |
|---|
| 339 | IF short[0] NE -1 THEN BEGIN |
|---|
| 340 | ilon = indexlon[short] |
|---|
| 341 | ; |
|---|
| 342 | delta = alon[ilon]-alon[ilon-1] |
|---|
| 343 | IF max(abs(delta-delta[0])) GE 1.e-6 THEN stop |
|---|
| 344 | delta = delta[0] |
|---|
| 345 | ; |
|---|
| 346 | d0 = (alon[ilon-1]-olon[short])/delta |
|---|
| 347 | IF min(d0, max = ma) LE -2 THEN stop |
|---|
| 348 | IF ma GT -1 THEN stop |
|---|
| 349 | wx0 = imoms3(temporary(d0)) |
|---|
| 350 | d1 = (alon[ilon ]-olon[short])/delta |
|---|
| 351 | IF min(d1, max = ma) LE -1 THEN stop |
|---|
| 352 | IF ma GT 0 THEN stop |
|---|
| 353 | wx1 = imoms3(temporary(d1)) |
|---|
| 354 | d2 = (alon[ilon+1]-olon[short])/delta |
|---|
| 355 | IF min(d2, max = ma) LE 0 THEN stop |
|---|
| 356 | IF ma GT 1 THEN stop |
|---|
| 357 | wx2 = imoms3(temporary(d2)) |
|---|
| 358 | d3 = (alon[ilon+2]-olon[short])/delta |
|---|
| 359 | IF min(d3, max = ma) LE 1 THEN stop |
|---|
| 360 | IF ma GT 2 THEN stop |
|---|
| 361 | wx3 = imoms3(temporary(d3)) |
|---|
| 362 | ; |
|---|
| 363 | mi = min(wx0+wx1+wx2+wx3, max = ma) |
|---|
| 364 | IF abs(mi-1) GE 1.e-6 THEN stop |
|---|
| 365 | IF abs(ma-1) GE 1.e-6 THEN stop |
|---|
| 366 | ; line 1 |
|---|
| 367 | xaddr[0, short] = ilon - 1L |
|---|
| 368 | xaddr[1, short] = ilon |
|---|
| 369 | xaddr[2, short] = ilon + 1L |
|---|
| 370 | xaddr[3, short] = ilon + 2L |
|---|
| 371 | yaddr[0:3, short] = 0 |
|---|
| 372 | weig[0, short] = wx0 |
|---|
| 373 | weig[1, short] = wx1 |
|---|
| 374 | weig[2, short] = wx2 |
|---|
| 375 | weig[3, short] = wx3 |
|---|
| 376 | ; |
|---|
| 377 | mi = min(total(weig[*, short], 1), max = ma) |
|---|
| 378 | IF abs(mi-1) GE 1.e-6 THEN stop |
|---|
| 379 | IF abs(ma-1) GE 1.e-6 THEN stop |
|---|
| 380 | ; |
|---|
| 381 | ENDIF |
|---|
| 382 | ; |
|---|
| 383 | ; for the ocean points located above jpia-1 |
|---|
| 384 | ; Interpolation only along the longitude |
|---|
| 385 | ; |
|---|
| 386 | short = where(indexlat EQ jpja-1) |
|---|
| 387 | IF short[0] NE -1 THEN BEGIN |
|---|
| 388 | ilon = indexlon[short] |
|---|
| 389 | ; |
|---|
| 390 | delta = alon[ilon]-alon[ilon-1] |
|---|
| 391 | IF max(abs(delta-delta[0])) GE 1.e-6 THEN stop |
|---|
| 392 | delta = delta[0] |
|---|
| 393 | ; |
|---|
| 394 | d0 = (alon[ilon-1]-olon[short])/delta |
|---|
| 395 | IF min(d0, max = ma) LE -2 THEN stop |
|---|
| 396 | IF ma GT -1 THEN stop |
|---|
| 397 | wx0 = imoms3(temporary(d0)) |
|---|
| 398 | d1 = (alon[ilon ]-olon[short])/delta |
|---|
| 399 | IF min(d1, max = ma) LE -1 THEN stop |
|---|
| 400 | IF ma GT 0 THEN stop |
|---|
| 401 | wx1 = imoms3(temporary(d1)) |
|---|
| 402 | d2 = (alon[ilon+1]-olon[short])/delta |
|---|
| 403 | IF min(d2, max = ma) LE 0 THEN stop |
|---|
| 404 | IF ma GT 1 THEN stop |
|---|
| 405 | wx2 = imoms3(temporary(d2)) |
|---|
| 406 | d3 = (alon[ilon+2]-olon[short])/delta |
|---|
| 407 | IF min(d3, max = ma) LE 1 THEN stop |
|---|
| 408 | IF ma GT 2 THEN stop |
|---|
| 409 | wx3 = imoms3(temporary(d3)) |
|---|
| 410 | ; |
|---|
| 411 | mi = min(wx0+wx1+wx2+wx3, max = ma) |
|---|
| 412 | IF abs(mi-1) GE 1.e-6 THEN stop |
|---|
| 413 | IF abs(ma-1) GE 1.e-6 THEN stop |
|---|
| 414 | ; line 1 |
|---|
| 415 | xaddr[0, short] = ilon-1L |
|---|
| 416 | xaddr[1, short] = ilon |
|---|
| 417 | xaddr[2, short] = ilon+1L |
|---|
| 418 | xaddr[3, short] = ilon+2L |
|---|
| 419 | yaddr[0:3, short] = jpja-1L |
|---|
| 420 | weig[0, short] = wx0 |
|---|
| 421 | weig[1, short] = wx1 |
|---|
| 422 | weig[2, short] = wx2 |
|---|
| 423 | weig[3, short] = wx3 |
|---|
| 424 | ; |
|---|
| 425 | mi = min(total(weig[*, short], 1), max = ma) |
|---|
| 426 | IF abs(mi-1) GE 1.e-6 THEN stop |
|---|
| 427 | IF abs(ma-1) GE 1.e-6 THEN stop |
|---|
| 428 | ; |
|---|
| 429 | ENDIF |
|---|
| 430 | ; |
|---|
| 431 | ; Come back to the original index of atm grid without longitudinal overlap. |
|---|
| 432 | ; |
|---|
| 433 | ; |
|---|
| 434 | xaddr = temporary(xaddr) - toadd |
|---|
| 435 | jpia = jpia - 2*toadd |
|---|
| 436 | ; make sure all values are ge 0 |
|---|
| 437 | xaddr = temporary(xaddr) + jpia |
|---|
| 438 | ; range the values between 0 and jpia-1 |
|---|
| 439 | xaddr = temporary(xaddr) mod jpia |
|---|
| 440 | ; |
|---|
| 441 | ; take into account shiftx if needed |
|---|
| 442 | IF shiftx NE 0 THEN xaddr = (temporary(xaddr) - shiftx) MOD jpia |
|---|
| 443 | ; take into account nosouthernline and nonorthernline |
|---|
| 444 | if keyword_set(nosouthernline) then BEGIN |
|---|
| 445 | yaddr = temporary(yaddr) + 1L |
|---|
| 446 | jpja = jpja + 1L |
|---|
| 447 | ENDIF |
|---|
| 448 | if keyword_set(nonorthernline) then jpja = jpja + 1L |
|---|
| 449 | ; take into account revy if needed |
|---|
| 450 | IF revy EQ 1 THEN yaddr = jpja - 1L - temporary(yaddr) |
|---|
| 451 | ; ; |
|---|
| 452 | addr = temporary(yaddr)*jpia+temporary(xaddr) |
|---|
| 453 | ; |
|---|
| 454 | RETURN |
|---|
| 455 | END |
|---|
