Index: /issm/trunk/README =================================================================== --- /issm/trunk/README (revision 22821) +++ /issm/trunk/README (revision 22822) @@ -16,5 +16,5 @@ -Copyright (c) 2002-2014, California Institute of Technology. +Copyright (c) 2002-2018, California Institute of Technology. All rights reserved. Based on Government Sponsored Research under contracts NAS7-1407 and/or NAS7-03001. Index: /issm/trunk/configure.ac =================================================================== --- /issm/trunk/configure.ac (revision 22821) +++ /issm/trunk/configure.ac (revision 22822) @@ -2,5 +2,5 @@ #AUTOCONF -AC_INIT([Ice Sheet System Model (ISSM)],[4.13],[issm@jpl.nasa.gov],[issm],[http://issm.jpl.nasa.gov]) #Initializing configure +AC_INIT([Ice Sheet System Model (ISSM)],[4.14],[issm@jpl.nasa.gov],[issm],[http://issm.jpl.nasa.gov]) #Initializing configure AC_CONFIG_AUX_DIR([./aux-config]) #Put config files in aux-config AC_CONFIG_MACRO_DIR([m4]) #m4 macros are located in m4 Index: /issm/trunk/examples/EsaGRACE/runme.m =================================================================== --- /issm/trunk/examples/EsaGRACE/runme.m (revision 22822) +++ /issm/trunk/examples/EsaGRACE/runme.m (revision 22822) @@ -0,0 +1,150 @@ + +clear all; +addpath('../Data','../Functions'); + +steps=[1]; % [1:5] + +if any(steps==1) + disp(' Step 1: Global mesh creation'); + + resolution=300; % [km] + radius = 6.371012*10^3; % [km] + + md=model; + md.mask=maskpsl(); + md.mesh=gmshplanet('radius',radius,'resolution',resolution); + + md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); + + save ./Models/EsaGRACE_Mesh md; + + plotmodel (md,'data',md.mask.ocean_levelset,'edgecolor','k'); + +end + +if any(steps==2) + disp(' Step 2: Define loads in meters of ice height equivalent'); + md = loadmodel('./Models/EsaGRACE_Mesh'); + + year_month = 2007+15/365; + time_range = [year_month year_month]; + + water_load = grace(md.mesh.elements,md.mesh.lat,md.mesh.long,time_range(1),time_range(2)); + + md.esa.deltathickness = water_load*md.materials.rho_freshwater/md.materials.rho_ice; % ice height equivalent + + save ./Models/EsaGRACE_Loads md; + + plotmodel (md,'data',md.esa.deltathickness,... + 'view',[90 -90],'caxis',[-.1 .1],... + 'title','Ice height equivalent [m]'); + +end + +if any(steps==3) + disp(' Step 3: Parameterization'); + md = loadmodel('./Models/EsaGRACE_Loads'); + + nlove=10001; + md.esa.love_h = love_numbers('h','CF'); md.esa.love_h(nlove+1:end)=[]; + md.esa.love_l = love_numbers('l','CF'); md.esa.love_l(nlove+1:end)=[]; + + md.mask.groundedice_levelset = ones(md.mesh.numberofvertices,1); + md.mask.ice_levelset = ones(md.mesh.numberofvertices,1); + pos=find(md.esa.deltathickness~=0); + md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; + md.mask.land_levelset = 1-md.mask.ocean_levelset; + + di=md.materials.rho_ice/md.materials.rho_water; + md.geometry.thickness=ones(md.mesh.numberofvertices,1); + md.geometry.surface=(1-di)*zeros(md.mesh.numberofvertices,1); + md.geometry.base=md.geometry.surface-md.geometry.thickness; + md.geometry.bed=md.geometry.base; + + md.initialization.temperature=273.25*ones(md.mesh.numberofvertices,1); + md.materials.rheology_B=paterson(md.initialization.temperature); + md.materials.rheology_n=3*ones(md.mesh.numberofelements,1); + + md.miscellaneous.name='EsaGRACE'; + + save ./Models/EsaGRACE_Parameterization md; + +end + +if any(steps==4) + disp(' Step 4: Solve Esa solver'); + md = loadmodel('./Models/EsaGRACE_Parameterization'); + + md.esa.requested_outputs = {'EsaUmotion','EsaNmotion','EsaEmotion'}; + + md.cluster=generic('name',oshostname(),'np',3); + md.verbose=verbose('111111111'); + + md=solve(md,'Esa'); + + save ./Models/EsaGRACE_Solution md; + +end + +if any(steps==5) + disp(' Step 5: Plot solutions'); + md = loadmodel('./Models/EsaGRACE_Solution'); + + sol1 = md.esa.deltathickness*100; % [cm] + sol2 = md.results.EsaSolution.EsaUmotion*1000; % [mm] + sol3 = md.results.EsaSolution.EsaNmotion*1000; % [mm] + sol4 = md.results.EsaSolution.EsaEmotion*1000; % [mm] + + sol_name={'Change in water equivalent height [cm]', 'Vertical displacement [mm]',... + 'Horizontal (NS) displacement [mm]', 'Horizontal (EW) displacement [mm]'}; + fig_name={'Fig_dH.pdf','Fig_vert.pdf','Fig_horzNS.pdf','Fig_horzEW.pdf'}; + + res = 1.0; % [degree] + + [lat_grid, lon_grid] = meshgrid(linspace(-90,90,180/res), linspace(-180,180,360/res)); + sol_grid = zeros(size(lat_grid)); + + for kk=1:4 + sol=eval(sprintf('sol%d',kk)); + + if length(sol)==md.mesh.numberofelements + for jj=1:md.mesh.numberofelements + ii=(jj-1)*3; + pp(ii+1:ii+3)=md.mesh.elements(jj,:); + end + for jj=1:md.mesh.numberofvertices + pos=ceil(find(pp==jj)/3); + temp(jj)=mean(sol(pos)); + end + sol=temp'; + end + + F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol); + F.Method = 'linear'; + F.ExtrapolationMethod = 'linear'; + + sol_grid = F(lat_grid, lon_grid); + sol_grid(isnan(sol_grid))=0; + sol_grid(lat_grid>85 & sol_grid==0)=NaN; + + set(0,'DefaultAxesFontSize',18,'DefaultAxesLineWidth',1,'DefaultTextFontSize',18,'DefaultLineMarkerSize',8) + figure1=figure('Position', [100, 100, 1000, 500]); + gcf; load coast; cla; + pcolor(lon_grid,lat_grid,sol_grid); shading flat; hold on; + if (kk==1) + geoshow(flipud(lat),flipud(long),'DisplayType','polygon','FaceColor','white'); + end + plot(long,lat,'k'); hold off; + c1=colorbar; + colormap('haxby'); + caxis([-min(abs(min(sol)),abs(max(sol))) min(abs(min(sol)),abs(max(sol)))]); + xlim([-180 180]); + ylim([-90 90]); + grid on; + title(sol_name(kk)); + set(gcf,'color','w'); + %export_fig(fig_name{kk}); + end + +end + Index: /issm/trunk/examples/EsaWahr/RoundDomain.exp =================================================================== --- /issm/trunk/examples/EsaWahr/RoundDomain.exp (revision 22822) +++ /issm/trunk/examples/EsaWahr/RoundDomain.exp (revision 22822) @@ -0,0 +1,49 @@ +## Name: +## Icon:0 +# Points Count Value +41 1 +# X pos Y pos +20000 0 +19754 3128.7 +19021 6180.3 +17820 9079.8 +16180 11756 +14142 14142 +11756 16180 +9079.8 17820 +6180.3 19021 +3128.7 19754 +0 20000 +-3128.7 19754 +-6180.3 19021 +-9079.8 17820 +-11756 16180 +-14142 14142 +-16180 11756 +-17820 9079.8 +-19021 6180.3 +-19754 3128.7 +-20000 0 +-19754 -3128.7 +-19021 -6180.3 +-17820 -9079.8 +-16180 -11756 +-14142 -14142 +-11756 -16180 +-9079.8 -17820 +-6180.3 -19021 +-3128.7 -19754 +-0 -20000 +3128.7 -19754 +6180.3 -19021 +9079.8 -17820 +11756 -16180 +14142 -14142 +16180 -11756 +17820 -9079.8 +19021 -6180.3 +19754 -3128.7 +20000 0 + + + Index: /issm/trunk/examples/EsaWahr/runme.m =================================================================== --- /issm/trunk/examples/EsaWahr/runme.m (revision 22822) +++ /issm/trunk/examples/EsaWahr/runme.m (revision 22822) @@ -0,0 +1,173 @@ + +clear all; +addpath('../Functions'); + +steps=[0]; % [0:6] + +if any(steps==0) + disp(' Step 0: Mesh creation'); + + md=roundmesh(model,100000,10000); % Domain radius and element size [m] + + save ./Models/EsaWahr_Mesh md; + + plotmodel(md,'data','mesh'); + +end + +if any(steps==1) + disp(' Step 1: Anisotropic mesh creation'); + + md=roundmesh(model,100000,1000); + + disc_radius=20*1000; + rad_dist=sqrt(md.mesh.x.^2+md.mesh.y.^2); + field = abs(rad_dist-disc_radius); + + md = bamg(md,'field',field,'err',50,'hmax',10000); + + save ./Models/EsaWahr_Mesh md; + + plotmodel (md,'data','mesh'); + +end + +if any(steps==2) + disp(' Step 2: Define loads'); + md = loadmodel('./Models/EsaWahr_Mesh'); + + rho_w_i=md.materials.rho_freshwater/md.materials.rho_ice; + + index=md.mesh.elements; + x_cent=mean(md.mesh.x(index),2); + y_cent=mean(md.mesh.y(index),2); + + md.esa.deltathickness = zeros(md.mesh.numberofelements,1); + disc_radius=20; % [km] + rad_dist=sqrt(x_cent.^2+y_cent.^2)/1000; + md.esa.deltathickness(rad_dist<=disc_radius) = -1.0*rho_w_i; + + save ./Models/EsaWahr_Loads md; + + plotmodel (md,'data',md.esa.deltathickness,'title','Ice height equivalent [m]'); + +end + +if any(steps==3) + disp(' Step 3: Parameterization'); + md = loadmodel('./Models/EsaWahr_Loads'); + + nlove=10001; + md.esa.love_h = love_numbers('h','CF'); md.esa.love_h(nlove+1:end)=[]; + md.esa.love_l = love_numbers('l','CF'); md.esa.love_l(nlove+1:end)=[]; + + md.mask.ice_levelset = -ones(md.mesh.numberofvertices,1); + md.mask.groundedice_levelset = ones(md.mesh.numberofvertices,1); + + di=md.materials.rho_ice/md.materials.rho_water; + md.geometry.thickness=ones(md.mesh.numberofvertices,1); + md.geometry.surface=(1-di)*zeros(md.mesh.numberofvertices,1); + md.geometry.base=md.geometry.surface-md.geometry.thickness; + md.geometry.bed=md.geometry.base; + + md.initialization.temperature=273.25*ones(md.mesh.numberofvertices,1); + md.materials.rheology_B=paterson(md.initialization.temperature); + md.materials.rheology_n=3*ones(md.mesh.numberofelements,1); + + md.miscellaneous.name='EsaWahr'; + + save ./Models/EsaWahr_Parameterization md; + +end + +if any(steps==4) + disp(' Step 4: Solve Esa solver'); + md = loadmodel('./Models/EsaWahr_Parameterization'); + + md.esa.requested_outputs = {'EsaUmotion','EsaXmotion','EsaYmotion'}; + + md.cluster=generic('name',oshostname(),'np',3); + md.verbose=verbose('111111111'); + + md=solve(md,'Esa'); + + save ./Models/EsaWahr_Solution md; + +end + +if any(steps==5) + disp(' Step 5: Plot solutions'); + md = loadmodel('./Models/EsaWahr_Solution'); + + vert = md.results.EsaSolution.EsaUmotion*1000; % [mm] + horz_n = md.results.EsaSolution.EsaYmotion*1000; % [mm] + horz_e = md.results.EsaSolution.EsaXmotion*1000; % [mm] + horz = sqrt(horz_n.^2+horz_e.^2); % [mm] + + set(0,'DefaultAxesFontSize',24,'DefaultAxesLineWidth',1,'DefaultTextFontSize',24,'DefaultLineMarkerSize',6); + figure('Position', [100, 100, 800, 600]); + plotmodel(md,'data',vert,... + 'xTickLabel#all',[],'yTickLabel#all',[],... + 'colormap#all','jet','colorbarcornerposition#all','south',... + 'expdisp#all','./RoundDomain.exp',... + 'gridded#all',1,... + 'axispos#1',[0.02 0.505 0.47 0.47],... + 'colorbarpos#1',[0.045,0.55,0.18,0.02],'colorbartitle#1','Vertical [mm]',... + 'caxis#1',[0 3.5],... + 'data',horz,... + 'axispos#2',[0.505 0.505 0.47 0.47],... + 'colorbarpos',[0.53,0.55,0.18,0.02],'colorbartitle#2','Horizontal [mm]',... + 'caxis#2',[0 0.5],... + 'data',horz_n,... + 'axispos',[0.02 0.02 0.47 0.47],... + 'colorbarpos',[0.045,0.065,0.18,0.02],'colorbartitle#3','North-south [mm]',... + 'data',horz_e,... + 'caxis#3-4',[-0.5 0.5],... + 'axispos',[0.505 0.02 0.47 0.47],... + 'colorbarpos',[0.53,0.065,0.18,0.02],'colorbartitle#4','East-west [mm]'); + %export_fig('Fig5.pdf'); + +end + +if any(steps==6) + disp(' Step 6: Compare results against Wahr semi-analytic solutions'); + md = loadmodel('./Models/EsaWahr_Solution'); + + vert = md.results.EsaSolution.EsaUmotion*1000; % [mm] + horz_n = md.results.EsaSolution.EsaYmotion*1000; % [mm] + horz_e = md.results.EsaSolution.EsaXmotion*1000; % [mm] + horz = sqrt(horz_n.^2+horz_e.^2); % [mm] + + xi=[0:500:100000]; % grid points [m] + yi=zeros(1,length(xi)); + vert_track=griddata(md.mesh.x,md.mesh.y,vert,xi,yi,'linear'); + horz_track=griddata(md.mesh.x,md.mesh.y,horz,xi,yi,'linear'); + + % semi-analytic solution (Wahr et al., 2013, JGR, Figure 1) + disc_radius = 20*1000; % [m] + [vert_wahr, horz_wahr] = wahr(disc_radius,xi,md.esa.love_h,md.esa.love_l); + + set(0,'DefaultAxesFontSize',16,'DefaultAxesLineWidth',1,'DefaultTextFontSize',16,'DefaultLineMarkerSize',6); + figure1=figure('Position', [100, 100, 700, 400]); + ylabel_1={'0',' ','1','','2','','3',''}; + axes1 = axes('Layer','top','Position',[0.1 0.15 0.8 0.8],... + 'XTick',[0:10:100],'xlim',[0 100],... + 'ylim',[0 3.5],'Ytick',[0:0.5:3.5],'yticklabel',ylabel_1); + box(axes1,'on'); hold(axes1,'all'); grid on; + xlabel(axes1,'Radial distance [km]'); + ylabel(axes1,'Displacement [mm]'); + plot([20 20],[0 3.5],'-k','linewidth',2,'parent',axes1); + % analytic soln + h3=plot(xi/1000,vert_wahr,'-r','linewidth',5,'parent',axes1); + h4=plot(xi/1000,horz_wahr,'-m','linewidth',5,'parent',axes1); + % ISSM soln + h1=plot(xi/1000,vert_track*917/1000,'-b','linewidth',3,'parent',axes1); + h2=plot(xi/1000,horz_track*917/1000,'-c','linewidth',3,'parent',axes1); + ag1 = gca; + leg1a = legend(ag1,[h3,h1,h4,h2],'Vertical (Wahr)','Vertical (ISSM)','Horizontal (Wahr)','Horizontal (ISSM)'); + set(leg1a,'location','east','Orientation','Vertical','Box','Off','FontSize',16); + set(gcf,'color','w'); + %export_fig('Fig6.pdf'); + +end + Index: /issm/trunk/examples/Functions/grace.m =================================================================== --- /issm/trunk/examples/Functions/grace.m (revision 22822) +++ /issm/trunk/examples/Functions/grace.m (revision 22822) @@ -0,0 +1,121 @@ +% map grace loads in meters [m] of water equivalent height +function water_load = grace(index,lat,lon,tmin,tmax); + + %compute centroids using (lat,lon) data {{{ + ne = length(index); % number of elements + % lat -> [0,180]; long -> [0,360] to compute centroids + lat=90-lat; + lon(lon<0)=180+(180+lon(lon<0)); + + ax_0=lat(index(:,1)); ay_0=lon(index(:,1)); + bx_0=lat(index(:,2)); by_0=lon(index(:,2)); + cx_0=lat(index(:,3)); cy_0=lon(index(:,3)); + % find whether long is 0 or 360! This is important to compute centroids as well as elemental area + for ii=1:ne + if (min([ay_0(ii),by_0(ii),cy_0(ii)])==0 && max([ay_0(ii),by_0(ii),cy_0(ii)])>180) + if ay_0(ii)==0 + ay_0(ii)=360; + end + if by_0(ii)==0 + by_0(ii)=360; + end + if cy_0(ii)==0 + cy_0(ii)=360; + end + end + end + % correction at the north pole + ay_0(ax_0==0)=(by_0(ax_0==0)+cy_0(ax_0==0))./2; + by_0(bx_0==0)=(cy_0(bx_0==0)+ay_0(bx_0==0))./2; + cy_0(cx_0==0)=(ay_0(cx_0==0)+by_0(cx_0==0))./2; + % correction at the south pole + ay_0(ax_0==180)=(by_0(ax_0==180)+cy_0(ax_0==180))./2; + by_0(bx_0==180)=(cy_0(bx_0==180)+ay_0(bx_0==180))./2; + cy_0(cx_0==180)=(ay_0(cx_0==180)+by_0(cx_0==180))./2; + % + lat_element=(ax_0+bx_0+cx_0)/3; + lon_element=(ay_0+by_0+cy_0)/3; + + % [lat,lon] \in [-90:90,-180,180]; + lat_element_0 = 90-lat_element; lon_element_0 = lon_element; + lon_element_0(lon_element>180) = (lon_element(lon_element>180)-180) - 180; + % }}} + + % Monthly GRACE data + filename=['GRCTellus.JPL.200204_201701.LND.RL05_1.DSTvSCS1411.nc']; + time_0=ncread(filename,'time'); % days since 2002-01-01 00:00:00 UTC + long_0=ncread(filename,'lon'); % longitudes: 0.27-359.75 + lati_0=ncread(filename,'lat'); % latitudes: -89.75:89.75 + rec=ncread(filename,'lwe_thickness');% rec_ensemble_mean [cm] + + time_yr = 2002+time_0/365; % [yr] + + [nn_0,mm_0] = size(squeeze(rec(:,:,1))); + for jj=1:mm_0 % chose a latitude + for kk=1:nn_0 + ii=nn_0*(jj-1)+kk; + lat_0(ii)=lati_0(jj); + lon_0(ii)=long_0(kk); + tws_monthly(:,ii) = rec(kk,jj,:); + end + end + % + %% translate variables as grace-related variables -- so I do not need to do too much editing + grace_monthly=tws_monthly; + grace_monthly(grace_monthly<-1344.6016 | grace_monthly>858.5046)=0; + + % detrend over the entire time period + grace_monthly = detrend(grace_monthly); % 159X64800 => remove trends from each column! + + % fill out the blanks {{{ + + lat_grace=lat_0; + lon_grace=lon_0; + num_org=length(lon_grace); + + qq=1; mm=1; + for jj=2:num_org-1 + if (lat_grace(jj)~=lat_grace(jj+1)) + lat_new(qq)=lat_grace(jj); + lon_new(qq)=lon_grace(jj)+(lon_grace(jj)-lon_grace(jj-1)); + load_new(:,qq)=grace_monthly(:,mm); + lat_new(qq+1)=lat_grace(jj); + lon_new(qq+1)=lon_grace(mm)-(lon_grace(jj)-lon_grace(jj-1)); + load_new(:,qq+1)=grace_monthly(:,jj); + qq=qq+2; + mm=jj+1; % to find out the value for monthly data + end + end + + num_add=length(lat_new); + num_plus=num_org+num_add; + + lat_grace_plus=zeros(num_plus,1); + lon_grace_plus=zeros(num_plus,1); + load_grace_plus=zeros(length(time_0),num_plus); + + lat_grace_plus(1:num_org)=lat_grace; + lat_grace_plus(1+num_org:num_plus)=lat_new; + lon_grace_plus(1:num_org)=lon_grace; + lon_grace_plus(1+num_org:num_plus)=lon_new; + load_grace_plus(:,1:num_org)=grace_monthly; + load_grace_plus(:,1+num_org:num_plus)=load_new; %(:,:); + % }}} + + % choose selected months ONLY + [diff1,pos1] = min(abs(tmin-time_yr)); + [diff2,pos2] = min(abs(tmax-time_yr)); + + time_yr=time_yr(pos1:pos2); + load_grace_plus=load_grace_plus(pos1:pos2,:); + num_yr=length(time_yr); + water_load_0=zeros(ne,num_yr); + + for jj=1:num_yr + water_load_0(:,jj) = griddata(lat_grace_plus,lon_grace_plus,load_grace_plus(jj,:),lat_element_0,lon_element); + disp([num2str(jj),' of ',num2str(num_yr),' months done!']); + end + + water_load=water_load_0/100; % cm -> meters of water + water_load(isnan(water_load))=0; + Index: /issm/trunk/examples/Functions/wahr.m =================================================================== --- /issm/trunk/examples/Functions/wahr.m (revision 22822) +++ /issm/trunk/examples/Functions/wahr.m (revision 22822) @@ -0,0 +1,46 @@ +% compute semi-analytic solutions for a disc load +function [vert, horz] = wahr(disc_rad,xi,love_h,love_l); + + disc_rad = disc_rad/1000; % km + % compute P(x), dP(x)/dx, d2P(x)/dx2 + %--------------------------------------------------------------------- + % compute p_value + theta=km2deg(xi/1000)'; + ang = theta/180*pi; + alpha=cos(ang); + m=length(alpha); + n=length(love_h)-1; + p_value = p_polynomial_value(m,n,alpha); + p_prime = p_polynomial_prime(m,n,alpha); + %--------------------------------------------------------------------- + nn=[0:n]; + nn_plus_1=nn+1; + + % disc radius in degree + disc_rad = km2deg(disc_rad)/180*pi; + tau=zeros(size(love_h)); + tau(1) = 0.5*(1-cos(disc_rad)); % tau_0 + p_value_disc = p_polynomial_value(1,n+1,cos(disc_rad)); + p_prime_disc = p_polynomial_prime(1,n,cos(disc_rad)); + for jj=2:n+1 + nnn = jj-1; + tau(jj) = 0.5 * (p_value_disc(jj-1) - p_value_disc(jj+1)); + end + + const=zeros(size(love_h)); + for jj=1:n+1 + const(jj) = 1/(2*(jj-1)+1); + end + + disc=sum(bsxfun(@times,p_value,tau'),2); + + g1 = -sum(bsxfun(@times,p_value,(tau.*love_h.*const)'),2); + g5 = -sum(bsxfun(@times,-sin(ang),bsxfun(@times,p_prime,(tau.*love_l.*const)')),2); + + % coeff + coeff = 1000*4*pi*(6.67408*10^-11)*(6.3781*10^6)/9.81; + + % vertical and horizontal solutions in mm + vert = g1*coeff*1000; % mm + horz = g5*coeff*1000; % mm + Index: /issm/trunk/examples/SlrFarrell/runme.m =================================================================== --- /issm/trunk/examples/SlrFarrell/runme.m (revision 22822) +++ /issm/trunk/examples/SlrFarrell/runme.m (revision 22822) @@ -0,0 +1,160 @@ + +clear all; + +steps=[1]; % [1:5] + +if any(steps==1) + disp(' Step 1: Global mesh creation'); + + numrefine=1; + resolution=150*1e3; % inital resolution [m] + radius = 6.371012*10^6; % mean radius of Earth, m + mindistance_coast=150*1e3; % coastal resolution [m] + mindistance_land=300*1e3; % resolution on the continents [m] + maxdistance=600*1e3; % max element size (on mid-oceans) [m] + + md=model; + md.mask=maskpsl(); + md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3); % attributes in [km] + + for i=1:numrefine, + + md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); + + distance=zeros(md.mesh.numberofvertices,1); + + pos=find(~md.mask.ocean_levelset); coaste.lat=md.mesh.lat(pos); coaste.long=md.mesh.long(pos); + pos=find(md.mask.ocean_levelset); coasto.lat=md.mesh.lat(pos); coasto.long=md.mesh.long(pos); + + for j=1:md.mesh.numberofvertices + phi1=md.mesh.lat(j)/180*pi; lambda1=md.mesh.long(j)/180*pi; + if md.mask.ocean_levelset(j), + phi2=coaste.lat/180*pi; lambda2=coaste.long/180*pi; + deltaphi=abs(phi2-phi1); deltalambda=abs(lambda2-lambda1); + d=radius*2*asin(sqrt(sin(deltaphi/2).^2+cos(phi1).*cos(phi2).*sin(deltalambda/2).^2)); + else + phi2=coasto.lat/180*pi; lambda2=coasto.long/180*pi; + deltaphi=abs(phi2-phi1); deltalambda=abs(lambda2-lambda1); + d=radius*2*asin(sqrt(sin(deltaphi/2).^2+cos(phi1).*cos(phi2).*sin(deltalambda/2).^2)); + end + distance(j)=min(d); + end + pos=find(distancemindistance_land); + distance(pos2)=mindistance_land; + + dist=min(maxdistance,distance); + md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3,'refine',md.mesh,'refinemetric',dist); + + end + + md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); + + save ./Models/SlrFarrell_Mesh md; + + plotmodel (md,'data',md.mask.ocean_levelset,'edgecolor','k'); + +end + +if any(steps==2) + disp(' Step 2: Define source as in Farrell, 1972, Figure 1'); + md = loadmodel('./Models/SlrFarrell_Mesh'); + + md.slr.sealevel=md.mask.ocean_levelset; % 1 m SLR everywhere + + md.slr.deltathickness=zeros(md.mesh.numberofelements,1); + md.slr.steric_rate=zeros(md.mesh.numberofvertices,1); + + save ./Models/SlrFarrell_Loads md; + + plotmodel (md,'data',md.slr.sealevel,'view',[90 90],'title#all','Initial sea-level [m]'); + +end + +if any(steps==3) + disp(' Step 3: Parameterization'); + md = loadmodel('./Models/SlrFarrell_Loads'); + + nlove=10001; + md.slr.love_h = love_numbers('h','CM'); md.slr.love_h(nlove+1:end)=[]; + md.slr.love_k = love_numbers('k','CM'); md.slr.love_k(nlove+1:end)=[]; + md.slr.love_l = love_numbers('l','CM'); md.slr.love_l(nlove+1:end)=[]; + + md.mask.land_levelset = 1-md.mask.ocean_levelset; + md.mask.ice_levelset = ones(md.mesh.numberofvertices,1); + md.mask.groundedice_levelset = -ones(md.mesh.numberofvertices,1); + pos=find(md.mesh.lat <-80); + md.mask.ice_levelset(pos(1))=-1; % ice yes! + md.mask.groundedice_levelset(pos(1))=1; % ice grounded! + + di=md.materials.rho_ice/md.materials.rho_water; + md.geometry.thickness=ones(md.mesh.numberofvertices,1); + md.geometry.surface=(1-di)*zeros(md.mesh.numberofvertices,1); + md.geometry.base=md.geometry.surface-md.geometry.thickness; + md.geometry.bed=md.geometry.base; + + md.initialization.temperature=273.25*ones(md.mesh.numberofvertices,1); + md.materials.rheology_B=paterson(md.initialization.temperature); + md.materials.rheology_n=3*ones(md.mesh.numberofelements,1); + + md.miscellaneous.name='SlrFarrell'; + + save ./Models/SlrFarrell_Parameterization md; + +end + +if any(steps==4) + disp(' Step 4: Solve Slr solver'); + md = loadmodel('./Models/SlrFarrell_Parameterization'); + + md.cluster=generic('name',oshostname(),'np',3); + md.verbose=verbose('111111111'); + + md.slr.reltol = 0.1/100; % per cent change in solution + + md=solve(md,'Slr'); + + save ./Models/SlrFarrell_Solution md; + +end + +if any(steps==5) + disp(' Step 5: Plot solutions'); + md = loadmodel('./Models/SlrFarrell_Solution'); + + sol = md.results.SealevelriseSolution.Sealevel*100; % per cent normalized by GMSL (which 1 m) + + res = 1; % [degree] + + [lat_grid, lon_grid] = meshgrid(linspace(-90,90,180/res), linspace(-180,180,360/res)); + sol_grid = zeros(size(lat_grid)); + + F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol); + F.Method = 'linear'; + F.ExtrapolationMethod = 'linear'; + + sol_grid = F(lat_grid, lon_grid); + sol_grid(isnan(sol_grid))=0; + sol_grid(lat_grid>85 & sol_grid==0) =NaN; + + set(0,'DefaultAxesFontSize',18,'DefaultAxesLineWidth',1,'DefaultTextFontSize',18,'DefaultLineMarkerSize',8) + figure1=figure('Position', [100, 100, 1000, 500]); + gcf; load coast; cla; + pcolor(lon_grid,lat_grid,sol_grid); shading flat; hold on; + [C,h]=contour(lon_grid,lat_grid,sol_grid,[96 98 100 102 104 105],'-k','linewidth',2); + clabel(C,h,'FontSize',18,'Color','red','LabelSpacing',500); + geoshow(lat,long,'DisplayType','polygon','FaceColor',[.82 .82 .82]); + plot(long,lat,'k'); hold off; + c1=colorbar; + colormap(flipud(haxby)); + caxis([96 105]); + xlim([-170 170]); + ylim([-85 85]); + grid on; + title('Relative sea-level [% of GMSL]'); + set(gcf,'color','w'); + %export_fig('Fig5.pdf'); + +end + Index: /issm/trunk/examples/SlrGRACE/runme.m =================================================================== --- /issm/trunk/examples/SlrGRACE/runme.m (revision 22822) +++ /issm/trunk/examples/SlrGRACE/runme.m (revision 22822) @@ -0,0 +1,305 @@ + +clear all; +addpath('../Data','../Functions'); + +steps=[1]; % [1:7] + +if any(steps==1) + disp(' Step 1: Global mesh creation'); + + numrefine=1; + resolution=150*1e3; % inital resolution [m] + radius = 6.371012*10^6; % mean radius of Earth, m + mindistance_coast=150*1e3; % coastal resolution [m] + mindistance_land=300*1e3; % resolution on the continents [m] + maxdistance=600*1e3; % max element size (on mid-oceans) [m] + + md=model; + md.mask=maskpsl(); + md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3); % attributes in [km] + + for i=1:numrefine, + + md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); + + distance=zeros(md.mesh.numberofvertices,1); + + pos=find(~md.mask.ocean_levelset); coaste.lat=md.mesh.lat(pos); coaste.long=md.mesh.long(pos); + pos=find(md.mask.ocean_levelset); coasto.lat=md.mesh.lat(pos); coasto.long=md.mesh.long(pos); + + for j=1:md.mesh.numberofvertices + phi1=md.mesh.lat(j)/180*pi; lambda1=md.mesh.long(j)/180*pi; + if md.mask.ocean_levelset(j), + phi2=coaste.lat/180*pi; lambda2=coaste.long/180*pi; + deltaphi=abs(phi2-phi1); deltalambda=abs(lambda2-lambda1); + d=radius*2*asin(sqrt(sin(deltaphi/2).^2+cos(phi1).*cos(phi2).*sin(deltalambda/2).^2)); + else + phi2=coasto.lat/180*pi; lambda2=coasto.long/180*pi; + deltaphi=abs(phi2-phi1); deltalambda=abs(lambda2-lambda1); + d=radius*2*asin(sqrt(sin(deltaphi/2).^2+cos(phi1).*cos(phi2).*sin(deltalambda/2).^2)); + end + distance(j)=min(d); + end + pos=find(distancemindistance_land); + distance(pos2)=mindistance_land; + + dist=min(maxdistance,distance); + md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3,'refine',md.mesh,'refinemetric',dist); + + end + + md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); + + save ./Models/SlrGRACE_Mesh md; + + plotmodel (md,'data',md.mask.ocean_levelset,'edgecolor','k'); + +end + +if any(steps==2) + disp(' Step 2: Define loads in meters of ice height equivalent'); + md = loadmodel('./Models/SlrGRACE_Mesh'); + + year_month = 2007+15/365; + time_range = [year_month year_month]; + + water_load = grace(md.mesh.elements,md.mesh.lat,md.mesh.long,time_range(1),time_range(2)); + + md.slr.deltathickness = water_load*md.materials.rho_freshwater/md.materials.rho_ice; + + save ./Models/SlrGRACE_Loads md; + + plotmodel (md,'data',md.slr.deltathickness,'view',[90 -90],'caxis',[-.1 .1],'title','Ice height equivalent [m]'); + +end + +if any(steps==3) + disp(' Step 3: Parameterization'); + md = loadmodel('./Models/SlrGRACE_Loads'); + + nlove=10001; + md.slr.love_h = love_numbers('h','CM'); md.slr.love_h(nlove+1:end)=[]; + md.slr.love_k = love_numbers('k','CM'); md.slr.love_k(nlove+1:end)=[]; + md.slr.love_l = love_numbers('l','CM'); md.slr.love_l(nlove+1:end)=[]; + + md.mask.groundedice_levelset = ones(md.mesh.numberofvertices,1); + md.mask.ice_levelset = ones(md.mesh.numberofvertices,1); + pos=find(md.slr.deltathickness~=0); + md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; + md.mask.land_levelset = 1-md.mask.ocean_levelset; + + md.slr.sealevel=zeros(md.mesh.numberofvertices,1); + md.slr.steric_rate=zeros(md.mesh.numberofvertices,1); + md.slr.ocean_area_scaling = 1; + + di=md.materials.rho_ice/md.materials.rho_water; + md.geometry.thickness=ones(md.mesh.numberofvertices,1); + md.geometry.surface=(1-di)*zeros(md.mesh.numberofvertices,1); + md.geometry.base=md.geometry.surface-md.geometry.thickness; + md.geometry.bed=md.geometry.base; + + md.initialization.temperature=273.25*ones(md.mesh.numberofvertices,1); + md.materials.rheology_B=paterson(md.initialization.temperature); + md.materials.rheology_n=3*ones(md.mesh.numberofelements,1); + + md.miscellaneous.name='SlrGRACE'; + + save ./Models/SlrGRACE_Parameterization md; + +end + +if any(steps==4) + disp(' Step 4: Solve Slr solver'); + md = loadmodel('./Models/SlrGRACE_Parameterization'); + + md.slr.rigid=1; + md.slr.elastic=1; + md.slr.rotation=1; + + md.cluster=generic('name',oshostname(),'np',3); + md.verbose=verbose('111111111'); + + md=solve(md,'Slr'); + + save ./Models/SlrGRACE_Solution md; + +end + +if any(steps==5) + disp(' Step 5: Plot solutions'); + md = loadmodel('./Models/SlrGRACE_Solution'); + + sol1 = md.slr.deltathickness*100; % [cm] + sol2 = md.results.SealevelriseSolution.Sealevel*1000; % [mm] + + sol_name={'Change in water equivalent height [cm]', 'Relative sea-level [mm]'}; + fig_name={'Fig_dH.pdf','Fig_slr.pdf'}; + + res = 1.0; + + [lat_grid, lon_grid] = meshgrid(linspace(-90,90,180/res), linspace(-180,180,360/res)); + sol_grid = zeros(size(lat_grid)); + + for kk=1:2 + sol=eval(sprintf('sol%d',kk)); + + if length(sol)==md.mesh.numberofelements + for jj=1:md.mesh.numberofelements + ii=(jj-1)*3; + pp(ii+1:ii+3)=md.mesh.elements(jj,:); + end + for jj=1:md.mesh.numberofvertices + pos=ceil(find(pp==jj)/3); + temp(jj)=mean(sol(pos)); + end + sol=temp'; + end + + F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol); + F.Method = 'linear'; + F.ExtrapolationMethod = 'linear'; + + sol_grid = F(lat_grid, lon_grid); + sol_grid(isnan(sol_grid))=0; + sol_grid(lat_grid>85 & sol_grid==0) = NaN; + + set(0,'DefaultAxesFontSize',18,'DefaultAxesLineWidth',1,'DefaultTextFontSize',18,'DefaultLineMarkerSize',8) + figure1=figure('Position', [100, 100, 1000, 500]); + gcf; load coast; cla; + pcolor(lon_grid,lat_grid,sol_grid); shading flat; hold on; + if (kk==1) + geoshow(flipud(lat),flipud(long),'DisplayType','polygon','FaceColor','white'); + else + geoshow(lat,long,'DisplayType','polygon','FaceColor','white'); + end + plot(long,lat,'k'); hold off; + c1=colorbar; + colormap('haxby'); + caxis([-floor(min(abs(min(sol)),abs(max(sol)))) floor(min(abs(min(sol)),abs(max(sol))))]); + xlim([-180 180]); + ylim([-90 90]); + grid on; + title(sol_name(kk)); + set(gcf,'color','w'); + %export_fig(fig_name{kk}); + end + +end + +if any(steps==6) + disp(' Step 6: Transient run'); + md = loadmodel('./Models/SlrGRACE_Parameterization'); + + disp('Projecting loads onto the mesh...'); + time_range = 2007 + [15 350]/365; + water_load = grace(md.mesh.elements,md.mesh.lat,md.mesh.long,time_range(1),time_range(2)); + + [ne,nt]=size(water_load); + md.slr.deltathickness = zeros(ne+1,nt); + md.slr.deltathickness(1:ne,:) = water_load*md.materials.rho_freshwater/md.materials.rho_ice; + md.slr.deltathickness(ne+1,:)=[1:nt]; % times + + md.transient.issmb=0; + md.transient.ismasstransport=0; + md.transient.isstressbalance=0; + md.transient.isthermal=0; + md.transient.isslr=1; + + md.timestepping.start_time=-1; + md.timestepping.final_time=nt; + md.timestepping.time_step=1; + md.timestepping.interp_forcings=0; + md.settings.output_frequency=1; + + md.cluster=generic('name',oshostname(),'np',3); + md.verbose=verbose('111111111'); + + md=solve(md,'tr'); + + save ./Models/SlrGRACE_Transient md; + +end + +if any(steps==7) + disp(' Step 7: Plot transient'); + md = loadmodel('./Models/SlrGRACE_Transient'); + + time = md.slr.deltathickness(end,:); + + for tt=1:length(time) + gmsl(tt) = md.results.TransientSolution(tt).SealevelEustatic*1000; % GMSL rate mm/yr + sol1(:,tt) = md.slr.deltathickness(1:end-1,tt)*100; % ice equivalent height [cm/yr] + sol2(:,tt) = md.results.TransientSolution(tt+1).Sealevel*1000; % mm/yr + end + sol_name = {'Change in water equivalent height [cm]', 'Relative sea-level [mm]'}; + movie_name = {'Movie_dH.avi','Movie_slr.avi'}; + + res = 1.0; + + [lat_grid, lon_grid] = meshgrid(linspace(-90,90,180/res), linspace(-180,180,360/res)); + sol_grid = zeros(size(lat_grid)); + + for kk=1:2 + sol=eval(sprintf('sol%d',kk)); + + if length(sol)==md.mesh.numberofelements + for jj=1:md.mesh.numberofelements + ii=(jj-1)*3; + pp(ii+1:ii+3)=md.mesh.elements(jj,:); + end + for jj=1:md.mesh.numberofvertices + pos=ceil(find(pp==jj)/3); + temp2(jj,:)=mean(sol(pos,:)); + end + sol=temp2; + end + + vidObj = VideoWriter(movie_name{kk}); + vidObj.FrameRate=2; % frames per second + open(vidObj); + + for jj=1:length(time) + F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol(:,jj)); + F.Method = 'linear'; + F.ExtrapolationMethod = 'linear'; + + sol_grid = F(lat_grid, lon_grid); + sol_grid(isnan(sol_grid))=0; + sol_grid(lat_grid>85 & sol_grid==0) = NaN; + + set(0,'DefaultAxesFontSize',18,'DefaultAxesLineWidth',1,'DefaultTextFontSize',18,'DefaultLineMarkerSize',8) + figure1=figure('Position', [100, 100, 1000, 500]); + gcf; load coast; cla; + pcolor(lon_grid,lat_grid,sol_grid); shading flat; hold on; + if (kk==1) + geoshow(flipud(lat),flipud(long),'DisplayType','polygon','FaceColor','white'); + else + geoshow(lat,long,'DisplayType','polygon','FaceColor','white'); + end + plot(long,lat,'k'); hold off; + c1=colorbar; + colormap('haxby'); + caxis([-floor(min(abs(min(min(sol))),abs(max(max(sol))))) floor(min(abs(min(min(sol))),abs(max(max(sol)))))]); + xlim([-180 180]); + ylim([-90 90]); + grid on; + title(sol_name(kk)); + set(gcf,'color','w'); + writeVideo(vidObj,getframe(gcf)); + close + end + close(vidObj); + end + !open *.avi; + + % plot GMSL time series + plot(time,gmsl,'-*','linewidth',3); grid on; + xlabel('# month'); + ylabel('GMSL [mm/yr]'); + set(gcf,'color','w'); + %export_fig('Fig7.pdf'); + +end + Index: /issm/trunk/jenkins/jenkins.sh =================================================================== --- /issm/trunk/jenkins/jenkins.sh (revision 22821) +++ /issm/trunk/jenkins/jenkins.sh (revision 22822) @@ -296,7 +296,7 @@ cd $ISSM_DIR/test/NightlyRun if [ "$OS" = "win" ]; then - $MATLAB_PATH/bin/matlab -nodisplay -nojvm -nosplash -nodesktop -r "addpath $ISSM_DIR_WIN/src/m/dev; devpath; addpath $ISSM_DIR_WIN/nightlylog/; matlab_run$i" -logfile $ISSM_DIR_WIN/nightlylog/matlab_log$i.log & - else - $MATLAB_PATH/bin/matlab -nodisplay -nojvm -nosplash -nodesktop -r "addpath $ISSM_DIR/src/m/dev; devpath; addpath $ISSM_DIR/nightlylog/; matlab_run$i" -logfile $ISSM_DIR/nightlylog/matlab_log$i.log & + $MATLAB_PATH/bin/matlab -nodisplay -nosplash -r "addpath $ISSM_DIR_WIN/src/m/dev; devpath; addpath $ISSM_DIR_WIN/nightlylog/; matlab_run$i" -logfile $ISSM_DIR_WIN/nightlylog/matlab_log$i.log & + else + $MATLAB_PATH/bin/matlab -nodisplay -nosplash -r "addpath $ISSM_DIR/src/m/dev; devpath; addpath $ISSM_DIR/nightlylog/; matlab_run$i" -logfile $ISSM_DIR/nightlylog/matlab_log$i.log & fi done @@ -470,5 +470,5 @@ echo 'end' >> $FILE - $MATLAB_PATH/bin/matlab -nodisplay -nojvm -nosplash -nodesktop -r "addpath $ISSM_DIR/src/m/dev; devpath; addpath $ISSM_DIR/nightlylog/; runme" -logfile $ISSM_DIR/nightlylog/$LOG_FILE + $MATLAB_PATH/bin/matlab -nodisplay -nosplash -r "addpath $ISSM_DIR/src/m/dev; devpath; addpath $ISSM_DIR/nightlylog/; runme" -logfile $ISSM_DIR/nightlylog/$LOG_FILE echo "starting: $(basename $dir)" >> $ISSM_DIR/nightlylog/matlab_log_examples.log cat $ISSM_DIR/nightlylog/$LOG_FILE >> $ISSM_DIR/nightlylog/matlab_log_examples.log Index: /issm/trunk/src/c/analyses/MasstransportAnalysis.cpp =================================================================== --- /issm/trunk/src/c/analyses/MasstransportAnalysis.cpp (revision 22821) +++ /issm/trunk/src/c/analyses/MasstransportAnalysis.cpp (revision 22822) @@ -179,8 +179,4 @@ } - if(isoceancoupling){ - iomodel->FetchDataToInput(elements,"md.mesh.lat",MeshLatEnum); - iomodel->FetchDataToInput(elements,"md.mesh.long",MeshLongEnum); - } if(!issmb){ iomodel->FetchDataToInput(elements,"md.smb.mass_balance",SmbMassBalanceEnum); Index: /issm/trunk/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk/src/c/classes/Elements/Element.cpp (revision 22821) +++ /issm/trunk/src/c/classes/Elements/Element.cpp (revision 22822) @@ -1400,5 +1400,9 @@ parameters->FindParam(&N,NULL,ControlInputSizeNEnum); parameters->FindParam(&M,NULL,ControlInputSizeMEnum); - if(control_index>0) for(int n=0;n0) { + for(int n=0;nObjectEnum()){ case TriaEnum: - values_input->AddTimeInput(new TriaInput(ControlInputValuesEnum,values,P1Enum)); + values_input->AddTimeInput(new TriaInput(input_enum,values,P1Enum)); mins_input->AddTimeInput(new TriaInput(ControlInputMinsEnum,values_min,P1Enum)); maxs_input->AddTimeInput(new TriaInput(ControlInputMaxsEnum,values_max,P1Enum)); break; case PentaEnum: - values_input->AddTimeInput(new PentaInput(ControlInputValuesEnum,values,P1Enum)); + values_input->AddTimeInput(new PentaInput(input_enum,values,P1Enum)); mins_input->AddTimeInput(new PentaInput(ControlInputMinsEnum,values_min,P1Enum)); maxs_input->AddTimeInput(new PentaInput(ControlInputMaxsEnum,values_max,P1Enum)); break; case TetraEnum: - values_input->AddTimeInput(new TetraInput(ControlInputValuesEnum,values,P1Enum)); + values_input->AddTimeInput(new TetraInput(input_enum,values,P1Enum)); mins_input->AddTimeInput(new TetraInput(ControlInputMinsEnum,values_min,P1Enum)); maxs_input->AddTimeInput(new TetraInput(ControlInputMaxsEnum,values_max,P1Enum)); Index: /issm/trunk/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/classes/Elements/Tria.cpp (revision 22821) +++ /issm/trunk/src/c/classes/Elements/Tria.cpp (revision 22822) @@ -2082,5 +2082,5 @@ /*Control Inputs*/ - if (control_analysis){ + if (control_analysis && !ad_analysis){ if(!ad_analysis)iomodel->FindConstant(&controls,NULL,"md.inversion.control_parameters"); if(ad_analysis)iomodel->FindConstant(&controls,NULL,"md.autodiff.independent_object_names"); Index: /issm/trunk/src/c/classes/ExternalResults/GenericExternalResult.h =================================================================== --- /issm/trunk/src/c/classes/ExternalResults/GenericExternalResult.h (revision 22821) +++ /issm/trunk/src/c/classes/ExternalResults/GenericExternalResult.h (revision 22822) @@ -179,4 +179,7 @@ } /*}}}*/ +void Transpose(void){ /*{{{*/ + _error_("not implemented yet"); +} /*}}}*/ char* GetResultName(void){ /*{{{*/ char* name = xNew(strlen(this->result_name)+1); @@ -594,183 +597,206 @@ } /*}}}*/ - -/*Specific instantiations for IssmDouble*: */ +template <> inline void GenericExternalResult::Transpose(void){/*{{{*/ + + + /*Perform transpose only if we have a matrix*/ + if(M>1 && N>1){ + IssmPDouble* temp=xNew(M*N); + for(int i=0;i(this->value); + this->value = temp; + } + + /*Switch dimensions*/ + int temp2 = this->N; + this->N = this->M; + this->M = temp2; + + + +} /*}}}*/ + + /*Specific instantiations for IssmDouble*: */ #if defined(_HAVE_ADOLC_) && !defined(_WRAPPERS_) //We hook off this specific specialization when not running ADOLC, otherwise we get a redeclaration with the next specialization. -template <> inline void GenericExternalResult::WriteData(FILE* fid,bool io_gather){ /*{{{*/ - - int i; - int my_rank; - int type; - int rows,cols; - char *name = NULL; - IssmPDouble passiveDouble; - IssmPDouble* passiveDoubles; - - /*recover my_rank:*/ - my_rank=IssmComm::GetRank(); - - if(io_gather){ - /*we are gathering the data on cpu 0, don't write on other cpus: */ - if(my_rank) return; - } - - /*First write enum: */ - int length=(strlen(this->result_name)+1)*sizeof(char); - fwrite(&length,sizeof(int),1,fid); - fwrite(this->result_name,length,1,fid); - - /*Now write time and step: */ - passiveDouble=reCast(time); - fwrite(&passiveDouble,sizeof(IssmPDouble),1,fid); - fwrite(&step,sizeof(int),1,fid); - - /*writing a IssmDouble array, type is 3:*/ - type=3; - fwrite(&type,sizeof(int),1,fid); - rows=this->M; - fwrite(&rows,sizeof(int),1,fid); - cols=this->N; - fwrite(&cols,sizeof(int),1,fid); - - passiveDoubles=xNew(this->M*this->N); - for (i=0;iM*this->N;i++)passiveDoubles[i]=reCast(value[i]); - fwrite(passiveDoubles,cols*rows*sizeof(IssmPDouble),1,fid); - xDelete(passiveDoubles); - -} -/*}}}*/ -#endif - -/*Specifics instantiations for Vector*/ -template <> inline GenericExternalResult*>::GenericExternalResult(int in_id, int in_enum_type,Vector* in_values, int in_M,int in_N,int in_step,IssmDouble in_time){/*{{{*/ - _error_("instanciation not correct"); -} -/*}}}*/ -template <> inline GenericExternalResult*>::GenericExternalResult(int in_id, int in_enum_type,Vector* in_value,int in_step, IssmDouble in_time){ /*{{{*/ - id = in_id; - M = 0; - N = 0; - - /*Convert enum to name*/ - EnumToStringx(&this->result_name,in_enum_type); - - step = in_step; - time = in_time; - - value = in_value; -} /*}}}*/ -template <> inline GenericExternalResult*>::~GenericExternalResult(){ /*{{{*/ - xDelete(this->result_name); - delete value; -} /*}}}*/ -template <> inline void GenericExternalResult*>::Echo(void){ /*{{{*/ - - _printf_("GenericExternalResult*>:\n"); - this->GenericEcho(); - this->value->Echo(); - -} /*}}}*/ -template <> inline void GenericExternalResult*>::DeepEcho(void){ /*{{{*/ - - this->Echo(); - -} /*}}}*/ -template <> inline Object* GenericExternalResult*>::copy(void){ /*{{{*/ - return new GenericExternalResult*>(this->id,StringToEnumx(this->result_name),this->value,this->step,this->time); -} /*}}}*/ -#if defined(_HAVE_ADOLC_) && !defined(_WRAPPERS_) //We hook off this specific specialization when not running ADOLC, otherwise we get a redeclaration with the next specialization. -template <> inline void GenericExternalResult*>::WriteData(FILE* fid,bool io_gather){ /*{{{*/ - - char *name = NULL; - int length,rows,cols=1; - - if(!io_gather){ - _error_("not supported yet"); - } - - /*Serialize vector*/ - IssmPDouble* serialvalues = this->value->ToMPISerial(); - this->value->GetSize(&rows); - - if(IssmComm::GetRank()==0){ - /*First write name: */ - length=(strlen(this->result_name)+1)*sizeof(char); + template <> inline void GenericExternalResult::WriteData(FILE* fid,bool io_gather){ /*{{{*/ + + int i; + int my_rank; + int type; + int rows,cols; + char *name = NULL; + IssmPDouble passiveDouble; + IssmPDouble* passiveDoubles; + + /*recover my_rank:*/ + my_rank=IssmComm::GetRank(); + + if(io_gather){ + /*we are gathering the data on cpu 0, don't write on other cpus: */ + if(my_rank) return; + } + + /*First write enum: */ + int length=(strlen(this->result_name)+1)*sizeof(char); fwrite(&length,sizeof(int),1,fid); fwrite(this->result_name,length,1,fid); /*Now write time and step: */ - IssmPDouble passiveDouble=reCast(time); + passiveDouble=reCast(time); fwrite(&passiveDouble,sizeof(IssmPDouble),1,fid); fwrite(&step,sizeof(int),1,fid); /*writing a IssmDouble array, type is 3:*/ - int type=3; + type=3; fwrite(&type,sizeof(int),1,fid); + rows=this->M; fwrite(&rows,sizeof(int),1,fid); + cols=this->N; fwrite(&cols,sizeof(int),1,fid); - fwrite(serialvalues,cols*rows*sizeof(IssmPDouble),1,fid); - } - - /*Clean up*/ - xDelete(serialvalues); - -} -/*}}}*/ + + passiveDoubles=xNew(this->M*this->N); + for (i=0;iM*this->N;i++)passiveDoubles[i]=reCast(value[i]); + fwrite(passiveDoubles,cols*rows*sizeof(IssmPDouble),1,fid); + xDelete(passiveDoubles); + + } + /*}}}*/ #endif -template <> inline int GenericExternalResult*>::ObjectEnum(void){ /*{{{*/ - return NoneEnum; - /*???? FIXME*/ -} /*}}}*/ - -/*Specifics instantiations for Vector*/ -template <> inline void GenericExternalResult*>::WriteData(FILE* fid,bool io_gather){ /*{{{*/ - - int i; - char *name = NULL; - int length,rows,cols=1; - IssmDouble* serialvalues = NULL; - IssmPDouble* pserialvalues = NULL; - - if(!io_gather){ - _error_("not supported yet"); - } - - /*Serialize vector*/ - serialvalues = this->value->ToMPISerial(); - this->value->GetSize(&rows); - - pserialvalues=xNew(rows); - for(i=0;i(serialvalues[i]); - - if(IssmComm::GetRank()==0){ - /*First write name: */ - length=(strlen(this->result_name)+1)*sizeof(char); - fwrite(&length,sizeof(int),1,fid); - fwrite(this->result_name,length,1,fid); - - /*Now write time and step: */ - IssmPDouble passiveDouble=reCast(time); - fwrite(&passiveDouble,sizeof(IssmPDouble),1,fid); - fwrite(&step,sizeof(int),1,fid); - - /*writing a IssmDouble array, type is 3:*/ - int type=3; - fwrite(&type,sizeof(int),1,fid); - fwrite(&rows,sizeof(int),1,fid); - fwrite(&cols,sizeof(int),1,fid); - fwrite(pserialvalues,cols*rows*sizeof(IssmPDouble),1,fid); - } - - /*Clean up*/ - xDelete(pserialvalues); - xDelete(serialvalues); - -} -/*}}}*/ -template <> inline void GenericExternalResult*>::Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){/*{{{*/ - - _error_("GenericExternalResult instantiated for type Vector* called " << result_name << " not implemented yet"); - -} /*}}}*/ + + /*Specifics instantiations for Vector*/ + template <> inline GenericExternalResult*>::GenericExternalResult(int in_id, int in_enum_type,Vector* in_values, int in_M,int in_N,int in_step,IssmDouble in_time){/*{{{*/ + _error_("instanciation not correct"); + } + /*}}}*/ + template <> inline GenericExternalResult*>::GenericExternalResult(int in_id, int in_enum_type,Vector* in_value,int in_step, IssmDouble in_time){ /*{{{*/ + id = in_id; + M = 0; + N = 0; + + /*Convert enum to name*/ + EnumToStringx(&this->result_name,in_enum_type); + + step = in_step; + time = in_time; + + value = in_value; + } /*}}}*/ + template <> inline GenericExternalResult*>::~GenericExternalResult(){ /*{{{*/ + xDelete(this->result_name); + delete value; + } /*}}}*/ + template <> inline void GenericExternalResult*>::Echo(void){ /*{{{*/ + + _printf_("GenericExternalResult*>:\n"); + this->GenericEcho(); + this->value->Echo(); + + } /*}}}*/ + template <> inline void GenericExternalResult*>::DeepEcho(void){ /*{{{*/ + + this->Echo(); + + } /*}}}*/ + template <> inline Object* GenericExternalResult*>::copy(void){ /*{{{*/ + return new GenericExternalResult*>(this->id,StringToEnumx(this->result_name),this->value,this->step,this->time); + } /*}}}*/ +#if defined(_HAVE_ADOLC_) && !defined(_WRAPPERS_) //We hook off this specific specialization when not running ADOLC, otherwise we get a redeclaration with the next specialization. + template <> inline void GenericExternalResult*>::WriteData(FILE* fid,bool io_gather){ /*{{{*/ + + char *name = NULL; + int length,rows,cols=1; + + if(!io_gather){ + _error_("not supported yet"); + } + + /*Serialize vector*/ + IssmPDouble* serialvalues = this->value->ToMPISerial(); + this->value->GetSize(&rows); + + if(IssmComm::GetRank()==0){ + /*First write name: */ + length=(strlen(this->result_name)+1)*sizeof(char); + fwrite(&length,sizeof(int),1,fid); + fwrite(this->result_name,length,1,fid); + + /*Now write time and step: */ + IssmPDouble passiveDouble=reCast(time); + fwrite(&passiveDouble,sizeof(IssmPDouble),1,fid); + fwrite(&step,sizeof(int),1,fid); + + /*writing a IssmDouble array, type is 3:*/ + int type=3; + fwrite(&type,sizeof(int),1,fid); + fwrite(&rows,sizeof(int),1,fid); + fwrite(&cols,sizeof(int),1,fid); + fwrite(serialvalues,cols*rows*sizeof(IssmPDouble),1,fid); + } + + /*Clean up*/ + xDelete(serialvalues); + + } + /*}}}*/ +#endif + template <> inline int GenericExternalResult*>::ObjectEnum(void){ /*{{{*/ + return NoneEnum; + /*???? FIXME*/ + } /*}}}*/ + + /*Specifics instantiations for Vector*/ + template <> inline void GenericExternalResult*>::WriteData(FILE* fid,bool io_gather){ /*{{{*/ + + int i; + char *name = NULL; + int length,rows,cols=1; + IssmDouble* serialvalues = NULL; + IssmPDouble* pserialvalues = NULL; + + if(!io_gather){ + _error_("not supported yet"); + } + + /*Serialize vector*/ + serialvalues = this->value->ToMPISerial(); + this->value->GetSize(&rows); + + pserialvalues=xNew(rows); + for(i=0;i(serialvalues[i]); + + if(IssmComm::GetRank()==0){ + /*First write name: */ + length=(strlen(this->result_name)+1)*sizeof(char); + fwrite(&length,sizeof(int),1,fid); + fwrite(this->result_name,length,1,fid); + + /*Now write time and step: */ + IssmPDouble passiveDouble=reCast(time); + fwrite(&passiveDouble,sizeof(IssmPDouble),1,fid); + fwrite(&step,sizeof(int),1,fid); + + /*writing a IssmDouble array, type is 3:*/ + int type=3; + fwrite(&type,sizeof(int),1,fid); + fwrite(&rows,sizeof(int),1,fid); + fwrite(&cols,sizeof(int),1,fid); + fwrite(pserialvalues,cols*rows*sizeof(IssmPDouble),1,fid); + } + + /*Clean up*/ + xDelete(pserialvalues); + xDelete(serialvalues); + + } + /*}}}*/ + template <> inline void GenericExternalResult*>::Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){/*{{{*/ + + _error_("GenericExternalResult instantiated for type Vector* called " << result_name << " not implemented yet"); + + } /*}}}*/ #endif /* _EXTERNAL_RESULTOBJECT_H */ Index: /issm/trunk/src/c/classes/Inputs/ControlInput.cpp =================================================================== --- /issm/trunk/src/c/classes/Inputs/ControlInput.cpp (revision 22821) +++ /issm/trunk/src/c/classes/Inputs/ControlInput.cpp (revision 22822) @@ -330,5 +330,5 @@ TransientInput* transient_input = xDynamicCast(input); IssmDouble time = transient_input->GetTimeByOffset(timeoffset); - TransientInput* new_trans_input = new TransientInput(ControlInputValuesEnum); + TransientInput* new_trans_input = new TransientInput(this->enum_type); for(int i=0;inumtimesteps;i++){ if(transient_input->timesteps[i]==time) new_trans_input->AddTimeInput(xDynamicCast(in_input),time); Index: /issm/trunk/src/c/classes/Vertex.cpp =================================================================== --- /issm/trunk/src/c/classes/Vertex.cpp (revision 22821) +++ /issm/trunk/src/c/classes/Vertex.cpp (revision 22822) @@ -32,4 +32,9 @@ this->domaintype = iomodel->domaintype; + if(iomodel->Data("md.mesh.lat") && iomodel->Data("md.mesh.long")){ + this->latitute = iomodel->Data("md.mesh.lat")[i]; + this->longitude = iomodel->Data("md.mesh.long")[i]; + } + switch(iomodel->domaintype){ case Domain3DEnum: @@ -38,4 +43,5 @@ break; case Domain3DsurfaceEnum: + _assert_(iomodel->Data("md.mesh.lat") && iomodel->Data("md.mesh.long") && iomodel->Data("md.mesh.r")); this->latitute = iomodel->Data("md.mesh.lat")[i]; this->longitude = iomodel->Data("md.mesh.long")[i]; @@ -195,21 +201,4 @@ /*}}}*/ int Vertex::Sid(void){ return sid; }/*{{{*/ -/*}}}*/ -void Vertex::ToXYZ(Matrix* matrix){/*{{{*/ - - IssmDouble xyz[3]; - int indices[3]; - - if (this->clone==true) return; - - xyz[0]=x; - xyz[1]=y; - xyz[2]=z; - indices[0]=0; - indices[1]=1; - indices[2]=2; - - matrix->SetValues(1,&sid,3,&indices[0],&xyz[0],INS_VAL); -} /*}}}*/ void Vertex::UpdateClonePids(int* alltruepids){/*{{{*/ Index: /issm/trunk/src/c/classes/Vertex.h =================================================================== --- /issm/trunk/src/c/classes/Vertex.h (revision 22821) +++ /issm/trunk/src/c/classes/Vertex.h (revision 22822) @@ -62,5 +62,4 @@ void ShowTruePids(int* borderpids); int Sid(void); - void ToXYZ(Matrix* matrix); void UpdateClonePids(int* allborderpids); void UpdatePosition(Vector* vx,Vector* vy,Vector* vz,Parameters* parameters,IssmDouble* thickness,IssmDouble* bed); Index: /issm/trunk/src/c/classes/Vertices.cpp =================================================================== --- /issm/trunk/src/c/classes/Vertices.cpp (revision 22821) +++ /issm/trunk/src/c/classes/Vertices.cpp (revision 22822) @@ -180,44 +180,39 @@ } /*}}}*/ -IssmDouble* Vertices::ToXYZ(void){/*{{{*/ - - /*intermediary: */ - int i; - int my_rank; - int num_vertices; +void Vertices::LatLonList(IssmDouble** plat,IssmDouble** plon){/*{{{*/ /*output: */ - Matrix* xyz = NULL; IssmDouble* xyz_serial=NULL; /*recover my_rank:*/ - my_rank=IssmComm::GetRank(); + int my_rank=IssmComm::GetRank(); /*First, figure out number of vertices: */ - num_vertices=this->NumberOfVertices(); - - /*Now, allocate matrix to hold all the vertices x,y and z values: */ - xyz= new Matrix(num_vertices,3); + int num_vertices=this->NumberOfVertices(); + + /*Now, allocate vectors*/ + Vector* lat = new Vector(num_vertices); + Vector* lon = new Vector(num_vertices); /*Go through vertices, and for each vertex, object, report it cpu: */ - for(i=0;iSize();i++){ - - /*let vertex fill matrix: */ - Vertex* vertex=xDynamicCast(this->GetObjectByOffset(i)); - vertex->ToXYZ(xyz); + for(int i=0;iSize();i++){ + Vertex* vertex=xDynamicCast(this->GetObjectByOffset(i)); + lat->SetValue(vertex->sid,vertex->GetLatitude() ,INS_VAL); + lon->SetValue(vertex->sid,vertex->GetLongitude(),INS_VAL); } /*Assemble:*/ - xyz->Assemble(); + lat->Assemble(); + lon->Assemble(); /*gather on cpu 0: */ - xyz_serial=xyz->ToSerial(); + IssmDouble* lat_serial=lat->ToMPISerial(); + IssmDouble* lon_serial=lon->ToMPISerial(); /*free ressources: */ - delete xyz; - if(my_rank!=0)delete xyz_serial; - - /*return matrix: */ - return xyz_serial; -} -/*}}}*/ + *plat = lat_serial; + *plon = lon_serial; + delete lat; + delete lon; +} +/*}}}*/ Index: /issm/trunk/src/c/classes/Vertices.h =================================================================== --- /issm/trunk/src/c/classes/Vertices.h (revision 22821) +++ /issm/trunk/src/c/classes/Vertices.h (revision 22822) @@ -25,5 +25,5 @@ int NumberOfVertices(void); void Ranks(int* ranks); - IssmDouble* ToXYZ(void); + void LatLonList(IssmDouble** lat,IssmDouble** lon); }; Index: /issm/trunk/src/c/cores/controladm1qn3_core.cpp =================================================================== --- /issm/trunk/src/c/cores/controladm1qn3_core.cpp (revision 22821) +++ /issm/trunk/src/c/cores/controladm1qn3_core.cpp (revision 22822) @@ -71,4 +71,5 @@ int* N = NULL; int N_add = 0; + int* control_enum = NULL; if (solution_type == TransientSolutionEnum){ @@ -88,4 +89,5 @@ femmodel->parameters->FindParam(&scaling_factors,NULL,InversionControlScalingFactorsEnum); femmodel->parameters->FindParam(&N,NULL,ControlInputSizeNEnum); + femmodel->parameters->FindParam(&control_enum,NULL,InversionControlParametersEnum); numberofvertices=femmodel->vertices->NumberOfVertices(); @@ -104,4 +106,5 @@ if(X[index]parameters->FindParamAndMakePassive(>tol,InversionGttolEnum); femmodel->parameters->FindParam(&scaling_factors,NULL,InversionControlScalingFactorsEnum); + femmodel->parameters->FindParam(&control_enum,NULL,InversionControlParametersEnum); femmodel->parameters->SetParam(false,SaveResultsEnum); femmodel->parameters->FindParam(&N,NULL,ControlInputSizeNEnum); @@ -478,8 +484,23 @@ int step = 1; femmodel->parameters->SetParam(step,StepEnum); - femmodel->OutputControlsx(&femmodel->results); femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,JEnum,mystruct.Jlist,(*mystruct.i),mystruct.N,1,0)); - femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,Gradient1Enum,G,numberofvertices,intn/numberofvertices,1,0)); - femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,InversionControlParametersEnum,X,numberofvertices,intn/numberofvertices,1,0)); + + int offset = 0; + for(int i=0;i* G_output = new GenericExternalResult(femmodel->results->Size()+1,Gradient1Enum+i,&G[offset],N[i],numberofvertices,1,0.); + GenericExternalResult* X_output = new GenericExternalResult(femmodel->results->Size()+1,control_enum[i],&X[offset],N[i],numberofvertices,1,0.); + + /*transpose for consistency with MATLAB's formating*/ + G_output->Transpose(); + X_output->Transpose(); + + /*Add to results*/ + femmodel->results->AddObject(G_output); + femmodel->results->AddObject(X_output); + + offset += N[i]*numberofvertices; + } } else{ Index: /issm/trunk/src/c/cores/sealevelrise_core.cpp =================================================================== --- /issm/trunk/src/c/cores/sealevelrise_core.cpp (revision 22821) +++ /issm/trunk/src/c/cores/sealevelrise_core.cpp (revision 22822) @@ -15,4 +15,5 @@ Vector *Sg_absolute = NULL; Vector *Sg_eustatic = NULL; + Vector *slr_initial = NULL; Vector *steric_rate_g = NULL; Vector *U_radial = NULL; @@ -80,4 +81,11 @@ if(isslr){ Sg_eustatic=sealevelrise_core_eustatic(femmodel); //generalized eustatic (Farrel and Clark, Eq 4, 1st, 3rd and 4rd terms on the RHS. + + /* The following is for reproducing Farrell&Clark1976 Fig. 1, and aimed for the workshop! + * md.slr.sealevel is considered as the "initial sea-level" where 1 m slr is distributed uniformly over the ocean + * similar strategy may work for computin SAL due to "internal mass distribution" associated with ocean dynamics + * if it breaks the code, it will be reverted. And, we will strategize how we want to accomodate */ + GetVectorFromInputsx(&slr_initial,femmodel,SealevelEnum,VertexSIdEnum); + Sg_eustatic->AXPY(slr_initial,1); Sg=sealevelrise_core_noneustatic(femmodel,Sg_eustatic); //ocean loading tems (2nd and 5th terms on the RHS of Farrel and Clark) Index: /issm/trunk/src/c/cores/transient_core.cpp =================================================================== --- /issm/trunk/src/c/cores/transient_core.cpp (revision 22821) +++ /issm/trunk/src/c/cores/transient_core.cpp (revision 22822) @@ -21,5 +21,5 @@ /*parameters: */ IssmDouble finaltime,dt,yts,starttime; - bool isstressbalance,ismasstransport,issmb,isFS,isthermal,isgroundingline,isgia,isslr,iscoupler,ismovingfront,isdamageevolution,ishydrology,isoceancoupling,iscontrol,isautodiff; + bool isstressbalance,ismasstransport,issmb,isFS,isthermal,isgroundingline,isgia,isesa,isslr,iscoupler,ismovingfront,isdamageevolution,ishydrology,isoceancoupling,iscontrol,isautodiff; bool save_results,dakota_analysis; int timestepping; @@ -54,4 +54,5 @@ femmodel->parameters->FindParam(&isthermal,TransientIsthermalEnum); femmodel->parameters->FindParam(&isgia,TransientIsgiaEnum); + femmodel->parameters->FindParam(&isesa,TransientIsesaEnum); femmodel->parameters->FindParam(&isslr,TransientIsslrEnum); femmodel->parameters->FindParam(&iscoupler,TransientIscouplerEnum); @@ -93,4 +94,5 @@ tomitgcmcomm=parcom->GetParameterValue(); if(my_rank==0){ + /*Recover fixed parameters and store them*/ ISSM_MPI_Send(&coupling_time,1,ISSM_MPI_DOUBLE,0,10001000,tomitgcmcomm); ISSM_MPI_Recv(&oceangridnxsize,1,ISSM_MPI_INT,0,10001003,tomitgcmcomm,&status); @@ -102,5 +104,8 @@ oceangridy = xNew(oceangridnxsize*oceangridnysize); ISSM_MPI_Recv(oceangridy,oceangridnxsize*oceangridnysize,ISSM_MPI_DOUBLE,0,10001006,tomitgcmcomm,&status); - + femmodel->parameters->SetParam(oceangridx,oceangridnxsize*oceangridnysize,OceanGridXEnum); + femmodel->parameters->SetParam(oceangridy,oceangridnxsize*oceangridnysize,OceanGridYEnum); + + /*Exchange varying parameters for the initialization*/ ISSM_MPI_Send(&time,1,ISSM_MPI_DOUBLE,0,10001001,tomitgcmcomm); ISSM_MPI_Recv(&oceantime,1,ISSM_MPI_DOUBLE,0,10001002,tomitgcmcomm,&status); @@ -162,11 +167,8 @@ femmodel->parameters->SetParam(save_results,SaveResultsEnum); - if(isoceancoupling){ {{{ + if(isoceancoupling){ /*{{{*/ + #ifndef _HAVE_ADOLC_ if(VerboseSolution()) _printf0_(" ocean coupling: exchanging information\n"); int my_rank; - int oceangridnxsize,oceangridnysize; - IssmDouble *oceanmelt; - IssmDouble *icebase; - IssmDouble oceantime; ISSM_MPI_Comm tomitgcmcomm; ISSM_MPI_Status status; @@ -177,23 +179,75 @@ tomitgcmcomm=parcom->GetParameterValue(); if(my_rank==0){ + int ngrids_ocean, nels_ocean; + IssmDouble oceantime; + IssmDouble rho_ice; + IssmDouble *oceanmelt; + IssmDouble *base_oceangrid; + IssmDouble *oceangridx; + IssmDouble *oceangridy; + IssmDouble* x_ice = NULL; + IssmDouble* y_ice = NULL; + IssmDouble* lat_ice = NULL; + IssmDouble* lon_ice = NULL; + IssmDouble* z_ice = NULL; + IssmDouble* icebase= NULL; + IssmDouble* melt_mesh = NULL; + int* index_ice= NULL; + int* index_ocean = NULL; + int ngrids_ice=femmodel->vertices->NumberOfVertices(); + int nels_ice=femmodel->elements->NumberOfElements(); + + /*Recover mesh and inputs needed*/ + femmodel->parameters->FindParam(&rho_ice,MaterialsRhoIceEnum); + femmodel->GetMesh(femmodel->vertices,femmodel->elements,&x_ice,&y_ice,&z_ice,&index_ice); + femmodel->parameters->FindParam(&oceangridx,&ngrids_ocean,OceanGridXEnum); + femmodel->parameters->FindParam(&oceangridy,&ngrids_ocean,OceanGridYEnum); + BamgTriangulatex(&index_ocean,&nels_ocean,oceangridx,oceangridy,ngrids_ocean); + + femmodel->vertices->LatLonList(&lat_ice,&lon_ice); + + /*Interpolate ice base onto ocean grid*/ + GetVectorFromInputsx(&icebase,femmodel,BaseEnum,VertexSIdEnum); + InterpFromMeshToMesh2dx(&base_oceangrid,index_ice,lon_ice,lat_ice,ngrids_ice,nels_ice, + icebase,ngrids_ice,1, + oceangridx,oceangridy,ngrids_ocean,NULL); + + /*Send and receive data*/ ISSM_MPI_Send(&time,1,ISSM_MPI_DOUBLE,0,10001001,tomitgcmcomm); + oceanmelt = xNew(ngrids_ocean); ISSM_MPI_Recv(&oceantime,1,ISSM_MPI_DOUBLE,0,10001002,tomitgcmcomm,&status); if((oceantime - time > 0.1*yts) & (oceantime - time < -0.1*yts)) _error_("Ocean and ice time are starting to diverge"); - femmodel->parameters->FindParam(&oceangridnxsize,OceanGridNxEnum); - femmodel->parameters->FindParam(&oceangridnysize,OceanGridNyEnum); - oceanmelt = xNew(oceangridnxsize*oceangridnysize); - ISSM_MPI_Recv(oceanmelt,oceangridnxsize*oceangridnysize,ISSM_MPI_DOUBLE,0,10001007,tomitgcmcomm,&status); - icebase = xNew(oceangridnxsize*oceangridnysize); - for(int i=0;i(index_ice); + xDelete(index_ocean); + xDelete(lat_ice); + xDelete(lon_ice); + xDelete(x_ice); + xDelete(y_ice); + xDelete(z_ice); + xDelete(melt_mesh); + xDelete(oceangridx); + xDelete(oceangridy); xDelete(oceanmelt); xDelete(icebase); - } - } - }}} + xDelete(base_oceangrid); + } + #else + _error_("not supported"); + #endif + } + /*}}}*/ if(isthermal && domaintype==Domain3DEnum){ @@ -258,4 +312,7 @@ } + /*esa: */ + if(isesa) esa_core(femmodel); + /*Sea level rise: */ if(isslr || iscoupler) sealevelrise_core(femmodel); Index: /issm/trunk/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp (revision 22821) +++ /issm/trunk/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp (revision 22822) @@ -16,4 +16,5 @@ bool dakota_analysis; bool adolc_analysis; + bool isoceancoupling; int nnat,dummy; int* nature=NULL; @@ -24,4 +25,5 @@ iomodel->FindConstant(&materials_type,"md.materials.type"); iomodel->FindConstant(&adolc_analysis,"md.autodiff.isautodiff"); + iomodel->FindConstant(&isoceancoupling,"md.transient.isoceancoupling"); /*Did we already create the elements? : */ @@ -232,4 +234,5 @@ if (iomodel->domaintype == Domain3DsurfaceEnum) iomodel->FetchData(3,"md.mesh.lat","md.mesh.long","md.mesh.r"); else iomodel->FetchDataToInput(elements,"md.mesh.scale_factor",MeshScaleFactorEnum,1.); + if (isoceancoupling) iomodel->FetchData(2,"md.mesh.lat","md.mesh.long"); CreateNumberNodeToElementConnectivity(iomodel,solution_type); @@ -242,3 +245,4 @@ iomodel->DeleteData(6,"md.mesh.x","md.mesh.y","md.mesh.z","md.geometry.base","md.geometry.thickness","md.mask.ice_levelset"); if (iomodel->domaintype == Domain3DsurfaceEnum) iomodel->DeleteData(3,"md.mesh.lat","md.mesh.long","md.mesh.r"); + if (isoceancoupling) iomodel->DeleteData(2,"md.mesh.lat","md.mesh.long"); } Index: /issm/trunk/src/c/modules/ModelProcessorx/CreateOutputDefinitions.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/CreateOutputDefinitions.cpp (revision 22821) +++ /issm/trunk/src/c/modules/ModelProcessorx/CreateOutputDefinitions.cpp (revision 22822) @@ -341,5 +341,4 @@ iomodel->FetchMultipleData(&cfsurfacelogvel_vyobs_string,&num_cfsurfacelogvels, "md.cfsurfacelogvel.vyobs_string"); iomodel->FetchMultipleData(&cfsurfacelogvel_weights,&cfsurfacelogvel_weights_M,&cfsurfacelogvel_weights_N,&num_cfsurfacelogvels, "md.cfsurfacelogvel.weights"); iomodel->FetchMultipleData(&cfsurfacelogvel_weightstring,&num_cfsurfacelogvels, "md.cfsurfacelogvel.weights_string"); - _printf_("Num with weight string: "<FetchMultipleData(&cfsurfacelogvel_datatime,&num_cfsurfacelogvels, "md.cfsurfacelogvel.datatime"); Index: /issm/trunk/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 22821) +++ /issm/trunk/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 22822) @@ -291,4 +291,7 @@ iomodel->FindConstant(&materialtype,"md.materials.type"); + if(materialtype==MatdamageiceEnum || materialtype==MaticeEnum || materialtype==MatenhancediceEnum){ + parameters->AddObject(iomodel->CopyConstantObject("md.materials.rho_ice",MaterialsRhoIceEnum)); + } if(materialtype==MatdamageiceEnum){ iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.damage.requested_outputs"); Index: /issm/trunk/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk/src/c/shared/Enum/EnumDefinitions.h (revision 22821) +++ /issm/trunk/src/c/shared/Enum/EnumDefinitions.h (revision 22822) @@ -200,4 +200,5 @@ MasstransportRequestedOutputsEnum, MasstransportStabilizationEnum, + MaterialsRhoIceEnum, MeltingOffsetEnum, MeshAverageVertexConnectivityEnum, @@ -834,5 +835,4 @@ MaterialsRheologyLawEnum, MaterialsRhoFreshwaterEnum, - MaterialsRhoIceEnum, MaterialsRhoSeawaterEnum, MaterialsTemperateiceconductivityEnum, @@ -855,6 +855,4 @@ MeltingAnalysisEnum, MeshElementsEnum, - MeshLatEnum, - MeshLongEnum, MeshXEnum, MeshYEnum, Index: /issm/trunk/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk/src/c/shared/Enum/EnumToStringx.cpp (revision 22821) +++ /issm/trunk/src/c/shared/Enum/EnumToStringx.cpp (revision 22822) @@ -208,4 +208,5 @@ case MasstransportRequestedOutputsEnum : return "MasstransportRequestedOutputs"; case MasstransportStabilizationEnum : return "MasstransportStabilization"; + case MaterialsRhoIceEnum : return "MaterialsRhoIce"; case MeltingOffsetEnum : return "MeltingOffset"; case MeshAverageVertexConnectivityEnum : return "MeshAverageVertexConnectivity"; @@ -838,5 +839,4 @@ case MaterialsRheologyLawEnum : return "MaterialsRheologyLaw"; case MaterialsRhoFreshwaterEnum : return "MaterialsRhoFreshwater"; - case MaterialsRhoIceEnum : return "MaterialsRhoIce"; case MaterialsRhoSeawaterEnum : return "MaterialsRhoSeawater"; case MaterialsTemperateiceconductivityEnum : return "MaterialsTemperateiceconductivity"; @@ -859,6 +859,4 @@ case MeltingAnalysisEnum : return "MeltingAnalysis"; case MeshElementsEnum : return "MeshElements"; - case MeshLatEnum : return "MeshLat"; - case MeshLongEnum : return "MeshLong"; case MeshXEnum : return "MeshX"; case MeshYEnum : return "MeshY"; Index: /issm/trunk/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk/src/c/shared/Enum/StringToEnumx.cpp (revision 22821) +++ /issm/trunk/src/c/shared/Enum/StringToEnumx.cpp (revision 22822) @@ -211,4 +211,5 @@ else if (strcmp(name,"MasstransportRequestedOutputs")==0) return MasstransportRequestedOutputsEnum; else if (strcmp(name,"MasstransportStabilization")==0) return MasstransportStabilizationEnum; + else if (strcmp(name,"MaterialsRhoIce")==0) return MaterialsRhoIceEnum; else if (strcmp(name,"MeltingOffset")==0) return MeltingOffsetEnum; else if (strcmp(name,"MeshAverageVertexConnectivity")==0) return MeshAverageVertexConnectivityEnum; @@ -259,9 +260,9 @@ else if (strcmp(name,"SealevelriseTidalLoveK")==0) return SealevelriseTidalLoveKEnum; else if (strcmp(name,"SealevelriseTransitions")==0) return SealevelriseTransitionsEnum; - else if (strcmp(name,"SealevelriseUElastic")==0) return SealevelriseUElasticEnum; else stage=3; } if(stage==3){ - if (strcmp(name,"SettingsIoGather")==0) return SettingsIoGatherEnum; + if (strcmp(name,"SealevelriseUElastic")==0) return SealevelriseUElasticEnum; + else if (strcmp(name,"SettingsIoGather")==0) return SettingsIoGatherEnum; else if (strcmp(name,"SettingsOutputFrequency")==0) return SettingsOutputFrequencyEnum; else if (strcmp(name,"SettingsRecordingFrequency")==0) return SettingsRecordingFrequencyEnum; @@ -382,9 +383,9 @@ else if (strcmp(name,"BasalforcingsPicoBasinId")==0) return BasalforcingsPicoBasinIdEnum; else if (strcmp(name,"BasalforcingsPicoBoxId")==0) return BasalforcingsPicoBoxIdEnum; - else if (strcmp(name,"BasalforcingsPicoSubShelfOceanOverturning")==0) return BasalforcingsPicoSubShelfOceanOverturningEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"BasalforcingsPicoSubShelfOceanSalinity")==0) return BasalforcingsPicoSubShelfOceanSalinityEnum; + if (strcmp(name,"BasalforcingsPicoSubShelfOceanOverturning")==0) return BasalforcingsPicoSubShelfOceanOverturningEnum; + else if (strcmp(name,"BasalforcingsPicoSubShelfOceanSalinity")==0) return BasalforcingsPicoSubShelfOceanSalinityEnum; else if (strcmp(name,"BasalforcingsPicoSubShelfOceanTemp")==0) return BasalforcingsPicoSubShelfOceanTempEnum; else if (strcmp(name,"Base")==0) return BaseEnum; @@ -505,9 +506,9 @@ else if (strcmp(name,"Sealevel")==0) return SealevelEnum; else if (strcmp(name,"SealevelriseDeltathickness")==0) return SealevelriseDeltathicknessEnum; - else if (strcmp(name,"SedimentHead")==0) return SedimentHeadEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"SedimentHeadOld")==0) return SedimentHeadOldEnum; + if (strcmp(name,"SedimentHead")==0) return SedimentHeadEnum; + else if (strcmp(name,"SedimentHeadOld")==0) return SedimentHeadOldEnum; else if (strcmp(name,"SedimentHeadResidual")==0) return SedimentHeadResidualEnum; else if (strcmp(name,"SigmaNN")==0) return SigmaNNEnum; @@ -628,9 +629,9 @@ else if (strcmp(name,"WaterfractionDrainageIntegrated")==0) return WaterfractionDrainageIntegratedEnum; else if (strcmp(name,"Waterfraction")==0) return WaterfractionEnum; - else if (strcmp(name,"Waterheight")==0) return WaterheightEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"WeightsSurfaceObservation")==0) return WeightsSurfaceObservationEnum; + if (strcmp(name,"Waterheight")==0) return WaterheightEnum; + else if (strcmp(name,"WeightsSurfaceObservation")==0) return WeightsSurfaceObservationEnum; else if (strcmp(name,"InputsEND")==0) return InputsENDEnum; else if (strcmp(name,"Absolute")==0) return AbsoluteEnum; @@ -751,9 +752,9 @@ else if (strcmp(name,"GaussTetra")==0) return GaussTetraEnum; else if (strcmp(name,"GaussTria")==0) return GaussTriaEnum; - else if (strcmp(name,"GenericOption")==0) return GenericOptionEnum; else stage=7; } if(stage==7){ - if (strcmp(name,"GenericParam")==0) return GenericParamEnum; + if (strcmp(name,"GenericOption")==0) return GenericOptionEnum; + else if (strcmp(name,"GenericParam")==0) return GenericParamEnum; else if (strcmp(name,"GiadWdt")==0) return GiadWdtEnum; else if (strcmp(name,"GiaIvinsAnalysis")==0) return GiaIvinsAnalysisEnum; @@ -856,5 +857,4 @@ else if (strcmp(name,"MaterialsRheologyLaw")==0) return MaterialsRheologyLawEnum; else if (strcmp(name,"MaterialsRhoFreshwater")==0) return MaterialsRhoFreshwaterEnum; - else if (strcmp(name,"MaterialsRhoIce")==0) return MaterialsRhoIceEnum; else if (strcmp(name,"MaterialsRhoSeawater")==0) return MaterialsRhoSeawaterEnum; else if (strcmp(name,"MaterialsTemperateiceconductivity")==0) return MaterialsTemperateiceconductivityEnum; @@ -880,6 +880,4 @@ if (strcmp(name,"MeltingAnalysis")==0) return MeltingAnalysisEnum; else if (strcmp(name,"MeshElements")==0) return MeshElementsEnum; - else if (strcmp(name,"MeshLat")==0) return MeshLatEnum; - else if (strcmp(name,"MeshLong")==0) return MeshLongEnum; else if (strcmp(name,"MeshX")==0) return MeshXEnum; else if (strcmp(name,"MeshY")==0) return MeshYEnum; @@ -998,10 +996,10 @@ else if (strcmp(name,"Outputdefinition90")==0) return Outputdefinition90Enum; else if (strcmp(name,"Outputdefinition91")==0) return Outputdefinition91Enum; + else if (strcmp(name,"Outputdefinition92")==0) return Outputdefinition92Enum; + else if (strcmp(name,"Outputdefinition93")==0) return Outputdefinition93Enum; else stage=9; } if(stage==9){ - if (strcmp(name,"Outputdefinition92")==0) return Outputdefinition92Enum; - else if (strcmp(name,"Outputdefinition93")==0) return Outputdefinition93Enum; - else if (strcmp(name,"Outputdefinition94")==0) return Outputdefinition94Enum; + if (strcmp(name,"Outputdefinition94")==0) return Outputdefinition94Enum; else if (strcmp(name,"Outputdefinition95")==0) return Outputdefinition95Enum; else if (strcmp(name,"Outputdefinition96")==0) return Outputdefinition96Enum; @@ -1121,10 +1119,10 @@ else if (strcmp(name,"TransientParam")==0) return TransientParamEnum; else if (strcmp(name,"TransientSolution")==0) return TransientSolutionEnum; + else if (strcmp(name,"Tria")==0) return TriaEnum; + else if (strcmp(name,"TriaInput")==0) return TriaInputEnum; else stage=10; } if(stage==10){ - if (strcmp(name,"Tria")==0) return TriaEnum; - else if (strcmp(name,"TriaInput")==0) return TriaInputEnum; - else if (strcmp(name,"UzawaPressureAnalysis")==0) return UzawaPressureAnalysisEnum; + if (strcmp(name,"UzawaPressureAnalysis")==0) return UzawaPressureAnalysisEnum; else if (strcmp(name,"VectorParam")==0) return VectorParamEnum; else if (strcmp(name,"Vertex")==0) return VertexEnum; Index: /issm/trunk/src/c/shared/io/Marshalling/IoCodeConversions.cpp =================================================================== --- /issm/trunk/src/c/shared/io/Marshalling/IoCodeConversions.cpp (revision 22821) +++ /issm/trunk/src/c/shared/io/Marshalling/IoCodeConversions.cpp (revision 22822) @@ -16,4 +16,9 @@ const char* field = "md.geometry.thickness"; input_enum = ThicknessEnum; + fieldname=xNew((strlen(field)+1)); xMemCpy(fieldname,field,(strlen(field)+1)); + } + else if(strcmp(string_in,"MaterialsRheologyBbar")==0){ + const char* field = "md.materials.rheology_B"; + input_enum = MaterialsRheologyBbarEnum; fieldname=xNew((strlen(field)+1)); xMemCpy(fieldname,field,(strlen(field)+1)); } @@ -141,4 +146,9 @@ const char* field = "md.calving.stress_threshold_groundedice"; input_enum = CalvingStressThresholdGroundediceEnum; + fieldname=xNew((strlen(field)+1)); xMemCpy(fieldname,field,(strlen(field)+1)); + } + else if(strcmp(string_in,"DamageDbar")==0){ + const char* field = "md.damage.D"; + input_enum = DamageDbarEnum; fieldname=xNew((strlen(field)+1)); xMemCpy(fieldname,field,(strlen(field)+1)); } Index: /issm/trunk/src/m/classes/clusters/generic.m =================================================================== --- /issm/trunk/src/m/classes/clusters/generic.m (revision 22821) +++ /issm/trunk/src/m/classes/clusters/generic.m (revision 22822) @@ -219,5 +219,11 @@ fid=fopen([modelname '.queue'],'w'); fprintf(fid,'#!%s\n',cluster.shell); - fprintf(fid,'mpiexec -np %i %s/%s %s %s %s : -np %i ./mitgcmuv\n',cluster.np,cluster.codepath,executable,solution,cluster.executionpath,modelname,cluster.npocean); + if ~isvalgrind, + fprintf(fid,'mpiexec -np %i %s/%s %s %s %s : -np %i ./mitgcmuv\n',cluster.np,cluster.codepath,executable,solution,cluster.executionpath,modelname,cluster.npocean); + + else + fprintf(fid,'mpiexec -np %i %s --leak-check=full --error-limit=no --dsymutil=yes --suppressions=%s %s/%s %s %s %s : -np %i ./mitgcmuv\n',... + cluster.np,cluster.valgrind,cluster.valgrindsup,cluster.codepath,executable,solution,cluster.executionpath,modelname,cluster.npocean); + end fclose(fid); Index: /issm/trunk/src/m/classes/geometry.js =================================================================== --- /issm/trunk/src/m/classes/geometry.js (revision 22821) +++ /issm/trunk/src/m/classes/geometry.js (revision 22822) @@ -37,5 +37,5 @@ checkfield(md,'fieldname','geometry.thickness','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices, 1],'>',0); for(var i=0;iMath.pow(10,9)){ + if (Math.abs(md.geometry.thickness[i]-md.geometry.surface[i]+md.geometry.base[i])>Math.pow(10,9)){ md = checkmessage(md,'equality thickness=surface-base violated'); break; Index: /issm/trunk/src/m/classes/slr.js =================================================================== --- /issm/trunk/src/m/classes/slr.js (revision 22821) +++ /issm/trunk/src/m/classes/slr.js (revision 22822) @@ -88,8 +88,8 @@ console.log(sprintf(' Sealevelrise solution parameters:')); - fielddisplay(this,'deltathickness','thickness change (main loading of the slr solution core [m]'); + fielddisplay(this,'deltathickness','thickness change: ice height equivalent [m]'); fielddisplay(this,'sealevel','current sea level (prior to computation) [m]'); fielddisplay(this,'reltol','sea level rise relative convergence criterion, (NaN: not applied)'); - fielddisplay(this,'abstol','sea level rise absolute convergence criterion, (default, NaN: not applied'); + fielddisplay(this,'abstol','sea level rise absolute convergence criterion, (default, NaN: not applied)'); fielddisplay(this,'maxiter','maximum number of nonlinear iterations'); fielddisplay(this,'love_h','load Love number for radial displacement'); @@ -105,6 +105,6 @@ fielddisplay(this,'elastic','elastic earth graviational potential perturbation'); fielddisplay(this,'rotation','rotational earth potential perturbation'); - fielddisplay(this,'ocean_area_scaling','correction for model representation of ocean area [default: No correction]'); - fielddisplay(this,'steric_rate','rate of steric ocean expansion (in mm/yr)'); + fielddisplay(this,'ocean_area_scaling','correction for model representation of ocean area (default: No correction)'); + fielddisplay(this,'steric_rate','rate of steric ocean expansion [mm/yr]'); fielddisplay(this,'degacc',"accuracy (default .01 deg) for numerical discretization of the Green's functions"); fielddisplay(this,'transitions','indices into parts of the mesh that will be icecaps'); Index: /issm/trunk/src/m/classes/slr.m =================================================================== --- /issm/trunk/src/m/classes/slr.m (revision 22821) +++ /issm/trunk/src/m/classes/slr.m (revision 22822) @@ -122,8 +122,8 @@ disp(sprintf(' slr parameters:')); - fielddisplay(self,'deltathickness','thickness change (main loading of the slr solution core [m]'); + fielddisplay(self,'deltathickness','thickness change: ice height equivalent [m]'); fielddisplay(self,'sealevel','current sea level (prior to computation) [m]'); fielddisplay(self,'reltol','sea level rise relative convergence criterion, (NaN: not applied)'); - fielddisplay(self,'abstol','sea level rise absolute convergence criterion, (default, NaN: not applied'); + fielddisplay(self,'abstol','sea level rise absolute convergence criterion, (default, NaN: not applied)'); fielddisplay(self,'maxiter','maximum number of nonlinear iterations'); fielddisplay(self,'love_h','load Love number for radial displacement'); @@ -139,6 +139,6 @@ fielddisplay(self,'elastic','elastic earth graviational potential perturbation'); fielddisplay(self,'rotation','earth rotational potential perturbation'); - fielddisplay(self,'ocean_area_scaling','correction for model representation of ocean area [default: No correction]'); - fielddisplay(self,'steric_rate','rate of steric ocean expansion (in mm/yr)'); + fielddisplay(self,'ocean_area_scaling','correction for model representation of ocean area (default: No correction)'); + fielddisplay(self,'steric_rate','rate of steric ocean expansion [mm/yr]'); fielddisplay(self,'degacc','accuracy (default .01 deg) for numerical discretization of the Green''s functions'); fielddisplay(self,'transitions','indices into parts of the mesh that will be icecaps'); Index: /issm/trunk/src/m/classes/slr.py =================================================================== --- /issm/trunk/src/m/classes/slr.py (revision 22821) +++ /issm/trunk/src/m/classes/slr.py (revision 22822) @@ -43,8 +43,8 @@ def __repr__(self): # {{{ string=' slr parameters:' - string="%s\n%s"%(string,fielddisplay(self,'deltathickness','thickness change (main loading of the slr solution core [m]')) + string="%s\n%s"%(string,fielddisplay(self,'deltathickness','thickness change: ice height equivalent [m]')) string="%s\n%s"%(string,fielddisplay(self,'sealevel','current sea level (prior to computation) [m]')) string="%s\n%s"%(string,fielddisplay(self,'reltol','sea level rise relative convergence criterion, (NaN: not applied)')) - string="%s\n%s"%(string,fielddisplay(self,'abstol','sea level rise absolute convergence criterion, (default, NaN: not applied')) + string="%s\n%s"%(string,fielddisplay(self,'abstol','sea level rise absolute convergence criterion, (default, NaN: not applied)')) string="%s\n%s"%(string,fielddisplay(self,'maxiter','maximum number of nonlinear iterations')) string="%s\n%s"%(string,fielddisplay(self,'love_h','load Love number for radial displacement')) @@ -60,6 +60,6 @@ string="%s\n%s"%(string,fielddisplay(self,'elastic','elastic earth graviational potential perturbation')) string="%s\n%s"%(string,fielddisplay(self,'rotation','earth rotational potential perturbation')) - string="%s\n%s"%(string,fielddisplay(self,'ocean_area_scaling','correction for model representation of ocean area [default: No correction]')) - string="%s\n%s"%(string,fielddisplay(self,'steric_rate','rate of steric ocean expansion (in mm/yr)')) + string="%s\n%s"%(string,fielddisplay(self,'ocean_area_scaling','correction for model representation of ocean area (default: No correction)')) + string="%s\n%s"%(string,fielddisplay(self,'steric_rate','rate of steric ocean expansion [mm/yr]')) string="%s\n%s"%(string,fielddisplay(self,'degacc','accuracy (default .01 deg) for numerical discretization of the Green''s functions')) string="%s\n%s"%(string,fielddisplay(self,'transitions','indices into parts of the mesh that will be icecaps')) Index: /issm/trunk/src/m/dev/issmversion.m =================================================================== --- /issm/trunk/src/m/dev/issmversion.m (revision 22821) +++ /issm/trunk/src/m/dev/issmversion.m (revision 22822) @@ -17,5 +17,5 @@ disp(['Compiled on ' IssmConfig('HOST_VENDOR') ' ' IssmConfig('HOST_OS') ' ' IssmConfig('HOST_ARCH') ' by ' IssmConfig('USER_NAME')]); disp([' ']); -disp(['Copyright (c) 2009-2017 California Institute of Technology']); +disp(['Copyright (c) 2009-2018 California Institute of Technology']); disp([' ']); disp([' to get started type: issmdoc']); Index: /issm/trunk/src/m/dev/issmversion.py =================================================================== --- /issm/trunk/src/m/dev/issmversion.py (revision 22821) +++ /issm/trunk/src/m/dev/issmversion.py (revision 22822) @@ -15,5 +15,5 @@ print ' ' print 'Build date: '+IssmConfig('PACKAGE_BUILD_DATE')[0] -print 'Copyright (c) 2009-2017 California Institute of Technology' +print 'Copyright (c) 2009-2018 California Institute of Technology' print ' ' print ' to get started type: issmdoc' Index: /issm/trunk/src/m/mesh/labelconnectedregions.m =================================================================== --- /issm/trunk/src/m/mesh/labelconnectedregions.m (revision 22822) +++ /issm/trunk/src/m/mesh/labelconnectedregions.m (revision 22822) @@ -0,0 +1,26 @@ +function labels = labelconnectedregions(md) +%LABELCONNECTEDREGIONS - label connected components of a mesh +% +% Usage: +% labels = labelconnectedregions(md) + +if size(md.mesh.elements,2)~=3, + error('not suppored yet (but easy to extend :)'); +end + +disp('Generate adjacency matrix'); +pairs = [ +md.mesh.elements(:,[1 2]) +md.mesh.elements(:,[2 1]) +md.mesh.elements(:,[2 3]) +md.mesh.elements(:,[3 2]) +md.mesh.elements(:,[3 1]) +md.mesh.elements(:,[1 3]) +]; +A = sparse(pairs(:,1), pairs(:,2), 1); + +disp('Construct graph'); +G = graph(A); + +disp('Label connected pieces'); +labels = conncomp(G); Index: /issm/trunk/src/m/mesh/planet/gmsh/gmshplanet.m =================================================================== --- /issm/trunk/src/m/mesh/planet/gmsh/gmshplanet.m (revision 22821) +++ /issm/trunk/src/m/mesh/planet/gmsh/gmshplanet.m (revision 22822) @@ -58,19 +58,19 @@ fprintf(fid,'Circle(12) = {6,1,2};\n'); fprintf(fid,'Line Loop(13) = {2,8,-10};\n'); - fprintf(fid,'Ruled Surface(14) = {13};\n'); + fprintf(fid,'Surface(14) = {13};\n'); fprintf(fid,'Line Loop(15) = {10,3,7};\n'); - fprintf(fid,'Ruled Surface(16) = {15};\n'); + fprintf(fid,'Surface(16) = {15};\n'); fprintf(fid,'Line Loop(17) = {-8,-9,1};\n'); - fprintf(fid,'Ruled Surface(18) = {17};\n'); + fprintf(fid,'Surface(18) = {17};\n'); fprintf(fid,'Line Loop(19) = {-11,-2,5};\n'); - fprintf(fid,'Ruled Surface(20) = {19};\n'); + fprintf(fid,'Surface(20) = {19};\n'); fprintf(fid,'Line Loop(21) = {-5,-12,-1};\n'); - fprintf(fid,'Ruled Surface(22) = {21};\n'); + fprintf(fid,'Surface(22) = {21};\n'); fprintf(fid,'Line Loop(23) = {-3,11,6};\n'); - fprintf(fid,'Ruled Surface(24) = {23};\n'); + fprintf(fid,'Surface(24) = {23};\n'); fprintf(fid,'Line Loop(25) = {-7,4,9};\n'); - fprintf(fid,'Ruled Surface(26) = {25};\n'); + fprintf(fid,'Surface(26) = {25};\n'); fprintf(fid,'Line Loop(27) = {-4,12,-6};\n'); - fprintf(fid,'Ruled Surface(28) = {27};\n'); + fprintf(fid,'Surface(28) = {27};\n'); fprintf(fid,'Surface Loop(29) = {28,26,16,14,20,24,22,18};\n'); fprintf(fid,'Volume(30) = {29};\n'); @@ -107,6 +107,5 @@ for i=paths if exist(i{1},'file'), - gmshpath = i{1} - break; + gmshpath = i{1}; end end Index: /issm/trunk/src/m/mesh/planet/gmsh/gmshplanet.py =================================================================== --- /issm/trunk/src/m/mesh/planet/gmsh/gmshplanet.py (revision 22821) +++ /issm/trunk/src/m/mesh/planet/gmsh/gmshplanet.py (revision 22822) @@ -64,19 +64,19 @@ fid.write('Circle(12) = {6,1,2};\n') fid.write('Line Loop(13) = {2,8,-10};\n') - fid.write('Ruled Surface(14) = {13};\n') + fid.write('Surface(14) = {13};\n') fid.write('Line Loop(15) = {10,3,7};\n') - fid.write('Ruled Surface(16) = {15};\n') + fid.write('Surface(16) = {15};\n') fid.write('Line Loop(17) = {-8,-9,1};\n') - fid.write('Ruled Surface(18) = {17};\n') + fid.write('Surface(18) = {17};\n') fid.write('Line Loop(19) = {-11,-2,5};\n') - fid.write('Ruled Surface(20) = {19};\n') + fid.write('Surface(20) = {19};\n') fid.write('Line Loop(21) = {-5,-12,-1};\n') - fid.write('Ruled Surface(22) = {21};\n') + fid.write('Surface(22) = {21};\n') fid.write('Line Loop(23) = {-3,11,6};\n') - fid.write('Ruled Surface(24) = {23};\n') + fid.write('Surface(24) = {23};\n') fid.write('Line Loop(25) = {-7,4,9};\n') - fid.write('Ruled Surface(26) = {25};\n') + fid.write('Surface(26) = {25};\n') fid.write('Line Loop(27) = {-4,12,-6};\n') - fid.write('Ruled Surface(28) = {27};\n') + fid.write('Surface(28) = {27};\n') fid.write('Surface Loop(29) = {28,26,16,14,20,24,22,18};\n') fid.write('Volume(30) = {29};\n') Index: /issm/trunk/test/MITgcm/code/cpl_issm.F =================================================================== --- /issm/trunk/test/MITgcm/code/cpl_issm.F (revision 22821) +++ /issm/trunk/test/MITgcm/code/cpl_issm.F (revision 22822) @@ -159,5 +159,5 @@ DO j=1,sNy DO i=1,sNx - local(i,j,bi,bj) = shelficeHeatFlux(i,j,bi,bj) + local(i,j,bi,bj)=shelficeFreshWaterFlux(i,j,bi,bj) ENDDO ENDDO @@ -175,5 +175,5 @@ ENDIF CALL BAR2( myThid ) - print*,'Done Sending shelficeHeatFlux array.' + print*,'Done Sending shelficeFreshWaterFlux array.' ENDIF