Index: /issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp (revision 25750) +++ /issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp (revision 25751) @@ -177,4 +177,5 @@ bool rigid=false; bool elastic=false; + bool rotation=false; int numoutputs; @@ -203,4 +204,5 @@ parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.computesealevelchange",SolidearthSettingsComputesealevelchangeEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.planetradius",SolidearthPlanetRadiusEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.glfraction",SolidearthSettingsGlfractionEnum)); parameters->AddObject(new DoubleParam(CumBslrEnum,0.0)); parameters->AddObject(new DoubleParam(CumBslrIceEnum,0.0)); @@ -212,4 +214,27 @@ planetarea=4*PI*planetradius*planetradius; parameters->AddObject(new DoubleParam(SolidearthPlanetAreaEnum,planetarea)); + + /*Deal with partition of the barystatic contribution:*/ + iomodel->FetchData(&npartice,"md.solidearth.npartice"); + parameters->AddObject(new IntParam(SolidearthNpartIceEnum,npartice)); + if(npartice){ + iomodel->FetchData(&partitionice,&nel,NULL,"md.solidearth.partitionice"); + parameters->AddObject(new DoubleMatParam(SolidearthPartitionIceEnum,partitionice,nel,1)); + parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.npartice",SolidearthNpartIceEnum)); + bslrice_partition=xNewZeroInit(npartice); + parameters->AddObject(new DoubleMatParam(CumBslrIcePartitionEnum,bslrice_partition,npartice,1)); + xDelete(partitionice); + } + iomodel->FetchData(&nparthydro,"md.solidearth.nparthydro"); + parameters->AddObject(new IntParam(SolidearthNpartHydroEnum,nparthydro)); + if(nparthydro){ + iomodel->FetchData(&partitionhydro,&nel,NULL,"md.solidearth.partitionhydro"); + parameters->AddObject(new DoubleMatParam(SolidearthPartitionHydroEnum,partitionhydro,nel,1)); + parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.nparthydro",SolidearthNpartHydroEnum)); + bslrhydro_partition=xNewZeroInit(nparthydro); + parameters->AddObject(new DoubleMatParam(CumBslrHydroPartitionEnum,bslrhydro_partition,nparthydro,1)); + xDelete(partitionhydro); + } + /*Deal with dsl multi-model ensembles: {{{*/ @@ -234,134 +259,168 @@ iomodel->FetchData(&rigid,"md.solidearth.settings.rigid"); iomodel->FetchData(&elastic,"md.solidearth.settings.elastic"); - if (rigid) { - if (elastic) { - /*love numbers: */ - iomodel->FetchData(&love_h,&nl,NULL,"md.solidearth.lovenumbers.h"); - iomodel->FetchData(&love_k,&nl,NULL,"md.solidearth.lovenumbers.k"); - iomodel->FetchData(&love_l,&nl,NULL,"md.solidearth.lovenumbers.l"); - iomodel->FetchData(&love_th,&nl,NULL,"md.solidearth.lovenumbers.th"); - iomodel->FetchData(&love_tk,&nl,NULL,"md.solidearth.lovenumbers.tk"); - iomodel->FetchData(&love_tl,&nl,NULL,"md.solidearth.lovenumbers.tl"); - parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.lovenumbers.tk2secular",TidalLoveK2SecularEnum)); - - parameters->AddObject(new DoubleMatParam(LoadLoveHEnum,love_h,nl,1)); - parameters->AddObject(new DoubleMatParam(LoadLoveKEnum,love_k,nl,1)); - parameters->AddObject(new DoubleMatParam(LoadLoveLEnum,love_l,nl,1)); - parameters->AddObject(new DoubleMatParam(TidalLoveHEnum,love_th,nl,1)); - parameters->AddObject(new DoubleMatParam(TidalLoveKEnum,love_tk,nl,1)); - parameters->AddObject(new DoubleMatParam(TidalLoveLEnum,love_tl,nl,1)); - - /*compute elastic green function for a range of angles*/ - iomodel->FetchData(°acc,"md.solidearth.settings.degacc"); - M=reCast(180./degacc+1.); - - // AD performance is sensitive to calls to ensurecontiguous. - // // Providing "t" will cause ensurecontiguous to be called. - #ifdef _HAVE_AD_ - G_rigid=xNew(M,"t"); - G_elastic=xNew(M,"t"); - U_elastic=xNew(M,"t"); - H_elastic=xNew(M,"t"); - #else - G_rigid=xNew(M); - G_elastic=xNew(M); - U_elastic=xNew(M); - H_elastic=xNew(M); - #endif - - /*compute combined legendre + love number (elastic green function:*/ - m=DetermineLocalSize(M,IssmComm::GetComm()); - GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm()); - // AD performance is sensitive to calls to ensurecontiguous. - // // Providing "t" will cause ensurecontiguous to be called. - #ifdef _HAVE_AD_ - G_elastic_local=xNew(m,"t"); - G_rigid_local=xNew(m,"t"); - U_elastic_local=xNew(m,"t"); - H_elastic_local=xNew(m,"t"); - #else - G_elastic_local=xNew(m); - G_rigid_local=xNew(m); - U_elastic_local=xNew(m); - H_elastic_local=xNew(m); - #endif - - for(int i=lower_row;i(i)*degacc * PI / 180.0; - - G_rigid_local[i-lower_row]= .5/sin(alpha/2.0); - G_elastic_local[i-lower_row]= (love_k[nl-1]-love_h[nl-1])*G_rigid_local[i-lower_row]; - U_elastic_local[i-lower_row]= (love_h[nl-1])*G_rigid_local[i-lower_row]; - H_elastic_local[i-lower_row]= 0; - IssmDouble Pn = 0.; - IssmDouble Pn1 = 0.; - IssmDouble Pn2 = 0.; - IssmDouble Pn_p = 0.; - IssmDouble Pn_p1 = 0.; - IssmDouble Pn_p2 = 0.; - - for (int n=0;nFetchData(&rotation,"md.solidearth.settings.rotation"); + + if(elastic | rigid){ + /*compute green functions for a range of angles*/ + iomodel->FetchData(°acc,"md.solidearth.settings.degacc"); + M=reCast(180./degacc+1.); + } + + /*love numbers: */ + if(elastic){ + iomodel->FetchData(&love_h,&nl,NULL,"md.solidearth.lovenumbers.h"); + iomodel->FetchData(&love_k,&nl,NULL,"md.solidearth.lovenumbers.k"); + iomodel->FetchData(&love_l,&nl,NULL,"md.solidearth.lovenumbers.l"); + iomodel->FetchData(&love_th,&nl,NULL,"md.solidearth.lovenumbers.th"); + iomodel->FetchData(&love_tk,&nl,NULL,"md.solidearth.lovenumbers.tk"); + iomodel->FetchData(&love_tl,&nl,NULL,"md.solidearth.lovenumbers.tl"); + + parameters->AddObject(new DoubleMatParam(LoadLoveHEnum,love_h,nl,1)); + parameters->AddObject(new DoubleMatParam(LoadLoveKEnum,love_k,nl,1)); + parameters->AddObject(new DoubleMatParam(LoadLoveLEnum,love_l,nl,1)); + parameters->AddObject(new DoubleMatParam(TidalLoveHEnum,love_th,nl,1)); + parameters->AddObject(new DoubleMatParam(TidalLoveKEnum,love_tk,nl,1)); + parameters->AddObject(new DoubleMatParam(TidalLoveLEnum,love_tl,nl,1)); + + // AD performance is sensitive to calls to ensurecontiguous. + // // Providing "t" will cause ensurecontiguous to be called. + #ifdef _HAVE_AD_ + G_elastic=xNew(M,"t"); + U_elastic=xNew(M,"t"); + H_elastic=xNew(M,"t"); + #else + G_elastic=xNew(M); + U_elastic=xNew(M); + H_elastic=xNew(M); + #endif + } + if(rigid){ + #ifdef _HAVE_AD_ + G_rigid=xNew(M,"t"); + #else + G_rigid=xNew(M); + #endif + } + + if(rotation)parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.lovenumbers.tk2secular",TidalLoveK2SecularEnum)); + + if(rigid | elastic){ + + /*compute combined legendre + love number (elastic green function:*/ + m=DetermineLocalSize(M,IssmComm::GetComm()); + GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm()); + } + if(elastic){ + #ifdef _HAVE_AD_ + G_elastic_local=xNew(m,"t"); + U_elastic_local=xNew(m,"t"); + H_elastic_local=xNew(m,"t"); + #else + G_elastic_local=xNew(m); + U_elastic_local=xNew(m); + H_elastic_local=xNew(m); + #endif + } + if(rigid){ + #ifdef _HAVE_AD_ + G_rigid_local=xNew(m,"t"); + #else + G_rigid_local=xNew(m); + #endif + } + + if(rigid){ + for(int i=lower_row;i(i)*degacc * PI / 180.0; + G_rigid_local[i-lower_row]= .5/sin(alpha/2.0); + } + } + if(elastic){ + for(int i=lower_row;i(i)*degacc * PI / 180.0; + + G_elastic_local[i-lower_row]= (love_k[nl-1]-love_h[nl-1])*G_rigid_local[i-lower_row]; + U_elastic_local[i-lower_row]= (love_h[nl-1])*G_rigid_local[i-lower_row]; + H_elastic_local[i-lower_row]= 0; + IssmDouble Pn = 0.; + IssmDouble Pn1 = 0.; + IssmDouble Pn2 = 0.; + IssmDouble Pn_p = 0.; + IssmDouble Pn_p1 = 0.; + IssmDouble Pn_p2 = 0.; + + for (int n=0;n(IssmComm::GetSize()); - int* displs=xNew(IssmComm::GetSize()); - - //recvcounts: - ISSM_MPI_Allgather(&m,1,ISSM_MPI_INT,recvcounts,1,ISSM_MPI_INT,IssmComm::GetComm()); - - /*displs: */ - ISSM_MPI_Allgather(&lower_row,1,ISSM_MPI_INT,displs,1,ISSM_MPI_INT,IssmComm::GetComm()); - - /*All gather:*/ - ISSM_MPI_Allgatherv(G_rigid_local, m, ISSM_MPI_DOUBLE, G_rigid, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm()); + } + } + if(rigid){ + + /*merge G_elastic_local into G_elastic; U_elastic_local into U_elastic; H_elastic_local to H_elastic:{{{*/ + int* recvcounts=xNew(IssmComm::GetSize()); + int* displs=xNew(IssmComm::GetSize()); + + //recvcounts: + ISSM_MPI_Allgather(&m,1,ISSM_MPI_INT,recvcounts,1,ISSM_MPI_INT,IssmComm::GetComm()); + + /*displs: */ + ISSM_MPI_Allgather(&lower_row,1,ISSM_MPI_INT,displs,1,ISSM_MPI_INT,IssmComm::GetComm()); + + /*All gather:*/ + ISSM_MPI_Allgatherv(G_rigid_local, m, ISSM_MPI_DOUBLE, G_rigid, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm()); + if(elastic){ ISSM_MPI_Allgatherv(G_elastic_local, m, ISSM_MPI_DOUBLE, G_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm()); ISSM_MPI_Allgatherv(U_elastic_local, m, ISSM_MPI_DOUBLE, U_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm()); ISSM_MPI_Allgatherv(H_elastic_local, m, ISSM_MPI_DOUBLE, H_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm()); - - /*free resources: */ - xDelete(recvcounts); - xDelete(displs); - - /*}}}*/ - - /*Avoid singularity at 0: */ - G_rigid[0]=G_rigid[1]; - parameters->AddObject(new DoubleVecParam(SealevelriseGRigidEnum,G_rigid,M)); - G_elastic[0]=G_elastic[1]; + } + + /*free resources: */ + xDelete(recvcounts); + xDelete(displs); + + /*Avoid singularity at 0: */ + G_rigid[0]=G_rigid[1]; + G_elastic[0]=G_elastic[1]; + U_elastic[0]=U_elastic[1]; + H_elastic[0]=H_elastic[1]; + + /*Save our precomputed tables into parameters*/ + parameters->AddObject(new DoubleVecParam(SealevelriseGRigidEnum,G_rigid,M)); + if(elastic){ parameters->AddObject(new DoubleVecParam(SealevelriseGElasticEnum,G_elastic,M)); - U_elastic[0]=U_elastic[1]; parameters->AddObject(new DoubleVecParam(SealevelriseUElasticEnum,U_elastic,M)); - H_elastic[0]=H_elastic[1]; parameters->AddObject(new DoubleVecParam(SealevelriseHElasticEnum,H_elastic,M)); - - /*free resources: */ + } + + /*free resources: */ + xDelete(G_rigid); + xDelete(G_rigid_local); + if(elastic){ xDelete(love_h); xDelete(love_k); @@ -370,6 +429,4 @@ xDelete(love_tk); xDelete(love_tl); - xDelete(G_rigid); - xDelete(G_rigid_local); xDelete(G_elastic); xDelete(G_elastic_local); @@ -378,67 +435,8 @@ xDelete(H_elastic); xDelete(H_elastic_local); - } else { /*elastic==false*/ - /*compute elastic green function for a range of angles*/ - iomodel->FetchData(°acc,"md.solidearth.settings.degacc"); - M=reCast(180./degacc+1.); - - // AD performance is sensitive to calls to ensurecontiguous. - // // Providing "t" will cause ensurecontiguous to be called. - #ifdef _HAVE_AD_ - G_rigid=xNew(M,"t"); - #else - G_rigid=xNew(M); - #endif - - /*compute combined legendre + love number (elastic green function:*/ - m=DetermineLocalSize(M,IssmComm::GetComm()); - GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm()); - // AD performance is sensitive to calls to ensurecontiguous. - // // Providing "t" will cause ensurecontiguous to be called. - #ifdef _HAVE_AD_ - G_rigid_local=xNew(m,"t"); - #else - G_rigid_local=xNew(m); - #endif - - for(int i=lower_row;i(i)*degacc*PI/180.0; - G_rigid_local[i-lower_row]=.5/sin(alpha/2.0); - } - - /*merge G_elastic_local into G_elastic; U_elastic_local into U_elastic; H_elastic_local to H_elastic:{{{*/ - int* recvcounts=xNew(IssmComm::GetSize()); - int* displs=xNew(IssmComm::GetSize()); - - //recvcounts: - ISSM_MPI_Allgather(&m,1,ISSM_MPI_INT,recvcounts,1,ISSM_MPI_INT,IssmComm::GetComm()); - - /*displs: */ - ISSM_MPI_Allgather(&lower_row,1,ISSM_MPI_INT,displs,1,ISSM_MPI_INT,IssmComm::GetComm()); - - /*All gather:*/ - ISSM_MPI_Allgatherv(G_rigid_local, m, ISSM_MPI_DOUBLE, G_rigid, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm()); - - /*free resources: */ - xDelete(recvcounts); - xDelete(displs); - - /*}}}*/ - - /*Avoid singularity at 0: */ - G_rigid[0]=G_rigid[1]; - parameters->AddObject(new DoubleVecParam(SealevelriseGRigidEnum,G_rigid,M)); - - /*free resources: */ - xDelete(G_rigid); - xDelete(G_rigid_local); - } - } else { /*rigid==false*/ - if (elastic) { - _error_("Must set md.solidearth.settings.rigid=1 when setting md.solidearth.settings.elastic=1"); - } - } - + } + } + /*}}}*/ + /*Indicate we have not yet run the Geometry Core module: */ parameters->AddObject(new BoolParam(SealevelriseGeometryDoneEnum,false)); Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 25750) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 25751) @@ -5463,158 +5463,119 @@ this->inputs->SetDoubleInput(AreaEnum,this->lid,area); this->AddInput(SealevelAreaEnum,&area,P0Enum); + if(!computerigid)return; /*recover precomputed green function kernels:*/ - if(computerigid){ - DoubleVecParam* parameter = static_cast(this->parameters->FindParamObject(SealevelriseGRigidEnum)); _assert_(parameter); - parameter->GetParameterValueByPointer((IssmDouble**)&G_rigid_precomputed,&M); + DoubleVecParam* parameter = static_cast(this->parameters->FindParamObject(SealevelriseGRigidEnum)); _assert_(parameter); + parameter->GetParameterValueByPointer((IssmDouble**)&G_rigid_precomputed,&M); + + /*allocate indices:*/ + indices=xNew(gsize); + G=xNewZeroInit(gsize); + + if(computeelastic){ + parameter = static_cast(this->parameters->FindParamObject(SealevelriseGElasticEnum)); _assert_(parameter); + parameter->GetParameterValueByPointer((IssmDouble**)&G_elastic_precomputed,&M); + + parameter = static_cast(this->parameters->FindParamObject(SealevelriseUElasticEnum)); _assert_(parameter); + parameter->GetParameterValueByPointer((IssmDouble**)&U_elastic_precomputed,&M); + + parameter = static_cast(this->parameters->FindParamObject(SealevelriseHElasticEnum)); _assert_(parameter); + parameter->GetParameterValueByPointer((IssmDouble**)&H_elastic_precomputed,&M); /*allocate indices:*/ - indices=xNew(gsize); - G=xNewZeroInit(gsize); - + GU=xNewZeroInit(gsize); + if(horiz){ + GN=xNewZeroInit(gsize); + GE=xNewZeroInit(gsize); + } + } + + /* Where is the centroid of this element:{{{*/ + + /*retrieve coordinates: */ + ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + + /*figure out gravity center of our element (Cartesian): */ + x_element=(xyz_list[0][0]+xyz_list[1][0]+xyz_list[2][0])/3.0; + y_element=(xyz_list[0][1]+xyz_list[1][1]+xyz_list[2][1])/3.0; + z_element=(xyz_list[0][2]+xyz_list[1][2]+xyz_list[2][2])/3.0; + + /*compute gravity center in lat,long: */ + late= asin(z_element/planetradius); + longe = atan2(y_element,x_element); + /*}}}*/ + + constant=3/rho_earth*area/planetarea; + for(int i=0;iPI)delLambda=2*PI-delLambda; + alpha=2.*asin(sqrt(pow(sin(delPhi/2),2)+cos(lati)*cos(late)*pow(sin(delLambda/2),2))); + indices[i]=alpha/PI*reCast(M-1); + index=reCast(indices[i]); + + /*Rigid earth gravitational perturbation: */ + G[i]+=G_rigid_precomputed[index]; + if(elastic) G[i]+=G_elastic_precomputed[index]; + G[i]=G[i]*constant; + + /*Elastic components:*/ if(computeelastic){ - parameter = static_cast(this->parameters->FindParamObject(SealevelriseGElasticEnum)); _assert_(parameter); - parameter->GetParameterValueByPointer((IssmDouble**)&G_elastic_precomputed,&M); - - parameter = static_cast(this->parameters->FindParamObject(SealevelriseUElasticEnum)); _assert_(parameter); - parameter->GetParameterValueByPointer((IssmDouble**)&U_elastic_precomputed,&M); - - parameter = static_cast(this->parameters->FindParamObject(SealevelriseHElasticEnum)); _assert_(parameter); - parameter->GetParameterValueByPointer((IssmDouble**)&H_elastic_precomputed,&M); - - /*allocate indices:*/ - GU=xNewZeroInit(gsize); + GU[i] = constant * U_elastic_precomputed[index]; if(horiz){ - GN=xNewZeroInit(gsize); - GE=xNewZeroInit(gsize); - } - - /* Where is the centroid of this element:*/ - - /*retrieve coordinates: */ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - - /*figure out gravity center of our element (Cartesian): */ - x_element=(xyz_list[0][0]+xyz_list[1][0]+xyz_list[2][0])/3.0; - y_element=(xyz_list[0][1]+xyz_list[1][1]+xyz_list[2][1])/3.0; - z_element=(xyz_list[0][2]+xyz_list[1][2]+xyz_list[2][2])/3.0; - - /*compute gravity center in lat,long: */ - late= asin(z_element/planetradius); - longe = atan2(y_element,x_element); - - constant=3/rho_earth*area/planetarea; - - for(int i=0;iPI)delLambda=2*PI-delLambda; - alpha=2.*asin(sqrt(pow(sin(delPhi/2),2)+cos(lati)*cos(late)*pow(sin(delLambda/2),2))); - indices[i]=alpha/PI*reCast(M-1); - index=reCast(indices[i]); - - /*Rigid earth gravitational perturbation: */ - G[i]+=G_rigid_precomputed[index]; - G[i]+=G_elastic_precomputed[index]; - G[i]=G[i]*constant; - - /*Elastic components:*/ - GU[i] = constant * U_elastic_precomputed[index]; - if(horiz){ - /*Compute azimuths, both north and east components: */ - x = xx[i]; y = yy[i]; z = zz[i]; - if(latitude[i]==90){ - x=1e-12; y=1e-12; - } - if(latitude[i]==-90){ - x=1e-12; y=1e-12; - } - dx = x_element-x; dy = y_element-y; dz = z_element-z; - N_azim = (-z*x*dx-z*y*dy+(pow(x,2)+pow(y,2))*dz) /pow((pow(x,2)+pow(y,2))*(pow(x,2)+pow(y,2)+pow(z,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5); - E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5); - - GN[i] = constant*H_elastic_precomputed[index]*N_azim; - GE[i] = constant*H_elastic_precomputed[index]*E_azim; + /*Compute azimuths, both north and east components: */ + x = xx[i]; y = yy[i]; z = zz[i]; + if(latitude[i]==90){ + x=1e-12; y=1e-12; } - } - - /*Add in inputs:*/ - this->inputs->SetArrayInput(SealevelriseIndicesEnum,this->lid,indices,gsize); - this->inputs->SetArrayInput(SealevelriseGEnum,this->lid,G,gsize); - this->inputs->SetArrayInput(SealevelriseGUEnum,this->lid,GU,gsize); - if(horiz){ - this->inputs->SetArrayInput(SealevelriseGNEnum,this->lid,GN,gsize); - this->inputs->SetArrayInput(SealevelriseGEEnum,this->lid,GE,gsize); - } - this->inputs->SetDoubleInput(AreaEnum,this->lid,area); - this->AddInput(SealevelAreaEnum,&area,P0Enum); - - /*Free allocations:*/ - #ifdef _HAVE_RESTRICT_ - delete GU; - if(horiz){ - delete GN; - delete GE; - } - #else - xDelete(GU); - if(horiz){ - xDelete(GN); - xDelete(GE); - } - #endif - } else { /*computeelastic==false*/ - /* Where is the centroid of this element:*/ - - /*retrieve coordinates: */ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - - /*figure out gravity center of our element (Cartesian): */ - x_element=(xyz_list[0][0]+xyz_list[1][0]+xyz_list[2][0])/3.0; - y_element=(xyz_list[0][1]+xyz_list[1][1]+xyz_list[2][1])/3.0; - z_element=(xyz_list[0][2]+xyz_list[1][2]+xyz_list[2][2])/3.0; - - /*compute gravity center in lat,long: */ - late= asin(z_element/planetradius); - longe = atan2(y_element,x_element); - - constant=3/rho_earth*area/planetarea; - for(int i=0;iPI)delLambda=2*PI-delLambda; - alpha=2.*asin(sqrt(pow(sin(delPhi/2),2)+cos(lati)*cos(late)*pow(sin(delLambda/2),2))); - indices[i]=alpha/PI*reCast(M-1); - index=reCast(indices[i]); - - /*Rigid earth gravitational perturbation: */ - G[i]+=G_rigid_precomputed[index]; - G[i]=G[i]*constant; - } - - /*Add in inputs:*/ - this->inputs->SetArrayInput(SealevelriseIndicesEnum,this->lid,indices,gsize); - this->inputs->SetArrayInput(SealevelriseGEnum,this->lid,G,gsize); - this->inputs->SetDoubleInput(AreaEnum,this->lid,area); - this->AddInput(SealevelAreaEnum,&area,P0Enum); - } - - /*Free allocations:*/ - #ifdef _HAVE_RESTRICT_ - delete indices; - delete G; - #else - xDelete(indices); - xDelete(G); - #endif - } else { /*computerigid==false*/ - // do nothing - } + if(latitude[i]==-90){ + x=1e-12; y=1e-12; + } + dx = x_element-x; dy = y_element-y; dz = z_element-z; + N_azim = (-z*x*dx-z*y*dy+(pow(x,2)+pow(y,2))*dz) /pow((pow(x,2)+pow(y,2))*(pow(x,2)+pow(y,2)+pow(z,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5); + E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5); + + GN[i] = constant*H_elastic_precomputed[index]*N_azim; + GE[i] = constant*H_elastic_precomputed[index]*E_azim; + } + } + } + + /*Add in inputs:*/ + this->inputs->SetArrayInput(SealevelriseIndicesEnum,this->lid,indices,gsize); + this->inputs->SetArrayInput(SealevelriseGEnum,this->lid,G,gsize); + if(computeelastic){ + this->inputs->SetArrayInput(SealevelriseGUEnum,this->lid,GU,gsize); + if(horiz){ + this->inputs->SetArrayInput(SealevelriseGNEnum,this->lid,GN,gsize); + this->inputs->SetArrayInput(SealevelriseGEEnum,this->lid,GE,gsize); + } + } + this->inputs->SetDoubleInput(AreaEnum,this->lid,area); + this->AddInput(SealevelAreaEnum,&area,P0Enum); + + /*Free allocations:*/ + #ifdef _HAVE_RESTRICT_ + delete G; + if(computeelastic){ + delete GU; + if(horiz){ + delete GN; + delete GE; + } + } + #else + xDelete(G); + if(elastic){ + xDelete(GU); + if(horiz){ + xDelete(GN); + xDelete(GE); + } + } + #endif return; @@ -5631,4 +5592,5 @@ bool scaleoceanarea= false; bool computerigid= false; + int glfraction=1; /*output: */ @@ -5690,4 +5652,8 @@ phi=this->GetGroundedPortion(&xyz_list[0][0]); //watch out, this only works because of the Thales theorem! We are in 3D, but this routine is inherently for 2D trias + if(glfraction==0)phi=1; + #ifdef _ISSM_DEBUG_ + this->AddInput2(SealevelEustaticMaskEnum,&phi,P0Enum); + #endif } else phi=1.0; @@ -5737,4 +5703,9 @@ } + /*Plug bslrice into barystatic contribution vector:*/ + if(barystatic_contribution){ + id=partition[this->Sid()]+1; + barystatic_contribution->SetValue(id,bslrice,ADD_VAL); + } /*return :*/ return bslrice; @@ -5762,9 +5733,10 @@ IssmDouble bslrhydro = 0; - /*early return if we are on an ice cap:*/ - if(masks->isiceonly[this->lid]){ + /*early return if we are on an ice cap. Nop, we may have hydro + * and ice on the same masscon:*/ + /*if(masks->isiceonly[this->lid]){ bslrhydro=0; return bslrhydro; - } + }*/ /*early return if we are fully floating:*/ @@ -5807,4 +5779,9 @@ } + /*Plug bslrice into barystatic contribution vector:*/ + if(barystatic_contribution){ + id=partition[this->Sid()]+1; + barystatic_contribution->SetValue(id,bslrhydro,ADD_VAL); + } /*output:*/ return bslrhydro; Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 25750) +++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 25751) @@ -114,4 +114,6 @@ syn keyword cConstant CumBslrIceEnum syn keyword cConstant CumBslrHydroEnum +syn keyword cConstant CumBslrIcePartitionEnum +syn keyword cConstant CumBslrHydroPartitionEnum syn keyword cConstant CumGmtslrEnum syn keyword cConstant CumGmslrEnum @@ -328,4 +330,8 @@ syn keyword cConstant ModelnameEnum syn keyword cConstant SaveResultsEnum +syn keyword cConstant SolidearthPartitionIceEnum +syn keyword cConstant SolidearthPartitionHydroEnum +syn keyword cConstant SolidearthNpartIceEnum +syn keyword cConstant SolidearthNpartHydroEnum syn keyword cConstant SolidearthPlanetRadiusEnum syn keyword cConstant SolidearthPlanetAreaEnum @@ -345,4 +351,5 @@ syn keyword cConstant SealevelriseGElasticEnum syn keyword cConstant SolidearthSettingsComputesealevelchangeEnum +syn keyword cConstant SolidearthSettingsGlfractionEnum syn keyword cConstant SolidearthSettingsRunFrequencyEnum syn keyword cConstant SealevelriseHElasticEnum Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 25750) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 25751) @@ -108,4 +108,6 @@ CumBslrIceEnum, CumBslrHydroEnum, + CumBslrIcePartitionEnum, + CumBslrHydroPartitionEnum, CumGmtslrEnum, CumGmslrEnum, @@ -322,4 +324,8 @@ ModelnameEnum, SaveResultsEnum, + SolidearthPartitionIceEnum, + SolidearthPartitionHydroEnum, + SolidearthNpartIceEnum, + SolidearthNpartHydroEnum, SolidearthPlanetRadiusEnum, SolidearthPlanetAreaEnum, @@ -339,4 +345,5 @@ SealevelriseGElasticEnum, SolidearthSettingsComputesealevelchangeEnum, + SolidearthSettingsGlfractionEnum, SolidearthSettingsRunFrequencyEnum, SealevelriseHElasticEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 25750) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 25751) @@ -116,4 +116,6 @@ case CumBslrIceEnum : return "CumBslrIce"; case CumBslrHydroEnum : return "CumBslrHydro"; + case CumBslrIcePartitionEnum : return "CumBslrIcePartition"; + case CumBslrHydroPartitionEnum : return "CumBslrHydroPartition"; case CumGmtslrEnum : return "CumGmtslr"; case CumGmslrEnum : return "CumGmslr"; @@ -330,4 +332,8 @@ case ModelnameEnum : return "Modelname"; case SaveResultsEnum : return "SaveResults"; + case SolidearthPartitionIceEnum : return "SolidearthPartitionIce"; + case SolidearthPartitionHydroEnum : return "SolidearthPartitionHydro"; + case SolidearthNpartIceEnum : return "SolidearthNpartIce"; + case SolidearthNpartHydroEnum : return "SolidearthNpartHydro"; case SolidearthPlanetRadiusEnum : return "SolidearthPlanetRadius"; case SolidearthPlanetAreaEnum : return "SolidearthPlanetArea"; @@ -347,4 +353,5 @@ case SealevelriseGElasticEnum : return "SealevelriseGElastic"; case SolidearthSettingsComputesealevelchangeEnum : return "SolidearthSettingsComputesealevelchange"; + case SolidearthSettingsGlfractionEnum : return "SolidearthSettingsGlfraction"; case SolidearthSettingsRunFrequencyEnum : return "SolidearthSettingsRunFrequency"; case SealevelriseHElasticEnum : return "SealevelriseHElastic"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 25750) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 25751) @@ -116,4 +116,6 @@ else if (strcmp(name,"CumBslrIce")==0) return CumBslrIceEnum; else if (strcmp(name,"CumBslrHydro")==0) return CumBslrHydroEnum; + else if (strcmp(name,"CumBslrIcePartition")==0) return CumBslrIcePartitionEnum; + else if (strcmp(name,"CumBslrHydroPartition")==0) return CumBslrHydroPartitionEnum; else if (strcmp(name,"CumGmtslr")==0) return CumGmtslrEnum; else if (strcmp(name,"CumGmslr")==0) return CumGmslrEnum; @@ -135,10 +137,10 @@ else if (strcmp(name,"DebugProfiling")==0) return DebugProfilingEnum; else if (strcmp(name,"DomainDimension")==0) return DomainDimensionEnum; - else if (strcmp(name,"DomainType")==0) return DomainTypeEnum; - else if (strcmp(name,"DslModel")==0) return DslModelEnum; else stage=2; } if(stage==2){ - if (strcmp(name,"DslModelid")==0) return DslModelidEnum; + if (strcmp(name,"DomainType")==0) return DomainTypeEnum; + else if (strcmp(name,"DslModel")==0) return DslModelEnum; + else if (strcmp(name,"DslModelid")==0) return DslModelidEnum; else if (strcmp(name,"DslNummodels")==0) return DslNummodelsEnum; else if (strcmp(name,"DslComputeFingerprints")==0) return DslComputeFingerprintsEnum; @@ -258,10 +260,10 @@ else if (strcmp(name,"LoveMu0")==0) return LoveMu0Enum; else if (strcmp(name,"LoveNfreq")==0) return LoveNfreqEnum; - else if (strcmp(name,"LoveR0")==0) return LoveR0Enum; - else if (strcmp(name,"LoveShNmax")==0) return LoveShNmaxEnum; else stage=3; } if(stage==3){ - if (strcmp(name,"LoveShNmin")==0) return LoveShNminEnum; + if (strcmp(name,"LoveR0")==0) return LoveR0Enum; + else if (strcmp(name,"LoveShNmax")==0) return LoveShNmaxEnum; + else if (strcmp(name,"LoveShNmin")==0) return LoveShNminEnum; else if (strcmp(name,"MassFluxSegments")==0) return MassFluxSegmentsEnum; else if (strcmp(name,"MassFluxSegmentsPresent")==0) return MassFluxSegmentsPresentEnum; @@ -336,4 +338,8 @@ else if (strcmp(name,"Modelname")==0) return ModelnameEnum; else if (strcmp(name,"SaveResults")==0) return SaveResultsEnum; + else if (strcmp(name,"SolidearthPartitionIce")==0) return SolidearthPartitionIceEnum; + else if (strcmp(name,"SolidearthPartitionHydro")==0) return SolidearthPartitionHydroEnum; + else if (strcmp(name,"SolidearthNpartIce")==0) return SolidearthNpartIceEnum; + else if (strcmp(name,"SolidearthNpartHydro")==0) return SolidearthNpartHydroEnum; else if (strcmp(name,"SolidearthPlanetRadius")==0) return SolidearthPlanetRadiusEnum; else if (strcmp(name,"SolidearthPlanetArea")==0) return SolidearthPlanetAreaEnum; @@ -353,4 +359,5 @@ else if (strcmp(name,"SealevelriseGElastic")==0) return SealevelriseGElasticEnum; else if (strcmp(name,"SolidearthSettingsComputesealevelchange")==0) return SolidearthSettingsComputesealevelchangeEnum; + else if (strcmp(name,"SolidearthSettingsGlfraction")==0) return SolidearthSettingsGlfractionEnum; else if (strcmp(name,"SolidearthSettingsRunFrequency")==0) return SolidearthSettingsRunFrequencyEnum; else if (strcmp(name,"SealevelriseHElastic")==0) return SealevelriseHElasticEnum; @@ -376,5 +383,8 @@ else if (strcmp(name,"SmbAIce")==0) return SmbAIceEnum; else if (strcmp(name,"SmbAIdx")==0) return SmbAIdxEnum; - else if (strcmp(name,"SmbASnow")==0) return SmbASnowEnum; + else stage=4; + } + if(stage==4){ + if (strcmp(name,"SmbASnow")==0) return SmbASnowEnum; else if (strcmp(name,"SmbAccualti")==0) return SmbAccualtiEnum; else if (strcmp(name,"SmbAccugrad")==0) return SmbAccugradEnum; @@ -383,8 +393,5 @@ else if (strcmp(name,"SmbAveraging")==0) return SmbAveragingEnum; else if (strcmp(name,"SmbDesfac")==0) return SmbDesfacEnum; - else stage=4; - } - if(stage==4){ - if (strcmp(name,"SmbDpermil")==0) return SmbDpermilEnum; + else if (strcmp(name,"SmbDpermil")==0) return SmbDpermilEnum; else if (strcmp(name,"SmbDsnowIdx")==0) return SmbDsnowIdxEnum; else if (strcmp(name,"SmbCldFrac")==0) return SmbCldFracEnum; @@ -499,5 +506,8 @@ else if (strcmp(name,"InputsSTART")==0) return InputsSTARTEnum; else if (strcmp(name,"Adjoint")==0) return AdjointEnum; - else if (strcmp(name,"Adjointp")==0) return AdjointpEnum; + else stage=5; + } + if(stage==5){ + if (strcmp(name,"Adjointp")==0) return AdjointpEnum; else if (strcmp(name,"Adjointx")==0) return AdjointxEnum; else if (strcmp(name,"Adjointy")==0) return AdjointyEnum; @@ -506,8 +516,5 @@ else if (strcmp(name,"Approximation")==0) return ApproximationEnum; else if (strcmp(name,"BalancethicknessMisfit")==0) return BalancethicknessMisfitEnum; - else stage=5; - } - if(stage==5){ - if (strcmp(name,"BalancethicknessOmega0")==0) return BalancethicknessOmega0Enum; + else if (strcmp(name,"BalancethicknessOmega0")==0) return BalancethicknessOmega0Enum; else if (strcmp(name,"BalancethicknessOmega")==0) return BalancethicknessOmegaEnum; else if (strcmp(name,"BalancethicknessSpcthickness")==0) return BalancethicknessSpcthicknessEnum; @@ -622,5 +629,8 @@ else if (strcmp(name,"FrictionTillFrictionAngle")==0) return FrictionTillFrictionAngleEnum; else if (strcmp(name,"FrictionWaterLayer")==0) return FrictionWaterLayerEnum; - else if (strcmp(name,"Frictionf")==0) return FrictionfEnum; + else stage=6; + } + if(stage==6){ + if (strcmp(name,"Frictionf")==0) return FrictionfEnum; else if (strcmp(name,"FrontalForcingsBasinId")==0) return FrontalForcingsBasinIdEnum; else if (strcmp(name,"FrontalForcingsSubglacialDischarge")==0) return FrontalForcingsSubglacialDischargeEnum; @@ -629,8 +639,5 @@ else if (strcmp(name,"GiaLithosphereThickness")==0) return GiaLithosphereThicknessEnum; else if (strcmp(name,"GiaMantleViscosity")==0) return GiaMantleViscosityEnum; - else stage=6; - } - if(stage==6){ - if (strcmp(name,"NGia")==0) return NGiaEnum; + else if (strcmp(name,"NGia")==0) return NGiaEnum; else if (strcmp(name,"NGiaRate")==0) return NGiaRateEnum; else if (strcmp(name,"UGia")==0) return UGiaEnum; @@ -745,5 +752,8 @@ else if (strcmp(name,"SealevelriseGN")==0) return SealevelriseGNEnum; else if (strcmp(name,"SedimentHead")==0) return SedimentHeadEnum; - else if (strcmp(name,"SedimentHeadOld")==0) return SedimentHeadOldEnum; + else stage=7; + } + if(stage==7){ + if (strcmp(name,"SedimentHeadOld")==0) return SedimentHeadOldEnum; else if (strcmp(name,"SedimentHeadSubstep")==0) return SedimentHeadSubstepEnum; else if (strcmp(name,"SedimentHeadTransient")==0) return SedimentHeadTransientEnum; @@ -752,8 +762,5 @@ else if (strcmp(name,"SigmaNN")==0) return SigmaNNEnum; else if (strcmp(name,"SigmaVM")==0) return SigmaVMEnum; - else stage=7; - } - if(stage==7){ - if (strcmp(name,"SmbA")==0) return SmbAEnum; + else if (strcmp(name,"SmbA")==0) return SmbAEnum; else if (strcmp(name,"SmbAValue")==0) return SmbAValueEnum; else if (strcmp(name,"SmbAccumulation")==0) return SmbAccumulationEnum; @@ -868,5 +875,8 @@ else if (strcmp(name,"SurfaceAbsVelMisfit")==0) return SurfaceAbsVelMisfitEnum; else if (strcmp(name,"Area")==0) return AreaEnum; - else if (strcmp(name,"SealevelArea")==0) return SealevelAreaEnum; + else stage=8; + } + if(stage==8){ + if (strcmp(name,"SealevelArea")==0) return SealevelAreaEnum; else if (strcmp(name,"SurfaceArea")==0) return SurfaceAreaEnum; else if (strcmp(name,"SurfaceAverageVelMisfit")==0) return SurfaceAverageVelMisfitEnum; @@ -875,8 +885,5 @@ else if (strcmp(name,"SurfaceOld")==0) return SurfaceOldEnum; else if (strcmp(name,"SurfaceLogVelMisfit")==0) return SurfaceLogVelMisfitEnum; - else stage=8; - } - if(stage==8){ - if (strcmp(name,"SurfaceLogVxVyMisfit")==0) return SurfaceLogVxVyMisfitEnum; + else if (strcmp(name,"SurfaceLogVxVyMisfit")==0) return SurfaceLogVxVyMisfitEnum; else if (strcmp(name,"SurfaceObservation")==0) return SurfaceObservationEnum; else if (strcmp(name,"SurfaceRelVelMisfit")==0) return SurfaceRelVelMisfitEnum; @@ -991,5 +998,8 @@ else if (strcmp(name,"Outputdefinition76")==0) return Outputdefinition76Enum; else if (strcmp(name,"Outputdefinition77")==0) return Outputdefinition77Enum; - else if (strcmp(name,"Outputdefinition78")==0) return Outputdefinition78Enum; + else stage=9; + } + if(stage==9){ + if (strcmp(name,"Outputdefinition78")==0) return Outputdefinition78Enum; else if (strcmp(name,"Outputdefinition79")==0) return Outputdefinition79Enum; else if (strcmp(name,"Outputdefinition7")==0) return Outputdefinition7Enum; @@ -998,8 +1008,5 @@ else if (strcmp(name,"Outputdefinition82")==0) return Outputdefinition82Enum; else if (strcmp(name,"Outputdefinition83")==0) return Outputdefinition83Enum; - else stage=9; - } - if(stage==9){ - if (strcmp(name,"Outputdefinition84")==0) return Outputdefinition84Enum; + else if (strcmp(name,"Outputdefinition84")==0) return Outputdefinition84Enum; else if (strcmp(name,"Outputdefinition85")==0) return Outputdefinition85Enum; else if (strcmp(name,"Outputdefinition86")==0) return Outputdefinition86Enum; @@ -1114,5 +1121,8 @@ else if (strcmp(name,"FemModel")==0) return FemModelEnum; else if (strcmp(name,"FileParam")==0) return FileParamEnum; - else if (strcmp(name,"FixedTimestepping")==0) return FixedTimesteppingEnum; + else stage=10; + } + if(stage==10){ + if (strcmp(name,"FixedTimestepping")==0) return FixedTimesteppingEnum; else if (strcmp(name,"FloatingArea")==0) return FloatingAreaEnum; else if (strcmp(name,"FloatingAreaScaled")==0) return FloatingAreaScaledEnum; @@ -1121,8 +1131,5 @@ else if (strcmp(name,"FreeSurfaceBaseAnalysis")==0) return FreeSurfaceBaseAnalysisEnum; else if (strcmp(name,"FreeSurfaceTopAnalysis")==0) return FreeSurfaceTopAnalysisEnum; - else stage=10; - } - if(stage==10){ - if (strcmp(name,"FrontalForcingsDefault")==0) return FrontalForcingsDefaultEnum; + else if (strcmp(name,"FrontalForcingsDefault")==0) return FrontalForcingsDefaultEnum; else if (strcmp(name,"FrontalForcingsRignot")==0) return FrontalForcingsRignotEnum; else if (strcmp(name,"Fset")==0) return FsetEnum; @@ -1237,5 +1244,8 @@ else if (strcmp(name,"MeshY")==0) return MeshYEnum; else if (strcmp(name,"MinVel")==0) return MinVelEnum; - else if (strcmp(name,"MinVx")==0) return MinVxEnum; + else stage=11; + } + if(stage==11){ + if (strcmp(name,"MinVx")==0) return MinVxEnum; else if (strcmp(name,"MinVy")==0) return MinVyEnum; else if (strcmp(name,"MinVz")==0) return MinVzEnum; @@ -1244,8 +1254,5 @@ else if (strcmp(name,"MpiDense")==0) return MpiDenseEnum; else if (strcmp(name,"Mpi")==0) return MpiEnum; - else stage=11; - } - if(stage==11){ - if (strcmp(name,"MpiSparse")==0) return MpiSparseEnum; + else if (strcmp(name,"MpiSparse")==0) return MpiSparseEnum; else if (strcmp(name,"Mumps")==0) return MumpsEnum; else if (strcmp(name,"NoFrictionOnPartiallyFloating")==0) return NoFrictionOnPartiallyFloatingEnum; @@ -1360,5 +1367,8 @@ else if (strcmp(name,"ThicknessErrorEstimator")==0) return ThicknessErrorEstimatorEnum; else if (strcmp(name,"TotalCalvingFluxLevelset")==0) return TotalCalvingFluxLevelsetEnum; - else if (strcmp(name,"TotalCalvingMeltingFluxLevelset")==0) return TotalCalvingMeltingFluxLevelsetEnum; + else stage=12; + } + if(stage==12){ + if (strcmp(name,"TotalCalvingMeltingFluxLevelset")==0) return TotalCalvingMeltingFluxLevelsetEnum; else if (strcmp(name,"TotalFloatingBmb")==0) return TotalFloatingBmbEnum; else if (strcmp(name,"TotalFloatingBmbScaled")==0) return TotalFloatingBmbScaledEnum; @@ -1367,8 +1377,5 @@ else if (strcmp(name,"TotalSmb")==0) return TotalSmbEnum; else if (strcmp(name,"TotalSmbScaled")==0) return TotalSmbScaledEnum; - else stage=12; - } - if(stage==12){ - if (strcmp(name,"TransientArrayParam")==0) return TransientArrayParamEnum; + else if (strcmp(name,"TransientArrayParam")==0) return TransientArrayParamEnum; else if (strcmp(name,"TransientInput")==0) return TransientInputEnum; else if (strcmp(name,"TransientParam")==0) return TransientParamEnum; Index: /issm/trunk-jpl/src/m/classes/model.m =================================================================== --- /issm/trunk-jpl/src/m/classes/model.m (revision 25750) +++ /issm/trunk-jpl/src/m/classes/model.m (revision 25751) @@ -192,13 +192,13 @@ methods function md = model(varargin) % {{{ + switch nargin case 0 - md=setdefaultparameters(md); - case 1 - error('model constructor not supported yet'); - + md=setdefaultparameters(md,'earth'); otherwise - error('model constructor error message: 0 of 1 argument only in input.'); + options=pairoptions(varargin{:}); + planet=getfieldvalue(options,'planet','earth'); + md=setdefaultparameters(md,planet); end end @@ -1287,5 +1287,5 @@ end end% }}} - function md = setdefaultparameters(md) % {{{ + function md = setdefaultparameters(md,planet) % {{{ %initialize subclasses @@ -1299,5 +1299,5 @@ md.friction = friction(); md.rifts = rifts(); - md.solidearth = solidearth(); + md.solidearth = solidearth(planet); md.dsl = dsl(); md.timestepping = timestepping(); Index: /issm/trunk-jpl/src/m/classes/solidearth.m =================================================================== --- /issm/trunk-jpl/src/m/classes/solidearth.m (revision 25750) +++ /issm/trunk-jpl/src/m/classes/solidearth.m (revision 25751) @@ -14,4 +14,6 @@ requested_outputs = {}; transitions = {}; + partitionice = []; + partitionhydro = []; end methods @@ -19,10 +21,12 @@ switch nargin case 0 - self=setdefaultparameters(self); + self=setdefaultparameters(self,earth); + case 1 + self=setdefaultparameters(self,varargin{:}); otherwise - error('constructor not supported'); + error('solidearth constructor error message: zero or one argument only!'); end end % }}} - function self = setdefaultparameters(self) % {{{ + function self = setdefaultparameters(self,planet) % {{{ %output default: @@ -32,6 +36,10 @@ self.transitions={}; + %no partitions requested for barystatic contribution: + self.partitionice=[]; + self.partitionhydro=[]; + %earth radius - self.planetradius= planetradius('earth'); + self.planetradius= planetradius(planet); end % }}} @@ -62,4 +70,6 @@ fielddisplay(self,'transitions','indices into parts of the mesh that will be icecaps'); fielddisplay(self,'requested_outputs','additional outputs requested'); + fielddisplay(self,'partitionice','ice partition vector for barystatic contribution'); + fielddisplay(self,'partitionhydro','hydro partition vector for barystatic contribution'); self.settings.disp(); self.surfaceload.disp(); @@ -73,5 +83,22 @@ WriteData(fid,prefix,'object',self,'fieldname','planetradius','format','Double'); WriteData(fid,prefix,'object',self,'fieldname','transitions','format','MatArray'); + + if ~isempty(self.partitionice), + npartice=max(self.partitionice)+2; + else + npartice=0; + end + if ~isempty(self.partitionhydro), + nparthydro=max(self.partitionhydro)+2; + else + nparthydro=0; + end + + WriteData(fid,prefix,'object',self,'fieldname','partitionice','mattype',1,'format','DoubleMat'); + WriteData(fid,prefix,'data',npartice,'format','Integer','name','md.solidearth.npartice'); + WriteData(fid,prefix,'object',self,'fieldname','partitionhydro','mattype',1,'format','DoubleMat'); + WriteData(fid,prefix,'data',nparthydro,'format','Integer','name','md.solidearth.nparthydro'); + self.settings.marshall(prefix,md,fid); self.surfaceload.marshall(prefix,md,fid); @@ -98,4 +125,5 @@ writejscellstring(fid,[modelname '.solidearth.requested_outputs'],self.requested_outputs); writejscellarray(fid,[modelname '.solidearth.transitions'],self.transitions); + writejscellarray(fid,[modelname '.solidearth.partition'],self.partition); end % }}} function self = extrude(self,md) % {{{ Index: /issm/trunk-jpl/src/m/classes/solidearthsettings.m =================================================================== --- /issm/trunk-jpl/src/m/classes/solidearthsettings.m (revision 25750) +++ /issm/trunk-jpl/src/m/classes/solidearthsettings.m (revision 25751) @@ -64,4 +64,9 @@ md = checkfield(md,'fieldname','solidearth.settings.horiz','NaN',1,'Inf',1,'values',[0 1]); + %checks on computational flags: + if self.elastic==1 & self.rigid==0, + error('solidearthsettings checkconsistency error message: need rigid on if elastic flag is set'); + end + %a coupler to a planet model is provided. if self.computesealevelchange, Index: /issm/trunk-jpl/src/m/geometry/planetradius.m =================================================================== --- /issm/trunk-jpl/src/m/geometry/planetradius.m (revision 25750) +++ /issm/trunk-jpl/src/m/geometry/planetradius.m (revision 25751) @@ -10,4 +10,6 @@ if strcmpi(planet,'earth'), radius=6.371012*10^6; +elseif strcmpi(planet,'europa'), + radius=1.5008*10^6; else error(['planet type ' planet ' not supported yet!']);