Index: /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 23796) +++ /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 23797) @@ -191,5 +191,5 @@ Element* element=xDynamicCast(elements->GetObjectByOffset(i)); for(int kk=0;kkDatasetInputAdd(BasalforcingsIsmp6TfEnum,array3d[kk],iomodel,Ms[kk],Ns[kk],1,BasalforcingsIsmp6TfEnum,7,BasalforcingsIsmp6TfEnum); + element->DatasetInputAdd(BasalforcingsIsmp6TfEnum,array3d[kk],iomodel,Ms[kk],Ns[kk],1,BasalforcingsIsmp6TfEnum,7,kk); } } @@ -536,9 +536,9 @@ element->FindParam(&melt_style,GroundinglineMeltInterpolationEnum); element->FindParam(&dt,TimesteppingTimeStepEnum); - Input* gmb_input = element->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(gmb_input); - Input* fmb_input = element->GetInput(BasalforcingsFloatingiceMeltingRateEnum); _assert_(fmb_input); - Input* gllevelset_input = element->GetInput(MaskGroundediceLevelsetEnum); _assert_(gllevelset_input); - Input* ms_input = element->GetInput(SmbMassBalanceEnum); _assert_(ms_input); - Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* gmb_input = element->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(gmb_input); + Input* fmb_input = element->GetInput(BasalforcingsFloatingiceMeltingRateEnum); _assert_(fmb_input); + Input* gllevelset_input = element->GetInput(MaskGroundediceLevelsetEnum); _assert_(gllevelset_input); + Input* ms_input = element->GetInput(SmbMassBalanceEnum); _assert_(ms_input); + Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input); /*Recover portion of element that is grounded*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 23796) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 23797) @@ -1901,9 +1901,9 @@ } /*}}}*/ -void Element::DatasetInputAdd(int enum_type,IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code,int input_enum){/*{{{*/ +void Element::DatasetInputAdd(int enum_type,IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int input_enum,int code,int input_id){/*{{{*/ /*enum_type: the name of the DatasetInput (eg Outputdefinition1) * vector: information being stored (eg observations) * vector_type: is if by element or by vertex - * vector_enum: is the name of the vector being stored + * input_enum: is the name of the vector being stored * code: what type of data is in the vector (booleans, ints, doubles) */ @@ -1942,7 +1942,7 @@ values[0]=vector[0]; switch(this->ObjectEnum()){ - case TriaEnum: datasetinput->AddInput(new TriaInput(vector_enum,values,P0Enum),input_enum); break; - case PentaEnum: datasetinput->AddInput(new PentaInput(vector_enum,values,P0Enum),input_enum); break; - case TetraEnum: datasetinput->AddInput(new TetraInput(vector_enum,values,P0Enum),input_enum); break; + case TriaEnum: datasetinput->AddInput(new TriaInput(input_enum,values,P0Enum),input_id); break; + case PentaEnum: datasetinput->AddInput(new PentaInput(input_enum,values,P0Enum),input_id); break; + case TetraEnum: datasetinput->AddInput(new TetraInput(input_enum,values,P0Enum),input_id); break; default: _error_("Not implemented yet"); } @@ -1951,7 +1951,7 @@ for(i=0;iObjectEnum()){ - case TriaEnum: datasetinput->AddInput(new TriaInput(vector_enum,values,P1Enum),input_enum); break; - case PentaEnum: datasetinput->AddInput(new PentaInput(vector_enum,values,P1Enum),input_enum); break; - case TetraEnum: datasetinput->AddInput(new TetraInput(vector_enum,values,P1Enum),input_enum); break; + case TriaEnum: datasetinput->AddInput(new TriaInput(input_enum,values,P1Enum),input_id); break; + case PentaEnum: datasetinput->AddInput(new PentaInput(input_enum,values,P1Enum),input_id); break; + case TetraEnum: datasetinput->AddInput(new TetraInput(input_enum,values,P1Enum),input_id); break; default: _error_("Not implemented yet"); } } @@ -1960,15 +1960,15 @@ IssmDouble* times = xNew(N); for(t=0;tObjectEnum()){ - case TriaEnum: transientinput->AddTimeInput(new TriaInput( vector_enum,values,P1Enum)); break; - case PentaEnum: transientinput->AddTimeInput(new PentaInput(vector_enum,values,P1Enum)); break; - case TetraEnum: transientinput->AddTimeInput(new TetraInput(vector_enum,values,P1Enum)); break; + case TriaEnum: transientinput->AddTimeInput(new TriaInput( input_enum,values,P1Enum)); break; + case PentaEnum: transientinput->AddTimeInput(new PentaInput(input_enum,values,P1Enum)); break; + case TetraEnum: transientinput->AddTimeInput(new TetraInput(input_enum,values,P1Enum)); break; default: _error_("Not implemented yet"); } } - datasetinput->AddInput(transientinput,input_enum); + datasetinput->AddInput(transientinput,input_id); xDelete(times); } @@ -1982,7 +1982,7 @@ if (N==this->GetNumberOfNodes(P1Enum) ){ switch(this->ObjectEnum()){ - case TriaEnum: datasetinput->AddInput(new TriaInput(vector_enum,values,P1Enum),input_enum); break; - case PentaEnum: datasetinput->AddInput(new PentaInput(vector_enum,values,P1Enum),input_enum); break; - case TetraEnum: datasetinput->AddInput(new TetraInput(vector_enum,values,P1Enum),input_enum); break; + case TriaEnum: datasetinput->AddInput(new TriaInput(input_enum,values,P1Enum),input_id); break; + case PentaEnum: datasetinput->AddInput(new PentaInput(input_enum,values,P1Enum),input_id); break; + case TetraEnum: datasetinput->AddInput(new TetraInput(input_enum,values,P1Enum),input_id); break; default: _error_("Not implemented yet"); } @@ -1990,7 +1990,7 @@ else if(N==this->GetNumberOfNodes(P0Enum) ){ switch(this->ObjectEnum()){ - case TriaEnum: datasetinput->AddInput(new TriaInput(vector_enum,values,P0Enum),input_enum); break; - case PentaEnum: datasetinput->AddInput(new PentaInput(vector_enum,values,P0Enum),input_enum); break; - case TetraEnum: datasetinput->AddInput(new TetraInput(vector_enum,values,P0Enum),input_enum); break; + case TriaEnum: datasetinput->AddInput(new TriaInput(input_enum,values,P0Enum),input_id); break; + case PentaEnum: datasetinput->AddInput(new PentaInput(input_enum,values,P0Enum),input_id); break; + case TetraEnum: datasetinput->AddInput(new TetraInput(input_enum,values,P0Enum),input_id); break; default: _error_("Not implemented yet"); } @@ -1998,7 +1998,7 @@ else if(N==this->GetNumberOfNodes(P1xP2Enum)){ switch(this->ObjectEnum()){ - case TriaEnum: datasetinput->AddInput(new TriaInput(vector_enum,values,P1xP2Enum),input_enum); break; - case PentaEnum: datasetinput->AddInput(new PentaInput(vector_enum,values,P1xP2Enum),input_enum); break; - case TetraEnum: datasetinput->AddInput(new TetraInput(vector_enum,values,P1xP2Enum),input_enum); break; + case TriaEnum: datasetinput->AddInput(new TriaInput(input_enum,values,P1xP2Enum),input_id); break; + case PentaEnum: datasetinput->AddInput(new PentaInput(input_enum,values,P1xP2Enum),input_id); break; + case TetraEnum: datasetinput->AddInput(new TetraInput(input_enum,values,P1xP2Enum),input_id); break; default: _error_("Not implemented yet"); } @@ -2006,7 +2006,7 @@ else if(N==this->GetNumberOfNodes(P1xP3Enum)) { switch(this->ObjectEnum()){ - case TriaEnum: datasetinput->AddInput(new TriaInput(vector_enum,values,P1xP3Enum),input_enum); break; - case PentaEnum: datasetinput->AddInput(new PentaInput(vector_enum,values,P1xP3Enum),input_enum); break; - case TetraEnum: datasetinput->AddInput(new TetraInput(vector_enum,values,P1xP3Enum),input_enum); break; + case TriaEnum: datasetinput->AddInput(new TriaInput(input_enum,values,P1xP3Enum),input_id); break; + case PentaEnum: datasetinput->AddInput(new PentaInput(input_enum,values,P1xP3Enum),input_id); break; + case TetraEnum: datasetinput->AddInput(new TetraInput(input_enum,values,P1xP3Enum),input_id); break; default: _error_("Not implemented yet"); } @@ -2016,5 +2016,5 @@ } else{ - _error_("nodal vector is either numberofvertices or numberofvertices+1 long. Field provided (" << EnumToStringx(vector_enum) << ") is " << M << " long"); + _error_("nodal vector is either numberofvertices or numberofvertices+1 long. Field provided (" << EnumToStringx(input_enum) << ") is " << M << " long"); } @@ -2029,11 +2029,11 @@ if(M==iomodel->numberofelements){ if (code==5){ //boolean - datasetinput->AddInput(new BoolInput(vector_enum,reCast(vector[this->Sid()])),input_enum); + datasetinput->AddInput(new BoolInput(input_enum,reCast(vector[this->Sid()])),input_id); } else if (code==6){ //integer - datasetinput->AddInput(new IntInput(vector_enum,reCast(vector[this->Sid()])),input_enum); + datasetinput->AddInput(new IntInput(input_enum,reCast(vector[this->Sid()])),input_id); } else if (code==7){ //IssmDouble - datasetinput->AddInput(new DoubleInput(vector_enum,vector[this->Sid()]),input_enum); + datasetinput->AddInput(new DoubleInput(input_enum,vector[this->Sid()]),input_id); } else _error_("could not recognize nature of vector from code " << code); @@ -2043,19 +2043,19 @@ IssmDouble* times = xNew(N); for(t=0;tSid()+t]; switch(this->ObjectEnum()){ - case TriaEnum: transientinput->AddTimeInput(new TriaInput( vector_enum,&value,P0Enum)); break; - case PentaEnum: transientinput->AddTimeInput(new PentaInput(vector_enum,&value,P0Enum)); break; - case TetraEnum: transientinput->AddTimeInput(new TetraInput(vector_enum,&value,P0Enum)); break; + case TriaEnum: transientinput->AddTimeInput(new TriaInput( input_enum,&value,P0Enum)); break; + case PentaEnum: transientinput->AddTimeInput(new PentaInput(input_enum,&value,P0Enum)); break; + case TetraEnum: transientinput->AddTimeInput(new TetraInput(input_enum,&value,P0Enum)); break; default: _error_("Not implemented yet"); } } - datasetinput->AddInput(transientinput,input_enum); + datasetinput->AddInput(transientinput,input_id); xDelete(times); } - else _error_("element vector is either numberofelements or numberofelements+1 long. Field provided (" << EnumToStringx(vector_enum) << ") is " << M << " long"); + else _error_("element vector is either numberofelements or numberofelements+1 long. Field provided (" << EnumToStringx(input_enum) << ") is " << M << " long"); } else if(vector_type==3){ //element vector @@ -2066,9 +2066,9 @@ IssmDouble* layers = xNewZeroInit(N);; for(t=0;tSid()+t]; - DoubleArrayInput* arrayinput=new DoubleArrayInput(vector_enum,layers,N); - datasetinput->AddInput(arrayinput,input_enum); + DoubleArrayInput* arrayinput=new DoubleArrayInput(input_enum,layers,N); + datasetinput->AddInput(arrayinput,input_id); xDelete(layers); } - else _error_("element vector is either numberofelements or numberofelements+1 long. Field provided (" << EnumToStringx(vector_enum) << ") is " << M << " long"); + else _error_("element vector is either numberofelements or numberofelements+1 long. Field provided (" << EnumToStringx(input_enum) << ") is " << M << " long"); } else _error_("Cannot add input for vector type " << vector_type << " (not supported)"); Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 23796) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 23797) @@ -2415,8 +2415,11 @@ /*Get variables*/ - IssmDouble rhoi = this->FindParam(MaterialsRhoIceEnum); - IssmDouble rhow = this->FindParam(MaterialsRhoSeawaterEnum)+5; - IssmDouble lf = this->FindParam(MaterialsLatentheatEnum); - IssmDouble cp = this->FindParam(MaterialsMixedLayerCapacityEnum); + IssmDouble rhoi = this->FindParam(MaterialsRhoIceEnum); + IssmDouble rhow = this->FindParam(MaterialsRhoSeawaterEnum); + IssmDouble lf = this->FindParam(MaterialsLatentheatEnum); + IssmDouble cp = this->FindParam(MaterialsMixedLayerCapacityEnum); + + /*Hard code sea water density to be consistent with ISMIP6 documentation*/ + rhow = 1028.; /* Get parameters and inputs */ @@ -2429,6 +2432,6 @@ delta_t_basin = delta_t[basinid]; mean_tf_basin = mean_tf[basinid]; - /* Compute melt rate */ - Gauss* gauss=this->NewGauss(); + /* Compute melt rate */ + Gauss* gauss=this->NewGauss(); for(int i=0;iGaussVertex(i); @@ -2439,9 +2442,6 @@ } - cout << "meltrate = " << basalmeltrate[0] << endl; - /*Return basal melt rate*/ this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,basalmeltrate,P1Enum); - _error_("STOP HERE"); /*Cleanup and return*/ @@ -2451,7 +2451,5 @@ xDelete(depths); - _error_("not implemented yet"); - - }/*}}}*/ +}/*}}}*/ bool Tria::IsNodeOnShelfFromFlags(IssmDouble* flags){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/Inputs/DatasetInput.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/DatasetInput.cpp (revision 23796) +++ /issm/trunk-jpl/src/c/classes/Inputs/DatasetInput.cpp (revision 23797) @@ -157,5 +157,10 @@ /*Object functions*/ void DatasetInput::Configure(Parameters* parameters){/*{{{*/ - /*do nothing: */ + + for(int i=0;inumids;i++){ + Input* input=xDynamicCast(this->inputs->GetObjectByOffset(i)); + input->Configure(parameters); + } + } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Inputs/TransientInput.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/TransientInput.cpp (revision 23796) +++ /issm/trunk-jpl/src/c/classes/Inputs/TransientInput.cpp (revision 23797) @@ -256,10 +256,9 @@ int TransientInput::GetTimeInputOffset(IssmDouble time){/*{{{*/ - int found; - int offset; + int offset; /*go through the timesteps, and figure out which interval we * *fall within. Then interpolate the values on this interval: */ - found=binary_search(&offset,time,this->timesteps,this->numtimesteps); + int found=binary_search(&offset,time,this->timesteps,this->numtimesteps); if(!found) _error_("Input not found (is TransientInput sorted ?)"); Index: /issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.cpp (revision 23796) +++ /issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.cpp (revision 23797) @@ -73,31 +73,76 @@ void FloatingiceMeltingRateIsmip6x(FemModel* femmodel){/*{{{*/ + int num_basins, basinid,num_depths; + IssmDouble area, tf, base, time; + IssmDouble* tf_depths = NULL; - int num_basins, basinid; - IssmDouble area, tf, base, time; - IssmDouble tf_test[3]; + femmodel->parameters->FindParam(&num_basins,BasalforcingsIsmp6NumBasinsEnum); + femmodel->parameters->FindParam(&tf_depths,&num_depths,BasalforcingsIsmp6TfDepthsEnum); _assert_(tf_depths); - femmodel->parameters->FindParam(&num_basins,BasalforcingsIsmp6NumBasinsEnum); + /*Binary search works for vectors that are sorted in increasing order only, make depths positive*/ + for(int i=0;i(num_basins); + IssmDouble* tf_weighted_avg = xNewZeroInit(num_basins); IssmDouble* tf_weighted_avg_cpu = xNewZeroInit(num_basins); - IssmDouble* areas_summed = xNewZeroInit(num_basins); - IssmDouble* areas_summed_cpu = xNewZeroInit(num_basins); - - /*Find and save TF at each ice shelf point*/ - // element->parameters->FindParam(&time,TimeEnum); - // Input* tf_input = element->GetInput(BasalforcingsIsmp6TfEnum); _assert_(tf_input); - // tf_input->Echo(); - // Input* base_input = element->GetInput(BaseEnum); _assert_(base_input); - // Gauss* gauss=element->NewGauss(1); gauss->GaussPoint(0); - // base_input->GetInputValue(&base,gauss); - // //tf_input->GetInputValue(&tf,gauss,BasalforcingsIsmp6TfEnum); + IssmDouble* areas_summed = xNewZeroInit(num_basins); + IssmDouble* areas_summed_cpu = xNewZeroInit(num_basins); /*TEST: Set tf=2 for all ice shelf elements*/ for(int i=0;ielements->Size();i++){ - Element* element=xDynamicCast(femmodel->elements->GetObjectByOffset(i)); - if(!element->IsIceInElement() || !element->IsFloating()) continue; - for(int k=0;k<3;k++) {tf_test[k] = 2.;} + Element* element = xDynamicCast(femmodel->elements->GetObjectByOffset(i)); + int numvertices = element->GetNumberOfVertices(); + + /*Set melt to 0 if non floating*/ + if(!element->IsIceInElement() || !element->IsFloating()){ + IssmDouble* values = xNewZeroInit(numvertices); + element->AddInput(BasalforcingsFloatingiceMeltingRateEnum,values,P1Enum); + element->AddInput(BasalforcingsIsmp6TfShelfEnum,values,P1Enum); + xDelete(values); + continue; + } + + /*Get TF on all vertices*/ + IssmDouble* tf_test = xNew(numvertices); + IssmDouble* depth_vertices = xNew(numvertices); + Input* tf_input = element->GetInput(BasalforcingsIsmp6TfEnum); _assert_(tf_input); + + element->GetInputListOnVertices(&depth_vertices[0],BaseEnum); + + Gauss* gauss=element->NewGauss(); + for(int iv=0;ivGaussVertex(iv); + + /*Find out where the ice shelf base is within tf_depths*/ + IssmDouble depth = -depth_vertices[iv]; /*NOTE: make sure we are dealing with depth>0*/ + int offset; + int found=binary_search(&offset,depth,tf_depths,num_depths); + if(!found) _error_("depth not found"); + + if (offset==-1){ + /*get values for the first depth: */ + _assert_(depth<=tf_depths[0]); + tf_input->GetInputValue(&tf_test[iv],gauss,0); + } + else if(offset==num_depths-1){ + /*get values for the last time: */ + _assert_(depth>=tf_depths[num_depths-1]); + tf_input->GetInputValue(&tf_test[iv],gauss,num_depths-1); + } + else { + /*get values between two times [offset:offset+1[, Interpolate linearly*/ + _assert_(depth>=tf_depths[offset] && depthGetInputValue(&tf1,gauss,offset); + tf_input->GetInputValue(&tf2,gauss,offset+1); + tf_test[iv] = alpha1*tf1 + alpha2*tf2; + } + } + element->AddInput(BasalforcingsIsmp6TfShelfEnum,tf_test,P1Enum); + xDelete(tf_test); + xDelete(depth_vertices); } @@ -129,4 +174,5 @@ xDelete(areas_summed); xDelete(areas_summed_cpu); + xDelete(tf_depths); /*Compute meltrates*/