Changeset 27936
- Timestamp:
- 10/03/23 09:26:58 (19 months ago)
- Location:
- issm/trunk-jpl/src
- Files:
-
- 9 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp
r27927 r27936 146 146 iomodel->FetchDataToInput(inputs,elements,"md.hydrology.bump_height",HydrologyBumpHeightEnum); 147 147 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);151 148 iomodel->FetchDataToInput(inputs,elements,"md.hydrology.channel_conductivity",HydrologyChannelConductivityEnum); 152 149 iomodel->FetchDataToInput(inputs,elements,"md.hydrology.neumannflux",HydrologyNeumannfluxEnum); … … 188 185 parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.channel_alpha",HydrologyChannelAlphaEnum)); 189 186 parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.channel_beta",HydrologyChannelBetaEnum)); 190 parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.islaminar",HydrologyIsLaminarEnum)); 187 parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.sheet_alpha",HydrologySheetAlphaEnum)); 188 parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.sheet_beta",HydrologySheetBetaEnum)); 189 parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.omega",HydrologyOmegaEnum)); 190 parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.istransition",HydrologyIsTransitionEnum)); 191 191 parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.englacial_void_ratio",HydrologyEnglacialVoidRatioEnum)); 192 192 … … 222 222 /*Intermediaries */ 223 223 IssmDouble Jdet,dphi[3],h,k; 224 IssmDouble alpha,beta,omega,h_r;224 IssmDouble h_r; 225 225 IssmDouble A,B,n,phi_old,phi,phi_0,H,b,v1; 226 226 IssmDouble* xyz_list = NULL; 227 228 /*Hard coded coefficients*/229 const IssmPDouble ALPHA = 5./4.;230 const IssmPDouble BETA = 3./2.;231 227 232 228 /*Fetch number of nodes and dof for this finite element*/ … … 242 238 243 239 /*Get all inputs and parameters*/ 244 bool islaminar; 245 element->FindParam(&islaminar,HydrologyIsLaminarEnum); 240 bool istransition; 241 element->FindParam(&istransition,HydrologyIsTransitionEnum); 242 IssmDouble alpha = element->FindParam(HydrologySheetAlphaEnum); 243 IssmDouble beta = element->FindParam(HydrologySheetBetaEnum); 244 IssmDouble omega = element->FindParam(HydrologyOmegaEnum); 246 245 IssmDouble dt = element->FindParam(TimesteppingTimeStepEnum); 247 246 IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum); … … 252 251 Input* hr_input = element->GetInput(HydrologyBumpHeightEnum);_assert_(hr_input); 253 252 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);257 253 Input* phi_input = element->GetInput(HydraulicPotentialEnum); _assert_(phi_input); 258 254 Input* h_input = element->GetInput(HydrologySheetThicknessEnum); _assert_(h_input); … … 274 270 h_input->GetInputValue(&h,gauss); 275 271 k_input->GetInputValue(&k,gauss); 276 omega_input->GetInputValue(&omega,gauss);277 alpha_input->GetInputValue(&alpha,gauss);278 beta_input->GetInputValue(&beta,gauss);279 272 B_input->GetInputValue(&B,gauss); 280 273 n_input->GetInputValue(&n,gauss); … … 290 283 if(normgradphi < AEPS) normgradphi = AEPS; 291 284 292 /*Use Laminar flowif specified*/285 /*Use transition model if specified*/ 293 286 IssmDouble nu = mu_water/rho_water; 294 287 IssmDouble coeff; 295 if(is laminar==1 && omega>=AEPS){288 if(istransition==1 && omega>=AEPS){ 296 289 IssmDouble hratio = fabs(h/h_r); 297 290 IssmDouble coarg = 1. + 4.*pow(hratio,3-2*alpha)*omega*k*pow(h,3)*normgradphi/nu; … … 300 293 else { 301 294 /*If omega is zero, use standard model, otherwise transition model*/ 302 coeff = k*pow(h, ALPHA)*pow(normgradphi,BETA-2.);295 coeff = k*pow(h,alpha)*pow(normgradphi,beta-2.); 303 296 } 304 297 … … 468 461 /*Intermediaries*/ 469 462 IssmDouble dphi[3],h,k,phi; 470 IssmDouble alpha,beta,omega,h_r;463 IssmDouble h_r; 471 464 IssmDouble oceanLS,iceLS; 472 465 IssmDouble* xyz_list = NULL; 473 474 /*Hard coded coefficients*/475 const IssmPDouble ALPHA = 5./4.;476 const IssmPDouble BETA = 3./2.;477 466 478 467 /*Fetch number vertices for this element*/ … … 495 484 496 485 /*Retrieve all inputs and parameters*/ 497 bool islaminar; 498 element->FindParam(&islaminar,HydrologyIsLaminarEnum); 486 bool istransition; 487 element->FindParam(&istransition,HydrologyIsTransitionEnum); 488 IssmDouble alpha = element->FindParam(HydrologySheetAlphaEnum); 489 IssmDouble beta = element->FindParam(HydrologySheetBetaEnum); 490 IssmDouble omega = element->FindParam(HydrologyOmegaEnum); 499 491 element->GetVerticesCoordinates(&xyz_list); 500 492 IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum); 501 493 IssmDouble mu_water = element->FindParam(MaterialsMuWaterEnum); 502 494 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_(beta_input);506 495 Input *phi_input = element->GetInput(HydraulicPotentialEnum); _assert_(phi_input); 507 496 Input *hr_input = element->GetInput(HydrologyBumpHeightEnum); _assert_(hr_input); … … 521 510 hr_input->GetInputValue(&h_r,gauss); 522 511 k_input->GetInputValue(&k,gauss); 523 omega_input->GetInputValue(&omega,gauss);524 alpha_input->GetInputValue(&alpha,gauss);525 beta_input->GetInputValue(&beta,gauss);526 512 oceanLS_input->GetInputValue(&oceanLS,gauss); 527 513 iceLS_input->GetInputValue(&iceLS,gauss); … … 541 527 IssmDouble nu = mu_water/rho_water; 542 528 IssmDouble coeff; 543 if(is laminar==1 && omega>=AEPS){529 if(istransition==1 && omega>=AEPS){ 544 530 IssmDouble hratio = fabs(h/h_r); 545 531 IssmDouble coarg = 1. + 4.*pow(hratio,3-2*alpha)*omega*k*pow(h,3)*normgradphi/nu; … … 548 534 } 549 535 else { 550 coeff = k*pow(h, ALPHA)*pow(normgradphi,BETA-2.)/max(AEPS,h); // divide by h to get speed instead of discharge536 coeff = k*pow(h,alpha)*pow(normgradphi,beta-2.)/max(AEPS,h); // divide by h to get speed instead of discharge 551 537 } 552 538 -
issm/trunk-jpl/src/c/classes/Loads/Channel.cpp
r27926 r27936 20 20 21 21 #define C_W 4.22e3 /*specific heat capacity of water (J/kg/K)*/ 22 #define ALPHA_C 5./4.23 #define BETA_C 3./2.24 /*Make sure these are the same as in HydrologyGlaDSAnalysis::CreateKMatrix*/25 #define ALPHA_S 5./4.26 #define BETA_S 3./2.27 22 #define AEPS 2.2204460492503131E-015 28 23 … … 367 362 IssmDouble Jdet,v1,qc,fFactor,Afactor,Bfactor,Xifactor; 368 363 IssmDouble A,B,n,phi_old,phi,phi_0,dPw,ks,kc,Ngrad; 369 IssmDouble alpha_s,beta_s,omega,h_r;364 IssmDouble h_r; 370 365 IssmDouble H,h,b,dphi[2],dphids,dphimds,db[2],dbds; 371 366 IssmDouble xyz_list[NUMVERTICES][3]; … … 383 378 GetVerticesCoordinates(&xyz_list_tria[0][0],tria->vertices,3); 384 379 385 bool is laminar;386 element->FindParam(&is laminar,HydrologyIsLaminarEnum);380 bool istransition; 381 element->FindParam(&istransition,HydrologyIsTransitionEnum); 387 382 IssmDouble L = element->FindParam(MaterialsLatentheatEnum); 388 383 IssmDouble mu_water = element->FindParam(MaterialsMuWaterEnum); … … 394 389 IssmDouble alpha_c = element->FindParam(HydrologyChannelAlphaEnum); 395 390 IssmDouble beta_c = element->FindParam(HydrologyChannelBetaEnum); 391 IssmDouble alpha_s = element->FindParam(HydrologySheetAlphaEnum); 392 IssmDouble beta_s = element->FindParam(HydrologySheetBetaEnum); 393 IssmDouble omega = element->FindParam(HydrologyOmegaEnum); 396 394 397 395 Input* h_input = element->GetInput(HydrologySheetThicknessEnum); _assert_(h_input); … … 402 400 Input* ks_input = element->GetInput(HydrologySheetConductivityEnum); _assert_(ks_input); 403 401 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 402 Input* hr_input = element->GetInput(HydrologyBumpHeightEnum);_assert_(hr_input); 408 403 Input* phi_input = element->GetInput(HydraulicPotentialEnum); _assert_(phi_input); … … 432 427 ks_input->GetInputValue(&ks,gauss); 433 428 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 429 hr_input->GetInputValue(&h_r,gauss); 438 430 B_input->GetInputValue(&B,gauss); … … 451 443 if(Ngrad<AEPS) Ngrad = AEPS; 452 444 453 /*Compute the effective conductivity Kc = k h^alpha |grad Phi|^{beta-2} (same for sheet) and use laminarif specified*/445 /*Compute the effective conductivity Kc = k h^alpha |grad Phi|^{beta-2} (same for sheet) and use transition model if specified*/ 454 446 IssmDouble Kc; 455 447 IssmDouble Ks; 456 448 IssmDouble nu = mu_water/rho_water; 457 if(is laminar==1 && omega>=AEPS){449 if(istransition==1 && omega>=AEPS){ 458 450 IssmDouble hratio = h/h_r; 459 451 IssmDouble coarg = 1. + 4.*omega*pow(hratio,3-2*alpha_s)*ks*pow(h,3)*Ngrad/nu; … … 462 454 } 463 455 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.);456 Ks = ks*pow(h,alpha_s)*pow(Ngrad,beta_s-2.); 457 Kc = kc * pow(this->S,alpha_c) * pow(Ngrad,beta_c-2.); 466 458 } 467 459 … … 529 521 IssmDouble A,B,n,phi_old,phi,phi_0,dphimds,dphi[2]; 530 522 IssmDouble H,h,b,db[2],dphids,qc,dPw,ks,kc,Ngrad; 531 IssmDouble alpha_s,beta_s,omega,h_r;523 IssmDouble h_r; 532 524 IssmDouble xyz_list[NUMVERTICES][3]; 533 525 IssmDouble xyz_list_tria[3][3]; … … 542 534 GetVerticesCoordinates(&xyz_list_tria[0][0],tria->vertices,3); 543 535 544 bool is laminar;545 element->FindParam(&is laminar,HydrologyIsLaminarEnum);536 bool istransition; 537 element->FindParam(&istransition,HydrologyIsTransitionEnum); 546 538 IssmDouble L = element->FindParam(MaterialsLatentheatEnum); 547 539 IssmDouble mu_water = element->FindParam(MaterialsMuWaterEnum); … … 551 543 IssmDouble lc = element->FindParam(HydrologyChannelSheetWidthEnum); 552 544 IssmDouble c_t = element->FindParam(HydrologyPressureMeltCoefficientEnum); 545 IssmDouble alpha_s = element->FindParam(HydrologySheetAlphaEnum); 546 IssmDouble beta_s = element->FindParam(HydrologySheetBetaEnum); 547 IssmDouble omega = element->FindParam(HydrologyOmegaEnum); 553 548 554 549 Input* h_input = element->GetInput(HydrologySheetThicknessEnum); _assert_(h_input); … … 559 554 Input* ks_input = element->GetInput(HydrologySheetConductivityEnum); _assert_(ks_input); 560 555 Input* kc_input = element->GetInput(HydrologyChannelConductivityEnum); _assert_(kc_input); 561 Input* omega_input = element->GetInput(HydrologyOmegaEnum); _assert_(omega_input);562 556 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 557 Input* hr_input = element->GetInput(HydrologyBumpHeightEnum);_assert_(hr_input); 566 558 … … 585 577 ks_input->GetInputValue(&ks,gauss); 586 578 kc_input->GetInputValue(&kc,gauss); 587 omega_input->GetInputValue(&omega,gauss);588 579 B_input->GetInputValue(&B,gauss); 589 580 n_input->GetInputValue(&n,gauss); … … 591 582 b_input->GetInputValue(&b,gauss); 592 583 H_input->GetInputValue(&H,gauss); 593 alpha_s_input->GetInputValue(&alpha_s,gauss);594 beta_s_input->GetInputValue(&beta_s,gauss);595 584 hr_input->GetInputValue(&h_r,gauss); 596 585 … … 606 595 607 596 608 /*Approx. discharge in the sheet flowing folwing in the direction of the channel ofver a width lc, use laminarif specified*/597 /*Approx. discharge in the sheet flowing folwing in the direction of the channel ofver a width lc, use transition model if specified*/ 609 598 IssmDouble Ks; 610 if (is laminar==1 && omega>=AEPS){599 if (istransition==1 && omega>=AEPS){ 611 600 IssmDouble hratio = h/h_r; 612 601 IssmDouble nu = mu_water/rho_water; … … 615 604 } 616 605 else { 617 Ks = ks * pow(h, ALPHA_S) * pow(Ngrad,BETA_S-2.);606 Ks = ks * pow(h,alpha_s) * pow(Ngrad,beta_s-2.); 618 607 } 619 608 … … 683 672 /*Intermediaries */ 684 673 IssmDouble A,B,n,phi,phi_0,ks,kc,Ngrad; 685 IssmDouble alpha_s,beta_s,omega,h_r;674 IssmDouble h_r; 686 675 IssmDouble H,h,b,dphi[2],dphids,dphimds,db[2],dbds; 687 676 IssmDouble xyz_list[NUMVERTICES][3]; … … 692 681 GetVerticesCoordinates(&xyz_list_tria[0][0],tria->vertices,3); 693 682 694 bool is laminar;695 element->FindParam(&is laminar,HydrologyIsLaminarEnum);683 bool istransition; 684 element->FindParam(&istransition,HydrologyIsTransitionEnum); 696 685 IssmDouble L = element->FindParam(MaterialsLatentheatEnum); 697 686 IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum); … … 704 693 IssmDouble alpha_c = element->FindParam(HydrologyChannelAlphaEnum); 705 694 IssmDouble beta_c = element->FindParam(HydrologyChannelBetaEnum); 695 IssmDouble alpha_s = element->FindParam(HydrologySheetAlphaEnum); 696 IssmDouble beta_s = element->FindParam(HydrologySheetBetaEnum); 697 IssmDouble omega = element->FindParam(HydrologyOmegaEnum); 706 698 707 699 Input* h_input = element->GetInput(HydrologySheetThicknessEnum); _assert_(h_input); … … 712 704 Input* ks_input = element->GetInput(HydrologySheetConductivityEnum); _assert_(ks_input); 713 705 Input* kc_input = element->GetInput(HydrologyChannelConductivityEnum); _assert_(kc_input); 714 Input* omega_input = element->GetInput(HydrologyOmegaEnum); _assert_(omega_input);715 706 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 707 Input* hr_input = element->GetInput(HydrologyBumpHeightEnum);_assert_(hr_input); 719 708 … … 731 720 ks_input->GetInputValue(&ks,gauss); 732 721 kc_input->GetInputValue(&kc,gauss); 733 omega_input->GetInputValue(&omega,gauss);734 722 B_input->GetInputValue(&B,gauss); 735 723 n_input->GetInputValue(&n,gauss); … … 738 726 H_input->GetInputValue(&H,gauss); 739 727 hr_input->GetInputValue(&h_r,gauss); 740 alpha_s_input->GetInputValue(&alpha_s,gauss);741 beta_s_input->GetInputValue(&beta_s,gauss);742 728 743 729 … … 755 741 IssmDouble dPw = dphids - dphimds; 756 742 757 /*Approx. discharge in the sheet flowing folwing in the direction of the channel ofver a width lc, use laminarif necessary*/743 /*Approx. discharge in the sheet flowing folwing in the direction of the channel ofver a width lc, use transition model if necessary*/ 758 744 IssmDouble qc; 759 if (is laminar==1 && omega>=AEPS){745 if (istransition==1 && omega>=AEPS){ 760 746 IssmDouble hratio = h/h_r; 761 747 IssmDouble nu = mu_water/rho_water; … … 764 750 } 765 751 else { 766 qc = - ks * pow(h, ALPHA_S) * pow(Ngrad,BETA_S-2.) * dphids;752 qc = - ks * pow(h,alpha_s) * pow(Ngrad,beta_s-2.) * dphids; 767 753 } 768 754 … … 771 757 772 758 IssmDouble C = C_W*c_t*rho_water; 773 IssmDouble Qprime = -kc * pow(Ngrad, BETA_C-2.)*dphids;759 IssmDouble Qprime = -kc * pow(Ngrad,beta_c-2.)*dphids; 774 760 IssmDouble N = phi_0 - phi; 775 761 … … 788 774 789 775 IssmDouble alpha = 1./(rho_ice*L)*( 790 fabs(Qprime*pow(Snew, ALPHA_C-1.)*dphids)791 + C*Qprime*pow(Snew, ALPHA_C-1.)*dPw776 fabs(Qprime*pow(Snew,alpha_c-1.)*dphids) 777 + C*Qprime*pow(Snew,alpha_c-1.)*dPw 792 778 ) - 2./pow(n,n)*A*pow(fabs(N),n-1.)*N; 793 779 … … 811 797 812 798 /*Compute new channel discharge for output only*/ 813 IssmDouble Kc = kc * pow(this->S, ALPHA_C) * pow(Ngrad,BETA_C-2.);799 IssmDouble Kc = kc * pow(this->S,alpha_c) * pow(Ngrad,beta_c-2.); 814 800 this->discharge = -Kc*dphids; 815 801 -
issm/trunk-jpl/src/c/shared/Enum/Enum.vim
r27925 r27936 267 267 syn keyword cConstant HydrologyEnglacialVoidRatioEnum 268 268 syn keyword cConstant HydrologyIschannelsEnum 269 syn keyword cConstant HydrologyIs LaminarEnum269 syn keyword cConstant HydrologyIsTransitionEnum 270 270 syn keyword cConstant HydrologyIsWaterPressureArmaEnum 271 271 syn keyword cConstant HydrologyMeltFlagEnum … … 273 273 syn keyword cConstant HydrologyNumBasinsEnum 274 274 syn keyword cConstant HydrologyNumRequestedOutputsEnum 275 syn keyword cConstant HydrologyOmegaEnum 275 276 syn keyword cConstant HydrologyPressureMeltCoefficientEnum 276 277 syn keyword cConstant HydrologyRelaxationEnum 277 278 syn keyword cConstant HydrologyRequestedOutputsEnum 278 279 syn keyword cConstant HydrologySedimentKmaxEnum 280 syn keyword cConstant HydrologySheetAlphaEnum 281 syn keyword cConstant HydrologySheetBetaEnum 279 282 syn keyword cConstant HydrologyStepsPerStepEnum 280 283 syn keyword cConstant HydrologyStorageEnum … … 898 901 syn keyword cConstant HydrologyMoulinInputEnum 899 902 syn keyword cConstant HydrologyNeumannfluxEnum 900 syn keyword cConstant HydrologyOmegaEnum901 903 syn keyword cConstant HydrologyReynoldsEnum 902 syn keyword cConstant HydrologySheetAlphaEnum903 syn keyword cConstant HydrologySheetBetaEnum904 904 syn keyword cConstant HydrologySheetConductivityEnum 905 905 syn keyword cConstant HydrologySheetThicknessEnum -
issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
r27925 r27936 261 261 HydrologyEnglacialVoidRatioEnum, 262 262 HydrologyIschannelsEnum, 263 HydrologyIs LaminarEnum,263 HydrologyIsTransitionEnum, 264 264 HydrologyIsWaterPressureArmaEnum, 265 265 HydrologyMeltFlagEnum, … … 267 267 HydrologyNumBasinsEnum, 268 268 HydrologyNumRequestedOutputsEnum, 269 HydrologyOmegaEnum, 269 270 HydrologyPressureMeltCoefficientEnum, 270 271 HydrologyRelaxationEnum, 271 272 HydrologyRequestedOutputsEnum, 272 273 HydrologySedimentKmaxEnum, 274 HydrologySheetAlphaEnum, 275 HydrologySheetBetaEnum, 273 276 HydrologyStepsPerStepEnum, 274 277 HydrologyStorageEnum, … … 894 897 HydrologyMoulinInputEnum, 895 898 HydrologyNeumannfluxEnum, 896 HydrologyOmegaEnum,897 899 HydrologyReynoldsEnum, 898 HydrologySheetAlphaEnum,899 HydrologySheetBetaEnum,900 900 HydrologySheetConductivityEnum, 901 901 HydrologySheetThicknessEnum, -
issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
r27925 r27936 269 269 case HydrologyEnglacialVoidRatioEnum : return "HydrologyEnglacialVoidRatio"; 270 270 case HydrologyIschannelsEnum : return "HydrologyIschannels"; 271 case HydrologyIs LaminarEnum : return "HydrologyIsLaminar";271 case HydrologyIsTransitionEnum : return "HydrologyIsTransition"; 272 272 case HydrologyIsWaterPressureArmaEnum : return "HydrologyIsWaterPressureArma"; 273 273 case HydrologyMeltFlagEnum : return "HydrologyMeltFlag"; … … 275 275 case HydrologyNumBasinsEnum : return "HydrologyNumBasins"; 276 276 case HydrologyNumRequestedOutputsEnum : return "HydrologyNumRequestedOutputs"; 277 case HydrologyOmegaEnum : return "HydrologyOmega"; 277 278 case HydrologyPressureMeltCoefficientEnum : return "HydrologyPressureMeltCoefficient"; 278 279 case HydrologyRelaxationEnum : return "HydrologyRelaxation"; 279 280 case HydrologyRequestedOutputsEnum : return "HydrologyRequestedOutputs"; 280 281 case HydrologySedimentKmaxEnum : return "HydrologySedimentKmax"; 282 case HydrologySheetAlphaEnum : return "HydrologySheetAlpha"; 283 case HydrologySheetBetaEnum : return "HydrologySheetBeta"; 281 284 case HydrologyStepsPerStepEnum : return "HydrologyStepsPerStep"; 282 285 case HydrologyStorageEnum : return "HydrologyStorage"; … … 900 903 case HydrologyMoulinInputEnum : return "HydrologyMoulinInput"; 901 904 case HydrologyNeumannfluxEnum : return "HydrologyNeumannflux"; 902 case HydrologyOmegaEnum : return "HydrologyOmega";903 905 case HydrologyReynoldsEnum : return "HydrologyReynolds"; 904 case HydrologySheetAlphaEnum : return "HydrologySheetAlpha";905 case HydrologySheetBetaEnum : return "HydrologySheetBeta";906 906 case HydrologySheetConductivityEnum : return "HydrologySheetConductivity"; 907 907 case HydrologySheetThicknessEnum : return "HydrologySheetThickness"; -
issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim
r27925 r27936 260 260 syn keyword juliaConstC HydrologyEnglacialVoidRatioEnum 261 261 syn keyword juliaConstC HydrologyIschannelsEnum 262 syn keyword juliaConstC HydrologyIs LaminarEnum262 syn keyword juliaConstC HydrologyIsTransitionEnum 263 263 syn keyword juliaConstC HydrologyIsWaterPressureArmaEnum 264 264 syn keyword juliaConstC HydrologyMeltFlagEnum … … 266 266 syn keyword juliaConstC HydrologyNumBasinsEnum 267 267 syn keyword juliaConstC HydrologyNumRequestedOutputsEnum 268 syn keyword juliaConstC HydrologyOmegaEnum 268 269 syn keyword juliaConstC HydrologyPressureMeltCoefficientEnum 269 270 syn keyword juliaConstC HydrologyRelaxationEnum 270 271 syn keyword juliaConstC HydrologyRequestedOutputsEnum 271 272 syn keyword juliaConstC HydrologySedimentKmaxEnum 273 syn keyword juliaConstC HydrologySheetAlphaEnum 274 syn keyword juliaConstC HydrologySheetBetaEnum 272 275 syn keyword juliaConstC HydrologyStepsPerStepEnum 273 276 syn keyword juliaConstC HydrologyStorageEnum … … 891 894 syn keyword juliaConstC HydrologyMoulinInputEnum 892 895 syn keyword juliaConstC HydrologyNeumannfluxEnum 893 syn keyword juliaConstC HydrologyOmegaEnum894 896 syn keyword juliaConstC HydrologyReynoldsEnum 895 syn keyword juliaConstC HydrologySheetAlphaEnum896 syn keyword juliaConstC HydrologySheetBetaEnum897 897 syn keyword juliaConstC HydrologySheetConductivityEnum 898 898 syn keyword juliaConstC HydrologySheetThicknessEnum -
issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
r27925 r27936 275 275 else if (strcmp(name,"HydrologyEnglacialVoidRatio")==0) return HydrologyEnglacialVoidRatioEnum; 276 276 else if (strcmp(name,"HydrologyIschannels")==0) return HydrologyIschannelsEnum; 277 else if (strcmp(name,"HydrologyIs Laminar")==0) return HydrologyIsLaminarEnum;277 else if (strcmp(name,"HydrologyIsTransition")==0) return HydrologyIsTransitionEnum; 278 278 else if (strcmp(name,"HydrologyIsWaterPressureArma")==0) return HydrologyIsWaterPressureArmaEnum; 279 279 else if (strcmp(name,"HydrologyMeltFlag")==0) return HydrologyMeltFlagEnum; … … 281 281 else if (strcmp(name,"HydrologyNumBasins")==0) return HydrologyNumBasinsEnum; 282 282 else if (strcmp(name,"HydrologyNumRequestedOutputs")==0) return HydrologyNumRequestedOutputsEnum; 283 else if (strcmp(name,"HydrologyOmega")==0) return HydrologyOmegaEnum; 283 284 else if (strcmp(name,"HydrologyPressureMeltCoefficient")==0) return HydrologyPressureMeltCoefficientEnum; 284 285 else if (strcmp(name,"HydrologyRelaxation")==0) return HydrologyRelaxationEnum; 285 286 else if (strcmp(name,"HydrologyRequestedOutputs")==0) return HydrologyRequestedOutputsEnum; 286 287 else if (strcmp(name,"HydrologySedimentKmax")==0) return HydrologySedimentKmaxEnum; 288 else if (strcmp(name,"HydrologySheetAlpha")==0) return HydrologySheetAlphaEnum; 289 else if (strcmp(name,"HydrologySheetBeta")==0) return HydrologySheetBetaEnum; 287 290 else if (strcmp(name,"HydrologyStepsPerStep")==0) return HydrologyStepsPerStepEnum; 288 291 else if (strcmp(name,"HydrologyStorage")==0) return HydrologyStorageEnum; … … 380 383 else if (strcmp(name,"MassFluxSegments")==0) return MassFluxSegmentsEnum; 381 384 else if (strcmp(name,"MassFluxSegmentsPresent")==0) return MassFluxSegmentsPresentEnum; 382 else if (strcmp(name,"MasstransportHydrostaticAdjustment")==0) return MasstransportHydrostaticAdjustmentEnum;383 else if (strcmp(name,"MasstransportIsfreesurface")==0) return MasstransportIsfreesurfaceEnum;384 else if (strcmp(name,"MasstransportMinThickness")==0) return MasstransportMinThicknessEnum;385 385 else stage=4; 386 386 } 387 387 if(stage==4){ 388 if (strcmp(name,"MasstransportNumRequestedOutputs")==0) return MasstransportNumRequestedOutputsEnum; 388 if (strcmp(name,"MasstransportHydrostaticAdjustment")==0) return MasstransportHydrostaticAdjustmentEnum; 389 else if (strcmp(name,"MasstransportIsfreesurface")==0) return MasstransportIsfreesurfaceEnum; 390 else if (strcmp(name,"MasstransportMinThickness")==0) return MasstransportMinThicknessEnum; 391 else if (strcmp(name,"MasstransportNumRequestedOutputs")==0) return MasstransportNumRequestedOutputsEnum; 389 392 else if (strcmp(name,"MasstransportPenaltyFactor")==0) return MasstransportPenaltyFactorEnum; 390 393 else if (strcmp(name,"MasstransportRequestedOutputs")==0) return MasstransportRequestedOutputsEnum; … … 503 506 else if (strcmp(name,"SealevelchangeTidalK2")==0) return SealevelchangeTidalK2Enum; 504 507 else if (strcmp(name,"SealevelchangeTidalH2")==0) return SealevelchangeTidalH2Enum; 505 else if (strcmp(name,"SealevelchangeTidalL2")==0) return SealevelchangeTidalL2Enum;506 else if (strcmp(name,"SolidearthSettingsSealevelLoading")==0) return SolidearthSettingsSealevelLoadingEnum;507 else if (strcmp(name,"SolidearthSettingsGRD")==0) return SolidearthSettingsGRDEnum;508 508 else stage=5; 509 509 } 510 510 if(stage==5){ 511 if (strcmp(name,"SolidearthSettingsRunFrequency")==0) return SolidearthSettingsRunFrequencyEnum; 511 if (strcmp(name,"SealevelchangeTidalL2")==0) return SealevelchangeTidalL2Enum; 512 else if (strcmp(name,"SolidearthSettingsSealevelLoading")==0) return SolidearthSettingsSealevelLoadingEnum; 513 else if (strcmp(name,"SolidearthSettingsGRD")==0) return SolidearthSettingsGRDEnum; 514 else if (strcmp(name,"SolidearthSettingsRunFrequency")==0) return SolidearthSettingsRunFrequencyEnum; 512 515 else if (strcmp(name,"SolidearthSettingsTimeAcc")==0) return SolidearthSettingsTimeAccEnum; 513 516 else if (strcmp(name,"SolidearthSettingsHoriz")==0) return SolidearthSettingsHorizEnum; … … 626 629 else if (strcmp(name,"SmbSemicTmin")==0) return SmbSemicTminEnum; 627 630 else if (strcmp(name,"SmbSemicTmid")==0) return SmbSemicTmidEnum; 628 else if (strcmp(name,"SmbSemicTmax")==0) return SmbSemicTmaxEnum;629 else if (strcmp(name,"SmbStepsPerStep")==0) return SmbStepsPerStepEnum;630 else if (strcmp(name,"SmbSwIdx")==0) return SmbSwIdxEnum;631 631 else stage=6; 632 632 } 633 633 if(stage==6){ 634 if (strcmp(name,"SmbSWgrad")==0) return SmbSWgradEnum; 634 if (strcmp(name,"SmbSemicTmax")==0) return SmbSemicTmaxEnum; 635 else if (strcmp(name,"SmbStepsPerStep")==0) return SmbStepsPerStepEnum; 636 else if (strcmp(name,"SmbSwIdx")==0) return SmbSwIdxEnum; 637 else if (strcmp(name,"SmbSWgrad")==0) return SmbSWgradEnum; 635 638 else if (strcmp(name,"SmbT0dry")==0) return SmbT0dryEnum; 636 639 else if (strcmp(name,"SmbT0wet")==0) return SmbT0wetEnum; … … 749 752 else if (strcmp(name,"BasalforcingsGroundediceMeltingRateObs")==0) return BasalforcingsGroundediceMeltingRateObsEnum; 750 753 else if (strcmp(name,"BasalforcingsLinearBasinId")==0) return BasalforcingsLinearBasinIdEnum; 751 else if (strcmp(name,"BasalforcingsPerturbationMeltingRate")==0) return BasalforcingsPerturbationMeltingRateEnum;752 else if (strcmp(name,"BasalforcingsSpatialDeepwaterElevation")==0) return BasalforcingsSpatialDeepwaterElevationEnum;753 else if (strcmp(name,"BasalforcingsSpatialDeepwaterMeltingRate")==0) return BasalforcingsSpatialDeepwaterMeltingRateEnum;754 754 else stage=7; 755 755 } 756 756 if(stage==7){ 757 if (strcmp(name,"BasalforcingsSpatialUpperwaterElevation")==0) return BasalforcingsSpatialUpperwaterElevationEnum; 757 if (strcmp(name,"BasalforcingsPerturbationMeltingRate")==0) return BasalforcingsPerturbationMeltingRateEnum; 758 else if (strcmp(name,"BasalforcingsSpatialDeepwaterElevation")==0) return BasalforcingsSpatialDeepwaterElevationEnum; 759 else if (strcmp(name,"BasalforcingsSpatialDeepwaterMeltingRate")==0) return BasalforcingsSpatialDeepwaterMeltingRateEnum; 760 else if (strcmp(name,"BasalforcingsSpatialUpperwaterElevation")==0) return BasalforcingsSpatialUpperwaterElevationEnum; 758 761 else if (strcmp(name,"BasalforcingsSpatialUpperwaterMeltingRate")==0) return BasalforcingsSpatialUpperwaterMeltingRateEnum; 759 762 else if (strcmp(name,"BasalforcingsIsmip6BasinId")==0) return BasalforcingsIsmip6BasinIdEnum; … … 872 875 else if (strcmp(name,"FrictionM")==0) return FrictionMEnum; 873 876 else if (strcmp(name,"FrictionP")==0) return FrictionPEnum; 874 else if (strcmp(name,"FrictionPressureAdjustedTemperature")==0) return FrictionPressureAdjustedTemperatureEnum;875 else if (strcmp(name,"FrictionQ")==0) return FrictionQEnum;876 else if (strcmp(name,"FrictionSedimentCompressibilityCoefficient")==0) return FrictionSedimentCompressibilityCoefficientEnum;877 877 else stage=8; 878 878 } 879 879 if(stage==8){ 880 if (strcmp(name,"FrictionTillFrictionAngle")==0) return FrictionTillFrictionAngleEnum; 880 if (strcmp(name,"FrictionPressureAdjustedTemperature")==0) return FrictionPressureAdjustedTemperatureEnum; 881 else if (strcmp(name,"FrictionQ")==0) return FrictionQEnum; 882 else if (strcmp(name,"FrictionSedimentCompressibilityCoefficient")==0) return FrictionSedimentCompressibilityCoefficientEnum; 883 else if (strcmp(name,"FrictionTillFrictionAngle")==0) return FrictionTillFrictionAngleEnum; 881 884 else if (strcmp(name,"FrictionWaterLayer")==0) return FrictionWaterLayerEnum; 882 885 else if (strcmp(name,"FrictionWaterPressure")==0) return FrictionWaterPressureEnum; … … 921 924 else if (strcmp(name,"HydrologyMoulinInput")==0) return HydrologyMoulinInputEnum; 922 925 else if (strcmp(name,"HydrologyNeumannflux")==0) return HydrologyNeumannfluxEnum; 923 else if (strcmp(name,"HydrologyOmega")==0) return HydrologyOmegaEnum;924 926 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;927 927 else if (strcmp(name,"HydrologySheetConductivity")==0) return HydrologySheetConductivityEnum; 928 928 else if (strcmp(name,"HydrologySheetThickness")==0) return HydrologySheetThicknessEnum; -
issm/trunk-jpl/src/m/classes/hydrologyglads.m
r27926 r27936 11 11 cavity_spacing = 0.; 12 12 bump_height = NaN; 13 omega = 0; % TH14 sheet_alpha = 5.0/4.0; % TH15 sheet_beta = 3.0/2.0; % TH13 omega = 0; 14 sheet_alpha = NaN; 15 sheet_beta = NaN; 16 16 17 17 %Channels … … 19 19 channel_conductivity = NaN; 20 20 channel_sheet_width = 0.; 21 channel_alpha = 5.0/4.0; % TH22 channel_beta = 3.0/2.0; % TH21 channel_alpha = NaN; 22 channel_beta = NaN; 23 23 24 24 %Other … … 29 29 requested_outputs = {}; 30 30 melt_flag = 0; 31 is laminar= 0;31 istransition = 0; 32 32 end 33 33 methods … … 51 51 self.pressure_melt_coefficient = 7.5e-8; %K/Pa (See table 1 in Erder et al. 2013) 52 52 self.cavity_spacing = 2.; %m 53 self.omega = 1./2000.; % TH 53 self.sheet_alpha = 5.0/4.0; 54 self.sheet_beta = 3.0/2.0; 55 self.omega = 1./2000.; 54 56 55 57 %Channel parameters … … 57 59 self.channel_conductivity = 5.e-2; %Dow's default, Table uses 0.1 58 60 self.channel_sheet_width = 2.; %m 61 self.channel_alpha = 5.0/4.0; 62 self.channel_beta = 3.0/2.0; 59 63 60 %Otherself.omega = 1./2000.; % TH64 %Otherself.omega = 1./2000.; 61 65 self.englacial_void_ratio = 1.e-5;% Dow's default, Table from Werder et al. uses 1e-3; 62 66 self.requested_outputs={'default'}; 63 67 self.melt_flag=0; 64 self.is laminar= 0; %by default use GlaDS default turbulent code68 self.istransition = 0; %by default use GlaDS default turbulent code 65 69 end % }}} 66 70 function md = checkconsistency(self,md,solution,analyses) % {{{ … … 76 80 md = checkfield(md,'fieldname','hydrology.cavity_spacing','numel',[1],'>',0); 77 81 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); % TH79 md = checkfield(md,'fieldname','hydrology.sheet_alpha', 'numel', [1], '>', 0); % TH80 md = checkfield(md,'fieldname','hydrology.sheet_beta', 'numel', [1], '>', 0); % TH82 md = checkfield(md,'fieldname','hydrology.omega', 'numel', [1], '>=', 0); 83 md = checkfield(md,'fieldname','hydrology.sheet_alpha', 'numel', [1], '>', 0); 84 md = checkfield(md,'fieldname','hydrology.sheet_beta', 'numel', [1], '>', 0); 81 85 82 86 %Channels … … 84 88 md = checkfield(md,'fieldname','hydrology.channel_conductivity','size',[md.mesh.numberofvertices 1],'>=',0,'NaN',1,'Inf',1); 85 89 md = checkfield(md,'fieldname','hydrology.channel_sheet_width','numel',[1],'>=',0); 86 md = checkfield(md,'fieldname','hydrology.channel_alpha', 'numel', [1], '>', 0); % TH87 md = checkfield(md,'fieldname','hydrology.channel_beta', 'numel', [1], '>', 0); % TH90 md = checkfield(md,'fieldname','hydrology.channel_alpha', 'numel', [1], '>', 0); 91 md = checkfield(md,'fieldname','hydrology.channel_beta', 'numel', [1], '>', 0); 88 92 89 93 %Other … … 94 98 md = checkfield(md,'fieldname','hydrology.requested_outputs','stringrow',1); 95 99 md = checkfield(md,'fieldname','hydrology.melt_flag','numel',[1],'values',[0 1 2]); 96 md = checkfield(md,'fieldname','hydrology.is laminar','numel',[1],'values',[0 1]);100 md = checkfield(md,'fieldname','hydrology.istransition','numel',[1],'values',[0 1]); 97 101 if self.melt_flag==1 || self.melt_flag==2 98 102 md = checkfield(md,'fieldname','basalforcings.groundedice_melting_rate','NaN',1,'Inf',1,'timeseries',1); … … 104 108 fielddisplay(self,'pressure_melt_coefficient','Pressure melt coefficient (c_t) [K Pa^-1]'); 105 109 fielddisplay(self,'sheet_conductivity','sheet conductivity (k) [m^(7/4) kg^(-1/2)]'); 106 fielddisplay(self,'sheet_alpha','First sheet-flow exponent (alpha_s) []'); % TH107 fielddisplay(self,'sheet_beta','Second sheet-flow exponent (beta_s) []'); % TH110 fielddisplay(self,'sheet_alpha','First sheet-flow exponent (alpha_s) []'); 111 fielddisplay(self,'sheet_beta','Second sheet-flow exponent (beta_s) []'); 108 112 fielddisplay(self,'cavity_spacing','cavity spacing (l_r) [m]'); 109 113 fielddisplay(self,'bump_height','typical bump height (h_r) [m]'); 110 fielddisplay(self,'omega','transition parameter (omega) []'); % TH114 fielddisplay(self,'omega','transition parameter (omega) []'); 111 115 disp(sprintf(' CHANNELS')); 112 116 fielddisplay(self,'ischannels','Do we allow for channels? 1: yes, 0: no'); 113 117 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) []'); % TH115 fielddisplay(self,'channel_beta','Second channel-flow exponent (beta_s) []'); % TH118 fielddisplay(self,'channel_alpha','First channel-flow exponent (alpha_s) []'); 119 fielddisplay(self,'channel_beta','Second channel-flow exponent (beta_s) []'); 116 120 fielddisplay(self,'channel_sheet_width','channel sheet width [m]'); 117 121 disp(sprintf(' OTHER')); … … 122 126 fielddisplay(self,'requested_outputs','additional outputs requested'); 123 127 fielddisplay(self,'melt_flag','User specified basal melt? 0: no (default), 1: use md.basalforcings.groundedice_melting_rate'); 124 fielddisplay(self,'is laminar','do we use laminar[1] or turbulent physics [0, default]');128 fielddisplay(self,'istransition','do we use transition [1] or turbulent physics [0, default]'); 125 129 end % }}} 126 130 function marshall(self,prefix,md,fid) % {{{ … … 131 135 WriteData(fid,prefix,'name','md.hydrology.model','data',5,'format','Integer'); 132 136 133 %Sheet137 %Sheet 134 138 WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','pressure_melt_coefficient','format','Double'); 135 139 WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','sheet_conductivity','format','DoubleMat','mattype',1); 136 140 WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','cavity_spacing','format','Double'); 137 141 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'); % TH139 WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','sheet_alpha','format','Double'); % TH140 WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','sheet_beta','format','Double'); % TH142 WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','omega','format','Double'); 143 WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','sheet_alpha','format','Double'); 144 WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','sheet_beta','format','Double'); 141 145 142 146 %Channels … … 144 148 WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','channel_conductivity','format','DoubleMat','mattype',1); 145 149 WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','channel_sheet_width','format','Double'); 146 WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','channel_alpha','format','Double'); % TH147 WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','channel_beta','format','Double'); % TH150 WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','channel_alpha','format','Double'); 151 WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','channel_beta','format','Double'); 148 152 149 153 %Others … … 153 157 WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','englacial_void_ratio','format','Double'); 154 158 WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','melt_flag','format','Integer'); 155 WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','is laminar','format','Boolean');159 WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','istransition','format','Boolean'); 156 160 outputs = self.requested_outputs; 157 161 pos = find(ismember(outputs,'default')); -
issm/trunk-jpl/src/m/classes/hydrologyglads.py
r27926 r27936 14 14 15 15 def __init__(self, *args): # {{{ 16 16 # Sheet 17 17 self.pressure_melt_coefficient = 0. 18 18 self.sheet_conductivity = np.nan 19 19 self.cavity_spacing = 0. 20 20 self.bump_height = np.nan 21 self.omega = 0.; #TH22 self.sheet_alpha = 5.0/4.0; # TH23 self.sheet_beta = 3.0/2.0; # TH21 self.omega = 0.; 22 self.sheet_alpha = np.nan; 23 self.sheet_beta = np.nan; 24 24 25 25 # Channels … … 27 27 self.channel_conductivity = np.nan 28 28 self.channel_sheet_width = 0. 29 self.channel_alpha = 5.0/4.0; # TH30 self.channel_beta = 3.0/2.0; # TH29 self.channel_alpha = np.nan; 30 self.channel_beta = np.nan; 31 31 32 32 # Other … … 37 37 self.requested_outputs = [] 38 38 self.melt_flag = 0 39 self.is laminar= 039 self.istransition = 0 40 40 41 41 nargs = len(args) … … 72 72 s += '{}\n'.format(fielddisplay(self, 'requested_outputs', 'additional outputs requested')) 73 73 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, 'is laminar','do we use laminar[1] or turbulent physics [0, default]'))74 s += '{}\n'.format(fielddisplay(self, 'istransition','do we use transition [1] or turbulent physics [0, default]')) 75 75 return string 76 76 # }}} … … 96 96 self.pressure_melt_coefficient = 7.5e-8 #K / Pa (See table 1 in Erder et al. 2013) 97 97 self.cavity_spacing = 2. #m 98 self.omega = 1./2000.; #TH 98 self.sheet_alpha = 5.0/4.0; 99 self.sheet_beta = 3.0/2.0; 100 self.omega = 1./2000.; 99 101 100 102 # Channel parameters … … 102 104 self.channel_conductivity = 5.e-2 #Dow's default, Table uses 0.1 103 105 self.channel_sheet_width = 2. #m 104 self.islaminar = 0 #by default use turbulent physics 106 self.channel_alpha = 5.0/4.0; 107 self.channel_beta = 3.0/2.0; 105 108 106 109 # Other … … 108 111 self.requested_outputs = ['default'] 109 112 self.melt_flag = 0 110 self.is laminar= 0 #by default use turbulent physics113 self.istransition = 0 #by default use turbulent physics 111 114 112 115 return self … … 123 126 md = checkfield(md, 'fieldname', 'hydrology.cavity_spacing', 'numel', [1], '>', 0) 124 127 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); # TH126 md = checkfield(md,'fieldname','hydrology.sheet_alpha', 'numel', [1], '>', 0); # TH127 md = checkfield(md,'fieldname','hydrology.sheet_beta', 'numel', [1], '>', 0); # TH128 md = checkfield(md,'fieldname','hydrology.omega', 'numel', [1], '>=', 0); 129 md = checkfield(md,'fieldname','hydrology.sheet_alpha', 'numel', [1], '>', 0); 130 md = checkfield(md,'fieldname','hydrology.sheet_beta', 'numel', [1], '>', 0); 128 131 129 132 # Channels … … 131 134 md = checkfield(md, 'fieldname', 'hydrology.channel_conductivity', 'size', [md.mesh.numberofvertices], '>', 0) 132 135 md = checkfield(md, 'fieldname', 'hydrology.channel_sheet_width', 'numel', [1], '>=', 0) 133 md = checkfield(md,'fieldname','hydrology.channel_alpha', 'numel', [1], '>', 0); # TH134 md = checkfield(md,'fieldname','hydrology.channel_beta', 'numel', [1], '>', 0); # TH136 md = checkfield(md,'fieldname','hydrology.channel_alpha', 'numel', [1], '>', 0); 137 md = checkfield(md,'fieldname','hydrology.channel_beta', 'numel', [1], '>', 0); 135 138 136 139 # Other … … 141 144 md = checkfield(md, 'fieldname', 'hydrology.requested_outputs', 'stringrow', 1) 142 145 md = checkfield(md, 'fieldname', 'hydrology.melt_flag', 'numel', [1], 'values', [0, 1]) 143 md = checkfield(md, 'fieldname', 'hydrology.is laminar', 'numel', [1], 'values', [0, 1])146 md = checkfield(md, 'fieldname', 'hydrology.istransition', 'numel', [1], 'values', [0, 1]) 144 147 if self.melt_flag == 1 or self.melt_flag == 2: 145 148 md = checkfield(md, 'fieldname', 'basalforcings.groundedice_melting_rate', 'NaN', 1, 'Inf', 1, 'timeseries', 1) … … 151 154 WriteData(fid, prefix, 'name', 'md.hydrology.model', 'data', 5, 'format', 'Integer') 152 155 153 156 # Sheet 154 157 WriteData(fid, prefix, 'object', self, 'class', 'hydrology', 'fieldname', 'pressure_melt_coefficient', 'format', 'Double') 155 158 WriteData(fid, prefix, 'object', self, 'class', 'hydrology', 'fieldname', 'sheet_conductivity', 'format', 'DoubleMat', 'mattype', 1) 156 159 WriteData(fid, prefix, 'object', self, 'class', 'hydrology', 'fieldname', 'cavity_spacing', 'format', 'Double') 157 160 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'); # TH159 WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','sheet_alpha','format','Double'); # TH160 WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','sheet_beta','format','Double'); # TH161 WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','omega','format','Double'); 162 WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','sheet_alpha','format','Double'); 163 WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','sheet_beta','format','Double'); 161 164 162 165 # Channels … … 164 167 WriteData(fid, prefix, 'object', self, 'class', 'hydrology', 'fieldname', 'channel_conductivity', 'format', 'DoubleMat', 'mattype', 1) 165 168 WriteData(fid, prefix, 'object', self, 'class', 'hydrology', 'fieldname', 'channel_sheet_width', 'format', 'Double') 166 WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','channel_alpha','format','Double'); # TH167 WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','channel_beta','format','Double'); # TH169 WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','channel_alpha','format','Double'); 170 WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','channel_beta','format','Double'); 168 171 169 172 # Others … … 173 176 WriteData(fid, prefix, 'object', self, 'class', 'hydrology', 'fieldname', 'englacial_void_ratio', 'format', 'Double') 174 177 WriteData(fid, prefix, 'object', self, 'class', 'hydrology', 'fieldname', 'melt_flag', 'format', 'Integer') 175 WriteData(fid, prefix, 'object', self, 'class', 'hydrology', 'fieldname', 'is laminar', 'format', 'Boolean')178 WriteData(fid, prefix, 'object', self, 'class', 'hydrology', 'fieldname', 'istransition', 'format', 'Boolean') 176 179 177 180 outputs = self.requested_outputs
Note:
See TracChangeset
for help on using the changeset viewer.