Index: /issm/trunk-jpl/src/c/cores/hydrology_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/hydrology_core.cpp (revision 23664) +++ /issm/trunk-jpl/src/c/cores/hydrology_core.cpp (revision 23665) @@ -19,7 +19,10 @@ int solution_type; int numoutputs = 0; + int smboutputs; bool save_results; bool modify_loads = true; + bool issmb; char **requested_outputs = NULL; + char **requested_smb_outputs = NULL; IssmDouble ThawedNodes; @@ -30,4 +33,6 @@ femmodel->parameters->FindParam(&numoutputs,HydrologyNumRequestedOutputsEnum); if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,HydrologyRequestedOutputsEnum); + femmodel->parameters->FindParam(&issmb,TransientIssmbEnum); + /*Using the Shreve based Model*/ if (hydrology_model==HydrologyshreveEnum){ @@ -40,5 +45,4 @@ femmodel->SetCurrentConfiguration(HydrologyShreveAnalysisEnum); solutionsequence_nonlinear(femmodel,modify_loads); - /*transfer water column thickness to old water column thickness: */ InputDuplicatex(femmodel,WatercolumnEnum,WaterColumnOldEnum); @@ -50,8 +54,11 @@ /*intermediary: */ bool isefficientlayer; - int hydrostep,hydroslices; - IssmDouble time,init_time,hydrotime,yts; + int hydrostep,hydroslices,numaveragedinput; + IssmDouble time,hydrotime,yts; IssmDouble dt,hydrodt; - + /*SMB related */ + int smb_model; + + /*recover parameters: */ femmodel->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum); femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum); @@ -60,54 +67,72 @@ femmodel->parameters->FindParam(&yts,ConstantsYtsEnum); + /*recover SMB related parameters: */ + if(issmb){ + femmodel->parameters->FindParam(&smb_model,SmbEnum); + femmodel->parameters->FindParam(&smboutputs,SmbNumRequestedOutputsEnum); + if(smboutputs) femmodel->parameters->FindParam(&requested_smb_outputs,&smboutputs,SmbRequestedOutputsEnum); + } + /*first we exclude frozen nodes of the solved nodes*/ femmodel->SetCurrentConfiguration(HydrologyDCInefficientAnalysisEnum); femmodel->HydrologyIDSupdateDomainx(&ThawedNodes); + if(ThawedNodes>0){ - init_time=time-dt; //getting the time back to the start of the timestep - hydrotime=init_time; + hydrotime=time-dt; //getting the time back to the start of the timestep hydrodt=dt/hydroslices; //computing hydro dt from dt and a divider hydrostep=0; femmodel->parameters->AddObject(new DoubleParam(HydrologydtEnum,hydrodt)); + if(hydroslices>1){ /*define which variable needs to be averaged on the sub-timestep and initialize as needed*/ + numaveragedinput = 2; + std::vector inputtostack = {EffectivePressureHydrostepEnum,SedimentHeadHydrostepEnum}; + std::vector stackedinput = {EffectivePressureStackedEnum,SedimentHeadStackedEnum}; + std::vector averagedinput = {EffectivePressureEnum,SedimentHeadEnum}; + if (isefficientlayer){ - int inputtostack[4]={EffectivePressureHydrostepEnum,SedimentHeadHydrostepEnum,EplHeadHydrostepEnum,HydrologydcEplThicknessHydrostepEnum}; - int stackedinput[4]={EffectivePressureStackedEnum,SedimentHeadStackedEnum,EplHeadStackedEnum,HydrologydcEplThicknessStackedEnum}; - int averagedinput[4]={EffectivePressureEnum,SedimentHeadEnum,EplHeadEnum,HydrologydcEplThicknessEnum}; - femmodel->InitMeanOutputx(&stackedinput[0],4); - while(hydrotimeInitMeanOutputx(&stackedinput[0],numaveragedinput); + + while(hydrotimeparameters->SetParam(hydrotime,TimeEnum); + if(VerboseSolution()) _printf0_("sub iteration " << hydrostep << "/" << hydroslices << " time [yr]: " << setprecision(4) << hydrotime/yts << " (time step: " << hydrodt/yts << ")\n"); + if(issmb){ + if(VerboseSolution()) _printf0_(" computing mass balance\n"); + SmbAnalysis* analysis = new SmbAnalysis(); + analysis->Core(femmodel); + delete analysis; + } + if(VerboseSolution()) _printf0_(" computing water heads\n"); + /*save preceding timestep*/ + InputDuplicatex(femmodel,SedimentHeadHydrostepEnum,SedimentHeadOldEnum); + if(isefficientlayer){ InputDuplicatex(femmodel,EplHeadHydrostepEnum,EplHeadOldEnum); InputDuplicatex(femmodel,HydrologydcEplThicknessHydrostepEnum,HydrologydcEplThicknessOldEnum); - /*Proceed now to heads computations*/ - solutionsequence_hydro_nonlinear(femmodel); - /*If we have a sub-timestep we stack the variables here*/ - femmodel->SumOutputx(&inputtostack[0],&stackedinput[0],4); } - femmodel->AverageSumOutputx(&stackedinput[0],&averagedinput[0],4); - } - else{ - int inputtostack[2]={EffectivePressureHydrostepEnum,SedimentHeadHydrostepEnum}; - int stackedinput[2]={EffectivePressureStackedEnum,SedimentHeadStackedEnum}; - int averagedinput[2]={EffectivePressureEnum,SedimentHeadEnum}; - femmodel->InitMeanOutputx(&stackedinput[0],2); - while(hydrotimeSumOutputx(&inputtostack[0],&stackedinput[0],2); - } - femmodel->AverageSumOutputx(&stackedinput[0],&averagedinput[0],2); - } + /*Proceed now to heads computations*/ + solutionsequence_hydro_nonlinear(femmodel); + /*If we have a sub-timestep we stack the variables here*/ + femmodel->SumOutputx(&inputtostack[0],&stackedinput[0],numaveragedinput); + } + /*Reseting to global time*/ + femmodel->parameters->SetParam(time,TimeEnum); + femmodel->AverageSumOutputx(&stackedinput[0],&averagedinput[0],numaveragedinput); } else{ + if(issmb){ + if(VerboseSolution()) _printf0_(" computing mass balance\n"); + SmbAnalysis* analysis = new SmbAnalysis(); + analysis->Core(femmodel); + delete analysis; + } InputDuplicatex(femmodel,SedimentHeadHydrostepEnum,SedimentHeadOldEnum); if (isefficientlayer){ @@ -116,9 +141,21 @@ } /*Proceed now to heads computations*/ + if(VerboseSolution()) _printf0_(" computing water heads\n"); solutionsequence_hydro_nonlinear(femmodel); } } - } - else if (hydrology_model==HydrologyshaktiEnum){ + else{ + /* If everything is frozen we still need smb */ + if(issmb){ + if(VerboseSolution()) _printf0_(" computing mass balance\n"); + SmbAnalysis* analysis = new SmbAnalysis(); + analysis->Core(femmodel); + delete analysis; + } + } + } + + /*Using the SHAKTI model*/ + else if (hydrology_model==HydrologyshaktiEnum){ femmodel->SetCurrentConfiguration(HydrologyShaktiAnalysisEnum); InputDuplicatex(femmodel,HydrologyHeadEnum,HydrologyHeadOldEnum); @@ -129,4 +166,6 @@ delete analysis; } + + /*Using the PISM hydrology model*/ else if (hydrology_model==HydrologypismEnum){ femmodel->SetCurrentConfiguration(HydrologyPismAnalysisEnum); @@ -143,4 +182,8 @@ if(hydrology_model==HydrologydcEnum && ThawedNodes==0){ if(VerboseSolution()) _printf0_(" No thawed node hydro is skiped \n");} + else if (hydrology_model==HydrologydcEnum && issmb){ + femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs); + femmodel->RequestedOutputsx(&femmodel->results,requested_smb_outputs,smboutputs); + } else{ femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs); @@ -154,5 +197,12 @@ xDelete(requested_outputs); } - + if(issmb){ + if(smboutputs){ + for(int i=0;i(requested_smb_outputs[i]); + } + xDelete(requested_smb_outputs); + } + } /*End profiler*/ femmodel->profiler->Stop(HYDROLOGYCORE); Index: /issm/trunk-jpl/src/c/cores/transient_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/transient_core.cpp (revision 23664) +++ /issm/trunk-jpl/src/c/cores/transient_core.cpp (revision 23665) @@ -74,5 +74,5 @@ #if defined(_HAVE_BAMG_) && !defined(_HAVE_AD_) if(amr_frequency){ - femmodel->parameters->FindParam(&amr_restart,AmrRestartEnum); + femmodel->parameters->FindParam(&amr_restart,AmrRestartEnum); if(amr_restart) femmodel->ReMesh(); } @@ -379,9 +379,14 @@ thermal_core(femmodel); } - - /*shifting smb position to have runoff value*/ - if(issmb) smb_core(femmodel); - - if(ishydrology) hydrology_core(femmodel); + /* Using Hydrology dc coupled we need to compute smb in the hydrology inner time loop*/ + if(ishydrology){ + int hydrology_model; + hydrology_core(femmodel); + femmodel->parameters->FindParam(&hydrology_model,HydrologyModelEnum); + if(hydrology_model!=HydrologydcEnum && issmb)smb_core(femmodel); + } + else{ + if(issmb) smb_core(femmodel); + } if(isstressbalance && (step%sb_coupling_frequency==0 || step==1) ) stressbalance_core(femmodel); Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 23664) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 23665) @@ -427,9 +427,4 @@ 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));