Changeset 28039
- Timestamp:
- 01/04/24 05:49:30 (15 months ago)
- Location:
- issm/trunk-jpl/src/c/classes/Elements
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
r27964 r28039 286 286 IssmDouble groundedice,bed,sealevel; 287 287 288 /*Depth average Bfor stress calculation*/288 /*Depth average velocity for stress calculation*/ 289 289 this->InputDepthAverageAtBase(VxEnum,VxAverageEnum); 290 290 this->InputDepthAverageAtBase(VyEnum,VyAverageEnum); … … 625 625 delete gauss; 626 626 } 627 } 628 /*}}}*/ 629 void Penta::CalvingRateCalvingMIP(){/*{{{*/ 630 631 if(!this->IsOnBase()) return; 632 633 IssmDouble calvingrate[NUMVERTICES]; 634 IssmDouble calvingratex[NUMVERTICES]; 635 IssmDouble calvingratey[NUMVERTICES]; 636 int experiment = 1; /* exp:1 by default */ 637 int dim, domaintype; 638 IssmDouble vx, vy, vel, c, wrate; 639 IssmDouble time, groundedice, yts; 640 641 /*Get problem dimension and whether there is moving front or not*/ 642 this->FindParam(&domaintype,DomainTypeEnum); 643 this->FindParam(&time,TimeEnum); 644 this->FindParam(&yts,ConstantsYtsEnum); 645 646 switch(domaintype){ 647 case Domain2DverticalEnum: dim = 1; break; 648 case Domain2DhorizontalEnum: dim = 2; break; 649 case Domain3DEnum: dim = 2; break; 650 default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet"); 651 } 652 if(dim==1) _error_("not implemented in 1D..."); 653 654 /*Depth average velocity for stress calculation*/ 655 this->InputDepthAverageAtBase(VxEnum,VxAverageEnum); 656 this->InputDepthAverageAtBase(VyEnum,VyAverageEnum); 657 658 /*Retrieve all inputs and parameters we will need*/ 659 Input *vx_input = this->GetInput(VxAverageEnum); _assert_(vx_input); 660 Input *vy_input = this->GetInput(VyAverageEnum); _assert_(vy_input); 661 Input *wrate_input = this->GetInput(CalvingAblationrateEnum); _assert_(wrate_input); 662 Input* gr_input = this->GetInput(MaskOceanLevelsetEnum); _assert_(gr_input); 663 664 /* Use which experiment: use existing Enum */ 665 this->FindParam(&experiment, CalvingUseParamEnum); 666 667 /* Start looping on the number of vertices: */ 668 GaussPenta gauss; 669 for(int iv=0;iv<3;iv++){ 670 gauss.GaussVertex(iv); 671 672 /*Get velocity components */ 673 vx_input->GetInputValue(&vx,&gauss); 674 vy_input->GetInputValue(&vy,&gauss); 675 vel=sqrt(vx*vx+vy*vy)+1.e-14; 676 677 /* no calving for grounded ice in EXP4 */ 678 gr_input->GetInputValue(&groundedice,&gauss); 679 680 switch (experiment) { 681 case 1: 682 case 3: 683 /* Exp 1 and 3: set c=v-wrate, wrate=0, so that w=0 */ 684 wrate = 0.0; 685 break; 686 case 2: 687 /* Exp 2: set c=v-wrate(given)*/ 688 wrate_input->GetInputValue(&wrate,&gauss); 689 break; 690 case 4: 691 /* Exp 4: set c=v-wrate(given), for the first 500 years, then c=0 for the second 500 years*/ 692 if((groundedice<0) && (time<=500.0*yts)) { 693 // wrate_input->GetInputValue(&wrate,&gauss); 694 wrate = -750*sin(2.0*M_PI*time/yts/1000)/yts; // m/a -> m/s 695 } 696 else { 697 /* no calving on the grounded ice*/ 698 wrate = vel; 699 } 700 break; 701 default: 702 _error_("The experiment is not supported yet!"); 703 } 704 705 calvingrate[iv] = vel - wrate; 706 calvingratex[iv] = vx - wrate*vx/vel; 707 calvingratey[iv] = vy - wrate*vy/vel; 708 } 709 /*Add input*/ 710 this->AddBasalInput(CalvingCalvingrateEnum,&calvingrate[0],P1DGEnum); 711 this->AddBasalInput(CalvingratexEnum,&calvingratex[0],P1DGEnum); 712 this->AddBasalInput(CalvingrateyEnum,&calvingratey[0],P1DGEnum); 713 714 /*Extrude*/ 715 this->InputExtrude(CalvingCalvingrateEnum,-1); 716 this->InputExtrude(CalvingratexEnum,-1); 717 this->InputExtrude(CalvingrateyEnum,-1); 627 718 } 628 719 /*}}}*/ … … 2810 2901 case CalvingVonmisesEnum: 2811 2902 case CalvingLevermannEnum: 2903 case CalvingCalvingMIPEnum: 2812 2904 calvingratex_input=this->GetInput(CalvingratexEnum); _assert_(calvingratex_input); 2813 2905 calvingratey_input=this->GetInput(CalvingrateyEnum); _assert_(calvingratey_input); … … 2847 2939 case CalvingVonmisesEnum: 2848 2940 case CalvingLevermannEnum: 2941 case CalvingCalvingMIPEnum: 2849 2942 calvingratex_input->GetInputValue(&c[0],&gauss); 2850 2943 calvingratey_input->GetInputValue(&c[1],&gauss); -
issm/trunk-jpl/src/c/classes/Elements/Penta.h
r27964 r28039 58 58 void CalvingFluxLevelset(); 59 59 void CalvingMeltingFluxLevelset(); 60 void CalvingRateCalvingMIP(); 60 61 IssmDouble CharacteristicLength(void){_error_("not implemented yet");}; 61 62 void ComputeBasalStress(void);
Note:
See TracChangeset
for help on using the changeset viewer.