;+ ; ; @file_comments ; Build the triangulation for a E-grid type ; ; @categories ; Graphics ; ; @param maskentree {in}{optional}{type=2d array} ; It is a 2d array which will serve to mask the field we will trace after ; with CONTOUR, ...TRIANGULATION=triangule(mask) ; If this argument is not specified, the function use tmask ; ; @keyword BASIC ; Specify that the mask is on a basic grid (use the triangulation for vertical ; cuts and hovmoellers) ; ; @keyword COINMONTE {type=array} ; To obtain the array of "ascending land corner" to be treated with ; completecointerre in the variable array instead of make it pass ; by the global variable twin_corners_up. ; ; @keyword COINDESCEND {type=array} ; See COINMONTE ; ; @keyword SHIFTED ; ; @uses ; common.pro ; ; @history ; Sebastien Masson (smasson\@lodyc.jussieu.fr) ; june 2001 ; ; @version ; $Id$ ; ; @todo ; seb L.152->153 je ne pense pas que ce soit ce que tu voulais dire mais ; c'est la traduction de ce qu'il y avait écrit. Correction si besoin. ; ;- FUNCTION triangule_e, maskentree $ , COINMONTE = coinmonte, COINDESCEND = coindescend $ , SHIFTED = shifted, BASIC = basic ; compile_opt idl2, strictarrsubs ; @cm_4mesh IF NOT keyword_set(key_forgetold) THEN BEGIN @updatenew ENDIF ;--------------------------------------------------------- tempsun = systime(1) ; For key_performance ;------------------------------------------------------------ ; Is the mask given or do we have to take tmask? ;------------------------------------------------------------ ; msk = maskentree sizem = size(msk) nx = sizem[1] ny = sizem[2] ;------------------------------------------------------------ if keyword_set(key_periodic)*(nx EQ jpi) $ AND NOT keyword_set(basic) then BEGIN msk = [msk, msk[0, *]] nx = nx+1 ENDIF ; ; we will find the diamond that must be cut in two triangle using the ; horizontal diagonal. ; index = lindgen(nx, ny) index = index[0:nx-2, 1:ny-2] if n_elements(shifted) EQ 0 then shifted = 1 oddeven = (index/nx+1-shifted) MOD 2 msk1 = msk[index] msk2 = msk[index+1] sum = msk[index-nx+oddeven]+msk[index+nx+oddeven] sum1 = msk2+sum sum2 = msk1+sum ; ; horizontal ; singularpoint = where((msk1 EQ 0 AND sum1 EQ 3) OR (msk1 EQ 1 AND sum1 EQ 0) $ OR (msk2 EQ 0 AND sum2 EQ 3) OR (msk2 EQ 1 AND sum2 EQ 0) $ OR (sum EQ 0 AND (msk1+msk2) EQ 2) ) if singularpoint[0] NE -1 then begin horizontal = index[singularpoint] triang = definetri_e(nx, ny, horizontal, SHIFTED = shifted) ENDIF ELSE triang = definetri_e(nx, ny, SHIFTED = shifted) ; coinmont = index[where(sum EQ 2 AND (msk1+msk2) EQ 0)] ; coindesc = index[where(sum EQ 0 AND (msk1+msk2) EQ 2)] ; ; we keep only the triangles which are outside the land ; but for some reasons we will in fact delete the land diamond ; ; allrecinland = where(sum1+msk1 EQ 0) ; indexallinland = index[allrecinland] ; otherrec = (lindgen(nx, ny))[0:nx-2, 1:ny-2] ; otherrec = different(otherrec, indexallinland) ; ; ; index = lindgen(nx, ny) ; index = index[0:nx-3, 2:ny-3] ; out = inter(index, indexallinland) ; IF out[0] NE -1 THEN begin ; out = inter(out+1, indexallinland) ; IF out[0] NE -1 THEN begin ; out = out-1 ; oddeven = (out/nx+1-shifted) MOD 2 ; out = inter(out-nx+oddeven, otherrec) ; IF out[0] NE -1 THEN begin ; out = inter(out+2*nx, otherrec) ; IF out[0] NE -1 THEN begin ; out = out-(nx+((out/nx+shifted) MOD 2)) ; endif ; endif ; endif ; ENDIF ; help, out ; ; ; index = lindgen(nx, ny) ; index = index[0:nx-3, 2:ny-3] ; out = inter(index, otherrec) ; IF out[0] NE -1 THEN begin ; out = inter(out+1, otherrec) ; IF out[0] NE -1 THEN begin ; out = out-1 ; oddeven = (out/nx+1-shifted) MOD 2 ; out = inter(out-nx+oddeven, indexallinland) ; IF out[0] NE -1 THEN begin ; out = inter(out+2*nx, indexallinland) ; IF out[0] NE -1 THEN begin ; out = out-(nx+((out/nx+shifted) MOD 2)) ; endif ; endif ; endif ; endif ; help, out ; ; ; IF out[0] EQ -1 THEN out = different(indexallinland, out) ELSE out = indexallinland ; triout = numtri(out, nx, ny) ; triout = [triout, triout+1] ; goodtri = lindgen(2*(nx-1)*(ny-1)) ; goodtri = different(goodtri, triout) ; triang = triang[*, temporary(goodtri)] ; ; ;------------------------------------------------------------ ; When key_periodic equal 1, triang is a list of indexes's array which ; have a surplus column. ; We have to put it back to the initial matrix by putting indexes of ; the last column equal to these of the last column... ;------------------------------------------------------------ tempdeux = systime(1) ; For key_performance =2 if keyword_set(key_periodic)*(nx-1 EQ jpi) $ AND NOT keyword_set(basic) then BEGIN indicey = triang/nx indicex = triang-indicey*nx nx = nx-1 liste = where(indicex EQ nx) if liste[0] NE -1 then indicex[liste] = 0 triang = indicex+nx*indicey nx = nx+1 ; if coinmont[0] NE -1 then begin ; indicey = coinmont/nx ; indicex = coinmont-indicey*nx ; nx = nx-1 ; liste = where(indicex EQ nx) ; if liste[0] NE -1 THEN indicex[liste] = 0 ; coinmont = indicex+nx*indicey ; nx = nx+1 ; endif ; if coindesc[0] NE -1 then begin ; indicey = coindesc/nx ; indicex = coindesc-indicey*nx ; nx = nx-1 ; liste = where(indicex EQ nx) ; if liste[0] NE -1 THEN indicex[liste] = 0 ; coindesc = indicex+nx*indicey ; nx = nx+1 ; endif endif IF testvar(var = key_performance) EQ 2 THEN $ print, 'temps triangule: finitions', systime(1)-tempdeux ;------------------------------------------------------------ if arg_present(coinmonte) THEN coinmonte = coinmont ELSE twin_corners_up = coinmont if arg_present(coindescend) THEN coindescend = coindesc ELSE twin_corners_dn = coindesc IF NOT keyword_set(key_forgetold) THEN BEGIN @updateold ENDIF ;------------------------------------------------------------ IF keyword_set(key_performance) THEN print, 'temps triangule', systime(1)-tempsun return, triang END