Index: /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 27788) +++ /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 27789) @@ -120,4 +120,6 @@ iomodel->FetchDataToInput(inputs,elements,"md.geometry.bed",BedEnum); break; + case CalvingCalvingMIPEnum: + break; /*"Discrete" calving laws (need to specify rate as 0 so that we can still solve the level set equation)*/ @@ -150,5 +152,5 @@ case FrontalForcingsDefaultEnum: iomodel->FetchDataToInput(inputs,elements,"md.frontalforcings.meltingrate",CalvingMeltingrateEnum); - if (calvinglaw == CalvingParameterizationEnum) { + if ((calvinglaw == CalvingParameterizationEnum) || (calvinglaw == CalvingCalvingMIPEnum)) { iomodel->FetchDataToInput(inputs,elements,"md.frontalforcings.ablationrate",CalvingAblationrateEnum); } @@ -256,4 +258,8 @@ case CalvingPollardEnum: parameters->AddObject(iomodel->CopyConstantObject("md.calving.rc",CalvingRcEnum)); + break; + case CalvingCalvingMIPEnum: + parameters->AddObject(iomodel->CopyConstantObject("md.calving.experiment",CalvingUseParamEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.calving.min_thickness",CalvingMinthicknessEnum)); break; default: @@ -748,5 +754,5 @@ /*Apply minimum thickness criterion*/ - if(calvinglaw==CalvingMinthicknessEnum || calvinglaw==CalvingVonmisesEnum || calvinglaw==CalvingParameterizationEnum || calvinglaw==CalvingVonmisesADEnum){ + if(calvinglaw==CalvingMinthicknessEnum || calvinglaw==CalvingVonmisesEnum || calvinglaw==CalvingParameterizationEnum || calvinglaw==CalvingVonmisesADEnum || calvinglaw==CalvingCalvingMIPEnum){ IssmDouble mig_max = femmodel->parameters->FindParam(MigrationMaxEnum); Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 27788) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 27789) @@ -3793,4 +3793,7 @@ case CalvingParameterizationEnum: this->CalvingRateParameterization(); + break; + case CalvingCalvingMIPEnum: + this->CalvingRateCalvingMIP(); break; case CalvingTestEnum: Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 27788) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 27789) @@ -238,4 +238,5 @@ virtual void BasalNodeIndices(int* pnumindices,int** pindices,int finiteelement){_error_("not implemented yet");}; virtual void CalvingRateParameterization(void){_error_("not implemented yet");}; + virtual void CalvingRateCalvingMIP(void){_error_("not implemented yet");}; virtual void CalvingRateVonmises(void){_error_("not implemented yet");}; virtual void CalvingRateVonmisesAD(void){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 27788) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 27789) @@ -1115,4 +1115,72 @@ this->AddInput(SigmaVMEnum,&sigma_vm[0],P1DGEnum); this->CalvingRateToVector(); + + /*Clean up and return*/ + delete gauss; +} +/*}}}*/ +void Tria::CalvingRateCalvingMIP(){/*{{{*/ + + IssmDouble calvingrate[NUMVERTICES]; + IssmDouble calvingratex[NUMVERTICES]; + IssmDouble calvingratey[NUMVERTICES]; + int experiment = 1; /* exp:1 by default */ + int dim, domaintype; + IssmDouble vx, vy, vel, c, wrate; + + /*Get problem dimension and whether there is moving front or not*/ + this->FindParam(&domaintype,DomainTypeEnum); + + switch(domaintype){ + case Domain2DverticalEnum: dim = 1; break; + case Domain2DhorizontalEnum: dim = 2; break; + case Domain3DEnum: dim = 2; break; + default: _error_("mesh "<GetInput(VxEnum); _assert_(vx_input); + Input *vy_input = this->GetInput(VyEnum); _assert_(vy_input); + Input *wrate_input = this->GetInput(CalvingAblationrateEnum); _assert_(wrate_input); + + /* Use which experiment: use existing Enum */ + this->FindParam(&experiment, CalvingUseParamEnum); + + /* Start looping on the number of vertices: */ + GaussTria* gauss=new GaussTria(); + for(int iv=0;ivGaussVertex(iv); + + /*Get velocity components */ + vx_input->GetInputValue(&vx,gauss); + vy_input->GetInputValue(&vy,gauss); + vel=sqrt(vx*vx+vy*vy)+1.e-14; + + + switch (experiment) { + case 1: + case 3: + /* Exp 1 and 3: set c=v-wrate, wrate=0, so that w=0 */ + wrate = 0.0; + break; + case 2: + case 4: + /* Exp 2 and 4: set c=v-wrate(given)*/ + /*Get wrate*/ + wrate_input->GetInputValue(&wrate,gauss); + break; + default: + _error_("The experiment is not supported yet!"); + } + + calvingrate[iv] = vel - wrate; + calvingratex[iv] = vx - wrate*vx/vel; + calvingratey[iv] = vy - wrate*vy/vel; + } + /*Add input*/ + this->AddInput(CalvingCalvingrateEnum,&calvingrate[0],P1DGEnum); + this->AddInput(CalvingratexEnum,&calvingratex[0],P1DGEnum); + this->AddInput(CalvingrateyEnum,&calvingratey[0],P1DGEnum); /*Clean up and return*/ @@ -4483,4 +4551,5 @@ case CalvingTestEnum: case CalvingParameterizationEnum: + case CalvingCalvingMIPEnum: calvingratex_input=this->GetInput(CalvingratexEnum); _assert_(calvingratex_input); calvingratey_input=this->GetInput(CalvingrateyEnum); _assert_(calvingratey_input); @@ -4528,4 +4597,5 @@ case CalvingLevermannEnum: case CalvingPollardEnum: + case CalvingCalvingMIPEnum: calvingratex_input->GetInputValue(&c[0],&gauss); calvingratey_input->GetInputValue(&c[1],&gauss); Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 27788) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 27789) @@ -63,4 +63,5 @@ void CalvingMeltingFluxLevelset(); void CalvingRateParameterization(); + void CalvingRateCalvingMIP(); IssmDouble CharacteristicLength(void); void ComputeBasalStress(void); Index: /issm/trunk-jpl/src/c/modules/Calvingx/Calvingx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/Calvingx/Calvingx.cpp (revision 27788) +++ /issm/trunk-jpl/src/c/modules/Calvingx/Calvingx.cpp (revision 27789) @@ -51,4 +51,7 @@ femmodel->ElementOperationx(&Element::CalvingPollard); break; + case CalvingCalvingMIPEnum: + femmodel->ElementOperationx(&Element::CalvingRateCalvingMIP); + break; default: _error_("Caving law "<=',0,'NaN',1,'Inf',1,'numel',1); + md = checkfield(md,'fieldname','calving.experiment','values',[1, 2, 3, 4, 5]); + end % }}} + function disp(self) % {{{ + disp(sprintf(' CalvingMIP parameters:')); + fielddisplay(self,'experiment','Experiment in CalvingMIP'); + fielddisplay(self,'min_thickness','minimum thickness below which no ice is allowed [m]'); + end % }}} + function marshall(self,prefix,md,fid) % {{{ + yts=md.constants.yts; + WriteData(fid,prefix,'name','md.calving.law','data',12,'format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','min_thickness','format','Double'); + WriteData(fid,prefix,'object',self,'fieldname','experiment','format','Integer'); + end % }}} + end +end Index: /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/expxy2shpll.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/expxy2shpll.m (revision 27789) +++ /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/expxy2shpll.m (revision 27789) @@ -0,0 +1,69 @@ +function expxy2shpll(expfilename,shpfilename,geometry) +%expxy2shpll +% Convert .exp to .shp file +% +% Usage: +% expxy2shpll('glacier.exp', 'glacier.shp') +% expxy2shpll('glacier.exp', 'glacier.shp', geometry) +% +% geometry (optional)- 'MultiPoint' : point clouds +% 'Point' : single point +% 'Line' : two points line +% 'polygon' : multiple points +% + +%check file extensions +[pathstr,name,ext] = fileparts(shpfilename); +if ~strcmp(ext,'.shp'), + error(['Shapefile ' shpfilename ' does not have an extension .shp']); +end + +[pathstr,name,ext] = fileparts(expfilename); +if ~strcmp(ext,'.exp'), + error(['Exp file ' expfilename ' does not have an extension .exp']); +end + +shp=expread(expfilename); +[lat, lon] = xy2ll(shp.x, shp.y, 1, 45, 70); +shp.x = lon; +shp.y = lat; + +%initialize number of profile +count=1; + +contours=struct([]); +for i=1:length(shp), + if nargin < 3 + + %TEMP + %if contains(shp(i).name,'_pointcloud'); + % continue; + %end + + if length(shp(i).x) == 0 + continue; + elseif contains(shp(i).name,'_pointcloud'); + geometry = 'MultiPoint'; + shp(i).name = erase(shp(i).name,'_pointcloud'); + elseif length(shp(i).x) == 1 + geometry = 'Point'; + elseif length(shp(i).x) < 3 + geometry = 'Line'; + else + if (shp(i).x(end)==shp(i).x(1) && shp(i).y(end)==shp(i).y(1)), + geometry = 'Polygon'; + else + geometry = 'Line'; + end + end + end + contours(count).Geometry=geometry; + contours(count).id=i; + contours(count).Name=shp(i).name; + contours(count).X=shp(i).x; + contours(count).Y=shp(i).y; + count = count+1; +end + +%Make sure it is one single geometry otherwise it will yell at you +shapewrite(contours,shpfilename); Index: sm/trunk-jpl/src/m/contrib/chenggong/expxy2shpll.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/chenggong/expxy2shpll.m (revision 27788) +++ (revision ) @@ -1,69 +1,0 @@ -function expxy2shpll(expfilename,shpfilename,geometry) -%expxy2shpll -% Convert .exp to .shp file -% -% Usage: -% expxy2shpll('glacier.exp', 'glacier.shp') -% expxy2shpll('glacier.exp', 'glacier.shp', geometry) -% -% geometry (optional)- 'MultiPoint' : point clouds -% 'Point' : single point -% 'Line' : two points line -% 'polygon' : multiple points -% - -%check file extensions -[pathstr,name,ext] = fileparts(shpfilename); -if ~strcmp(ext,'.shp'), - error(['Shapefile ' shpfilename ' does not have an extension .shp']); -end - -[pathstr,name,ext] = fileparts(expfilename); -if ~strcmp(ext,'.exp'), - error(['Exp file ' expfilename ' does not have an extension .exp']); -end - -shp=expread(expfilename); -[lat, lon] = xy2ll(shp.x, shp.y, 1, 45, 70); -shp.x = lon; -shp.y = lat; - -%initialize number of profile -count=1; - -contours=struct([]); -for i=1:length(shp), - if nargin < 3 - - %TEMP - %if contains(shp(i).name,'_pointcloud'); - % continue; - %end - - if length(shp(i).x) == 0 - continue; - elseif contains(shp(i).name,'_pointcloud'); - geometry = 'MultiPoint'; - shp(i).name = erase(shp(i).name,'_pointcloud'); - elseif length(shp(i).x) == 1 - geometry = 'Point'; - elseif length(shp(i).x) < 3 - geometry = 'Line'; - else - if (shp(i).x(end)==shp(i).x(1) && shp(i).y(end)==shp(i).y(1)), - geometry = 'Polygon'; - else - geometry = 'Line'; - end - end - end - contours(count).Geometry=geometry; - contours(count).id=i; - contours(count).Name=shp(i).name; - contours(count).X=shp(i).x; - contours(count).Y=shp(i).y; - count = count+1; -end - -%Make sure it is one single geometry otherwise it will yell at you -shapewrite(contours,shpfilename); Index: sm/trunk-jpl/src/m/contrib/chenggong/interpFromMEaSUREsGeotiff.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/chenggong/interpFromMEaSUREsGeotiff.m (revision 27788) +++ (revision ) @@ -1,68 +1,0 @@ -function dataout = interpFromMEaSUREsGeotiff(X,Y,Tstart,Tend,varargin) -%interpFromMEaSUREsGeotiff: -% This function calls src/m/contrib/morlighem/modeldata/interpFromGeotiff.m for multiple times to load all avaliable -% tif data in /totten_1/ModelData/Greenland/VelMEaSUREs/Jakobshavn_2008_2021/ within the given time period (in decimal years) -% For some reason, each .tif file in this folder contains two sets of data, only the first dataset is useful -% -% Usage: -% dataout = interpFromMEaSUREsGeotiff(X,Y,Tstart,Tend, varargin) -% -% X, Y are the coordinates of the mesh -% Tstart and Tend decimal year of the start and end time -% -% Example: -% obsData = interpFromMEaSUREsGeotiff(md.mesh.x,md.mesh.y, tstart, tend); -% -% Options: -% - 'glacier': which glacier to look for -options = pairoptions(varargin{:}); -glacier = getfieldvalue(options,'glacier','Jakobshavn'); - -if strcmp(glacier, 'Jakobshavn') - foldername = '/totten_1/ModelData/Greenland/VelMEaSUREs/Jakobshavn_2008_2021/'; -elseif strcmp(glacier, 'Kangerlussuaq') - foldername = '/totten_1/ModelData/Greenland/VelMEaSUREs/Kangerlussuaq_2006_2021/'; -elseif strcmp(glacier, 'Store') - foldername = '/totten_1/ModelData/Greenland/VelMEaSUREs/Store_2008_2021/'; -elseif strcmp(glacier, 'Rink') - foldername = '/totten_1/ModelData/Greenland/VelMEaSUREs/Rink_2008_2022/'; -elseif strcmp(glacier, 'Upernavik') - foldername = '/totten_1/ModelData/Greenland/VelMEaSUREs/Upernavik_2008_2022/'; -elseif strcmp(glacier, 'Helheim') - foldername = '/totten_1/ModelData/Greenland/VelMEaSUREs/Helheim_2008_2023/'; -else - error(['The velocity data for ', glacier, ' is not available, please download from NSIDC first.']); -end - -% get the time info from file names -templist = dir([foldername,'*.meta']); -Ndata = length(templist); -dataTstart = zeros(Ndata,1); -dataTend = zeros(Ndata,1); - -for i = 1:Ndata - tempConv = split(templist(i).name, '_'); - % follow the naming convention - dataPrefix(i) = join(tempConv(1:5), '_'); - dataTstart(i) = date2decyear(datenum(tempConv{3})); - dataTend(i) = date2decyear(datenum(tempConv{4})); -end -disp([' Found ', num2str(Ndata), ' records in ', foldername]); -disp([' from ', datestr(decyear2date(min(dataTstart)),'yyyy-mm-dd'), ' to ', datestr(decyear2date(max(dataTend)),'yyyy-mm-dd') ]); - - -% find all the data files with Tstart<=t<=Tend -dataInd = (dataTend>=Tstart) & (dataTstart<=Tend); -disp([' For the selected period: ', datestr(decyear2date((Tstart)),'yyyy-mm-dd'), ' to ', datestr(decyear2date((Tend)),'yyyy-mm-dd'), ', there are ', num2str(sum(dataInd)), ' records' ]); - -dataToLoad = dataPrefix(dataInd); -TstartToload = dataTstart(dataInd); -TendToload = dataTend(dataInd); - -for i = 1:length(dataToLoad) - dataout(i).vx = interpFromGeotiff([foldername, dataToLoad{i}, '_vx_v04.0.tif'], X, Y, 2e9); - dataout(i).vy = interpFromGeotiff([foldername, dataToLoad{i}, '_vy_v04.0.tif'], X, Y, 2e9); - dataout(i).vel = interpFromGeotiff([foldername, dataToLoad{i}, '_vv_v04.0.tif'], X, Y, -1); - dataout(i).Tstart = TstartToload(i); - dataout(i).Tend = TendToload(i); -end Index: /issm/trunk-jpl/src/m/contrib/chenggong/modeldata/interpFromAtlasDEM.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/chenggong/modeldata/interpFromAtlasDEM.m (revision 27789) +++ /issm/trunk-jpl/src/m/contrib/chenggong/modeldata/interpFromAtlasDEM.m (revision 27789) @@ -0,0 +1,106 @@ +function dataout = interpFromAtlasDEM(X,Y,Tstart,Tend,varargin) + %interpFromAtlasDEM: + % This function calls src/m/contrib/morlighem/modeldata/interpFromGeotiff.m for multiple times to load all avaliable + % + % Usage: + % dataout = interpFromAtlasDEM(X,Y,Tstart,Tend, varargin) + % + % X, Y are the coordinates of the mesh + % Tstart and Tend decimal year of the start and end time + % + % Example: + % obsData = interpFromAtlasDEM(md.mesh.x,md.mesh.y, tstart, tend); + % + % Options: + options = pairoptions(varargin{:}); + + foldername = '/totten_1/ModelData/Greenland/Helheim_ATLAS/'; + + % get the time info from file names + templist = dir([foldername,'*.tif']); + Ndata = length(templist); + dataTime = zeros(Ndata,1); + + for i = 1:Ndata + tempConv = split(templist(i).name, '-'); + % follow the naming convention + dataTime(i) = date2decyear(datenum(tempConv{1}, 'yymmdd_hhMMss')); + end + disp([' Found ', num2str(Ndata), ' records in ', foldername]); + disp([' from ', datestr(decyear2date(min(dataTime)),'yyyy-mm-dd'), ' to ', datestr(decyear2date(max(dataTime)),'yyyy-mm-dd') ]); + + + % find all the data files with Tstart<=t<=Tend + dataInd = (dataTime>=Tstart) & (dataTime<=Tend); + disp([' For the selected period: ', datestr(decyear2date((Tstart)),'yyyy-mm-dd'), ' to ', datestr(decyear2date((Tend)),'yyyy-mm-dd'), ', there are ', num2str(sum(dataInd)), ' records' ]); + + dataToLoad = {templist(dataInd).name}; + timeToload = dataTime(dataInd); + + for i = 1:length(dataToLoad) + tifdata= interpFromTif([foldername, dataToLoad{i}], X, Y, 2e9); + dataout(i).name = dataToLoad{i}; + dataout(i).surface = tifdata(:,:,3); + dataout(i).Time = timeToload(i); + end + +end + + function dataout = interpFromTif(tifname,X,Y,nanValue) % {{{ + + if nargin < 4 + nanValue = 10^30; + end + + usemap = 0; + + %Get image info + Tinfo = imfinfo(tifname); + N = Tinfo(1).Width; + M = Tinfo(1).Height; + dx = Tinfo(1).ModelPixelScaleTag(1); + dy = Tinfo(1).ModelPixelScaleTag(2); + minx = Tinfo(1).ModelTiepointTag(4); + maxy = Tinfo(1).ModelTiepointTag(5); + + %Generate vectors + xdata = minx + dx/2 + ((0:N-1).*dx); + ydata = maxy - dy/2 - ((M -1:-1:0).*dy); + + %Read image + assert(dx>0); assert(dy>0); + ydata = fliplr(ydata); + + %Get pixels we are interested in + offset=2; + xmin=min(X(:)); xmax=max(X(:)); + posx=find(xdata<=xmax); + id1x=max(1,find(xdata>=xmin,1)-offset); + id2x=min(numel(xdata),posx(end)+offset); + + ymin=min(Y(:)); ymax=max(Y(:)); + posy=find(ydata>=ymin); + id1y=max(1,find(ydata<=ymax,1)-offset); + id2y=min(numel(ydata),posy(end)+offset); + + data = double(imread(tifname,'PixelRegion',{[id1y,id2y],[id1x,id2x]})); + xdata=xdata(id1x:id2x); + ydata=ydata(id1y:id2y); + + if nanValue > 0 + data(find(abs(data)>=nanValue))=NaN; + else + data(find(data<=nanValue))=NaN; + end + + if ndims(data) == 2 + dataout = InterpFromGrid(xdata,ydata,data,X,Y); + elseif ndims(data) == 3 + for i = 1:size(data, 3) + dataout(:,:,i) = InterpFromGrid(xdata, ydata, data(:,:, i), X, Y); + end + else + error(['not implemented for data of ', num2str(ndims(data)), 'dimensions!']) + end + dataout(dataout==-9999)=NaN; + end %}}} Index: /issm/trunk-jpl/src/m/contrib/chenggong/modeldata/interpFromMEaSUREsGeotiff.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/chenggong/modeldata/interpFromMEaSUREsGeotiff.m (revision 27789) +++ /issm/trunk-jpl/src/m/contrib/chenggong/modeldata/interpFromMEaSUREsGeotiff.m (revision 27789) @@ -0,0 +1,68 @@ +function dataout = interpFromMEaSUREsGeotiff(X,Y,Tstart,Tend,varargin) +%interpFromMEaSUREsGeotiff: +% This function calls src/m/contrib/morlighem/modeldata/interpFromGeotiff.m for multiple times to load all avaliable +% tif data in /totten_1/ModelData/Greenland/VelMEaSUREs/Jakobshavn_2008_2021/ within the given time period (in decimal years) +% For some reason, each .tif file in this folder contains two sets of data, only the first dataset is useful +% +% Usage: +% dataout = interpFromMEaSUREsGeotiff(X,Y,Tstart,Tend, varargin) +% +% X, Y are the coordinates of the mesh +% Tstart and Tend decimal year of the start and end time +% +% Example: +% obsData = interpFromMEaSUREsGeotiff(md.mesh.x,md.mesh.y, tstart, tend); +% +% Options: +% - 'glacier': which glacier to look for +options = pairoptions(varargin{:}); +glacier = getfieldvalue(options,'glacier','Jakobshavn'); + +if strcmp(glacier, 'Jakobshavn') + foldername = '/totten_1/ModelData/Greenland/VelMEaSUREs/Jakobshavn_2008_2021/'; +elseif strcmp(glacier, 'Kangerlussuaq') + foldername = '/totten_1/ModelData/Greenland/VelMEaSUREs/Kangerlussuaq_2006_2021/'; +elseif strcmp(glacier, 'Store') + foldername = '/totten_1/ModelData/Greenland/VelMEaSUREs/Store_2008_2021/'; +elseif strcmp(glacier, 'Rink') + foldername = '/totten_1/ModelData/Greenland/VelMEaSUREs/Rink_2008_2022/'; +elseif strcmp(glacier, 'Upernavik') + foldername = '/totten_1/ModelData/Greenland/VelMEaSUREs/Upernavik_2008_2022/'; +elseif strcmp(glacier, 'Helheim') + foldername = '/totten_1/ModelData/Greenland/VelMEaSUREs/Helheim_2008_2023/'; +else + error(['The velocity data for ', glacier, ' is not available, please download from NSIDC first.']); +end + +% get the time info from file names +templist = dir([foldername,'*.meta']); +Ndata = length(templist); +dataTstart = zeros(Ndata,1); +dataTend = zeros(Ndata,1); + +for i = 1:Ndata + tempConv = split(templist(i).name, '_'); + % follow the naming convention + dataPrefix(i) = join(tempConv(1:5), '_'); + dataTstart(i) = date2decyear(datenum(tempConv{3})); + dataTend(i) = date2decyear(datenum(tempConv{4})); +end +disp([' Found ', num2str(Ndata), ' records in ', foldername]); +disp([' from ', datestr(decyear2date(min(dataTstart)),'yyyy-mm-dd'), ' to ', datestr(decyear2date(max(dataTend)),'yyyy-mm-dd') ]); + + +% find all the data files with Tstart<=t<=Tend +dataInd = (dataTend>=Tstart) & (dataTstart<=Tend); +disp([' For the selected period: ', datestr(decyear2date((Tstart)),'yyyy-mm-dd'), ' to ', datestr(decyear2date((Tend)),'yyyy-mm-dd'), ', there are ', num2str(sum(dataInd)), ' records' ]); + +dataToLoad = dataPrefix(dataInd); +TstartToload = dataTstart(dataInd); +TendToload = dataTend(dataInd); + +for i = 1:length(dataToLoad) + dataout(i).vx = interpFromGeotiff([foldername, dataToLoad{i}, '_vx_v04.0.tif'], X, Y, 2e9); + dataout(i).vy = interpFromGeotiff([foldername, dataToLoad{i}, '_vy_v04.0.tif'], X, Y, 2e9); + dataout(i).vel = interpFromGeotiff([foldername, dataToLoad{i}, '_vv_v04.0.tif'], X, Y, -1); + dataout(i).Tstart = TstartToload(i); + dataout(i).Tend = TendToload(i); +end Index: /issm/trunk-jpl/src/m/contrib/chenggong/visualization/plotCompareTransientFlowline.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/chenggong/visualization/plotCompareTransientFlowline.m (revision 27788) +++ /issm/trunk-jpl/src/m/contrib/chenggong/visualization/plotCompareTransientFlowline.m (revision 27789) @@ -16,5 +16,5 @@ % % Author: Cheng Gong -% Date: 2023-04-10 +% Date: 2021-12-06 N = length(velList); @@ -46,7 +46,3 @@ h=colorbar; title(h,'m/a') - % add a vertical line to indicate initial ice front position - id = max(find(~isnan(cumsum(vel_flowline{p}(:,1))))); - icefront = flowline.Xmain(id); - plot([icefront, icefront], [2007,2020], '-r') end Index: /issm/trunk-jpl/src/m/exp/isoline.m =================================================================== --- /issm/trunk-jpl/src/m/exp/isoline.m (revision 27788) +++ /issm/trunk-jpl/src/m/exp/isoline.m (revision 27789) @@ -258,4 +258,7 @@ elseif strcmp(outputformat,'struct') %nothing to do, this is the default +elseif strcmp(outputformat,'longest') + [~, mId] = max([contours.nods]); + contours = contours(mId); else disp('output format not supported, returning struct'); Index: /issm/trunk-jpl/src/m/modeldata/interpBedmachineGreenland.m =================================================================== --- /issm/trunk-jpl/src/m/modeldata/interpBedmachineGreenland.m (revision 27788) +++ /issm/trunk-jpl/src/m/modeldata/interpBedmachineGreenland.m (revision 27789) @@ -32,5 +32,8 @@ ncdate='2021-08-27'; ncdate='2022-03-17'; + ncdate='2022-05-18'; ncdate='2022-07-28'; + ncdate='v6.0'; + ncdate='v6.1'; end basename = 'BedMachineGreenland';