Changeset 27925


Ignore:
Timestamp:
09/27/23 13:18:06 (18 months ago)
Author:
koitf1
Message:

NEW: Implemented laminar flow when islaminar is called in GlaDS hydrology

Location:
issm/trunk-jpl/src
Files:
9 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp

    r27913 r27925  
    146146        iomodel->FetchDataToInput(inputs,elements,"md.hydrology.bump_height",HydrologyBumpHeightEnum);
    147147        iomodel->FetchDataToInput(inputs,elements,"md.hydrology.sheet_conductivity",HydrologySheetConductivityEnum);
     148        iomodel->FetchDataToInput(inputs,elements,"md.hydrology.omega",HydrologyOmegaEnum);
     149        iomodel->FetchDataToInput(inputs,elements,"md.hydrology.sheet_alpha",HydrologySheetAlphaEnum);
     150        iomodel->FetchDataToInput(inputs,elements,"md.hydrology.sheet_beta",HydrologySheetBetaEnum);
    148151        iomodel->FetchDataToInput(inputs,elements,"md.hydrology.channel_conductivity",HydrologyChannelConductivityEnum);
    149152        iomodel->FetchDataToInput(inputs,elements,"md.hydrology.neumannflux",HydrologyNeumannfluxEnum);
    … …  
    183186        parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.melt_flag",HydrologyMeltFlagEnum));
    184187        parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.channel_sheet_width",HydrologyChannelSheetWidthEnum));
     188        parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.channel_alpha",HydrologyChannelAlphaEnum));
     189        parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.channel_beta",HydrologyChannelBetaEnum));
    185190        parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.islaminar",HydrologyIsLaminarEnum));
    186191        parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.englacial_void_ratio",HydrologyEnglacialVoidRatioEnum));
    … …  
    217222        /*Intermediaries */
    218223        IssmDouble  Jdet,dphi[3],h,k;
     224        IssmDouble  alpha,beta,omega,h_r;
    219225        IssmDouble  A,B,n,phi_old,phi,phi_0,H,b,v1;
    220226        IssmDouble* xyz_list = NULL;
    221227
    222228        /*Hard coded coefficients*/
    223         const IssmPDouble alpha = 5./4.;
    224         const IssmPDouble beta  = 3./2.;
     229        const IssmPDouble ALPHA = 5./4.;
     230        const IssmPDouble BETA  = 3./2.;
    225231
    226232        /*Fetch number of nodes and dof for this finite element*/
    … …  
    236242
    237243        /*Get all inputs and parameters*/
     244        int islaminar;
     245        element->FindParam(&islaminar,HydrologyIsLaminarEnum);
    238246        IssmDouble dt        = element->FindParam(TimesteppingTimeStepEnum);
    239247        IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
     248        IssmDouble mu_water  = element->FindParam(MaterialsMuWaterEnum);
    240249        IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
    241250        IssmDouble g         = element->FindParam(ConstantsGEnum);
    242251        IssmDouble e_v       = element->FindParam(HydrologyEnglacialVoidRatioEnum);
     252        Input* hr_input      = element->GetInput(HydrologyBumpHeightEnum);_assert_(hr_input);
    243253        Input* k_input   = element->GetInput(HydrologySheetConductivityEnum);_assert_(k_input);
     254        Input* alpha_input   = element->GetInput(HydrologySheetAlphaEnum);_assert_(alpha_input);
     255        Input* beta_input    = element->GetInput(HydrologySheetBetaEnum);_assert_(beta_input);
     256        Input* omega_input   = element->GetInput(HydrologyOmegaEnum);_assert_(omega_input);
    244257        Input* phi_input = element->GetInput(HydraulicPotentialEnum);      _assert_(phi_input);
    245258        Input* h_input   = element->GetInput(HydrologySheetThicknessEnum); _assert_(h_input);
    … …  
    261274                h_input->GetInputValue(&h,gauss);
    262275                k_input->GetInputValue(&k,gauss);
     276                omega_input->GetInputValue(&omega,gauss);
     277                alpha_input->GetInputValue(&alpha,gauss);
     278                beta_input->GetInputValue(&beta,gauss);
    263279                B_input->GetInputValue(&B,gauss);
    264280                n_input->GetInputValue(&n,gauss);
     281                hr_input->GetInputValue(&h_r,gauss);
    265282                b_input->GetInputValue(&b,gauss);
    266283                H_input->GetInputValue(&H,gauss);
    … …  
    268285                /*Hard code B*/
    269286                B = Cuffey(273.15-2);
    270 
     287               
    271288                /*Get norm of gradient of hydraulic potential and make sure it is >0*/
    272289                IssmDouble normgradphi = sqrt(dphi[0]*dphi[0] + dphi[1]*dphi[1]);
    273290                if(normgradphi < AEPS) normgradphi = AEPS;
    274 
    275                 IssmDouble coeff = k*pow(h,alpha)*pow(normgradphi,beta-2.);
     291               
     292                /*Use Laminar flow if specified*/
     293                IssmDouble nu = mu_water/rho_water;
     294                IssmDouble coeff;
     295                if(islaminar==1 && omega>=AEPS){
     296                        IssmDouble hratio = fabs(h/h_r);
     297                        IssmDouble coarg = 1. + 4.*pow(hratio,3-2*alpha)*omega*k*pow(h,3)*normgradphi/nu;
     298                        coeff = nu/2./omega*pow(hratio,2*alpha-3) * (-1 + pow(coarg, 0.5))/normgradphi;
     299                }
     300                else {
     301                        /*If omega is zero, use standard model, otherwise transition model*/
     302                        coeff = k*pow(h,ALPHA)*pow(normgradphi,BETA-2.);
     303                }
    276304
    277305                /*Diffusive term*/
    … …  
    440468        /*Intermediaries*/
    441469   IssmDouble  dphi[3],h,k,phi;
     470        IssmDouble  alpha,beta,omega,h_r;
    442471        IssmDouble  oceanLS,iceLS;
    443472        IssmDouble* xyz_list = NULL;
    444473
    445474        /*Hard coded coefficients*/
    446         const IssmPDouble alpha = 5./4.;
    447         const IssmPDouble beta  = 3./2.;
     475        const IssmPDouble ALPHA = 5./4.;
     476        const IssmPDouble BETA  = 3./2.;
    448477
    449478        /*Fetch number vertices for this element*/
    … …  
    466495
    467496        /*Retrieve all inputs and parameters*/
     497        int islaminar;
     498        element->FindParam(&islaminar,HydrologyIsLaminarEnum);
    468499        element->GetVerticesCoordinates(&xyz_list);
     500        IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
     501        IssmDouble mu_water  = element->FindParam(MaterialsMuWaterEnum);
    469502        Input *k_input       = element->GetInput(HydrologySheetConductivityEnum); _assert_(k_input);
     503        Input *omega_input   = element->GetInput(HydrologyOmegaEnum);             _assert_(omega_input);
     504        Input *alpha_input   = element->GetInput(HydrologySheetAlphaEnum);        _assert_(alpha_input);
     505        Input *beta_input    = element->GetInput(HydrologySheetBetaEnum);         _assert_(bet_input);
    470506        Input *phi_input     = element->GetInput(HydraulicPotentialEnum);         _assert_(phi_input);
     507        Input *hr_input      = element->GetInput(HydrologyBumpHeightEnum);        _assert_(hr_input);
    471508        Input *h_input       = element->GetInput(HydrologySheetThicknessEnum);    _assert_(h_input);
    472509        Input *oceanLS_input = element->GetInput(MaskOceanLevelsetEnum);          _assert_(oceanLS_input);
    … …  
    482519      phi_input->GetInputValue(&phi,gauss);
    483520      h_input->GetInputValue(&h,gauss);
     521        hr_input->GetInputValue(&h_r,gauss);
    484522      k_input->GetInputValue(&k,gauss);
     523        omega_input->GetInputValue(&omega,gauss);
     524        alpha_input->GetInputValue(&alpha,gauss);
     525        beta_input->GetInputValue(&beta,gauss);
    485526                oceanLS_input->GetInputValue(&oceanLS,gauss);
    486527                iceLS_input->GetInputValue(&iceLS,gauss);
    … …  
    496537         IssmDouble normgradphi = sqrt(dphi[0]*dphi[0] + dphi[1]*dphi[1]);
    497538         if(normgradphi < AEPS) normgradphi = AEPS;
    498 
    499          IssmDouble coeff = k*pow(h,alpha)*pow(normgradphi,beta-2.)/max(AEPS,h); // divide by h to get speed instead of discharge
     539         
     540         /*If omega is zero, use standard model, otherwise transition model*/
     541         IssmDouble nu = mu_water/rho_water;
     542        IssmDouble coeff;
     543        if(islaminar==1 && omega>=AEPS){
     544                IssmDouble hratio = fabs(h/h_r);
     545                IssmDouble coarg = 1. + 4.*pow(hratio,3-2*alpha)*omega*k*pow(h,3)*normgradphi/nu;
     546                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
     547                coeff = k*pow(h,alpha)*pow(normgradphi,beta-2.)/max(AEPS,h);  // divide by h to get speed instead of discharge
     548        }
     549        else {
     550                coeff = k*pow(h,ALPHA)*pow(normgradphi,BETA-2.)/max(AEPS,h);  // divide by h to get speed instead of discharge
     551        }
     552
    500553
    501554                        vx[iv] = -coeff*dphi[0];

  • issm/trunk-jpl/src/c/classes/Loads/Channel.cpp

    r27545 r27925  
    367367        IssmDouble  Jdet,v1,qc,fFactor,Afactor,Bfactor,Xifactor;
    368368        IssmDouble  A,B,n,phi_old,phi,phi_0,dPw,ks,kc,Ngrad;
     369        IssmDouble  alpha_s,beta_s,omega,h_r;
    369370        IssmDouble  H,h,b,dphi[2],dphids,dphimds,db[2],dbds;
    370371        IssmDouble  xyz_list[NUMVERTICES][3];
    … …  
    382383        GetVerticesCoordinates(&xyz_list_tria[0][0],tria->vertices,3);
    383384
     385        int islaminar;
     386        element->FindParam(&islaminar,HydrologyIsLaminarEnum);
    384387        IssmDouble L         = element->FindParam(MaterialsLatentheatEnum);
     388        IssmDouble mu_water  = element->FindParam(MaterialsMuWaterEnum);
    385389        IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
    386390        IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
    … …  
    388392        IssmDouble lc        = element->FindParam(HydrologyChannelSheetWidthEnum);
    389393        IssmDouble c_t       = element->FindParam(HydrologyPressureMeltCoefficientEnum);
     394        IssmDouble alpha_c   = element->FindParam(HydrologyChannelAlphaEnum);
     395        IssmDouble beta_c    = element->FindParam(HydrologyChannelBetaEnum);
    390396
    391397        Input* h_input      = element->GetInput(HydrologySheetThicknessEnum);      _assert_(h_input);
    … …  
    396402        Input* ks_input     = element->GetInput(HydrologySheetConductivityEnum);   _assert_(ks_input);
    397403        Input* kc_input     = element->GetInput(HydrologyChannelConductivityEnum); _assert_(kc_input);
     404        Input* omega_input  = element->GetInput(HydrologyOmegaEnum); _assert_(ks_input);
     405        Input* alpha_s_input= element->GetInput(HydrologySheetAlphaEnum); _assert_(alpha_s_input);
     406        Input* beta_s_input = element->GetInput(HydrologySheetBetaEnum); _assert_(beta_s_input);
     407        Input* hr_input     = element->GetInput(HydrologyBumpHeightEnum);_assert_(hr_input);
    398408        Input* phi_input    = element->GetInput(HydraulicPotentialEnum);           _assert_(phi_input);
    399409
    … …  
    422432                ks_input->GetInputValue(&ks,gauss);
    423433                kc_input->GetInputValue(&kc,gauss);
     434                alpha_s_input->GetInputValue(&alpha_s,gauss);
     435                beta_s_input->GetInputValue(&beta_s,gauss);
     436                omega_input->GetInputValue(&omega,gauss);
     437                hr_input->GetInputValue(&h_r,gauss);
    424438                B_input->GetInputValue(&B,gauss);
    425439                n_input->GetInputValue(&n,gauss);
    … …  
    436450                Ngrad   = fabs(dphids);
    437451                if(Ngrad<AEPS) Ngrad = AEPS;
    438 
    439                 /*Compute the effective conductivity Kc = k h^alpha |grad Phi|^{beta-2} (same for sheet)*/
    440                 IssmDouble Kc = kc * pow(this->S,ALPHA_C) * pow(Ngrad,BETA_C-2.);
    441                 IssmDouble Ks = ks * pow(h      ,ALPHA_S) * pow(Ngrad,BETA_S-2.);
     452               
     453                /*Compute the effective conductivity Kc = k h^alpha |grad Phi|^{beta-2} (same for sheet) and use laminar if specified*/
     454                IssmDouble Kc;
     455                IssmDouble Ks;
     456                IssmDouble nu = mu_water/rho_water;
     457                if(islaminar==1 && omega>=AEPS){
     458                        IssmDouble hratio = h/h_r;
     459                        IssmDouble coarg = 1. + 4.*omega*pow(hratio,3-2*alpha_s)*ks*pow(h,3)*Ngrad/nu;
     460                        Ks = nu/2./omega*pow(hratio,2*alpha_s-3) * (-1 + pow(coarg, 0.5))/Ngrad;
     461                        Kc = kc * pow(this->S,alpha_c) * pow(Ngrad,beta_c-2.);
     462                }
     463                else {
     464                        Ks = ks*pow(h,ALPHA_S)*pow(Ngrad,BETA_S-2.);
     465                        Kc = kc * pow(this->S,ALPHA_C) * pow(Ngrad,BETA_C-2.);
     466                }
    442467
    443468                /*Approx. discharge in the sheet flowing folwing in the direction of the channel ofver a width lc*/
    … …  
    504529        IssmDouble  A,B,n,phi_old,phi,phi_0,dphimds,dphi[2];
    505530        IssmDouble  H,h,b,db[2],dphids,qc,dPw,ks,kc,Ngrad;
     531        IssmDouble  alpha_s,beta_s,omega,h_r;
    506532        IssmDouble  xyz_list[NUMVERTICES][3];
    507533        IssmDouble  xyz_list_tria[3][3];
    … …  
    516542        GetVerticesCoordinates(&xyz_list_tria[0][0],tria->vertices,3);
    517543
     544        int islaminar;
     545        element->FindParam(&islaminar,HydrologyIsLaminarEnum);
    518546        IssmDouble L         = element->FindParam(MaterialsLatentheatEnum);
     547        IssmDouble mu_water  = element->FindParam(MaterialsMuWaterEnum);
    519548        IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
    520549        IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
    … …  
    530559        Input* ks_input     = element->GetInput(HydrologySheetConductivityEnum);   _assert_(ks_input);
    531560        Input* kc_input     = element->GetInput(HydrologyChannelConductivityEnum); _assert_(kc_input);
     561        Input* omega_input  = element->GetInput(HydrologyOmegaEnum); _assert_(omega_input);
    532562        Input* phi_input    = element->GetInput(HydraulicPotentialEnum);           _assert_(phi_input);
     563        Input* alpha_s_input= element->GetInput(HydrologySheetAlphaEnum);_assert_(alpha_s_input);
     564        Input* beta_s_input = element->GetInput(HydrologySheetBetaEnum);_assert_(beta_s_input);
     565        Input* hr_input     = element->GetInput(HydrologyBumpHeightEnum);_assert_(hr_input);
    533566
    534567        /*Get tangent vector*/
    … …  
    552585                ks_input->GetInputValue(&ks,gauss);
    553586                kc_input->GetInputValue(&kc,gauss);
     587                omega_input->GetInputValue(&omega,gauss);
    554588                B_input->GetInputValue(&B,gauss);
    555589                n_input->GetInputValue(&n,gauss);
    … …  
    557591                b_input->GetInputValue(&b,gauss);
    558592                H_input->GetInputValue(&H,gauss);
     593                alpha_s_input->GetInputValue(&alpha_s,gauss);
     594                beta_s_input->GetInputValue(&beta_s,gauss);
     595                hr_input->GetInputValue(&h_r,gauss);
    559596
    560597                /*Hard code B*/
    … …  
    568605                if(Ngrad<AEPS) Ngrad = AEPS;
    569606
    570                 /*Compute the effective conductivity Ks = k h^alpha |grad Phi|^{beta-2} (same for sheet)*/
    571                 IssmDouble Ks = ks * pow(h,ALPHA_S) * pow(Ngrad,BETA_S-2.);
     607
     608                /*Approx. discharge in the sheet flowing folwing in the direction of the channel ofver a width lc, use laminar if specified*/
     609                IssmDouble Ks;
     610                if (islaminar==1 && omega>=AEPS){
     611                IssmDouble hratio = h/h_r;
     612                        IssmDouble nu = mu_water/rho_water;
     613                        IssmDouble coarg = 1. + 4.*omega*pow(hratio,3-2*alpha_s)*ks*pow(h,3)*Ngrad/nu;
     614                        Ks = nu/2./omega*pow(hratio,2*alpha_s-3) * (-1 + pow(coarg, 0.5))/Ngrad;
     615                }
     616                else {
     617                        Ks = ks * pow(h,ALPHA_S) * pow(Ngrad,BETA_S-2.);
     618                }
    572619
    573620                /*Approx. discharge in the sheet flowing folwing in the direction of the channel ofver a width lc*/
    … …  
    636683        /*Intermediaries */
    637684        IssmDouble  A,B,n,phi,phi_0,ks,kc,Ngrad;
     685        IssmDouble  alpha_s,beta_s,omega,h_r;
    638686        IssmDouble  H,h,b,dphi[2],dphids,dphimds,db[2],dbds;
    639687        IssmDouble  xyz_list[NUMVERTICES][3];
    … …  
    644692        GetVerticesCoordinates(&xyz_list_tria[0][0],tria->vertices,3);
    645693
     694        int islaminar;
     695        element->FindParam(&islaminar,HydrologyIsLaminarEnum);
    646696        IssmDouble L         = element->FindParam(MaterialsLatentheatEnum);
    647697        IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
    648698        IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
     699        IssmDouble mu_water  = element->FindParam(MaterialsMuWaterEnum);
    649700        IssmDouble g         = element->FindParam(ConstantsGEnum);
    650701        IssmDouble lc        = element->FindParam(HydrologyChannelSheetWidthEnum);
    651702        IssmDouble c_t       = element->FindParam(HydrologyPressureMeltCoefficientEnum);
    652703        IssmDouble dt        = element->FindParam(TimesteppingTimeStepEnum);
     704        IssmDouble alpha_c   = element->FindParam(HydrologyChannelAlphaEnum);
     705        IssmDouble beta_c    = element->FindParam(HydrologyChannelBetaEnum);
    653706
    654707        Input* h_input      = element->GetInput(HydrologySheetThicknessEnum);      _assert_(h_input);
    … …  
    659712        Input* ks_input     = element->GetInput(HydrologySheetConductivityEnum);   _assert_(ks_input);
    660713        Input* kc_input     = element->GetInput(HydrologyChannelConductivityEnum); _assert_(kc_input);
     714        Input* omega_input  = element->GetInput(HydrologyOmegaEnum); _assert_(omega_input);
    661715        Input* phi_input    = element->GetInput(HydraulicPotentialEnum);           _assert_(phi_input);
     716        Input* alpha_s_input= element->GetInput(HydrologySheetAlphaEnum);_assert_(alpha_s_input);
     717        Input* beta_s_input = element->GetInput(HydrologySheetBetaEnum); _assert_(beta_s_input);
     718        Input* hr_input     = element->GetInput(HydrologyBumpHeightEnum);_assert_(hr_input);
    662719
    663720        /*Get tangent vector*/
    … …  
    674731        ks_input->GetInputValue(&ks,gauss);
    675732        kc_input->GetInputValue(&kc,gauss);
     733        omega_input->GetInputValue(&omega,gauss);
    676734        B_input->GetInputValue(&B,gauss);
    677735        n_input->GetInputValue(&n,gauss);
    … …  
    679737        b_input->GetInputDerivativeValue(&db[0],&xyz_list_tria[0][0],gauss);
    680738        H_input->GetInputValue(&H,gauss);
     739        hr_input->GetInputValue(&h_r,gauss);
     740        alpha_s_input->GetInputValue(&alpha_s,gauss);
     741        beta_s_input->GetInputValue(&beta_s,gauss);
     742
    681743
    682744        /*Hard code B*/
    … …  
    693755        IssmDouble dPw = dphids - dphimds;
    694756
    695         /*Approx. discharge in the sheet flowing folwing in the direction of the channel ofver a width lc*/
    696         IssmDouble qc = - ks * pow(h,ALPHA_S) * pow(Ngrad,BETA_S-2.) * dphids;
     757        /*Approx. discharge in the sheet flowing folwing in the direction of the channel ofver a width lc, use laminar if necessary*/
     758        IssmDouble qc;
     759        if (islaminar==1 && omega>=AEPS){
     760        IssmDouble hratio = h/h_r;
     761                IssmDouble nu = mu_water/rho_water;
     762                IssmDouble coarg = 1. + 4.*omega*pow(hratio,3-2*alpha_s)*ks*pow(h,3)*fabs(Ngrad)/nu;
     763                qc = -nu/2./omega*pow(hratio,2*alpha_s-3) * (-1 + pow(coarg, 0.5))*dphids/Ngrad;
     764        }
     765        else {
     766                qc = - ks * pow(h,ALPHA_S) * pow(Ngrad,BETA_S-2.) * dphids;
     767        }
    697768
    698769        /*Ice rate factor*/

  • issm/trunk-jpl/src/c/shared/Enum/Enum.vim

    r27920 r27925  
    261261syn keyword cConstant HydrologyarmaTimestepEnum
    262262syn keyword cConstant HydrologyAveragingEnum
     263syn keyword cConstant HydrologyChannelAlphaEnum
     264syn keyword cConstant HydrologyChannelBetaEnum
    263265syn keyword cConstant HydrologyCavitySpacingEnum
    264266syn keyword cConstant HydrologyChannelSheetWidthEnum
    … …  
    896898syn keyword cConstant HydrologyMoulinInputEnum
    897899syn keyword cConstant HydrologyNeumannfluxEnum
     900syn keyword cConstant HydrologyOmegaEnum
    898901syn keyword cConstant HydrologyReynoldsEnum
     902syn keyword cConstant HydrologySheetAlphaEnum
     903syn keyword cConstant HydrologySheetBetaEnum
    899904syn keyword cConstant HydrologySheetConductivityEnum
    900905syn keyword cConstant HydrologySheetThicknessEnum

  • issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

    r27920 r27925  
    255255   HydrologyarmaTimestepEnum,
    256256        HydrologyAveragingEnum,
     257        HydrologyChannelAlphaEnum,
     258        HydrologyChannelBetaEnum,
    257259        HydrologyCavitySpacingEnum,     
    258260        HydrologyChannelSheetWidthEnum,
    … …  
    892894        HydrologyMoulinInputEnum,
    893895        HydrologyNeumannfluxEnum,
     896        HydrologyOmegaEnum,
    894897        HydrologyReynoldsEnum,
     898        HydrologySheetAlphaEnum,
     899        HydrologySheetBetaEnum,
    895900        HydrologySheetConductivityEnum,
    896901        HydrologySheetThicknessEnum,

  • issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

    r27920 r27925  
    263263                case HydrologyarmaTimestepEnum : return "HydrologyarmaTimestep";
    264264                case HydrologyAveragingEnum : return "HydrologyAveraging";
     265                case HydrologyChannelAlphaEnum : return "HydrologyChannelAlpha";
     266                case HydrologyChannelBetaEnum : return "HydrologyChannelBeta";
    265267                case HydrologyCavitySpacingEnum : return "HydrologyCavitySpacing";
    266268                case HydrologyChannelSheetWidthEnum : return "HydrologyChannelSheetWidth";
    … …  
    898900                case HydrologyMoulinInputEnum : return "HydrologyMoulinInput";
    899901                case HydrologyNeumannfluxEnum : return "HydrologyNeumannflux";
     902                case HydrologyOmegaEnum : return "HydrologyOmega";
    900903                case HydrologyReynoldsEnum : return "HydrologyReynolds";
     904                case HydrologySheetAlphaEnum : return "HydrologySheetAlpha";
     905                case HydrologySheetBetaEnum : return "HydrologySheetBeta";
    901906                case HydrologySheetConductivityEnum : return "HydrologySheetConductivity";
    902907                case HydrologySheetThicknessEnum : return "HydrologySheetThickness";

  • issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim

    r27920 r27925  
    254254syn keyword juliaConstC HydrologyarmaTimestepEnum
    255255syn keyword juliaConstC HydrologyAveragingEnum
     256syn keyword juliaConstC HydrologyChannelAlphaEnum
     257syn keyword juliaConstC HydrologyChannelBetaEnum
    256258syn keyword juliaConstC HydrologyCavitySpacingEnum
    257259syn keyword juliaConstC HydrologyChannelSheetWidthEnum
    … …  
    889891syn keyword juliaConstC HydrologyMoulinInputEnum
    890892syn keyword juliaConstC HydrologyNeumannfluxEnum
     893syn keyword juliaConstC HydrologyOmegaEnum
    891894syn keyword juliaConstC HydrologyReynoldsEnum
     895syn keyword juliaConstC HydrologySheetAlphaEnum
     896syn keyword juliaConstC HydrologySheetBetaEnum
    892897syn keyword juliaConstC HydrologySheetConductivityEnum
    893898syn keyword juliaConstC HydrologySheetThicknessEnum

  • issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

    r27920 r27925  
    269269              else if (strcmp(name,"HydrologyarmaTimestep")==0) return HydrologyarmaTimestepEnum;
    270270              else if (strcmp(name,"HydrologyAveraging")==0) return HydrologyAveragingEnum;
     271              else if (strcmp(name,"HydrologyChannelAlpha")==0) return HydrologyChannelAlphaEnum;
     272              else if (strcmp(name,"HydrologyChannelBeta")==0) return HydrologyChannelBetaEnum;
    271273              else if (strcmp(name,"HydrologyCavitySpacing")==0) return HydrologyCavitySpacingEnum;
    272274              else if (strcmp(name,"HydrologyChannelSheetWidth")==0) return HydrologyChannelSheetWidthEnum;
    … …  
    381383              else if (strcmp(name,"MasstransportIsfreesurface")==0) return MasstransportIsfreesurfaceEnum;
    382384              else if (strcmp(name,"MasstransportMinThickness")==0) return MasstransportMinThicknessEnum;
    383               else if (strcmp(name,"MasstransportNumRequestedOutputs")==0) return MasstransportNumRequestedOutputsEnum;
    384               else if (strcmp(name,"MasstransportPenaltyFactor")==0) return MasstransportPenaltyFactorEnum;
    385385         else stage=4;
    386386   }
    387387   if(stage==4){
    388               if (strcmp(name,"MasstransportRequestedOutputs")==0) return MasstransportRequestedOutputsEnum;
     388              if (strcmp(name,"MasstransportNumRequestedOutputs")==0) return MasstransportNumRequestedOutputsEnum;
     389              else if (strcmp(name,"MasstransportPenaltyFactor")==0) return MasstransportPenaltyFactorEnum;
     390              else if (strcmp(name,"MasstransportRequestedOutputs")==0) return MasstransportRequestedOutputsEnum;
    389391              else if (strcmp(name,"MasstransportStabilization")==0) return MasstransportStabilizationEnum;
    390392              else if (strcmp(name,"MaterialsBeta")==0) return MaterialsBetaEnum;
    … …  
    504506              else if (strcmp(name,"SolidearthSettingsSealevelLoading")==0) return SolidearthSettingsSealevelLoadingEnum;
    505507              else if (strcmp(name,"SolidearthSettingsGRD")==0) return SolidearthSettingsGRDEnum;
    506               else if (strcmp(name,"SolidearthSettingsRunFrequency")==0) return SolidearthSettingsRunFrequencyEnum;
    507               else if (strcmp(name,"SolidearthSettingsTimeAcc")==0) return SolidearthSettingsTimeAccEnum;
    508508         else stage=5;
    509509   }
    510510   if(stage==5){
    511               if (strcmp(name,"SolidearthSettingsHoriz")==0) return SolidearthSettingsHorizEnum;
     511              if (strcmp(name,"SolidearthSettingsRunFrequency")==0) return SolidearthSettingsRunFrequencyEnum;
     512              else if (strcmp(name,"SolidearthSettingsTimeAcc")==0) return SolidearthSettingsTimeAccEnum;
     513              else if (strcmp(name,"SolidearthSettingsHoriz")==0) return SolidearthSettingsHorizEnum;
    512514              else if (strcmp(name,"SolidearthSettingsMaxiter")==0) return SolidearthSettingsMaxiterEnum;
    513515              else if (strcmp(name,"SolidearthSettingsGrdOcean")==0) return SolidearthSettingsGrdOceanEnum;
    … …  
    627629              else if (strcmp(name,"SmbStepsPerStep")==0) return SmbStepsPerStepEnum;
    628630              else if (strcmp(name,"SmbSwIdx")==0) return SmbSwIdxEnum;
    629               else if (strcmp(name,"SmbSWgrad")==0) return SmbSWgradEnum;
    630               else if (strcmp(name,"SmbT0dry")==0) return SmbT0dryEnum;
    631631         else stage=6;
    632632   }
    633633   if(stage==6){
    634               if (strcmp(name,"SmbT0wet")==0) return SmbT0wetEnum;
     634              if (strcmp(name,"SmbSWgrad")==0) return SmbSWgradEnum;
     635              else if (strcmp(name,"SmbT0dry")==0) return SmbT0dryEnum;
     636              else if (strcmp(name,"SmbT0wet")==0) return SmbT0wetEnum;
    635637              else if (strcmp(name,"SmbTcIdx")==0) return SmbTcIdxEnum;
    636638              else if (strcmp(name,"SmbTeThresh")==0) return SmbTeThreshEnum;
    … …  
    750752              else if (strcmp(name,"BasalforcingsSpatialDeepwaterElevation")==0) return BasalforcingsSpatialDeepwaterElevationEnum;
    751753              else if (strcmp(name,"BasalforcingsSpatialDeepwaterMeltingRate")==0) return BasalforcingsSpatialDeepwaterMeltingRateEnum;
    752               else if (strcmp(name,"BasalforcingsSpatialUpperwaterElevation")==0) return BasalforcingsSpatialUpperwaterElevationEnum;
    753               else if (strcmp(name,"BasalforcingsSpatialUpperwaterMeltingRate")==0) return BasalforcingsSpatialUpperwaterMeltingRateEnum;
    754754         else stage=7;
    755755   }
    756756   if(stage==7){
    757               if (strcmp(name,"BasalforcingsIsmip6BasinId")==0) return BasalforcingsIsmip6BasinIdEnum;
     757              if (strcmp(name,"BasalforcingsSpatialUpperwaterElevation")==0) return BasalforcingsSpatialUpperwaterElevationEnum;
     758              else if (strcmp(name,"BasalforcingsSpatialUpperwaterMeltingRate")==0) return BasalforcingsSpatialUpperwaterMeltingRateEnum;
     759              else if (strcmp(name,"BasalforcingsIsmip6BasinId")==0) return BasalforcingsIsmip6BasinIdEnum;
    758760              else if (strcmp(name,"BasalforcingsIsmip6Tf")==0) return BasalforcingsIsmip6TfEnum;
    759761              else if (strcmp(name,"BasalforcingsIsmip6TfShelf")==0) return BasalforcingsIsmip6TfShelfEnum;
    … …  
    873875              else if (strcmp(name,"FrictionQ")==0) return FrictionQEnum;
    874876              else if (strcmp(name,"FrictionSedimentCompressibilityCoefficient")==0) return FrictionSedimentCompressibilityCoefficientEnum;
    875               else if (strcmp(name,"FrictionTillFrictionAngle")==0) return FrictionTillFrictionAngleEnum;
    876               else if (strcmp(name,"FrictionWaterLayer")==0) return FrictionWaterLayerEnum;
    877877         else stage=8;
    878878   }
    879879   if(stage==8){
    880               if (strcmp(name,"FrictionWaterPressure")==0) return FrictionWaterPressureEnum;
     880              if (strcmp(name,"FrictionTillFrictionAngle")==0) return FrictionTillFrictionAngleEnum;
     881              else if (strcmp(name,"FrictionWaterLayer")==0) return FrictionWaterLayerEnum;
     882              else if (strcmp(name,"FrictionWaterPressure")==0) return FrictionWaterPressureEnum;
    881883              else if (strcmp(name,"FrictionWaterPressureNoise")==0) return FrictionWaterPressureNoiseEnum;
    882884              else if (strcmp(name,"Frictionf")==0) return FrictionfEnum;
    … …  
    919921              else if (strcmp(name,"HydrologyMoulinInput")==0) return HydrologyMoulinInputEnum;
    920922              else if (strcmp(name,"HydrologyNeumannflux")==0) return HydrologyNeumannfluxEnum;
     923              else if (strcmp(name,"HydrologyOmega")==0) return HydrologyOmegaEnum;
    921924              else if (strcmp(name,"HydrologyReynolds")==0) return HydrologyReynoldsEnum;
     925              else if (strcmp(name,"HydrologySheetAlpha")==0) return HydrologySheetAlphaEnum;
     926              else if (strcmp(name,"HydrologySheetBeta")==0) return HydrologySheetBetaEnum;
    922927              else if (strcmp(name,"HydrologySheetConductivity")==0) return HydrologySheetConductivityEnum;
    923928              else if (strcmp(name,"HydrologySheetThickness")==0) return HydrologySheetThicknessEnum;
    … …  
    993998              else if (strcmp(name,"SamplingKappa")==0) return SamplingKappaEnum;
    994999              else if (strcmp(name,"SamplingPhi")==0) return SamplingPhiEnum;
    995               else if (strcmp(name,"SamplingTau")==0) return SamplingTauEnum;
     1000         else stage=9;
     1001   }
     1002   if(stage==9){
     1003              if (strcmp(name,"SamplingTau")==0) return SamplingTauEnum;
    9961004              else if (strcmp(name,"Sealevel")==0) return SealevelEnum;
    9971005              else if (strcmp(name,"SealevelGRD")==0) return SealevelGRDEnum;
    9981006              else if (strcmp(name,"SatGraviGRD")==0) return SatGraviGRDEnum;
    9991007              else if (strcmp(name,"SealevelBarystaticMask")==0) return SealevelBarystaticMaskEnum;
    1000          else stage=9;
    1001    }
    1002    if(stage==9){
    1003               if (strcmp(name,"SealevelBarystaticIceMask")==0) return SealevelBarystaticIceMaskEnum;
     1008              else if (strcmp(name,"SealevelBarystaticIceMask")==0) return SealevelBarystaticIceMaskEnum;
    10041009              else if (strcmp(name,"SealevelBarystaticIceWeights")==0) return SealevelBarystaticIceWeightsEnum;
    10051010              else if (strcmp(name,"SealevelBarystaticIceArea")==0) return SealevelBarystaticIceAreaEnum;
    … …  
    11161121              else if (strcmp(name,"SmbEAir")==0) return SmbEAirEnum;
    11171122              else if (strcmp(name,"SmbEC")==0) return SmbECEnum;
    1118               else if (strcmp(name,"SmbECDt")==0) return SmbECDtEnum;
     1123         else stage=10;
     1124   }
     1125   if(stage==10){
     1126              if (strcmp(name,"SmbECDt")==0) return SmbECDtEnum;
    11191127              else if (strcmp(name,"SmbECini")==0) return SmbECiniEnum;
    11201128              else if (strcmp(name,"SmbEla")==0) return SmbElaEnum;
    11211129              else if (strcmp(name,"SmbEvaporation")==0) return SmbEvaporationEnum;
    11221130              else if (strcmp(name,"SmbFAC")==0) return SmbFACEnum;
    1123          else stage=10;
    1124    }
    1125    if(stage==10){
    1126               if (strcmp(name,"SmbGdn")==0) return SmbGdnEnum;
     1131              else if (strcmp(name,"SmbGdn")==0) return SmbGdnEnum;
    11271132              else if (strcmp(name,"SmbGdnini")==0) return SmbGdniniEnum;
    11281133              else if (strcmp(name,"SmbGsp")==0) return SmbGspEnum;
    … …  
    12391244              else if (strcmp(name,"SurfaceArea")==0) return SurfaceAreaEnum;
    12401245              else if (strcmp(name,"SurfaceAverageVelMisfit")==0) return SurfaceAverageVelMisfitEnum;
    1241               else if (strcmp(name,"SurfaceCrevasse")==0) return SurfaceCrevasseEnum;
     1246         else stage=11;
     1247   }
     1248   if(stage==11){
     1249              if (strcmp(name,"SurfaceCrevasse")==0) return SurfaceCrevasseEnum;
    12421250              else if (strcmp(name,"Surface")==0) return SurfaceEnum;
    12431251              else if (strcmp(name,"SurfaceOld")==0) return SurfaceOldEnum;
    12441252              else if (strcmp(name,"SurfaceLogVelMisfit")==0) return SurfaceLogVelMisfitEnum;
    12451253              else if (strcmp(name,"SurfaceLogVxVyMisfit")==0) return SurfaceLogVxVyMisfitEnum;
    1246          else stage=11;
    1247    }
    1248    if(stage==11){
    1249               if (strcmp(name,"SurfaceObservation")==0) return SurfaceObservationEnum;
     1254              else if (strcmp(name,"SurfaceObservation")==0) return SurfaceObservationEnum;
    12501255              else if (strcmp(name,"SurfaceRelVelMisfit")==0) return SurfaceRelVelMisfitEnum;
    12511256              else if (strcmp(name,"SurfaceSlopeX")==0) return SurfaceSlopeXEnum;
    … …  
    13621367              else if (strcmp(name,"Outputdefinition58")==0) return Outputdefinition58Enum;
    13631368              else if (strcmp(name,"Outputdefinition59")==0) return Outputdefinition59Enum;
    1364               else if (strcmp(name,"Outputdefinition5")==0) return Outputdefinition5Enum;
     1369         else stage=12;
     1370   }
     1371   if(stage==12){
     1372              if (strcmp(name,"Outputdefinition5")==0) return Outputdefinition5Enum;
    13651373              else if (strcmp(name,"Outputdefinition60")==0) return Outputdefinition60Enum;
    13661374              else if (strcmp(name,"Outputdefinition61")==0) return Outputdefinition61Enum;
    13671375              else if (strcmp(name,"Outputdefinition62")==0) return Outputdefinition62Enum;
    13681376              else if (strcmp(name,"Outputdefinition63")==0) return Outputdefinition63Enum;
    1369          else stage=12;
    1370    }
    1371    if(stage==12){
    1372               if (strcmp(name,"Outputdefinition64")==0) return Outputdefinition64Enum;
     1377              else if (strcmp(name,"Outputdefinition64")==0) return Outputdefinition64Enum;
    13731378              else if (strcmp(name,"Outputdefinition65")==0) return Outputdefinition65Enum;
    13741379              else if (strcmp(name,"Outputdefinition66")==0) return Outputdefinition66Enum;
    … …  
    14851490              else if (strcmp(name,"DataSetParam")==0) return DataSetParamEnum;
    14861491              else if (strcmp(name,"DatasetInput")==0) return DatasetInputEnum;
    1487               else if (strcmp(name,"DebrisAnalysis")==0) return DebrisAnalysisEnum;
     1492         else stage=13;
     1493   }
     1494   if(stage==13){
     1495              if (strcmp(name,"DebrisAnalysis")==0) return DebrisAnalysisEnum;
    14881496              else if (strcmp(name,"DebrisSolution")==0) return DebrisSolutionEnum;
    14891497              else if (strcmp(name,"DefaultAnalysis")==0) return DefaultAnalysisEnum;
    … …  
    14911499              else if (strcmp(name,"Dense")==0) return DenseEnum;
    14921500              else if (strcmp(name,"DependentObject")==0) return DependentObjectEnum;
    1493          else stage=13;
    1494    }
    1495    if(stage==13){
    1496               if (strcmp(name,"DepthAverageAnalysis")==0) return DepthAverageAnalysisEnum;
     1501              else if (strcmp(name,"DepthAverageAnalysis")==0) return DepthAverageAnalysisEnum;
    14971502              else if (strcmp(name,"DeviatoricStressErrorEstimator")==0) return DeviatoricStressErrorEstimatorEnum;
    14981503              else if (strcmp(name,"Divergence")==0) return DivergenceEnum;
    … …  
    16081613              else if (strcmp(name,"Loads")==0) return LoadsEnum;
    16091614              else if (strcmp(name,"LoveAnalysis")==0) return LoveAnalysisEnum;
    1610               else if (strcmp(name,"LoveHf")==0) return LoveHfEnum;
     1615         else stage=14;
     1616   }
     1617   if(stage==14){
     1618              if (strcmp(name,"LoveHf")==0) return LoveHfEnum;
    16111619              else if (strcmp(name,"LoveHt")==0) return LoveHtEnum;
    16121620              else if (strcmp(name,"LoveKernelsImag")==0) return LoveKernelsImagEnum;
    … …  
    16141622              else if (strcmp(name,"LoveKf")==0) return LoveKfEnum;
    16151623              else if (strcmp(name,"LoveKt")==0) return LoveKtEnum;
    1616          else stage=14;
    1617    }
    1618    if(stage==14){
    1619               if (strcmp(name,"LoveLf")==0) return LoveLfEnum;
     1624              else if (strcmp(name,"LoveLf")==0) return LoveLfEnum;
    16201625              else if (strcmp(name,"LoveLt")==0) return LoveLtEnum;
    16211626              else if (strcmp(name,"LoveTidalHt")==0) return LoveTidalHtEnum;
    … …  
    17311736              else if (strcmp(name,"SMBpdd")==0) return SMBpddEnum;
    17321737              else if (strcmp(name,"SMBpddSicopolis")==0) return SMBpddSicopolisEnum;
    1733               else if (strcmp(name,"SMBsemic")==0) return SMBsemicEnum;
     1738         else stage=15;
     1739   }
     1740   if(stage==15){
     1741              if (strcmp(name,"SMBsemic")==0) return SMBsemicEnum;
    17341742              else if (strcmp(name,"SSAApproximation")==0) return SSAApproximationEnum;
    17351743              else if (strcmp(name,"SSAFSApproximation")==0) return SSAFSApproximationEnum;
    … …  
    17371745              else if (strcmp(name,"Scaled")==0) return ScaledEnum;
    17381746              else if (strcmp(name,"SealevelAbsolute")==0) return SealevelAbsoluteEnum;
    1739          else stage=15;
    1740    }
    1741    if(stage==15){
    1742               if (strcmp(name,"SealevelEmotion")==0) return SealevelEmotionEnum;
     1747              else if (strcmp(name,"SealevelEmotion")==0) return SealevelEmotionEnum;
    17431748              else if (strcmp(name,"SealevelchangePolarMotionX")==0) return SealevelchangePolarMotionXEnum;
    17441749              else if (strcmp(name,"SealevelchangePolarMotionY")==0) return SealevelchangePolarMotionYEnum;

  • issm/trunk-jpl/src/m/classes/hydrologyglads.m

    r27913 r27925  
    1111                cavity_spacing            = 0.;
    1212                bump_height               = NaN;
     13                omega                     = 0; % TH
     14                sheet_alpha               = 5.0/4.0; % TH
     15                sheet_beta                = 3.0/2.0; % TH
    1316
    1417                %Channels
    … …  
    1619                channel_conductivity = NaN;
    1720                channel_sheet_width  = 0.;
    18                 islaminar            = 0;
     21                channel_alpha        = 5.0/4.0; % TH
     22                channel_beta         = 3.0/2.0; % TH
    1923
    2024                %Other
    … …  
    2529                requested_outputs    = {};
    2630                melt_flag            = 0;
     31                islaminar            = 0;
    2732        end
    2833        methods
    … …  
    4651                        self.pressure_melt_coefficient = 7.5e-8; %K/Pa (See table 1 in Erder et al. 2013)
    4752                        self.cavity_spacing = 2.; %m
     53                        self.omega = 1./2000.; % TH
    4854
    4955                        %Channel parameters
    … …  
    5157                        self.channel_conductivity = 5.e-2; %Dow's default, Table uses 0.1
    5258                        self.channel_sheet_width = 2.; %m
    53                         self.islaminar = 0; %by default use GlaDS default turbulent code
    5459
    55                         %Other
     60                        %Otherself.omega = 1./2000.; % TH
    5661                        self.englacial_void_ratio = 1.e-5;% Dow's default, Table from Werder et al. uses 1e-3;
    5762                        self.requested_outputs={'default'};
    5863                        self.melt_flag=0;
     64                        self.islaminar = 0; %by default use GlaDS default turbulent code
    5965                end % }}}
    6066                function md = checkconsistency(self,md,solution,analyses) % {{{
    … …  
    7076                        md = checkfield(md,'fieldname','hydrology.cavity_spacing','numel',[1],'>',0);
    7177                        md = checkfield(md,'fieldname','hydrology.bump_height','size',[md.mesh.numberofvertices 1],'>=',0,'NaN',1,'Inf',1);
     78                        md = checkfield(md,'fieldname','hydrology.omega', 'numel', [1], '>=', 0); % TH
     79                        md = checkfield(md,'fieldname','hydrology.sheet_alpha', 'numel', [1], '>', 0); % TH
     80                        md = checkfield(md,'fieldname','hydrology.sheet_beta', 'numel', [1], '>', 0); % TH
    7281
    7382                        %Channels
    … …  
    7584                        md = checkfield(md,'fieldname','hydrology.channel_conductivity','size',[md.mesh.numberofvertices 1],'>=',0,'NaN',1,'Inf',1);
    7685                        md = checkfield(md,'fieldname','hydrology.channel_sheet_width','numel',[1],'>=',0);
    77                         md = checkfield(md,'fieldname','hydrology.islaminar','numel',[1],'values',[0 1]);
     86                        md = checkfield(md,'fieldname','hydrology.channel_alpha', 'numel', [1], '>', 0); % TH
     87                        md = checkfield(md,'fieldname','hydrology.channel_beta', 'numel', [1], '>', 0); % TH
    7888
    7989                        %Other
    … …  
    8494                        md = checkfield(md,'fieldname','hydrology.requested_outputs','stringrow',1);
    8595                        md = checkfield(md,'fieldname','hydrology.melt_flag','numel',[1],'values',[0 1 2]);
     96                        md = checkfield(md,'fieldname','hydrology.islaminar','numel',[1],'values',[0 1]);
    8697                        if self.melt_flag==1 || self.melt_flag==2
    8798                                md = checkfield(md,'fieldname','basalforcings.groundedice_melting_rate','NaN',1,'Inf',1,'timeseries',1);
    … …  
    93104                        fielddisplay(self,'pressure_melt_coefficient','Pressure melt coefficient (c_t) [K Pa^-1]');
    94105                        fielddisplay(self,'sheet_conductivity','sheet conductivity (k) [m^(7/4) kg^(-1/2)]');
     106                        fielddisplay(self,'sheet_alpha','First sheet-flow exponent (alpha_s) []'); % TH
     107                        fielddisplay(self,'sheet_beta','Second sheet-flow exponent (beta_s) []'); % TH
    95108                        fielddisplay(self,'cavity_spacing','cavity spacing (l_r) [m]');
    96109                        fielddisplay(self,'bump_height','typical bump height (h_r) [m]');
     110                        fielddisplay(self,'omega','transition parameter (omega) []'); % TH
    97111                        disp(sprintf('      CHANNELS'));
    98112                        fielddisplay(self,'ischannels','Do we allow for channels? 1: yes, 0: no');
    99113                        fielddisplay(self,'channel_conductivity','channel conductivity (k_c) [m^(3/2) kg^(-1/2)]');
     114                        fielddisplay(self,'channel_alpha','First channel-flow exponent (alpha_s) []'); % TH
     115                        fielddisplay(self,'channel_beta','Second channel-flow exponent (beta_s) []'); % TH
    100116                        fielddisplay(self,'channel_sheet_width','channel sheet width [m]');
    101                         fielddisplay(self,'islaminar','do we use laminar [1] or turbulent physics [0, default]');
    102117                        disp(sprintf('      OTHER'));
    103118                        fielddisplay(self,'spcphi','Hydraulic potential Dirichlet constraints [Pa]');
    … …  
    107122                        fielddisplay(self,'requested_outputs','additional outputs requested');
    108123                        fielddisplay(self,'melt_flag','User specified basal melt? 0: no (default), 1: use md.basalforcings.groundedice_melting_rate');
     124                        fielddisplay(self,'islaminar','do we use laminar [1] or turbulent physics [0, default]');
    109125                end % }}}
    110126                function marshall(self,prefix,md,fid) % {{{
    … …  
    120136                        WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','cavity_spacing','format','Double');
    121137                        WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','bump_height','format','DoubleMat','mattype',1);
     138                        WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','omega','format','Double'); % TH
     139                        WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','sheet_alpha','format','Double'); % TH
     140                        WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','sheet_beta','format','Double'); % TH
    122141
    123142                        %Channels
    … …  
    125144                        WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','channel_conductivity','format','DoubleMat','mattype',1);
    126145                        WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','channel_sheet_width','format','Double');
    127                         WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','islaminar','format','Boolean');
     146                        WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','channel_alpha','format','Double'); % TH
     147                        WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','channel_beta','format','Double'); % TH
    128148
    129149                        %Others
    … …  
    133153                        WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','englacial_void_ratio','format','Double');
    134154                        WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','melt_flag','format','Integer');
     155                        WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','islaminar','format','Integer');
    135156                        outputs = self.requested_outputs;
    136157                        pos  = find(ismember(outputs,'default'));

  • issm/trunk-jpl/src/m/classes/hydrologyglads.py

    r27913 r27925  
    1919        self.cavity_spacing = 0.
    2020        self.bump_height = np.nan
     21        self.omega = 0.; #TH
     22        self.sheet_alpha = 5.0/4.0; # TH
     23        self.sheet_beta = 3.0/2.0; # TH
    2124
    2225        # Channels
    … …  
    2427        self.channel_conductivity = np.nan
    2528        self.channel_sheet_width = 0.
    26         self.islaminar = 0
     29        self.channel_alpha = 5.0/4.0; # TH
     30        self.channel_beta = 3.0/2.0; # TH
    2731
    2832        # Other
    … …  
    3337        self.requested_outputs = []
    3438        self.melt_flag = 0
     39        self.islaminar = 0
    3540
    3641        nargs = len(args)
    … …  
    4954        s += '{}\n'.format(fielddisplay(self, 'pressure_melt_coefficient', 'Pressure melt coefficient (c_t) [K Pa^ - 1]'))
    5055        s += '{}\n'.format(fielddisplay(self, 'sheet_conductivity', 'sheet conductivity (k) [m^(7 / 4) kg^(- 1 / 2)]'))
     56        s += '{}\n'.format(fielddisplay(self, 'sheet_alpha', 'First sheet-flow exponent (alpha_s) []')) #TH
     57        s += '{}\n'.format(fielddisplay(self, 'sheet_beta', 'Second sheet-flow exponent (beta_s) []')) #TH
    5158        s += '{}\n'.format(fielddisplay(self, 'cavity_spacing', 'cavity spacing (l_r) [m]'))
    5259        s += '{}\n'.format(fielddisplay(self, 'bump_height', 'typical bump height (h_r) [m]'))
     60        s += '{}\n'.format(fielddisplay(self, 'omega', 'transition parameter (omega) []')) #TH
    5361        s = '\t--CHANNELS\n'
    5462        s += '{}\n'.format(fielddisplay(self, 'ischannels', 'Do we allow for channels? 1: yes, 0: no'))
    5563        s += '{}\n'.format(fielddisplay(self, 'channel_conductivity', 'channel conductivity (k_c) [m^(3 / 2) kg^(- 1 / 2)]'))
    5664        s += '{}\n'.format(fielddisplay(self, 'channel_sheet_width', 'channel sheet width [m]'))
    57         s += '{}\n'.format(fielddisplay(self, 'islaminar','do we use laminar [1] or turbulent physics [0, default]'))
     65        s += '{}\n'.format(fielddisplay(self, 'channel_alpha', 'First channel-flow exponent (alpha_s) []')) #TH
     66        s += '{}\n'.format(fielddisplay(self, 'channel_beta', 'Second channel-flow exponent (beta_s) []')) #TH
    5867        s = '\t--OTHER\n'
    5968        s += '{}\n'.format(fielddisplay(self, 'spcphi', 'Hydraulic potential Dirichlet constraints [Pa]'))
    … …  
    6372        s += '{}\n'.format(fielddisplay(self, 'requested_outputs', 'additional outputs requested'))
    6473        s += '{}\n'.format(fielddisplay(self, 'melt_flag', 'User specified basal melt? 0: no (default), 1: use md.basalforcings.groundedice_melting_rate'))
     74        s += '{}\n'.format(fielddisplay(self, 'islaminar','do we use laminar [1] or turbulent physics [0, default]'))
    6575        return string
    6676    # }}}
    … …  
    8696        self.pressure_melt_coefficient = 7.5e-8  #K / Pa (See table 1 in Erder et al. 2013)
    8797        self.cavity_spacing = 2.  #m
     98        self.omega = 1./2000.; #TH
    8899
    89100        # Channel parameters
    … …  
    97108        self.requested_outputs = ['default']
    98109        self.melt_flag = 0
     110        self.islaminar = 0  #by default use turbulent physics
    99111
    100112        return self
    … …  
    111123        md = checkfield(md, 'fieldname', 'hydrology.cavity_spacing', 'numel', [1], '>', 0)
    112124        md = checkfield(md, 'fieldname', 'hydrology.bump_height', 'size', [md.mesh.numberofvertices], '>=', 0, 'np.nan', 1, 'Inf', 1)
     125        md = checkfield(md,'fieldname','hydrology.omega', 'numel', [1], '>=', 0); # TH
     126        md = checkfield(md,'fieldname','hydrology.sheet_alpha', 'numel', [1], '>', 0); # TH
     127        md = checkfield(md,'fieldname','hydrology.sheet_beta', 'numel', [1], '>', 0); # TH
    113128
    114129        # Channels
    … …  
    116131        md = checkfield(md, 'fieldname', 'hydrology.channel_conductivity', 'size', [md.mesh.numberofvertices], '>', 0)
    117132        md = checkfield(md, 'fieldname', 'hydrology.channel_sheet_width', 'numel', [1], '>=', 0)
    118         md = checkfield(md, 'fieldname', 'hydrology.islaminar', 'numel', [1], 'values', [0, 1])
     133        md = checkfield(md,'fieldname','hydrology.channel_alpha', 'numel', [1], '>', 0); # TH
     134        md = checkfield(md,'fieldname','hydrology.channel_beta', 'numel', [1], '>', 0); # TH
    119135
    120136        # Other
    … …  
    125141        md = checkfield(md, 'fieldname', 'hydrology.requested_outputs', 'stringrow', 1)
    126142        md = checkfield(md, 'fieldname', 'hydrology.melt_flag', 'numel', [1], 'values', [0, 1])
     143        md = checkfield(md, 'fieldname', 'hydrology.islaminar', 'numel', [1], 'values', [0, 1])
    127144        if self.melt_flag == 1 or self.melt_flag == 2:
    128145            md = checkfield(md, 'fieldname', 'basalforcings.groundedice_melting_rate', 'NaN', 1, 'Inf', 1, 'timeseries', 1)
    … …  
    139156        WriteData(fid, prefix, 'object', self, 'class', 'hydrology', 'fieldname', 'cavity_spacing', 'format', 'Double')
    140157        WriteData(fid, prefix, 'object', self, 'class', 'hydrology', 'fieldname', 'bump_height', 'format', 'DoubleMat', 'mattype', 1)
     158        WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','omega','format','Double'); # TH
     159        WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','sheet_alpha','format','Double'); # TH
     160        WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','sheet_beta','format','Double'); # TH
    141161
    142162        # Channels
    … …  
    144164        WriteData(fid, prefix, 'object', self, 'class', 'hydrology', 'fieldname', 'channel_conductivity', 'format', 'DoubleMat', 'mattype', 1)
    145165        WriteData(fid, prefix, 'object', self, 'class', 'hydrology', 'fieldname', 'channel_sheet_width', 'format', 'Double')
    146         WriteData(fid, prefix, 'object', self, 'class', 'hydrology', 'fieldname', 'islaminar', 'format', 'Boolean')
     166        WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','channel_alpha','format','Double'); # TH
     167        WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','channel_beta','format','Double'); # TH
    147168
    148169        # Others
    … …  
    152173        WriteData(fid, prefix, 'object', self, 'class', 'hydrology', 'fieldname', 'englacial_void_ratio', 'format', 'Double')
    153174        WriteData(fid, prefix, 'object', self, 'class', 'hydrology', 'fieldname', 'melt_flag', 'format', 'Integer')
     175        WriteData(fid, prefix, 'object', self, 'class', 'hydrology', 'fieldname', 'islaminar', 'format', 'Integer')
    154176
    155177        outputs = self.requested_outputs
Note: See TracChangeset for help on using the changeset viewer.