;+ ; ; @file_comments ; Choose the grid which must be used to do the graph in function of ; vargrid and send back corresponding parameters calculated in ; domdef and reduced at the domain defined by ; domdef (contrarily to ; grandegrille) ; BEWARE!! The choice of the grid is made from the value of the ; global variable vargrid, which can be equal to 'T', 'U', 'V', 'W' ou 'F'. ; ; @categories ; Grid ; ; @keyword TRI ; This keyword serve to obtain the triangulation which ; refer to the grid but only on the part of the zoom. This array of triangulation ; is passed in the variable we have equate at TRI. ; For example: grille,...,tri=triangulation_reduite. ; This keyword is used in plt. ; ; @keyword WDEPTH ; To specify that the field is at W depth instead of T ; depth (automatically activated if vargrid eq 'W') ; ; @keyword FORPLT ; In plt, we want that land points, ; glam and gphi, be equal to glamt and ; gphit regardless of the grid. ; ; @keyword NOTRI ; Useful only when TRI is activated. In this case, grille send back -1 in the ; variable tri even if the variable of the common triangles_list is defined ; and different of -1 ; ; @keyword _EXTRA ; Used to pass keywords ; ; @keyword TOUT ; ; @param MASK {out}{optional} ; For the definition, see domdef and the management of subdomains on the web. ; ; @param GLAM {out}{optional} ; For the definition, see domdef and the management of subdomains on the web. ; ; @param GPHI {out}{optional} ; For the definition, see domdef and the management of subdomains on the web. ; ; @param GDEP {out}{optional} ; For the definition, see domdef and the management of subdomains on the web. ; ; @param NX {out}{optional} ; For the definition, see domdef and the management of subdomains on the web. ; ; @param NY {out}{optional} ; For the definition, see domdef and the management of subdomains on the web. ; ; @param NZ {out}{optional} ; For the definition, see domdef and the management of subdomains on the web. ; ; @param FIRSTX {out}{optional} ; For the definition, see domdef and the management of subdomains on the web. ; ; @param FIRSTY {out}{optional} ; For the definition, see domdef and the management of subdomains on the web. ; ; @param FIRSTZ {out}{optional} ; For the definition, see domdef and the management of subdomains on the web. ; ; @param LASTX {out}{optional} ; For the definition, see domdef and the management of subdomains on the web. ; ; @param LASTY {out}{optional} ; For the definition, see domdef and the management of subdomains on the web. ; ; @param LASTZ {out}{optional} ; For the definition, see domdef and the management of subdomains on the web. ; ; @param E1 {out}{optional} ; For the definition, see domdef and the management of subdomains on the web. ; ; @param E2 {out}{optional} ; For the definition, see domdef and the management of subdomains on the web. ; ; @param E3 {out}{optional} ; For the definition, see domdef and the management of subdomains on the web. ; ; @uses ; cm_4mesh ; cm_4data ; ; @restrictions ; Use the variable vargrid ; ; @restrictions ; Vargrid must be 'T', 'W', 'U', 'V' or 'F' ; ; @history ; Sebastien Masson (smasson\@lodyc.jussieu.fr) ; 12/2/1999 ; 10/11/1999 /forplt ; ; @version ; $Id$ ; ; @todo ; Comment ecrire la remarque sur les inputs? ; ;- PRO grille, mask, glam, gphi, gdep, nx, ny, nz $ , firstx, firsty, firstz, lastx, lasty, lastz $ , e1, e2, e3 $ , TRI=tri, NOTRI=notri, TOUT=tout, FORPLT=forplt $ , TYPE=type, WDEPTH=wdepth, _EXTRA=ex ; compile_opt idl2, strictarrsubs ; @cm_4mesh @cm_4data @cm_4cal ; for jpt IF NOT keyword_set(key_forgetold) THEN BEGIN @updatenew ENDIF ;--------------------- tempsun = systime(1) ; For key_performance ;------------------------------------------------------------ vargrid = strupcase(strmid(vargrid,0,/reverse_offset)) ; if vargrid eq 'W' then wdepth = 1 if keyword_set(tout) then begin savedbox = 1b saveboxparam, 'boxparam4grille.dat' domdef, gridtype = vargrid, _EXTRA = ex endif tempdeux = systime(1) ; For key_performance =2 ;------------------------------------------------------------ ;------------------------------------------------------------ IF keyword_set(wdepth) THEN BEGIN firstz = firstzw lastz = lastzw nz = nzw ENDIF ELSE BEGIN firstz = firstzt lastz = lastzt nz = nzt ENDELSE ;------------------------------------------------------------ ;------------------------------------------------------------ CASE 1 OF ;------------------------------------------------------------ ; grid T and W ;------------------------------------------------------------ vargrid eq 'T' OR vargrid eq 'W' : begin ;scalars nx = nxt ny = nyt firstx = firstxt firsty = firstyt lastx = lastxt lasty = lastyt ;2d vectors IF arg_present(glam) THEN glam = glamt[firstx:lastx, firsty:lasty] IF arg_present(gphi) THEN gphi = gphit[firstx:lastx, firsty:lasty] IF arg_present(e1) THEN e1 = e1t[firstx:lastx, firsty:lasty] IF arg_present(e2) THEN e2 = e2t[firstx:lastx, firsty:lasty] ;3d vectors IF keyword_set(forplt) THEN mask = tmask[firstx:lastx, firsty:lasty, firstz] $ ELSE IF arg_present(mask) THEN mask = tmask[firstx:lastx, firsty:lasty, firstz:lastz] end ;------------------------------------------------------------ ; grid U ;------------------------------------------------------------ vargrid eq 'U': begin ;scalars nx = nxu ny = nyu firstx = firstxu firsty = firstyu lastx = lastxu lasty = lastyu ;2d vectors IF arg_present(glam) THEN glam = glamu[firstx:lastx, firsty:lasty] IF arg_present(gphi) THEN gphi = gphiu[firstx:lastx, firsty:lasty] if keyword_set(forplt) then BEGIN mask = 1b-tmask[firstx:lastx, firsty:lasty, firstz] eastboarder = mask-shift(mask, 1, 0)*mask westboarder = mask-shift(mask, -1, 0)*mask if key_periodic NE 1 OR nx NE jpi then westboarder[nx-1, *] = 0b tmp1 = shift(eastboarder, 0, 1) tmp1[*, 0] = 0b tmp2 = shift(eastboarder, 0, -1) tmp2[*, ny-1] = 0b add = (temporary(tmp1)+temporary(tmp2))*(1b-eastboarder)*(1b-temporary(westboarder)) eastboarder = temporary(eastboarder)+temporary(add) tmp1 = (mask+shift(mask, 0, -1)+shift(mask, 0, 1)) NE 1b tmp1[*, ny-1] = 1b tmp1[*, 0] = 1b tmp2 = (mask+shift(mask, -1, 0)+shift(mask, 1, 0)) NE 1b if key_periodic NE 1 OR nx NE jpi then begin tmp2[nx-1, *] = 1b tmp2[0, *] = 0b endif no1 = temporary(tmp1)*temporary(tmp2) tmp = temporary(eastboarder)*temporary(no1)*mask mask[0:nx-2, *] = 0b tmp = temporary(tmp)+temporary(mask) tmp = where(tmp GE 1) if tmp[0] NE -1 then begin glam[tmp] = (glamt[firstx:lastx, firsty:lasty])[tmp] gphi[tmp] = (gphit[firstx:lastx, firsty:lasty])[tmp] endif ENDIF IF arg_present(e1) THEN e1 = e1u[firstx:lastx, firsty:lasty] IF arg_present(e2) THEN e2 = e2u[firstx:lastx, firsty:lasty] ;3d vectors IF keyword_set(forplt) THEN mask = (umask())[firstx:lastx, firsty:lasty, firstz] $ ELSE IF arg_present(mask) THEN mask = (umask())[firstx:lastx, firsty:lasty, firstz:lastz] end ;------------------------------------------------------------ ; grid V ;------------------------------------------------------------ vargrid eq 'OPAPTDHV' or vargrid eq 'OPAPT3DV' $ or vargrid eq 'V': begin ;scalars nx = nxv ny = nyv firstx = firstxv firsty = firstyv lastx = lastxv lasty = lastyv ;2d vectors IF arg_present(glam) THEN glam = glamv[firstx:lastx, firsty:lasty] IF arg_present(gphi) THEN gphi = gphiv[firstx:lastx, firsty:lasty] if keyword_set(forplt) then BEGIN mask = 1b-tmask[firstx:lastx, firsty:lasty, firstz] northboarder = mask-shift(mask, 0, 1)*mask southboarder = mask-shift(mask, 0, -1)*mask southboarder[*, ny-1] = 0b tmp1 = shift(northboarder, -1, 0) if key_periodic NE 1 OR nx NE jpi then tmp1[nx-1, *] = 0b tmp2 = shift(northboarder, 1, 0) if key_periodic NE 1 OR nx NE jpi then tmp2[0, *] = 0b add = (temporary(tmp1)+temporary(tmp2))*(1b-northboarder)*(1b-southboarder) northboarder = temporary(northboarder)+temporary(add) tmp1 = (mask+shift(mask, 0, -1)+shift(mask, 0, 1)) NE 1b tmp1[*, ny-1] = 1b tmp1[*, 0] = 0b tmp2 = (mask+shift(mask, -1, 0)+shift(mask, 1, 0)) NE 1b if key_periodic NE 1 OR nx NE jpi then begin tmp2[nx-1, *] = 1b tmp2[0, *] = 1b endif no1 = temporary(tmp1)*temporary(tmp2) tmp = temporary(northboarder)*mask*temporary(no1) mask[*, 0:ny-2] = 0b tmp = temporary(tmp)+temporary(mask) tmp = where(tmp GE 1) if tmp[0] NE -1 then begin glam[tmp] = (glamt[firstx:lastx, firsty:lasty])[tmp] gphi[tmp] = (gphit[firstx:lastx, firsty:lasty])[tmp] endif ENDIF IF arg_present(e1) THEN e1 = e1v[firstx:lastx, firsty:lasty] IF arg_present(e2) THEN e2 = e2v[firstx:lastx, firsty:lasty] ;3d vecteurs IF keyword_set(forplt) THEN mask = (vmask())[firstx:lastx, firsty:lasty, firstz] $ ELSE IF arg_present(mask) THEN mask = (vmask())[firstx:lastx, firsty:lasty, firstz:lastz] end ;------------------------------------------------------------ ; grid F ;------------------------------------------------------------ vargrid eq 'OPAPTDHF' or vargrid eq 'OPAPT3DF' $ or vargrid eq 'F': begin ;scalars nx = nxf ny = nyf firstx = firstxf firsty = firstyf lastx = lastxf lasty = lastyf ;2d vectors IF arg_present(glam) THEN glam = glamf[firstx:lastx, firsty:lasty] IF arg_present(gphi) THEN gphi = gphif[firstx:lastx, firsty:lasty] if keyword_set(forplt) then BEGIN mask = 1b-tmask[firstx:lastx, firsty:lasty, firstz] eastboarder = mask-shift(mask, 1, 0)*mask westboarder = mask-shift(mask, -1, 0)*mask westboarder[nx-1, *] = 0b northboarder = mask-shift(mask, 0, 1)*mask southboarder = mask-shift(mask, 0, -1)*mask southboarder[*, ny-1] = 0b tmp1 = shift(northboarder, -1, 0) if key_periodic NE 1 OR nx NE jpi then tmp1[nx-1, *] = 0b tmp2 = shift(northboarder, 1, 0) if key_periodic NE 1 OR nx NE jpi then tmp2[0, *] = 0b add = (temporary(tmp1)+temporary(tmp2))*(1b-northboarder)*(1b-southboarder) northboarder = temporary(northboarder)+temporary(add) tmp1 = shift(eastboarder, 0, 1) tmp1[*, 0] = 0b tmp2 = shift(eastboarder, 0, -1) tmp2[*, ny-1] = 0b add = (temporary(tmp1)+temporary(tmp2))*(1b-eastboarder)*(1b-temporary(westboarder)) eastboarder = temporary(eastboarder)+temporary(add) tmp1 = (mask+shift(mask, 0, -1)+shift(mask, 0, 1)) NE 1b tmp1[*, ny-1] = 1b tmp1[*, 0] = 1b tmp2 = (mask+shift(mask, -1, 0)+shift(mask, 1, 0)) NE 1b if key_periodic NE 1 OR nx NE jpi then begin tmp2[nx-1, *] = 1b tmp2[0, *] = 1b endif no1 = temporary(tmp1)*temporary(tmp2) tmp = (temporary(northboarder)+temporary(eastboarder))*mask*temporary(no1) mask[0:nx-2, *] = 0b mask[*, 0:ny-2] = 0b tmp = temporary(tmp)+temporary(mask) tmp = where(tmp GE 1) if tmp[0] NE -1 then begin glam[tmp] = (glamt[firstx:lastx, firsty:lasty])[tmp] gphi[tmp] = (gphit[firstx:lastx, firsty:lasty])[tmp] endif ENDIF IF arg_present(e1) THEN e1 = e1f[firstx:lastx, firsty:lasty] IF arg_present(e2) THEN e2 = e2f[firstx:lastx, firsty:lasty] ;3d vectors IF keyword_set(forplt) THEN mask = (fmask())[firstx:lastx, firsty:lasty, firstz] $ ELSE IF arg_present(mask) THEN mask = (fmask())[firstx:lastx, firsty:lasty, firstz:lastz] END ;------------------------------------------------------------ ELSE:BEGIN ras = report('Wrong definition of vargrid = '+vargrid+'. Only T, U, V, W or F are acceptable') stop END ENDCASE IF testvar(var = key_performance) EQ 2 THEN $ print, 'temps grille: attribution des scalaires, vecteurs et tableaux ', systime(1)-tempdeux ; ;------------------------------------------------------------ ;------------------------------------------------------------ ;------------------------------------------------------------ ; Variables refering to the vertical dimension ;------------------------------------------------------------ ;------------------------------------------------------------ ;------------------------------------------------------------ ; ; tempdeux = systime(1) ; For key_performance =2 if keyword_set(wdepth) then begin gdep = gdepw[firstz:lastz] e3 = e3w[firstz:lastz] endif else begin gdep = gdept[firstz:lastz] e3 = e3t[firstz:lastz] ENDELSE ; for the vertical sections with partial steps IF keyword_set(type) AND keyword_set(key_partialstep) THEN BEGIN CASE 1 OF type EQ 'xz' AND ny EQ 1:BEGIN bottom = total(tmask[firstx:lastx, firsty:lasty, firstz:lastz], 3) good = where(bottom NE 0 AND bottom NE nz+keyword_set(wdepth)) bottom = lindgen(nx)+(bottom-1l+keyword_set(wdepth))*nx IF good[0] NE -1 THEN BEGIN bottom = bottom[good] IF lastz EQ jpk-1 THEN gdep[nz-1] = max(hdepw) gdep = replicate(1, nx)#gdep if keyword_set(wdepth) THEN $ truegdep = hdepw[firstx:lastx, firsty:lasty] $ ELSE truegdep = hdept[firstx:lastx, firsty:lasty] gdep[bottom] = truegdep[good] ENDIF END type EQ 'yz' AND nx EQ 1:BEGIN bottom = total(tmask[firstx:lastx, firsty:lasty, firstz:lastz], 3) good = where(bottom NE 0 AND bottom NE nz+keyword_set(wdepth)) bottom = lindgen(ny)+(bottom-1l+keyword_set(wdepth))*ny IF good[0] NE -1 THEN BEGIN bottom = bottom[good] IF lastz EQ jpk-1 THEN gdep[nz-1] = max(hdepw) gdep = replicate(1, ny)#gdep if keyword_set(wdepth) THEN $ truegdep = hdepw[firstx:lastx, firsty:lasty] $ ELSE truegdep = hdept[firstx:lastx, firsty:lasty] gdep[bottom] = truegdep[good] ENDIF END ELSE: ENDCASE ENDIF ; for the vertical sections with roms IF keyword_set(type) AND n_elements(romszinfos) NE 0 THEN BEGIN romsdp = romsdepth() IF romsdp[0] NE -1 THEN BEGIN IF jpt EQ 1 THEN BEGIN CASE type OF 'z':gdep = moyenne(temporary(romsdp), 'xy') 'xz':gdep = moyenne(temporary(romsdp), 'y') 'yz':gdep = moyenne(temporary(romsdp), 'x') 'zt':gdep = moyenne(temporary(romsdp), 'xy') ELSE: ENDCASE ENDIF ELSE BEGIN CASE type OF 'z':gdep = moyenne(temporary(romsdp), 'xyt') 'xz':gdep = grossemoyenne(temporary(romsdp), 'yt') 'yz':gdep = grossemoyenne(temporary(romsdp), 'xt') 'zt':gdep = grossemoyenne(temporary(romsdp), 'xy') ELSE: ENDCASE ENDELSE IF size(gdep, /n_dimensions) EQ 2 THEN $ gdep = remplit(gdep, niter = 32000, mask = gdep LT 1.E19, /basique, /fillxdir) ENDIF ENDIF ; IF testvar(var = key_performance) EQ 2 THEN $ print, 'temps grille: Variables se rapportant a la dimension verticale ', systime(1)-tempdeux ;------------------------------------------------------------ ; Triangulation vector when TRI is activated. ;------------------------------------------------------------ if arg_present(TRI) then $ if triangles_list[0] EQ -1 OR keyword_set(notri) then tri = -1 ELSE BEGIN tempdeux = systime(1) ; pour key_performance =2 msk = bytarr(jpi, jpj) msk[firstx:lastx, firsty:lasty] = 1 ind = where( msk[triangles_list[0, *]]*msk[triangles_list[1, *]]*msk[triangles_list[2, *]] EQ 1 ) tri = triangles_list[*, ind]-(firstx+firsty*jpi) y = tri/jpi x = tri-y*jpi tri = x+y*nx IF testvar(var = key_performance) EQ 2 THEN $ print, 'temps grille: decoupage de la triangulation ', systime(1)-tempdeux ENDELSE ;------------------------------------------------------------------ ; To make sure there is not any degenerated dimension (=1) ;------------------------------------------------------------------ ; mask=reform(mask, /over) ; glam=reform(glam, /over) ; gphi=reform(gphi, /over) ; gdep=reform(gdep, /over) ; e1=reform(e1, /over) ; e2=reform(e2, /over) ; e3=reform(e3, /over) if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4grille.dat' if keyword_set(key_performance) THEN print, 'temps grille', systime(1)-tempsun ;------------------------------------------------------------ IF NOT keyword_set(key_forgetold) THEN BEGIN @updateold ENDIF ;--------------------- return end