Index: /issm/oecreview/Archive/22755-22818/Date.tex =================================================================== --- /issm/oecreview/Archive/22755-22818/Date.tex (revision 22819) +++ /issm/oecreview/Archive/22755-22818/Date.tex (revision 22819) @@ -0,0 +1,1 @@ +May-31-2018 Index: /issm/oecreview/Archive/22755-22818/ISSM-22759-22760.diff =================================================================== --- /issm/oecreview/Archive/22755-22818/ISSM-22759-22760.diff (revision 22819) +++ /issm/oecreview/Archive/22755-22818/ISSM-22759-22760.diff (revision 22819) @@ -0,0 +1,52 @@ +Index: ../trunk-jpl/src/m/dev/issmversion.py +=================================================================== +--- ../trunk-jpl/src/m/dev/issmversion.py (revision 22759) ++++ ../trunk-jpl/src/m/dev/issmversion.py (revision 22760) +@@ -14,7 +14,7 @@ + print '(website: '+IssmConfig('PACKAGE_URL')[0]+' contact: '+IssmConfig('PACKAGE_BUGREPORT')[0]+')' + 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' + print ' ' +Index: ../trunk-jpl/src/m/dev/issmversion.m +=================================================================== +--- ../trunk-jpl/src/m/dev/issmversion.m (revision 22759) ++++ ../trunk-jpl/src/m/dev/issmversion.m (revision 22760) +@@ -16,7 +16,7 @@ + disp(['Build date: ' IssmConfig('PACKAGE_BUILD_DATE')]); + 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']); + disp([' ']); +Index: ../trunk-jpl/README +=================================================================== +--- ../trunk-jpl/README (revision 22759) ++++ ../trunk-jpl/README (revision 22760) +@@ -14,7 +14,7 @@ + ISSM is California Institute of Technology Copyright, Distributed under BSD Three Clause License. + + +-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: ../trunk-jpl/configure.ac +=================================================================== +--- ../trunk-jpl/configure.ac (revision 22759) ++++ ../trunk-jpl/configure.ac (revision 22760) +@@ -1,7 +1,7 @@ + # Process this file with autoconf to produce a configure script. + + #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 + m4_include([m4/issm_options.m4]) Index: /issm/oecreview/Archive/22755-22818/ISSM-22760-22761.diff =================================================================== --- /issm/oecreview/Archive/22755-22818/ISSM-22760-22761.diff (revision 22819) +++ /issm/oecreview/Archive/22755-22818/ISSM-22760-22761.diff (revision 22819) @@ -0,0 +1,25 @@ +Index: ../trunk-jpl/jenkins/jenkins.sh +=================================================================== +--- ../trunk-jpl/jenkins/jenkins.sh (revision 22760) ++++ ../trunk-jpl/jenkins/jenkins.sh (revision 22761) +@@ -295,9 +295,9 @@ + EOF + 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 & ++ $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 -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/src/m/dev; devpath; addpath $ISSM_DIR/nightlylog/; matlab_run$i" -logfile $ISSM_DIR/nightlylog/matlab_log$i.log & + fi + done + +@@ -469,7 +469,7 @@ + echo "disp('FAILURE');" >> $FILE + 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 + echo "finished: $(basename $dir)" >> $ISSM_DIR/nightlylog/matlab_log_examples.log Index: /issm/oecreview/Archive/22755-22818/ISSM-22761-22762.diff =================================================================== --- /issm/oecreview/Archive/22755-22818/ISSM-22761-22762.diff (revision 22819) +++ /issm/oecreview/Archive/22755-22818/ISSM-22761-22762.diff (revision 22819) @@ -0,0 +1,21 @@ +Index: ../trunk-jpl/jenkins/macosx_pine-island_static +=================================================================== +--- ../trunk-jpl/jenkins/macosx_pine-island_static (revision 22761) ++++ ../trunk-jpl/jenkins/macosx_pine-island_static (revision 22762) +@@ -21,6 +21,7 @@ + --with-mumps-dir=$ISSM_DIR/externalpackages/petsc/install \ + --with-metis-dir=$ISSM_DIR/externalpackages/petsc/install \ + --with-m1qn3-dir=$ISSM_DIR/externalpackages/m1qn3/install \ ++ --with-math77-dir=$ISSM_DIR/externalpackages/math77/install \ + --with-fortran-lib="/usr/local/gfortran/lib/libgfortran.a /usr/local/gfortran/lib/libquadmath.a /usr/local/gfortran/lib/gcc/x86_64-apple-darwin14/5.2.0/libgcc.a" \ + --with-numthreads=4' + +@@ -39,6 +40,8 @@ + m1qn3 install.sh + petsc install-3.6-macosx64-static.sh + triangle install-macosx64.sh ++ math77 install.sh ++ gmsh install.sh + shell2junit install.sh" + + #-----------------# Index: /issm/oecreview/Archive/22755-22818/ISSM-22769-22770.diff =================================================================== --- /issm/oecreview/Archive/22755-22818/ISSM-22769-22770.diff (revision 22819) +++ /issm/oecreview/Archive/22755-22818/ISSM-22769-22770.diff (revision 22819) @@ -0,0 +1,129 @@ +Index: ../trunk-jpl/src/m/plot/subplotmodel.m +=================================================================== +--- ../trunk-jpl/src/m/plot/subplotmodel.m (revision 22769) ++++ ../trunk-jpl/src/m/plot/subplotmodel.m (revision 22770) +@@ -29,7 +29,7 @@ + for j = 1:ncols + if(((i-1)*ncols+j)==num) + ha = axes('Units','normalized', ... +- 'Position',[xmin ymin width height]);%,'XTickLabel','','YTickLabel','','Visible','off'); ++ 'Position',[xmin ymin width height],'XTickLabel','','YTickLabel','','Visible','off'); + return + end + xmin = xmin+width+gap(2); +Index: ../trunk-jpl/src/m/classes/clusters/pfe.m +=================================================================== +--- ../trunk-jpl/src/m/classes/clusters/pfe.m (revision 22769) ++++ ../trunk-jpl/src/m/classes/clusters/pfe.m (revision 22770) +@@ -10,7 +10,7 @@ + % {{{ + name = 'pfe' + login = ''; +- modules = {'comp-intel/2018.0.128' 'mpi-sgi/mpt'}; ++ modules = {'comp-intel/2016.2.181' 'mpi-sgi/mpt'}; + numnodes = 20; + cpuspernode = 8; + port = 1025; +@@ -162,12 +162,10 @@ + fprintf(fid,'#PBS -o %s.outlog \n',[cluster.executionpath '/' dirname '/' modelname]); + fprintf(fid,'#PBS -e %s.errlog \n\n',[cluster.executionpath '/' dirname '/' modelname]); + fprintf(fid,'. /usr/share/modules/init/bash\n\n'); +- for i=1:numel(cluster.modules), +- fprintf(fid,['module load ' cluster.modules{i} '\n']); +- end ++ fprintf(fid,'module load comp-intel/2016.2.181\n'); ++ fprintf(fid,'module load mpi-sgi/mpt\n'); + fprintf(fid,'export PATH="$PATH:."\n\n'); + fprintf(fid,'export MPI_GROUP_MAX=64\n\n'); +- fprintf(fid,'export MKL_NUM_THREADS=2\n\n'); + fprintf(fid,'export ISSM_DIR="%s/../"\n',cluster.codepath); %FIXME + fprintf(fid,'source $ISSM_DIR/etc/environment.sh\n'); %FIXME + fprintf(fid,'cd %s/%s/\n\n',cluster.executionpath,dirname); +@@ -227,9 +225,8 @@ + fprintf(fid,'#PBS -o %s.outlog \n',[cluster.executionpath '/' dirname '/' modelname]); + fprintf(fid,'#PBS -e %s.errlog \n\n',[cluster.executionpath '/' dirname '/' modelname]); + fprintf(fid,'. /usr/share/modules/init/bash\n\n'); +- for i=1:numel(cluster.modules), +- fprintf(fid,['module load ' cluster.modules{i} '\n']); +- end ++ fprintf(fid,'module load comp-intel/2016.2.181\n'); ++ fprintf(fid,'module load mpi-sgi/mpt\n'); + fprintf(fid,'export PATH="$PATH:."\n\n'); + fprintf(fid,'export MPI_GROUP_MAX=64\n\n'); + fprintf(fid,'export ISSM_DIR="%s/../"\n',cluster.codepath); %FIXME +@@ -355,9 +352,6 @@ + fprintf(fid,'#PBS -o %s.outlog \n',modelname); + fprintf(fid,'#PBS -e %s.errlog \n\n',modelname); + fprintf(fid,'. /usr/share/modules/init/bash\n\n'); +- %for i=1:numel(cluster.modules), +- % fprintf(fid,['module load ' cluster.modules{i} '\n']); +- %end + fprintf(fid,'module load comp-intel/2016.2.181\n'); + fprintf(fid,'module load netcdf/4.4.1.1_mpt\n'); + fprintf(fid,'module load mpi-sgi/mpt.2.15r20\n'); +Index: ../trunk-jpl/src/m/classes/clusters/generic.m +=================================================================== +--- ../trunk-jpl/src/m/classes/clusters/generic.m (revision 22769) ++++ ../trunk-jpl/src/m/classes/clusters/generic.m (revision 22770) +@@ -110,10 +110,10 @@ + fprintf(fid,'\n gprof %s/issm.exe gmon.out > %s.performance',cluster.codepath,modelname); + else + %Add --gen-suppressions=all to get suppression lines +- %fprintf(fid,'LD_PRELOAD=%s \\\n',cluster.valgrindlib); ++ fprintf(fid,'LD_PRELOAD=%s \\\n',cluster.valgrindlib); + if ismac, + if IssmConfig('_HAVE_MPI_'), +- fprintf(fid,'mpiexec -np %i %s --leak-check=full --gen-suppressions=all --error-limit=no --dsymutil=yes --suppressions=%s %s/%s %s %s %s 2> %s.errlog >%s.outlog ',... ++ fprintf(fid,'mpiexec -np %i %s --leak-check=full --error-limit=no --dsymutil=yes --suppressions=%s %s/%s %s %s %s 2> %s.errlog >%s.outlog ',... + cluster.np,cluster.valgrind,cluster.valgrindsup,cluster.codepath,executable,solution,[cluster.executionpath '/' dirname], modelname,modelname,modelname); + else + fprintf(fid,'%s --leak-check=full --dsymutil=yes --error-limit=no --suppressions=%s %s/%s %s %s %s 2> %s.errlog >%s.outlog ',... +@@ -121,7 +121,7 @@ + end + else + if IssmConfig('_HAVE_MPI_'), +- fprintf(fid,'mpiexec -np %i %s --leak-check=full --gen-suppressions=all --error-limit=no --suppressions=%s %s/%s %s %s %s 2> %s.errlog >%s.outlog ',... ++ fprintf(fid,'mpiexec -np %i %s --leak-check=full --error-limit=no --suppressions=%s %s/%s %s %s %s 2> %s.errlog >%s.outlog ',... + cluster.np,cluster.valgrind,cluster.valgrindsup,cluster.codepath,executable,solution,[cluster.executionpath '/' dirname],modelname,modelname,modelname); + else + fprintf(fid,'%s --leak-check=full --error-limit=no --suppressions=%s %s/%s %s %s %s 2> %s.errlog >%s.outlog ',... +Index: ../trunk-jpl/src/m/solvers/mumpsoptions.m +=================================================================== +--- ../trunk-jpl/src/m/solvers/mumpsoptions.m (revision 22769) ++++ ../trunk-jpl/src/m/solvers/mumpsoptions.m (revision 22770) +@@ -26,11 +26,9 @@ + mumps.pc_type=getfieldvalue(options,'pc_type','lu'); + mumps.pc_factor_mat_solver_package=getfieldvalue(options,'pc_factor_mat_solver_package','mumps'); + mumps.mat_mumps_icntl_14=getfieldvalue(options,'mat_mumps_icntl_14',120); +- +- %Seems like this one is not needed anymore + mumps.pc_factor_shift_positive_definite=getfieldvalue(options,'pc_factor_shift_positive_definite','true'); + + %These 2 lines make raijin break (ptwgts error during solver with PETSc 3.3) +- %mumps.mat_mumps_icntl_28=2; %1=serial, 2=parallel +- %mumps.mat_mumps_icntl_29=2; %parallel ordering 1 = ptscotch, 2 = parmetis ++ mumps.mat_mumps_icntl_28=2; %1=serial, 2=parallel ++ mumps.mat_mumps_icntl_29=2; %parallel ordering 1 = ptscotch, 2 = parmetis + end +Index: ../trunk-jpl/src/m/inversions/supportedcontrols.m +=================================================================== +--- ../trunk-jpl/src/m/inversions/supportedcontrols.m (revision 22769) ++++ ../trunk-jpl/src/m/inversions/supportedcontrols.m (revision 22770) +@@ -1,7 +1,6 @@ + function list = supportedcontrols(), + + list = {... +- 'BalancethicknessSpcthickness',... + 'BalancethicknessThickeningRate',... + 'FrictionCoefficient',... + 'FrictionAs',... +Index: ../trunk-jpl/src/m/inversions/marshallcostfunctions.m +=================================================================== +--- ../trunk-jpl/src/m/inversions/marshallcostfunctions.m (revision 22769) ++++ ../trunk-jpl/src/m/inversions/marshallcostfunctions.m (revision 22770) +@@ -14,5 +14,3 @@ + pos=find(cost_functions==507); data(pos) = {'RheologyBAbsGradient'}; + pos=find(cost_functions==510); data(pos) = {'ThicknessPositive'}; + pos=find(cost_functions==601); data(pos) = {'SurfaceAbsMisfit'}; +- pos=find(cost_functions==602); data(pos) = {'OmegaAbsGradient'}; +- pos=find(cost_functions==603); data(pos) = {'EtaDiff'}; Index: /issm/oecreview/Archive/22755-22818/ISSM-22777-22778.diff =================================================================== --- /issm/oecreview/Archive/22755-22818/ISSM-22777-22778.diff (revision 22819) +++ /issm/oecreview/Archive/22755-22818/ISSM-22777-22778.diff (revision 22819) @@ -0,0 +1,270 @@ +Index: ../trunk-jpl/src/c/classes/Vertex.cpp +=================================================================== +--- ../trunk-jpl/src/c/classes/Vertex.cpp (revision 22777) ++++ ../trunk-jpl/src/c/classes/Vertex.cpp (revision 22778) +@@ -31,6 +31,9 @@ + this->z = iomodel->Data("md.mesh.z")[i]; + this->domaintype = iomodel->domaintype; + ++ this->latitute = iomodel->Data("md.mesh.lat")[i]; ++ this->longitude = iomodel->Data("md.mesh.long")[i]; ++ + switch(iomodel->domaintype){ + case Domain3DEnum: + _assert_(iomodel->Data("md.geometry.base") && iomodel->Data("md.geometry.thickness")); +@@ -195,23 +198,6 @@ + /*}}}*/ + 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){/*{{{*/ + + /*If we are not a clone, don't update, we already have pids: */ +Index: ../trunk-jpl/src/c/classes/Vertices.h +=================================================================== +--- ../trunk-jpl/src/c/classes/Vertices.h (revision 22777) ++++ ../trunk-jpl/src/c/classes/Vertices.h (revision 22778) +@@ -24,7 +24,7 @@ + void FlagClones(int numberofnodes); + int NumberOfVertices(void); + void Ranks(int* ranks); +- IssmDouble* ToXYZ(void); ++ void LatLonList(IssmDouble** lat,IssmDouble** lon); + }; + + #endif //ifndef _VERTICES_H_ +Index: ../trunk-jpl/src/c/classes/Vertex.h +=================================================================== +--- ../trunk-jpl/src/c/classes/Vertex.h (revision 22777) ++++ ../trunk-jpl/src/c/classes/Vertex.h (revision 22778) +@@ -61,7 +61,6 @@ + void SetClone(int* minranks); + 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); + void VertexCoordinates(Vector* vx,Vector* vy,Vector* vz,bool spherical=false); +Index: ../trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp +=================================================================== +--- ../trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp (revision 22777) ++++ ../trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp (revision 22778) +@@ -15,6 +15,7 @@ + bool control_analysis; + bool dakota_analysis; + bool adolc_analysis; ++ bool isoceancoupling; + int nnat,dummy; + int* nature=NULL; + +@@ -23,6 +24,7 @@ + iomodel->FindConstant(&dakota_analysis,"md.qmu.isdakota"); + 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? : */ + _assert_(elements->Size()==0); +@@ -231,6 +233,7 @@ + iomodel->FetchData(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->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); + +@@ -241,4 +244,5 @@ + /*Free data: */ + 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: ../trunk-jpl/src/c/cores/transient_core.cpp +=================================================================== +--- ../trunk-jpl/src/c/cores/transient_core.cpp (revision 22777) ++++ ../trunk-jpl/src/c/cores/transient_core.cpp (revision 22778) +@@ -92,6 +92,7 @@ + if(!parcom)_error_("TransferForcing error message: could not find ToMITgcmCommEnum communicator"); + 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); + ISSM_MPI_Recv(&oceangridnysize,1,ISSM_MPI_INT,0,10001004,tomitgcmcomm,&status); +@@ -101,7 +102,10 @@ + ISSM_MPI_Recv(oceangridx,oceangridnxsize*oceangridnysize,ISSM_MPI_DOUBLE,0,10001005,tomitgcmcomm,&status); + 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); + icebase = xNew(oceangridnxsize*oceangridnysize); +@@ -161,12 +165,14 @@ + save_results=false; + femmodel->parameters->SetParam(save_results,SaveResultsEnum); + +- if(isoceancoupling){ {{{ ++ if(isoceancoupling){ /*{{{*/ + if(VerboseSolution()) _printf0_(" ocean coupling: exchanging information\n"); + int my_rank; +- int oceangridnxsize,oceangridnysize; ++ int oceangridnxsize,oceangridnysize,ngrids_ocean; + IssmDouble *oceanmelt; + IssmDouble *icebase; ++ IssmDouble *oceangridx; ++ IssmDouble *oceangridy; + IssmDouble oceantime; + ISSM_MPI_Comm tomitgcmcomm; + ISSM_MPI_Status status; +@@ -179,22 +185,59 @@ + ISSM_MPI_Send(&time,1,ISSM_MPI_DOUBLE,0,10001001,tomitgcmcomm); + 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"); ++ ++ /*Recover inputs needed*/ + 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); ++ femmodel->parameters->FindParam(&oceangridx,&ngrids_ocean,OceanGridXEnum); ++ femmodel->parameters->FindParam(&oceangridy,&ngrids_ocean,OceanGridYEnum); ++ ++ /*Exchange information*/ ++ oceanmelt = xNew(ngrids_ocean); ++ ISSM_MPI_Recv(oceanmelt,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001007,tomitgcmcomm,&status); ++ icebase = xNew(ngrids_ocean); + for(int i=0;ivertices->NumberOfVertices(); ++ lat_ice= xNew(ngrids_ice); ++ lon_ice= xNew(ngrids_ice); ++ melt_mesh= xNew(ngrids_ice); ++ femmodel->GetMesh(femmodel->vertices,femmodel->elements,&x_ice,&y_ice,&z_ice,&index_ice); ++ BamgTriangulatex(&index_ocean,&nels_ocean,oceangridx,oceangridy,ngrids_ocean); ++ femmodel->vertices->LatLonList(&lat_ice,&lon_ice); ++ InterpFromMeshToMesh2dx(&melt_mesh,index_ocean,oceangridx,oceangridy,ngrids_ocean,nels_ocean, ++ oceanmelt,ngrids_ocean,1, ++ lat_ice,lon_ice,ngrids_ice,NULL); /*FIXME: may be wrong: y=lat x=lon*/ ++ InputUpdateFromVectorx(femmodel,melt_mesh,BasalforcingsFloatingiceMeltingRateEnum,VertexSIdEnum); ++ xDelete(lat_ice); ++ xDelete(lon_ice); ++ xDelete(x_ice); ++ xDelete(y_ice); ++ xDelete(z_ice); ++ xDelete(melt_mesh); ++ xDelete(index_ice); ++ xDelete(oceangridx); ++ xDelete(oceangridy); + xDelete(oceanmelt); + xDelete(icebase); + } + } +- }}} ++ /*}}}*/ + + if(isthermal && domaintype==Domain3DEnum){ + if(issmb){ +Index: ../trunk-jpl/src/c/classes/Vertices.cpp +=================================================================== +--- ../trunk-jpl/src/c/classes/Vertices.cpp (revision 22777) ++++ ../trunk-jpl/src/c/classes/Vertices.cpp (revision 22778) +@@ -179,45 +179,41 @@ + } + } + /*}}}*/ +-IssmDouble* Vertices::ToXYZ(void){/*{{{*/ ++void Vertices::LatLonList(IssmDouble** plat,IssmDouble** plon){/*{{{*/ + +- /*intermediary: */ +- int i; +- int my_rank; +- int num_vertices; +- + /*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(); ++ int num_vertices=this->NumberOfVertices(); + +- /*Now, allocate matrix to hold all the vertices x,y and z values: */ +- xyz= new Matrix(num_vertices,3); ++ /*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: */ ++ for(int i=0;iSize();i++){ + Vertex* vertex=xDynamicCast(this->GetObjectByOffset(i)); +- vertex->ToXYZ(xyz); ++ 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(); + ++ printf("lat_serial %g %g %g \n",lat_serial[0],lat_serial[2], lat_serial[799]); + /*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/oecreview/Archive/22755-22818/ISSM-22778-22779.diff =================================================================== --- /issm/oecreview/Archive/22755-22818/ISSM-22778-22779.diff (revision 22819) +++ /issm/oecreview/Archive/22755-22818/ISSM-22778-22779.diff (revision 22819) @@ -0,0 +1,14 @@ +Index: ../trunk-jpl/src/c/classes/Vertices.cpp +=================================================================== +--- ../trunk-jpl/src/c/classes/Vertices.cpp (revision 22778) ++++ ../trunk-jpl/src/c/classes/Vertices.cpp (revision 22779) +@@ -213,7 +213,7 @@ + /*free ressources: */ + *plat = lat_serial; + *plon = lon_serial; +- //delete lat; +- //delete lon; ++ delete lat; ++ delete lon; + } + /*}}}*/ Index: /issm/oecreview/Archive/22755-22818/ISSM-22779-22780.diff =================================================================== --- /issm/oecreview/Archive/22755-22818/ISSM-22779-22780.diff (revision 22819) +++ /issm/oecreview/Archive/22755-22818/ISSM-22779-22780.diff (revision 22819) @@ -0,0 +1,12 @@ +Index: ../trunk-jpl/src/c/classes/Vertices.cpp +=================================================================== +--- ../trunk-jpl/src/c/classes/Vertices.cpp (revision 22779) ++++ ../trunk-jpl/src/c/classes/Vertices.cpp (revision 22780) +@@ -209,7 +209,6 @@ + IssmDouble* lat_serial=lat->ToMPISerial(); + IssmDouble* lon_serial=lon->ToMPISerial(); + +- printf("lat_serial %g %g %g \n",lat_serial[0],lat_serial[2], lat_serial[799]); + /*free ressources: */ + *plat = lat_serial; + *plon = lon_serial; Index: /issm/oecreview/Archive/22755-22818/ISSM-22786-22787.diff =================================================================== --- /issm/oecreview/Archive/22755-22818/ISSM-22786-22787.diff (revision 22819) +++ /issm/oecreview/Archive/22755-22818/ISSM-22786-22787.diff (revision 22819) @@ -0,0 +1,86 @@ +Index: ../trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp +=================================================================== +--- ../trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 22786) ++++ ../trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 22787) +@@ -178,10 +178,6 @@ + _error_("Basal forcing model "<FetchDataToInput(elements,"md.mesh.lat",MeshLatEnum); +- iomodel->FetchDataToInput(elements,"md.mesh.long",MeshLongEnum); +- } + if(!issmb){ + iomodel->FetchDataToInput(elements,"md.smb.mass_balance",SmbMassBalanceEnum); + } +Index: ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h +=================================================================== +--- ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 22786) ++++ ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 22787) +@@ -854,8 +854,6 @@ + MelangeEnum, + MeltingAnalysisEnum, + MeshElementsEnum, +- MeshLatEnum, +- MeshLongEnum, + MeshXEnum, + MeshYEnum, + MINIcondensedEnum, +Index: ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp +=================================================================== +--- ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 22786) ++++ ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 22787) +@@ -858,8 +858,6 @@ + case MelangeEnum : return "Melange"; + case MeltingAnalysisEnum : return "MeltingAnalysis"; + case MeshElementsEnum : return "MeshElements"; +- case MeshLatEnum : return "MeshLat"; +- case MeshLongEnum : return "MeshLong"; + case MeshXEnum : return "MeshX"; + case MeshYEnum : return "MeshY"; + case MINIcondensedEnum : return "MINIcondensed"; +Index: ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp +=================================================================== +--- ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 22786) ++++ ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 22787) +@@ -879,8 +879,6 @@ + if(stage==8){ + 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; + else if (strcmp(name,"MINIcondensed")==0) return MINIcondensedEnum; +@@ -997,12 +995,12 @@ + else if (strcmp(name,"Outputdefinition8")==0) return Outputdefinition8Enum; + 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; + else if (strcmp(name,"Outputdefinition97")==0) return Outputdefinition97Enum; +@@ -1120,12 +1118,12 @@ + else if (strcmp(name,"TransientInput")==0) return TransientInputEnum; + 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; + else if (strcmp(name,"VertexPId")==0) return VertexPIdEnum; Index: /issm/oecreview/Archive/22755-22818/ISSM-22793-22794.diff =================================================================== --- /issm/oecreview/Archive/22755-22818/ISSM-22793-22794.diff (revision 22819) +++ /issm/oecreview/Archive/22755-22818/ISSM-22793-22794.diff (revision 22819) @@ -0,0 +1,71 @@ +Index: ../trunk-jpl/src/c/cores/controladm1qn3_core.cpp +=================================================================== +--- ../trunk-jpl/src/c/cores/controladm1qn3_core.cpp (revision 22793) ++++ ../trunk-jpl/src/c/cores/controladm1qn3_core.cpp (revision 22794) +@@ -105,6 +105,7 @@ + if(X[index]>XU[index]) X[index]=XU[index]; + if(X[index](scaling_factors[c]); + Gnorm += G[index]*G[index]; + } ++ N_add+=N[c]; + } + Gnorm = sqrt(Gnorm); + +Index: ../trunk-jpl/src/c/classes/Elements/Element.cpp +=================================================================== +--- ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 22793) ++++ ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 22794) +@@ -1661,7 +1661,7 @@ + IssmDouble* times = xNew(N); + for(int t=0;tObjectEnum()){ + 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; +Index: ../trunk-jpl/src/c/classes/Inputs/ControlInput.cpp +=================================================================== +--- ../trunk-jpl/src/c/classes/Inputs/ControlInput.cpp (revision 22793) ++++ ../trunk-jpl/src/c/classes/Inputs/ControlInput.cpp (revision 22794) +@@ -329,7 +329,7 @@ + if(input->ObjectEnum()!=TransientInputEnum)_error_("cannot have timeoffset argument if not TransientInput Control"); + 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); + else { +Index: ../trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp +=================================================================== +--- ../trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp (revision 22793) ++++ ../trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp (revision 22794) +@@ -142,6 +142,11 @@ + 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)); ++ } + else{ + _error_("Field \""<',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/oecreview/Archive/22755-22818/ISSM-22795-22796.diff =================================================================== --- /issm/oecreview/Archive/22755-22818/ISSM-22795-22796.diff (revision 22819) +++ /issm/oecreview/Archive/22755-22818/ISSM-22795-22796.diff (revision 22819) @@ -0,0 +1,309 @@ +Index: ../trunk-jpl/examples/EsaGRACE/runme.m +=================================================================== +--- ../trunk-jpl/examples/EsaGRACE/runme.m (nonexistent) ++++ ../trunk-jpl/examples/EsaGRACE/runme.m (revision 22796) +@@ -0,0 +1,171 @@ ++ ++clear all; ++steps=[5]; % [1:6]; ++ ++if any(steps==0) % Download GRACE land_mass data {{{ ++ disp(' Step 0: Download GRACE land_mass data'); ++ ++ % Connect to ftp server and download ++ f = ftp('podaac-ftp.jpl.nasa.gov'); ++ %dir(f) ++ cd(f,'allData/tellus/L3/land_mass/RL05/netcdf/'); ++ mget(f,'GRCTellus.JPL.200204_201701.LND.RL05_1.DSTvSCS1411.nc'); ++ ++ % display the content of the netcdf file. ++ %ncdisp('GRCTellus.JPL.200204_201701.LND.RL05_1.DSTvSCS1411.nc') ++ ++end % }}} ++ ++if any(steps==1) % Global mesh creation {{{ ++ disp(' Step 1: Global mesh creation'); ++ ++ resolution=300; % km. uniform meshing... ++ radius = 6.371012*10^3; % mean radius of Earth, km ++ ++ %mesh earth: ++ md=model; ++ md.mask=maskpsl(); % use maskpsl class (instead of mask) to store the ocean function as a ocean_levelset ++ md.mesh=gmshplanet('radius',radius,'resolution',resolution); ++ ++ %figure out mask: ++ 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'); ++ %export_fig('Fig1.pdf'); ++ ++end % }}} ++ ++if any(steps==2) % Define loads {{{ ++ disp(' Step 2: Define loads in meters of ice height equivalent'); ++ md = loadmodel('./Models/EsaGRACE.Mesh'); ++ ++ year_month = 2007+15/365; ++ time_min=year_month; ++ time_max=year_month; ++ ++ % map GRACE water load on to the mesh for the seleted month(s) ++ water_load = grace(md.mesh.elements,md.mesh.lat,md.mesh.long,time_min,time_max); ++ ++ 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]'); ++ %export_fig('Fig2.pdf'); ++ ++end % }}} ++ ++if any(steps==3) % Parameterization {{{ ++ disp(' Step 3: Parameterization'); ++ md = loadmodel('./Models/EsaGRACE.Loads'); ++ ++ % Love numbers and reference frame: CF or CM (choose one!) ++ nlove=10001; % up to 10,000 degree ++ 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)=[]; ++ ++ % Mask: for computational efficiency only those elements that have loads are convolved! ++ md.mask.groundedice_levelset = ones(md.mesh.numberofvertices,1); % 1 = ice is grounnded ++ md.mask.ice_levelset = ones(md.mesh.numberofvertices,1); ++ pos=find(md.esa.deltathickness~=0); ++ md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; % -1 = ice loads ++ md.mask.land_levelset = 1-md.mask.ocean_levelset; ++ ++ %% IGNORE BUT DO NOT DELETE %% {{{ ++ % Geometry: Important only when you want to couple with Ice Flow Model ++ 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; ++ % Materials: ++ 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); ++ % Miscellaneous: ++ md.miscellaneous.name='EsaGRACE'; ++ %% IGNORE BUT DO NOT DELETE %% }}} ++ ++ save ./Models/EsaGRACE.Parameterization md; ++ ++end % }}} ++ ++if any(steps==4) % Solve {{{ ++ disp(' Step 4: Solve Esa solver'); ++ md = loadmodel('./Models/EsaGRACE.Parameterization'); ++ ++ % Request outputs ++ md.esa.requested_outputs = {'EsaUmotion','EsaNmotion','EsaEmotion'}; ++ ++ % Cluster info ++ md.cluster=generic('name',oshostname(),'np',3); ++ md.verbose=verbose('111111111'); ++ ++ % Solve ++ md=solve(md,'Esa'); ++ ++ save ./Models/EsaGRACE.Solution md; ++ ++end % }}} ++ ++if any(steps==5) % Plot solutions {{{ ++ disp(' Step 5: Plot solutions'); ++ md = loadmodel('./Models/EsaGRACE.Solution'); ++ ++ % solutions. ++ ur = md.results.EsaSolution.EsaUmotion*1000; % [mm] ++ un = md.results.EsaSolution.EsaNmotion*1000; % [mm] ++ ue = md.results.EsaSolution.EsaEmotion*1000; % [mm] ++ ++ res = 1.0; % degree ++ ++ % Make a grid of lats and lons, based on the min and max of the original vectors ++ [lat_grid, lon_grid] = meshgrid(linspace(-90,90,180/res), linspace(-180,180,360/res)); ++ sol_grid = zeros(size(lat_grid)); ++ ++ sol = ue; ++ ++ % Make a interpolation object ++ F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol); ++ F.Method = 'linear'; ++ F.ExtrapolationMethod = 'linear'; ++ ++ % Do the interpolation to get gridded solutions... ++ sol_grid = F(lat_grid, lon_grid); ++ sol_grid(isnan(sol_grid))=0; ++ ++ % set polar unphysical 0s to Nan ++ 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; ++ %caxis([-0.3 0.3]) ++ hold on ++ plot(long,lat,'k'); ++ %geoshow(lat,long,'DisplayType','polygon','FaceColor','white'); ++ hold off; ++ c1=colorbar; ++ colormap(jet); ++ xlim([-180 180]); ++ ylim([-90 90]); ++ grid on; ++ title(['Average change in relative sea-level [mm/yr]']); ++ set(gcf,'color','w'); ++ ++ %export_fig('Fig5.pdf'); ++ ++end % }}} ++ ++if any(steps==6) % {{{ Transient ++ disp(' Step 6: Transient run'); ++ ++end % }}} ++ +Index: ../trunk-jpl/examples/EsaGRACE/grace.m +=================================================================== +--- ../trunk-jpl/examples/EsaGRACE/grace.m (nonexistent) ++++ ../trunk-jpl/examples/EsaGRACE/grace.m (revision 22796) +@@ -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; ++ + +Property changes on: ../trunk-jpl/examples/EsaGRACE/grace.m +___________________________________________________________________ +Added: svn:executable +## -0,0 +1 ## ++* +\ No newline at end of property Index: /issm/oecreview/Archive/22755-22818/ISSM-22796-22797.diff =================================================================== --- /issm/oecreview/Archive/22755-22818/ISSM-22796-22797.diff (revision 22819) +++ /issm/oecreview/Archive/22755-22818/ISSM-22796-22797.diff (revision 22819) @@ -0,0 +1,22 @@ +Index: ../trunk-jpl/test/MITgcm/code/cpl_issm.F +=================================================================== +--- ../trunk-jpl/test/MITgcm/code/cpl_issm.F (revision 22796) ++++ ../trunk-jpl/test/MITgcm/code/cpl_issm.F (revision 22797) +@@ -158,7 +158,7 @@ + DO bi=1,nSx + 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 + ENDDO +@@ -174,7 +174,7 @@ + _END_MASTER( myThid ) + ENDIF + CALL BAR2( myThid ) +- print*,'Done Sending shelficeHeatFlux array.' ++ print*,'Done Sending shelficeFreshWaterFlux array.' + + ENDIF + C End recurring step C3. Index: /issm/oecreview/Archive/22755-22818/ISSM-22797-22798.diff =================================================================== --- /issm/oecreview/Archive/22755-22818/ISSM-22797-22798.diff (revision 22819) +++ /issm/oecreview/Archive/22755-22818/ISSM-22797-22798.diff (revision 22819) @@ -0,0 +1,253 @@ +Index: ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h +=================================================================== +--- ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 22797) ++++ ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 22798) +@@ -199,6 +199,7 @@ + MasstransportPenaltyFactorEnum, + MasstransportRequestedOutputsEnum, + MasstransportStabilizationEnum, ++ MaterialsRhoIceEnum, + MeltingOffsetEnum, + MeshAverageVertexConnectivityEnum, + MeshElementtypeEnum, +@@ -833,7 +834,6 @@ + MaterialsMuWaterEnum, + MaterialsRheologyLawEnum, + MaterialsRhoFreshwaterEnum, +- MaterialsRhoIceEnum, + MaterialsRhoSeawaterEnum, + MaterialsTemperateiceconductivityEnum, + MaterialsThermalconductivityEnum, +Index: ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp +=================================================================== +--- ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 22797) ++++ ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 22798) +@@ -207,6 +207,7 @@ + case MasstransportPenaltyFactorEnum : return "MasstransportPenaltyFactor"; + case MasstransportRequestedOutputsEnum : return "MasstransportRequestedOutputs"; + case MasstransportStabilizationEnum : return "MasstransportStabilization"; ++ case MaterialsRhoIceEnum : return "MaterialsRhoIce"; + case MeltingOffsetEnum : return "MeltingOffset"; + case MeshAverageVertexConnectivityEnum : return "MeshAverageVertexConnectivity"; + case MeshElementtypeEnum : return "MeshElementtype"; +@@ -837,7 +838,6 @@ + case MaterialsMuWaterEnum : return "MaterialsMuWater"; + case MaterialsRheologyLawEnum : return "MaterialsRheologyLaw"; + case MaterialsRhoFreshwaterEnum : return "MaterialsRhoFreshwater"; +- case MaterialsRhoIceEnum : return "MaterialsRhoIce"; + case MaterialsRhoSeawaterEnum : return "MaterialsRhoSeawater"; + case MaterialsTemperateiceconductivityEnum : return "MaterialsTemperateiceconductivity"; + case MaterialsThermalconductivityEnum : return "MaterialsThermalconductivity"; +Index: ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp +=================================================================== +--- ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 22797) ++++ ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 22798) +@@ -210,6 +210,7 @@ + else if (strcmp(name,"MasstransportPenaltyFactor")==0) return MasstransportPenaltyFactorEnum; + 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; + else if (strcmp(name,"MeshElementtype")==0) return MeshElementtypeEnum; +@@ -258,11 +259,11 @@ + else if (strcmp(name,"SealevelriseTidalLoveH")==0) return SealevelriseTidalLoveHEnum; + 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; + else if (strcmp(name,"SettingsResultsOnNodes")==0) return SettingsResultsOnNodesEnum; +@@ -381,11 +382,11 @@ + else if (strcmp(name,"BasalforcingsGroundediceMeltingRate")==0) return BasalforcingsGroundediceMeltingRateEnum; + 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; + else if (strcmp(name,"Bed")==0) return BedEnum; +@@ -504,11 +505,11 @@ + else if (strcmp(name,"RheologyBbarAbsGradient")==0) return RheologyBbarAbsGradientEnum; + 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; + else if (strcmp(name,"SmbAccumulation")==0) return SmbAccumulationEnum; +@@ -627,11 +628,11 @@ + else if (strcmp(name,"WaterfractionDrainage")==0) return WaterfractionDrainageEnum; + 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; + else if (strcmp(name,"AdaptiveTimestepping")==0) return AdaptiveTimesteppingEnum; +@@ -750,11 +751,11 @@ + else if (strcmp(name,"GaussSeg")==0) return GaussSegEnum; + 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; + else if (strcmp(name,"GiaSolution")==0) return GiaSolutionEnum; +@@ -855,7 +856,6 @@ + else if (strcmp(name,"MaterialsMuWater")==0) return MaterialsMuWaterEnum; + 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; + else if (strcmp(name,"MaterialsThermalconductivity")==0) return MaterialsThermalconductivityEnum; +Index: ../trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp +=================================================================== +--- ../trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 22797) ++++ ../trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 22798) +@@ -104,6 +104,7 @@ + parameters->AddObject(iomodel->CopyConstantObject("md.transient.isslr",TransientIsslrEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.transient.iscoupler",TransientIscouplerEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.transient.isoceancoupling",TransientIsoceancouplingEnum)); ++ parameters->AddObject(iomodel->CopyConstantObject("md.materials.rho_ice",MaterialsRhoIceEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.transient.amr_frequency",TransientAmrFrequencyEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.gia.cross_section_shape",GiaCrossSectionShapeEnum)); + +Index: ../trunk-jpl/src/c/cores/transient_core.cpp +=================================================================== +--- ../trunk-jpl/src/c/cores/transient_core.cpp (revision 22797) ++++ ../trunk-jpl/src/c/cores/transient_core.cpp (revision 22798) +@@ -169,12 +169,6 @@ + #ifndef _HAVE_ADOLC_ + if(VerboseSolution()) _printf0_(" ocean coupling: exchanging information\n"); + int my_rank; +- int oceangridnxsize,oceangridnysize,ngrids_ocean; +- IssmDouble *oceanmelt; +- IssmDouble *icebase; +- IssmDouble *oceangridx; +- IssmDouble *oceangridy; +- IssmDouble oceantime; + ISSM_MPI_Comm tomitgcmcomm; + ISSM_MPI_Status status; + +@@ -183,48 +177,63 @@ + if(!parcom)_error_("TransferForcing error message: could not find ToMITgcmCommEnum communicator"); + tomitgcmcomm=parcom->GetParameterValue(); + if(my_rank==0){ +- ISSM_MPI_Send(&time,1,ISSM_MPI_DOUBLE,0,10001001,tomitgcmcomm); +- 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"); +- +- /*Recover inputs needed*/ +- femmodel->parameters->FindParam(&oceangridnxsize,OceanGridNxEnum); +- femmodel->parameters->FindParam(&oceangridnysize,OceanGridNyEnum); +- femmodel->parameters->FindParam(&oceangridx,&ngrids_ocean,OceanGridXEnum); +- femmodel->parameters->FindParam(&oceangridy,&ngrids_ocean,OceanGridYEnum); +- +- /*Exchange information*/ +- oceanmelt = xNew(ngrids_ocean); +- ISSM_MPI_Recv(oceanmelt,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001007,tomitgcmcomm,&status); +- icebase = xNew(ngrids_ocean); +- for(int i=0;ivertices->NumberOfVertices(); ++ int nels_ice=femmodel->elements->NumberOfElements(); + lat_ice= xNew(ngrids_ice); + lon_ice= xNew(ngrids_ice); ++ icebase= xNew(ngrids_ice); + melt_mesh= xNew(ngrids_ice); ++ ++ /*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); ++ base_oceangrid= xNew(ngrids_ocean); ++ oceanmelt = xNew(ngrids_ocean); ++ + femmodel->vertices->LatLonList(&lat_ice,&lon_ice); ++ ++ /*Interpolate ice base onto ocean grid*/ ++// InputUpdateFromVectorx(femmodel,melt_mesh,BasalforcingsFloatingiceMeltingRateEnum,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); ++ 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"); ++ ISSM_MPI_Recv(oceanmelt,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001007,tomitgcmcomm,&status); ++ ISSM_MPI_Send(base_oceangrid,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001008,tomitgcmcomm); ++ ++ /*Interp melt onto ice grid*/ + InterpFromMeshToMesh2dx(&melt_mesh,index_ocean,oceangridx,oceangridy,ngrids_ocean,nels_ocean, + oceanmelt,ngrids_ocean,1, + lon_ice,lat_ice,ngrids_ice,NULL); ++ ++ for(int i=0;i(lat_ice); + xDelete(lon_ice); + xDelete(x_ice); +@@ -236,6 +245,7 @@ + xDelete(oceangridy); + xDelete(oceanmelt); + xDelete(icebase); ++ xDelete(base_oceangrid); + } + #else + _error_("not supported"); Index: /issm/oecreview/Archive/22755-22818/ISSM-22798-22799.diff =================================================================== --- /issm/oecreview/Archive/22755-22818/ISSM-22798-22799.diff (revision 22819) +++ /issm/oecreview/Archive/22755-22818/ISSM-22798-22799.diff (revision 22819) @@ -0,0 +1,30 @@ +Index: ../trunk-jpl/src/c/cores/transient_core.cpp +=================================================================== +--- ../trunk-jpl/src/c/cores/transient_core.cpp (revision 22798) ++++ ../trunk-jpl/src/c/cores/transient_core.cpp (revision 22799) +@@ -177,7 +177,7 @@ + if(!parcom)_error_("TransferForcing error message: could not find ToMITgcmCommEnum communicator"); + tomitgcmcomm=parcom->GetParameterValue(); + if(my_rank==0){ +- int ngrids_ocean; ++ int ngrids_ocean, nels_ocean; + IssmDouble oceantime; + IssmDouble rho_ice; + IssmDouble *oceanmelt; +@@ -193,7 +193,6 @@ + IssmDouble* melt_mesh = NULL; + int* index_ice= NULL; + int* index_ocean = NULL; +- int nels_ocean; + int ngrids_ice=femmodel->vertices->NumberOfVertices(); + int nels_ice=femmodel->elements->NumberOfElements(); + lat_ice= xNew(ngrids_ice); +@@ -213,7 +212,7 @@ + femmodel->vertices->LatLonList(&lat_ice,&lon_ice); + + /*Interpolate ice base onto ocean grid*/ +-// InputUpdateFromVectorx(femmodel,melt_mesh,BasalforcingsFloatingiceMeltingRateEnum,VertexSIdEnum); ++ 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); Index: /issm/oecreview/Archive/22755-22818/ISSM-22799-22800.diff =================================================================== --- /issm/oecreview/Archive/22755-22818/ISSM-22799-22800.diff (revision 22819) +++ /issm/oecreview/Archive/22755-22818/ISSM-22799-22800.diff (revision 22819) @@ -0,0 +1,53 @@ +Index: ../trunk-jpl/src/c/cores/transient_core.cpp +=================================================================== +--- ../trunk-jpl/src/c/cores/transient_core.cpp (revision 22799) ++++ ../trunk-jpl/src/c/cores/transient_core.cpp (revision 22800) +@@ -195,10 +195,6 @@ + int* index_ocean = NULL; + int ngrids_ice=femmodel->vertices->NumberOfVertices(); + int nels_ice=femmodel->elements->NumberOfElements(); +- lat_ice= xNew(ngrids_ice); +- lon_ice= xNew(ngrids_ice); +- icebase= xNew(ngrids_ice); +- melt_mesh= xNew(ngrids_ice); + + /*Recover mesh and inputs needed*/ + femmodel->parameters->FindParam(&rho_ice,MaterialsRhoIceEnum); +@@ -206,8 +202,6 @@ + femmodel->parameters->FindParam(&oceangridx,&ngrids_ocean,OceanGridXEnum); + femmodel->parameters->FindParam(&oceangridy,&ngrids_ocean,OceanGridYEnum); + BamgTriangulatex(&index_ocean,&nels_ocean,oceangridx,oceangridy,ngrids_ocean); +- base_oceangrid= xNew(ngrids_ocean); +- oceanmelt = xNew(ngrids_ocean); + + femmodel->vertices->LatLonList(&lat_ice,&lon_ice); + +@@ -219,9 +213,11 @@ + + /*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"); + ISSM_MPI_Recv(oceanmelt,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001007,tomitgcmcomm,&status); ++ for(int i=0;i(index_ice); ++ xDelete(index_ocean); + xDelete(lat_ice); + xDelete(lon_ice); + xDelete(x_ice); +@@ -239,7 +237,6 @@ + xDelete(y_ice); + xDelete(z_ice); + xDelete(melt_mesh); +- xDelete(index_ice); + xDelete(oceangridx); + xDelete(oceangridy); + xDelete(oceanmelt); Index: /issm/oecreview/Archive/22755-22818/ISSM-22800-22801.diff =================================================================== --- /issm/oecreview/Archive/22755-22818/ISSM-22800-22801.diff (revision 22819) +++ /issm/oecreview/Archive/22755-22818/ISSM-22800-22801.diff (revision 22819) @@ -0,0 +1,83 @@ +Index: ../trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp +=================================================================== +--- ../trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp (revision 22800) ++++ ../trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp (revision 22801) +@@ -17,6 +17,11 @@ + 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)); ++ } + else if(strcmp(string_in,"MaterialsRheologyB")==0){ + const char* field = "md.materials.rheology_B"; + input_enum = MaterialsRheologyBEnum; +Index: ../trunk-jpl/src/c/cores/controladm1qn3_core.cpp +=================================================================== +--- ../trunk-jpl/src/c/cores/controladm1qn3_core.cpp (revision 22800) ++++ ../trunk-jpl/src/c/cores/controladm1qn3_core.cpp (revision 22801) +@@ -485,7 +485,6 @@ + femmodel->parameters->SetParam(step,StepEnum); + femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,JEnum,mystruct.Jlist,(*mystruct.i),mystruct.N,1,0)); + +- /*TEMP START*/ + int offset = 0; + for(int i=0;iresults->AddObject(new GenericExternalResult(femmodel->results->Size()+1,GradientEnum,G,numberofvertices,intn/numberofvertices,1,0)); +- +- femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,InversionControlParametersEnum,X,numberofvertices,intn/numberofvertices,1,0)); +- +- //femmodel->OutputControlsx(&femmodel->results); + } + else{ + femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,JEnum,mystruct.Jlist,(*mystruct.i),mystruct.N,0,0)); +Index: ../trunk-jpl/src/c/classes/Elements/Element.cpp +=================================================================== +--- ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 22800) ++++ ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 22801) +@@ -1399,7 +1399,11 @@ + int start = 0; + 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(input_enum,values,P1Enum)); +@@ -1678,12 +1682,12 @@ + 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)); + break; Index: /issm/oecreview/Archive/22755-22818/ISSM-22801-22802.diff =================================================================== --- /issm/oecreview/Archive/22755-22818/ISSM-22801-22802.diff (revision 22819) +++ /issm/oecreview/Archive/22755-22818/ISSM-22801-22802.diff (revision 22819) @@ -0,0 +1,22 @@ +Index: ../trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp +=================================================================== +--- ../trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 22801) ++++ ../trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 22802) +@@ -104,7 +104,6 @@ + parameters->AddObject(iomodel->CopyConstantObject("md.transient.isslr",TransientIsslrEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.transient.iscoupler",TransientIscouplerEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.transient.isoceancoupling",TransientIsoceancouplingEnum)); +- parameters->AddObject(iomodel->CopyConstantObject("md.materials.rho_ice",MaterialsRhoIceEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.transient.amr_frequency",TransientAmrFrequencyEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.gia.cross_section_shape",GiaCrossSectionShapeEnum)); + +@@ -291,6 +290,9 @@ + iomodel->DeleteData(&requestedoutputs,numoutputs,"md.steadystate.requested_outputs"); + + 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"); + parameters->AddObject(new IntParam(DamageEvolutionNumRequestedOutputsEnum,numoutputs)); Index: /issm/oecreview/Archive/22755-22818/ISSM-22802-22803.diff =================================================================== --- /issm/oecreview/Archive/22755-22818/ISSM-22802-22803.diff (revision 22819) +++ /issm/oecreview/Archive/22755-22818/ISSM-22802-22803.diff (revision 22819) @@ -0,0 +1,133 @@ +Index: ../trunk-jpl/examples/EsaGRACE/runme.m +=================================================================== +--- ../trunk-jpl/examples/EsaGRACE/runme.m (revision 22802) ++++ ../trunk-jpl/examples/EsaGRACE/runme.m (revision 22803) +@@ -1,6 +1,6 @@ + + clear all; +-steps=[5]; % [1:6]; ++steps=[0]; % [1:5]; + + if any(steps==0) % Download GRACE land_mass data {{{ + disp(' Step 0: Download GRACE land_mass data'); +@@ -41,12 +41,12 @@ + disp(' Step 2: Define loads in meters of ice height equivalent'); + md = loadmodel('./Models/EsaGRACE.Mesh'); + ++ % define time interval for analysis + year_month = 2007+15/365; +- time_min=year_month; +- time_max=year_month; +- ++ time_range = [year_month year_month]; ++ + % map GRACE water load on to the mesh for the seleted month(s) +- water_load = grace(md.mesh.elements,md.mesh.lat,md.mesh.long,time_min,time_max); ++ 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 + +@@ -116,10 +116,13 @@ + disp(' Step 5: Plot solutions'); + md = loadmodel('./Models/EsaGRACE.Solution'); + +- % solutions. +- ur = md.results.EsaSolution.EsaUmotion*1000; % [mm] +- un = md.results.EsaSolution.EsaNmotion*1000; % [mm] +- ue = md.results.EsaSolution.EsaEmotion*1000; % [mm] ++ % loads and solutions. ++ sol1 = md.esa.deltathickness*100; % WEH 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]'}; + + res = 1.0; % degree + +@@ -127,45 +130,50 @@ + [lat_grid, lon_grid] = meshgrid(linspace(-90,90,180/res), linspace(-180,180,360/res)); + sol_grid = zeros(size(lat_grid)); + +- sol = ue; ++ for kk=1:4 ++ sol=eval(sprintf('sol%d',kk)); ++ ++ % if data are on elements, map those on to the vertices {{{ ++ if length(sol)==md.mesh.numberofelements ++ % map on to the vertices ++ 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 % }}} + +- % Make a interpolation object +- F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol); +- F.Method = 'linear'; +- F.ExtrapolationMethod = 'linear'; ++ % Make a interpolation object ++ F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol); ++ F.Method = 'linear'; ++ F.ExtrapolationMethod = 'linear'; + +- % Do the interpolation to get gridded solutions... +- sol_grid = F(lat_grid, lon_grid); +- sol_grid(isnan(sol_grid))=0; ++ % Do the interpolation to get gridded solutions... ++ sol_grid = F(lat_grid, lon_grid); ++ sol_grid(isnan(sol_grid))=0; ++ sol_grid(lat_grid>85 & sol_grid==0) =NaN; % set polar unphysical 0s to Nan + +- % set polar unphysical 0s to Nan +- 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; ++ plot(long,lat,'k'); hold off; ++ c1=colorbar; ++ colormap(jet); ++ xlim([-180 180]); ++ ylim([-90 90]); ++ grid on; ++ title(sol_name(kk)); ++ set(gcf,'color','w'); + +- 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; +- %caxis([-0.3 0.3]) +- hold on +- plot(long,lat,'k'); +- %geoshow(lat,long,'DisplayType','polygon','FaceColor','white'); +- hold off; +- c1=colorbar; +- colormap(jet); +- xlim([-180 180]); +- ylim([-90 90]); +- grid on; +- title(['Average change in relative sea-level [mm/yr]']); +- set(gcf,'color','w'); ++ %export_fig('Fig5.pdf'); ++ end + +- %export_fig('Fig5.pdf'); +- + end % }}} + +-if any(steps==6) % {{{ Transient +- disp(' Step 6: Transient run'); +- +-end % }}} +- Index: /issm/oecreview/Archive/22755-22818/ISSM-22803-22804.diff =================================================================== --- /issm/oecreview/Archive/22755-22818/ISSM-22803-22804.diff (revision 22819) +++ /issm/oecreview/Archive/22755-22818/ISSM-22803-22804.diff (revision 22819) @@ -0,0 +1,279 @@ +Index: ../trunk-jpl/examples/SlrGRACE/runme.m +=================================================================== +--- ../trunk-jpl/examples/SlrGRACE/runme.m (nonexistent) ++++ ../trunk-jpl/examples/SlrGRACE/runme.m (revision 22804) +@@ -0,0 +1,274 @@ ++ ++clear all; ++steps=[5]; % [1:6]; ++ ++if any(steps==0) % Download GRACE land_mass data {{{ ++ disp(' Step 0: Download GRACE land_mass data'); ++ ++ % Connect to ftp server and download ++ f = ftp('podaac-ftp.jpl.nasa.gov'); ++ %dir(f) ++ cd(f,'allData/tellus/L3/land_mass/RL05/netcdf/'); ++ mget(f,'GRCTellus.JPL.200204_201701.LND.RL05_1.DSTvSCS1411.nc'); ++ ++ % display the content of the netcdf file. ++ %ncdisp('GRCTellus.JPL.200204_201701.LND.RL05_1.DSTvSCS1411.nc') ++ ++end % }}} ++ ++if any(steps==1) % Global mesh creation {{{ ++ disp(' Step 1: Global mesh creation'); ++ ++ numrefine=1; ++ resolution=150*1e3; % inital resolution [m]. It determines, e.g., whether we capture small islands. ++ 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] ++ ++ %mesh earth: ++ md=model; ++ md.mask=maskpsl(); % use maskpsl class (instead of mask) to store the ocean function as a ocean_levelset ++ md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3); % attributes should be in km. ++ ++ for i=1:numrefine, ++ ++ %figure out mask: ++ md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); ++ ++ %figure out distance to the coastline, in lat,long (not x,y,z): ++ 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 ++ %figure out nearest coastline (using the great circle distance) ++ 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); % max size 1000 km ++ %use distance to the coastline to refine mesh: ++ md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3,'refine',md.mesh,'refinemetric',dist); ++ end ++ ++ %figure out mask: ++ 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'); ++ %export_fig('Fig1.pdf'); ++ ++end % }}} ++ ++if any(steps==2) % Define loads {{{ ++ disp(' Step 2: Define loads in meters of ice height equivalent'); ++ md = loadmodel('./Models/SlrGRACE.Mesh'); ++ ++ % define time interval for analysis ++ year_month = 2007+15/365; ++ time_range = [year_month year_month]; ++ ++ % map GRACE water load on to the mesh for the seleted month(s) ++ 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; % ice height equivalent ++ ++ save ./Models/SlrGRACE.Loads md; ++ ++ plotmodel (md,'data',md.slr.deltathickness,... ++ 'view',[90 -90],'caxis',[-.1 .1],... ++ 'title','Ice height equivalent [m]'); ++ %export_fig('Fig2.pdf'); ++ ++end % }}} ++ ++if any(steps==3) % Parameterization {{{ ++ disp(' Step 3: Parameterization'); ++ md = loadmodel('./Models/SlrGRACE.Loads'); ++ ++ % Love numbers and reference frame: CF or CM (choose one!) ++ nlove=10001; % up to 10,000 degree ++ 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)=[]; ++ ++ % Mask: for computational efficiency only those elements that have loads are convolved! ++ md.mask.groundedice_levelset = ones(md.mesh.numberofvertices,1); % 1 = ice is grounnded ++ md.mask.ice_levelset = ones(md.mesh.numberofvertices,1); ++ pos=find(md.slr.deltathickness~=0); ++ md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; % -1 = ice loads ++ md.mask.land_levelset = 1-md.mask.ocean_levelset; ++ ++ % initialize ++ md.slr.sealevel=zeros(md.mesh.numberofvertices,1); ++ % steric rate ++ md.slr.steric_rate=zeros(md.mesh.numberofvertices,1); ++ % solution is scaled according to "exact" area of the oceans ++ md.slr.ocean_area_scaling = 1; ++ ++ % convergence criteria ++ %md.slr.reltol=NaN; ++ %md.slr.abstol=1e-4; ++ ++ %% IGNORE BUT DO NOT DELETE %% {{{ ++ % Geometry: Important only when you want to couple with Ice Flow Model ++ 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; ++ % Materials: ++ 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); ++ % Miscellaneous: ++ md.miscellaneous.name='SlrGRACE'; ++ %% IGNORE BUT DO NOT DELETE %% }}} ++ ++ save ./Models/SlrGRACE.Parameterization md; ++ ++end % }}} ++ ++if any(steps==4) % Solve {{{ ++ disp(' Step 4: Solve Slr solver'); ++ md = loadmodel('./Models/SlrGRACE.Parameterization'); ++ ++ % What physics to capture? ++ md.slr.rigid=1; ++ md.slr.elastic=1; ++ md.slr.rotation=1; ++ ++ % Request outputs ++% md.slr.requested_outputs = {'SlrUmotion','SlrNmotion','SlrEmotion'}; ++ ++ % Cluster info ++ md.cluster=generic('name',oshostname(),'np',3); ++ md.verbose=verbose('111111111'); ++ ++ % Solve ++ md=solve(md,'Slr'); ++ ++ save ./Models/SlrGRACE.Solution md; ++ ++end % }}} ++ ++if any(steps==5) % Plot solutions {{{ ++ disp(' Step 5: Plot solutions'); ++ md = loadmodel('./Models/SlrGRACE.Solution'); ++ ++ % loads and solutions. ++ sol1 = md.slr.deltathickness*100; % WEH cm ++ sol2 = md.results.SealevelriseSolution.Sealevel*1000; % [mm] ++ sol_name={'Change in water equivalent height [cm]', 'Relative sea-level [mm]'}; ++ ++ res = 1.0; % degree ++ ++ % Make a grid of lats and lons, based on the min and max of the original vectors ++ [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 data are on elements, map those on to the vertices {{{ ++ if length(sol)==md.mesh.numberofelements ++ % map on to the vertices ++ 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 % }}} ++ ++ % Make a interpolation object ++ F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol); ++ F.Method = 'linear'; ++ F.ExtrapolationMethod = 'linear'; ++ ++ % Do the interpolation to get gridded solutions... ++ sol_grid = F(lat_grid, lon_grid); ++ sol_grid(isnan(sol_grid))=0; ++ sol_grid(lat_grid>85 & sol_grid==0) =NaN; % set polar unphysical 0s to 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; ++ plot(long,lat,'k'); hold off; ++ c1=colorbar; ++ colormap(jet); ++ xlim([-180 180]); ++ ylim([-90 90]); ++ grid on; ++ title(sol_name(kk)); ++ set(gcf,'color','w'); ++ ++ %export_fig('Fig5.pdf'); ++ end ++ ++end % }}} ++ ++if any(steps==6) % Transient {{{ ++ disp(' Step 6: Transient run'); ++ md = loadmodel('./Models/SlrGRACE.Parameterization'); ++ ++ % read GRACE data during the chosen time period << steps=2 >> ++ disp('Projecting grace 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)); ++ ++ % define ice load ++ [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; % ice height equivalent ++ md.slr.deltathickness(ne+1,:)=[1:nt]; % times ++ ++ % define transient run ++ md.transient.issmb=0; ++ md.transient.ismasstransport=0; ++ md.transient.isstressbalance=0; ++ md.transient.isthermal=0; ++ md.transient.isslr=1; ++ ++ % time stepping, output frequency etc. ++ 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; ++ ++ % Request outputs ++ md.slr.requested_outputs = {'SlrUmotion','SlrNmotion','SlrEmotion'}; ++ ++ % Cluster info ++ md.cluster=generic('name',oshostname(),'np',3); ++ md.verbose=verbose('111111111'); ++ ++ % Solve ++ md=solve(md,'tr'); ++ ++end % }}} ++ Index: /issm/oecreview/Archive/22755-22818/ISSM-22804-22805.diff =================================================================== --- /issm/oecreview/Archive/22755-22818/ISSM-22804-22805.diff (revision 22819) +++ /issm/oecreview/Archive/22755-22818/ISSM-22804-22805.diff (revision 22819) @@ -0,0 +1,205 @@ +Index: ../trunk-jpl/examples/SlrFarrell/runme.m +=================================================================== +--- ../trunk-jpl/examples/SlrFarrell/runme.m (nonexistent) ++++ ../trunk-jpl/examples/SlrFarrell/runme.m (revision 22805) +@@ -0,0 +1,200 @@ ++ ++clear all; ++steps=[2]; % ++ ++if any(steps==1) % Global mesh creation {{{ ++ disp(' Step 1: Global mesh creation'); ++ ++ numrefine=1; ++ resolution=150*1e3; % inital resolution [m]. It determines, e.g., whether we capture small islands. ++ 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] ++ ++ %mesh earth: ++ md=model; ++ md.mask=maskpsl(); % use maskpsl class (instead of mask) to store the ocean function as a ocean_levelset ++ md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3); % attributes should be in km. ++ ++ for i=1:numrefine, ++ ++ %figure out mask: ++ md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); ++ ++ %figure out distance to the coastline, in lat,long (not x,y,z): ++ 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 ++ %figure out nearest coastline (using the great circle distance) ++ 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); % max size 1000 km ++ %use distance to the coastline to refine mesh: ++ md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3,'refine',md.mesh,'refinemetric',dist); ++ end ++ ++ %figure out mask: ++ 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'); ++ %export_fig('Fig1.pdf'); ++ ++end % }}} ++ ++if any(steps==2) % Define source {{{ ++ disp(' Step 2: Define source as in Farrell, 1972, Figure 1'); ++ md = loadmodel('./Models/SlrFarrell.Mesh'); ++ ++ % initial sea-level: 1 m RSL everywhere. ++ md.slr.sealevel=md.mask.ocean_levelset; ++ ++ save ./Models/SlrFarrell.Loads md; ++ ++ plotmodel (md,'data',md.slr.sealevel,'view',[90 90],... ++ 'title#all','Initial sea-level [m]'); ++ %export_fig('Fig2.pdf'); ++ ++end % }}} ++ ++if any(steps==3) % Parameterization {{{ ++ disp(' Step 3: Parameterization'); ++ md = loadmodel('./Models/SlrFarrell.Loads'); ++ ++ % Love numbers and reference frame: CF or CM (choose one!) ++ nlove=10001; % up to 10,000 degree ++ md.slr.love_h = love_numbers('h','CM'); md.slr.love_h(nlov+1:end)=[]; ++ md.slr.love_k = love_numbers('k','CM'); md.slr.love_k(nlov+1:end)=[]; ++ md.slr.love_l = love_numbers('l','CM'); md.slr.love_l(nlov+1:end)=[]; ++ ++ % Mask: for computational efficiency only those elements that have loads are convolved! ++ md.mask.groundedice_levelset = ones(md.mesh.numberofvertices,1); % 1 = ice is grounnded ++ md.mask.ice_levelset = ones(md.mesh.numberofvertices,1); ++ pos=find(md.slr.deltathickness~=0); ++ md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; % -1 = ice loads ++ md.mask.land_levelset = 1-md.mask.ocean_levelset; ++ ++ %% IGNORE BUT DO NOT DELETE %% {{{ ++ % Geometry: Important only when you want to couple with Ice Flow Model ++ 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; ++ % Materials: ++ 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); ++ % Miscellaneous: ++ md.miscellaneous.name='SlrFarrell'; ++ %% IGNORE BUT DO NOT DELETE %% }}} ++ ++ save ./Models/SlrFarrell.Parameterization md; ++ ++end % }}} ++ ++if any(steps==4) % Solve {{{ ++ disp(' Step 4: Solve Slr solver'); ++ md = loadmodel('./Models/SlrFarrell.Parameterization'); ++ ++ % Request outputs ++ md.slr.requested_outputs = {'SlrUmotion','SlrNmotion','SlrEmotion'}; ++ ++ % Cluster info ++ md.cluster=generic('name',oshostname(),'np',3); ++ md.verbose=verbose('111111111'); ++ ++ % Solve ++ md=solve(md,'Slr'); ++ ++ save ./Models/SlrFarrell.Solution md; ++ ++end % }}} ++ ++if any(steps==5) % Plot solutions {{{ ++ disp(' Step 5: Plot solutions'); ++ md = loadmodel('./Models/SlrFarrell.Solution'); ++ ++ % loads and solutions. ++ sol1 = md.slr.deltathickness*100; % WEH cm ++ sol2 = md.results.SlrSolution.SlrUmotion*1000; % [mm] ++ sol3 = md.results.SlrSolution.SlrNmotion*1000; % [mm] ++ sol4 = md.results.SlrSolution.SlrEmotion*1000; % [mm] ++ sol_name={'Change in water equivalent height [cm]', 'Vertical displacement [mm]',... ++ 'Horizontal (NS) displacement [mm]', 'Horizontal (EW) displacement [mm]'}; ++ ++ res = 1.0; % degree ++ ++ % Make a grid of lats and lons, based on the min and max of the original vectors ++ [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 data are on elements, map those on to the vertices {{{ ++ if length(sol)==md.mesh.numberofelements ++ % map on to the vertices ++ 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 % }}} ++ ++ % Make a interpolation object ++ F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol); ++ F.Method = 'linear'; ++ F.ExtrapolationMethod = 'linear'; ++ ++ % Do the interpolation to get gridded solutions... ++ sol_grid = F(lat_grid, lon_grid); ++ sol_grid(isnan(sol_grid))=0; ++ sol_grid(lat_grid>85 & sol_grid==0) =NaN; % set polar unphysical 0s to 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; ++ plot(long,lat,'k'); hold off; ++ c1=colorbar; ++ colormap(jet); ++ xlim([-180 180]); ++ ylim([-90 90]); ++ grid on; ++ title(sol_name(kk)); ++ set(gcf,'color','w'); ++ ++ %export_fig('Fig5.pdf'); ++ end ++ ++end % }}} ++ ++ Index: /issm/oecreview/Archive/22755-22818/ISSM-22805-22806.diff =================================================================== --- /issm/oecreview/Archive/22755-22818/ISSM-22805-22806.diff (revision 22819) +++ /issm/oecreview/Archive/22755-22818/ISSM-22805-22806.diff (revision 22819) @@ -0,0 +1,68 @@ +Index: ../trunk-jpl/src/m/mesh/planet/gmsh/gmshplanet.py +=================================================================== +--- ../trunk-jpl/src/m/mesh/planet/gmsh/gmshplanet.py (revision 22805) ++++ ../trunk-jpl/src/m/mesh/planet/gmsh/gmshplanet.py (revision 22806) +@@ -63,21 +63,21 @@ + fid.write('Circle(11) = {4,1,6};\n') + 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') + fid.write('Physical Surface(1) = {28,26,16,14,20,24,22,18};\n') +Index: ../trunk-jpl/src/m/mesh/planet/gmsh/gmshplanet.m +=================================================================== +--- ../trunk-jpl/src/m/mesh/planet/gmsh/gmshplanet.m (revision 22805) ++++ ../trunk-jpl/src/m/mesh/planet/gmsh/gmshplanet.m (revision 22806) +@@ -57,21 +57,21 @@ + fprintf(fid,'Circle(11) = {4,1,6};\n'); + 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'); + fprintf(fid,'Physical Surface(1) = {28,26,16,14,20,24,22,18};\n'); Index: /issm/oecreview/Archive/22755-22818/ISSM-22806-22807.diff =================================================================== --- /issm/oecreview/Archive/22755-22818/ISSM-22806-22807.diff (revision 22819) +++ /issm/oecreview/Archive/22755-22818/ISSM-22806-22807.diff (revision 22819) @@ -0,0 +1,13 @@ +Index: ../trunk-jpl/src/m/mesh/planet/gmsh/gmshplanet.m +=================================================================== +--- ../trunk-jpl/src/m/mesh/planet/gmsh/gmshplanet.m (revision 22806) ++++ ../trunk-jpl/src/m/mesh/planet/gmsh/gmshplanet.m (revision 22807) +@@ -106,7 +106,7 @@ + gmshpath = ''; + for i=paths + if exist(i{1},'file'), +- gmshpath = i{1} ++ gmshpath = i{1}; + end + end + if isempty(gmshpath), Index: /issm/oecreview/Archive/22755-22818/ISSM-22807-22808.diff =================================================================== --- /issm/oecreview/Archive/22755-22818/ISSM-22807-22808.diff (revision 22819) +++ /issm/oecreview/Archive/22755-22818/ISSM-22807-22808.diff (revision 22819) @@ -0,0 +1,84 @@ +Index: ../trunk-jpl/src/m/classes/slr.m +=================================================================== +--- ../trunk-jpl/src/m/classes/slr.m (revision 22807) ++++ ../trunk-jpl/src/m/classes/slr.m (revision 22808) +@@ -121,10 +121,10 @@ + function disp(self) % {{{ + 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'); + fielddisplay(self,'love_k','load Love number for gravitational potential perturbation'); +@@ -138,8 +138,8 @@ + fielddisplay(self,'rigid','rigid earth graviational potential perturbation'); + 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'); + fielddisplay(self,'requested_outputs','additional outputs requested'); +Index: ../trunk-jpl/src/m/classes/slr.py +=================================================================== +--- ../trunk-jpl/src/m/classes/slr.py (revision 22807) ++++ ../trunk-jpl/src/m/classes/slr.py (revision 22808) +@@ -42,10 +42,10 @@ + #}}} + 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')) + string="%s\n%s"%(string,fielddisplay(self,'love_k','load Love number for gravitational potential perturbation')) +@@ -59,8 +59,8 @@ + string="%s\n%s"%(string,fielddisplay(self,'rigid','rigid earth graviational potential perturbation')) + 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')) + string="%s\n%s"%(string,fielddisplay(self,'requested_outputs','additional outputs requested')) +Index: ../trunk-jpl/src/m/classes/slr.js +=================================================================== +--- ../trunk-jpl/src/m/classes/slr.js (revision 22807) ++++ ../trunk-jpl/src/m/classes/slr.js (revision 22808) +@@ -87,10 +87,10 @@ + + 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'); + fielddisplay(this,'love_k','load Love number for gravitational potential perturbation'); +@@ -104,8 +104,8 @@ + fielddisplay(this,'rigid','rigid earth graviational potential perturbation'); + 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'); + fielddisplay(this,'requested_outputs','additional outputs requested'); Index: /issm/oecreview/Archive/22755-22818/ISSM-22808-22809.diff =================================================================== --- /issm/oecreview/Archive/22755-22818/ISSM-22808-22809.diff (revision 22819) +++ /issm/oecreview/Archive/22755-22818/ISSM-22808-22809.diff (revision 22819) @@ -0,0 +1,413 @@ +Index: ../trunk-jpl/examples/Functions/wahr.m +=================================================================== +--- ../trunk-jpl/examples/Functions/wahr.m (nonexistent) ++++ ../trunk-jpl/examples/Functions/wahr.m (revision 22809) +@@ -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 ++ + +Property changes on: ../trunk-jpl/examples/Functions/wahr.m +___________________________________________________________________ +Added: svn:executable +## -0,0 +1 ## ++* +\ No newline at end of property +Index: ../trunk-jpl/examples/SlrGRACE/runme.m +=================================================================== +--- ../trunk-jpl/examples/SlrGRACE/runme.m (revision 22808) ++++ ../trunk-jpl/examples/SlrGRACE/runme.m (revision 22809) +@@ -1,7 +1,9 @@ + + clear all; +-steps=[5]; % [1:6]; ++addpath('../Data','../Functions'); + ++steps=[7]; % [0:7]; ++ + if any(steps==0) % Download GRACE land_mass data {{{ + disp(' Step 0: Download GRACE land_mass data'); + +@@ -11,6 +13,8 @@ + cd(f,'allData/tellus/L3/land_mass/RL05/netcdf/'); + mget(f,'GRCTellus.JPL.200204_201701.LND.RL05_1.DSTvSCS1411.nc'); + ++ !mv *.nc '../Data/' ++ + % display the content of the netcdf file. + %ncdisp('GRCTellus.JPL.200204_201701.LND.RL05_1.DSTvSCS1411.nc') + +@@ -31,7 +35,7 @@ + md.mask=maskpsl(); % use maskpsl class (instead of mask) to store the ocean function as a ocean_levelset + md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3); % attributes should be in km. + +- for i=1:numrefine, ++ for i=1:numrefine; % refine mesh... {{{ + + %figure out mask: + md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); +@@ -65,7 +69,7 @@ + dist=min(maxdistance,distance); % max size 1000 km + %use distance to the coastline to refine mesh: + md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3,'refine',md.mesh,'refinemetric',dist); +- end ++ end % }}} + + %figure out mask: + md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); +@@ -155,9 +159,6 @@ + md.slr.elastic=1; + md.slr.rotation=1; + +- % Request outputs +-% md.slr.requested_outputs = {'SlrUmotion','SlrNmotion','SlrEmotion'}; +- + % Cluster info + md.cluster=generic('name',oshostname(),'np',3); + md.verbose=verbose('111111111'); +@@ -205,7 +206,7 @@ + F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol); + F.Method = 'linear'; + F.ExtrapolationMethod = 'linear'; +- ++ + % Do the interpolation to get gridded solutions... + sol_grid = F(lat_grid, lon_grid); + sol_grid(isnan(sol_grid))=0; +@@ -213,15 +214,22 @@ + + 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; ++ gcf; load coast; cla; ++ pcolor(lon_grid,lat_grid,sol_grid); shading flat; hold on; ++ % mask out land or oceans {{{ ++ 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; ++ % define colormap, caxis, xlim etc {{{ + c1=colorbar; +- colormap(jet); ++ 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'); +@@ -236,7 +244,7 @@ + md = loadmodel('./Models/SlrGRACE.Parameterization'); + + % read GRACE data during the chosen time period << steps=2 >> +- disp('Projecting grace loads onto the mesh...'); ++ 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)); + +@@ -260,9 +268,6 @@ + md.timestepping.interp_forcings=0; + md.settings.output_frequency=1; + +- % Request outputs +- md.slr.requested_outputs = {'SlrUmotion','SlrNmotion','SlrEmotion'}; +- + % Cluster info + md.cluster=generic('name',oshostname(),'np',3); + md.verbose=verbose('111111111'); +@@ -270,5 +275,101 @@ + % Solve + md=solve(md,'tr'); + ++ save ./Models/SlrGRACE.Transient md; ++ + end % }}} + ++if any(steps==7) % Plot transient {{{ ++ disp(' Step 7: Plot transient'); ++ md = loadmodel('./Models/SlrGRACE.Transient'); ++ ++ time = md.slr.deltathickness(end,:); % year ++ ++ % loads and solutions. ++ 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] ++ %% something weired happened here!!! first month solution is duplicated. Use offset of 1. Will fix it asap! ++ 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_loads.avi','movie_fingerprints.avi'}; ++ ++ res = 1.0; % degree ++ ++ % Make a grid of lats and lons, based on the min and max of the original vectors ++ [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 data are on elements, map those on to the vertices {{{ ++ if length(sol)==md.mesh.numberofelements ++ % map on to the vertices ++ 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) ++ % Make a interpolation object ++ F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol(:,jj)); ++ F.Method = 'linear'; ++ F.ExtrapolationMethod = 'linear'; ++ ++ % Do the interpolation to get gridded solutions... ++ sol_grid = F(lat_grid, lon_grid); ++ sol_grid(isnan(sol_grid))=0; ++ sol_grid(lat_grid>85 & sol_grid==0) =NaN; % set polar unphysical 0s to 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; ++ % mask out land or oceans {{{ ++ 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; ++ % define colormap, caxis, xlim etc {{{ ++ 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 % close current figure ++ end ++ ++ % Close the file. ++ 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: ../trunk-jpl/examples/Functions/grace.m +=================================================================== +--- ../trunk-jpl/examples/Functions/grace.m (nonexistent) ++++ ../trunk-jpl/examples/Functions/grace.m (revision 22809) +@@ -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; ++ + +Property changes on: ../trunk-jpl/examples/Functions/grace.m +___________________________________________________________________ +Added: svn:executable +## -0,0 +1 ## ++* +\ No newline at end of property +Index: ../trunk-jpl/examples/EsaGRACE/runme.m +=================================================================== +--- ../trunk-jpl/examples/EsaGRACE/runme.m (revision 22808) ++++ ../trunk-jpl/examples/EsaGRACE/runme.m (revision 22809) +@@ -1,6 +1,6 @@ + + clear all; +-steps=[0]; % [1:5]; ++steps=[0:5]; % [1:5]; + + if any(steps==0) % Download GRACE land_mass data {{{ + disp(' Step 0: Download GRACE land_mass data'); Index: /issm/oecreview/Archive/22755-22818/ISSM-22809-22810.diff =================================================================== --- /issm/oecreview/Archive/22755-22818/ISSM-22809-22810.diff (revision 22819) +++ /issm/oecreview/Archive/22755-22818/ISSM-22809-22810.diff (revision 22819) @@ -0,0 +1,191 @@ +Index: ../trunk-jpl/examples/EsaGRACE/grace.m +=================================================================== +--- ../trunk-jpl/examples/EsaGRACE/grace.m (revision 22809) ++++ ../trunk-jpl/examples/EsaGRACE/grace.m (nonexistent) +@@ -1,121 +0,0 @@ +-% 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; +- + +Property changes on: ../trunk-jpl/examples/EsaGRACE/grace.m +___________________________________________________________________ +Deleted: svn:executable +## -1 +0,0 ## +-* +\ No newline at end of property +Index: ../trunk-jpl/examples/EsaWahr/wahr.m +=================================================================== +--- ../trunk-jpl/examples/EsaWahr/wahr.m (revision 22809) ++++ ../trunk-jpl/examples/EsaWahr/wahr.m (nonexistent) +@@ -1,46 +0,0 @@ +-% 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 +- + +Property changes on: ../trunk-jpl/examples/EsaWahr/wahr.m +___________________________________________________________________ +Deleted: svn:executable +## -1 +0,0 ## +-* +\ No newline at end of property Index: /issm/oecreview/Archive/22755-22818/ISSM-22810-22811.diff =================================================================== --- /issm/oecreview/Archive/22755-22818/ISSM-22810-22811.diff (revision 22819) +++ /issm/oecreview/Archive/22755-22818/ISSM-22810-22811.diff (revision 22819) @@ -0,0 +1,119 @@ +Index: ../trunk-jpl/examples/EsaWahr/runme.m +=================================================================== +--- ../trunk-jpl/examples/EsaWahr/runme.m (revision 22810) ++++ ../trunk-jpl/examples/EsaWahr/runme.m (revision 22811) +@@ -1,8 +1,10 @@ + + clear all; +-steps=[1:6]; % [1:6]; ++addpath('../Functions'); + +-if any(steps==0) ++steps=[0:6]; % [0:6]; ++ ++if any(steps==0) % Simple mesh creation {{{ + disp(' Step 0: Mesh creation'); + + % Generate initial uniform mesh +@@ -13,9 +15,9 @@ + plotmodel(md,'data','mesh'); + %export_fig('Fig0.pdf'); + +-end ++end % }}} + +-if any(steps==1) ++if any(steps==1) % {{{ Anisotropic mesh creation + disp(' Step 1: Anisotropic mesh creation'); + + % Generate initial uniform mesh +@@ -35,9 +37,9 @@ + plotmodel (md,'data','mesh'); + %export_fig('Fig1.pdf'); + +-end ++end % }}} + +-if any(steps==2) ++if any(steps==2) % Define loads {{{ + disp(' Step 2: Define loads'); + md = loadmodel('./Models/EsaWahr.Mesh'); + +@@ -60,9 +62,9 @@ + plotmodel (md,'data',md.esa.deltathickness,'title','Ice height equivalent [m]'); + %export_fig('Fig2.pdf'); + +-end ++end % }}} + +-if any(steps==3) ++if any(steps==3) % Parameterization {{{ + disp(' Step 3: Parameterization'); + md = loadmodel('./Models/EsaWahr.Loads'); + +@@ -75,7 +77,7 @@ + md.mask.ice_levelset = -ones(md.mesh.numberofvertices,1); % -1 = ice loads + md.mask.groundedice_levelset = ones(md.mesh.numberofvertices,1); % 1 = ice is grounnded + +- %% IGNORE BUT DO NOT DELETE %% ++ %% IGNORE BUT DO NOT DELETE %% {{{ + % Geometry: Important only when you want to couple with Ice Flow Model + di=md.materials.rho_ice/md.materials.rho_water; + md.geometry.thickness=ones(md.mesh.numberofvertices,1); +@@ -88,13 +90,13 @@ + md.materials.rheology_n=3*ones(md.mesh.numberofelements,1); + % Miscellaneous: + md.miscellaneous.name='EsaWahr'; +- %% IGNORE BUT DO NOT DELETE %% ++ %% IGNORE BUT DO NOT DELETE %% }}} + + save ./Models/EsaWahr.Parameterization md; + +-end ++end % }}} + +-if any(steps==4) ++if any(steps==4) % Solve {{{ + disp(' Step 4: Solve Esa solver'); + md = loadmodel('./Models/EsaWahr.Parameterization'); + +@@ -110,9 +112,9 @@ + + save ./Models/EsaWahr.Solution md; + +-end ++end % }}} + +-if any(steps==5) ++if any(steps==5) % Plot solutions {{{ + disp(' Step 5: Plot solutions'); + md = loadmodel('./Models/EsaWahr.Solution'); + +@@ -147,9 +149,9 @@ + + %export_fig('Fig5.pdf'); + +-end ++end % }}} + +-if any(steps==6) ++if any(steps==6) % Compare results against semi-analytic solutions {{{ + disp(' Step 6: Compare results against Wahr semi-analytic solutions'); + md = loadmodel('./Models/EsaWahr.Solution'); + +@@ -188,7 +190,7 @@ + h2=plot(xi/1000,horz_track*917/1000,'-c','linewidth',3,'parent',axes1); + % first box + ag1 = gca; +- leg1a = legend(ag1,[h3,h1,h4,h2],'Vertical (Wahr)','Vertical (ISSM)','Horizontal (Wahr)','Horizontal (ISSM)',1); ++ 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'); +@@ -195,5 +197,5 @@ + + %export_fig('Fig6.pdf'); + +-end ++end % }}} + Index: /issm/oecreview/Archive/22755-22818/ISSM-22811-22812.diff =================================================================== --- /issm/oecreview/Archive/22755-22818/ISSM-22811-22812.diff (revision 22819) +++ /issm/oecreview/Archive/22755-22818/ISSM-22811-22812.diff (revision 22819) @@ -0,0 +1,59 @@ +Index: ../trunk-jpl/examples/SlrGRACE/runme.m +=================================================================== +--- ../trunk-jpl/examples/SlrGRACE/runme.m (revision 22811) ++++ ../trunk-jpl/examples/SlrGRACE/runme.m (revision 22812) +@@ -2,7 +2,7 @@ + clear all; + addpath('../Data','../Functions'); + +-steps=[7]; % [0:7]; ++steps=[0]; % [0:7]; + + if any(steps==0) % Download GRACE land_mass data {{{ + disp(' Step 0: Download GRACE land_mass data'); +@@ -177,7 +177,9 @@ + % loads and solutions. + sol1 = md.slr.deltathickness*100; % WEH 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; % degree + +@@ -226,7 +228,7 @@ + % define colormap, caxis, xlim etc {{{ + c1=colorbar; + colormap('haxby'); +- caxis([-floor(min(abs(min(sol)),abs(max(sol)))) floor(min(abs(min(sol)),abs(max(sol))))]) ++ caxis([-floor(min(abs(min(sol)),abs(max(sol)))) floor(min(abs(min(sol)),abs(max(sol))))]); + xlim([-180 180]); + ylim([-90 90]); + % }}} +@@ -234,7 +236,7 @@ + title(sol_name(kk)); + set(gcf,'color','w'); + +- %export_fig('Fig5.pdf'); ++ %export_fig(fig_name{kk}); + end + + end % }}} +@@ -293,7 +295,7 @@ + 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_loads.avi','movie_fingerprints.avi'}; ++ movie_name = {'Movie_dH.avi','Movie_slr.avi'}; + + res = 1.0; % degree + +@@ -347,7 +349,7 @@ + % define colormap, caxis, xlim etc {{{ + 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)))))]) ++ 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]); + % }}} Index: /issm/oecreview/Archive/22755-22818/ISSM-22812-22813.diff =================================================================== --- /issm/oecreview/Archive/22755-22818/ISSM-22812-22813.diff (revision 22819) +++ /issm/oecreview/Archive/22755-22818/ISSM-22812-22813.diff (revision 22819) @@ -0,0 +1,154 @@ +Index: ../trunk-jpl/examples/SlrFarrell/runme.m +=================================================================== +--- ../trunk-jpl/examples/SlrFarrell/runme.m (revision 22812) ++++ ../trunk-jpl/examples/SlrFarrell/runme.m (revision 22813) +@@ -1,6 +1,6 @@ + + clear all; +-steps=[2]; % ++steps=[5]; % + + if any(steps==1) % Global mesh creation {{{ + disp(' Step 1: Global mesh creation'); +@@ -69,6 +69,9 @@ + + % initial sea-level: 1 m RSL everywhere. + md.slr.sealevel=md.mask.ocean_levelset; ++ ++ md.slr.deltathickness=zeros(md.mesh.numberofelements,1); ++ md.slr.steric_rate=zeros(md.mesh.numberofvertices,1); + + save ./Models/SlrFarrell.Loads md; + +@@ -84,16 +87,18 @@ + + % Love numbers and reference frame: CF or CM (choose one!) + nlove=10001; % up to 10,000 degree +- md.slr.love_h = love_numbers('h','CM'); md.slr.love_h(nlov+1:end)=[]; +- md.slr.love_k = love_numbers('k','CM'); md.slr.love_k(nlov+1:end)=[]; +- md.slr.love_l = love_numbers('l','CM'); md.slr.love_l(nlov+1:end)=[]; ++ 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)=[]; + + % Mask: for computational efficiency only those elements that have loads are convolved! +- md.mask.groundedice_levelset = ones(md.mesh.numberofvertices,1); % 1 = ice is grounnded +- md.mask.ice_levelset = ones(md.mesh.numberofvertices,1); +- pos=find(md.slr.deltathickness~=0); +- md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; % -1 = ice loads + md.mask.land_levelset = 1-md.mask.ocean_levelset; ++ % fake ice load in one element! ++ md.mask.ice_levelset = ones(md.mesh.numberofvertices,1); % no ice ++ md.mask.groundedice_levelset = -ones(md.mesh.numberofvertices,1); % floated... ++ pos=find(md.mesh.lat <-80); ++ md.mask.ice_levelset(pos(1))=-1; % ice yes! ++ md.mask.groundedice_levelset(pos(1))=1; % ice grounded! + + %% IGNORE BUT DO NOT DELETE %% {{{ + % Geometry: Important only when you want to couple with Ice Flow Model +@@ -118,9 +123,6 @@ + disp(' Step 4: Solve Slr solver'); + md = loadmodel('./Models/SlrFarrell.Parameterization'); + +- % Request outputs +- md.slr.requested_outputs = {'SlrUmotion','SlrNmotion','SlrEmotion'}; +- + % Cluster info + md.cluster=generic('name',oshostname(),'np',3); + md.verbose=verbose('111111111'); +@@ -136,14 +138,9 @@ + disp(' Step 5: Plot solutions'); + md = loadmodel('./Models/SlrFarrell.Solution'); + +- % loads and solutions. +- sol1 = md.slr.deltathickness*100; % WEH cm +- sol2 = md.results.SlrSolution.SlrUmotion*1000; % [mm] +- sol3 = md.results.SlrSolution.SlrNmotion*1000; % [mm] +- sol4 = md.results.SlrSolution.SlrEmotion*1000; % [mm] +- sol_name={'Change in water equivalent height [cm]', 'Vertical displacement [mm]',... +- 'Horizontal (NS) displacement [mm]', 'Horizontal (EW) displacement [mm]'}; +- ++ % solutions. ++ sol = md.results.SealevelriseSolution.Sealevel*100; % per cent normalized by GMSL (which 1 m) ++ + res = 1.0; % degree + + % Make a grid of lats and lons, based on the min and max of the original vectors +@@ -150,51 +147,35 @@ + [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)); ++ % Make a interpolation object ++ F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol); ++ F.Method = 'linear'; ++ F.ExtrapolationMethod = 'linear'; + +- % if data are on elements, map those on to the vertices {{{ +- if length(sol)==md.mesh.numberofelements +- % map on to the vertices +- 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 % }}} ++ % Do the interpolation to get gridded solutions... ++ sol_grid = F(lat_grid, lon_grid); ++ sol_grid(isnan(sol_grid))=0; ++ sol_grid(lat_grid>85 & sol_grid==0) =NaN; % set polar unphysical 0s to Nan + +- % Make a interpolation object +- F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol); +- F.Method = 'linear'; +- F.ExtrapolationMethod = 'linear'; ++ 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; ++ contourf(lon_grid,lat_grid,sol_grid,[96 98 100 102 104 105]); shading flat; hold on; ++ geoshow(lat,long,'DisplayType','polygon','FaceColor','white'); ++ plot(long,lat,'k'); hold off; ++ % define colormap, caxis, xlim etc {{{ ++ c1=colorbar; ++ %colormap('haxby'); ++ %caxis([96 104]); ++ xlim([-180 180]); ++ ylim([-90 90]); ++ % }}} ++ grid on; ++ title('Relative sea-level [mm]'); ++ set(gcf,'color','w'); + +- % Do the interpolation to get gridded solutions... +- sol_grid = F(lat_grid, lon_grid); +- sol_grid(isnan(sol_grid))=0; +- sol_grid(lat_grid>85 & sol_grid==0) =NaN; % set polar unphysical 0s to Nan ++ %export_fig('Fig5.pdf'); + +- 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; +- plot(long,lat,'k'); hold off; +- c1=colorbar; +- colormap(jet); +- xlim([-180 180]); +- ylim([-90 90]); +- grid on; +- title(sol_name(kk)); +- set(gcf,'color','w'); +- +- %export_fig('Fig5.pdf'); +- end +- + end % }}} + +- Index: /issm/oecreview/Archive/22755-22818/ISSM-22813-22814.diff =================================================================== --- /issm/oecreview/Archive/22755-22818/ISSM-22813-22814.diff (revision 22819) +++ /issm/oecreview/Archive/22755-22818/ISSM-22813-22814.diff (revision 22819) @@ -0,0 +1,15 @@ +Index: ../trunk-jpl/examples/SlrFarrell/Farrell_1976_Fig1.pdf +=================================================================== +Cannot display: file marked as a binary type. +svn:mime-type = application/octet-stream +Index: ../trunk-jpl/examples/SlrFarrell/Farrell_1976_Fig1.pdf +=================================================================== +--- ../trunk-jpl/examples/SlrFarrell/Farrell_1976_Fig1.pdf (nonexistent) ++++ ../trunk-jpl/examples/SlrFarrell/Farrell_1976_Fig1.pdf (revision 22814) + +Property changes on: ../trunk-jpl/examples/SlrFarrell/Farrell_1976_Fig1.pdf +___________________________________________________________________ +Added: svn:mime-type +## -0,0 +1 ## ++application/octet-stream +\ No newline at end of property Index: /issm/oecreview/Archive/22755-22818/ISSM-22814-22815.diff =================================================================== --- /issm/oecreview/Archive/22755-22818/ISSM-22814-22815.diff (revision 22819) +++ /issm/oecreview/Archive/22755-22818/ISSM-22814-22815.diff (revision 22819) @@ -0,0 +1,156 @@ +Index: ../trunk-jpl/examples/EsaWahr/runme.m +=================================================================== +--- ../trunk-jpl/examples/EsaWahr/runme.m (revision 22814) ++++ ../trunk-jpl/examples/EsaWahr/runme.m (revision 22815) +@@ -2,7 +2,7 @@ + clear all; + addpath('../Functions'); + +-steps=[0:6]; % [0:6]; ++steps=[0]; % [0:6]; + + if any(steps==0) % Simple mesh creation {{{ + disp(' Step 0: Mesh creation'); +Index: ../trunk-jpl/examples/SlrFarrell/runme.m +=================================================================== +--- ../trunk-jpl/examples/SlrFarrell/runme.m (revision 22814) ++++ ../trunk-jpl/examples/SlrFarrell/runme.m (revision 22815) +@@ -1,7 +1,8 @@ + + clear all; +-steps=[5]; % + ++steps=[1]; % [1:5]; ++ + if any(steps==1) % Global mesh creation {{{ + disp(' Step 1: Global mesh creation'); + +@@ -127,6 +128,9 @@ + md.cluster=generic('name',oshostname(),'np',3); + md.verbose=verbose('111111111'); + ++ % Choose different convergence threshold. [10% 1% 0.1%] to match Farrell 3 panels in Fig. 1 ++ md.slr.reltol = 0.1/100; % per cent change in solution ++ + % Solve + md=solve(md,'Slr'); + +@@ -140,9 +144,9 @@ + + % solutions. + sol = md.results.SealevelriseSolution.Sealevel*100; % per cent normalized by GMSL (which 1 m) +- +- res = 1.0; % degree + ++ res = 1; % degree ++ + % Make a grid of lats and lons, based on the min and max of the original vectors + [lat_grid, lon_grid] = meshgrid(linspace(-90,90,180/res), linspace(-180,180,360/res)); + sol_grid = zeros(size(lat_grid)); +@@ -149,8 +153,8 @@ + + % Make a interpolation object + F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol); +- F.Method = 'linear'; +- F.ExtrapolationMethod = 'linear'; ++ F.Method = 'natural'; % for smooth contour ++ F.ExtrapolationMethod = 'none'; + + % Do the interpolation to get gridded solutions... + sol_grid = F(lat_grid, lon_grid); +@@ -160,19 +164,20 @@ + 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; +- contourf(lon_grid,lat_grid,sol_grid,[96 98 100 102 104 105]); shading flat; hold on; +- geoshow(lat,long,'DisplayType','polygon','FaceColor','white'); ++ 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; + % define colormap, caxis, xlim etc {{{ + c1=colorbar; +- %colormap('haxby'); +- %caxis([96 104]); +- xlim([-180 180]); +- ylim([-90 90]); ++ colormap(flipud(haxby)); ++ caxis([96 105]); ++ xlim([-170 170]); ++ ylim([-85 85]); + % }}} + grid on; +- title('Relative sea-level [mm]'); ++ title('Relative sea-level [% of GMSL]'); + set(gcf,'color','w'); + + %export_fig('Fig5.pdf'); +Index: ../trunk-jpl/examples/EsaGRACE/runme.m +=================================================================== +--- ../trunk-jpl/examples/EsaGRACE/runme.m (revision 22814) ++++ ../trunk-jpl/examples/EsaGRACE/runme.m (revision 22815) +@@ -1,7 +1,9 @@ + + clear all; +-steps=[0:5]; % [1:5]; ++addpath('../Data','../Functions'); + ++steps=[0]; % [0:5]; ++ + if any(steps==0) % Download GRACE land_mass data {{{ + disp(' Step 0: Download GRACE land_mass data'); + +@@ -11,6 +13,8 @@ + cd(f,'allData/tellus/L3/land_mass/RL05/netcdf/'); + mget(f,'GRCTellus.JPL.200204_201701.LND.RL05_1.DSTvSCS1411.nc'); + ++ !mv *.nc '../Data/' ++ + % display the content of the netcdf file. + %ncdisp('GRCTellus.JPL.200204_201701.LND.RL05_1.DSTvSCS1411.nc') + +@@ -121,8 +125,10 @@ + 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 + +@@ -159,20 +165,25 @@ + + set(0,'DefaultAxesFontSize',18,'DefaultAxesLineWidth',1,'DefaultTextFontSize',18,'DefaultLineMarkerSize',8) + figure1=figure('Position', [100, 100, 1000, 500]); +- gcf; +- load coast; +- cla; ++ gcf; load coast; cla; + pcolor(lon_grid,lat_grid,sol_grid); shading flat; hold on; ++ % mask out oceans {{{ ++ if (kk==1) ++ geoshow(flipud(lat),flipud(long),'DisplayType','polygon','FaceColor','white'); ++ end % }}} + plot(long,lat,'k'); hold off; ++ % colormap, xlim etc {{{ + c1=colorbar; +- colormap(jet); ++ 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('Fig5.pdf'); ++ ++ %export_fig(fig_name{kk}); + end + + end % }}} Index: /issm/oecreview/Archive/22755-22818/ISSM-22815-22816.diff =================================================================== --- /issm/oecreview/Archive/22755-22818/ISSM-22815-22816.diff (revision 22819) +++ /issm/oecreview/Archive/22755-22818/ISSM-22815-22816.diff (revision 22819) @@ -0,0 +1,31 @@ +Index: ../trunk-jpl/src/c/cores/transient_core.cpp +=================================================================== +--- ../trunk-jpl/src/c/cores/transient_core.cpp (revision 22815) ++++ ../trunk-jpl/src/c/cores/transient_core.cpp (revision 22816) +@@ -20,7 +20,7 @@ + + /*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; + int output_frequency; +@@ -53,6 +53,7 @@ + femmodel->parameters->FindParam(&issmb,TransientIssmbEnum); + 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); + femmodel->parameters->FindParam(&isgroundingline,TransientIsgroundinglineEnum); +@@ -310,6 +311,9 @@ + #endif + } + ++ /*esa: */ ++ if(isesa) esa_core(femmodel); ++ + /*Sea level rise: */ + if(isslr || iscoupler) sealevelrise_core(femmodel); + Index: /issm/oecreview/Archive/22755-22818/ISSM-22816-22817.diff =================================================================== --- /issm/oecreview/Archive/22755-22818/ISSM-22816-22817.diff (revision 22819) +++ /issm/oecreview/Archive/22755-22818/ISSM-22816-22817.diff (revision 22819) @@ -0,0 +1,26 @@ +Index: ../trunk-jpl/src/c/cores/sealevelrise_core.cpp +=================================================================== +--- ../trunk-jpl/src/c/cores/sealevelrise_core.cpp (revision 22816) ++++ ../trunk-jpl/src/c/cores/sealevelrise_core.cpp (revision 22817) +@@ -14,6 +14,7 @@ + Vector *Sg = NULL; + Vector *Sg_absolute = NULL; + Vector *Sg_eustatic = NULL; ++ Vector *slr_initial = NULL; + Vector *steric_rate_g = NULL; + Vector *U_radial = NULL; + Vector *U_north = NULL; +@@ -79,6 +80,13 @@ + /*call sea-level rise sub cores:*/ + 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/oecreview/Archive/22755-22818/ISSM-22817-22818.diff =================================================================== --- /issm/oecreview/Archive/22755-22818/ISSM-22817-22818.diff (revision 22819) +++ /issm/oecreview/Archive/22755-22818/ISSM-22817-22818.diff (revision 22819) @@ -0,0 +1,1116 @@ +Index: ../trunk-jpl/examples/EsaWahr/runme.m +=================================================================== +--- ../trunk-jpl/examples/EsaWahr/runme.m (revision 22817) ++++ ../trunk-jpl/examples/EsaWahr/runme.m (revision 22818) +@@ -2,129 +2,108 @@ + clear all; + addpath('../Functions'); + +-steps=[0]; % [0:6]; ++steps=[0]; % [0:6] + +-if any(steps==0) % Simple mesh creation {{{ ++if any(steps==0) + disp(' Step 0: Mesh creation'); + +- % Generate initial uniform mesh +- md=roundmesh(model,100000,10000); % Domain radius and element size (meters) ++ md=roundmesh(model,100000,10000); % Domain radius and element size [m] + +- save ./Models/EsaWahr.Mesh md; ++ save ./Models/EsaWahr_Mesh md; + + plotmodel(md,'data','mesh'); +- %export_fig('Fig0.pdf'); + +-end % }}} ++end + +-if any(steps==1) % {{{ Anisotropic mesh creation ++if any(steps==1) + disp(' Step 1: Anisotropic mesh creation'); + +- % Generate initial uniform mesh +- md=roundmesh(model,100000,1000); % Domain radius and element size (meters) +- %plotmodel(md,'data','mesh'); ++ md=roundmesh(model,100000,1000); + +- % Define a fake field to have anisotropic meshing +- disc_radius=20*1000; % km => m +- rad_dist=sqrt(md.mesh.x.^2+md.mesh.y.^2); % radial distance [m] ++ disc_radius=20*1000; ++ rad_dist=sqrt(md.mesh.x.^2+md.mesh.y.^2); + field = abs(rad_dist-disc_radius); +- %plotmodel(md,'data',field); + +- md = bamg(md,'field',field,'err',50,'hmax',10000); % error in field in m ++ md = bamg(md,'field',field,'err',50,'hmax',10000); + +- save ./Models/EsaWahr.Mesh md; ++ save ./Models/EsaWahr_Mesh md; + + plotmodel (md,'data','mesh'); +- %export_fig('Fig1.pdf'); + +-end % }}} ++end + +-if any(steps==2) % Define loads {{{ ++if any(steps==2) + disp(' Step 2: Define loads'); +- md = loadmodel('./Models/EsaWahr.Mesh'); ++ md = loadmodel('./Models/EsaWahr_Mesh'); + +- % primer + rho_w_i=md.materials.rho_freshwater/md.materials.rho_ice; + +- % elemental centroids + index=md.mesh.elements; + x_cent=mean(md.mesh.x(index),2); + y_cent=mean(md.mesh.y(index),2); + +- % Define a 20-km radius disc and unload it by 1 meter water height equivalent uniformly + md.esa.deltathickness = zeros(md.mesh.numberofelements,1); +- disc_radius=20; % km +- rad_dist=sqrt(x_cent.^2+y_cent.^2)/1000; % radial distance in km +- md.esa.deltathickness(rad_dist<=disc_radius) = -1.0*rho_w_i; % 1 m water withdrawl ++ 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; ++ save ./Models/EsaWahr_Loads md; + + plotmodel (md,'data',md.esa.deltathickness,'title','Ice height equivalent [m]'); +- %export_fig('Fig2.pdf'); + +-end % }}} ++end + +-if any(steps==3) % Parameterization {{{ ++if any(steps==3) + disp(' Step 3: Parameterization'); +- md = loadmodel('./Models/EsaWahr.Loads'); ++ md = loadmodel('./Models/EsaWahr_Loads'); + +- % Love numbers and reference frame: CF or CM (choose one!) +- nlove=10001; % up to 10,000 degree ++ 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)=[]; + +- % Mask: for computational efficiency only those elements that have loads are convolved! +- md.mask.ice_levelset = -ones(md.mesh.numberofvertices,1); % -1 = ice loads +- md.mask.groundedice_levelset = ones(md.mesh.numberofvertices,1); % 1 = ice is grounnded ++ md.mask.ice_levelset = -ones(md.mesh.numberofvertices,1); ++ md.mask.groundedice_levelset = ones(md.mesh.numberofvertices,1); + +- %% IGNORE BUT DO NOT DELETE %% {{{ +- % Geometry: Important only when you want to couple with Ice Flow Model + 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; +- % Materials: ++ + 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); +- % Miscellaneous: ++ + md.miscellaneous.name='EsaWahr'; +- %% IGNORE BUT DO NOT DELETE %% }}} + +- save ./Models/EsaWahr.Parameterization md; ++ save ./Models/EsaWahr_Parameterization md; + +-end % }}} ++end + +-if any(steps==4) % Solve {{{ ++if any(steps==4) + disp(' Step 4: Solve Esa solver'); +- md = loadmodel('./Models/EsaWahr.Parameterization'); ++ md = loadmodel('./Models/EsaWahr_Parameterization'); + +- % Request outputs + md.esa.requested_outputs = {'EsaUmotion','EsaXmotion','EsaYmotion'}; + +- % Cluster info + md.cluster=generic('name',oshostname(),'np',3); + md.verbose=verbose('111111111'); + +- % Solve + md=solve(md,'Esa'); + +- save ./Models/EsaWahr.Solution md; ++ save ./Models/EsaWahr_Solution md; + +-end % }}} ++end + +-if any(steps==5) % Plot solutions {{{ ++if any(steps==5) + disp(' Step 5: Plot solutions'); +- md = loadmodel('./Models/EsaWahr.Solution'); ++ md = loadmodel('./Models/EsaWahr_Solution'); + +- % m => mm +- 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); ++ 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] + +- % plot + set(0,'DefaultAxesFontSize',24,'DefaultAxesLineWidth',1,'DefaultTextFontSize',24,'DefaultLineMarkerSize',6); + figure('Position', [100, 100, 800, 600]); + plotmodel(md,'data',vert,... +@@ -146,32 +125,28 @@ + '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 % }}} ++end + +-if any(steps==6) % Compare results against semi-analytic solutions {{{ ++if any(steps==6) + disp(' Step 6: Compare results against Wahr semi-analytic solutions'); +- md = loadmodel('./Models/EsaWahr.Solution'); ++ md = loadmodel('./Models/EsaWahr_Solution'); + +- % numerical solutions +- % m => mm +- 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); +- % grid data from the center +- xi=[0:500:100000]; % 500 m resolution ++ 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; % 20 km ++ disc_radius = 20*1000; % [m] + [vert_wahr, horz_wahr] = wahr(disc_radius,xi,md.esa.love_h,md.esa.love_l); + +- % plot + set(0,'DefaultAxesFontSize',16,'DefaultAxesLineWidth',1,'DefaultTextFontSize',16,'DefaultLineMarkerSize',6); + figure1=figure('Position', [100, 100, 700, 400]); + ylabel_1={'0',' ','1','','2','','3',''}; +@@ -185,17 +160,14 @@ + % 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 ++ % 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); +- % first box + 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'); +- ++ set(gcf,'color','w'); + %export_fig('Fig6.pdf'); + +-end % }}} ++end + +Index: ../trunk-jpl/examples/SlrGRACE/runme.m +=================================================================== +--- ../trunk-jpl/examples/SlrGRACE/runme.m (revision 22817) ++++ ../trunk-jpl/examples/SlrGRACE/runme.m (revision 22818) +@@ -2,45 +2,26 @@ + clear all; + addpath('../Data','../Functions'); + +-steps=[0]; % [0:7]; ++steps=[1]; % [1:7] + +-if any(steps==0) % Download GRACE land_mass data {{{ +- disp(' Step 0: Download GRACE land_mass data'); +- +- % Connect to ftp server and download +- f = ftp('podaac-ftp.jpl.nasa.gov'); +- %dir(f) +- cd(f,'allData/tellus/L3/land_mass/RL05/netcdf/'); +- mget(f,'GRCTellus.JPL.200204_201701.LND.RL05_1.DSTvSCS1411.nc'); +- +- !mv *.nc '../Data/' +- +- % display the content of the netcdf file. +- %ncdisp('GRCTellus.JPL.200204_201701.LND.RL05_1.DSTvSCS1411.nc') +- +-end % }}} +- +-if any(steps==1) % Global mesh creation {{{ ++if any(steps==1) + disp(' Step 1: Global mesh creation'); + + numrefine=1; +- resolution=150*1e3; % inital resolution [m]. It determines, e.g., whether we capture small islands. +- 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] ++ 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] + +- %mesh earth: + md=model; +- md.mask=maskpsl(); % use maskpsl class (instead of mask) to store the ocean function as a ocean_levelset +- md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3); % attributes should be in km. ++ md.mask=maskpsl(); ++ md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3); % attributes in [km] + +- for i=1:numrefine; % refine mesh... {{{ ++ for i=1:numrefine, + +- %figure out mask: + md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); + +- %figure out distance to the coastline, in lat,long (not x,y,z): + distance=zeros(md.mesh.numberofvertices,1); + + pos=find(~md.mask.ocean_levelset); coaste.lat=md.mesh.lat(pos); coaste.long=md.mesh.long(pos); +@@ -47,7 +28,6 @@ + pos=find(md.mask.ocean_levelset); coasto.lat=md.mesh.lat(pos); coasto.long=md.mesh.long(pos); + + for j=1:md.mesh.numberofvertices +- %figure out nearest coastline (using the great circle distance) + 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; +@@ -62,128 +42,103 @@ + end + pos=find(distancemindistance_land); + distance(pos2)=mindistance_land; + +- dist=min(maxdistance,distance); % max size 1000 km +- %use distance to the coastline to refine mesh: ++ dist=min(maxdistance,distance); + md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3,'refine',md.mesh,'refinemetric',dist); +- end % }}} + +- %figure out mask: ++ end ++ + md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); + +- save ./Models/SlrGRACE.Mesh md; ++ save ./Models/SlrGRACE_Mesh md; + + plotmodel (md,'data',md.mask.ocean_levelset,'edgecolor','k'); +- %export_fig('Fig1.pdf'); + +-end % }}} ++end + +-if any(steps==2) % Define loads {{{ ++if any(steps==2) + disp(' Step 2: Define loads in meters of ice height equivalent'); +- md = loadmodel('./Models/SlrGRACE.Mesh'); ++ md = loadmodel('./Models/SlrGRACE_Mesh'); + +- % define time interval for analysis + year_month = 2007+15/365; + time_range = [year_month year_month]; + +- % map GRACE water load on to the mesh for the seleted month(s) + 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; % ice height equivalent ++ md.slr.deltathickness = water_load*md.materials.rho_freshwater/md.materials.rho_ice; + +- save ./Models/SlrGRACE.Loads md; ++ save ./Models/SlrGRACE_Loads md; + +- plotmodel (md,'data',md.slr.deltathickness,... +- 'view',[90 -90],'caxis',[-.1 .1],... +- 'title','Ice height equivalent [m]'); +- %export_fig('Fig2.pdf'); ++ plotmodel (md,'data',md.slr.deltathickness,'view',[90 -90],'caxis',[-.1 .1],'title','Ice height equivalent [m]'); + +-end % }}} ++end + +-if any(steps==3) % Parameterization {{{ ++if any(steps==3) + disp(' Step 3: Parameterization'); +- md = loadmodel('./Models/SlrGRACE.Loads'); ++ md = loadmodel('./Models/SlrGRACE_Loads'); + +- % Love numbers and reference frame: CF or CM (choose one!) +- nlove=10001; % up to 10,000 degree ++ 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)=[]; + +- % Mask: for computational efficiency only those elements that have loads are convolved! +- md.mask.groundedice_levelset = ones(md.mesh.numberofvertices,1); % 1 = ice is grounnded ++ 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; % -1 = ice loads ++ md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; + md.mask.land_levelset = 1-md.mask.ocean_levelset; + +- % initialize + md.slr.sealevel=zeros(md.mesh.numberofvertices,1); +- % steric rate + md.slr.steric_rate=zeros(md.mesh.numberofvertices,1); +- % solution is scaled according to "exact" area of the oceans + md.slr.ocean_area_scaling = 1; + +- % convergence criteria +- %md.slr.reltol=NaN; +- %md.slr.abstol=1e-4; +- +- %% IGNORE BUT DO NOT DELETE %% {{{ +- % Geometry: Important only when you want to couple with Ice Flow Model + 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; +- % Materials: ++ + 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); +- % Miscellaneous: ++ + md.miscellaneous.name='SlrGRACE'; +- %% IGNORE BUT DO NOT DELETE %% }}} + +- save ./Models/SlrGRACE.Parameterization md; ++ save ./Models/SlrGRACE_Parameterization md; + +-end % }}} ++end + +-if any(steps==4) % Solve {{{ ++if any(steps==4) + disp(' Step 4: Solve Slr solver'); +- md = loadmodel('./Models/SlrGRACE.Parameterization'); ++ md = loadmodel('./Models/SlrGRACE_Parameterization'); + +- % What physics to capture? + md.slr.rigid=1; + md.slr.elastic=1; + md.slr.rotation=1; + +- % Cluster info + md.cluster=generic('name',oshostname(),'np',3); + md.verbose=verbose('111111111'); + +- % Solve + md=solve(md,'Slr'); + +- save ./Models/SlrGRACE.Solution md; ++ save ./Models/SlrGRACE_Solution md; + +-end % }}} ++end + +-if any(steps==5) % Plot solutions {{{ ++if any(steps==5) + disp(' Step 5: Plot solutions'); +- md = loadmodel('./Models/SlrGRACE.Solution'); ++ md = loadmodel('./Models/SlrGRACE_Solution'); + +- % loads and solutions. +- sol1 = md.slr.deltathickness*100; % WEH cm ++ 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; % degree ++ res = 1.0; + +- % Make a grid of lats and lons, based on the min and max of the original vectors + [lat_grid, lon_grid] = meshgrid(linspace(-90,90,180/res), linspace(-180,180,360/res)); + sol_grid = zeros(size(lat_grid)); + +@@ -190,9 +145,7 @@ + for kk=1:2 + sol=eval(sprintf('sol%d',kk)); + +- % if data are on elements, map those on to the vertices {{{ + if length(sol)==md.mesh.numberofelements +- % map on to the vertices + for jj=1:md.mesh.numberofelements + ii=(jj-1)*3; + pp(ii+1:ii+3)=md.mesh.elements(jj,:); +@@ -202,61 +155,52 @@ + temp(jj)=mean(sol(pos)); + end + sol=temp'; +- end % }}} ++ end + +- % Make a interpolation object + F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol); + F.Method = 'linear'; + F.ExtrapolationMethod = 'linear'; + +- % Do the interpolation to get gridded solutions... + sol_grid = F(lat_grid, lon_grid); + sol_grid(isnan(sol_grid))=0; +- sol_grid(lat_grid>85 & sol_grid==0) =NaN; % set polar unphysical 0s to Nan ++ 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; +- % mask out land or oceans {{{ + if (kk==1) + geoshow(flipud(lat),flipud(long),'DisplayType','polygon','FaceColor','white'); + else + geoshow(lat,long,'DisplayType','polygon','FaceColor','white'); +- end % }}} ++ end + plot(long,lat,'k'); hold off; +- % define colormap, caxis, xlim etc {{{ + 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 % }}} ++end + +-if any(steps==6) % Transient {{{ ++if any(steps==6) + disp(' Step 6: Transient run'); +- md = loadmodel('./Models/SlrGRACE.Parameterization'); ++ md = loadmodel('./Models/SlrGRACE_Parameterization'); + +- % read GRACE data during the chosen time period << steps=2 >> + 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)); + +- % define ice load + [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; % ice height equivalent ++ md.slr.deltathickness(1:ne,:) = water_load*md.materials.rho_freshwater/md.materials.rho_ice; + md.slr.deltathickness(ne+1,:)=[1:nt]; % times + +- % define transient run + md.transient.issmb=0; + md.transient.ismasstransport=0; + md.transient.isstressbalance=0; +@@ -263,7 +207,6 @@ + md.transient.isthermal=0; + md.transient.isslr=1; + +- % time stepping, output frequency etc. + md.timestepping.start_time=-1; + md.timestepping.final_time=nt; + md.timestepping.time_step=1; +@@ -270,36 +213,31 @@ + md.timestepping.interp_forcings=0; + md.settings.output_frequency=1; + +- % Cluster info + md.cluster=generic('name',oshostname(),'np',3); + md.verbose=verbose('111111111'); + +- % Solve + md=solve(md,'tr'); + +- save ./Models/SlrGRACE.Transient md; ++ save ./Models/SlrGRACE_Transient md; + +-end % }}} ++end + +-if any(steps==7) % Plot transient {{{ ++if any(steps==7) + disp(' Step 7: Plot transient'); +- md = loadmodel('./Models/SlrGRACE.Transient'); ++ md = loadmodel('./Models/SlrGRACE_Transient'); + +- time = md.slr.deltathickness(end,:); % year ++ time = md.slr.deltathickness(end,:); + +- % loads and solutions. + 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] +- %% something weired happened here!!! first month solution is duplicated. Use offset of 1. Will fix it asap! + 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; % degree ++ res = 1.0; + +- % Make a grid of lats and lons, based on the min and max of the original vectors + [lat_grid, lon_grid] = meshgrid(linspace(-90,90,180/res), linspace(-180,180,360/res)); + sol_grid = zeros(size(lat_grid)); + +@@ -306,9 +244,7 @@ + for kk=1:2 + sol=eval(sprintf('sol%d',kk)); + +- % if data are on elements, map those on to the vertices {{{ + if length(sol)==md.mesh.numberofelements +- % map on to the vertices + for jj=1:md.mesh.numberofelements + ii=(jj-1)*3; + pp(ii+1:ii+3)=md.mesh.elements(jj,:); +@@ -318,7 +254,7 @@ + temp2(jj,:)=mean(sol(pos,:)); + end + sol=temp2; +- end % }}} ++ end + + vidObj = VideoWriter(movie_name{kk}); + vidObj.FrameRate=2; % frames per second +@@ -325,42 +261,35 @@ + open(vidObj); + + for jj=1:length(time) +- % Make a interpolation object + F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol(:,jj)); + F.Method = 'linear'; + F.ExtrapolationMethod = 'linear'; + +- % Do the interpolation to get gridded solutions... + sol_grid = F(lat_grid, lon_grid); + sol_grid(isnan(sol_grid))=0; +- sol_grid(lat_grid>85 & sol_grid==0) =NaN; % set polar unphysical 0s to Nan ++ 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; +- % mask out land or oceans {{{ + if (kk==1) + geoshow(flipud(lat),flipud(long),'DisplayType','polygon','FaceColor','white'); + else + geoshow(lat,long,'DisplayType','polygon','FaceColor','white'); +- end % }}} ++ end + plot(long,lat,'k'); hold off; +- % define colormap, caxis, xlim etc {{{ + 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 % close current figure ++ close + end +- +- % Close the file. + close(vidObj); + end + !open *.avi; +@@ -370,8 +299,7 @@ + xlabel('# month'); + ylabel('GMSL [mm/yr]'); + set(gcf,'color','w'); +- + %export_fig('Fig7.pdf'); + +-end % }}} ++end + +Index: ../trunk-jpl/examples/SlrFarrell/runme.m +=================================================================== +--- ../trunk-jpl/examples/SlrFarrell/runme.m (revision 22817) ++++ ../trunk-jpl/examples/SlrFarrell/runme.m (revision 22818) +@@ -1,29 +1,26 @@ + + clear all; + +-steps=[1]; % [1:5]; ++steps=[1]; % [1:5] + +-if any(steps==1) % Global mesh creation {{{ ++if any(steps==1) + disp(' Step 1: Global mesh creation'); + + numrefine=1; +- resolution=150*1e3; % inital resolution [m]. It determines, e.g., whether we capture small islands. +- 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] ++ 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] + +- %mesh earth: + md=model; +- md.mask=maskpsl(); % use maskpsl class (instead of mask) to store the ocean function as a ocean_levelset +- md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3); % attributes should be in km. ++ md.mask=maskpsl(); ++ md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3); % attributes in [km] + + for i=1:numrefine, + +- %figure out mask: + md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); + +- %figure out distance to the coastline, in lat,long (not x,y,z): + distance=zeros(md.mesh.numberofvertices,1); + + pos=find(~md.mask.ocean_levelset); coaste.lat=md.mesh.lat(pos); coaste.long=md.mesh.long(pos); +@@ -30,7 +27,6 @@ + pos=find(md.mask.ocean_levelset); coasto.lat=md.mesh.lat(pos); coasto.long=md.mesh.long(pos); + + for j=1:md.mesh.numberofvertices +- %figure out nearest coastline (using the great circle distance) + 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; +@@ -45,121 +41,102 @@ + end + pos=find(distancemindistance_land); + distance(pos2)=mindistance_land; + +- dist=min(maxdistance,distance); % max size 1000 km +- %use distance to the coastline to refine mesh: ++ dist=min(maxdistance,distance); + md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3,'refine',md.mesh,'refinemetric',dist); ++ + end + +- %figure out mask: + md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); + +- save ./Models/SlrFarrell.Mesh md; ++ save ./Models/SlrFarrell_Mesh md; + + plotmodel (md,'data',md.mask.ocean_levelset,'edgecolor','k'); +- %export_fig('Fig1.pdf'); + +-end % }}} ++end + +-if any(steps==2) % Define source {{{ ++if any(steps==2) + disp(' Step 2: Define source as in Farrell, 1972, Figure 1'); +- md = loadmodel('./Models/SlrFarrell.Mesh'); ++ md = loadmodel('./Models/SlrFarrell_Mesh'); + +- % initial sea-level: 1 m RSL everywhere. +- md.slr.sealevel=md.mask.ocean_levelset; ++ 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; ++ save ./Models/SlrFarrell_Loads md; + +- plotmodel (md,'data',md.slr.sealevel,'view',[90 90],... +- 'title#all','Initial sea-level [m]'); +- %export_fig('Fig2.pdf'); ++ plotmodel (md,'data',md.slr.sealevel,'view',[90 90],'title#all','Initial sea-level [m]'); + +-end % }}} ++end + +-if any(steps==3) % Parameterization {{{ ++if any(steps==3) + disp(' Step 3: Parameterization'); +- md = loadmodel('./Models/SlrFarrell.Loads'); ++ md = loadmodel('./Models/SlrFarrell_Loads'); + +- % Love numbers and reference frame: CF or CM (choose one!) +- nlove=10001; % up to 10,000 degree ++ 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)=[]; + +- % Mask: for computational efficiency only those elements that have loads are convolved! + md.mask.land_levelset = 1-md.mask.ocean_levelset; +- % fake ice load in one element! +- md.mask.ice_levelset = ones(md.mesh.numberofvertices,1); % no ice +- md.mask.groundedice_levelset = -ones(md.mesh.numberofvertices,1); % floated... ++ 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! + +- %% IGNORE BUT DO NOT DELETE %% {{{ +- % Geometry: Important only when you want to couple with Ice Flow Model + 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; +- % Materials: ++ + 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); +- % Miscellaneous: ++ + md.miscellaneous.name='SlrFarrell'; +- %% IGNORE BUT DO NOT DELETE %% }}} + +- save ./Models/SlrFarrell.Parameterization md; ++ save ./Models/SlrFarrell_Parameterization md; + +-end % }}} ++end + +-if any(steps==4) % Solve {{{ ++if any(steps==4) + disp(' Step 4: Solve Slr solver'); +- md = loadmodel('./Models/SlrFarrell.Parameterization'); ++ md = loadmodel('./Models/SlrFarrell_Parameterization'); + +- % Cluster info + md.cluster=generic('name',oshostname(),'np',3); + md.verbose=verbose('111111111'); + +- % Choose different convergence threshold. [10% 1% 0.1%] to match Farrell 3 panels in Fig. 1 + md.slr.reltol = 0.1/100; % per cent change in solution + +- % Solve + md=solve(md,'Slr'); + +- save ./Models/SlrFarrell.Solution md; ++ save ./Models/SlrFarrell_Solution md; + +-end % }}} ++end + +-if any(steps==5) % Plot solutions {{{ ++if any(steps==5) + disp(' Step 5: Plot solutions'); +- md = loadmodel('./Models/SlrFarrell.Solution'); ++ md = loadmodel('./Models/SlrFarrell_Solution'); + +- % solutions. + sol = md.results.SealevelriseSolution.Sealevel*100; % per cent normalized by GMSL (which 1 m) + +- res = 1; % degree ++ res = 1; % [degree] + +- % Make a grid of lats and lons, based on the min and max of the original vectors + [lat_grid, lon_grid] = meshgrid(linspace(-90,90,180/res), linspace(-180,180,360/res)); + sol_grid = zeros(size(lat_grid)); + +- % Make a interpolation object + F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol); +- F.Method = 'natural'; % for smooth contour +- F.ExtrapolationMethod = 'none'; ++ F.Method = 'linear'; ++ F.ExtrapolationMethod = 'linear'; + +- % Do the interpolation to get gridded solutions... + sol_grid = F(lat_grid, lon_grid); + sol_grid(isnan(sol_grid))=0; +- sol_grid(lat_grid>85 & sol_grid==0) =NaN; % set polar unphysical 0s to Nan ++ 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]); +@@ -169,18 +146,15 @@ + clabel(C,h,'FontSize',18,'Color','red','LabelSpacing',500); + geoshow(lat,long,'DisplayType','polygon','FaceColor',[.82 .82 .82]); + plot(long,lat,'k'); hold off; +- % define colormap, caxis, xlim etc {{{ + 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 % }}} ++end + +Index: ../trunk-jpl/examples/EsaGRACE/runme.m +=================================================================== +--- ../trunk-jpl/examples/EsaGRACE/runme.m (revision 22817) ++++ ../trunk-jpl/examples/EsaGRACE/runme.m (revision 22818) +@@ -2,137 +2,105 @@ + clear all; + addpath('../Data','../Functions'); + +-steps=[0]; % [0:5]; ++steps=[1]; % [1:5] + +-if any(steps==0) % Download GRACE land_mass data {{{ +- disp(' Step 0: Download GRACE land_mass data'); +- +- % Connect to ftp server and download +- f = ftp('podaac-ftp.jpl.nasa.gov'); +- %dir(f) +- cd(f,'allData/tellus/L3/land_mass/RL05/netcdf/'); +- mget(f,'GRCTellus.JPL.200204_201701.LND.RL05_1.DSTvSCS1411.nc'); +- +- !mv *.nc '../Data/' +- +- % display the content of the netcdf file. +- %ncdisp('GRCTellus.JPL.200204_201701.LND.RL05_1.DSTvSCS1411.nc') +- +-end % }}} +- +-if any(steps==1) % Global mesh creation {{{ ++if any(steps==1) + disp(' Step 1: Global mesh creation'); + +- resolution=300; % km. uniform meshing... +- radius = 6.371012*10^3; % mean radius of Earth, km ++ resolution=300; % [km] ++ radius = 6.371012*10^3; % [km] + +- %mesh earth: + md=model; +- md.mask=maskpsl(); % use maskpsl class (instead of mask) to store the ocean function as a ocean_levelset ++ md.mask=maskpsl(); + md.mesh=gmshplanet('radius',radius,'resolution',resolution); + +- %figure out mask: + md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); + +- save ./Models/EsaGRACE.Mesh md; ++ save ./Models/EsaGRACE_Mesh md; + + plotmodel (md,'data',md.mask.ocean_levelset,'edgecolor','k'); +- %export_fig('Fig1.pdf'); + +-end % }}} ++end + +-if any(steps==2) % Define loads {{{ ++if any(steps==2) + disp(' Step 2: Define loads in meters of ice height equivalent'); +- md = loadmodel('./Models/EsaGRACE.Mesh'); ++ md = loadmodel('./Models/EsaGRACE_Mesh'); + +- % define time interval for analysis + year_month = 2007+15/365; + time_range = [year_month year_month]; + +- % map GRACE water load on to the mesh for the seleted month(s) + 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; ++ save ./Models/EsaGRACE_Loads md; + + plotmodel (md,'data',md.esa.deltathickness,... + 'view',[90 -90],'caxis',[-.1 .1],... + 'title','Ice height equivalent [m]'); +- %export_fig('Fig2.pdf'); + +-end % }}} ++end + +-if any(steps==3) % Parameterization {{{ ++if any(steps==3) + disp(' Step 3: Parameterization'); +- md = loadmodel('./Models/EsaGRACE.Loads'); ++ md = loadmodel('./Models/EsaGRACE_Loads'); + +- % Love numbers and reference frame: CF or CM (choose one!) +- nlove=10001; % up to 10,000 degree ++ 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)=[]; + +- % Mask: for computational efficiency only those elements that have loads are convolved! +- md.mask.groundedice_levelset = ones(md.mesh.numberofvertices,1); % 1 = ice is grounnded ++ 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; % -1 = ice loads ++ md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; + md.mask.land_levelset = 1-md.mask.ocean_levelset; + +- %% IGNORE BUT DO NOT DELETE %% {{{ +- % Geometry: Important only when you want to couple with Ice Flow Model + 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; +- % Materials: ++ + 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); +- % Miscellaneous: ++ + md.miscellaneous.name='EsaGRACE'; +- %% IGNORE BUT DO NOT DELETE %% }}} + +- save ./Models/EsaGRACE.Parameterization md; ++ save ./Models/EsaGRACE_Parameterization md; + +-end % }}} ++end + +-if any(steps==4) % Solve {{{ ++if any(steps==4) + disp(' Step 4: Solve Esa solver'); +- md = loadmodel('./Models/EsaGRACE.Parameterization'); ++ md = loadmodel('./Models/EsaGRACE_Parameterization'); + +- % Request outputs + md.esa.requested_outputs = {'EsaUmotion','EsaNmotion','EsaEmotion'}; + +- % Cluster info + md.cluster=generic('name',oshostname(),'np',3); + md.verbose=verbose('111111111'); + +- % Solve + md=solve(md,'Esa'); + +- save ./Models/EsaGRACE.Solution md; ++ save ./Models/EsaGRACE_Solution md; + +-end % }}} ++end + +-if any(steps==5) % Plot solutions {{{ ++if any(steps==5) + disp(' Step 5: Plot solutions'); +- md = loadmodel('./Models/EsaGRACE.Solution'); ++ md = loadmodel('./Models/EsaGRACE_Solution'); + +- % loads and solutions. +- sol1 = md.esa.deltathickness*100; % WEH cm +- sol2 = md.results.EsaSolution.EsaUmotion*1000; % [mm] +- sol3 = md.results.EsaSolution.EsaNmotion*1000; % [mm] +- sol4 = md.results.EsaSolution.EsaEmotion*1000; % [mm] ++ 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 ++ res = 1.0; % [degree] + +- % Make a grid of lats and lons, based on the min and max of the original vectors + [lat_grid, lon_grid] = meshgrid(linspace(-90,90,180/res), linspace(-180,180,360/res)); + sol_grid = zeros(size(lat_grid)); + +@@ -139,9 +107,7 @@ + for kk=1:4 + sol=eval(sprintf('sol%d',kk)); + +- % if data are on elements, map those on to the vertices {{{ + if length(sol)==md.mesh.numberofelements +- % map on to the vertices + for jj=1:md.mesh.numberofelements + ii=(jj-1)*3; + pp(ii+1:ii+3)=md.mesh.elements(jj,:); +@@ -151,40 +117,34 @@ + temp(jj)=mean(sol(pos)); + end + sol=temp'; +- end % }}} ++ end + +- % Make a interpolation object + F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol); + F.Method = 'linear'; + F.ExtrapolationMethod = 'linear'; + +- % Do the interpolation to get gridded solutions... + sol_grid = F(lat_grid, lon_grid); + sol_grid(isnan(sol_grid))=0; +- sol_grid(lat_grid>85 & sol_grid==0) =NaN; % set polar unphysical 0s to Nan ++ 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; +- % mask out oceans {{{ + if (kk==1) + geoshow(flipud(lat),flipud(long),'DisplayType','polygon','FaceColor','white'); +- end % }}} ++ end + plot(long,lat,'k'); hold off; +- % colormap, xlim etc {{{ + 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 % }}} ++end + Index: /issm/oecreview/Archive/22755-22818/ISSM-DocReview-22755-22818.tex =================================================================== --- /issm/oecreview/Archive/22755-22818/ISSM-DocReview-22755-22818.tex (revision 22819) +++ /issm/oecreview/Archive/22755-22818/ISSM-DocReview-22755-22818.tex (revision 22819) @@ -0,0 +1,67 @@ +\documentclass[]{report} % list options between brackets +\usepackage{graphicx} % list packages between braces + +% type user-defined commands here + +\begin{document} + +\title{JPL EXPORT ADMINISTRATION \\ DOCUMENT REVIEW RECORD} % type title between braces +\author{Tom Scavo} % type author(s) between braces +\date{October 27, 1995} % type date between braces +%\maketitle + + +\begin{center} +\begin{tabular}{ |c|c| } +\hline +JPL EXPORT ADMINISTRATION \\ +\textit{ DOCUMENT REVIEW RECORD} \\ +\hline +\end{tabular} +\end{center} + +\hfill Log \#: \underline{\input{LogNumber}} + +\vspace{1cm} +\noindent (Note: This form and process do not replace the procedures described in JPL Policy relating to review and approval of proposals and contractual +documents. This process is intended to document the review and coordination of requests to ascertain the export control ramifications relating to specific +documents. Export Administration signature does not convey authority to export or release the "Exporter of Record" from any export laws or regulations.)\\ + +\noindent \textbf{Program:} \underline{ISSM: Ice Sheet System Model} \\ \\ +\noindent \textbf{Person Requesting or Initiating Export}: \underline{Dr. Eric Larour}\\ \\ +\noindent \textbf{Date Received}: \underline{\input{Date}}\\ \\ +\noindent \textbf{Document Title/Description}: ISSM changes from revision \input{r1} to revision \input{r2} \\ \\ +\noindent \textbf{Release to:} \underline{http://issm.ess.uci.edu/svn/issm/issm/trunk on ISSM svn repository}\\ \\ +\noindent \textbf{JPL Intranet:} \underline{murdo.jpl.nasa.gov/proj/ice/larour/issm-uci/trunk-jpl}\\ \\ + +\noindent \textbf{Disposition: 6 } +Does not contain export-controlled information. May be released/disclosed as requested subject to +Company guidelines on protection of proprietary information (if applicable). \\ + +\noindent \textbf{Comments:} see table of changes below. \\ \\ +\noindent \textbf{Reviewed by ISSM Export Transfer Liaison:} Dr. Eric Larour \hfill \textbf{Date:} \input{Date} \\ +\includegraphics[scale=1]{signature} + +\noindent JPL Export Administration Form TBS – June 29, 2011 + +\begin{center} +\line(1,0){250} +\end{center} + +\noindent \textbf{Disposition:} \\ +1: Public Domain Information (Ref ITAR Section 120.11) \\ +2: Qualifies for ITAR Exemption \\ +3: Covered by Department of State License/Agreement Number \\ +4: Covered by Department of Commerce validated license or exception \\ +5: New License Required \\ +6: Does not contain export-controlled information. May be released/disclosed as requested subject + to Company guidelines on protection of proprietary information (if applicable). \\ +7: Other (specify) +\begin{center} +\line(1,0){250} +\end{center} + + +\input{log} + +\end{document} Index: /issm/oecreview/Archive/22755-22818/LogNumber.tex =================================================================== --- /issm/oecreview/Archive/22755-22818/LogNumber.tex (revision 22819) +++ /issm/oecreview/Archive/22755-22818/LogNumber.tex (revision 22819) @@ -0,0 +1,1 @@ +22755-22818 Index: /issm/oecreview/Archive/22755-22818/Makefile =================================================================== --- /issm/oecreview/Archive/22755-22818/Makefile (revision 22819) +++ /issm/oecreview/Archive/22755-22818/Makefile (revision 22819) @@ -0,0 +1,8 @@ +TARGET=ISSM-DocReview-22755-22818 + +all: + pdflatex -interaction=errorstopmode -file-line-error -halt-on-error $(TARGET).tex + rm -rf *.log *.aux + +clean: + rm -rf *.log *.aux Index: /issm/oecreview/Archive/22755-22818/log.tex =================================================================== --- /issm/oecreview/Archive/22755-22818/log.tex (revision 22819) +++ /issm/oecreview/Archive/22755-22818/log.tex (revision 22819) @@ -0,0 +1,301 @@ +\noindent \textbf{Change \#1} with diff file ISSM-22757-22758.diff: \\ +Function name: \\ +M /issm/trunk-jpl M /issm/trunk-jpl/externalpackages/mpich/install-3.0-linux64-static.sh M /issm/trunk-jpl/externalpackages/mpich/install-3.2-linux64-static.sh M /issm/trunk-jpl/externalpackages/petsc/install-3.6-linux64-static.sh M /issm/trunk-jpl/externalpackages/petsc/install-3.7-linux64-static.sh M /issm/trunk-jpl/externalpackages/petsc/install-3.7-macosx64-static.sh M /issm/trunk-jpl/jenkins/jenkins.sh M /issm/trunk-jpl/jenkins/macosx\_pine-island\_dakota\_static M /issm/trunk-jpl/jenkins/macosx\_pine-island\_static M /issm/trunk-jpl/packagers/macosx/package.sh M /issm/trunk-jpl/packagers/macosx-dakota/package.sh M /issm/trunk-jpl/packagers/ubuntu/package.sh M /issm/trunk-jpl/src\\ +Export determination: 6. \\ +Rationale: merged trunk and trunk-jpl\\ +\vspace{3em} + +\noindent \textbf{Change \#2} with diff file ISSM-22759-22760.diff: \\ +Function name: \\ +M /issm/trunk-jpl M /issm/trunk-jpl/src M /issm/trunk-jpl/test\\ +Export determination: 6. \\ +Rationale: Block revision 22758 from being merged into trunk-jpl\\ +\vspace{3em} + +\noindent \textbf{Change \#3} with diff file ISSM-22760-22761.diff: \\ +Function name: \\ +M /issm/trunk-jpl/README M /issm/trunk-jpl/configure.ac M /issm/trunk-jpl/src/m/dev/issmversion.m M /issm/trunk-jpl/src/m/dev/issmversion.py\\ +Export determination: 6. \\ +Rationale: CHG: new version number\\ +\vspace{3em} + +\noindent \textbf{Change \#4} with diff file ISSM-22761-22762.diff: \\ +Function name: \\ +M /issm/trunk-jpl/jenkins/jenkins.sh\\ +Export determination: 6. \\ +Rationale: CHG: matlab options were overkill\\ +\vspace{3em} + +\noindent \textbf{Change \#5} with diff file ISSM-22762-22763.diff: \\ +Function name: \\ +M /issm/trunk-jpl/jenkins/macosx\_pine-island\_static\\ +Export determination: 6. \\ +Rationale: CHG: now require to install math77 and gmsh\\ +\vspace{3em} + +\noindent \textbf{Change \#6} with diff file ISSM-22764-22765.diff: \\ +Function name: \\ +M /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp M /issm/trunk-jpl/src/c/classes/ExternalResults/GenericExternalResult.h M /issm/trunk-jpl/src/c/cores/controladm1qn3\_core.cpp M /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateOutputDefinitions.cpp\\ +Export determination: 6. \\ +Rationale: CHG: clean up code and enable transpose of ExternalResults\\ +\vspace{3em} + +\noindent \textbf{Change \#7} with diff file ISSM-22769-22770.diff: \\ +Function name: \\ +M /issm/trunk-jpl/src/m/classes/clusters/generic.m M /issm/trunk-jpl/src/m/classes/clusters/pfe.m M /issm/trunk-jpl/src/m/inversions/marshallcostfunctions.m M /issm/trunk-jpl/src/m/inversions/supportedcontrols.m A /issm/trunk-jpl/src/m/mesh/labelconnectedregions.m M /issm/trunk-jpl/src/m/mesh/planet/gmsh/gmshplanet.m M /issm/trunk-jpl/src/m/plot/subplotmodel.m M /issm/trunk-jpl/src/m/solvers/mumpsoptions.m\\ +Export determination: 6. \\ +Rationale: CHG: fixing call to gmsh for workshop binaries\\ +\vspace{3em} + +\noindent \textbf{Change \#8} with diff file ISSM-22770-22771.diff: \\ +Function name: \\ +M /issm/trunk-jpl/src/m/classes/clusters/generic.m M /issm/trunk-jpl/src/m/classes/clusters/pfe.m M /issm/trunk-jpl/src/m/inversions/marshallcostfunctions.m M /issm/trunk-jpl/src/m/inversions/supportedcontrols.m M /issm/trunk-jpl/src/m/plot/subplotmodel.m M /issm/trunk-jpl/src/m/solvers/mumpsoptions.m\\ +Export determination: 6. \\ +Rationale: CHG: reverting changes\\ +\vspace{3em} + +\noindent \textbf{Change \#9} with diff file ISSM-22777-22778.diff: \\ +Function name: \\ +M /issm/trunk-jpl/src/m/classes/clusters/generic.m\\ +Export determination: 6. \\ +Rationale: NEW: added valgind in coupled simulations\\ +\vspace{3em} + +\noindent \textbf{Change \#10} with diff file ISSM-22778-22779.diff: \\ +Function name: \\ +M /issm/trunk-jpl/src/c/classes/Vertex.cpp M /issm/trunk-jpl/src/c/classes/Vertex.h M /issm/trunk-jpl/src/c/classes/Vertices.cpp M /issm/trunk-jpl/src/c/classes/Vertices.h M /issm/trunk-jpl/src/c/cores/transient\_core.cpp M /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp\\ +Export determination: 6. \\ +Rationale: NEW: starting to interpolate data between ice and ocean grids\\ +\vspace{3em} + +\noindent \textbf{Change \#11} with diff file ISSM-22779-22780.diff: \\ +Function name: \\ +M /issm/trunk-jpl/src/c/classes/Vertices.cpp\\ +Export determination: 6. \\ +Rationale: FIX: forgot to delete variables\\ +\vspace{3em} + +\noindent \textbf{Change \#12} with diff file ISSM-22780-22781.diff: \\ +Function name: \\ +M /issm/trunk-jpl/src/c/classes/Vertices.cpp\\ +Export determination: 6. \\ +Rationale: FIX: no need to print lat/long\\ +\vspace{3em} + +\noindent \textbf{Change \#13} with diff file ISSM-22782-22783.diff: \\ +Function name: \\ +M /issm/trunk-jpl/src/c/cores/transient\_core.cpp\\ +Export determination: 6. \\ +Rationale: FIX: fixed problem with lat/long\\ +\vspace{3em} + +\noindent \textbf{Change \#14} with diff file ISSM-22784-22785.diff: \\ +Function name: \\ +M /issm/trunk-jpl/src/c/cores/transient\_core.cpp\\ +Export determination: 6. \\ +Rationale: CHG: fixed AD\\ +\vspace{3em} + +\noindent \textbf{Change \#15} with diff file ISSM-22786-22787.diff: \\ +Function name: \\ +M /issm/trunk-jpl/src/c/classes/Vertex.cpp\\ +Export determination: 6. \\ +Rationale: FIX: lat/long only if they exist\\ +\vspace{3em} + +\noindent \textbf{Change \#16} with diff file ISSM-22787-22788.diff: \\ +Function name: \\ +M /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp M /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h M /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp M /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp\\ +Export determination: 6. \\ +Rationale: CHG: no need to have lat/long enums\\ +\vspace{3em} + +\noindent \textbf{Change \#17} with diff file ISSM-22789-22790.diff: \\ +Function name: \\ +M /issm/trunk-jpl/src/c/classes/Vertex.cpp\\ +Export determination: 6. \\ +Rationale: CHG: added some asserts\\ +\vspace{3em} + +\noindent \textbf{Change \#18} with diff file ISSM-22791-22792.diff: \\ +Function name: \\ +A /issm/trunk-jpl/examples/EsaGRACE A /issm/trunk-jpl/examples/EsaGRACE/Models A /issm/trunk-jpl/examples/EsaWahr A /issm/trunk-jpl/examples/EsaWahr/Models A /issm/trunk-jpl/examples/EsaWahr/runme.m A /issm/trunk-jpl/examples/SlrFarrell A /issm/trunk-jpl/examples/SlrFarrell/Models A /issm/trunk-jpl/examples/SlrGRACE A /issm/trunk-jpl/examples/SlrGRACE/Models\\ +Export determination: 6. \\ +Rationale: NEW: four tutorials for Hawaii workshop\\ +\vspace{3em} + +\noindent \textbf{Change \#19} with diff file ISSM-22793-22794.diff: \\ +Function name: \\ +A /issm/trunk-jpl/examples/EsaWahr/RoundDomain.exp M /issm/trunk-jpl/examples/EsaWahr/runme.m A /issm/trunk-jpl/examples/EsaWahr/wahr.m\\ +Export determination: 6. \\ +Rationale: CHG: tutorial 1 completed\\ +\vspace{3em} + +\noindent \textbf{Change \#20} with diff file ISSM-22794-22795.diff: \\ +Function name: \\ +M /issm/trunk-jpl/src/c/classes/Elements/Element.cpp M /issm/trunk-jpl/src/c/classes/Inputs/ControlInput.cpp M /issm/trunk-jpl/src/c/cores/controladm1qn3\_core.cpp M /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp\\ +Export determination: 6. \\ +Rationale: CHG: added basal melt to IoCodeConversions, added better names to transient inputs for transient controls, fixed bug in controladm1qn3\\ +\vspace{3em} + +\noindent \textbf{Change \#21} with diff file ISSM-22795-22796.diff: \\ +Function name: \\ +M /issm/trunk-jpl/src/m/classes/geometry.js\\ +Export determination: 6. \\ +Rationale: CHG (JS): Fixing Matlab-JS semantic translation error.\\ +\vspace{3em} + +\noindent \textbf{Change \#22} with diff file ISSM-22796-22797.diff: \\ +Function name: \\ +A /issm/trunk-jpl/examples/EsaGRACE/grace.m A /issm/trunk-jpl/examples/EsaGRACE/runme.m\\ +Export determination: 6. \\ +Rationale: CHG: esa tutorial for GRACE loads\\ +\vspace{3em} + +\noindent \textbf{Change \#23} with diff file ISSM-22797-22798.diff: \\ +Function name: \\ +M /issm/trunk-jpl/test/MITgcm/code/cpl\_issm.F\\ +Export determination: 6. \\ +Rationale: MITgcm sends shelficeFreshWaterFlux instead of shelficeHeatFlux to ISSM\\ +\vspace{3em} + +\noindent \textbf{Change \#24} with diff file ISSM-22798-22799.diff: \\ +Function name: \\ +M /issm/trunk-jpl/src/c/cores/transient\_core.cpp M /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp M /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h M /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp M /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp\\ +Export determination: 6. \\ +Rationale: CHG: fixed ocean heat flux used in ISSM for ocean coupling\\ +\vspace{3em} + +\noindent \textbf{Change \#25} with diff file ISSM-22799-22800.diff: \\ +Function name: \\ +M /issm/trunk-jpl/src/c/cores/transient\_core.cpp\\ +Export determination: 6. \\ +Rationale: CHG: send the right base to the ocean model\\ +\vspace{3em} + +\noindent \textbf{Change \#26} with diff file ISSM-22800-22801.diff: \\ +Function name: \\ +M /issm/trunk-jpl/src/c/cores/transient\_core.cpp\\ +Export determination: 6. \\ +Rationale: BUG: fixed leaks in coupling\\ +\vspace{3em} + +\noindent \textbf{Change \#27} with diff file ISSM-22801-22802.diff: \\ +Function name: \\ +M /issm/trunk-jpl/src/c/classes/Elements/Element.cpp M /issm/trunk-jpl/src/c/cores/controladm1qn3\_core.cpp M /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp\\ +Export determination: 6. \\ +Rationale: CHG: fixing non transient adolc controls and adding MaterialsRheologyBbar to IoCodeConversions\\ +\vspace{3em} + +\noindent \textbf{Change \#28} with diff file ISSM-22802-22803.diff: \\ +Function name: \\ +M /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp\\ +Export determination: 6. \\ +Rationale: BUG: fixed problem with materials in solid earth tests\\ +\vspace{3em} + +\noindent \textbf{Change \#29} with diff file ISSM-22803-22804.diff: \\ +Function name: \\ +M /issm/trunk-jpl/examples/EsaGRACE/runme.m\\ +Export determination: 6. \\ +Rationale: CHG: added plotting scripts\\ +\vspace{3em} + +\noindent \textbf{Change \#30} with diff file ISSM-22804-22805.diff: \\ +Function name: \\ +A /issm/trunk-jpl/examples/SlrGRACE/runme.m\\ +Export determination: 6. \\ +Rationale: NEW: last tutorial on GRACE sea-level fingerprint\\ +\vspace{3em} + +\noindent \textbf{Change \#31} with diff file ISSM-22805-22806.diff: \\ +Function name: \\ +A /issm/trunk-jpl/examples/SlrFarrell/runme.m\\ +Export determination: 6. \\ +Rationale: CHG: third tutorial place holder\\ +\vspace{3em} + +\noindent \textbf{Change \#32} with diff file ISSM-22806-22807.diff: \\ +Function name: \\ +M /issm/trunk-jpl/src/m/mesh/planet/gmsh/gmshplanet.m M /issm/trunk-jpl/src/m/mesh/planet/gmsh/gmshplanet.py\\ +Export determination: 6. \\ +Rationale: CHG: ruled surface is changed to surface\\ +\vspace{3em} + +\noindent \textbf{Change \#33} with diff file ISSM-22807-22808.diff: \\ +Function name: \\ +M /issm/trunk-jpl/src/m/mesh/planet/gmsh/gmshplanet.m\\ +Export determination: 6. \\ +Rationale: minor\\ +\vspace{3em} + +\noindent \textbf{Change \#34} with diff file ISSM-22808-22809.diff: \\ +Function name: \\ +M /issm/trunk-jpl/src/m/classes/slr.js M /issm/trunk-jpl/src/m/classes/slr.m M /issm/trunk-jpl/src/m/classes/slr.py\\ +Export determination: 6. \\ +Rationale: minor\\ +\vspace{3em} + +\noindent \textbf{Change \#35} with diff file ISSM-22809-22810.diff: \\ +Function name: \\ +M /issm/trunk-jpl/examples/EsaGRACE/runme.m A /issm/trunk-jpl/examples/Functions A /issm/trunk-jpl/examples/Functions/grace.m A /issm/trunk-jpl/examples/Functions/wahr.m M /issm/trunk-jpl/examples/SlrGRACE/runme.m\\ +Export determination: 6. \\ +Rationale: CHG: functions are located in a seperate directory etc\\ +\vspace{3em} + +\noindent \textbf{Change \#36} with diff file ISSM-22810-22811.diff: \\ +Function name: \\ +D /issm/trunk-jpl/examples/EsaGRACE/grace.m D /issm/trunk-jpl/examples/EsaWahr/wahr.m\\ +Export determination: 6. \\ +Rationale: CHG: duplicate functions removed from svn\\ +\vspace{3em} + +\noindent \textbf{Change \#37} with diff file ISSM-22811-22812.diff: \\ +Function name: \\ +M /issm/trunk-jpl/examples/EsaWahr/runme.m\\ +Export determination: 6. \\ +Rationale: minor\\ +\vspace{3em} + +\noindent \textbf{Change \#38} with diff file ISSM-22812-22813.diff: \\ +Function name: \\ +M /issm/trunk-jpl/examples/SlrGRACE/runme.m\\ +Export determination: 6. \\ +Rationale: CHG: minor\\ +\vspace{3em} + +\noindent \textbf{Change \#39} with diff file ISSM-22813-22814.diff: \\ +Function name: \\ +M /issm/trunk-jpl/examples/SlrFarrell/runme.m\\ +Export determination: 6. \\ +Rationale: CHG: Farrell solutions recovered\\ +\vspace{3em} + +\noindent \textbf{Change \#40} with diff file ISSM-22814-22815.diff: \\ +Function name: \\ +A /issm/trunk-jpl/examples/SlrFarrell/Farrell\_1976\_Fig1.pdf\\ +Export determination: 6. \\ +Rationale: NEW: Farrell solutions for benchmarking\\ +\vspace{3em} + +\noindent \textbf{Change \#41} with diff file ISSM-22815-22816.diff: \\ +Function name: \\ +M /issm/trunk-jpl/examples/EsaGRACE/runme.m M /issm/trunk-jpl/examples/EsaWahr/runme.m M /issm/trunk-jpl/examples/SlrFarrell/runme.m\\ +Export determination: 6. \\ +Rationale: CHG: final minor changes to all tutorials\\ +\vspace{3em} + +\noindent \textbf{Change \#42} with diff file ISSM-22816-22817.diff: \\ +Function name: \\ +M /issm/trunk-jpl/src/c/cores/transient\_core.cpp\\ +Export determination: 6. \\ +Rationale: CHG: esa added in transient core\\ +\vspace{3em} + +\noindent \textbf{Change \#43} with diff file ISSM-22817-22818.diff: \\ +Function name: \\ +M /issm/trunk-jpl/src/c/cores/sealevelrise\_core.cpp\\ +Export determination: 6. \\ +Rationale: CHG: minor adjustment to let intial slr relax according to SAL\\ +\vspace{3em} + Index: /issm/oecreview/Archive/22755-22818/r1.tex =================================================================== --- /issm/oecreview/Archive/22755-22818/r1.tex (revision 22819) +++ /issm/oecreview/Archive/22755-22818/r1.tex (revision 22819) @@ -0,0 +1,1 @@ +22755 Index: /issm/oecreview/Archive/22755-22818/r2.tex =================================================================== --- /issm/oecreview/Archive/22755-22818/r2.tex (revision 22819) +++ /issm/oecreview/Archive/22755-22818/r2.tex (revision 22819) @@ -0,0 +1,1 @@ +22818