; pro defzonaldata,line,data,linedata field = line.field fieldlevs = *field.levs nfieldlevs = n_elements(fieldlevs) nlon = field.nlon & nlat = field.nlat & nlev = field.nlev lons = *line.lons iLat = line.ilat ilon = line.ilon iLev = line.ilev ; ; For WACCM histories, when picking vertical level, need to do special handling of ; fields on interface levels ; IF line.ftype EQ 'WACCM' THEN BEGIN iFieldLev = 0 IF iLev GT 0 THEN BEGIN nLevs = N_ELEMENTS(fieldlevs) newLevs = *line.levs newNLevs = N_ELEMENTS(newLevs) newLevSelected = newLevs[iLev] iFieldLev = IXFIND_NEAREST(fieldlevs, newlevSelected) levSelected = fieldlevs[iFieldLev] IF iFieldLev GT 0 THEN levSelectedLower = fieldlevs[iFieldLev-1] IF iFieldLev LT nLevs-1 THEN levSelectedUpper = fieldlevs[iFieldLev+1] IF newLevSelected LT levSelected AND iFieldLev GT 0 THEN BEGIN levDiffLower = levSelected - newLevSelected levDiffUpper = newLevSelected - levSelectedLower levDiffAll = levSelected - levSelectedLower levRatLower = levDiffLower / levDiffAll levRatUpper = levDiffUpper / levDiffAll ENDIF ELSE IF newLevSelected GT levSelected AND iFieldLev LT nLevs-1 THEN BEGIN levDiffLower = newLevSelected - levSelected levDiffUpper = levSelectedUpper - newLevSelected levDiffAll = levSelectedUpper - levSelected levRatLower = levDiffLower / levDiffAll levRatUpper = levDiffUpper / levDiffAll ENDIF ENDIF ENDIF linedata = fltarr(nlon) linedata[*] = 0. ; ; Pick out data based on zonal type of line plot ; ; if line.latmin lt line.latmax then begin ; kk = 0 ; for k=0,nlon-1 do $ ; if lons[k] ge line.lonmin and lons[k] le line.lonmax then begin ; if line.meridionalmean eq 0 then begin ; linedata[kk] = float(data[k,iLat,iLev,line.imtime]) ; endif else begin ; plot meridional mean ; linedata[kk] = float(mean(data[k,*,iLev,line.imtime])) ; endelse ; kk = kk+1 ; endif ; endif else begin ; print, ' Loading data2 for time ', line.imtime ; ; Handle WACCM special case ; IF line.ftype EQ 'WACCM' THEN BEGIN ; ; No mean values needed for trace ; if line.meridionalmean eq 0 then begin linedata[*] = float(data[*,iLat,iFieldLev,line.imtime]) ; ; Need to handle fields on interfaces levels(instead of midpoints) ; IF iLev GT 0 THEN BEGIN IF iFieldLev GT 0 AND iFieldLev LT nLevs-1 THEN BEGIN IF isdensity(line.field.name) THEN BEGIN dataLower = ALOG10(data[*,iLat,IFieldLev-1,line.imtime]) dataLev = ALOG10(data[*,iLat,IFieldLev,line.imtime]) dataUpper = ALOG10(data[*,iLat,IFieldLev+1,line.imtime]) ENDIF ELSE BEGIN dataLower = data[*,iLat,IFieldLev-1,line.imtime] dataLev = data[*,iLat,IFieldLev,line.imtime] dataUpper = data[*,iLat,IFieldLev+1,line.imtime] ENDELSE IF newLevSelected LT levSelected THEN linedata = float(dataLower)*levRatUpper + $ float(dataLev)*levRatLower IF newLevSelected GT levSelected THEN linedata = float(dataLev)*levRatUpper + $ float(dataUpper)*levRatLower IF isdensity(line.field.name) THEN linedata = 10^(linedata) ENDIF ENDIF ; ; Handle wave amplitude calculation here ; IF line.amplitude NE 0 THEN BEGIN ampdata = fltarr(nlon) ampmin = 1.0E+24 ampdata = ABS(FFT(linedata, -1)) * 2. FOR imin = 0, nlon-1 DO BEGIN if ampdata[imin] LT ampmin AND ampdata[imin] GT 0.0 then ampmin = ampdata[imin] ENDFOR zeroind = WHERE(ampdata[*] EQ 0.0) if zeroind[0] NE -1 then ampdata[zeroind] = ampmin linedata = ampdata IF line.amplitude EQ 2 THEN linedata = ALOG10(ampdata) fmin = min(linedata) & fmax = max(linedata) line.fmin = fmin line.fmax = fmax line.lonmin = 0 line.lonmax = nlon ENDIF endif else begin ; meridional mean for k=0,nlon-1 do begin linedata[k] = float(mean(data[k,*,iFieldLev,line.imtime])) endfor ; ; Need to handle fields on interfaces levels(instead of midpoints) ; IF iLev GT 0 THEN BEGIN dataLower = linedata dataLev = linedata dataUpper = linedata IF iFieldLev GT 0 AND iFieldLev LT nLevs-1 THEN BEGIN IF isdensity(line.field.name) THEN BEGIN FOR k=0,nlon-1 do begin dataLower[k] = ALOG10(mean(data[k,*,IFieldLev-1,line.imtime])) dataLev[k] = ALOG10(mean(data[k,*,IFieldLev,line.imtime])) dataUpper[k] = ALOG10(mean(data[k,*,IFieldLev+1,line.imtime])) ENDFOR ENDIF ELSE BEGIN FOR k=0,nlon-1 do begin dataLower[k] = mean(data[k,*,IFieldLev-1,line.imtime]) dataLev[k] = mean(data[k,*,IFieldLev,line.imtime]) dataUpper[k] = mean(data[k,*,IFieldLev+1,line.imtime]) ENDFOR ENDELSE IF newLevSelected LT levSelected THEN linedata = float(dataLower)*levRatUpper + $ float(dataLev)*levRatLower IF newLevSelected GT levSelected THEN linedata = float(dataLev)*levRatUpper + $ float(dataUpper)*levRatLower IF isdensity(line.field.name) THEN linedata = 10^(linedata) ENDIF ENDIF endelse ; meridional mean ENDIF ELSE BEGIN ; Not WACCM if line.meridionalmean eq 0 then begin linedata[*] = float(data[*,iLat,iLev,line.imtime]) ; ; Handle wave amplitude calculation here ; IF line.amplitude NE 0 THEN BEGIN ampdata = fltarr(nlon) ampmin = 1.0E+24 ampdata = ABS(FFT(linedata, -1)) * 2. FOR imin = 0, nlon-1 DO BEGIN if ampdata[imin] LT ampmin AND ampdata[imin] GT 0.0 then ampmin = ampdata[imin] ENDFOR zeroind = WHERE(ampdata[*] EQ 0.0) if zeroind[0] NE -1 then ampdata[zeroind] = ampmin linedata = ampdata IF line.amplitude EQ 2 THEN linedata = ALOG10(ampdata) fmin = min(linedata) & fmax = max(linedata) line.fmin = fmin line.fmax = fmax line.lonmin = 0 line.lonmax = nlon ENDIF endif else begin ; meridional mean for k=0,nlon-1 do begin linedata[k] = float(mean(data[k,*,iLev,line.imtime])) endfor endelse ; meridional mean ENDELSE ;WACCM or not ; ; Take log10 of densities: ; if (line.log10 eq 'density fields' and isdensity(field.name) gt 0) or $ (line.log10 eq 'this frame only') or $ (line.log10 eq 'all fields') AND (line.amplitude EQ 0) then begin ; print,'defzonaldata: taking log10 of field ',field.name,'...' missing_value = float(line.missing_value) log10f,linedata,missing_value endif ; ; Update line structure: ; line.data = ptr_new(linedata) end