Index: /issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp (revision 24504) +++ /issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp (revision 24505) @@ -58,5 +58,5 @@ /*dynamic sea level: */ iomodel->FetchData(&dslmodel,"md.dsl.model"); - if (dslmodel==1){ + if (dslmodel==1){ /*standard dsl model:{{{*/ /*deal with global mean steric rate: */ @@ -94,10 +94,183 @@ /*deal with dynamic sea level fields: */ - iomodel->FetchDataToInput(inputs2,elements,"md.dsl.sea_surface_height_change_above_geoid", DslSeaSurfaceHeightChangeAboveGeoidEnum); iomodel->FetchDataToInput(inputs2,elements,"md.dsl.sea_water_pressure_change_at_sea_floor", DslSeaWaterPressureChangeAtSeaFloor); - //TransientInput2* transientinput=inputs2->SetDatasetTransientInput(DslGlobalAverageThermostericSeaLevelChangeEnum,counter, times,N); - } + } /*}}}*/ + else if (dslmodel==2){ /*multi-model ensemble dsl model:{{{*/ + + /*variables:*/ + int nummodels; + IssmDouble** pstr=NULL; + IssmDouble* str=NULL; + IssmDouble* times = NULL; + int* pM = NULL; + int* pN = NULL; + int M,N; + + /*deal with dsl.sea_surface_height_change_above_geoid {{{*/ + iomodel->FetchData(&pstr,&pM,&pN,&nummodels,"md.dsl.global_average_thermosteric_sea_level_change"); + + /*go through the mat array and create a dataset of transient inputs:*/ + for (int i=0;i(N); + for(int t=0;tSetDatasetTransientInput(DslGlobalAverageThermostericSeaLevelChangeEnum,i, times,N); + + for(int j=0;jSize();j++){ + Element* element=xDynamicCast(elements->GetObjectByOffset(j)); + + for(int t=0;tObjectEnum()){ + case TriaEnum: transientinput->AddTriaTimeInput( t,1,&element->lid,&str[t],P0Enum); break; + case PentaEnum: transientinput->AddPentaTimeInput(t,1,&element->lid,&str[t],P0Enum); break; + default: _error_("Not implemented yet"); + } + } + } + xDelete(times); + } + /*Delete data:*/ + for(int i=0;i(str); + } + xDelete(pstr); + xDelete(pM); + xDelete(pN); + /*}}}*/ + /*now do the same with the dsl.sea_surface_height_change_above_geoid field:{{{ */ + iomodel->FetchData(&pstr,&pM,&pN,&nummodels,"md.dsl.sea_surface_height_change_above_geoid"); + + for (int i=0;i(N); + for(int t=0;tSetDatasetTransientInput(DslSeaSurfaceHeightChangeAboveGeoidEnum,i, times,N); + + for(int j=0;jSize();j++){ + + /*Get the right transient input*/ + Element* element=xDynamicCast(elements->GetObjectByOffset(j)); + + /*Get values and lid list*/ + const int numvertices = element->GetNumberOfVertices(); + int *vertexlids = xNew(numvertices); + int *vertexsids = xNew(numvertices); + + + /*Recover vertices ids needed to initialize inputs*/ + _assert_(iomodel->elements); + for(int k=0;k(iomodel->elements[numvertices*element->Sid()+k]); //ids for vertices are in the elements array from Matlab + vertexlids[k]=iomodel->my_vertices_lids[vertexsids[k]-1]; + } + + //element->GetVerticesLidList(vertexlids); + //element->GetVerticesSidList(vertexsids); + IssmDouble* values=xNew(numvertices); + + for(int t=0;tObjectEnum()){ + case TriaEnum: transientinput->AddTriaTimeInput( t,numvertices,vertexlids,values,P1Enum); break; + case PentaEnum: transientinput->AddPentaTimeInput(t,numvertices,vertexlids,values,P1Enum); break; + default: _error_("Not implemented yet"); + } + } + xDelete(values); + xDelete(vertexlids); + xDelete(vertexsids); + } + + xDelete(times); + } + + /*Delete data:*/ + for(int i=0;i(str); + } + xDelete(pstr); + xDelete(pM); + xDelete(pN); + /*}}}*/ + /*now do the same with the dsl.sea_water_pressure_change_at_sea_floor field:{{{ */ + iomodel->FetchData(&pstr,&pM,&pN,&nummodels,"md.dsl.sea_water_pressure_change_at_sea_floor"); + + for (int i=0;i(N); + for(int t=0;tSetDatasetTransientInput(DslSeaWaterPressureChangeAtSeaFloor,i, times,N); + + for(int j=0;jSize();j++){ + + /*Get the right transient input*/ + Element* element=xDynamicCast(elements->GetObjectByOffset(j)); + + /*Get values and lid list*/ + const int numvertices = element->GetNumberOfVertices(); + int *vertexlids = xNew(numvertices); + int *vertexsids = xNew(numvertices); + + + /*Recover vertices ids needed to initialize inputs*/ + _assert_(iomodel->elements); + for(int k=0;k(iomodel->elements[numvertices*element->Sid()+k]); //ids for vertices are in the elements array from Matlab + vertexlids[k]=iomodel->my_vertices_lids[vertexsids[k]-1]; + } + //element->GetVerticesLidList(vertexlids); + //element->GetVerticesSidList(vertexsids); + + IssmDouble* values=xNew(numvertices); + + for(int t=0;tObjectEnum()){ + case TriaEnum: transientinput->AddTriaTimeInput( t,numvertices,vertexlids,values,P1Enum); break; + case PentaEnum: transientinput->AddPentaTimeInput(t,numvertices,vertexlids,values,P1Enum); break; + default: _error_("Not implemented yet"); + } + } + xDelete(values); + xDelete(vertexlids); + xDelete(vertexsids); + } + xDelete(times); + } + + /*Delete data:*/ + for(int i=0;i(str); + } + xDelete(pstr); + xDelete(pM); + xDelete(pN); + /*}}}*/ + + } /*}}}*/ else _error_("Dsl model " << dslmodel << " not implemented yet!"); @@ -157,5 +330,20 @@ parameters->AddObject(iomodel->CopyConstantObject("md.slr.geodetic",SealevelriseGeodeticEnum)); - + /*Deal with dsl multi-model ensembles: {{{*/ + iomodel->FetchData(&dslmodel,"md.dsl.model"); + parameters->AddObject(iomodel->CopyConstantObject("md.dsl.compute_fingerprints",DslComputeFingerprintsEnum)); + if(dslmodel==2){ + int modelid; + int nummodels; + + parameters->AddObject(iomodel->CopyConstantObject("md.dsl.modelid",DslModelidEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.dsl.nummodels",DslNummodelsEnum)); + iomodel->FetchData(&modelid,"md.dsl.modelid"); + iomodel->FetchData(&nummodels,"md.dsl.nummodels"); + + /*quick checks: */ + if(nummodels<=0)_error_("dslmme object in md.dsl field should contain at least 1 ensemble model!"); + if(modelid<=0 || modelid>nummodels)_error_("modelid field in dslmme object of md.dsl field should be between 1 and the number of ensemble runs!"); + } /*}}}*/ /*Deal with elasticity {{{*/ iomodel->FetchData(&elastic,"md.slr.elastic"); Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 24504) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 24505) @@ -5555,4 +5555,24 @@ void Tria::SealevelriseEustatic(Vector* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/ + /*Computational flags:*/ + int bp_compute_fingerprints= 0; + + /*some paramters first: */ + this->parameters->FindParam(&bp_compute_fingerprints,DslComputeFingerprintsEnum); + + if(!IsOceanInElement()){ + /*ok, there is ocean in this element, we should compute eustatic loads for the ocean if we have requested + *bottom pressure fingerprints:*/ + if(bp_compute_fingerprints)this->SealevelriseEustaticBottomPressure(pSgi,peustatic,latitude,longitude,radius,oceanarea,eartharea); + } + if(!IsIceInElement()){ + /*there is ice in this eleemnt, let's compute the eustatic response for ice changes:*/ + this->SealevelriseEustaticIce(pSgi,peustatic,latitude,longitude,radius,oceanarea,eartharea); + } + +} +/*}}}*/ +void Tria::SealevelriseEustaticIce(Vector* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/ + /*diverse:*/ int gsize; @@ -5583,5 +5603,5 @@ bool computeelastic= true; bool scaleoceanarea= false; - + /*early return if we are not on an ice cap:*/ if(!IsIceOnlyInElement()){ @@ -5591,5 +5611,5 @@ } - /*early return if we are fully floating: */ + /*early return if we are fully floating:*/ Input2* gr_input=this->GetInput2(MaskGroundediceLevelsetEnum); _assert_(gr_input); if (gr_input->GetInputMax()<=0){ @@ -5599,5 +5619,5 @@ } - /*If we are here, we are on ice that is fully grounded or half-way to floating: */ + /*If we are here, we are on ice that is fully grounded or half-way to floating:*/ if ((gr_input->GetInputMin())<0){ notfullygrounded=true; //used later on. @@ -5677,5 +5697,5 @@ } - /*Compute ice thickness change: */ + /*Compute ice thickness: */ Input2* deltathickness_input=this->GetInput2(SealevelriseDeltathicknessEnum); if (!deltathickness_input)_error_("delta thickness input needed to compute sea level rise!"); @@ -5740,4 +5760,158 @@ /*Add all components to the pSgi or pSgo solution vectors:*/ values[i]=3*rho_ice/rho_earth*area/eartharea*I*(G_rigid+G_elastic); + } + pSgi->SetValues(gsize,indices,values,ADD_VAL); + + /*free ressources:*/ + xDelete(values); + xDelete(indices); + } + + /*Assign output pointer:*/ + _assert_(!xIsNan(eustatic)); + _assert_(!xIsInf(eustatic)); + *peustatic=eustatic; + return; +} +/*}}}*/ +void Tria::SealevelriseEustaticBottomPressure(Vector* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/ + + /*diverse:*/ + int gsize; + bool spherical=true; + IssmDouble llr_list[NUMVERTICES][3]; + IssmDouble area; + IssmDouble I; //change in ice thickness or water level(Farrel and Clarke, Equ. 4) + IssmDouble rho; + IssmDouble late,longe,re; + IssmDouble lati,longi,ri; + + /*elastic green function:*/ + IssmDouble* G_elastic_precomputed=NULL; + int M; + + /*ice properties: */ + IssmDouble rho_water,rho_earth; + + /*constants:*/ + IssmDouble constant=0; + + /*Initialize eustatic component: do not skip this step :):*/ + IssmDouble eustatic = 0.; + + /*Computational flags:*/ + bool computerigid = true; + bool computeelastic= true; + bool scaleoceanarea= false; + bool bp_compute_fingerprints= false; + + /*we are here to compute fingerprints originating fromn bottom pressure loads:*/ + + /*Inform mask: */ + constant=1; this->AddInput2(SealevelEustaticMaskEnum,&constant,P0Enum); + + /*recover material parameters: */ + rho_water=FindParam(MaterialsRhoFreshwaterEnum); + rho_earth=FindParam(MaterialsEarthDensityEnum); + + /*recover love numbers and computational flags: */ + this->parameters->FindParam(&computerigid,SealevelriseRigidEnum); + this->parameters->FindParam(&computeelastic,SealevelriseElasticEnum); + this->parameters->FindParam(&scaleoceanarea,SealevelriseOceanAreaScalingEnum); + + /*recover elastic green function:*/ + if(computeelastic){ + DoubleVecParam* parameter = static_cast(this->parameters->FindParamObject(SealevelriseGElasticEnum)); + _assert_(parameter); + parameter->GetParameterValueByPointer(&G_elastic_precomputed,&M); + } + + /*how many dofs are we working with here? */ + this->parameters->FindParam(&gsize,MeshNumberofverticesEnum); + + /* Where is the centroid of this element?:{{{*/ + + /*retrieve coordinates: */ + ::GetVerticesCoordinates(&llr_list[0][0],this->vertices,NUMVERTICES,spherical); + + IssmDouble minlong=400; + IssmDouble maxlong=-20; + for (int i=0;imaxlong)maxlong=llr_list[i][1]; + if(llr_list[i][1]180){ + if (llr_list[0][1]==0)llr_list[0][1]=360; + if (llr_list[1][1]==0)llr_list[1][1]=360; + if (llr_list[2][1]==0)llr_list[2][1]=360; + } + + // correction at the north pole + if(llr_list[0][0]==0)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0; + if(llr_list[1][0]==0)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0; + if(llr_list[2][0]==0)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0; + + //correction at the south pole + if(llr_list[0][0]==180)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0; + if(llr_list[1][0]==180)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0; + if(llr_list[2][0]==180)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0; + + late=(llr_list[0][0]+llr_list[1][0]+llr_list[2][0])/3.0; + longe=(llr_list[0][1]+llr_list[1][1]+llr_list[2][1])/3.0; + + late=90-late; + if(longe>180)longe=(longe-180)-180; + + late=late/180*PI; + longe=longe/180*PI; + /*}}}*/ + + /*Compute area of element. For now, we dont do partially grounded elements:*/ + area=GetAreaSpherical(); + + /*Compute bottom pressure change: */ + Input2* bottompressure_change_input=this->GetInput2(DslSeaWaterPressureChangeAtSeaFloor); + if (!bottompressure_change_input)_error_("bottom pressure input needed to compute sea level rise fingerprint!"); + + /*If we are fully grounded, take the average over the element: */ + bottompressure_change_input->GetInputAverage(&I); + + /*Compute eustatic compoent:*/ + _assert_(oceanarea>0.); + if(scaleoceanarea) oceanarea=3.619e+14; // use true ocean area, m^2 + + /*We do not need to add the bottom pressure component to the eustatic value: */ + eustatic += 0; + + if(computeelastic | computerigid){ + int* indices=xNew(gsize); + IssmDouble* values=xNew(gsize); + IssmDouble alpha; + IssmDouble delPhi,delLambda; + for(int i=0;i(alpha/PI*reCast(M-1)); + G_elastic += G_elastic_precomputed[index]; + } + + /*Add all components to the pSgi or pSgo solution vectors:*/ + values[i]=3*rho_water/rho_earth*area/eartharea*I*(G_rigid+G_elastic); } pSgi->SetValues(gsize,indices,values,ADD_VAL); Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 24504) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 24505) @@ -167,4 +167,6 @@ void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea); void SealevelriseEustatic(Vector* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea); + void SealevelriseEustaticIce(Vector* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea); + void SealevelriseEustaticBottomPressure(Vector* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea); void SealevelriseNonEustatic(Vector* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea); void SealevelriseGeodetic(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz); Index: /issm/trunk-jpl/src/c/cores/sealevelrise_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/sealevelrise_core.cpp (revision 24504) +++ /issm/trunk-jpl/src/c/cores/sealevelrise_core.cpp (revision 24505) @@ -8,4 +8,5 @@ #include "../classes/Inputs2/TriaInput2.h" #include "../classes/Inputs2/TransientInput2.h" +#include "../classes/Inputs2/DatasetInput2.h" #include "../shared/shared.h" #include "../modules/modules.h" @@ -583,4 +584,23 @@ femmodel->inputs2->AddInput(tria_input_copy); } + else if(dslmodel==2){ + + int modelid; + + /*Recover modelid:*/ + femmodel->parameters->FindParam(&modelid,DslModelidEnum); + modelid--; //from matlab. + + /*find the DslSeaSurfaceHeightChangeAboveGeoidEnum dataset of transient inputs:*/ + DatasetInput2* dataset_input=femmodel->inputs2->GetDatasetInput2(DslSeaSurfaceHeightChangeAboveGeoidEnum); + + /*Go find the modelid'th transient input:*/ + TriaInput2* tria_input=dataset_input->GetTriaInputByOffset(modelid); + + /*Plug into DslDynamicRate input: */ + Input2* tria_input_copy=tria_input->copy(); + tria_input_copy->ChangeEnum(DslDynamicRateEnum); + femmodel->inputs2->AddInput(tria_input_copy); + } else _error_("not implemented yet"); @@ -603,4 +623,23 @@ TransientInput2* transient_input = femmodel->inputs2->GetTransientInput(DslGlobalAverageThermostericSeaLevelChangeEnum); TriaInput2* tria_input=transient_input->GetTriaInput(time); + Input2* tria_input_copy=tria_input->copy(); + tria_input_copy->ChangeEnum(DslStericRateEnum); + femmodel->inputs2->AddInput(tria_input_copy); + } + else if (dslmodel==2){ + int modelid; + + /*Recover modelid:*/ + femmodel->parameters->FindParam(&modelid,DslModelidEnum); + + modelid--; //from matlab. + + /*find the DslGlobalAverageThermostericSeaLevelChangeEnum dataset of transient inputs:*/ + DatasetInput2* dataset_input=femmodel->inputs2->GetDatasetInput2(DslGlobalAverageThermostericSeaLevelChangeEnum); + + /*Go find the modelid'th transient input:*/ + TriaInput2* tria_input=dataset_input->GetTriaInputByOffset(modelid); + + /*Plug into DslStericRate input: */ Input2* tria_input_copy=tria_input->copy(); tria_input_copy->ChangeEnum(DslStericRateEnum); Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 24504) +++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 24505) @@ -130,4 +130,7 @@ syn keyword cConstant DomainTypeEnum syn keyword cConstant DslModelEnum +syn keyword cConstant DslModelidEnum +syn keyword cConstant DslNummodelsEnum +syn keyword cConstant DslComputeFingerprintsEnum syn keyword cConstant EarthIdEnum syn keyword cConstant EplZigZagCounterEnum Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 24504) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 24505) @@ -124,4 +124,7 @@ DomainTypeEnum, DslModelEnum, + DslModelidEnum, + DslNummodelsEnum, + DslComputeFingerprintsEnum, EarthIdEnum, EplZigZagCounterEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 24504) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 24505) @@ -132,4 +132,7 @@ case DomainTypeEnum : return "DomainType"; case DslModelEnum : return "DslModel"; + case DslModelidEnum : return "DslModelid"; + case DslNummodelsEnum : return "DslNummodels"; + case DslComputeFingerprintsEnum : return "DslComputeFingerprints"; case EarthIdEnum : return "EarthId"; case EplZigZagCounterEnum : return "EplZigZagCounter"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 24504) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 24505) @@ -132,13 +132,16 @@ else 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; else if (strcmp(name,"EarthId")==0) return EarthIdEnum; else if (strcmp(name,"EplZigZagCounter")==0) return EplZigZagCounterEnum; - else if (strcmp(name,"EsaHElastic")==0) return EsaHElasticEnum; - else if (strcmp(name,"EsaHemisphere")==0) return EsaHemisphereEnum; - else if (strcmp(name,"EsaRequestedOutputs")==0) return EsaRequestedOutputsEnum; else stage=2; } if(stage==2){ - if (strcmp(name,"EsaUElastic")==0) return EsaUElasticEnum; + if (strcmp(name,"EsaHElastic")==0) return EsaHElasticEnum; + else if (strcmp(name,"EsaHemisphere")==0) return EsaHemisphereEnum; + else if (strcmp(name,"EsaRequestedOutputs")==0) return EsaRequestedOutputsEnum; + else if (strcmp(name,"EsaUElastic")==0) return EsaUElasticEnum; else if (strcmp(name,"ExtrapolationVariable")==0) return ExtrapolationVariableEnum; else if (strcmp(name,"FemModelComm")==0) return FemModelCommEnum; @@ -257,11 +260,11 @@ else if (strcmp(name,"MaterialsLatentheat")==0) return MaterialsLatentheatEnum; else if (strcmp(name,"MaterialsLithosphereDensity")==0) return MaterialsLithosphereDensityEnum; - else if (strcmp(name,"MaterialsLithosphereShearModulus")==0) return MaterialsLithosphereShearModulusEnum; - else if (strcmp(name,"MaterialsMantleDensity")==0) return MaterialsMantleDensityEnum; - else if (strcmp(name,"MaterialsMantleShearModulus")==0) return MaterialsMantleShearModulusEnum; else stage=3; } if(stage==3){ - if (strcmp(name,"MaterialsMeltingpoint")==0) return MaterialsMeltingpointEnum; + if (strcmp(name,"MaterialsLithosphereShearModulus")==0) return MaterialsLithosphereShearModulusEnum; + else if (strcmp(name,"MaterialsMantleDensity")==0) return MaterialsMantleDensityEnum; + else if (strcmp(name,"MaterialsMantleShearModulus")==0) return MaterialsMantleShearModulusEnum; + else if (strcmp(name,"MaterialsMeltingpoint")==0) return MaterialsMeltingpointEnum; else if (strcmp(name,"MaterialsMixedLayerCapacity")==0) return MaterialsMixedLayerCapacityEnum; else if (strcmp(name,"MaterialsMuWater")==0) return MaterialsMuWaterEnum; @@ -380,11 +383,11 @@ else if (strcmp(name,"SmbRunoffalti")==0) return SmbRunoffaltiEnum; else if (strcmp(name,"SmbRunoffgrad")==0) return SmbRunoffgradEnum; - else if (strcmp(name,"SmbRunoffref")==0) return SmbRunoffrefEnum; - else if (strcmp(name,"SmbSealev")==0) return SmbSealevEnum; - else if (strcmp(name,"SmbStepsPerStep")==0) return SmbStepsPerStepEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"SmbSwIdx")==0) return SmbSwIdxEnum; + if (strcmp(name,"SmbRunoffref")==0) return SmbRunoffrefEnum; + else if (strcmp(name,"SmbSealev")==0) return SmbSealevEnum; + else if (strcmp(name,"SmbStepsPerStep")==0) return SmbStepsPerStepEnum; + else if (strcmp(name,"SmbSwIdx")==0) return SmbSwIdxEnum; else if (strcmp(name,"SmbT0dry")==0) return SmbT0dryEnum; else if (strcmp(name,"SmbT0wet")==0) return SmbT0wetEnum; @@ -503,11 +506,11 @@ else if (strcmp(name,"CalvingrateyAverage")==0) return CalvingrateyAverageEnum; else if (strcmp(name,"Calvingratey")==0) return CalvingrateyEnum; - else if (strcmp(name,"CalvingFluxLevelset")==0) return CalvingFluxLevelsetEnum; - else if (strcmp(name,"CalvingMeltingFluxLevelset")==0) return CalvingMeltingFluxLevelsetEnum; - else if (strcmp(name,"Converged")==0) return ConvergedEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"CrevasseDepth")==0) return CrevasseDepthEnum; + if (strcmp(name,"CalvingFluxLevelset")==0) return CalvingFluxLevelsetEnum; + else if (strcmp(name,"CalvingMeltingFluxLevelset")==0) return CalvingMeltingFluxLevelsetEnum; + else if (strcmp(name,"Converged")==0) return ConvergedEnum; + else if (strcmp(name,"CrevasseDepth")==0) return CrevasseDepthEnum; else if (strcmp(name,"DamageD")==0) return DamageDEnum; else if (strcmp(name,"DamageDOld")==0) return DamageDOldEnum; @@ -626,11 +629,11 @@ else if (strcmp(name,"LevelsetfunctionSlopeX")==0) return LevelsetfunctionSlopeXEnum; else if (strcmp(name,"LevelsetfunctionSlopeY")==0) return LevelsetfunctionSlopeYEnum; - else if (strcmp(name,"LoadingforceX")==0) return LoadingforceXEnum; - else if (strcmp(name,"LoadingforceY")==0) return LoadingforceYEnum; - else if (strcmp(name,"LoadingforceZ")==0) return LoadingforceZEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"MaskGroundediceLevelset")==0) return MaskGroundediceLevelsetEnum; + if (strcmp(name,"LoadingforceX")==0) return LoadingforceXEnum; + else if (strcmp(name,"LoadingforceY")==0) return LoadingforceYEnum; + else if (strcmp(name,"LoadingforceZ")==0) return LoadingforceZEnum; + else if (strcmp(name,"MaskGroundediceLevelset")==0) return MaskGroundediceLevelsetEnum; else if (strcmp(name,"MaskIceLevelset")==0) return MaskIceLevelsetEnum; else if (strcmp(name,"MaskLandLevelset")==0) return MaskLandLevelsetEnum; @@ -749,11 +752,11 @@ else if (strcmp(name,"SmbMonthlytemperatures")==0) return SmbMonthlytemperaturesEnum; else if (strcmp(name,"SmbNetLW")==0) return SmbNetLWEnum; - else if (strcmp(name,"SmbNetSW")==0) return SmbNetSWEnum; - else if (strcmp(name,"SmbPAir")==0) return SmbPAirEnum; - else if (strcmp(name,"SmbP")==0) return SmbPEnum; else stage=7; } if(stage==7){ - if (strcmp(name,"SmbPddfacIce")==0) return SmbPddfacIceEnum; + if (strcmp(name,"SmbNetSW")==0) return SmbNetSWEnum; + else if (strcmp(name,"SmbPAir")==0) return SmbPAirEnum; + else if (strcmp(name,"SmbP")==0) return SmbPEnum; + else if (strcmp(name,"SmbPddfacIce")==0) return SmbPddfacIceEnum; else if (strcmp(name,"SmbPddfacSnow")==0) return SmbPddfacSnowEnum; else if (strcmp(name,"SmbPrecipitation")==0) return SmbPrecipitationEnum; @@ -872,11 +875,11 @@ else if (strcmp(name,"Outputdefinition23")==0) return Outputdefinition23Enum; else if (strcmp(name,"Outputdefinition24")==0) return Outputdefinition24Enum; - else if (strcmp(name,"Outputdefinition25")==0) return Outputdefinition25Enum; - else if (strcmp(name,"Outputdefinition26")==0) return Outputdefinition26Enum; - else if (strcmp(name,"Outputdefinition27")==0) return Outputdefinition27Enum; else stage=8; } if(stage==8){ - if (strcmp(name,"Outputdefinition28")==0) return Outputdefinition28Enum; + if (strcmp(name,"Outputdefinition25")==0) return Outputdefinition25Enum; + else if (strcmp(name,"Outputdefinition26")==0) return Outputdefinition26Enum; + else if (strcmp(name,"Outputdefinition27")==0) return Outputdefinition27Enum; + else if (strcmp(name,"Outputdefinition28")==0) return Outputdefinition28Enum; else if (strcmp(name,"Outputdefinition29")==0) return Outputdefinition29Enum; else if (strcmp(name,"Outputdefinition2")==0) return Outputdefinition2Enum; @@ -995,11 +998,11 @@ else if (strcmp(name,"CalvingHab")==0) return CalvingHabEnum; else if (strcmp(name,"CalvingLevermann")==0) return CalvingLevermannEnum; - else if (strcmp(name,"CalvingVonmises")==0) return CalvingVonmisesEnum; - else if (strcmp(name,"Cfdragcoeffabsgrad")==0) return CfdragcoeffabsgradEnum; - else if (strcmp(name,"Cfsurfacelogvel")==0) return CfsurfacelogvelEnum; else stage=9; } if(stage==9){ - if (strcmp(name,"Cfsurfacesquare")==0) return CfsurfacesquareEnum; + if (strcmp(name,"CalvingVonmises")==0) return CalvingVonmisesEnum; + else if (strcmp(name,"Cfdragcoeffabsgrad")==0) return CfdragcoeffabsgradEnum; + else if (strcmp(name,"Cfsurfacelogvel")==0) return CfsurfacelogvelEnum; + else if (strcmp(name,"Cfsurfacesquare")==0) return CfsurfacesquareEnum; else if (strcmp(name,"Channel")==0) return ChannelEnum; else if (strcmp(name,"ChannelArea")==0) return ChannelAreaEnum; @@ -1118,11 +1121,11 @@ else if (strcmp(name,"ElementInput2")==0) return ElementInput2Enum; else if (strcmp(name,"SegInput2")==0) return SegInput2Enum; - else if (strcmp(name,"TriaInput2")==0) return TriaInput2Enum; - else if (strcmp(name,"PentaInput2")==0) return PentaInput2Enum; - else if (strcmp(name,"IntMatExternalResult")==0) return IntMatExternalResultEnum; else stage=10; } if(stage==10){ - if (strcmp(name,"IntMatParam")==0) return IntMatParamEnum; + if (strcmp(name,"TriaInput2")==0) return TriaInput2Enum; + else if (strcmp(name,"PentaInput2")==0) return PentaInput2Enum; + else if (strcmp(name,"IntMatExternalResult")==0) return IntMatExternalResultEnum; + else if (strcmp(name,"IntMatParam")==0) return IntMatParamEnum; else if (strcmp(name,"IntParam")==0) return IntParamEnum; else if (strcmp(name,"IntVecParam")==0) return IntVecParamEnum; @@ -1241,11 +1244,11 @@ else if (strcmp(name,"SMBd18opdd")==0) return SMBd18opddEnum; else if (strcmp(name,"SMBforcing")==0) return SMBforcingEnum; - else if (strcmp(name,"SMBgcm")==0) return SMBgcmEnum; - else if (strcmp(name,"SMBgemb")==0) return SMBgembEnum; - else if (strcmp(name,"SMBgradients")==0) return SMBgradientsEnum; else stage=11; } if(stage==11){ - if (strcmp(name,"SMBgradientscomponents")==0) return SMBgradientscomponentsEnum; + if (strcmp(name,"SMBgcm")==0) return SMBgcmEnum; + else if (strcmp(name,"SMBgemb")==0) return SMBgembEnum; + else if (strcmp(name,"SMBgradients")==0) return SMBgradientsEnum; + else if (strcmp(name,"SMBgradientscomponents")==0) return SMBgradientscomponentsEnum; else if (strcmp(name,"SMBgradientsela")==0) return SMBgradientselaEnum; else if (strcmp(name,"SMBhenning")==0) return SMBhenningEnum; Index: /issm/trunk-jpl/src/m/classes/dsl.m =================================================================== --- /issm/trunk-jpl/src/m/classes/dsl.m (revision 24504) +++ /issm/trunk-jpl/src/m/classes/dsl.m (revision 24505) @@ -10,4 +10,5 @@ sea_surface_height_change_above_geoid; %corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable rate (mm/yr) sea_water_pressure_change_at_sea_floor; %corresponds to bpo field in CMIP5 archives. Specified as a spatio-temporally variable rate (in mm/yr equivalent, not in Pa/yr!) for each ensemble + compute_fingerprints; %do we use the sea water pressure change to compute fingerprints and correct sea_surface_height_change_above_geoid end @@ -32,4 +33,5 @@ self.sea_surface_height_change_above_geoid=NaN; self.sea_water_pressure_change_at_sea_floor=NaN; + self.compute_fingerprints=0; end % }}} @@ -38,7 +40,15 @@ %Early return if ~ismember('SealevelriseAnalysis',analyses), return; end + if (strcmp(solution,'TransientSolution') & md.transient.isslr == 0), return; end md = checkfield(md,'fieldname','dsl.global_average_thermosteric_sea_level_change','NaN',1,'Inf',1); md = checkfield(md,'fieldname','dsl.sea_surface_height_change_above_geoid','NaN',1,'Inf',1,'timeseries',1); md = checkfield(md,'fieldname','dsl.sea_water_pressure_change_at_sea_floor','NaN',1,'Inf',1,'timeseries',1); + md = checkfield(md,'fieldname','dsl.compute_fingerprints','NaN',1,'Inf',1,'values',[0,1]); + if self.compute_fingerprints, + %check geodetic flag of slr is on: + if md.slr.geodetic==0, + error('DSL checkconsistency error message: if bottom pressure fingerprints computations are requested, slr class should have geodetic flag on'); + end + end end % }}} @@ -49,4 +59,5 @@ fielddisplay(self,'sea_surface_height_change_above_geoid','corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable rate (mm/yr)'); fielddisplay(self,'sea_water_pressure_change_at_sea_floor','corresponds to bpo field in CMIP5 archives. Specified as a spatio-temporally variable rate (in Pa/yr)'); + fielddisplay(self,'compute_fingerprints','%do we use the sea water pressure change to compute fingerprints and correct sea_surface_height_change_above_geoid'); end % }}} @@ -54,4 +65,5 @@ WriteData(fid,prefix,'name','md.dsl.model','data',1,'format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','compute_fingerprints','format','Integer'); WriteData(fid,prefix,'object',self,'fieldname','global_average_thermosteric_sea_level_change','format','DoubleMat','mattype',1,'timeseries',1,'timeserieslength',2,'yts',md.constants.yts,'scale',1e-3/md.constants.yts); WriteData(fid,prefix,'object',self,'fieldname','sea_water_pressure_change_at_sea_floor','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts,'scale',1e-3/md.constants.yts); @@ -62,4 +74,5 @@ writejs1Darray(fid,[modelname '.dsl.global_average_thermosteric_sea_level_change'],self.global_average_thermosteric_sea_level_change); + writejs1Darray(fid,[modelname '.dsl.compute_fingerprints'],self.compute_fingerprints); writejs1Darray(fid,[modelname '.dsl.sea_surface_height_change_above_geoid'],self.sea_surface_height_change_above_geoid); writejs1Darray(fid,[modelname '.dsl.sea_water_pressure_change_at_sea_floor'],self.sea_water_pressure_change_at_sea_floor); Index: /issm/trunk-jpl/src/m/classes/dslmme.m =================================================================== --- /issm/trunk-jpl/src/m/classes/dslmme.m (revision 24505) +++ /issm/trunk-jpl/src/m/classes/dslmme.m (revision 24505) @@ -0,0 +1,92 @@ +%DSLMME class definition +% +% Usage: +% dsl=dslmme(); %dynamic sea level class based on a multi-model ensemble of CMIP5 outputs + +classdef dslmme + properties (SetAccess=public) + + modelid; %index into the multi-model ensemble, determine which field will be used. + global_average_thermosteric_sea_level_change; %corresponds to zostoga fields in CMIP5 archives. Specified as a temporally variable global rate (mm/yr) for each ensemble. + sea_surface_height_change_above_geoid; %corresponds to zos fields in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable rate (mm/yr) for each ensemble + sea_water_pressure_change_at_sea_floor; %corresponds to bpo fields in CMIP5 archives. Specified as a spatio-temporally variable rate (in mm/yr equivalent, not in Pa/yr!) for each ensemble + compute_fingerprints; %do we use the sea water pressure change to compute fingerprints and correct sea_surface_height_change_above_geoid + + end + methods + function self = extrude(self,md) % {{{ + for i=1:length(self.global_average_thermosteric_sea_level_change), + self.sea_surface_height_change_above_geoid{i}=project3d(md,'vector',self.sea_surface_height_change_above_geoid{i},'type','node','layer',1); + self.sea_water_pressure_change_at_sea_floor{i}=project3d(md,'vector',self.sea_water_pressure_change_at_sea_floor{i},'type','node','layer',1); + end + end % }}} + function self = dslmme(varargin) % {{{ + switch nargin + case 0 + self=setdefaultparameters(self); + case 1 + self=structtoobj(dsl(),varargin{1}); + otherwise + error('constructor not supported'); + end + end % }}} + function self = setdefaultparameters(self) % {{{ + + self.modelid=0; + self.global_average_thermosteric_sea_level_change={}; + self.sea_surface_height_change_above_geoid={}; + self.sea_water_pressure_change_at_sea_floor={}; + self.compute_fingerprints=0; + + end % }}} + function md = checkconsistency(self,md,solution,analyses) % {{{ + + %Early return + if ~ismember('SealevelriseAnalysis',analyses), return; end + if (strcmp(solution,'TransientSolution') & md.transient.isslr == 0), return; end + for i=1:length(self.global_average_thermosteric_sea_level_change), + md = checkfield(md,'field',self.global_average_thermosteric_sea_level_change{i},'NaN',1,'Inf',1); + md = checkfield(md,'field',self.sea_surface_height_change_above_geoid{i},'NaN',1,'Inf',1,'timeseries',1); + md = checkfield(md,'field',self.sea_water_pressure_change_at_sea_floor{i},'NaN',1,'Inf',1,'timeseries',1); + end + md = checkfield(md,'field',self.modelid,'NaN',1,'Inf',1,'>=',1,'<=',length(self.global_average_thermosteric_sea_level_change)); + if self.compute_fingerprints, + %check geodetic flag of slr is on: + if md.slr.geodetic==0, + error('DSL checkconsistency error message: if bottom pressure fingerprints computations are requested, slr class should have geodetic flag on'); + end + end + + + end % }}} + function disp(self) % {{{ + + disp(sprintf(' dsl mme parameters:')); + fielddisplay(self,'modelid','index into the multi-model ensemble, determine which field will be used.'); + fielddisplay(self,'global_average_thermosteric_sea_level_change','corresponds to zostoga fields in CMIP5 archives. Specified as a temporally variable global rate (mm/yr) for each ensemble.'); + fielddisplay(self,'sea_surface_height_change_above_geoid','corresponds to zos fields in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable rate (mm/yr) for each ensemble.'); + fielddisplay(self,'sea_water_pressure_change_at_sea_floor','corresponds to bpo fields in CMIP5 archives. Specified as a spatio-temporally variable rate (in Pa/yr) for each ensemble.'); + fielddisplay(self,'compute_fingerprints','%do we use the sea water pressure change to compute fingerprints and correct sea_surface_height_change_above_geoid'); + end % }}} + function marshall(self,prefix,md,fid) % {{{ + + WriteData(fid,prefix,'name','md.dsl.model','data',2,'format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','compute_fingerprints','format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','modelid','format','Integer'); + WriteData(fid,prefix,'name','md.dsl.nummodels','data',length(self.global_average_thermosteric_sea_level_change),'format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','global_average_thermosteric_sea_level_change','format','MatArray','timeseries',1,'timeserieslength',2,'yts',md.constants.yts,'scale',1e-3/md.constants.yts); + WriteData(fid,prefix,'object',self,'fieldname','sea_water_pressure_change_at_sea_floor','format','MatArray','timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts,'scale',1e-3/md.constants.yts); + WriteData(fid,prefix,'object',self,'fieldname','sea_surface_height_change_above_geoid','format','MatArray','timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts,'scale',1e-3/md.constants.yts); + + end % }}} + function savemodeljs(self,fid,modelname) % {{{ + + writejsdouble(fid,[modelname '.dsl.modelid'],self.modelid); + writejscellarray(fid,[modelname '.dsl.global_average_thermosteric_sea_level_change'],self.global_average_thermosteric_sea_level_change); + writejscellarray(fid,[modelname '.dsl.sea_surface_height_change_above_geoid'],self.sea_surface_height_change_above_geoid); + writejs1Darray(fid,[modelname '.dsl.compute_fingerprints'],self.compute_fingerprints); + writejscellarray(fid,[modelname '.dsl.sea_water_pressure_change_at_sea_floor'],self.sea_water_pressure_change_at_sea_floor); + + end % }}} + end +end