Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 8560) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 8561) @@ -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,hydro_CR,hydro_n,rho_water,g; + double hydro_p,hydro_q,hydro_CR,hydro_n,rho_water,g,gamma; double DL_scalar1; double DL_scalar2; @@ -1005,13 +1005,11 @@ /*retrieve material parameters: */ -// 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)); + g =matpar->GetG(); + rho_water=matpar->GetRhoWater(); + hydro_p =matpar->GetHydroP(); + hydro_q =matpar->GetHydroQ(); + hydro_CR =matpar->GetHydroCR(); + hydro_n =matpar->GetHydroN(); + gamma=1/(hydro_n*pow(hydro_CR,hydro_p)*pow((rho_water*g),hydro_q)); /* Start looping on the number of gaussian points: */ @@ -1024,17 +1022,17 @@ GetNodalFunctions(&L[0],gauss); GetNodalFunctionsDerivatives(&dL[0][0],&xyz_list[0][0],gauss); - //printf("NUMVERTICES=%d, Jdettria=%f\n",NUMVERTICES,Jdettria); //BK + //printf("NUMVERTICES=%d, Jdettria=%f\n",NUMVERTICES,Jdettria); //BK for(i=0;iweight*Jdettria*hydro_gamma*K[0]*dt; //don't forget the - - DL_scalar3=-gauss->weight*Jdettria*hydro_gamma*K[1]*dt; //don't forget the - + DL_scalar2=-gauss->weight*Jdettria*gamma*K[0]*dt; //don't forget the - + DL_scalar3=-gauss->weight*Jdettria*gamma*K[1]*dt; //don't forget the - } else{ - DL_scalar2=-gauss->weight*Jdettria*hydro_gamma*K[0]; //don't forget the - - DL_scalar3=-gauss->weight*Jdettria*hydro_gamma*K[1]; //don't forget the - + DL_scalar2=-gauss->weight*Jdettria*gamma*K[0]; //don't forget the - + DL_scalar3=-gauss->weight*Jdettria*gamma*K[1]; //don't forget the - } @@ -2105,6 +2103,6 @@ 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 + //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]; @@ -2715,5 +2713,4 @@ watercolumn_input->GetParameterValue(&watercolumn,gauss); values[i]=watercolumn; - // if((values[i])<0)values[i]=0; } @@ -2803,13 +2800,10 @@ //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);} - - -} - + //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*/ @@ -3940,5 +3934,5 @@ values[i]=solution[doflist[i]]; if(isnan(values[i])) _error_("NaN found in solution vector"); - // if (values[i]numforcings,iomodel_handle,"numforcings"); - /*!Get hydrology parameters: */ - // IoModelFetchData(&this->hydro_gamma,iomodel_handle,"hydro_gamma"); + /*!Get hydrology parameters: */ IoModelFetchData(&this->hydro_kn,iomodel_handle,"hydro_kn"); IoModelFetchData(&this->hydro_p,iomodel_handle,"hydro_p"); @@ -233,4 +232,5 @@ 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 8560) +++ /issm/trunk/src/c/objects/IoModel.h (revision 8561) @@ -227,11 +227,10 @@ 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; + /*hydrology parameters: */ + double hydro_kn; + double hydro_p; + double hydro_q; + double hydro_CR; + double hydro_n; /*}}}*/ /*Methods: {{{1*/ Index: /issm/trunk/src/c/objects/Materials/Matpar.cpp =================================================================== --- /issm/trunk/src/c/objects/Materials/Matpar.cpp (revision 8560) +++ /issm/trunk/src/c/objects/Materials/Matpar.cpp (revision 8561) @@ -36,5 +36,4 @@ double matpar_thermal_exchange_velocity; double matpar_g; - double matpar_gamma; double matpar_kn; @@ -50,5 +49,4 @@ matpar_mixed_layer_capacity=iomodel->mixed_layer_capacity; matpar_thermal_exchange_velocity=iomodel->thermal_exchange_velocity; -// matpar_gamma=iomodel->hydro_gamma; matpar_kn=iomodel->hydro_kn; @@ -65,5 +63,4 @@ this->thermal_exchange_velocity=matpar_thermal_exchange_velocity; this->g=matpar_g; - this->gamma=matpar_gamma; this->kn=matpar_kn; this->hydro_p=iomodel->hydro_p; @@ -96,5 +93,4 @@ printf(" thermal_exchange_velocity: %g\n",thermal_exchange_velocity); printf(" g: %g\n",g); - printf(" gamma: %g\n",gamma); printf(" kn: %g\n",kn); return; @@ -117,5 +113,4 @@ printf(" thermal_exchange_velocity: %g\n",thermal_exchange_velocity); printf(" g: %g\n",g); - printf(" gamma: %g\n",gamma); printf(" kn: %g\n",kn); return; @@ -159,5 +154,4 @@ memcpy(marshalled_dataset,&thermal_exchange_velocity,sizeof(thermal_exchange_velocity));marshalled_dataset+=sizeof(thermal_exchange_velocity); memcpy(marshalled_dataset,&g,sizeof(g));marshalled_dataset+=sizeof(g); - memcpy(marshalled_dataset,&gamma,sizeof(gamma));marshalled_dataset+=sizeof(gamma); memcpy(marshalled_dataset,&kn,sizeof(kn));marshalled_dataset+=sizeof(kn); @@ -181,5 +175,4 @@ sizeof(thermal_exchange_velocity)+ sizeof(g)+ - sizeof(gamma)+ sizeof(kn)+ sizeof(int); //sizeof(int) for enum type @@ -209,5 +202,4 @@ memcpy(&thermal_exchange_velocity,marshalled_dataset,sizeof(thermal_exchange_velocity));marshalled_dataset+=sizeof(thermal_exchange_velocity); memcpy(&g,marshalled_dataset,sizeof(g));marshalled_dataset+=sizeof(g); - memcpy(&gamma,marshalled_dataset,sizeof(gamma));marshalled_dataset+=sizeof(gamma); memcpy(&kn,marshalled_dataset,sizeof(kn));marshalled_dataset+=sizeof(kn); @@ -393,9 +385,4 @@ double Matpar::GetThermalExchangeVelocity(){ return thermal_exchange_velocity; -} -/*}}}1*/ -/*FUNCTION Matpar::GetGamma {{{1*/ -double Matpar::GetGamma(){ - return gamma; } /*}}}1*/ Index: /issm/trunk/src/c/objects/Materials/Matpar.h =================================================================== --- /issm/trunk/src/c/objects/Materials/Matpar.h (revision 8560) +++ /issm/trunk/src/c/objects/Materials/Matpar.h (revision 8561) @@ -30,13 +30,11 @@ /*hydrology: */ double kn; - double gamma; double hydro_p; double hydro_q; - double hydro_CR; + double hydro_CR; double hydro_n; + public: - Matpar(); - Matpar(int matpar_id, IoModel* iomodel); ~Matpar(); @@ -83,10 +81,9 @@ double GetMeltingPoint(); double GetReferenceTemperature(); - double GetGamma(); double GetKn(); double GetHydroP(); double GetHydroQ(); - double GetHydroCR(); - double GetHydroN(); + 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 8560) +++ /issm/trunk/src/c/solutions/hydrology_core.cpp (revision 8561) @@ -46,5 +46,4 @@ time=(i+1)*dt; -// printf("time step=%d=============================\n",i+1);//bk /*call hydrology_core step: */ hydrology_core_step(femmodel,i,time);