Index: /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp (revision 27924) +++ /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp (revision 27925) @@ -146,4 +146,7 @@ iomodel->FetchDataToInput(inputs,elements,"md.hydrology.bump_height",HydrologyBumpHeightEnum); iomodel->FetchDataToInput(inputs,elements,"md.hydrology.sheet_conductivity",HydrologySheetConductivityEnum); + iomodel->FetchDataToInput(inputs,elements,"md.hydrology.omega",HydrologyOmegaEnum); + iomodel->FetchDataToInput(inputs,elements,"md.hydrology.sheet_alpha",HydrologySheetAlphaEnum); + iomodel->FetchDataToInput(inputs,elements,"md.hydrology.sheet_beta",HydrologySheetBetaEnum); iomodel->FetchDataToInput(inputs,elements,"md.hydrology.channel_conductivity",HydrologyChannelConductivityEnum); iomodel->FetchDataToInput(inputs,elements,"md.hydrology.neumannflux",HydrologyNeumannfluxEnum); @@ -183,4 +186,6 @@ parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.melt_flag",HydrologyMeltFlagEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.channel_sheet_width",HydrologyChannelSheetWidthEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.channel_alpha",HydrologyChannelAlphaEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.channel_beta",HydrologyChannelBetaEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.islaminar",HydrologyIsLaminarEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.englacial_void_ratio",HydrologyEnglacialVoidRatioEnum)); @@ -217,10 +222,11 @@ /*Intermediaries */ IssmDouble Jdet,dphi[3],h,k; + IssmDouble alpha,beta,omega,h_r; IssmDouble A,B,n,phi_old,phi,phi_0,H,b,v1; IssmDouble* xyz_list = NULL; /*Hard coded coefficients*/ - const IssmPDouble alpha = 5./4.; - const IssmPDouble beta = 3./2.; + const IssmPDouble ALPHA = 5./4.; + const IssmPDouble BETA = 3./2.; /*Fetch number of nodes and dof for this finite element*/ @@ -236,10 +242,17 @@ /*Get all inputs and parameters*/ + int islaminar; + element->FindParam(&islaminar,HydrologyIsLaminarEnum); IssmDouble dt = element->FindParam(TimesteppingTimeStepEnum); IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum); + IssmDouble mu_water = element->FindParam(MaterialsMuWaterEnum); IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum); IssmDouble g = element->FindParam(ConstantsGEnum); IssmDouble e_v = element->FindParam(HydrologyEnglacialVoidRatioEnum); + Input* hr_input = element->GetInput(HydrologyBumpHeightEnum);_assert_(hr_input); Input* k_input = element->GetInput(HydrologySheetConductivityEnum);_assert_(k_input); + Input* alpha_input = element->GetInput(HydrologySheetAlphaEnum);_assert_(alpha_input); + Input* beta_input = element->GetInput(HydrologySheetBetaEnum);_assert_(beta_input); + Input* omega_input = element->GetInput(HydrologyOmegaEnum);_assert_(omega_input); Input* phi_input = element->GetInput(HydraulicPotentialEnum); _assert_(phi_input); Input* h_input = element->GetInput(HydrologySheetThicknessEnum); _assert_(h_input); @@ -261,6 +274,10 @@ h_input->GetInputValue(&h,gauss); k_input->GetInputValue(&k,gauss); + omega_input->GetInputValue(&omega,gauss); + alpha_input->GetInputValue(&alpha,gauss); + beta_input->GetInputValue(&beta,gauss); B_input->GetInputValue(&B,gauss); n_input->GetInputValue(&n,gauss); + hr_input->GetInputValue(&h_r,gauss); b_input->GetInputValue(&b,gauss); H_input->GetInputValue(&H,gauss); @@ -268,10 +285,21 @@ /*Hard code B*/ B = Cuffey(273.15-2); - + /*Get norm of gradient of hydraulic potential and make sure it is >0*/ IssmDouble normgradphi = sqrt(dphi[0]*dphi[0] + dphi[1]*dphi[1]); if(normgradphi < AEPS) normgradphi = AEPS; - - IssmDouble coeff = k*pow(h,alpha)*pow(normgradphi,beta-2.); + + /*Use Laminar flow if specified*/ + IssmDouble nu = mu_water/rho_water; + IssmDouble coeff; + if(islaminar==1 && omega>=AEPS){ + IssmDouble hratio = fabs(h/h_r); + IssmDouble coarg = 1. + 4.*pow(hratio,3-2*alpha)*omega*k*pow(h,3)*normgradphi/nu; + coeff = nu/2./omega*pow(hratio,2*alpha-3) * (-1 + pow(coarg, 0.5))/normgradphi; + } + else { + /*If omega is zero, use standard model, otherwise transition model*/ + coeff = k*pow(h,ALPHA)*pow(normgradphi,BETA-2.); + } /*Diffusive term*/ @@ -440,10 +468,11 @@ /*Intermediaries*/ IssmDouble dphi[3],h,k,phi; + IssmDouble alpha,beta,omega,h_r; IssmDouble oceanLS,iceLS; IssmDouble* xyz_list = NULL; /*Hard coded coefficients*/ - const IssmPDouble alpha = 5./4.; - const IssmPDouble beta = 3./2.; + const IssmPDouble ALPHA = 5./4.; + const IssmPDouble BETA = 3./2.; /*Fetch number vertices for this element*/ @@ -466,7 +495,15 @@ /*Retrieve all inputs and parameters*/ + int islaminar; + element->FindParam(&islaminar,HydrologyIsLaminarEnum); element->GetVerticesCoordinates(&xyz_list); + IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum); + IssmDouble mu_water = element->FindParam(MaterialsMuWaterEnum); Input *k_input = element->GetInput(HydrologySheetConductivityEnum); _assert_(k_input); + Input *omega_input = element->GetInput(HydrologyOmegaEnum); _assert_(omega_input); + Input *alpha_input = element->GetInput(HydrologySheetAlphaEnum); _assert_(alpha_input); + Input *beta_input = element->GetInput(HydrologySheetBetaEnum); _assert_(bet_input); Input *phi_input = element->GetInput(HydraulicPotentialEnum); _assert_(phi_input); + Input *hr_input = element->GetInput(HydrologyBumpHeightEnum); _assert_(hr_input); Input *h_input = element->GetInput(HydrologySheetThicknessEnum); _assert_(h_input); Input *oceanLS_input = element->GetInput(MaskOceanLevelsetEnum); _assert_(oceanLS_input); @@ -482,5 +519,9 @@ phi_input->GetInputValue(&phi,gauss); h_input->GetInputValue(&h,gauss); + hr_input->GetInputValue(&h_r,gauss); k_input->GetInputValue(&k,gauss); + omega_input->GetInputValue(&omega,gauss); + alpha_input->GetInputValue(&alpha,gauss); + beta_input->GetInputValue(&beta,gauss); oceanLS_input->GetInputValue(&oceanLS,gauss); iceLS_input->GetInputValue(&iceLS,gauss); @@ -496,6 +537,18 @@ IssmDouble normgradphi = sqrt(dphi[0]*dphi[0] + dphi[1]*dphi[1]); if(normgradphi < AEPS) normgradphi = AEPS; - - IssmDouble coeff = k*pow(h,alpha)*pow(normgradphi,beta-2.)/max(AEPS,h); // divide by h to get speed instead of discharge + + /*If omega is zero, use standard model, otherwise transition model*/ + IssmDouble nu = mu_water/rho_water; + IssmDouble coeff; + if(islaminar==1 && omega>=AEPS){ + IssmDouble hratio = fabs(h/h_r); + IssmDouble coarg = 1. + 4.*pow(hratio,3-2*alpha)*omega*k*pow(h,3)*normgradphi/nu; + coeff = nu/2./omega*pow(hratio,2*alpha-3) * (-1 + pow(coarg, 0.5))/normgradphi/max(AEPS,h); // divide by h to get speed instead of discharge + coeff = k*pow(h,alpha)*pow(normgradphi,beta-2.)/max(AEPS,h); // divide by h to get speed instead of discharge + } + else { + coeff = k*pow(h,ALPHA)*pow(normgradphi,BETA-2.)/max(AEPS,h); // divide by h to get speed instead of discharge + } + vx[iv] = -coeff*dphi[0]; Index: /issm/trunk-jpl/src/c/classes/Loads/Channel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Channel.cpp (revision 27924) +++ /issm/trunk-jpl/src/c/classes/Loads/Channel.cpp (revision 27925) @@ -367,4 +367,5 @@ IssmDouble Jdet,v1,qc,fFactor,Afactor,Bfactor,Xifactor; IssmDouble A,B,n,phi_old,phi,phi_0,dPw,ks,kc,Ngrad; + IssmDouble alpha_s,beta_s,omega,h_r; IssmDouble H,h,b,dphi[2],dphids,dphimds,db[2],dbds; IssmDouble xyz_list[NUMVERTICES][3]; @@ -382,5 +383,8 @@ GetVerticesCoordinates(&xyz_list_tria[0][0],tria->vertices,3); + int islaminar; + element->FindParam(&islaminar,HydrologyIsLaminarEnum); IssmDouble L = element->FindParam(MaterialsLatentheatEnum); + IssmDouble mu_water = element->FindParam(MaterialsMuWaterEnum); IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum); IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum); @@ -388,4 +392,6 @@ IssmDouble lc = element->FindParam(HydrologyChannelSheetWidthEnum); IssmDouble c_t = element->FindParam(HydrologyPressureMeltCoefficientEnum); + IssmDouble alpha_c = element->FindParam(HydrologyChannelAlphaEnum); + IssmDouble beta_c = element->FindParam(HydrologyChannelBetaEnum); Input* h_input = element->GetInput(HydrologySheetThicknessEnum); _assert_(h_input); @@ -396,4 +402,8 @@ Input* ks_input = element->GetInput(HydrologySheetConductivityEnum); _assert_(ks_input); Input* kc_input = element->GetInput(HydrologyChannelConductivityEnum); _assert_(kc_input); + Input* omega_input = element->GetInput(HydrologyOmegaEnum); _assert_(ks_input); + Input* alpha_s_input= element->GetInput(HydrologySheetAlphaEnum); _assert_(alpha_s_input); + Input* beta_s_input = element->GetInput(HydrologySheetBetaEnum); _assert_(beta_s_input); + Input* hr_input = element->GetInput(HydrologyBumpHeightEnum);_assert_(hr_input); Input* phi_input = element->GetInput(HydraulicPotentialEnum); _assert_(phi_input); @@ -422,4 +432,8 @@ ks_input->GetInputValue(&ks,gauss); kc_input->GetInputValue(&kc,gauss); + alpha_s_input->GetInputValue(&alpha_s,gauss); + beta_s_input->GetInputValue(&beta_s,gauss); + omega_input->GetInputValue(&omega,gauss); + hr_input->GetInputValue(&h_r,gauss); B_input->GetInputValue(&B,gauss); n_input->GetInputValue(&n,gauss); @@ -436,8 +450,19 @@ Ngrad = fabs(dphids); if(NgradS,ALPHA_C) * pow(Ngrad,BETA_C-2.); - IssmDouble Ks = ks * pow(h ,ALPHA_S) * pow(Ngrad,BETA_S-2.); + + /*Compute the effective conductivity Kc = k h^alpha |grad Phi|^{beta-2} (same for sheet) and use laminar if specified*/ + IssmDouble Kc; + IssmDouble Ks; + IssmDouble nu = mu_water/rho_water; + if(islaminar==1 && omega>=AEPS){ + IssmDouble hratio = h/h_r; + IssmDouble coarg = 1. + 4.*omega*pow(hratio,3-2*alpha_s)*ks*pow(h,3)*Ngrad/nu; + Ks = nu/2./omega*pow(hratio,2*alpha_s-3) * (-1 + pow(coarg, 0.5))/Ngrad; + Kc = kc * pow(this->S,alpha_c) * pow(Ngrad,beta_c-2.); + } + else { + Ks = ks*pow(h,ALPHA_S)*pow(Ngrad,BETA_S-2.); + Kc = kc * pow(this->S,ALPHA_C) * pow(Ngrad,BETA_C-2.); + } /*Approx. discharge in the sheet flowing folwing in the direction of the channel ofver a width lc*/ @@ -504,4 +529,5 @@ IssmDouble A,B,n,phi_old,phi,phi_0,dphimds,dphi[2]; IssmDouble H,h,b,db[2],dphids,qc,dPw,ks,kc,Ngrad; + IssmDouble alpha_s,beta_s,omega,h_r; IssmDouble xyz_list[NUMVERTICES][3]; IssmDouble xyz_list_tria[3][3]; @@ -516,5 +542,8 @@ GetVerticesCoordinates(&xyz_list_tria[0][0],tria->vertices,3); + int islaminar; + element->FindParam(&islaminar,HydrologyIsLaminarEnum); IssmDouble L = element->FindParam(MaterialsLatentheatEnum); + IssmDouble mu_water = element->FindParam(MaterialsMuWaterEnum); IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum); IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum); @@ -530,5 +559,9 @@ Input* ks_input = element->GetInput(HydrologySheetConductivityEnum); _assert_(ks_input); Input* kc_input = element->GetInput(HydrologyChannelConductivityEnum); _assert_(kc_input); + Input* omega_input = element->GetInput(HydrologyOmegaEnum); _assert_(omega_input); Input* phi_input = element->GetInput(HydraulicPotentialEnum); _assert_(phi_input); + Input* alpha_s_input= element->GetInput(HydrologySheetAlphaEnum);_assert_(alpha_s_input); + Input* beta_s_input = element->GetInput(HydrologySheetBetaEnum);_assert_(beta_s_input); + Input* hr_input = element->GetInput(HydrologyBumpHeightEnum);_assert_(hr_input); /*Get tangent vector*/ @@ -552,4 +585,5 @@ ks_input->GetInputValue(&ks,gauss); kc_input->GetInputValue(&kc,gauss); + omega_input->GetInputValue(&omega,gauss); B_input->GetInputValue(&B,gauss); n_input->GetInputValue(&n,gauss); @@ -557,4 +591,7 @@ b_input->GetInputValue(&b,gauss); H_input->GetInputValue(&H,gauss); + alpha_s_input->GetInputValue(&alpha_s,gauss); + beta_s_input->GetInputValue(&beta_s,gauss); + hr_input->GetInputValue(&h_r,gauss); /*Hard code B*/ @@ -568,6 +605,16 @@ if(Ngrad=AEPS){ + IssmDouble hratio = h/h_r; + IssmDouble nu = mu_water/rho_water; + IssmDouble coarg = 1. + 4.*omega*pow(hratio,3-2*alpha_s)*ks*pow(h,3)*Ngrad/nu; + Ks = nu/2./omega*pow(hratio,2*alpha_s-3) * (-1 + pow(coarg, 0.5))/Ngrad; + } + else { + Ks = ks * pow(h,ALPHA_S) * pow(Ngrad,BETA_S-2.); + } /*Approx. discharge in the sheet flowing folwing in the direction of the channel ofver a width lc*/ @@ -636,4 +683,5 @@ /*Intermediaries */ IssmDouble A,B,n,phi,phi_0,ks,kc,Ngrad; + IssmDouble alpha_s,beta_s,omega,h_r; IssmDouble H,h,b,dphi[2],dphids,dphimds,db[2],dbds; IssmDouble xyz_list[NUMVERTICES][3]; @@ -644,11 +692,16 @@ GetVerticesCoordinates(&xyz_list_tria[0][0],tria->vertices,3); + int islaminar; + element->FindParam(&islaminar,HydrologyIsLaminarEnum); IssmDouble L = element->FindParam(MaterialsLatentheatEnum); IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum); IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum); + IssmDouble mu_water = element->FindParam(MaterialsMuWaterEnum); IssmDouble g = element->FindParam(ConstantsGEnum); IssmDouble lc = element->FindParam(HydrologyChannelSheetWidthEnum); IssmDouble c_t = element->FindParam(HydrologyPressureMeltCoefficientEnum); IssmDouble dt = element->FindParam(TimesteppingTimeStepEnum); + IssmDouble alpha_c = element->FindParam(HydrologyChannelAlphaEnum); + IssmDouble beta_c = element->FindParam(HydrologyChannelBetaEnum); Input* h_input = element->GetInput(HydrologySheetThicknessEnum); _assert_(h_input); @@ -659,5 +712,9 @@ Input* ks_input = element->GetInput(HydrologySheetConductivityEnum); _assert_(ks_input); Input* kc_input = element->GetInput(HydrologyChannelConductivityEnum); _assert_(kc_input); + Input* omega_input = element->GetInput(HydrologyOmegaEnum); _assert_(omega_input); Input* phi_input = element->GetInput(HydraulicPotentialEnum); _assert_(phi_input); + Input* alpha_s_input= element->GetInput(HydrologySheetAlphaEnum);_assert_(alpha_s_input); + Input* beta_s_input = element->GetInput(HydrologySheetBetaEnum); _assert_(beta_s_input); + Input* hr_input = element->GetInput(HydrologyBumpHeightEnum);_assert_(hr_input); /*Get tangent vector*/ @@ -674,4 +731,5 @@ ks_input->GetInputValue(&ks,gauss); kc_input->GetInputValue(&kc,gauss); + omega_input->GetInputValue(&omega,gauss); B_input->GetInputValue(&B,gauss); n_input->GetInputValue(&n,gauss); @@ -679,4 +737,8 @@ b_input->GetInputDerivativeValue(&db[0],&xyz_list_tria[0][0],gauss); H_input->GetInputValue(&H,gauss); + hr_input->GetInputValue(&h_r,gauss); + alpha_s_input->GetInputValue(&alpha_s,gauss); + beta_s_input->GetInputValue(&beta_s,gauss); + /*Hard code B*/ @@ -693,6 +755,15 @@ IssmDouble dPw = dphids - dphimds; - /*Approx. discharge in the sheet flowing folwing in the direction of the channel ofver a width lc*/ - IssmDouble qc = - ks * pow(h,ALPHA_S) * pow(Ngrad,BETA_S-2.) * dphids; + /*Approx. discharge in the sheet flowing folwing in the direction of the channel ofver a width lc, use laminar if necessary*/ + IssmDouble qc; + if (islaminar==1 && omega>=AEPS){ + IssmDouble hratio = h/h_r; + IssmDouble nu = mu_water/rho_water; + IssmDouble coarg = 1. + 4.*omega*pow(hratio,3-2*alpha_s)*ks*pow(h,3)*fabs(Ngrad)/nu; + qc = -nu/2./omega*pow(hratio,2*alpha_s-3) * (-1 + pow(coarg, 0.5))*dphids/Ngrad; + } + else { + qc = - ks * pow(h,ALPHA_S) * pow(Ngrad,BETA_S-2.) * dphids; + } /*Ice rate factor*/ Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 27924) +++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 27925) @@ -261,4 +261,6 @@ syn keyword cConstant HydrologyarmaTimestepEnum syn keyword cConstant HydrologyAveragingEnum +syn keyword cConstant HydrologyChannelAlphaEnum +syn keyword cConstant HydrologyChannelBetaEnum syn keyword cConstant HydrologyCavitySpacingEnum syn keyword cConstant HydrologyChannelSheetWidthEnum @@ -896,5 +898,8 @@ syn keyword cConstant HydrologyMoulinInputEnum syn keyword cConstant HydrologyNeumannfluxEnum +syn keyword cConstant HydrologyOmegaEnum syn keyword cConstant HydrologyReynoldsEnum +syn keyword cConstant HydrologySheetAlphaEnum +syn keyword cConstant HydrologySheetBetaEnum syn keyword cConstant HydrologySheetConductivityEnum syn keyword cConstant HydrologySheetThicknessEnum Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 27924) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 27925) @@ -255,4 +255,6 @@ HydrologyarmaTimestepEnum, HydrologyAveragingEnum, + HydrologyChannelAlphaEnum, + HydrologyChannelBetaEnum, HydrologyCavitySpacingEnum, HydrologyChannelSheetWidthEnum, @@ -892,5 +894,8 @@ HydrologyMoulinInputEnum, HydrologyNeumannfluxEnum, + HydrologyOmegaEnum, HydrologyReynoldsEnum, + HydrologySheetAlphaEnum, + HydrologySheetBetaEnum, HydrologySheetConductivityEnum, HydrologySheetThicknessEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 27924) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 27925) @@ -263,4 +263,6 @@ case HydrologyarmaTimestepEnum : return "HydrologyarmaTimestep"; case HydrologyAveragingEnum : return "HydrologyAveraging"; + case HydrologyChannelAlphaEnum : return "HydrologyChannelAlpha"; + case HydrologyChannelBetaEnum : return "HydrologyChannelBeta"; case HydrologyCavitySpacingEnum : return "HydrologyCavitySpacing"; case HydrologyChannelSheetWidthEnum : return "HydrologyChannelSheetWidth"; @@ -898,5 +900,8 @@ case HydrologyMoulinInputEnum : return "HydrologyMoulinInput"; case HydrologyNeumannfluxEnum : return "HydrologyNeumannflux"; + case HydrologyOmegaEnum : return "HydrologyOmega"; case HydrologyReynoldsEnum : return "HydrologyReynolds"; + case HydrologySheetAlphaEnum : return "HydrologySheetAlpha"; + case HydrologySheetBetaEnum : return "HydrologySheetBeta"; case HydrologySheetConductivityEnum : return "HydrologySheetConductivity"; case HydrologySheetThicknessEnum : return "HydrologySheetThickness"; Index: /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim (revision 27924) +++ /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim (revision 27925) @@ -254,4 +254,6 @@ syn keyword juliaConstC HydrologyarmaTimestepEnum syn keyword juliaConstC HydrologyAveragingEnum +syn keyword juliaConstC HydrologyChannelAlphaEnum +syn keyword juliaConstC HydrologyChannelBetaEnum syn keyword juliaConstC HydrologyCavitySpacingEnum syn keyword juliaConstC HydrologyChannelSheetWidthEnum @@ -889,5 +891,8 @@ syn keyword juliaConstC HydrologyMoulinInputEnum syn keyword juliaConstC HydrologyNeumannfluxEnum +syn keyword juliaConstC HydrologyOmegaEnum syn keyword juliaConstC HydrologyReynoldsEnum +syn keyword juliaConstC HydrologySheetAlphaEnum +syn keyword juliaConstC HydrologySheetBetaEnum syn keyword juliaConstC HydrologySheetConductivityEnum syn keyword juliaConstC HydrologySheetThicknessEnum Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 27924) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 27925) @@ -269,4 +269,6 @@ else if (strcmp(name,"HydrologyarmaTimestep")==0) return HydrologyarmaTimestepEnum; else if (strcmp(name,"HydrologyAveraging")==0) return HydrologyAveragingEnum; + else if (strcmp(name,"HydrologyChannelAlpha")==0) return HydrologyChannelAlphaEnum; + else if (strcmp(name,"HydrologyChannelBeta")==0) return HydrologyChannelBetaEnum; else if (strcmp(name,"HydrologyCavitySpacing")==0) return HydrologyCavitySpacingEnum; else if (strcmp(name,"HydrologyChannelSheetWidth")==0) return HydrologyChannelSheetWidthEnum; @@ -381,10 +383,10 @@ else if (strcmp(name,"MasstransportIsfreesurface")==0) return MasstransportIsfreesurfaceEnum; else if (strcmp(name,"MasstransportMinThickness")==0) return MasstransportMinThicknessEnum; - else if (strcmp(name,"MasstransportNumRequestedOutputs")==0) return MasstransportNumRequestedOutputsEnum; - else if (strcmp(name,"MasstransportPenaltyFactor")==0) return MasstransportPenaltyFactorEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"MasstransportRequestedOutputs")==0) return MasstransportRequestedOutputsEnum; + if (strcmp(name,"MasstransportNumRequestedOutputs")==0) return MasstransportNumRequestedOutputsEnum; + else if (strcmp(name,"MasstransportPenaltyFactor")==0) return MasstransportPenaltyFactorEnum; + else if (strcmp(name,"MasstransportRequestedOutputs")==0) return MasstransportRequestedOutputsEnum; else if (strcmp(name,"MasstransportStabilization")==0) return MasstransportStabilizationEnum; else if (strcmp(name,"MaterialsBeta")==0) return MaterialsBetaEnum; @@ -504,10 +506,10 @@ else if (strcmp(name,"SolidearthSettingsSealevelLoading")==0) return SolidearthSettingsSealevelLoadingEnum; else if (strcmp(name,"SolidearthSettingsGRD")==0) return SolidearthSettingsGRDEnum; - else if (strcmp(name,"SolidearthSettingsRunFrequency")==0) return SolidearthSettingsRunFrequencyEnum; - else if (strcmp(name,"SolidearthSettingsTimeAcc")==0) return SolidearthSettingsTimeAccEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"SolidearthSettingsHoriz")==0) return SolidearthSettingsHorizEnum; + if (strcmp(name,"SolidearthSettingsRunFrequency")==0) return SolidearthSettingsRunFrequencyEnum; + else if (strcmp(name,"SolidearthSettingsTimeAcc")==0) return SolidearthSettingsTimeAccEnum; + else if (strcmp(name,"SolidearthSettingsHoriz")==0) return SolidearthSettingsHorizEnum; else if (strcmp(name,"SolidearthSettingsMaxiter")==0) return SolidearthSettingsMaxiterEnum; else if (strcmp(name,"SolidearthSettingsGrdOcean")==0) return SolidearthSettingsGrdOceanEnum; @@ -627,10 +629,10 @@ else if (strcmp(name,"SmbStepsPerStep")==0) return SmbStepsPerStepEnum; else if (strcmp(name,"SmbSwIdx")==0) return SmbSwIdxEnum; - else if (strcmp(name,"SmbSWgrad")==0) return SmbSWgradEnum; - else if (strcmp(name,"SmbT0dry")==0) return SmbT0dryEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"SmbT0wet")==0) return SmbT0wetEnum; + if (strcmp(name,"SmbSWgrad")==0) return SmbSWgradEnum; + else if (strcmp(name,"SmbT0dry")==0) return SmbT0dryEnum; + else if (strcmp(name,"SmbT0wet")==0) return SmbT0wetEnum; else if (strcmp(name,"SmbTcIdx")==0) return SmbTcIdxEnum; else if (strcmp(name,"SmbTeThresh")==0) return SmbTeThreshEnum; @@ -750,10 +752,10 @@ else if (strcmp(name,"BasalforcingsSpatialDeepwaterElevation")==0) return BasalforcingsSpatialDeepwaterElevationEnum; else if (strcmp(name,"BasalforcingsSpatialDeepwaterMeltingRate")==0) return BasalforcingsSpatialDeepwaterMeltingRateEnum; - else if (strcmp(name,"BasalforcingsSpatialUpperwaterElevation")==0) return BasalforcingsSpatialUpperwaterElevationEnum; - else if (strcmp(name,"BasalforcingsSpatialUpperwaterMeltingRate")==0) return BasalforcingsSpatialUpperwaterMeltingRateEnum; else stage=7; } if(stage==7){ - if (strcmp(name,"BasalforcingsIsmip6BasinId")==0) return BasalforcingsIsmip6BasinIdEnum; + if (strcmp(name,"BasalforcingsSpatialUpperwaterElevation")==0) return BasalforcingsSpatialUpperwaterElevationEnum; + else if (strcmp(name,"BasalforcingsSpatialUpperwaterMeltingRate")==0) return BasalforcingsSpatialUpperwaterMeltingRateEnum; + else if (strcmp(name,"BasalforcingsIsmip6BasinId")==0) return BasalforcingsIsmip6BasinIdEnum; else if (strcmp(name,"BasalforcingsIsmip6Tf")==0) return BasalforcingsIsmip6TfEnum; else if (strcmp(name,"BasalforcingsIsmip6TfShelf")==0) return BasalforcingsIsmip6TfShelfEnum; @@ -873,10 +875,10 @@ else if (strcmp(name,"FrictionQ")==0) return FrictionQEnum; else if (strcmp(name,"FrictionSedimentCompressibilityCoefficient")==0) return FrictionSedimentCompressibilityCoefficientEnum; - else if (strcmp(name,"FrictionTillFrictionAngle")==0) return FrictionTillFrictionAngleEnum; - else if (strcmp(name,"FrictionWaterLayer")==0) return FrictionWaterLayerEnum; else stage=8; } if(stage==8){ - if (strcmp(name,"FrictionWaterPressure")==0) return FrictionWaterPressureEnum; + if (strcmp(name,"FrictionTillFrictionAngle")==0) return FrictionTillFrictionAngleEnum; + else if (strcmp(name,"FrictionWaterLayer")==0) return FrictionWaterLayerEnum; + else if (strcmp(name,"FrictionWaterPressure")==0) return FrictionWaterPressureEnum; else if (strcmp(name,"FrictionWaterPressureNoise")==0) return FrictionWaterPressureNoiseEnum; else if (strcmp(name,"Frictionf")==0) return FrictionfEnum; @@ -919,5 +921,8 @@ else if (strcmp(name,"HydrologyMoulinInput")==0) return HydrologyMoulinInputEnum; else if (strcmp(name,"HydrologyNeumannflux")==0) return HydrologyNeumannfluxEnum; + else if (strcmp(name,"HydrologyOmega")==0) return HydrologyOmegaEnum; else if (strcmp(name,"HydrologyReynolds")==0) return HydrologyReynoldsEnum; + else if (strcmp(name,"HydrologySheetAlpha")==0) return HydrologySheetAlphaEnum; + else if (strcmp(name,"HydrologySheetBeta")==0) return HydrologySheetBetaEnum; else if (strcmp(name,"HydrologySheetConductivity")==0) return HydrologySheetConductivityEnum; else if (strcmp(name,"HydrologySheetThickness")==0) return HydrologySheetThicknessEnum; @@ -993,13 +998,13 @@ else if (strcmp(name,"SamplingKappa")==0) return SamplingKappaEnum; else if (strcmp(name,"SamplingPhi")==0) return SamplingPhiEnum; - else if (strcmp(name,"SamplingTau")==0) return SamplingTauEnum; + else stage=9; + } + if(stage==9){ + if (strcmp(name,"SamplingTau")==0) return SamplingTauEnum; else if (strcmp(name,"Sealevel")==0) return SealevelEnum; else if (strcmp(name,"SealevelGRD")==0) return SealevelGRDEnum; else if (strcmp(name,"SatGraviGRD")==0) return SatGraviGRDEnum; else if (strcmp(name,"SealevelBarystaticMask")==0) return SealevelBarystaticMaskEnum; - else stage=9; - } - if(stage==9){ - if (strcmp(name,"SealevelBarystaticIceMask")==0) return SealevelBarystaticIceMaskEnum; + else if (strcmp(name,"SealevelBarystaticIceMask")==0) return SealevelBarystaticIceMaskEnum; else if (strcmp(name,"SealevelBarystaticIceWeights")==0) return SealevelBarystaticIceWeightsEnum; else if (strcmp(name,"SealevelBarystaticIceArea")==0) return SealevelBarystaticIceAreaEnum; @@ -1116,13 +1121,13 @@ else if (strcmp(name,"SmbEAir")==0) return SmbEAirEnum; else if (strcmp(name,"SmbEC")==0) return SmbECEnum; - else if (strcmp(name,"SmbECDt")==0) return SmbECDtEnum; + else stage=10; + } + if(stage==10){ + if (strcmp(name,"SmbECDt")==0) return SmbECDtEnum; else if (strcmp(name,"SmbECini")==0) return SmbECiniEnum; else if (strcmp(name,"SmbEla")==0) return SmbElaEnum; else if (strcmp(name,"SmbEvaporation")==0) return SmbEvaporationEnum; else if (strcmp(name,"SmbFAC")==0) return SmbFACEnum; - else stage=10; - } - if(stage==10){ - if (strcmp(name,"SmbGdn")==0) return SmbGdnEnum; + else if (strcmp(name,"SmbGdn")==0) return SmbGdnEnum; else if (strcmp(name,"SmbGdnini")==0) return SmbGdniniEnum; else if (strcmp(name,"SmbGsp")==0) return SmbGspEnum; @@ -1239,13 +1244,13 @@ else if (strcmp(name,"SurfaceArea")==0) return SurfaceAreaEnum; else if (strcmp(name,"SurfaceAverageVelMisfit")==0) return SurfaceAverageVelMisfitEnum; - else if (strcmp(name,"SurfaceCrevasse")==0) return SurfaceCrevasseEnum; + else stage=11; + } + if(stage==11){ + if (strcmp(name,"SurfaceCrevasse")==0) return SurfaceCrevasseEnum; else if (strcmp(name,"Surface")==0) return SurfaceEnum; else if (strcmp(name,"SurfaceOld")==0) return SurfaceOldEnum; else if (strcmp(name,"SurfaceLogVelMisfit")==0) return SurfaceLogVelMisfitEnum; else if (strcmp(name,"SurfaceLogVxVyMisfit")==0) return SurfaceLogVxVyMisfitEnum; - else stage=11; - } - if(stage==11){ - if (strcmp(name,"SurfaceObservation")==0) return SurfaceObservationEnum; + else if (strcmp(name,"SurfaceObservation")==0) return SurfaceObservationEnum; else if (strcmp(name,"SurfaceRelVelMisfit")==0) return SurfaceRelVelMisfitEnum; else if (strcmp(name,"SurfaceSlopeX")==0) return SurfaceSlopeXEnum; @@ -1362,13 +1367,13 @@ else if (strcmp(name,"Outputdefinition58")==0) return Outputdefinition58Enum; else if (strcmp(name,"Outputdefinition59")==0) return Outputdefinition59Enum; - else if (strcmp(name,"Outputdefinition5")==0) return Outputdefinition5Enum; + else stage=12; + } + if(stage==12){ + if (strcmp(name,"Outputdefinition5")==0) return Outputdefinition5Enum; else if (strcmp(name,"Outputdefinition60")==0) return Outputdefinition60Enum; else if (strcmp(name,"Outputdefinition61")==0) return Outputdefinition61Enum; else if (strcmp(name,"Outputdefinition62")==0) return Outputdefinition62Enum; else if (strcmp(name,"Outputdefinition63")==0) return Outputdefinition63Enum; - else stage=12; - } - if(stage==12){ - if (strcmp(name,"Outputdefinition64")==0) return Outputdefinition64Enum; + else if (strcmp(name,"Outputdefinition64")==0) return Outputdefinition64Enum; else if (strcmp(name,"Outputdefinition65")==0) return Outputdefinition65Enum; else if (strcmp(name,"Outputdefinition66")==0) return Outputdefinition66Enum; @@ -1485,5 +1490,8 @@ else if (strcmp(name,"DataSetParam")==0) return DataSetParamEnum; else if (strcmp(name,"DatasetInput")==0) return DatasetInputEnum; - else if (strcmp(name,"DebrisAnalysis")==0) return DebrisAnalysisEnum; + else stage=13; + } + if(stage==13){ + if (strcmp(name,"DebrisAnalysis")==0) return DebrisAnalysisEnum; else if (strcmp(name,"DebrisSolution")==0) return DebrisSolutionEnum; else if (strcmp(name,"DefaultAnalysis")==0) return DefaultAnalysisEnum; @@ -1491,8 +1499,5 @@ else if (strcmp(name,"Dense")==0) return DenseEnum; else if (strcmp(name,"DependentObject")==0) return DependentObjectEnum; - else stage=13; - } - if(stage==13){ - if (strcmp(name,"DepthAverageAnalysis")==0) return DepthAverageAnalysisEnum; + else if (strcmp(name,"DepthAverageAnalysis")==0) return DepthAverageAnalysisEnum; else if (strcmp(name,"DeviatoricStressErrorEstimator")==0) return DeviatoricStressErrorEstimatorEnum; else if (strcmp(name,"Divergence")==0) return DivergenceEnum; @@ -1608,5 +1613,8 @@ else if (strcmp(name,"Loads")==0) return LoadsEnum; else if (strcmp(name,"LoveAnalysis")==0) return LoveAnalysisEnum; - else if (strcmp(name,"LoveHf")==0) return LoveHfEnum; + else stage=14; + } + if(stage==14){ + if (strcmp(name,"LoveHf")==0) return LoveHfEnum; else if (strcmp(name,"LoveHt")==0) return LoveHtEnum; else if (strcmp(name,"LoveKernelsImag")==0) return LoveKernelsImagEnum; @@ -1614,8 +1622,5 @@ else if (strcmp(name,"LoveKf")==0) return LoveKfEnum; else if (strcmp(name,"LoveKt")==0) return LoveKtEnum; - else stage=14; - } - if(stage==14){ - if (strcmp(name,"LoveLf")==0) return LoveLfEnum; + else if (strcmp(name,"LoveLf")==0) return LoveLfEnum; else if (strcmp(name,"LoveLt")==0) return LoveLtEnum; else if (strcmp(name,"LoveTidalHt")==0) return LoveTidalHtEnum; @@ -1731,5 +1736,8 @@ else if (strcmp(name,"SMBpdd")==0) return SMBpddEnum; else if (strcmp(name,"SMBpddSicopolis")==0) return SMBpddSicopolisEnum; - else if (strcmp(name,"SMBsemic")==0) return SMBsemicEnum; + else stage=15; + } + if(stage==15){ + if (strcmp(name,"SMBsemic")==0) return SMBsemicEnum; else if (strcmp(name,"SSAApproximation")==0) return SSAApproximationEnum; else if (strcmp(name,"SSAFSApproximation")==0) return SSAFSApproximationEnum; @@ -1737,8 +1745,5 @@ else if (strcmp(name,"Scaled")==0) return ScaledEnum; else if (strcmp(name,"SealevelAbsolute")==0) return SealevelAbsoluteEnum; - else stage=15; - } - if(stage==15){ - if (strcmp(name,"SealevelEmotion")==0) return SealevelEmotionEnum; + else if (strcmp(name,"SealevelEmotion")==0) return SealevelEmotionEnum; else if (strcmp(name,"SealevelchangePolarMotionX")==0) return SealevelchangePolarMotionXEnum; else if (strcmp(name,"SealevelchangePolarMotionY")==0) return SealevelchangePolarMotionYEnum; Index: /issm/trunk-jpl/src/m/classes/hydrologyglads.m =================================================================== --- /issm/trunk-jpl/src/m/classes/hydrologyglads.m (revision 27924) +++ /issm/trunk-jpl/src/m/classes/hydrologyglads.m (revision 27925) @@ -11,4 +11,7 @@ cavity_spacing = 0.; bump_height = NaN; + omega = 0; % TH + sheet_alpha = 5.0/4.0; % TH + sheet_beta = 3.0/2.0; % TH %Channels @@ -16,5 +19,6 @@ channel_conductivity = NaN; channel_sheet_width = 0.; - islaminar = 0; + channel_alpha = 5.0/4.0; % TH + channel_beta = 3.0/2.0; % TH %Other @@ -25,4 +29,5 @@ requested_outputs = {}; melt_flag = 0; + islaminar = 0; end methods @@ -46,4 +51,5 @@ self.pressure_melt_coefficient = 7.5e-8; %K/Pa (See table 1 in Erder et al. 2013) self.cavity_spacing = 2.; %m + self.omega = 1./2000.; % TH %Channel parameters @@ -51,10 +57,10 @@ self.channel_conductivity = 5.e-2; %Dow's default, Table uses 0.1 self.channel_sheet_width = 2.; %m - self.islaminar = 0; %by default use GlaDS default turbulent code - %Other + %Otherself.omega = 1./2000.; % TH self.englacial_void_ratio = 1.e-5;% Dow's default, Table from Werder et al. uses 1e-3; self.requested_outputs={'default'}; self.melt_flag=0; + self.islaminar = 0; %by default use GlaDS default turbulent code end % }}} function md = checkconsistency(self,md,solution,analyses) % {{{ @@ -70,4 +76,7 @@ md = checkfield(md,'fieldname','hydrology.cavity_spacing','numel',[1],'>',0); md = checkfield(md,'fieldname','hydrology.bump_height','size',[md.mesh.numberofvertices 1],'>=',0,'NaN',1,'Inf',1); + md = checkfield(md,'fieldname','hydrology.omega', 'numel', [1], '>=', 0); % TH + md = checkfield(md,'fieldname','hydrology.sheet_alpha', 'numel', [1], '>', 0); % TH + md = checkfield(md,'fieldname','hydrology.sheet_beta', 'numel', [1], '>', 0); % TH %Channels @@ -75,5 +84,6 @@ md = checkfield(md,'fieldname','hydrology.channel_conductivity','size',[md.mesh.numberofvertices 1],'>=',0,'NaN',1,'Inf',1); md = checkfield(md,'fieldname','hydrology.channel_sheet_width','numel',[1],'>=',0); - md = checkfield(md,'fieldname','hydrology.islaminar','numel',[1],'values',[0 1]); + md = checkfield(md,'fieldname','hydrology.channel_alpha', 'numel', [1], '>', 0); % TH + md = checkfield(md,'fieldname','hydrology.channel_beta', 'numel', [1], '>', 0); % TH %Other @@ -84,4 +94,5 @@ md = checkfield(md,'fieldname','hydrology.requested_outputs','stringrow',1); md = checkfield(md,'fieldname','hydrology.melt_flag','numel',[1],'values',[0 1 2]); + md = checkfield(md,'fieldname','hydrology.islaminar','numel',[1],'values',[0 1]); if self.melt_flag==1 || self.melt_flag==2 md = checkfield(md,'fieldname','basalforcings.groundedice_melting_rate','NaN',1,'Inf',1,'timeseries',1); @@ -93,11 +104,15 @@ fielddisplay(self,'pressure_melt_coefficient','Pressure melt coefficient (c_t) [K Pa^-1]'); fielddisplay(self,'sheet_conductivity','sheet conductivity (k) [m^(7/4) kg^(-1/2)]'); + fielddisplay(self,'sheet_alpha','First sheet-flow exponent (alpha_s) []'); % TH + fielddisplay(self,'sheet_beta','Second sheet-flow exponent (beta_s) []'); % TH fielddisplay(self,'cavity_spacing','cavity spacing (l_r) [m]'); fielddisplay(self,'bump_height','typical bump height (h_r) [m]'); + fielddisplay(self,'omega','transition parameter (omega) []'); % TH disp(sprintf(' CHANNELS')); fielddisplay(self,'ischannels','Do we allow for channels? 1: yes, 0: no'); fielddisplay(self,'channel_conductivity','channel conductivity (k_c) [m^(3/2) kg^(-1/2)]'); + fielddisplay(self,'channel_alpha','First channel-flow exponent (alpha_s) []'); % TH + fielddisplay(self,'channel_beta','Second channel-flow exponent (beta_s) []'); % TH fielddisplay(self,'channel_sheet_width','channel sheet width [m]'); - fielddisplay(self,'islaminar','do we use laminar [1] or turbulent physics [0, default]'); disp(sprintf(' OTHER')); fielddisplay(self,'spcphi','Hydraulic potential Dirichlet constraints [Pa]'); @@ -107,4 +122,5 @@ fielddisplay(self,'requested_outputs','additional outputs requested'); fielddisplay(self,'melt_flag','User specified basal melt? 0: no (default), 1: use md.basalforcings.groundedice_melting_rate'); + fielddisplay(self,'islaminar','do we use laminar [1] or turbulent physics [0, default]'); end % }}} function marshall(self,prefix,md,fid) % {{{ @@ -120,4 +136,7 @@ WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','cavity_spacing','format','Double'); WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','bump_height','format','DoubleMat','mattype',1); + WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','omega','format','Double'); % TH + WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','sheet_alpha','format','Double'); % TH + WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','sheet_beta','format','Double'); % TH %Channels @@ -125,5 +144,6 @@ WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','channel_conductivity','format','DoubleMat','mattype',1); WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','channel_sheet_width','format','Double'); - WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','islaminar','format','Boolean'); + WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','channel_alpha','format','Double'); % TH + WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','channel_beta','format','Double'); % TH %Others @@ -133,4 +153,5 @@ WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','englacial_void_ratio','format','Double'); WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','melt_flag','format','Integer'); + WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','islaminar','format','Integer'); outputs = self.requested_outputs; pos = find(ismember(outputs,'default')); Index: /issm/trunk-jpl/src/m/classes/hydrologyglads.py =================================================================== --- /issm/trunk-jpl/src/m/classes/hydrologyglads.py (revision 27924) +++ /issm/trunk-jpl/src/m/classes/hydrologyglads.py (revision 27925) @@ -19,4 +19,7 @@ self.cavity_spacing = 0. self.bump_height = np.nan + self.omega = 0.; #TH + self.sheet_alpha = 5.0/4.0; # TH + self.sheet_beta = 3.0/2.0; # TH # Channels @@ -24,5 +27,6 @@ self.channel_conductivity = np.nan self.channel_sheet_width = 0. - self.islaminar = 0 + self.channel_alpha = 5.0/4.0; # TH + self.channel_beta = 3.0/2.0; # TH # Other @@ -33,4 +37,5 @@ self.requested_outputs = [] self.melt_flag = 0 + self.islaminar = 0 nargs = len(args) @@ -49,11 +54,15 @@ s += '{}\n'.format(fielddisplay(self, 'pressure_melt_coefficient', 'Pressure melt coefficient (c_t) [K Pa^ - 1]')) s += '{}\n'.format(fielddisplay(self, 'sheet_conductivity', 'sheet conductivity (k) [m^(7 / 4) kg^(- 1 / 2)]')) + s += '{}\n'.format(fielddisplay(self, 'sheet_alpha', 'First sheet-flow exponent (alpha_s) []')) #TH + s += '{}\n'.format(fielddisplay(self, 'sheet_beta', 'Second sheet-flow exponent (beta_s) []')) #TH s += '{}\n'.format(fielddisplay(self, 'cavity_spacing', 'cavity spacing (l_r) [m]')) s += '{}\n'.format(fielddisplay(self, 'bump_height', 'typical bump height (h_r) [m]')) + s += '{}\n'.format(fielddisplay(self, 'omega', 'transition parameter (omega) []')) #TH s = '\t--CHANNELS\n' s += '{}\n'.format(fielddisplay(self, 'ischannels', 'Do we allow for channels? 1: yes, 0: no')) s += '{}\n'.format(fielddisplay(self, 'channel_conductivity', 'channel conductivity (k_c) [m^(3 / 2) kg^(- 1 / 2)]')) s += '{}\n'.format(fielddisplay(self, 'channel_sheet_width', 'channel sheet width [m]')) - s += '{}\n'.format(fielddisplay(self, 'islaminar','do we use laminar [1] or turbulent physics [0, default]')) + s += '{}\n'.format(fielddisplay(self, 'channel_alpha', 'First channel-flow exponent (alpha_s) []')) #TH + s += '{}\n'.format(fielddisplay(self, 'channel_beta', 'Second channel-flow exponent (beta_s) []')) #TH s = '\t--OTHER\n' s += '{}\n'.format(fielddisplay(self, 'spcphi', 'Hydraulic potential Dirichlet constraints [Pa]')) @@ -63,4 +72,5 @@ s += '{}\n'.format(fielddisplay(self, 'requested_outputs', 'additional outputs requested')) s += '{}\n'.format(fielddisplay(self, 'melt_flag', 'User specified basal melt? 0: no (default), 1: use md.basalforcings.groundedice_melting_rate')) + s += '{}\n'.format(fielddisplay(self, 'islaminar','do we use laminar [1] or turbulent physics [0, default]')) return string # }}} @@ -86,4 +96,5 @@ self.pressure_melt_coefficient = 7.5e-8 #K / Pa (See table 1 in Erder et al. 2013) self.cavity_spacing = 2. #m + self.omega = 1./2000.; #TH # Channel parameters @@ -97,4 +108,5 @@ self.requested_outputs = ['default'] self.melt_flag = 0 + self.islaminar = 0 #by default use turbulent physics return self @@ -111,4 +123,7 @@ md = checkfield(md, 'fieldname', 'hydrology.cavity_spacing', 'numel', [1], '>', 0) md = checkfield(md, 'fieldname', 'hydrology.bump_height', 'size', [md.mesh.numberofvertices], '>=', 0, 'np.nan', 1, 'Inf', 1) + md = checkfield(md,'fieldname','hydrology.omega', 'numel', [1], '>=', 0); # TH + md = checkfield(md,'fieldname','hydrology.sheet_alpha', 'numel', [1], '>', 0); # TH + md = checkfield(md,'fieldname','hydrology.sheet_beta', 'numel', [1], '>', 0); # TH # Channels @@ -116,5 +131,6 @@ md = checkfield(md, 'fieldname', 'hydrology.channel_conductivity', 'size', [md.mesh.numberofvertices], '>', 0) md = checkfield(md, 'fieldname', 'hydrology.channel_sheet_width', 'numel', [1], '>=', 0) - md = checkfield(md, 'fieldname', 'hydrology.islaminar', 'numel', [1], 'values', [0, 1]) + md = checkfield(md,'fieldname','hydrology.channel_alpha', 'numel', [1], '>', 0); # TH + md = checkfield(md,'fieldname','hydrology.channel_beta', 'numel', [1], '>', 0); # TH # Other @@ -125,4 +141,5 @@ md = checkfield(md, 'fieldname', 'hydrology.requested_outputs', 'stringrow', 1) md = checkfield(md, 'fieldname', 'hydrology.melt_flag', 'numel', [1], 'values', [0, 1]) + md = checkfield(md, 'fieldname', 'hydrology.islaminar', 'numel', [1], 'values', [0, 1]) if self.melt_flag == 1 or self.melt_flag == 2: md = checkfield(md, 'fieldname', 'basalforcings.groundedice_melting_rate', 'NaN', 1, 'Inf', 1, 'timeseries', 1) @@ -139,4 +156,7 @@ WriteData(fid, prefix, 'object', self, 'class', 'hydrology', 'fieldname', 'cavity_spacing', 'format', 'Double') WriteData(fid, prefix, 'object', self, 'class', 'hydrology', 'fieldname', 'bump_height', 'format', 'DoubleMat', 'mattype', 1) + WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','omega','format','Double'); # TH + WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','sheet_alpha','format','Double'); # TH + WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','sheet_beta','format','Double'); # TH # Channels @@ -144,5 +164,6 @@ WriteData(fid, prefix, 'object', self, 'class', 'hydrology', 'fieldname', 'channel_conductivity', 'format', 'DoubleMat', 'mattype', 1) WriteData(fid, prefix, 'object', self, 'class', 'hydrology', 'fieldname', 'channel_sheet_width', 'format', 'Double') - WriteData(fid, prefix, 'object', self, 'class', 'hydrology', 'fieldname', 'islaminar', 'format', 'Boolean') + WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','channel_alpha','format','Double'); # TH + WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','channel_beta','format','Double'); # TH # Others @@ -152,4 +173,5 @@ WriteData(fid, prefix, 'object', self, 'class', 'hydrology', 'fieldname', 'englacial_void_ratio', 'format', 'Double') WriteData(fid, prefix, 'object', self, 'class', 'hydrology', 'fieldname', 'melt_flag', 'format', 'Integer') + WriteData(fid, prefix, 'object', self, 'class', 'hydrology', 'fieldname', 'islaminar', 'format', 'Integer') outputs = self.requested_outputs