Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 23365) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp (revision 23366) @@ -313,12 +313,14 @@ } /*Intermediaries */ + int smb_model; IssmDouble dt,scalar,water_head; - IssmDouble water_load,transfer; + IssmDouble water_load,transfer,runoff_value; IssmDouble epl_storing,epl_transmitivity; IssmDouble Jdet; IssmDouble residual,connectivity; - IssmDouble *xyz_list = NULL; - Input* old_wh_input = NULL; + IssmDouble *xyz_list = NULL; + Input* old_wh_input = NULL; + Input* surface_runoff_input = NULL; /*Fetch number of nodes and dof for this finite element*/ @@ -334,15 +336,20 @@ //basalelement->FindParam(&dt,TimesteppingTimeStepEnum); basalelement ->FindParam(&dt,HydrologydtEnum); - - Input* epl_thick_input = basalelement->GetInput(HydrologydcEplThicknessHydrostepEnum); _assert_(epl_thick_input); - Input* sed_head_input = basalelement->GetInput(SedimentHeadHydrostepEnum); _assert_(sed_head_input); - Input* epl_head_input = basalelement->GetInput(EplHeadHydrostepEnum); _assert_(epl_head_input); - Input* water_input = basalelement->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(water_input); - Input* residual_input = basalelement->GetInput(SedimentHeadResidualEnum); _assert_(residual_input); - Input* base_input = basalelement->GetInput(BaseEnum); _assert_(base_input); + basalelement ->FindParam(&smb_model,SmbEnum); + + Input* epl_thick_input = basalelement->GetInput(HydrologydcEplThicknessHydrostepEnum); _assert_(epl_thick_input); + Input* sed_head_input = basalelement->GetInput(SedimentHeadHydrostepEnum); _assert_(sed_head_input); + Input* epl_head_input = basalelement->GetInput(EplHeadHydrostepEnum); _assert_(epl_head_input); + Input* basal_melt_input = basalelement->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(basal_melt_input); + Input* residual_input = basalelement->GetInput(SedimentHeadResidualEnum); _assert_(residual_input); + Input* base_input = basalelement->GetInput(BaseEnum); _assert_(base_input); if(dt!= 0.){ old_wh_input = basalelement->GetInput(EplHeadOldEnum); _assert_(old_wh_input); } + if(smb_model==SMBgradientscomponentsEnum){ + surface_runoff_input = basalelement->GetInput(SmbRunoffEnum); + } + /* Start looping on the number of gaussian points: */ Gauss* gauss = basalelement->NewGauss(2); @@ -355,6 +362,7 @@ /*Loading term*/ - water_input->GetInputValue(&water_load,gauss); - scalar = Jdet*gauss->weight*(water_load); + basal_melt_input->GetInputValue(&water_load,gauss); + surface_runoff_input->GetInputValue(&runoff_value,gauss); + scalar = Jdet*gauss->weight*(water_load+runoff_value); //scalar = Jdet*gauss->weight*(water_load)/epl_transmitivity; if(dt!=0.) scalar = scalar*dt; @@ -524,4 +532,5 @@ IssmDouble EPLgrad2; IssmDouble EPL_N; + IssmDouble opening,closing; femmodel->parameters->FindParam(&domaintype,DomainTypeEnum); @@ -589,7 +598,20 @@ EPLgrad2 = (epl_slopeX[i]*epl_slopeX[i])+(epl_slopeY[i]*epl_slopeY[i]); /*And proceed to the real thing*/ - thickness[i] = old_thickness[i]/(1.0 - -((rho_water*gravity*epl_conductivity*EPLgrad2*dt)/(rho_ice*latentheat)) - +((2.0*A*dt*pow(EPL_N,n))/(pow(n,n)))); + opening=(rho_water*gravity*epl_conductivity*EPLgrad2*dt)/(rho_ice*latentheat); + closing=(2.0*A*dt*pow(EPL_N,n))/(pow(n,n)); + /*implicit*/ + thickness[i] = old_thickness[i]/(1.0-opening+closing); + /*explicit*/ + //thickness[i] = old_thickness[i]*(1.0+opening-closing); + /*centered*/ + //thickness[i] = old_thickness[i]*(1.0+opening-closing)/(1.0-opening+closing); + + if(element->nodes[i]->Sid()==2299){ + printf("for node %i \n",element->nodes[i]->Sid()); + printf(" old thickness is %g \n",old_thickness[i]); + printf(" new thickness is %g \n",thickness[i]); + printf(" closing is %g \n",((2.0*A*dt*pow(EPL_N,n))/(pow(n,n)))); + printf(" opening is %g \n",((rho_water*gravity*epl_conductivity*EPLgrad2*dt)/(rho_ice*latentheat))); + } /*Take care of otherthikening*/ if(thickness[i]>max_thick){ @@ -689,5 +711,9 @@ vec_mask->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL); /* If epl thickness gets under colapse thickness, close the layer */ + if(element->nodes[i]->Sid()==2299){ + printf("thickness for maksing for node %i is %g \n",element->nodes[i]->Sid(),epl_thickness[i]);} if(epl_thickness[i]nodes[i]->Sid()==2299){ + printf("maksing node %i\n",element->nodes[i]->Sid());} vec_mask->SetValue(basalelement->nodes[i]->Sid(),0.,INS_VAL); recurence->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL); Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 23365) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 23366) @@ -105,5 +105,4 @@ iomodel->FetchDataToInput(elements,"md.mesh.vertexonsurface",MeshVertexonsurfaceEnum); } - if(isefficientlayer){ iomodel->FetchDataToInput(elements,"md.hydrology.mask_eplactive_node",HydrologydcMaskEplactiveNodeEnum); @@ -342,12 +341,14 @@ /*Intermediaries */ bool active_element,isefficientlayer; + int smb_model; IssmDouble dt,scalar,sediment_storing; IssmDouble water_head,sediment_transmitivity; - IssmDouble water_load,transfer; + IssmDouble water_load,runoff_value,transfer; IssmDouble Jdet; IssmDouble *xyz_list = NULL; Input* active_element_input = NULL; - Input* old_wh_input = NULL; + Input* old_wh_input = NULL; + Input* surface_runoff_input = NULL; /*Fetch number of nodes and dof for this finite element*/ @@ -361,16 +362,21 @@ basalelement->GetVerticesCoordinates(&xyz_list); //basalelement->FindParam(&dt,TimesteppingTimeStepEnum); - basalelement ->FindParam(&dt,HydrologydtEnum); + basalelement->FindParam(&dt,HydrologydtEnum); basalelement->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum); - - Input* sed_head_input = basalelement->GetInput(SedimentHeadHydrostepEnum); - Input* epl_head_input = basalelement->GetInput(EplHeadHydrostepEnum); - Input* base_input = basalelement->GetInput(BaseEnum); - Input* water_input = basalelement->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(water_input); - Input* SedTrans_input = basalelement->GetInput(HydrologydcSedimentTransmitivityEnum); _assert_(SedTrans_input); + basalelement->FindParam(&smb_model,SmbEnum); + + Input* sed_head_input = basalelement->GetInput(SedimentHeadHydrostepEnum); + Input* epl_head_input = basalelement->GetInput(EplHeadHydrostepEnum); + Input* base_input = basalelement->GetInput(BaseEnum); + Input* basal_melt_input = basalelement->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(basal_melt_input); + Input* SedTrans_input = basalelement->GetInput(HydrologydcSedimentTransmitivityEnum); _assert_(SedTrans_input); if(dt!= 0.){ old_wh_input = basalelement->GetInput(SedimentHeadOldEnum); _assert_(old_wh_input); } + if(smb_model==SMBgradientscomponentsEnum){ + surface_runoff_input = basalelement->GetInput(SmbRunoffEnum); + } + /*Transfer related Inputs*/ if(isefficientlayer){ @@ -388,6 +394,7 @@ /*Loading term*/ if(!isefficientlayer){ - water_input->GetInputValue(&water_load,gauss); - scalar = Jdet*gauss->weight*(water_load); + basal_melt_input->GetInputValue(&water_load,gauss); + surface_runoff_input->GetInputValue(&runoff_value,gauss); + scalar = Jdet*gauss->weight*(water_load+runoff_value); //scalar = Jdet*gauss->weight*(water_load)/sediment_transmitivity; if(dt!=0.) scalar = scalar*dt; @@ -400,6 +407,7 @@ active_element_input->GetInputValue(&active_element); if(!active_element){ - water_input->GetInputValue(&water_load,gauss); - scalar = Jdet*gauss->weight*(water_load); + basal_melt_input->GetInputValue(&water_load,gauss); + surface_runoff_input->GetInputValue(&runoff_value,gauss); + scalar = Jdet*gauss->weight*(water_load+runoff_value); //scalar = Jdet*gauss->weight*(water_load)/sediment_transmitivity; if(dt!=0.) scalar = scalar*dt; Index: /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp (revision 23365) +++ /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp (revision 23366) @@ -38,5 +38,4 @@ /*Figure out smb model: */ iomodel->FindConstant(&smb_model,"md.smb.model"); - switch(smb_model){ case SMBforcingEnum: @@ -152,4 +151,10 @@ iomodel->FetchDataToInput(elements,"md.smb.refreeze",SmbRefreezeEnum,0.); break; + case SMBgradientscomponentsEnum: + iomodel->FetchDataToInput(elements,"md.smb.accualti",SmbAccualtiEnum); + iomodel->FetchDataToInput(elements,"md.smb.accugrad",SmbAccugradEnum); + iomodel->FetchDataToInput(elements,"md.smb.runoffalti",SmbRunoffaltiEnum); + iomodel->FetchDataToInput(elements,"md.smb.runoffgrad",SmbRunoffgradEnum); + break; default: _error_("Surface mass balance model "<AddObject(iomodel->CopyConstantObject("md.smb.accualti",SmbAccualtiEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.smb.accugrad",SmbAccugradEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.smb.runoffalti",SmbRunoffaltiEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.smb.runoffgrad",SmbRunoffgradEnum)); + + iomodel->FetchData(&temp,&N,&M,"md.smb.accuref"); _assert_(N==2); + parameters->AddObject(new TransientParam(SmbAccurefEnum,&temp[0],&temp[M],interp,M)); + iomodel->DeleteData(temp,"md.smb.accuref"); + + iomodel->FetchData(&temp,&N,&M,"md.smb.runoffref"); _assert_(N==2); + parameters->AddObject(new TransientParam(SmbRunoffrefEnum,&temp[0],&temp[M],interp,M)); + iomodel->DeleteData(temp,"md.smb.runoffref"); + break; default: _error_("Surface mass balance model "<(tmp); +} +/*}}}*/ +void Element::SmbGradCompParameterization(void){/*{{{*/ + + /*Are we on the base? If not, return*/ + if(!IsOnBase()) return; + int numvertices = this->GetNumberOfVertices(); + + int i; + IssmDouble accuref, runoffref; //reference values at given altitude + IssmDouble accualti, runoffalti; //reference altitudes + IssmDouble accugrad, runoffgrad; //gradients from reference altitude + IssmDouble rho_water, rho_ice; + IssmDouble time; + + IssmDouble* smb = xNew(numvertices); + IssmDouble* surf = xNew(numvertices); + IssmDouble* accu = xNew(numvertices); + IssmDouble* runoff = xNew(numvertices); + + /*Get material parameters :*/ + rho_water=this->matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum); + rho_ice=this->matpar->GetMaterialParameter(MaterialsRhoIceEnum); + + /*Recover parameters*/ + parameters->FindParam(&time,TimeEnum); + parameters->FindParam(&accualti,SmbAccualtiEnum); + parameters->FindParam(&accugrad,SmbAccugradEnum); + parameters->FindParam(&runoffalti,SmbRunoffaltiEnum); + parameters->FindParam(&runoffgrad,SmbRunoffgradEnum); + + /*Recover reference values at current time*/ + parameters->FindParam(&accuref,SmbAccurefEnum,time); + parameters->FindParam(&runoffref,SmbRunoffrefEnum,time); + + /*Recover surface elevation*/ + GetInputListOnVertices(&surf[0],SurfaceEnum); + + /*Compute the temperature and precipitation*/ + for(int iv=0;ivObjectEnum()){ + case TriaEnum: + this->inputs->AddInput(new TriaInput(SmbMassBalanceEnum,&smb[0],P1Enum)); + this->inputs->AddInput(new TriaInput(SmbRunoffEnum,&runoff[0],P1Enum)); + break; + case PentaEnum: + this->inputs->AddInput(new PentaInput(SmbMassBalanceEnum,&smb[0],P1Enum)); + this->inputs->AddInput(new PentaInput(SmbRunoffEnum,&runoff[0],P1Enum)); + this->InputExtrude(SmbMassBalanceEnum,-1); + this->InputExtrude(SmbRunoffEnum,-1); + break; + case TetraEnum: + this->inputs->AddInput(new TetraInput(SmbMassBalanceEnum,&smb[0],P1Enum)); + this->inputs->AddInput(new TetraInput(SmbRunoffEnum,&runoff[0],P1Enum)); + this->InputExtrude(SmbMassBalanceEnum,-1); + this->InputExtrude(SmbRunoffEnum,-1); + break; + default: _error_("Not implemented yet"); + } + /* this->AddInput(SmbMassBalanceEnum,smb,P1Enum); */ + /* this->AddInput(SmbRunoffEnum,runoff,P1Enum); */ + /*clean-up*/ + xDelete(surf); + xDelete(accu); + xDelete(runoff); + xDelete(smb); } /*}}}*/ @@ -2588,9 +2658,9 @@ rho_water=this->matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum); rho_ice=this->matpar->GetMaterialParameter(MaterialsRhoIceEnum); - + /*Get parameters for height corrections*/ desfac=this->matpar->GetMaterialParameter(SmbDesfacEnum); rlaps=this->matpar->GetMaterialParameter(SmbRlapsEnum); - + /*Recover monthly temperatures and precipitation*/ Input* input=this->inputs->GetInput(SmbMonthlytemperaturesEnum); _assert_(input); @@ -2598,5 +2668,5 @@ /*Recover smb correction term */ Input* input3=this->inputs->GetInput(SmbSmbCorrEnum); _assert_(input3); - + /* Get time */ this->parameters->FindParam(&time,TimeEnum); @@ -2607,5 +2677,5 @@ IssmDouble MU_0 = 9.7155; //Firn-warming correction, in (d*deg C)/(mm WE) IssmDouble mu = MU_0*(1000.0*86400.0)*(rho_ice/rho_water); // (d*deg C)/(mm WE) --> (s*deg C)/(m IE) - + /*loop over vertices: */ Gauss* gauss=this->NewGauss(); @@ -2638,6 +2708,6 @@ /*Get yearlytemperatures */ for(int month=0;month<12;month++) yearlytemperatures[iv]=yearlytemperatures[iv]+((monthlytemperatures[iv*12+month]+273.15)+t_ampl[iv])*inv_twelve; // Has to be in Kelvin - - if(isfirnwarming){ + + if(isfirnwarming){ if(melt_star[iv]>=melt[iv]){ yearlytemperatures[iv]= yearlytemperatures[iv]+mu*(melt_star[iv]-melt[iv]); @@ -2720,6 +2790,6 @@ xDelete(s0t); xDelete(s0p); - xDelete(t_ampl); - xDelete(p_ampl); + xDelete(t_ampl); + xDelete(p_ampl); xDelete(smbcorr); xDelete(melt_star); @@ -2733,6 +2803,6 @@ /*Some intputs need to be computed, even if they are already in inputs, they might not be up to date!*/ switch(output_enum){ - case ViscousHeatingEnum: this->ViscousHeatingCreateInput(); break; - case StressMaxPrincipalEnum: this->StressMaxPrincipalCreateInput(); break; + case ViscousHeatingEnum: this->ViscousHeatingCreateInput(); break; + case StressMaxPrincipalEnum: this->StressMaxPrincipalCreateInput(); break; case StressTensorxxEnum: case StressTensorxyEnum: @@ -2747,5 +2817,5 @@ case StrainRateyzEnum: case StrainRatezzEnum: - case StrainRateeffectiveEnum: this->ComputeStrainRate(); break; + case StrainRateeffectiveEnum: this->ComputeStrainRate(); break; case DeviatoricStressxxEnum: case DeviatoricStressxyEnum: Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 23365) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 23366) @@ -1,7 +1,7 @@ /*!\file: Element.h * \brief abstract class for Element object - * This class is a place holder for the Tria and the Penta elements. + * This class is a place holder for the Tria and the Penta elements. * It is derived from Element, so DataSets can contain them. - */ + */ #ifndef _ELEMENT_H_ @@ -50,5 +50,5 @@ int element_type; - public: + public: /*Constructors/Destructores*/ Element(); @@ -68,4 +68,5 @@ void Delta18oParameterization(void); void Delta18opdParameterization(void); + void SmbGradCompParameterization(void); IssmDouble Divergence(void); void dViscositydBFS(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input); @@ -130,14 +131,14 @@ void InputUpdateFromConstant(int constant, int name); void InputUpdateFromConstant(bool constant, int name); - bool IsFloating(); + bool IsFloating(); bool IsIceInElement(); bool IsLandInElement(); bool IsWaterInElement(); - void LinearFloatingiceMeltingRate(); - void SpatialLinearFloatingiceMeltingRate(); - void MantlePlumeGeothermalFlux(); + void LinearFloatingiceMeltingRate(); + void SpatialLinearFloatingiceMeltingRate(); + void MantlePlumeGeothermalFlux(); void MarshallElement(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction,int numanalyses); void MigrateGroundingLine(IssmDouble* sheet_ungrounding); - void MismipFloatingiceMeltingRate(); + void MismipFloatingiceMeltingRate(); void MungsmtpParameterization(void); ElementMatrix* NewElementMatrix(int approximation_enum=NoneApproximationEnum); @@ -241,5 +242,5 @@ virtual bool IsFaceOnBoundary(void)=0; virtual bool IsIcefront(void)=0; - virtual bool IsNodeOnShelfFromFlags(IssmDouble* flags)=0; + virtual bool IsNodeOnShelfFromFlags(IssmDouble* flags)=0; virtual bool IsOnBase()=0; virtual bool IsOnSurface()=0; @@ -331,5 +332,5 @@ virtual IssmDouble OceanAverage(IssmDouble* Sg)=0; virtual IssmDouble OceanArea(void)=0; - virtual void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea)=0; + virtual void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea)=0; virtual void SealevelriseEustatic(Vector* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea)=0; virtual void SealevelriseNonEustatic(Vector* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea)=0; Index: /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp (revision 23365) +++ /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp (revision 23366) @@ -115,5 +115,5 @@ iomodel->FindConstant(&this->rlaps,"md.smb.rlaps"); iomodel->FindConstant(&this->rlapslgm,"md.smb.rlapslgm"); - iomodel->FindConstant(&this->dpermil,"md.smb.dpermil"); + iomodel->FindConstant(&this->dpermil,"md.smb.dpermil"); case SMBgradientsEnum: /*Nothing to add*/ @@ -129,4 +129,7 @@ break; case SMBmeltcomponentsEnum: + /*Nothing to add*/ + break; + case SMBgradientscomponentsEnum: /*Nothing to add*/ break; @@ -174,9 +177,9 @@ break; case MaterialsEnum: - //we have several types of materials. Retrieve this info first: + //we have several types of materials. Retrieve this info first: iomodel->FetchData(&nature,&nnat,&dummy,"md.materials.nature"); - //go through list of materials, and create constant parameters accordingly: - for(int i=0;iFindConstant(&this->rlaps,"md.smb.rlaps"); iomodel->FindConstant(&this->rlapslgm,"md.smb.rlapslgm"); - iomodel->FindConstant(&this->dpermil,"md.smb.dpermil"); + iomodel->FindConstant(&this->dpermil,"md.smb.dpermil"); case SMBgradientsEnum: /*Nothing to add*/ @@ -241,5 +244,5 @@ default: _error_("Surface mass balance model "<FindConstant(&this->sediment_compressibility,"md.hydrology.sediment_compressibility"); @@ -358,5 +361,5 @@ this->Echo(); -} +} /*}}}*/ void Matpar::Echo(void){/*{{{*/ @@ -389,16 +392,16 @@ _printf_(" sediment_porosity: " << sediment_porosity << "\n"); _printf_(" sediment_thickness: " << sediment_thickness << "\n"); - _printf_(" water_compressibility: " << water_compressibility << "\n"); - _printf_(" epl_compressibility: " << epl_compressibility << "\n"); - _printf_(" epl_porosity: " << epl_porosity << "\n"); - _printf_(" epl_init_thickness: " << epl_init_thickness << "\n"); - _printf_(" epl_colapse_thickness: " << epl_colapse_thickness << "\n"); - _printf_(" epl_max_thickness: " << epl_max_thickness << "\n"); - _printf_(" epl_conductivity: " << epl_conductivity << "\n"); - _printf_(" lithosphere_shear_modulus: " << lithosphere_shear_modulus << "\n"); - _printf_(" lithosphere_density: " << lithosphere_density << "\n"); - _printf_(" mantle_shear_modulus: " << mantle_shear_modulus << "\n"); - _printf_(" mantle_density: " << mantle_density << "\n"); - _printf_(" earth_density: " << earth_density << "\n"); + _printf_(" water_compressibility: " << water_compressibility << "\n"); + _printf_(" epl_compressibility: " << epl_compressibility << "\n"); + _printf_(" epl_porosity: " << epl_porosity << "\n"); + _printf_(" epl_init_thickness: " << epl_init_thickness << "\n"); + _printf_(" epl_colapse_thickness: " << epl_colapse_thickness << "\n"); + _printf_(" epl_max_thickness: " << epl_max_thickness << "\n"); + _printf_(" epl_conductivity: " << epl_conductivity << "\n"); + _printf_(" lithosphere_shear_modulus: " << lithosphere_shear_modulus << "\n"); + _printf_(" lithosphere_density: " << lithosphere_density << "\n"); + _printf_(" mantle_shear_modulus: " << mantle_shear_modulus << "\n"); + _printf_(" mantle_density: " << mantle_density << "\n"); + _printf_(" earth_density: " << earth_density << "\n"); return; } @@ -520,5 +523,5 @@ this->dpermil=constant; break; - default: + default: break; } Index: /issm/trunk-jpl/src/c/cores/transient_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/transient_core.cpp (revision 23365) +++ /issm/trunk-jpl/src/c/cores/transient_core.cpp (revision 23366) @@ -180,5 +180,5 @@ for(int i=0;i0.) icebase_oceangrid[i]=+9999.; xDelete(icemask_oceangrid); - + if(my_rank==0){ ISSM_MPI_Send(icebase_oceangrid,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001008,tomitgcmcomm); @@ -384,4 +384,7 @@ } + /*duplicating smb position to have runoff value*/ + if(issmb) smb_core(femmodel); + if(ishydrology) hydrology_core(femmodel); @@ -393,5 +396,4 @@ /* from here on, prepare geometry for next time step*/ - if(issmb) smb_core(femmodel); @@ -403,8 +405,8 @@ if(isgroundingline){ - + /*Start profiler*/ femmodel->profiler->Start(GROUNDINGLINECORE); - + if(VerboseSolution()) _printf0_(" computing new grounding line position\n"); GroundinglineMigrationx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); @@ -416,5 +418,5 @@ femmodel->parameters->SetParam(SurfaceEnum,InputToExtrudeEnum); extrudefrombase_core(femmodel); - + /*Stop profiler*/ femmodel->profiler->Stop(GROUNDINGLINECORE); Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 23365) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 23366) @@ -1,4 +1,4 @@ /*!\file SurfaceMassBalancex - * \brief: calculates SMB + * \brief: calculates SMB */ @@ -105,5 +105,5 @@ for(v=0;vela[v]){ + if(s[v]>ela[v]){ smb[v]=b_pos[v]*(s[v]-ela[v]); } @@ -167,5 +167,5 @@ // INPUT: surface elevation (m): hd(NA) // monthly mean surface sealevel temperature (degrees C): vTempsea(NA - // ,NTIME) + // ,NTIME) // monthly mean precip rate (m/yr water equivalent): vPrec(NA,NTIME) // OUTPUT: mass-balance (m/yr ice): agd(NA) @@ -174,5 +174,5 @@ int i, it, jj, itm; IssmDouble DT = 0.02, sigfac, snormfac; - IssmDouble signorm = 5.5; // signorm : sigma of the temperature distribution for a normal day + IssmDouble signorm = 5.5; // signorm : sigma of the temperature distribution for a normal day IssmDouble siglim; // sigma limit for the integration which is equal to 2.5 sigmanorm IssmDouble signormc = signorm - 0.5; // sigma of the temperature distribution for cloudy day @@ -180,5 +180,5 @@ IssmDouble tstep, tsint, tint, tstepc; int NPDMAX = 1504, NPDCMAX = 1454; - //IssmDouble pdds[NPDMAX]={0}; + //IssmDouble pdds[NPDMAX]={0}; //IssmDouble pds[NPDCMAX]={0}; IssmDouble pddt, pd ; // pd : snow/precip fraction, precipitation falling as snow @@ -193,6 +193,6 @@ Element *element = NULL; - pdds=xNew(NPDMAX+1); - pds=xNew(NPDCMAX+1); + pdds=xNew(NPDMAX+1); + pds=xNew(NPDCMAX+1); // Get ismungsm parameter @@ -216,5 +216,5 @@ if(itm >= NPDMAX) _error_("increase NPDMAX in massBalance.cpp"); - for(it = 0; it < itm; it++){ + for(it = 0; it < itm; it++){ // tstar = REAL(it)*DT-siglim; tstar = it*DT-siglim; @@ -263,8 +263,8 @@ }/*}}}*/ void PositiveDegreeDaySicopolisx(FemModel* femmodel){/*{{{*/ - + bool isfirnwarming; femmodel->parameters->FindParam(&isfirnwarming,SmbIsfirnwarmingEnum); - + for(int i=0;ielements->Size();i++){ Element* element=xDynamicCast(femmodel->elements->GetObjectByOffset(i)); @@ -357,7 +357,7 @@ /*Allocate all arrays*/ int numvertices = element->GetNumberOfVertices(); - IssmDouble* acc = xNew(numvertices); + IssmDouble* acc = xNew(numvertices); IssmDouble* evap = xNew(numvertices); - IssmDouble* runoff = xNew(numvertices); + IssmDouble* runoff = xNew(numvertices); IssmDouble* smb = xNew(numvertices); @@ -399,5 +399,5 @@ int numvertices = element->GetNumberOfVertices(); IssmDouble* acc = xNew(numvertices); - IssmDouble* evap = xNew(numvertices); + IssmDouble* evap = xNew(numvertices); IssmDouble* melt = xNew(numvertices); IssmDouble* refreeze = xNew(numvertices); @@ -425,2 +425,15 @@ }/*}}}*/ +void SmbGradientsComponentsx(FemModel* femmodel){/*{{{*/ + + // void SurfaceMassBalancex(hd,agd,ni){ + // INPUT parameters: ni: working size of arrays + // INPUT: surface elevation (m): hd(NA) + // OUTPUT: mass-balance (m/yr ice): agd(NA) + + for(int i=0;ielements->Size();i++){ + Element* element=xDynamicCast(femmodel->elements->GetObjectByOffset(i)); + element->SmbGradCompParameterization(); + } + +}/*}}}*/ Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h (revision 23365) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h (revision 23366) @@ -1,5 +1,5 @@ /*!\file: SurfaceMassBalancex.h * \brief header file for SMB - */ + */ #ifndef _SurfaceMassBalancex_H @@ -19,9 +19,9 @@ void SmbHenningx(FemModel* femmodel); void SmbComponentsx(FemModel* femmodel); -void SmbMeltComponentsx(FemModel* femmodel); - +void SmbMeltComponentsx(FemModel* femmodel); +void SmbGradientsComponentsx(FemModel* femmodel); /*GEMB: */ void Gembx(FemModel* femmodel); -void GembgridInitialize(IssmDouble** pdz, int* psize, IssmDouble zTop, IssmDouble dzTop, IssmDouble zMax, IssmDouble zY); +void GembgridInitialize(IssmDouble** pdz, int* psize, IssmDouble zTop, IssmDouble dzTop, IssmDouble zMax, IssmDouble zY); IssmDouble Marbouty(IssmDouble T, IssmDouble d, IssmDouble dT); void grainGrowth(IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, IssmDouble* T,IssmDouble* dz,IssmDouble* d, IssmDouble* W,IssmDouble smb_dt,int m,int aIdx, int sid); @@ -29,5 +29,5 @@ void shortwave(IssmDouble** pswf, int swIdx, int aIdx, IssmDouble dsw, IssmDouble as, IssmDouble* d, IssmDouble* dz, IssmDouble* re, IssmDouble dIce, int m, int sid); void thermo(IssmDouble* pEC, IssmDouble** T, IssmDouble* dz, IssmDouble* d, IssmDouble* swf, IssmDouble dlw, IssmDouble Ta, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble teValue, IssmDouble Ws, IssmDouble dt0, int m, IssmDouble Vz, IssmDouble Tz, IssmDouble thermo_scaling, IssmDouble dIce, int sid); -void accumulation(IssmDouble** pT, IssmDouble** pdz, IssmDouble** pd, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pm, int aIdx,IssmDouble Ta, IssmDouble P, IssmDouble dzMin, IssmDouble aSnow, IssmDouble dIce, int sid); +void accumulation(IssmDouble** pT, IssmDouble** pdz, IssmDouble** pd, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pm, int aIdx,IssmDouble Ta, IssmDouble P, IssmDouble dzMin, IssmDouble aSnow, IssmDouble dIce, int sid); void melt(IssmDouble* pM, IssmDouble* pR, IssmDouble* pmAdd, IssmDouble* pdz_add, IssmDouble** pT, IssmDouble** pd, IssmDouble** pdz, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pn, IssmDouble dzMin, IssmDouble zMax, IssmDouble zMin, IssmDouble zTop, IssmDouble dIce, int sid); void densification(IssmDouble** pd,IssmDouble** pdz, IssmDouble* T, IssmDouble* re, int denIdx, IssmDouble C, IssmDouble dt, IssmDouble Tmean, IssmDouble dIce, int m, int sid); Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 23365) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 23366) @@ -273,4 +273,7 @@ SettingsSolverResidueThresholdEnum, SettingsWaitonlockEnum, + SmbAccualtiEnum, + SmbAccugradEnum, + SmbAccurefEnum, SmbAdThreshEnum, SmbAIceEnum, @@ -304,4 +307,7 @@ SmbPfacEnum, SmbRequestedOutputsEnum, + SmbRunoffaltiEnum, + SmbRunoffgradEnum, + SmbRunoffrefEnum, SmbSealevEnum, SmbSwIdxEnum, @@ -528,7 +534,7 @@ SealevelNEsaRateEnum, SealevelRSLRateEnum, - SealevelUEsaEnum, + SealevelUEsaEnum, SealevelRSLEustaticRateEnum, - SealevelNEsaEnum, + SealevelNEsaEnum, SealevelUGiaEnum, SealevelNGiaEnum, @@ -1101,4 +1107,5 @@ SMBpddEnum, SMBpddSicopolisEnum, + SMBgradientscomponentsEnum, SmbRlapsEnum, SmbRlapslgmEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 23365) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 23366) @@ -281,4 +281,7 @@ case SettingsSolverResidueThresholdEnum : return "SettingsSolverResidueThreshold"; case SettingsWaitonlockEnum : return "SettingsWaitonlock"; + case SmbAccualtiEnum : return "SmbAccualti"; + case SmbAccugradEnum : return "SmbAccugrad"; + case SmbAccurefEnum : return "SmbAccuref"; case SmbAdThreshEnum : return "SmbAdThresh"; case SmbAIceEnum : return "SmbAIce"; @@ -312,4 +315,7 @@ case SmbPfacEnum : return "SmbPfac"; case SmbRequestedOutputsEnum : return "SmbRequestedOutputs"; + case SmbRunoffaltiEnum : return "SmbRunoffalti"; + case SmbRunoffgradEnum : return "SmbRunoffgrad"; + case SmbRunoffrefEnum : return "SmbRunoffref"; case SmbSealevEnum : return "SmbSealev"; case SmbSwIdxEnum : return "SmbSwIdx"; @@ -1105,4 +1111,5 @@ case SMBpddEnum : return "SMBpdd"; case SMBpddSicopolisEnum : return "SMBpddSicopolis"; + case SMBgradientscomponentsEnum : return "SMBgradientscomponents"; case SmbRlapsEnum : return "SmbRlaps"; case SmbRlapslgmEnum : return "SmbRlapslgm"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 23365) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 23366) @@ -287,4 +287,7 @@ else if (strcmp(name,"SettingsSolverResidueThreshold")==0) return SettingsSolverResidueThresholdEnum; else if (strcmp(name,"SettingsWaitonlock")==0) return SettingsWaitonlockEnum; + else if (strcmp(name,"SmbAccualti")==0) return SmbAccualtiEnum; + else if (strcmp(name,"SmbAccugrad")==0) return SmbAccugradEnum; + else if (strcmp(name,"SmbAccuref")==0) return SmbAccurefEnum; else if (strcmp(name,"SmbAdThresh")==0) return SmbAdThreshEnum; else if (strcmp(name,"SmbAIce")==0) return SmbAIceEnum; @@ -318,4 +321,7 @@ else if (strcmp(name,"SmbPfac")==0) return SmbPfacEnum; else if (strcmp(name,"SmbRequestedOutputs")==0) return SmbRequestedOutputsEnum; + 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,"SmbSwIdx")==0) return SmbSwIdxEnum; @@ -377,5 +383,8 @@ else if (strcmp(name,"TransientIsmovingfront")==0) return TransientIsmovingfrontEnum; else if (strcmp(name,"TransientIsoceancoupling")==0) return TransientIsoceancouplingEnum; - else if (strcmp(name,"TransientIsslr")==0) return TransientIsslrEnum; + else stage=4; + } + if(stage==4){ + if (strcmp(name,"TransientIsslr")==0) return TransientIsslrEnum; else if (strcmp(name,"TransientIssmb")==0) return TransientIssmbEnum; else if (strcmp(name,"TransientIsstressbalance")==0) return TransientIsstressbalanceEnum; @@ -383,8 +392,5 @@ else if (strcmp(name,"TransientNumRequestedOutputs")==0) return TransientNumRequestedOutputsEnum; else if (strcmp(name,"TransientRequestedOutputs")==0) return TransientRequestedOutputsEnum; - else stage=4; - } - if(stage==4){ - if (strcmp(name,"Velocity")==0) return VelocityEnum; + else if (strcmp(name,"Velocity")==0) return VelocityEnum; else if (strcmp(name,"WorldComm")==0) return WorldCommEnum; else if (strcmp(name,"ParametersEND")==0) return ParametersENDEnum; @@ -500,5 +506,8 @@ else if (strcmp(name,"Ice")==0) return IceEnum; else if (strcmp(name,"IceMaskNodeActivation")==0) return IceMaskNodeActivationEnum; - else if (strcmp(name,"Input")==0) return InputEnum; + else stage=5; + } + if(stage==5){ + if (strcmp(name,"Input")==0) return InputEnum; else if (strcmp(name,"InversionCostFunctionsCoefficients")==0) return InversionCostFunctionsCoefficientsEnum; else if (strcmp(name,"InversionSurfaceObs")==0) return InversionSurfaceObsEnum; @@ -506,8 +515,5 @@ else if (strcmp(name,"InversionVelObs")==0) return InversionVelObsEnum; else if (strcmp(name,"InversionVxObs")==0) return InversionVxObsEnum; - else stage=5; - } - if(stage==5){ - if (strcmp(name,"InversionVyObs")==0) return InversionVyObsEnum; + else if (strcmp(name,"InversionVyObs")==0) return InversionVyObsEnum; else if (strcmp(name,"LevelsetfunctionSlopeX")==0) return LevelsetfunctionSlopeXEnum; else if (strcmp(name,"LevelsetfunctionSlopeY")==0) return LevelsetfunctionSlopeYEnum; @@ -623,5 +629,8 @@ else if (strcmp(name,"SmbVz")==0) return SmbVzEnum; else if (strcmp(name,"SmbW")==0) return SmbWEnum; - else if (strcmp(name,"SmbWini")==0) return SmbWiniEnum; + else stage=6; + } + if(stage==6){ + if (strcmp(name,"SmbWini")==0) return SmbWiniEnum; else if (strcmp(name,"SmbZMax")==0) return SmbZMaxEnum; else if (strcmp(name,"SmbZMin")==0) return SmbZMinEnum; @@ -629,8 +638,5 @@ else if (strcmp(name,"SmbZY")==0) return SmbZYEnum; else if (strcmp(name,"StrainRateeffective")==0) return StrainRateeffectiveEnum; - else stage=6; - } - if(stage==6){ - if (strcmp(name,"StrainRateparallel")==0) return StrainRateparallelEnum; + else if (strcmp(name,"StrainRateparallel")==0) return StrainRateparallelEnum; else if (strcmp(name,"StrainRateperpendicular")==0) return StrainRateperpendicularEnum; else if (strcmp(name,"StrainRatexx")==0) return StrainRatexxEnum; @@ -746,5 +752,8 @@ else if (strcmp(name,"DatasetInput")==0) return DatasetInputEnum; else if (strcmp(name,"DataSetParam")==0) return DataSetParamEnum; - else if (strcmp(name,"DefaultAnalysis")==0) return DefaultAnalysisEnum; + else stage=7; + } + if(stage==7){ + if (strcmp(name,"DefaultAnalysis")==0) return DefaultAnalysisEnum; else if (strcmp(name,"DefaultCalving")==0) return DefaultCalvingEnum; else if (strcmp(name,"DegreeOfChannelization")==0) return DegreeOfChannelizationEnum; @@ -752,8 +761,5 @@ else if (strcmp(name,"DependentObject")==0) return DependentObjectEnum; else if (strcmp(name,"DepthAverageAnalysis")==0) return DepthAverageAnalysisEnum; - else stage=7; - } - if(stage==7){ - if (strcmp(name,"DeviatoricStresseffective")==0) return DeviatoricStresseffectiveEnum; + else if (strcmp(name,"DeviatoricStresseffective")==0) return DeviatoricStresseffectiveEnum; else if (strcmp(name,"DeviatoricStressErrorEstimator")==0) return DeviatoricStressErrorEstimatorEnum; else if (strcmp(name,"Divergence")==0) return DivergenceEnum; @@ -869,5 +875,8 @@ else if (strcmp(name,"LambdaS")==0) return LambdaSEnum; else if (strcmp(name,"LATaylorHood")==0) return LATaylorHoodEnum; - else if (strcmp(name,"LevelsetAnalysis")==0) return LevelsetAnalysisEnum; + else stage=8; + } + if(stage==8){ + if (strcmp(name,"LevelsetAnalysis")==0) return LevelsetAnalysisEnum; else if (strcmp(name,"LevelsetfunctionPicard")==0) return LevelsetfunctionPicardEnum; else if (strcmp(name,"LinearFloatingMeltRate")==0) return LinearFloatingMeltRateEnum; @@ -875,8 +884,5 @@ else if (strcmp(name,"LliboutryDuval")==0) return LliboutryDuvalEnum; else if (strcmp(name,"Loads")==0) return LoadsEnum; - else stage=8; - } - if(stage==8){ - if (strcmp(name,"LoveAnalysis")==0) return LoveAnalysisEnum; + else if (strcmp(name,"LoveAnalysis")==0) return LoveAnalysisEnum; else if (strcmp(name,"LoveHi")==0) return LoveHiEnum; else if (strcmp(name,"LoveHr")==0) return LoveHrEnum; @@ -992,5 +998,8 @@ else if (strcmp(name,"Outputdefinition37")==0) return Outputdefinition37Enum; else if (strcmp(name,"Outputdefinition38")==0) return Outputdefinition38Enum; - else if (strcmp(name,"Outputdefinition39")==0) return Outputdefinition39Enum; + else stage=9; + } + if(stage==9){ + if (strcmp(name,"Outputdefinition39")==0) return Outputdefinition39Enum; else if (strcmp(name,"Outputdefinition3")==0) return Outputdefinition3Enum; else if (strcmp(name,"Outputdefinition40")==0) return Outputdefinition40Enum; @@ -998,8 +1007,5 @@ else if (strcmp(name,"Outputdefinition42")==0) return Outputdefinition42Enum; else if (strcmp(name,"Outputdefinition43")==0) return Outputdefinition43Enum; - else stage=9; - } - if(stage==9){ - if (strcmp(name,"Outputdefinition44")==0) return Outputdefinition44Enum; + else if (strcmp(name,"Outputdefinition44")==0) return Outputdefinition44Enum; else if (strcmp(name,"Outputdefinition45")==0) return Outputdefinition45Enum; else if (strcmp(name,"Outputdefinition46")==0) return Outputdefinition46Enum; @@ -1115,5 +1121,8 @@ else if (strcmp(name,"SMBcomponents")==0) return SMBcomponentsEnum; else if (strcmp(name,"SMBd18opdd")==0) return SMBd18opddEnum; - else if (strcmp(name,"SmbDesfac")==0) return SmbDesfacEnum; + else stage=10; + } + if(stage==10){ + if (strcmp(name,"SmbDesfac")==0) return SmbDesfacEnum; else if (strcmp(name,"SmbDpermil")==0) return SmbDpermilEnum; else if (strcmp(name,"SmbDzAdd")==0) return SmbDzAddEnum; @@ -1121,8 +1130,5 @@ else if (strcmp(name,"SMBforcing")==0) return SMBforcingEnum; else if (strcmp(name,"SMBgcm")==0) return SMBgcmEnum; - else stage=10; - } - if(stage==10){ - if (strcmp(name,"SMBgemb")==0) return SMBgembEnum; + else if (strcmp(name,"SMBgemb")==0) return SMBgembEnum; else if (strcmp(name,"SMBgradientsela")==0) return SMBgradientselaEnum; else if (strcmp(name,"SMBgradients")==0) return SMBgradientsEnum; @@ -1132,4 +1138,5 @@ else if (strcmp(name,"SMBpdd")==0) return SMBpddEnum; else if (strcmp(name,"SMBpddSicopolis")==0) return SMBpddSicopolisEnum; + else if (strcmp(name,"SMBgradientscomponents")==0) return SMBgradientscomponentsEnum; else if (strcmp(name,"SmbRlaps")==0) return SmbRlapsEnum; else if (strcmp(name,"SmbRlapslgm")==0) return SmbRlapslgmEnum; Index: /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp (revision 23365) +++ /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp (revision 23366) @@ -174,4 +174,5 @@ case 9: return SMBgradientselaEnum; case 10: return SMBpddSicopolisEnum; + case 11: return SMBgradientscomponentsEnum; default: _error_("Marshalled SMB code \""<