PRO plt_map, cmd, iplot, win, iover, landscape @common @com_eg ;; ;;------------------------------------- ;; ;; make window plot ;; ;; EG 19-2-99 ;;------------------------------------- ; ; ; cmd = command line of window ; win = position of window (petit[,,,]) ; iover = overlay index ; ; empty window case for on=2 ; ioverchk = iover cmd_wrk = cmd IF cmd.on EQ 2 THEN return ; ; incompatible options of plt_def and post_it ; trend_typp = cmd.trend ; CASE cmd.timave OF ; '1mm': cmd.trend = '0' ; '3mm': cmd.trend = '0' ; ELSE: ; ENDCASE ; ; ============= ; 1.read data ; ============= ; ; type of plot ; time series of rms deviation on a sigma surface? st_rms = 0 IF strpos(cmd.plt, 'st') NE -1 AND strpos(cmd.var, '@s') GT 4 THEN BEGIN st_rms = 1 sig_surface = strmid(cmd.var, strpos(cmd.var, '@s')+2, strlen(cmd.var)) cmd.var = strmid(cmd.var, 0, strpos(cmd.var, '@s')) legendsurf = ' on sigma='+sig_surface ENDIF CASE cmd.plt OF 'x': cmd.plt = cmd.plt+'_global' 'y': cmd.plt = cmd.plt+'_global' 'z': cmd.plt = cmd.plt+'_global' 's': cmd.plt = cmd.plt+'_global' 'yz': cmd.plt = cmd.plt+'_global' ELSE: ENDCASE ; Definition of plttyp, hotyp, dimplot datyp = def_dptyp(cmd) plttyp = datyp.plttyp hotyp = datyp.hotyp dimplot = datyp.dimplot pltztyp = datyp.pltztyp plt1dtyp = datyp.plt1dtyp splot = datyp.splot hotypchk = hotyp ; stat computations: CASE cmd.plt OF 'xyt': BEGIN & plttyp = 'plt' & hotyp = 'xyt' & dimplot = 2 & END ELSE: ENDCASE ; xy plot on density surface (cmd.var = @s legsurf = '' sig_surf = 0. IF strpos(cmd.var, '@s') GT 4 THEN BEGIN splot = 1 sig_surf = strmid(cmd.var, strpos(cmd.var, '@s')+2, strlen(cmd.var)) cmd.var = strmid(cmd.var, 0, strpos(cmd.var, '@s')) legsurf = ' on sigma='+sig_surf ENDIF ; contour of coasts CASE splot OF 1: c_cont = 0 0: c_cont = (!d.n_colors-1) < 255 ENDCASE ; scatter plot y=f(x) IF strpos(cmd.var, '=f(') GE 1 THEN BEGIN CASE cmd.plt OF 'xy': cmd.plt = cmd.plt+'_global' ELSE: ENDCASE plttyp = 'yfx' dimplot = 1 ENDIF ; vector plot IF strpos(cmd.var, '[') EQ 0 THEN BEGIN vecplot = 1 ; sampling of vector IF vector_sample GE 2 THEN BEGIN vect_sample = ',UNVECTSUR=strtrim(string(vector_sample),2)' ENDIF ELSE vect_sample = '' ENDIF ELSE BEGIN vecplot = 0 ENDELSE ; spectrum spectplt = 0 IF strpos(cmd.plt, '@s') GE 1 THEN BEGIN IF hotyp EQ 't' THEN spectplt = 1 ELSE BEGIN print, ' *** spectrum [@s] applies to 1D time serie only (t_)' stop ENDELSE cmd.plt = strmid(cmd.plt, 0, strpos(cmd.plt, '@s')) ENDIF ; wavelet waveplt = 0 IF strpos(cmd.plt, '@w') GE 1 THEN BEGIN IF hotyp EQ 't' THEN waveplt = 1 ELSE BEGIN print, ' *** wavelet [@w] applies to 1D time serie only (t_)' stop ENDELSE cmd.plt = strmid(cmd.plt, 0, strpos(cmd.plt, '@w')) ENDIF ; running std dev run_stddev = 0 IF strpos(cmd.plt, '@r') GE 1 THEN BEGIN IF hotyp NE 't' THEN BEGIN print, ' *** running [@r] std dev applies to 1D time serie only (t_ & 1m@t412)' stop ENDIF run_stddev = strmid(cmd.plt, strpos(cmd.plt, '@r')+2, strlen(cmd.plt)) cmd.plt = strmid(cmd.plt, 0, strpos(cmd.plt, '@r')) ENDIF ; define grids ; read if new type + triangule ; sw_diffg = 0 def_grid, cmd IF debug_w THEN print, 'def_grid,cmd dans plt_map ' IF debug_w THEN print, 'cmd.grid : ', cmd.grid IF debug_w THEN print, 'vargrid : ', vargrid ; vertical average if 3D field vert_switch = 0 IF plttyp EQ 'plt' AND vert_type ne '0' THEN vert_switch = 1 IF plttyp EQ 'pltt' AND vert_type ne '0' THEN vert_switch = 1 IF plttyp EQ 'plt1d' AND vert_type ne '0' THEN vert_switch = 1 IF plt1dtyp EQ 'z' THEN vert_switch = 0 ; read all data if netcdf output IF cmd.out NE 'cdf' THEN BEGIN all_data_str = '' ENDIF ELSE BEGIN all_data_str = ',/ALL_DATA' END ; =========== ; Read data ; =========== really_1m_st = 1 IF cmd.timave EQ '1y' AND strmid(cmd.plt, 0, 2) EQ 'st' THEN BEGIN @densit_pltmap_read ENDIF ELSE BEGIN pltcmd = 'field = data_read(cmd,'''+hotyp+''','''+plttyp+''','+string(dimplot)+','+string(iover)+all_data_str+', ZMTYP = '''+cmd.plt+''')' printf, nulhis, strcompress(pltcmd) printf, nulhis, ' ' IF execute(pltcmd) EQ 0 THEN stop ENDELSE ; data specific treatments IF cmd.var EQ 'sla' THEN field.origin = 'diff' IF cmd.var EQ 'sosstsst' AND cmd.exp EQ 'HadISST1' THEN BEGIN ; limit data to -1.8 C idx = where(field.data EQ valmask) field.data = field.data > (-1.8) IF idx[0] NE -1 THEN field.data(idx) = valmask ENDIF IF n_elements(field.data) EQ 1 THEN return ; change cmd.var if macro CASE STRMID(cmd.var, 0, 2) OF '@@': BEGIN cmd.var = field.name cmd.var = STRMID(cmd.var, 2, 15) IF debug_w THEN print, 'cmd.var : ', cmd.var END ELSE: ENDCASE IF vecplot EQ 1 THEN BEGIN ; vectors case CASE STRMID(cmd.var2, 0, 2) OF '@@': begin cmd.var2 = STRMID(cmd.var2, 2, 15) END ELSE: ENDCASE ENDIF ; spectum plot IF spectplt EQ 1 THEN plttyp = 'spec' ; wavelet plot IF waveplt EQ 1 THEN plttyp = 'wavelet' xlogax = '' ; density pojection plot (splot=1) IF splot EQ 1 THEN BEGIN ; keep grid attributes izminmesh_b = izminmesh izmaxmesh_b = izmaxmesh jpk_b = jpk jpkglo_b = jpkglo gdept_b = gdept gdepw_b = gdepw e3t_b = e3t e3w_b = e3w tmask_b = tmask ; IF cmd.var EQ 'vosigvol' THEN xlogax = ',xlog = 1' ; modify vertical grid -> sigma n_sig = (sig_max - sig_min)/sig_del + 1 z_sig = sig_min+findgen(n_sig)*sig_del izminmesh = 0 izmaxmesh = n_sig-1 jpk = long(izmaxmesh-izminmesh+1) jpkglo = jpk gdept = z_sig gdepw = gdept e3t = shift(gdept, 1)-gdept e3t[0] = e3t[1] e3w = e3t prof1 = sig_min prof2 = sig_max ; grille,mask,glam,gphi,gdep,nx,ny,nz,premierx,premiery,premierz,dernierx,derniery,dernierz ; added PDW 12/5/04 ; tmask = mask ; added PDW 12/5/04 tmask = reform(reform(tmask[*, *, 0], jpi*jpj)#replicate(1, jpk), jpi, jpj, jpk) domdef END ; ; ====================== ; 2. window attributes ; ====================== ; ; find field attributes ; --------------------- fldatt = {name:'', assos:'', min: 0.0, max: 0.0, int: 0.0, mult: 0.0, add: 0.0, unit: '', mid: 0.0, homin:0.0, homax:0.0, min1d:0.0, max1d:0.0, dmax:0.0} ; extremum fldext = fld_pltext(cmd.var, plttyp, dimplot, hotyp) fldatt.name = fldext.name fldatt.min = fldext.min fldatt.max = fldext.max fldatt.homin = fldext.homin fldatt.homax = fldext.homax fldatt.min1d = fldext.min1d fldatt.max1d = fldext.max1d fldatt.dmax = fldext.dmax fldatt.assos = fldext.assos ; interval + change of unit fldint = fld_pltint(cmd.var, plttyp, dimplot, hotyp) ; help, fldint, /struct fldatt.int = fldint.int fldatt.mult = fldint.mult fldatt.add = fldint.add fldatt.unit = fldint.unit fldatt.mid = fldint.mid ; print, 'plt_map 1 :', fldatt.mult, fldatt.add, fldatt.min, fldatt.max IF fldatt.mult NE 1. OR fldatt.add NE 0. THEN field.units = fldatt.unit fld = field.data*fldatt.mult + fldatt.add IF (where(field.data EQ valmask))[0] NE -1 THEN $ fld[where(field.data EQ valmask)] = valmask ; field 2 if needed IF plttyp EQ 'yfx' OR vecplot EQ 1 THEN BEGIN fldatt2 = fldatt fldext2 = fld_pltext(cmd.var2, plttyp, dimplot, hotyp) fldatt2.name = fldext2.name fldatt2.min = fldext2.min fldatt2.max = fldext2.max fldatt2.homin = fldext2.homin fldatt2.homax = fldext2.homax fldatt2.min1d = fldext2.min1d fldatt2.max1d = fldext2.max1d fldatt2.dmax = fldext2.dmax ; interval + change of unit fldint2 = fld_pltint(cmd.var2, plttyp, dimplot, hotyp) fldatt2.int = fldint2.int fldatt2.mult = fldint2.mult fldatt2.add = fldint2.add fldatt2.unit = fldint2.unit fldatt2.mid = fldint2.mid IF fldatt2.mult NE 1. OR fldatt2.add NE 0. THEN field.units2 = fldatt2.unit fld2 = field.data2*fldatt2.mult + fldatt2.add IF (where(field.data2 EQ valmask))[0] NE -1 THEN $ fld2[where(field.data2 EQ valmask)] = valmask ENDIF ; If running std dev change min/max IF run_stddev GT 0 THEN BEGIN fldatt.homin = 0. fldatt.homax = fldext.dmax ENDIF ; legend text ; ----------- ; IF plttyp EQ "pltt" THEN BEGIN CASE strmid(cmd.trend, 0, 1) OF '1': trd_txt = ' trend' '2': trd_txt = ' drift' '3': trd_txt = ' inverse trend' '4': trd_txt = ' anomaly' '5': trd_txt = ' filter' '7': trd_txt = ' month sel' ELSE: trd_txt = '' ENDCASE ; ENDIF ELSE trd_txt = '' IF field_int EQ 1 AND dimplot EQ 1 THEN BEGIN trd_txt = trd_txt+' integral' ENDIF varlegend = field.legend+trd_txt+' ('+field.units+')' CASE plttyp OF 'yfx': varlegend2 = field.legend2+trd_txt+' ('+field.units2+')' ELSE: BEGIN IF run_stddev GT 0 THEN trd_txt = trd_txt+' running Std Dev ['+string(run_stddev, format = '(I3)')+']' END ENDCASE ; date text CASE strmid(cmd.timave, 0, 2) OF '1m': BEGIN mn = def_month(cmd.timave, cmd.date1) IF strmid(cmd.timave, 0, 3) EQ '1mm' THEN $ date_txt = mn+' ('+strmid(cmd.date1, 3, strlen(cmd.date1)-3)+')' $ ELSE date_txt = mn+' '+strmid(cmd.date1, 0, strlen(cmd.date1)-2) END '3m': BEGIN mn = def_month(cmd.timave, cmd.date1) IF strmid(cmd.timave, 0, 3) EQ '3mm' THEN $ date_txt = mn+' ('+strmid(cmd.date1, 3, strlen(cmd.date1)-3)+')' $ ELSE date_txt = mn+' '+strmid(cmd.date1, 0, strlen(cmd.date1)-2) END ELSE:date_txt = cmd.timave+' '+cmd.date1 ENDCASE CASE plttyp OF 'pltt':BEGIN date_txt = cmd.timave IF strmid(cmd.timave, 1, 2) EQ 'mm' THEN $ date_txt = cmd.timave+' ('+strmid(cmd.date1, 3, strlen(cmd.date1)-3)+')' END 'yfx': BEGIN IF hotyp NE '-' THEN date_txt = cmd.timave+' '+cmd.date1+' '+cmd.spec END 'wavelet': BEGIN date_txt = cmd.timave+' Wavelet' IF strmid(cmd.timave, 1, 2) EQ 'mm' THEN $ date_txt = cmd.timave+' Wavelet ('+strmid(cmd.date1, 3, strlen(cmd.date1)-3)+')' END ELSE: ENDCASE ; min/max/int ; defined if fraction of x/y.range is added to plot domain ; IF free_1d_minmax EQ 'no' THEN fraction = 0. ELSE fraction = 1.0 fraction = 0. CASE dimplot OF 1: BEGIN CASE plttyp OF 'pltt': BEGIN & minc = fldatt.homin & maxc = fldatt.homax & END 'wavelet': BEGIN & minc = fldatt.homin & maxc = fldatt.homax & END 'plt1d': BEGIN & minc = fldatt.min1d & maxc = fldatt.max1d & END 'yfx': BEGIN ; y=f(x) scatter plot IF long(strmid(cmd.trend, 0, 1)) GT 0 THEN BEGIN IF run_stddev EQ 0 THEN BEGIN minc = -fldatt.dmax & maxc = fldatt.dmax minc2 = -fldatt2.dmax & maxc2 = fldatt2.dmax ENDIF ELSE BEGIN minc = 0. & maxc = fldatt.dmax minc2 = 0. & maxc2 = fldatt2.dmax ENDELSE ENDIF ELSE BEGIN minc = fldatt.min1d & maxc = fldatt.max1d minc2 = fldatt2.min1d & maxc2 = fldatt2.max1d ENDELSE END ELSE: BEGIN & minc = 0. & maxc = 0. & END ENDCASE fldatt.int = !VALUES.F_NAN END 2: BEGIN CASE plttyp OF 'pltt': BEGIN & minc = fldatt.homin & maxc = fldatt.homax & END ELSE: BEGIN & minc = fldatt.min & maxc = fldatt.max & END ENDCASE END ENDCASE IF cmd.var NE fld_prev OR cmd.trend NE '0' THEN BEGIN ; field.origin EQ 'diff' OR print, '' CASE dimplot OF 2: BEGIN print, ' '+cmd.var, ' plot attributes ('+plttyp+') : min/max/int=', $ minc, maxc, fldatt.int IF ( fldatt.mult NE 1.0 OR fldatt.add NE 0.0 ) THEN $ print, ' --> Modify field using : ', 'mult/add=' , $ fldatt.mult, fldatt.add, ' to obtain : ', fldatt.unit END 1: BEGIN print, ' '+cmd.var, ' plot attributes ('+plttyp+') : min/max=',minc, maxc IF ( fldatt.mult NE 1.0 OR fldatt.add NE 0.0 ) THEN $ print, ' --> Modify field using : ', 'mult/add=' , $ fldatt.mult, fldatt.add, ' to obtain : ', fldatt.unit END 0: BEGIN print, ' '+cmd.var, ' plot attributes ('+plttyp+') : min/max=',minc, maxc IF ( fldatt.mult NE 1.0 OR fldatt.add NE 0.0 ) THEN $ print, ' --> Modify field using : ', 'mult/add=' , $ fldatt.mult, fldatt.add, ' to obtain : ', fldatt.unit print, ' '+cmd.var2, ' plot attributes ('+plttyp+') : min/max=',minc2, maxc2 IF ( fldatt2.mult NE 1.0 OR fldatt2.add NE 0.0 ) THEN $ print, ' --> Modify field using : ', 'mult/add=' , $ fldatt2.mult, fldatt2.add, ' to obtain : ', fldatt2.unit END ELSE: ENDCASE print, '' ENDIF IF ( fldatt.mult EQ -1.0) THEN BEGIN temp_m = minc minc = -maxc maxc = -temp_m ENDIF IF finite(fldatt.int) EQ 0 THEN BEGIN intcmd = '' colbarfor = '' fmt = '(f6.1)' ENDIF ELSE BEGIN intcmd = ',int = '+string(fldatt.int) IF fldatt.int LT 0.1 THEN BEGIN fmt = '(f6.2)' ENDIF ELSE IF fldatt.int LT 1 THEN BEGIN fmt = '(f5.1)' ENDIF ELSE BEGIN fmt = '(f5.0)' ENDELSE IF long(fldatt.int) NE 0 THEN BEGIN IF fldatt.int/long(fldatt.int) NE 1 THEN fmt = '(f5.1)' ENDIF colbarfor = ', format = '''+fmt+'''' ENDELSE ; if window > 1 or overlay > 1 noerase IF win[2] NE 1 OR iover GT 1 THEN BEGIN noerase = 1 ENDIF ELSE BEGIN noerase = 0 ENDELSE ; choose window number CASE multi_win OF 1: window_number = ', window='+string(iplot) ELSE: window_number = '' ENDCASE ; color label control labstr = '' IF iover EQ 1 THEN BEGIN nofill = 1-shading ENDIF ELSE BEGIN nofill = 1 ENDELSE nocoltxt = '' IF nofill EQ 1 THEN nocoltxt = ',NOCOULEUR=1' IF dimplot NE 2 THEN nofill = 1 nocolorbar = 0 IF col_palette EQ 'no' THEN nocolorbar = 1 IF pal_type EQ '2dom' THEN nocolorbar = 1 IF iover GT 1 THEN nocolorbar = 1 IF dimplot NE 2 THEN nocolorbar = 1 IF cmd.out NE 'cdf' THEN BEGIN readpal = 0 IF nofill EQ 0 AND dimplot GT 1 THEN readpal = 1 IF waveplt EQ 1 THEN readpal = 2 IF field.origin EQ 'div' AND plttyp NE 'pltt' THEN readpal = 3 ;; Necessary for correlation plots (overlay of 2D fields) IF dimplot GT 1 AND nover GT 1 AND iover GT 1 THEN readpal = 1 IF readpal GE 1 THEN BEGIN lec_pal_gmt, cmd.var, c_anot_str, fmt, found, readpal colbarfor = ', format = '''+fmt+'''' IF found EQ 1 THEN BEGIN labstr = ',label=3' ; if mincmd/naxcmd not defined then use one from palette ; if define take min of 2 mins and max of 2 maxs min_palette = min(levels_gmt) max_palette = max(levels_gmt) IF finite(minc) NE 0 THEN minc = min([minc, min_palette]) ELSE minc = min_palette IF finite(maxc) NE 0 THEN maxc = max([maxc, max_palette]) ELSE maxc = max_palette print, ' '+cmd.var, ' plot attributes ('+plttyp+') modified via color palette: min/max=', $ minc, maxc ENDIF ENDIF ENDIF colbar = colbarfor mincmd = '' maxcmd = '' IF finite(minc) NE 0 THEN BEGIN mincmd = ','+string(minc) maxcmd = ','+string(maxc) ENDIF ; contour annotation + other color bar controls IF n_elements(c_anot_str) EQ 0 THEN BEGIN c_anot = '' ENDIF ELSE BEGIN c_anot = ',c_annotation=c_anot_str' ; sampling of colorbar annotation IF col_palette EQ 'yes' THEN BEGIN sampling = win[0]*win[1] ; idx_cb = lindgen((ncont_gmt-2)/sampling)*sampling+1 idx_cb0 = lindgen(ncont_gmt-1)+1 idx_cb = lindgen(ncont_gmt-2)+2 idx_cb1 = lindgen(ncont_gmt-2)+1 colbar = colbar+',discret=coul_gmt[idx_cb1],div=n_elements(idx_cb0)-1,cb_label=levels_gmt[idx_cb0]' ; c_anot = ',c_annotation=levels_gmt[idx_cb0]' ENDIF ENDELSE ; ; vertical grid type ; IF debug_w THEN print, 'splot=', splot CASE mesh_type OF 'oce': BEGIN IF splot EQ 1 THEN BEGIN type_yz = ',type_yz =''sigma''' zoom_txt = ', zoom='+string(sig_max) ENDIF ELSE BEGIN type_yz = ',type_yz=''m''' zoom_txt = ', zoom='+string(zoom_z) ENDELSE END 'atm': BEGIN & type_yz = ',type_yz=''hPa''' & zoom_txt = ', zoom='+string(hpa_max) & END ELSE:type_yz = '' ENDCASE ; continents color, thickness ; c_cont = 100 ; continent fill ; real continent drawn strcont ='' IF mesh_type NE 'oce' THEN BEGIN IF cont_fill EQ 0 THEN strcont = ',/CONT_NOFILL, nite=0' IF cont_real GE 1 THEN strcont = strcont+', REALCONT = '+string(strtrim(cont_real, 2))+', COAST_THICK = 2' ENDIF ; calendar type cal_typ = '' IF calendar_type GT 1 THEN cal_typ = ',ndayspm=calendar_type' ; titles options ; By default, SAXO puts a title and a subtitle for the first plot ; Then, if you overlay another field, title and subtitle are set to '' CASE title_type OF "T": tit_str = ',subtitle='''','+marge_option "S": tit_str = ',title='''','+marge_option "TS": tit_str = ','+marge_option "off":tit_str = ',title='''',subtitle='''','+marge_option ENDCASE ; fill_space option filltxt = '' IF fill_space EQ 1 THEN filltxt = ',/rempli' ; common command string to plots com_strplt = ',petit = ['+string(win[0])+','+string(win[1])+','+string(win[2])+']'+filltxt+nocoltxt+', LANDSCAPE = '+string(landscape)+', NOCOLORBAR = '+string(nocolorbar)+', NOERASE = '+string(noerase)+look+labstr+c_anot+colbar+cal_typ+',cont_thick=2'+type_yz+window_number+tit_str+xlogax+', /MEMEINDICES' ; add contour_options not overlay IF iover EQ 1 THEN com_strplt = com_strplt+contour_options ; run specific fixes IF (strpos(cmd.exp, 'MIMR') NE -1 OR strpos(cmd.exp, 'MIHR') NE -1) AND cmd.var EQ 'tauu' AND strpos(cmd.timave, '1m') NE -1 AND strpos(cmd.timave, '1mm') EQ -1 THEN BEGIN fld = -fld print, ' *** WARNING: Multiply tauu by -1 ', cmd.exp, cmd.var, cmd.timave ENDIF ; ============================ ; 3. make output (data/plot) ; ============================ ; CASE cmd.out OF 'data': write_data = long((iplot-1)*10+win[2]) ; write 1D data ascii in trends.pro 'tcdf': write_data = -long((iplot-1)*10+win[2]) ; write 1D data netcdf in trends.pro 'cdf': write_data = 0 ELSE: write_data = 0 ENDCASE ; make output CASE cmd.out OF 'tv': BEGIN ; $$$$$$$$$$$$$$$$$$$$$$ tvnplot $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ erase CASE (size(fld))[0] OF 2: fldtv = fld 3: BEGIN niveau = xquestion('3d data : which level ?', '1', /chkwidget) fldtv = fld[*, *, niveau] END 4: BEGIN niveau = xquestion('3d data : which level ?', '1', /chkwidget) timeplot = xquestion('4d data : which date ?', '1', /chkwidget) fldtv = fld[*, *, niveau, timeplot] END ELSE: ENDCASE CASE cmd.grid OF ; 'T': tvnplot, fldtv*tmask 'T': tvplus, fldtv, min = fldatt.min, max = fldatt.max 'U': tvplus, fldtv, min = fldatt.min, max = fldatt.max 'V': tvplus, fldtv, min = fldatt.min, max = fldatt.max 'W': tvplus, fldtv, min = fldatt.min, max = fldatt.max ELSE: tvplus, fldtv, min = fldatt.min, max = fldatt.max ENDCASE END 'cdf': BEGIN ; $$$$$$$$$$$$$$$$$$$$$$ netCDF output $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ ; build structures containing grids cdf_description = nc_build(cmd, field, field.direc, vargrid) ; build netCDF file (here if whole data) IF hotyp EQ '-' OR hotyp EQ 'xyt' THEN BEGIN fldcdf = {data:field.data, units:field.units, short_name:cmd.var, long_name:field.legend, missing_value:valmask, direc:field.direc} file_out = cmd.var+'_'+strtrim(string(FORMAT = '(I2.2)', iplot), 2)+'.nc' pltcmd ='nc_put, fldcdf, file_out, ncdf_db = hom_idl'+cdf_description printf, nulhis, strcompress(pltcmd) res = execute(pltcmd) ENDIF END 'tcdf': BEGIN ; $$$$$$$$$$$$$$$$$$$$$$ 1D netCDF output $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ mask_z, fld, cmd, boite_plt1d, dimplot fld = checkfield(fld, 'pltt', type = datyp.hotyp, boite = boite_plt1d) IF debug_w THEN print, 'size(fld)', size(fld) IF debug_w THEN print, 'boite_plt1d', boite_plt1d IF cmd.trend GT 0 THEN BEGIN fld = trends(fld, cmd.trend, datyp.hotyp) add_txt = trd_txt ENDIF ELSE print, 'You need to have a trend to make 1D netcdf output' END ELSE: BEGIN ; $$$$$$$$$$$$$$$$$$$$$$$$$$ make plot $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ CASE plttyp OF 'plt': BEGIN ; ; Map plot ; -------- ; print, ' '+cmd.var+' data to plot min and max : ', $ min(fld(where (fld NE -valmask)), /nan), max(fld(where (fld NE valmask)), /nan) print, '' map_cmd = '' niveau = 1 IF vecplot EQ 0 THEN BEGIN ; eddy field ? (if so, remove zonal mean) IF strpos(cmd.plt, '@z') GT 1 THEN BEGIN fldzm = moyenne(fld, 'x', boite=[20,380,-90,90], /NAN) fldzm = replicate(1, (size(fld))(1))#fldzm fld = fld-fldzm filleg = ' Eddy ' ENDIF ELSE BEGIN filleg = '' ENDELSE IF hotyp EQ 'xyt' THEN date_txt = 'monthly' pltcmd = 'plt,{a:fld, d:'''+date_txt+''', n:'''+filleg+varlegend+' '+legbox+legsurf+''', u:'''+field.units+''', g:'''+vargrid+'''}'+mincmd+maxcmd+intcmd+com_strplt+strcont+map_cmd ;; For 2D plots with overlays (correlation and proba for instance) ;; Scratch the titles and define contours style for the second plot IF nover GT 1 AND iover GT 1 THEN BEGIN contour_style = ',c_thick=2,c_linestyle=2,c_charthick=2' CASE title_type OF 'T': pltcmd = pltcmd+',title='''''+contour_style 'S': pltcmd = pltcmd+',subtitle='''''+contour_style 'TS': pltcmd = pltcmd+',title='''',subtitle='''''+contour_style 'off': pltcmd = pltcmd+contour_style END ENDIF ENDIF ELSE BEGIN fld = {a:fld, g:strmid(cmd.grid, 0, 1)} fld2 = {a:fld2, g:strmid(cmd.grid, 1, 1)} fldv = {data1: fld, data2: fld2} vargrid = strmid(cmd.grid, 2, 1) pltcmd = 'ajoutvect,fldv'+vect_sample+com_strplt+strcont+map_cmd ENDELSE IF debug_w THEN print, strcompress(pltcmd) printf, nulhis, strcompress(pltcmd) res = execute(pltcmd(0)) END 'pltz': BEGIN ; ; Vertical section/mean plot ; -------------------------- ; g = gphit t = tmask ; if already 1D, reform fld sz = size(fld) IF sz[0] EQ 2 THEN BEGIN IF sz[1] EQ 1 OR sz[2] EQ 1 THEN fld = reform(fld, sz[1]*sz[2]) ENDIF IF sz[0] EQ 3 THEN BEGIN IF sz[1] EQ 1 THEN fld = reform(fld, sz[2], sz[3]) IF sz[2] EQ 1 THEN fld = reform(fld, sz[1], sz[3]) IF sz[3] EQ 1 THEN fld = reform(fld, sz[1], sz[2]) ENDIF ; mask fld mask_z, fld, cmd, boite_pltz, dimplot, legz printf, nulhis, ' boite_pltz = ', boite_pltz print, ' '+cmd.var+' data to plot min and max : ', $ min(fld(where (fld NE valmask)), /NAN), max(fld(where (fld NE valmask)), /NAN) print, '' IF vecplot EQ 0 THEN BEGIN IF cmd.var EQ 'vozonbsf' THEN xindstr = ', /xindex' ELSE xindstr = '' pltcmd = 'pltz,{a:fld, d:'''+date_txt+''', n:'''+varlegend+' '+legz+''', u:'''+field.units+''', g:'''+vargrid+'''}'+mincmd+maxcmd+intcmd+',/'+pltztyp+', boite=boite_pltz'+zoom_txt+com_strplt+xindstr ENDIF ELSE BEGIN fld = {a:fld, g:strmid(cmd.grid, 0, 1)} fld2 = {a:fld2, g:strmid(cmd.grid, 1, 1)} fldv = {data1: fld, data2: fld2} vargrid = strmid(cmd.grid, 2, 1) pltcmd = 'ajoutvectz,fldv'+com_strplt+strcont+',type='''+strmid(cmd.plt, 0, 2)+''', boite=boite_pltz' ENDELSE IF debug_w THEN print, pltcmd printf, nulhis, strcompress(pltcmd) res = execute(pltcmd(0)) ; overlay bowl if sig_bowl=1 IF sig_bowl EQ 1 THEN begin ; define line color, thickness and type iover = 2 overc = overlay_type(iover, dimplot) ; type of latitude axis IF strmid(cmd.plt, 0, 1) EQ 'y' AND lat_axis EQ 'sin' THEN $ sin = ',/sin' ELSE sin = '' plt1dtyp = strmid(pltztyp, 0, 1) noerase = 1 com_strplt = ', NOERASE = '+string(noerase) pltcmd = 'plt1d,field.bowl,/'+plt1dtyp+', boite=boite_pltz'+overc+sin+com_strplt IF debug_w THEN print, pltcmd res = execute(pltcmd(0)) ENDIF gphit = g tmask = t END 'plt1d': BEGIN ; ; 1D, x, y, z plot ; ---------------- ; g = gphit t = tmask ; if already 1D, reform fld sz = size(fld) IF sz[0] EQ 2 THEN BEGIN IF sz[1] EQ 1 OR sz[2] EQ 1 THEN fld = reform(fld, sz[1]*sz[2]) ENDIF IF sz[0] EQ 3 THEN BEGIN IF sz[1] EQ 1 THEN fld = reform(fld, sz[2], sz[3]) IF sz[2] EQ 1 THEN fld = reform(fld, sz[1], sz[3]) IF sz[3] EQ 1 THEN fld = reform(fld, sz[1], sz[2]) ENDIF ; mask fld mask_z, fld, cmd, boite_plt1d, dimplot, legz IF debug_w THEN print, boite_plt1d niveau = 1 ; define line color, thickness and type overc = overlay_type(iover, dimplot) ; type of latitude axis IF strmid(cmd.plt, 0, 1) EQ 'y' AND lat_axis EQ 'sin' THEN $ sin = ',/sin' ELSE sin = '' print, ' '+cmd.var+' data to plot min and max : ', $ min(fld(where (fld NE valmask))), max(fld(where (fld NE valmask))) print, '' take_log = ',take_log=0' IF cmd.var EQ 'vosigvol' AND n_elements(str_sep(cmd.exp,'-')) EQ 1 THEN BEGIN ; don't take log if difference plot take_log = ',take_log=1' ENDIF pltcmd = 'plt1d,{a:fld, d:'''+date_txt+''', n:'''+varlegend+' '+legz+''', u:'''+field.units+''', g:'''+vargrid+'''},'''+plt1dtyp+''''+mincmd+maxcmd+intcmd+', boite=boite_plt1d'+overc+sin+com_strplt+take_log printf, nulhis, 'boite_plt1d=',boite_plt1d printf, nulhis, strcompress(pltcmd) IF debug_w THEN print, pltcmd res = execute(pltcmd(0)) @legend_overlay gphit = g tmask = t END 'pltt': BEGIN ; ; hovmoeller ; ----------- ; g = gphit t = tmask ; mask fld mask_z, fld, cmd, boite_pltt, dimplot, legz ; define line color, thickness and type overc = overlay_type(iover, dimplot) ; define additional text for pltt @add_txt_pltt ; repeat for seasonal mean CASE cmd.timave OF '1mm': rep_txt = ',repeat_c=nb_cycles' ELSE: rep_txt = '' ENDCASE print, ' '+cmd.var+' data to plot min and max : ', $ min(fld(where (fld NE valmask))), max(fld(where (fld NE valmask))) print, '' vardate = 'toto' grille, mask, glam, gphi, gdep, nx, ny,nz IF st_rms EQ 1 THEN BEGIN ; time series of rms deviation on a sigma surface @densit_pltmap_plot ENDIF ELSE BEGIN pltcmd = 'pltt,{a:fld, n:'''+date_txt+' '+varlegend+' '+legbox+''', u:'''+field.units+filleg+''', g:'''+vargrid+'''}, '''+hotyp+''''+mincmd+maxcmd+intcmd+',boite=boite_pltt'+com_strplt+overc+filtxt+',TREND_TYPE='+cmd.trend+rep_txt ; +',/integration' ENDELSE printf, nulhis, 'boite_pltt=',boite_pltt printf, nulhis, strcompress(pltcmd) IF debug_w THEN print, pltcmd res = execute(pltcmd(0)) IF iover EQ 4 THEN print, 'There might be a pb with the legends !' ; legend if plot=t IF hotyp EQ 't' THEN BEGIN ; positions of the legend depend on nb_cycles ; it is different from 1 only with 1mm case IF cmd.timave NE '1mm' THEN nb_cycles = 1 @legend_overlay ENDIF gphit = g tmask = t END 'yfx': BEGIN ; ; Scatter plot y=f(x) ; ------------------- ; t = tmask boxx = '' boxy = '' add_txt = '' ; arrange/average data if needed ; space and time plots fld IF strpos(cmd.plt, 'xy_') NE -1 THEN BEGIN mask_z, fld, cmd, boite_plt, dimplot, boxx boxx = ' ['+boxx+']' print, ' Averaging (space) '+cmd.var+' in'+boxx fld = checkfield(fld, 'plt', type = 'xy', boite = boite_plt) ENDIF IF datyp.hotyp NE '-' THEN BEGIN mask_z, fld, cmd, boite_plt1d, dimplot, boxx boxx = ' ['+boxx+']' vargrid = vargrid1 IF debug_w THEN print, 'domdef and vargrid : ', boite_plt1d[0:3], vargrid domdef, boite_plt1d[0:3] print, ' Averaging (time serie plot) cmd.var '+cmd.var+' in'+boxx IF debug_w THEN print, 'fld size', size(fld) fld = checkfield(fld, 'plt1d', type = datyp.hotyp, boite = boite_plt1d) IF cmd.trend GT 0 THEN BEGIN fld_flag = 1 fld = trends(fld, cmd.trend, datyp.hotyp) add_txt = trd_txt ENDIF ENDIF ; space and time plots fld2 IF sw_diffg EQ 1 THEN BEGIN jptb = jpt def_grid, cmd2 IF debug_w THEN print, 'def_grid, cmd2 dans plt_map ' jpt = jptb ENDIF IF strpos(cmd2.plt, 'xy_') NE -1 THEN BEGIN mask_z, fld2, cmd2, boite_plt2, dimplot, boxy boxy = ' ['+boxy+']' print, ' Averaging (space) '+cmd2.var+' in'+boxy fld2 = checkfield(fld2, 'plt', type = 'xy', boite = boite_plt2) ENDIF IF datyp2.hotyp NE '-' THEN BEGIN run_stddev = 0 IF strpos(cmd2.plt, '@r') GE 1 THEN BEGIN run_stddev = strmid(cmd2.plt, strpos(cmd2.plt, '@r')+2, strlen(cmd2.plt)) cmd2.plt = strmid(cmd2.plt, 0, strpos(cmd2.plt, '@r')) ENDIF mask_z, fld2, cmd2, boite_plt1d2, dimplot, boxy boxy = ' ['+boxy+']' vargrid = vargrid2 IF debug_w THEN print, 'domdef and vargrid : ', boite_plt1d2[0:3], vargrid domdef, boite_plt1d2[0:3] print, ' Averaging (time serie plot) cmd2.var '+cmd2.var+' in'+boxy IF debug_w THEN print, 'fld size', size(fld2) fld2 = checkfield(fld2, 'plt1d', type = datyp2.hotyp, direc = 'xy', boite = boite_plt1d2, /timearray) IF cmd2.trend GT 0 THEN BEGIN fld_flag = 2 fld2 = trends(fld2, cmd2.trend, datyp2.hotyp) add_txt = trd_txt ENDIF ENDIF IF sw_diffg EQ 1 THEN BEGIN sw_diffg = 0 jptb = jpt def_grid, cmd IF debug_w THEN print, 'def_grid, cmd dans plt_map ' jpt = jptb ENDIF coeff = linfit(fld2, fld, CHISQ = linerr) print, ' Linfit coef + error (full serie)=', coeff(0), coeff(1), linerr ; plot domain boxyfx = def_box(cmd.plt, dimplot, legz, time_stride) ; define line color and type overc = overlay_type(iover, dimplot) vardate = date_txt varname = varlegend varunit = ' '+cmd.var+boxx+' =f('+cmd.var2+boxy+') '+add_txt pltcmdstd = 'pltsc,fld,fld2,minc,maxc,minc2,maxc2,varlegend2, boite=boxyfx'+overc+com_strplt IF hotyp EQ 't' THEN BEGIN ; time scatter plot IF mean_sc_only EQ 0 OR mean_sc_only EQ 4 THEN BEGIN ; number of time colors IF strpos(symbol_families, 'x') NE -1 THEN BEGIN coding = str_sep(symbol_families, 'x') ncolors = long(coding(0)) ntimes = long(coding(1)) ENDIF ELSE BEGIN ncolors = long(symbol_families) ntimes = 1 ENDELSE IF ncolors GE 2 THEN BEGIN ; time color coding ic = 0 jpto = jpt jpt = jpt/ncolors lincoef = fltarr(ncolors) linerro = fltarr(ncolors) linprob = fltarr(ncolors) WHILE ic LE ncolors-1 DO BEGIN idx0 = (floor(findgen(jpto/(ncolors*ntimes)))*ncolors*ntimes)+ic*(ntimes) idx = idx0 jl = 1 WHILE jl LE ntimes-1 DO BEGIN idx = [idx, idx0+jl] jl = jl+1 ENDWHILE fldp = fld(idx) fldp2 = fld2(idx) coeff = linfit(fldp2(where (fldp2(*) GT 0.)), fldp(where (fldp2(*) GT 0.)), CHISQ = errlin, PROB = prblin) linerro(ic) = errlin linprob(ic) = prblin lincoef(ic) = coeff(1) ; print, ' Period '+strtrim(string(ic), 2)+' Linear fit slope, chisq, prob =', coeff(1)*1000., errlin*1000., prblin IF mean_sc_only EQ 0 THEN BEGIN pltcmd = 'pltsc,fldp,fldp2,minc,maxc,minc2,maxc2,varlegend2, boite=boxyfx'+com_strplt+',ov1d='+string(ic)+',COLOR='+string(symbol_color(ic))+', STY1D='+string(symbol_style(ic)-1) printf, nulhis, strcompress(pltcmd) IF debug_w THEN print, pltcmd res = execute(pltcmd) ENDIF IF mean_sc_only EQ 4 AND ic EQ ncolors-1 THEN BEGIN print, ' Linfit slope + error (June-Nov)=', mean(lincoef(5:10)), mean(linerro(5:10)) print, ' Linfit slope + error (Dec-May)=', mean([lincoef(0:4), lincoef(11)]), mean([linerro(0:4), linerro(11)]) vardate = 'toto' varname = 'Coupling strength '+cmd.date1+' '+cmd.spec jpt = ncolors time = lindgen(12)*30*1+ $ julday(1,1,01, ndayspm = calendar_type) hotyp = 't' minmax = ',0,25' pltcmd = 'pltt,lincoef*1000., '''+hotyp+''''+minmax+com_strplt+',ov1d=0,COLOR=1, thick=4, STY1D=0' printf, nulhis, strcompress(pltcmd) IF debug_w THEN print, pltcmd res = execute(pltcmd) pltcmd = 'pltt,(lincoef-linerro)*1000., '''+hotyp+''''+minmax+com_strplt+',ov1d=1,COLOR=1, thick=1, STY1D=0' printf, nulhis, strcompress(pltcmd) IF debug_w THEN print, pltcmd res = execute(pltcmd) pltcmd = 'pltt,(lincoef+linerro)*1000., '''+hotyp+''''+minmax+com_strplt+',ov1d=1,COLOR=1, thick=1, STY1D=0' printf, nulhis, strcompress(pltcmd) IF debug_w THEN print, pltcmd res = execute(pltcmd) ENDIF ic = ic + 1 ENDWHILE ENDIF ELSE BEGIN ; one color IF debug_w THEN print, pltcmdstd res = execute(pltcmdstd(0)) ENDELSE ENDIF ENDIF ELSE BEGIN ; standard scatter plot IF debug_w THEN print, pltcmdstd res = execute(pltcmdstd(0)) ENDELSE ; add mean SC in 4x3 plots IF symbol_families EQ '4x3' THEN BEGIN IF cmd.trend EQ '412' THEN BEGIN fldrem_t1 = fldrem_t1-mean(fldrem_t1) fldrem_t2 = fldrem_t2-mean(fldrem_t2) ENDIF ELSE BEGIN running = 12L lenght = (size(fld))[1] fldrem_t1 = fltarr(running) fldrem_t2 = fltarr(running) FOR t1 = 0, running-1 DO BEGIN fldrem_t1(t1) = mean(fld(long(findgen(lenght/running)*running+t1))) fldrem_t2(t1) = mean(fld2(long(findgen(lenght/running)*running+t1))) ENDFOR ENDELSE print, ' Seasonal cycle var/stddev fld1 ', (moment(fldrem_t1))[1], sqrt((moment(fldrem_t1))[1]) print, ' Seasonal cycle var/stddev fld2 ', (moment(fldrem_t2))[1], sqrt((moment(fldrem_t2))[1]) fldrem_t1 = [fldrem_t1, fldrem_t1(0)] fldrem_t2 = [fldrem_t2, fldrem_t2(0)] sw_ov1d = mean_sc_only IF mean_sc_only EQ 2 AND cmd.trend EQ '412' THEN BEGIN ; SC of std dev stat = spectra(fld, jpt/12, 20, 15, 0) stat2 = spectra(fld2, jpt/12, 20, 15, 0) print, ' Stddev of monthly stddevs fld', sqrt((moment(stat.sc_std-(moment(stat.sc_std))[0]))[1]) print, ' Stddev of monthly stddevs fld2', sqrt((moment(stat2.sc_std-(moment(stat2.sc_std))[0]))[1]) stdsc = [stat.sc_std, stat.sc_std(0)] stdsc2 = [stat2.sc_std, stat2.sc_std(0)] pltcmd = 'pltsc,stdsc,stdsc2,0. ,maxc,0.,maxc2,varlegend2, boite=boxyfx'+com_strplt+',ov1d=1-min(1,sw_ov1d), STY1D=-1, THICKN=5, SYMSIZE=1, FRACTION=fraction' printf, nulhis, strcompress(pltcmd) IF debug_w THEN print, pltcmd res = execute(pltcmd(0)) pltcmd = 'xyouts,stdsc2-(maxc2)/(.3*win(0)*win(1)),stdsc,string(long(findgen(12))+1),charsize=1.3,alignment=0.5' printf, nulhis, strcompress(pltcmd) IF debug_w THEN print, pltcmd res = execute(pltcmd(0)) ENDIF ELSE BEGIN ov1d_val = 0 IF mean_sc_only EQ 0 THEN ov1d_val = 1 pltcmd = 'pltsc,fldrem_t1,fldrem_t2,minc,maxc,minc2,maxc2,varlegend2, boite=boxyfx'+com_strplt+',ov1d=ov1d_val, STY1D=-1, THICKN=5, SYMSIZE=1, FRACTION=fraction' printf, nulhis, strcompress(pltcmd) IF debug_w THEN print, pltcmd res = execute(pltcmd(0)) IF mean_sc_only EQ 1 THEN BEGIN ; add month info win(0)*win(1) pltcmd = 'xyouts,fldrem_t2-(maxc2-minc2)/(0.5*win(0)*win(1)),fldrem_t1,string(long(findgen(12))+1),charsize=1.3,alignment=0.5' printf, nulhis, strcompress(pltcmd) IF debug_w THEN print, pltcmd res = execute(pltcmd(0)) ENDIF ENDELSE ENDIF ; overplot sigma contours in T-S plane IF (cmd.var EQ 'sosstsst' and cmd.var2 EQ 'sosaline') OR (cmd.var EQ 'votemper' and cmd.var2 EQ 'vosaline') THEN BEGIN IF noerase EQ 0 THEN BEGIN min_t = minc-2.0 max_t = maxc+2.0 min_s = minc2-2.0 max_s = maxc2+2.0 delta_t = 0.5 delta_s = 0.05 t_vec = findgen((max_t-min_t)/delta_t+1)*delta_t+min_t s_vec = findgen((max_s-min_s)/delta_s+1)*delta_s+min_s t_ = t ; since t is already defined t = t_vec FOR i = 1, (max_s-min_s)/delta_s DO t = [[t], [t_vec]] s = s_vec FOR i = 1, (max_t-min_t)/delta_t DO s = [[s], [s_vec]] s = transpose(s) sr=sqrt(abs(s)) r1=((((6.536332E-9*t-1.120083E-6)*t+1.001685E-4)*t $ -9.095290E-3)*t+6.793952E-2)*t+999.842594 r2=(((5.3875E-9*t-8.2467E-7)*t+7.6438E-5)*t-4.0899E-3)*t+8.24493E-1 r3=(-1.6546E-6*t+1.0227E-4)*t-5.72466E-3 rhopn = ( ( 4.8314E-4*s + r3*sr +r2)*s +r1) -1000.0 min_sig = round(min(rhopn)+0.5)-1. max_sig = round(max(rhopn)-0.5)+1. delta_sig = 1. ; define contour interval here! levels = findgen((max_sig-min_sig)/delta_sig+1)*delta_sig+min_sig labels = levels labels(*) = 1 contour, rhopn, s, t, /overplot, levels = levels, c_labels = labels, c_linestyle = 0 t = t_ ENDIF ENDIF printf, nulhis, 'boite_yfx=',boxyfx printf, nulhis, strcompress(pltcmd) tmask = t ; force read grid for next window cmd_prev.grid = 'toto' END 'spec': BEGIN ; ; Spectrum plot y=fft(x) ; --------------------- ; t = tmask ; plot domain mask_z, fld, cmd, boite_pltspec, dimplot, legspec ; define line color and type overc = overlay_type(iover, dimplot) vardate = date_txt varname = varlegend varunit = cmd.exp+' '+cmd.timave+' '+cmd.date1+'-'+cmd.spec+' '+varlegend+' '+legspec ; make mean in box fld = checkfield(fld, 'pltt', type = 't', boite = boite_pltspec, _extra = ex) ; apply trends IF long(cmd.trend) GT 0 THEN fld = trends(fld, long(cmd.trend), 't') ; make spectrum time_interval = time[1]-time[0] print, ' Max plot spectrum in plt_def : days/month/year ', max_spec, max_spec/30, max_spec/360 print, ' lenght of time serie (years) = ', long(time(jpt-1)-time(0)+time_interval)/360 print, ' spectrum window (years)/chunks = ', spec_win/360, long(time(jpt-1)-time(0)+time_interval)/spec_win print, ' ' ; number of time windows¨ nt_win = long((time(jpt-1)-time(0)+time_interval)/spec_win) idx_win_size = long(spec_win/time_interval) IF idx_win_size*nt_win NE jpt THEN BEGIN print, ' *** Warning : spec_win must divide lenght of time serie ', idx_win_size, jpt return ENDIF ; mean of idx_win_size chunks spect = spectrum(fld[0:idx_win_size-1], time_interval) FOR it = 2, nt_win DO BEGIN idx1 = idx_win_size*(it-1) idx2 = idx1 + idx_win_size - 1 specc = spectrum(fld[idx1:idx2], time_interval) spect = spect + specc ENDFOR spect = spect/nt_win ; build new time array idx = where (spect[0, *] LE max_spec) jpt = n_elements(idx) fld = reverse(reform(spect[1, idx], jpt)) time = findgen(jpt) time = reverse(reform(spect[0, idx], jpt)) IF max(time) GT 20*360 THEN BEGIN time = time/360 xtitle = 'Years' ENDIF ELSE IF max(time) GT 5*360 THEN BEGIN time = time/30 xtitle = 'Months' ENDIF ELSE xtitle = 'Days' ; time/space filter ? IF strpos(cmd.plt, '@f') GT 1 THEN BEGIN filter = long(strmid(cmd.plt, strpos(cmd.plt, '@f')+3, strlen(cmd.plt)-strpos(cmd.plt, '@f')-3)) fildirec = strmid(cmd.plt, strpos(cmd.plt, '@f')+2, 1) IF fildirec EQ 't' THEN BEGIN xtitle = xtitle+' (filter '+strtrim(string(filter), 2)+' warning : discrete filter)' fld = smooth(fld, filter) fld[0:filter/2-1] = 0. fld[(size(fld))[1]-filter/2-1:(size(fld))[1]-1] = 0. ENDIF ENDIF ; plot min/max min1 = min(fld) max1 = max(fld) min2 = 0 max2 = max(time) !x.range = [min2-abs(max2-min2)/5.,max2+abs(max2-min2)/5.] !y.range = [min1-abs(max1-min1)/5.,max1+abs(max1-min1)/5.] ; plot ytitle = 'Power spectrum (window='+strtrim(string(spec_win/360), 2)+'y)' pltcmd = 'splot,time,fld,xstyle=1,ystyle=1,title=varunit,xtitle=xtitle,ytitle=ytitle'+overc+com_strplt printf, nulhis, 'boite_pltspec=',boite_pltspec printf, nulhis, strcompress(pltcmd) IF debug_w THEN print, pltcmd res = execute(pltcmd(0)) tmask = t END 'wavelet': BEGIN ; ; Wavelet plot ; ------------ ; t = tmask ; plot domain mask_z, fld, cmd, boite_pltspec, dimplot, legspec ; define line color and type vardate = date_txt varname = varlegend + ' '+date_txt varunit = cmd.exp+' '+cmd.timave+' '+cmd.date1+'-'+cmd.spec+' '+varlegend+' '+legspec ; make mean in box fld = checkfield(fld, 'pltt', type = 't', boite = boite_pltspec, _extra = ex) ; apply trends IF long(cmd.trend) GT 0 THEN fld = trends(fld, long(cmd.trend), 't') ; make spectrum time_interval = time[1]-time[0] print, ' Max plot spectrum in plt_def : days/month/year ', max_spec, max_spec/30, max_spec/360 print, ' lenght of time serie (years) = ', long(time(jpt-1)-time(0)+time_interval)/360 print, ' ' ; make wavelet wave = wavelet(fld,time_interval,period=period,coi=coi,/pad,signif=signif) power = abs(wave)^2 nscale = n_elements(period) period = period/365 ; to have period in years ; compute significance signif = rebin(transpose(signif),jpt,nscale) ; plot wavelet mincmd = ','+string(min(power)) maxcmd = ','+string(max(power)) boite_pltspec = [0, 0, min(period), max(period)] ; domdef, boite_pltspec lat1r = lat1 lat2r = lat2 lat1 = 0 lat2 = max_spec/360 key_onearth = 0 pltcmd = 'plttg,power,period'+mincmd+maxcmd+intcmd+',boite=boite_pltspec,reverse_y=1,nocontour=1'+com_strplt printf, nulhis, 'boite_pltspec=',boite_pltspec printf, nulhis, strcompress(pltcmd) IF debug_w THEN print, pltcmd res = execute(pltcmd(0)) contour,abs(wave)^2/signif,time,period, /overplot,level=1.0,c_annot='95%' plot,time,coi/365, noclip = 0, /noerase tmask = t lat1 = lat1r lat2 = lat2r key_onearth = 1 END ELSE: BEGIN print, ' unknown projection plot ', cmd.plt, ' / plot type = ', plttyp END ENDCASE END ENDCASE fld_prev = cmd.var ; ; reset incompatible options of plt_def cmd.trend = trend_typp ; ; reset vertical grid after density bining IF splot EQ 1 THEN BEGIN izminmesh = izminmesh_b izmaxmesh = izmaxmesh_b jpk = jpk_b jpkglo = jpkglo_b gdept = gdept_b gdepw = gdepw_b e3t = e3t_b e3w = e3w_b tmask = tmask_b ENDIF END