Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 25307) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 25308) @@ -184,5 +184,4 @@ basalelement->GetInputValue(&active_element,HydrologydcMaskEplactiveEltEnum); - /*Check that all nodes are active, else return empty matrix*/ if(!active_element) { @@ -213,21 +212,18 @@ basalelement->FindParam(&dt,TimesteppingTimeStepEnum); Input2* epl_thick_input = basalelement->GetInput2(HydrologydcEplThicknessSubstepEnum); _assert_(epl_thick_input); - Input2* epl_head_input = basalelement->GetInput2(EplHeadSubstepEnum); _assert_(epl_head_input); - Input2* base_input = basalelement->GetInput2(BaseEnum); _assert_(base_input); /* Start looping on the number of gaussian points: */ Gauss* gauss = basalelement->NewGauss(2); for(int ig=gauss->begin();igend();ig++){ - gauss ->GaussPoint(ig); - basalelement ->JacobianDeterminant(&Jdet,xyz_list,gauss); + gauss ->GaussPoint(ig); + basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss); basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); basalelement->NodalFunctions(basis,gauss); - epl_transmitivity = EplTransmitivity(basalelement,gauss,epl_thick_input,epl_head_input,base_input); - epl_storing = EplStoring(basalelement,gauss,epl_thick_input,epl_head_input,base_input); + epl_transmitivity = EplTransmitivity(basalelement,gauss,epl_thick_input); + epl_storing = EplStoring(basalelement,gauss,epl_thick_input); /*Diffusivity*/ D_scalar=epl_transmitivity*gauss->weight*Jdet; - //D_scalar=gauss->weight*Jdet; if(dt!=0.) D_scalar=D_scalar*dt; for(int i=0;iNodalFunctions(&basis[0],gauss); D_scalar=epl_storing*gauss->weight*Jdet; - //D_scalar=(epl_storing/epl_transmitivity)*gauss->weight*Jdet; for(int i=0;ivalues[i*numnodes+j] += D_scalar*basis[j]*basis[i]; - /*Transfer EPL part*/ transfer=GetHydrologyKMatrixTransfer(basalelement); D_scalar=dt*transfer*gauss->weight*Jdet; - //D_scalar=dt*(transfer/epl_transmitivity)*gauss->weight*Jdet; for(int i=0;ivalues[i*numnodes+j] += D_scalar*basis[j]*basis[i]; } @@ -291,12 +282,10 @@ } /*Intermediaries */ - int smb_model; - int smb_averaging; - int smbsubstepping; - int hydrologysubstepping; + int smb_model,smb_averaging; + int smbsubstepping,hydrologysubstepping; IssmDouble dt,scalar,water_head; IssmDouble water_load,transfer,runoff_value; IssmDouble epl_storing,epl_transmitivity; - IssmDouble Jdet; + IssmDouble Jdet,time; IssmDouble residual,connectivity; @@ -318,15 +307,9 @@ basalelement->FindParam(&dt,TimesteppingTimeStepEnum); basalelement ->FindParam(&smb_model,SmbEnum); - basalelement->FindParam(&smb_averaging,SmbAveragingEnum); Input2* epl_thick_input = basalelement->GetInput2(HydrologydcEplThicknessSubstepEnum); _assert_(epl_thick_input); Input2* sed_head_input = basalelement->GetInput2(SedimentHeadSubstepEnum); _assert_(sed_head_input); - Input2* epl_head_input = basalelement->GetInput2(EplHeadSubstepEnum); _assert_(epl_head_input); Input2* basal_melt_input = basalelement->GetInput2(BasalforcingsGroundediceMeltingRateEnum); _assert_(basal_melt_input); Input2* residual_input = basalelement->GetInput2(SedimentHeadResidualEnum); _assert_(residual_input); - Input2* base_input = basalelement->GetInput2(BaseEnum); _assert_(base_input); - - IssmDouble time; - basalelement->FindParam(&time,TimeEnum); if(dt!= 0.){ @@ -334,4 +317,5 @@ } if(smb_model==SMBgradientscomponentsEnum){ + basalelement->FindParam(&time,TimeEnum); basalelement->FindParam(&smbsubstepping,SmbStepsPerStepEnum); basalelement->FindParam(&hydrologysubstepping,HydrologyStepsPerStepEnum); @@ -347,4 +331,5 @@ else{ //finer stepping in smb, we average the runoff from transient input + basalelement->FindParam(&smb_averaging,SmbAveragingEnum); dummy_input = basalelement->GetInput2(SmbRunoffTransientEnum,time-dt,time,smb_averaging); _assert_(dummy_input); } @@ -358,6 +343,6 @@ basalelement ->JacobianDeterminant(&Jdet,xyz_list,gauss); basalelement ->NodalFunctions(basis,gauss); - epl_storing = EplStoring(basalelement,gauss,epl_thick_input,epl_head_input,base_input); - epl_transmitivity = EplTransmitivity(basalelement,gauss,epl_thick_input,epl_head_input,base_input); + epl_storing = EplStoring(basalelement,gauss,epl_thick_input); + epl_transmitivity = EplTransmitivity(basalelement,gauss,epl_thick_input); /*Loading term*/ @@ -366,5 +351,4 @@ else runoff_value = 0.; scalar = Jdet*gauss->weight*(water_load+runoff_value); - //scalar = Jdet*gauss->weight*(water_load)/epl_transmitivity; if(dt!=0.) scalar = scalar*dt; for(int i=0;ivalues[i]+=scalar*basis[i]; @@ -376,5 +360,4 @@ transfer=GetHydrologyPVectorTransfer(basalelement,gauss,sed_head_input); scalar = Jdet*gauss->weight*((water_head*epl_storing)+(dt*transfer)); - //scalar = Jdet*gauss->weight*((water_head*epl_storing)+(dt*transfer))/epl_transmitivity; for(int i=0;ivalues[i]+=scalar*basis[i]; } @@ -386,9 +369,8 @@ for(int iv=0;ivGaussVertex(iv); - epl_transmitivity = EplTransmitivity(basalelement,gauss,epl_thick_input,epl_head_input,base_input); + epl_transmitivity = EplTransmitivity(basalelement,gauss,epl_thick_input); connectivity = IssmDouble(basalelement->VertexConnectivity(iv)); residual_input->GetInputValue(&residual,gauss); pe->values[iv]+=residual/connectivity; - //pe->values[iv]+=residual/(epl_transmitivity*connectivity); } /*Clean up and return*/ @@ -453,47 +435,22 @@ /*Intermediaries*/ -IssmDouble HydrologyDCEfficientAnalysis::EplStoring(Element* element,Gauss* gauss, Input2* epl_thick_input, Input2* epl_head_input, Input2* base_input){/*{{{*/ +IssmDouble HydrologyDCEfficientAnalysis::EplStoring(Element* element,Gauss* gauss, Input2* epl_thick_input){/*{{{*/ IssmDouble epl_storing; - IssmDouble water_sheet,storing; - IssmDouble epl_thickness,prestep_head,base_elev; - IssmDouble rho_freshwater = element->FindParam(MaterialsRhoFreshwaterEnum); - IssmDouble g = element->FindParam(ConstantsGEnum); - IssmDouble epl_porosity = element->FindParam(HydrologydcEplPorosityEnum); - IssmDouble epl_compressibility = element->FindParam(HydrologydcEplCompressibilityEnum); - IssmDouble water_compressibility = element->FindParam(HydrologydcWaterCompressibilityEnum); + IssmDouble epl_thickness; + IssmDouble porewater_mass = element->FindParam(HydrologydcEplPoreWaterMassEnum); + IssmDouble layer_compressibility = element->FindParam(HydrologydcEplLayerCompressibilityEnum); epl_thick_input->GetInputValue(&epl_thickness,gauss); - epl_head_input->GetInputValue(&prestep_head,gauss); - base_input->GetInputValue(&base_elev,gauss); - water_sheet=max(0.0,(prestep_head-base_elev)); - storing=rho_freshwater*g*epl_porosity*epl_thickness*(water_compressibility+(epl_compressibility/epl_porosity)); - - /* //porosity for unconfined region */ - /* if (water_sheet<=0.9*epl_thickness){ */ - /* epl_storing=epl_porosity; */ - /* } */ - /* //continuity ramp */ - /* else if((water_sheet0.9*epl_thickness)){ */ - /* epl_storing=(epl_thickness-water_sheet)*(epl_porosity-storing)/(0.1*epl_thickness)+storing; */ - /* } */ - /* //storing coefficient for confined */ - /* else{ */ - /* epl_storing=storing; */ - /* } */ - /* return epl_storing; */ - return storing; -}/*}}}*/ -IssmDouble HydrologyDCEfficientAnalysis::EplTransmitivity(Element* element,Gauss* gauss, Input2* epl_thick_input, Input2* epl_head_input, Input2* base_input){/*{{{*/ + epl_storing=porewater_mass*epl_thickness*layer_compressibility; + + return epl_storing; +}/*}}}*/ +IssmDouble HydrologyDCEfficientAnalysis::EplTransmitivity(Element* element,Gauss* gauss, Input2* epl_thick_input){/*{{{*/ IssmDouble epl_transmitivity; - IssmDouble water_sheet; - IssmDouble epl_thickness,base_elev,prestep_head; + IssmDouble epl_thickness; IssmDouble epl_conductivity = element->FindParam(HydrologydcEplConductivityEnum); epl_thick_input->GetInputValue(&epl_thickness,gauss); - epl_head_input->GetInputValue(&prestep_head,gauss); - base_input->GetInputValue(&base_elev,gauss); - - water_sheet=max(0.0,(prestep_head-base_elev)); + epl_transmitivity=epl_conductivity*epl_thickness; - //epl_transmitivity=max(1.0e-6,(epl_conductivity*min(water_sheet,epl_thickness))); return epl_transmitivity; }/*}}}*/ @@ -584,8 +541,4 @@ int iseplthickcomp; int domaintype; - IssmDouble dt,A; - IssmDouble EPLgrad2; - IssmDouble EPL_N; - IssmDouble opening,closing; femmodel->parameters->FindParam(&domaintype,DomainTypeEnum); @@ -612,30 +565,23 @@ element->GetInputValue(&active_element,HydrologydcMaskEplactiveEltEnum); - element->FindParam(&dt,TimesteppingTimeStepEnum); - - /*For now, assuming just one way to compute EPL thickness*/ - IssmDouble gravity = element->FindParam(ConstantsGEnum); - IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum); - IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum); - IssmDouble latentheat = element->FindParam(MaterialsLatentheatEnum); - IssmDouble epl_conductivity = element->FindParam(HydrologydcEplConductivityEnum); - IssmDouble init_thick = element->FindParam(HydrologydcEplInitialThicknessEnum); - IssmDouble max_thick = element->FindParam(HydrologydcEplMaxThicknessEnum); - - switch(domaintype){ - case Domain2DhorizontalEnum: element->GetInputListOnVertices(&B[0],MaterialsRheologyBbarEnum); break; - case Domain3DEnum: element->GetInputListOnVertices(&B[0],MaterialsRheologyBEnum); break; - default: _error_("not Implemented Yet"); - } - - element->GetInputListOnVertices(&eplhead[0],EplHeadSubstepEnum); - element->GetInputListOnVertices(&epl_slopeX[0],EplHeadSlopeXEnum); - element->GetInputListOnVertices(&epl_slopeY[0],EplHeadSlopeYEnum); - element->GetInputListOnVertices(&old_thickness[0],HydrologydcEplThicknessOldEnum); - element->GetInputListOnVertices(&ice_thickness[0],ThicknessEnum); - element->GetInputListOnVertices(&bed[0],BaseEnum); - element->GetInputListOnVertices(&n[0],MaterialsRheologyNEnum); + + /*parameters*/ + IssmDouble gravity; + IssmDouble rho_water; + IssmDouble rho_ice; + IssmDouble latentheat; + IssmDouble epl_conductivity; + IssmDouble init_thick; + IssmDouble max_thick; + IssmDouble dt; + + /* Intermiedaries */ + IssmDouble A; + IssmDouble EPLgrad2; + IssmDouble EPL_N; + IssmDouble opening,closing; if(!active_element){ + init_thick = element->FindParam(HydrologydcEplInitialThicknessEnum); /*Keeping thickness to initial value if EPL is not active*/ for(int i=0;iGetInputListOnVertices(&B[0],MaterialsRheologyBbarEnum); break; + case Domain3DEnum: element->GetInputListOnVertices(&B[0],MaterialsRheologyBEnum); break; + default: _error_("not Implemented Yet"); + } + + element->FindParam(&max_thick,HydrologydcEplMaxThicknessEnum); + element->FindParam(&epl_conductivity,HydrologydcEplConductivityEnum); + element->FindParam(&latentheat,MaterialsLatentheatEnum); + element->FindParam(&rho_ice,MaterialsRhoIceEnum); + element->FindParam(&rho_water,MaterialsRhoFreshwaterEnum); + element->FindParam(&gravity,ConstantsGEnum); + element->FindParam(&dt,TimesteppingTimeStepEnum); + + element->GetInputListOnVertices(&eplhead[0],EplHeadSubstepEnum); + element->GetInputListOnVertices(&epl_slopeX[0],EplHeadSlopeXEnum); + element->GetInputListOnVertices(&epl_slopeY[0],EplHeadSlopeYEnum); + element->GetInputListOnVertices(&old_thickness[0],HydrologydcEplThicknessOldEnum); + element->GetInputListOnVertices(&ice_thickness[0],ThicknessEnum); + element->GetInputListOnVertices(&bed[0],BaseEnum); + element->GetInputListOnVertices(&n[0],MaterialsRheologyNEnum); + for(int i=0;i(numnodes); IssmDouble* residual =xNew(numnodes); - IssmDouble* base =xNew(numnodes); IssmDouble init_thick =basalelement->FindParam(HydrologydcEplInitialThicknessEnum); @@ -712,10 +679,9 @@ basalelement->GetInputValue(&active_element,HydrologydcMaskEplactiveEltEnum); + basalelement-> GetInputListOnVertices(&sedhead[0],SedimentHeadSubstepEnum); basalelement-> GetInputListOnVertices(&old_active[0],HydrologydcMaskEplactiveNodeEnum); + basalelement-> GetInputListOnVertices(&residual[0],SedimentHeadResidualEnum); basalelement-> GetInputListOnVertices(&epl_thickness[0],HydrologydcEplThicknessSubstepEnum); - basalelement-> GetInputListOnVertices(&sedhead[0],SedimentHeadSubstepEnum); basalelement-> GetInputListOnVertices(&eplhead[0],EplHeadSubstepEnum); - basalelement-> GetInputListOnVertices(&residual[0],SedimentHeadResidualEnum); - basalelement-> GetInputListOnVertices(&base[0],BaseEnum); /*Get minimum sediment head of the element*/ @@ -769,5 +735,4 @@ xDelete(eplhead); xDelete(residual); - xDelete(base); } /*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h (revision 25307) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h (revision 25308) @@ -35,6 +35,6 @@ /*Intermediaries*/ - IssmDouble EplStoring(Element* element,Gauss* gauss, Input2* epl_thick_input, Input2* epl_head_input, Input2* base_input); - IssmDouble EplTransmitivity(Element* element,Gauss* gauss, Input2* epl_thick_input, Input2* epl_head_input, Input2* base_input); + IssmDouble EplStoring(Element* element,Gauss* gauss, Input2* epl_thick_input); + IssmDouble EplTransmitivity(Element* element,Gauss* gauss, Input2* epl_thick_input); void GetHydrologyDCInefficientHmax(IssmDouble* ph_max,Element* element, Node* innode); IssmDouble GetHydrologyKMatrixTransfer(Element* element); Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 25307) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 25308) @@ -238,5 +238,4 @@ Input2* sed_head_input = basalelement->GetInput2(SedimentHeadSubstepEnum); Input2* base_input = basalelement->GetInput2(BaseEnum); - Input2* old_wh_input = basalelement->GetInput2(SedimentHeadOldEnum); _assert_(old_wh_input); /*Transfer related Inputs*/ @@ -244,11 +243,10 @@ basalelement->GetInput2Value(&active_element,HydrologydcMaskEplactiveEltEnum); } - /* Start looping on the number of gaussian points: */ Gauss* gauss=basalelement->NewGauss(2); for(int ig=gauss -> begin();igend();ig++){ - gauss -> GaussPoint(ig); - basalelement -> JacobianDeterminant(&Jdet,xyz_list,gauss); + gauss ->GaussPoint(ig); + basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss); basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); basalelement->NodalFunctions(basis,gauss); @@ -259,5 +257,4 @@ /*Diffusivity*/ D_scalar=sediment_transmitivity*gauss->weight*Jdet; - //D_scalar=gauss->weight*Jdet; if(dt!=0.) D_scalar=D_scalar*dt; for(int i=0;iweight*Jdet; - //D_scalar=(sediment_storing/sediment_transmitivity)*gauss->weight*Jdet; for(int i=0;ivalues[i*numnodes+j] += D_scalar*basis[j]*basis[i]; - /*Transfer EPL part*/ if(isefficientlayer){ @@ -278,5 +272,4 @@ transfer=GetHydrologyKMatrixTransfer(basalelement); D_scalar=dt*transfer*gauss->weight*Jdet; - //D_scalar=dt*(transfer/sediment_transmitivity)*gauss->weight*Jdet; for(int i=0;ivalues[i*numnodes+j] += D_scalar*basis[j]*basis[i]; } @@ -298,6 +291,6 @@ /*Intermediaries*/ - bool thawed_element; - int domaintype; + bool thawed_element; + int domaintype; Element* basalelement; @@ -328,12 +321,10 @@ /*Intermediaries */ bool active_element,isefficientlayer; - int smb_model; - int smbsubstepping; - int hydrologysubstepping; - int smb_averaging; + int smb_model,smbsubstepping; + int hydrologysubstepping,smb_averaging; IssmDouble dt,scalar,sediment_storing; IssmDouble water_head,sediment_transmitivity; IssmDouble water_load,runoff_value,transfer; - IssmDouble Jdet; + IssmDouble Jdet,time; IssmDouble *xyz_list = NULL; @@ -355,5 +346,4 @@ basalelement->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum); basalelement->FindParam(&smb_model,SmbEnum); - basalelement->FindParam(&smb_averaging,SmbAveragingEnum); Input2* sed_head_input = basalelement->GetInput2(SedimentHeadSubstepEnum); @@ -363,11 +353,9 @@ Input2* SedTrans_input = basalelement->GetInput2(HydrologydcSedimentTransmitivityEnum); _assert_(SedTrans_input); - IssmDouble time; - basalelement->FindParam(&time,TimeEnum); - if(dt!= 0.){ old_wh_input = basalelement->GetInput2(SedimentHeadOldEnum); _assert_(old_wh_input); } if(smb_model==SMBgradientscomponentsEnum){ + basalelement->FindParam(&time,TimeEnum); basalelement->FindParam(&smbsubstepping,SmbStepsPerStepEnum); basalelement->FindParam(&hydrologysubstepping,HydrologyStepsPerStepEnum); @@ -383,4 +371,5 @@ else{ //finer stepping in smb, we average the runoff from transient input + basalelement->FindParam(&smb_averaging,SmbAveragingEnum); dummy_input = basalelement->GetInput2(SmbRunoffTransientEnum,time-dt,time,smb_averaging); _assert_(dummy_input); } @@ -395,7 +384,4 @@ /* Start looping on the number of gaussian points: */ Gauss* gauss=basalelement->NewGauss(2); - - IssmDouble yts; - basalelement->FindParam(&yts,ConstantsYtsEnum); for(int ig=gauss->begin();igend();ig++){ @@ -411,5 +397,4 @@ else runoff_value = 0.; scalar = Jdet*gauss->weight*(water_load+runoff_value); - //scalar = Jdet*gauss->weight*(water_load)/sediment_transmitivity; if(dt!=0.) scalar = scalar*dt; for(int i=0;iweight*(water_load+runoff_value); - //scalar = Jdet*gauss->weight*(water_load)/sediment_transmitivity; if(dt!=0.) scalar = scalar*dt; for(int i=0;iweight*((water_head*sediment_storing)+(dt*transfer)); - //scalar = Jdet*gauss->weight*((water_head*sediment_storing)+(dt*transfer))/sediment_transmitivity; for(int i=0;ivalues[i]+=scalar*basis[i]; } else{ scalar = Jdet*gauss->weight*(water_head*sediment_storing); - //scalar = Jdet*gauss->weight*(water_head*sediment_storing)/sediment_transmitivity; for(int i=0;ivalues[i]+=scalar*basis[i]; } @@ -571,30 +553,21 @@ IssmDouble storing,yield; IssmDouble base_elev,prestep_head,water_sheet; - IssmDouble rho_freshwater = element->FindParam(MaterialsRhoFreshwaterEnum); - IssmDouble g = element->FindParam(ConstantsGEnum); - IssmDouble sediment_porosity = element->FindParam(HydrologydcSedimentPorosityEnum); + IssmDouble porewater_mass = element->FindParam(HydrologydcSedimentPoreWaterMassEnum); + IssmDouble layer_compressibility = element->FindParam(HydrologydcSedimentLayerCompressibilityEnum); IssmDouble sediment_thickness = element->FindParam(HydrologydcSedimentThicknessEnum); - IssmDouble sediment_compressibility = element->FindParam(HydrologydcSedimentCompressibilityEnum); - IssmDouble water_compressibility = element->FindParam(HydrologydcWaterCompressibilityEnum); element->FindParam(&unconf_scheme,HydrologydcUnconfinedFlagEnum); switch(unconf_scheme){ case 0: - sediment_storing=rho_freshwater*g*sediment_porosity*sediment_thickness*(water_compressibility+(sediment_compressibility/sediment_porosity)); + sediment_storing=porewater_mass*sediment_thickness*layer_compressibility; break; case 1: + yield = element->FindParam(HydrologydcSedimentPorosityEnum); base_input->GetInputValue(&base_elev,gauss); sed_head_input->GetInputValue(&prestep_head,gauss); + water_sheet=max(0.0,(prestep_head-(base_elev-sediment_thickness))); - - /* if (water_sheetFindParam(HydrologydcSedimentThicknessEnum); Index: /issm/trunk-jpl/src/c/classes/Inputs2/TransientInput2.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs2/TransientInput2.cpp (revision 25307) +++ /issm/trunk-jpl/src/c/classes/Inputs2/TransientInput2.cpp (revision 25308) @@ -420,5 +420,4 @@ /*If already processed return*/ if(fabs(this->current_step-this_step)<1.e-5) return; - // if(this->current_step>this_step-1.e-5 && this->current_step(start_offset)+0.5*timespan; - /*If already processed return, we set step in the middle of the interval*/ if(fabs(this->current_step-mid_step)<1.e-5) return; Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 25307) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 25308) @@ -1,5 +1,5 @@ /*!\file: CreateParameters.cpp * \brief general driver for creating parameters dataset - */ + */ #ifdef HAVE_CONFIG_H @@ -63,5 +63,5 @@ parameters->AddObject(iomodel->CopyConstantObject("md.calving.law",CalvingLawEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.frontalforcings.parameterization",FrontalForcingsParamEnum)); - parameters->AddObject(new IntParam(SealevelriseRunCountEnum,1)); + parameters->AddObject(new IntParam(SealevelriseRunCountEnum,1)); {/*This is specific to ice...*/ @@ -98,5 +98,5 @@ } - /*amr properties*/ + /*amr properties*/ int amrtype,amr_frequency; iomodel->FindConstant(&amr_frequency,"md.transient.amr_frequency"); @@ -156,5 +156,5 @@ case LinearFloatingMeltRateEnum: iomodel->FindConstant(&interp,"md.timestepping.interp_forcings"); - iomodel->FetchData(&transparam,&N,&M,"md.basalforcings.deepwater_melting_rate"); + iomodel->FetchData(&transparam,&N,&M,"md.basalforcings.deepwater_melting_rate"); if(N==1){ _assert_(M==1); @@ -166,5 +166,5 @@ } xDelete(transparam); - iomodel->FetchData(&transparam,&N,&M,"md.basalforcings.upperwater_melting_rate"); + iomodel->FetchData(&transparam,&N,&M,"md.basalforcings.upperwater_melting_rate"); if(N==1){ _assert_(M==1); @@ -176,5 +176,5 @@ } xDelete(transparam); - iomodel->FetchData(&transparam,&N,&M,"md.basalforcings.deepwater_elevation"); + iomodel->FetchData(&transparam,&N,&M,"md.basalforcings.deepwater_elevation"); if(N==1){ _assert_(M==1); @@ -183,8 +183,8 @@ else{ _assert_(N==2); - parameters->AddObject(new TransientParam(BasalforcingsDeepwaterElevationEnum,&transparam[0],&transparam[M],interp,M)); - } - xDelete(transparam); - iomodel->FetchData(&transparam,&N,&M,"md.basalforcings.upperwater_elevation"); + parameters->AddObject(new TransientParam(BasalforcingsDeepwaterElevationEnum,&transparam[0],&transparam[M],interp,M)); + } + xDelete(transparam); + iomodel->FetchData(&transparam,&N,&M,"md.basalforcings.upperwater_elevation"); if(N==1){ _assert_(M==1); @@ -226,9 +226,9 @@ parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.isplume",BasalforcingsPicoIsplumeEnum)); iomodel->FetchData(&transparam,&M,&N,"md.basalforcings.farocean_temperature"); - _assert_(M>=1 && N>=1); + _assert_(M>=1 && N>=1); parameters->AddObject(new TransientArrayParam(BasalforcingsPicoFarOceantemperatureEnum,transparam,&transparam[N*(M-1)],interp,N,M)); xDelete(transparam); iomodel->FetchData(&transparam,&M,&N,"md.basalforcings.farocean_salinity"); - _assert_(M>=1 && N>=1); + _assert_(M>=1 && N>=1); parameters->AddObject(new TransientArrayParam(BasalforcingsPicoFarOceansalinityEnum,transparam,&transparam[N*(M-1)],interp,N,M)); xDelete(transparam); @@ -237,5 +237,5 @@ parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.num_basins",BasalforcingsIsmip6NumBasinsEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.gamma_0",BasalforcingsIsmip6Gamma0Enum)); - parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.islocal",BasalforcingsIsmip6IsLocalEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.islocal",BasalforcingsIsmip6IsLocalEnum)); iomodel->FetchData(&transparam,&M,&N,"md.basalforcings.delta_t"); parameters->AddObject(new DoubleVecParam(BasalforcingsIsmip6DeltaTEnum,transparam,N)); @@ -279,6 +279,6 @@ } iomodel->FindConstant(&time,"md.timestepping.start_time"); - parameters->AddObject(new DoubleParam(TimeEnum,time)); - parameters->AddObject(new IntParam(StepEnum,0)); + parameters->AddObject(new DoubleParam(TimeEnum,time)); + parameters->AddObject(new IntParam(StepEnum,0)); /*By default, save all results*/ @@ -421,20 +421,34 @@ iomodel->FindConstant(&hydrology_model,"md.hydrology.model"); if(hydrology_model==HydrologydcEnum){ - /*FIXME: this cshould go to Analysis!!!*/ - parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.sediment_compressibility",HydrologydcSedimentCompressibilityEnum)); + IssmDouble sedcomp, sedporo, watcomp, rhofresh, g; + + /*FIXME: this should go to Analysis!!!*/ parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.sediment_porosity",HydrologydcSedimentPorosityEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.sediment_thickness",HydrologydcSedimentThicknessEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.water_compressibility",HydrologydcWaterCompressibilityEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.isefficientlayer",HydrologydcIsefficientlayerEnum)); + + iomodel->FindConstant(&sedcomp,"md.hydrology.sediment_compressibility"); + iomodel->FindConstant(&sedporo,"md.hydrology.sediment_porosity"); + iomodel->FindConstant(&watcomp,"md.hydrology.water_compressibility"); + iomodel->FindConstant(&rhofresh,"md.materials.rho_freshwater"); + iomodel->FindConstant(&g,"md.constants.g"); + + parameters->AddObject(new DoubleParam(HydrologydcSedimentLayerCompressibilityEnum,(watcomp + sedcomp/sedporo))); + parameters->AddObject(new DoubleParam(HydrologydcSedimentPoreWaterMassEnum,(rhofresh*g*sedporo))); bool isefficientlayer; iomodel->FindConstant(&isefficientlayer,"md.hydrology.isefficientlayer"); if(isefficientlayer){ - parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.epl_compressibility",HydrologydcEplCompressibilityEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.epl_porosity",HydrologydcEplPorosityEnum)); + IssmDouble eplcomp, eplporo; parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.epl_initial_thickness",HydrologydcEplInitialThicknessEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.epl_colapse_thickness",HydrologydcEplColapseThicknessEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.epl_max_thickness",HydrologydcEplMaxThicknessEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.epl_conductivity",HydrologydcEplConductivityEnum)); + + iomodel->FindConstant(&eplcomp,"md.hydrology.epl_compressibility"); + iomodel->FindConstant(&eplporo,"md.hydrology.epl_porosity"); + parameters->AddObject(new DoubleParam(HydrologydcEplLayerCompressibilityEnum,(watcomp + eplcomp/eplporo))); + parameters->AddObject(new DoubleParam(HydrologydcEplPoreWaterMassEnum,(rhofresh*g*eplporo))); + } } @@ -472,5 +486,5 @@ if(mass_flux_present){ - /*Fetch the mass flux segments necessary to compute the mass fluxes. Build a DoubleMatArrayParam object out of them: */ + /*Fetch the mass flux segments necessary to compute the mass fluxes. Build a DoubleMatArrayParam object out of them: */ iomodel->FetchData(&array,&mdims_array,&ndims_array,&mass_flux_num_profiles,"md.qmu.mass_flux_segments"); if(mass_flux_num_profiles==0)_error_("mass_flux_num_profiles is 0, when MassFlux computations were requested!"); @@ -518,5 +532,5 @@ xDelete(matrix); } - xDelete(mdims_array); + xDelete(mdims_array); xDelete(ndims_array); xDelete(array); Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 25307) +++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 25308) @@ -182,9 +182,9 @@ syn keyword cConstant HydrologyStorageEnum syn keyword cConstant HydrologydcEplColapseThicknessEnum -syn keyword cConstant HydrologydcEplCompressibilityEnum syn keyword cConstant HydrologydcEplConductivityEnum syn keyword cConstant HydrologydcEplInitialThicknessEnum +syn keyword cConstant HydrologydcEplLayerCompressibilityEnum syn keyword cConstant HydrologydcEplMaxThicknessEnum -syn keyword cConstant HydrologydcEplPorosityEnum +syn keyword cConstant HydrologydcEplPoreWaterMassEnum syn keyword cConstant HydrologydcEplThickCompEnum syn keyword cConstant HydrologydcEplflipLockEnum @@ -195,12 +195,12 @@ syn keyword cConstant HydrologydcPenaltyLockEnum syn keyword cConstant HydrologydcRelTolEnum -syn keyword cConstant HydrologydcSedimentCompressibilityEnum syn keyword cConstant HydrologydcSedimentlimitEnum syn keyword cConstant HydrologydcSedimentlimitFlagEnum +syn keyword cConstant HydrologydcSedimentLayerCompressibilityEnum +syn keyword cConstant HydrologydcSedimentPoreWaterMassEnum syn keyword cConstant HydrologydcSedimentPorosityEnum syn keyword cConstant HydrologydcSedimentThicknessEnum syn keyword cConstant HydrologydcTransferFlagEnum syn keyword cConstant HydrologydcUnconfinedFlagEnum -syn keyword cConstant HydrologydcWaterCompressibilityEnum syn keyword cConstant HydrologyshreveStabilizationEnum syn keyword cConstant IcecapToEarthCommEnum @@ -1361,4 +1361,5 @@ syn keyword cType Cfsurfacesquare syn keyword cType Channel +syn keyword cType classes syn keyword cType Constraint syn keyword cType Constraints @@ -1367,6 +1368,6 @@ syn keyword cType ControlInput2 syn keyword cType Covertree +syn keyword cType DatasetInput2 syn keyword cType DataSetParam -syn keyword cType DatasetInput2 syn keyword cType Definition syn keyword cType DependentObject @@ -1381,6 +1382,6 @@ syn keyword cType ElementInput2 syn keyword cType ElementMatrix +syn keyword cType Elements syn keyword cType ElementVector -syn keyword cType Elements syn keyword cType ExponentialVariogram syn keyword cType ExternalResult @@ -1389,9 +1390,10 @@ syn keyword cType Friction syn keyword cType Gauss +syn keyword cType GaussianVariogram +syn keyword cType gaussobjects syn keyword cType GaussPenta syn keyword cType GaussSeg syn keyword cType GaussTetra syn keyword cType GaussTria -syn keyword cType GaussianVariogram syn keyword cType GenericExternalResult syn keyword cType GenericOption @@ -1408,4 +1410,5 @@ syn keyword cType IssmDirectApplicInterface syn keyword cType IssmParallelDirectApplicInterface +syn keyword cType krigingobjects syn keyword cType Load syn keyword cType Loads @@ -1418,4 +1421,5 @@ syn keyword cType Matice syn keyword cType Matlitho +syn keyword cType matrixobjects syn keyword cType MatrixParam syn keyword cType Misfit @@ -1430,6 +1434,6 @@ syn keyword cType Observations syn keyword cType Option +syn keyword cType Options syn keyword cType OptionUtilities -syn keyword cType Options syn keyword cType Param syn keyword cType Parameters @@ -1445,11 +1449,11 @@ syn keyword cType Regionaloutput syn keyword cType Results +syn keyword cType Riftfront syn keyword cType RiftStruct -syn keyword cType Riftfront syn keyword cType SealevelMasks syn keyword cType Seg syn keyword cType SegInput2 +syn keyword cType Segment syn keyword cType SegRef -syn keyword cType Segment syn keyword cType SpcDynamic syn keyword cType SpcStatic @@ -1470,8 +1474,4 @@ syn keyword cType Vertex syn keyword cType Vertices -syn keyword cType classes -syn keyword cType gaussobjects -syn keyword cType krigingobjects -syn keyword cType matrixobjects syn keyword cType AdjointBalancethickness2Analysis syn keyword cType AdjointBalancethicknessAnalysis @@ -1492,6 +1492,6 @@ syn keyword cType FreeSurfaceBaseAnalysis syn keyword cType FreeSurfaceTopAnalysis +syn keyword cType GiaAnalysis syn keyword cType GLheightadvectionAnalysis -syn keyword cType GiaAnalysis syn keyword cType HydrologyDCEfficientAnalysis syn keyword cType HydrologyDCInefficientAnalysis Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 25307) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 25308) @@ -176,9 +176,9 @@ HydrologyStorageEnum, HydrologydcEplColapseThicknessEnum, - HydrologydcEplCompressibilityEnum, HydrologydcEplConductivityEnum, HydrologydcEplInitialThicknessEnum, + HydrologydcEplLayerCompressibilityEnum, HydrologydcEplMaxThicknessEnum, - HydrologydcEplPorosityEnum, + HydrologydcEplPoreWaterMassEnum, HydrologydcEplThickCompEnum, HydrologydcEplflipLockEnum, @@ -189,12 +189,12 @@ HydrologydcPenaltyLockEnum, HydrologydcRelTolEnum, - HydrologydcSedimentCompressibilityEnum, HydrologydcSedimentlimitEnum, HydrologydcSedimentlimitFlagEnum, + HydrologydcSedimentLayerCompressibilityEnum, + HydrologydcSedimentPoreWaterMassEnum, HydrologydcSedimentPorosityEnum, HydrologydcSedimentThicknessEnum, HydrologydcTransferFlagEnum, HydrologydcUnconfinedFlagEnum, - HydrologydcWaterCompressibilityEnum, HydrologyshreveStabilizationEnum, IcecapToEarthCommEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 25307) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 25308) @@ -184,9 +184,9 @@ case HydrologyStorageEnum : return "HydrologyStorage"; case HydrologydcEplColapseThicknessEnum : return "HydrologydcEplColapseThickness"; - case HydrologydcEplCompressibilityEnum : return "HydrologydcEplCompressibility"; case HydrologydcEplConductivityEnum : return "HydrologydcEplConductivity"; case HydrologydcEplInitialThicknessEnum : return "HydrologydcEplInitialThickness"; + case HydrologydcEplLayerCompressibilityEnum : return "HydrologydcEplLayerCompressibility"; case HydrologydcEplMaxThicknessEnum : return "HydrologydcEplMaxThickness"; - case HydrologydcEplPorosityEnum : return "HydrologydcEplPorosity"; + case HydrologydcEplPoreWaterMassEnum : return "HydrologydcEplPoreWaterMass"; case HydrologydcEplThickCompEnum : return "HydrologydcEplThickComp"; case HydrologydcEplflipLockEnum : return "HydrologydcEplflipLock"; @@ -197,12 +197,12 @@ case HydrologydcPenaltyLockEnum : return "HydrologydcPenaltyLock"; case HydrologydcRelTolEnum : return "HydrologydcRelTol"; - case HydrologydcSedimentCompressibilityEnum : return "HydrologydcSedimentCompressibility"; case HydrologydcSedimentlimitEnum : return "HydrologydcSedimentlimit"; case HydrologydcSedimentlimitFlagEnum : return "HydrologydcSedimentlimitFlag"; + case HydrologydcSedimentLayerCompressibilityEnum : return "HydrologydcSedimentLayerCompressibility"; + case HydrologydcSedimentPoreWaterMassEnum : return "HydrologydcSedimentPoreWaterMass"; case HydrologydcSedimentPorosityEnum : return "HydrologydcSedimentPorosity"; case HydrologydcSedimentThicknessEnum : return "HydrologydcSedimentThickness"; case HydrologydcTransferFlagEnum : return "HydrologydcTransferFlag"; case HydrologydcUnconfinedFlagEnum : return "HydrologydcUnconfinedFlag"; - case HydrologydcWaterCompressibilityEnum : return "HydrologydcWaterCompressibility"; case HydrologyshreveStabilizationEnum : return "HydrologyshreveStabilization"; case IcecapToEarthCommEnum : return "IcecapToEarthComm"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 25307) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 25308) @@ -187,9 +187,9 @@ else if (strcmp(name,"HydrologyStorage")==0) return HydrologyStorageEnum; else if (strcmp(name,"HydrologydcEplColapseThickness")==0) return HydrologydcEplColapseThicknessEnum; - else if (strcmp(name,"HydrologydcEplCompressibility")==0) return HydrologydcEplCompressibilityEnum; else if (strcmp(name,"HydrologydcEplConductivity")==0) return HydrologydcEplConductivityEnum; else if (strcmp(name,"HydrologydcEplInitialThickness")==0) return HydrologydcEplInitialThicknessEnum; + else if (strcmp(name,"HydrologydcEplLayerCompressibility")==0) return HydrologydcEplLayerCompressibilityEnum; else if (strcmp(name,"HydrologydcEplMaxThickness")==0) return HydrologydcEplMaxThicknessEnum; - else if (strcmp(name,"HydrologydcEplPorosity")==0) return HydrologydcEplPorosityEnum; + else if (strcmp(name,"HydrologydcEplPoreWaterMass")==0) return HydrologydcEplPoreWaterMassEnum; else if (strcmp(name,"HydrologydcEplThickComp")==0) return HydrologydcEplThickCompEnum; else if (strcmp(name,"HydrologydcEplflipLock")==0) return HydrologydcEplflipLockEnum; @@ -200,12 +200,12 @@ else if (strcmp(name,"HydrologydcPenaltyLock")==0) return HydrologydcPenaltyLockEnum; else if (strcmp(name,"HydrologydcRelTol")==0) return HydrologydcRelTolEnum; - else if (strcmp(name,"HydrologydcSedimentCompressibility")==0) return HydrologydcSedimentCompressibilityEnum; else if (strcmp(name,"HydrologydcSedimentlimit")==0) return HydrologydcSedimentlimitEnum; else if (strcmp(name,"HydrologydcSedimentlimitFlag")==0) return HydrologydcSedimentlimitFlagEnum; + else if (strcmp(name,"HydrologydcSedimentLayerCompressibility")==0) return HydrologydcSedimentLayerCompressibilityEnum; + else if (strcmp(name,"HydrologydcSedimentPoreWaterMass")==0) return HydrologydcSedimentPoreWaterMassEnum; else if (strcmp(name,"HydrologydcSedimentPorosity")==0) return HydrologydcSedimentPorosityEnum; else if (strcmp(name,"HydrologydcSedimentThickness")==0) return HydrologydcSedimentThicknessEnum; else if (strcmp(name,"HydrologydcTransferFlag")==0) return HydrologydcTransferFlagEnum; else if (strcmp(name,"HydrologydcUnconfinedFlag")==0) return HydrologydcUnconfinedFlagEnum; - else if (strcmp(name,"HydrologydcWaterCompressibility")==0) return HydrologydcWaterCompressibilityEnum; else if (strcmp(name,"HydrologyshreveStabilization")==0) return HydrologyshreveStabilizationEnum; else if (strcmp(name,"IcecapToEarthComm")==0) return IcecapToEarthCommEnum;