Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h (revision 25246) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h (revision 25247) @@ -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 25246) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 25247) @@ -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 25246) +++ /issm/trunk-jpl/src/c/classes/Inputs2/TransientInput2.cpp (revision 25247) @@ -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_stepAddObject(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,38 @@ iomodel->FindConstant(&hydrology_model,"md.hydrology.model"); if(hydrology_model==HydrologydcEnum){ + IssmDouble sedcomp, sedporo, watcomp, rhofresh, g; + /*FIXME: this cshould go to Analysis!!!*/ - parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.sediment_compressibility",HydrologydcSedimentCompressibilityEnum)); + /* parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.sediment_compressibility",HydrologydcSedimentCompressibilityEnum)); */ + /* parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.water_compressibility",HydrologydcWaterCompressibilityEnum)); */ 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_compressibility",HydrologydcEplCompressibilityEnum)); */ + /* parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.epl_porosity",HydrologydcEplPorosityEnum)); */ 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 +490,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 +536,5 @@ xDelete(matrix); } - xDelete(mdims_array); + xDelete(mdims_array); xDelete(ndims_array); xDelete(array); Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 25246) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 25247) @@ -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,