Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 8548) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 8549) @@ -991,5 +991,5 @@ double Ke_gg_gaussian2[numdof][numdof] ={0.0}; double Ke_gg_gaussian3[numdof][numdof] ={0.0}; - double hydro_gamma,hydro_p,hydro_q; + double hydro_gamma,hydro_p,hydro_q,hydro_CR,hydro_n,rho_water,g; double DL_scalar1; double DL_scalar2; @@ -1005,8 +1005,14 @@ /*retrieve material parameters: */ - hydro_gamma=matpar->GetGamma(); +// hydro_gamma=matpar->GetGamma(); + g =matpar->GetG(); + rho_water =matpar->GetRhoWater(); hydro_p =matpar->GetHydroP(); hydro_q =matpar->GetHydroQ(); - + hydro_CR =matpar->GetHydroCR(); + hydro_n =matpar->GetHydroN(); + + hydro_gamma=1/(hydro_n*pow(hydro_CR,hydro_p)*pow((rho_water*g),hydro_q)); + /* Start looping on the number of gaussian points: */ gauss=new GaussTria(2); @@ -1018,11 +1024,18 @@ GetNodalFunctions(&L[0],gauss); GetNodalFunctionsDerivatives(&dL[0][0],&xyz_list[0][0],gauss); - + //printf("NUMVERTICES=%d, Jdettria=%f\n",NUMVERTICES,Jdettria); //BK for(i=0;iweight*Jdettria; @@ -1049,4 +1062,5 @@ for(i=0;ivalues[i*numdof+j]+=Ke_gg_gaussian2[i][j]; for(i=0;ivalues[i*numdof+j]+=Ke_gg_gaussian3[i][j]; + } @@ -2091,5 +2105,7 @@ basal_melting_input->GetParameterValue(&basal_melting_g,gauss); old_watercolumn_input->GetParameterValue(&old_watercolumn_g,gauss); - +// printf("ig=%d, old_watercolumn_g=%f,Jdettria=%f\n",ig,old_watercolumn_g,Jdettria);//bk +// if(old_watercolumn_g<0) printf("old_watercolumn_g=%g,Jdettria=%f\n",old_watercolumn_g,Jdettria);//bk + if(dt)for(i=0;ivalues[i]+=Jdettria*gauss->weight*(old_watercolumn_g+dt*basal_melting_g)*L[i]; else for(i=0;ivalues[i]+=Jdettria*gauss->weight*basal_melting_g*L[i]; @@ -2699,5 +2715,5 @@ watercolumn_input->GetParameterValue(&watercolumn,gauss); values[i]=watercolumn; - if((values[i])<0)values[i]=0; + // if((values[i])<0)values[i]=0; } @@ -2785,6 +2801,15 @@ K[0]=pow(w,2)*(rho_ice*g*dsdx+(rho_water/rho_ice-1)*rho_ice*g*dbdx) - rho_ice * g * kn* w * (dsdx - dbdx ) * surface_slope; K[1]=pow(w,2)*(rho_ice*g*dsdy+(rho_water/rho_ice-1)*rho_ice*g*dbdy) - rho_ice * g * kn *w * (dsdy - dbdy ) * surface_slope; - -} + + //bk + // K[0]=fabs(pow(w,2)*(rho_ice*g*dsdx+(rho_water/rho_ice-1)*rho_ice*g*dbdx) - rho_ice * g * kn* w * (dsdx - dbdx ) * surface_slope); + // K[1]=fabs(pow(w,2)*(rho_ice*g*dsdy+(rho_water/rho_ice-1)*rho_ice*g*dbdy) - rho_ice * g * kn *w * (dsdy - dbdy ) * surface_slope); +// printf("K[0]=%g,K[1]=%g,w=%g\n",K[0],K[1],w); //bk +// if (w<0) {printf("negative w!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n"); +// printf("K[0]=%g,K[1]=%g,w=%g\n",K[0],K[1],w);} + + +} + /*}}}*/ /*FUNCTION Tria::Gradj {{{1*/ @@ -3915,5 +3940,6 @@ values[i]=solution[doflist[i]]; if(isnan(values[i])) _error_("NaN found in solution vector"); - //if (values[i]beta,iomodel_handle,"beta"); - IoModelFetchData(&this->hydro_gamma,iomodel_handle,"hydro_gamma"); - IoModelFetchData(&this->hydro_kn,iomodel_handle,"hydro_kn"); - IoModelFetchData(&this->hydro_p,iomodel_handle,"hydro_p"); - IoModelFetchData(&this->hydro_q,iomodel_handle,"hydro_q"); IoModelFetchData(&this->meltingpoint,iomodel_handle,"meltingpoint"); IoModelFetchData(&this->referencetemperature,iomodel_handle,"referencetemperature"); @@ -230,4 +226,11 @@ IoModelFetchData(&this->numforcings,iomodel_handle,"numforcings"); + /*!Get hydrology parameters: */ +// IoModelFetchData(&this->hydro_gamma,iomodel_handle,"hydro_gamma"); + IoModelFetchData(&this->hydro_kn,iomodel_handle,"hydro_kn"); + IoModelFetchData(&this->hydro_p,iomodel_handle,"hydro_p"); + IoModelFetchData(&this->hydro_q,iomodel_handle,"hydro_q"); + IoModelFetchData(&this->hydro_CR,iomodel_handle,"hydro_CR"); + IoModelFetchData(&this->hydro_n,iomodel_handle,"hydro_n"); /*qmu: */ if(this->qmu_analysis){ Index: /issm/trunk/src/c/objects/IoModel.h =================================================================== --- /issm/trunk/src/c/objects/IoModel.h (revision 8548) +++ /issm/trunk/src/c/objects/IoModel.h (revision 8549) @@ -170,8 +170,4 @@ /*thermal parameters: */ double beta; - double hydro_gamma; - double hydro_kn; - double hydro_p; - double hydro_q; double meltingpoint; double referencetemperature; @@ -230,5 +226,13 @@ int loadcounter; //keep track of how many loads are being created in each analysis type int constraintcounter; //keep track of how many constraints are being created in each analysis type - /*}}}*/ + + /*hydrology parameters: */ +// double hydro_gamma; + double hydro_kn; + double hydro_p; + double hydro_q; + double hydro_CR; + double hydro_n; + /*}}}*/ /*Methods: {{{1*/ ~IoModel(); Index: /issm/trunk/src/c/objects/Materials/Matpar.cpp =================================================================== --- /issm/trunk/src/c/objects/Materials/Matpar.cpp (revision 8548) +++ /issm/trunk/src/c/objects/Materials/Matpar.cpp (revision 8549) @@ -50,5 +50,5 @@ matpar_mixed_layer_capacity=iomodel->mixed_layer_capacity; matpar_thermal_exchange_velocity=iomodel->thermal_exchange_velocity; - matpar_gamma=iomodel->hydro_gamma; +// matpar_gamma=iomodel->hydro_gamma; matpar_kn=iomodel->hydro_kn; @@ -69,5 +69,6 @@ this->hydro_p=iomodel->hydro_p; this->hydro_q=iomodel->hydro_q; - + this->hydro_CR=iomodel->hydro_CR; + this->hydro_n=iomodel->hydro_n; } /*}}}1*/ @@ -414,4 +415,14 @@ } /*}}}1*/ +/*FUNCTION Matpar::GetHydroCR {{{1*/ +double Matpar::GetHydroCR(){ + return hydro_CR; +} +/*}}}1*/ +/*FUNCTION Matpar::GetHydroN {{{1*/ +double Matpar::GetHydroN(){ + return hydro_n; +} +/*}}}1*/ /*FUNCTION Matpar::TMeltingPoint {{{1*/ double Matpar::TMeltingPoint(double pressure){ Index: /issm/trunk/src/c/objects/Materials/Matpar.h =================================================================== --- /issm/trunk/src/c/objects/Materials/Matpar.h (revision 8548) +++ /issm/trunk/src/c/objects/Materials/Matpar.h (revision 8549) @@ -33,5 +33,6 @@ double hydro_p; double hydro_q; - + double hydro_CR; + double hydro_n; public: @@ -86,4 +87,6 @@ double GetHydroP(); double GetHydroQ(); + double GetHydroCR(); + double GetHydroN(); double TMeltingPoint(double pressure); double PureIceEnthalpy(double pressure); Index: /issm/trunk/src/c/solutions/hydrology_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/hydrology_core.cpp (revision 8548) +++ /issm/trunk/src/c/solutions/hydrology_core.cpp (revision 8549) @@ -43,7 +43,8 @@ for(i=0;i