![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 20500 for issm/trunk/src/c/classes/Elements/Element.cpp
- Timestamp:
- 04/12/16 21:32:01 (9 years ago)
- Location:
- issm/trunk
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk
- Property svn:ignore
-
old new 1 build-js 2 build-esmf 3 build-gcm 1 4 build-fw 2 5 build-ad
-
- Property svn:mergeinfo changed
/issm/trunk-jpl merged: 19104,19106-19126,19128-19134,19136-19170,19172-19299,19302,19306-19405,19407-19604,19606-19668,19670-20496
- Property svn:ignore
-
issm/trunk/src
- Property svn:mergeinfo changed
-
issm/trunk/src/c/classes/Elements/Element.cpp
r19105 r20500 9 9 #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" 10 10 #endif 11 #include <math.h> 11 12 #include <stdio.h> 12 13 #include <string.h> 13 14 #include "../classes.h" 14 15 #include "../../shared/shared.h" 16 #include "../../modules/SurfaceMassBalancex/SurfaceMassBalancex.h" 15 17 /*}}}*/ 16 18 … … 46 48 IssmDouble eps_xx,eps_xy,eps_yy,eps_xz,eps_yz,eps_zz,eps_eff; 47 49 IssmDouble epsmin=1.e-27; 48 IssmDouble eps_0, eps_f,sigma_0,B,D,n;50 IssmDouble eps_0,kappa,sigma_0,B,D,n,envelopeD; 49 51 int dim,counter=0; 50 52 IssmDouble k1,k2,threshold=1.e-12; … … 54 56 this->ComputeStrainRate(); 55 57 parameters->FindParam(&dim,DomainDimensionEnum); 56 parameters->FindParam(& eps_f,DamageC1Enum);58 parameters->FindParam(&kappa,DamageKappaEnum); 57 59 parameters->FindParam(&sigma_0,DamageStressThresholdEnum); 58 60 … … 113 115 /* Compute threshold strain rate from threshold stress */ 114 116 eps_0=pow(sigma_0/B,n); 115 _assert_(eps_f>eps_0); 116 117 /* Compute kappa (k) from pre-existing level of damage, using Newton-Raphson */ 118 /* provide a reasonable initial guess */ 119 if(D==0){ 120 k1=eps_0; 121 } 122 else{ 123 k1=exp(n*eps_0/(eps_f-eps_0)-n*log(1.-D)+log(eps_0)); /* initial guess */ 124 } 125 126 counter=0; 127 while(true){ 128 /*Newton step k2=k1-f(k1)/f'(k1) */ 129 k2=k1-(k1+(eps_f-eps_0)/n*log(k1)-eps_0+(eps_f-eps_0)*(log(1.-D)-1./n*log(eps_0)))/(1.+(eps_f-eps_0)/n/k1); 130 131 if( fabs(k2-k1)/(fabs(k2))<threshold ){ 132 break; 133 } 134 else{ 135 k1=k2; 136 counter++; 137 } 138 139 if(counter>50) break; 140 } 141 142 if(eps_eff>k2){ 143 newD[i]=1.-pow(eps_0/eps_eff,1./n)*exp(-(eps_eff-eps_0)/(eps_f-eps_0)); 117 118 if(eps_eff>eps_0){ 119 /* Compute damage on envelope curve for existing level of effective strain rate */ 120 envelopeD=1.-pow(eps_0/eps_eff,1./n)*exp(-(eps_eff-eps_0)/(eps_0*(kappa-1.))); 121 122 if(envelopeD>D){ 123 newD[i]=envelopeD; 124 } 125 else newD[i]=D; 144 126 } 145 127 else newD[i]=D; … … 235 217 xDelete<IssmDouble>(eps_xz); 236 218 xDelete<IssmDouble>(eps_yz); 219 xDelete<IssmDouble>(eps_ef); 237 220 238 221 } … … 394 377 395 378 /*Clean up and return*/ 379 xDelete<IssmDouble>(xyz_list); 396 380 delete gauss; 397 381 return divergence; … … 503 487 *pdmudB=dmudB; 504 488 489 } 490 /*}}}*/ 491 void Element::Delta18oParameterization(void){/*{{{*/ 492 493 /*Are we on the base? If not, return*/ 494 if(!IsOnBase()) return; 495 496 int numvertices = this->GetNumberOfVertices(); 497 498 int i; 499 IssmDouble* monthlytemperatures=xNew<IssmDouble>(12*numvertices); 500 IssmDouble* monthlyprec=xNew<IssmDouble>(12*numvertices); 501 IssmDouble* TemperaturesPresentday=xNew<IssmDouble>(12*numvertices); 502 IssmDouble* TemperaturesLgm=xNew<IssmDouble>(12*numvertices); 503 IssmDouble* PrecipitationsPresentday=xNew<IssmDouble>(12*numvertices); 504 IssmDouble* tmp=xNew<IssmDouble>(numvertices); 505 IssmDouble Delta18oPresent,Delta18oLgm,Delta18oTime; 506 IssmDouble Delta18oSurfacePresent,Delta18oSurfaceLgm,Delta18oSurfaceTime; 507 IssmDouble time,yts,finaltime,time_yr; 508 509 /*Recover parameters*/ 510 this->parameters->FindParam(&time,TimeEnum); 511 this->parameters->FindParam(&yts,ConstantsYtsEnum); 512 this->parameters->FindParam(&finaltime,TimesteppingFinalTimeEnum); 513 time_yr=floor(time/yts)*yts; 514 515 /*Recover present day temperature and precipitation*/ 516 Input* input=this->inputs->GetInput(SmbTemperaturesPresentdayEnum); _assert_(input); 517 Input* input2=this->inputs->GetInput(SmbTemperaturesLgmEnum); _assert_(input2); 518 Input* input3=this->inputs->GetInput(SmbPrecipitationsPresentdayEnum); _assert_(input3); 519 /*loop over vertices: */ 520 Gauss* gauss=this->NewGauss(); 521 for(int month=0;month<12;month++){ 522 for(int iv=0;iv<numvertices;iv++){ 523 gauss->GaussVertex(iv); 524 input->GetInputValue(&TemperaturesPresentday[iv*12+month],gauss,month/12.*yts); 525 input2->GetInputValue(&TemperaturesLgm[iv*12+month],gauss,month/12.*yts); 526 input3->GetInputValue(&PrecipitationsPresentday[iv*12+month],gauss,month/12.*yts); 527 528 PrecipitationsPresentday[iv*12+month]=PrecipitationsPresentday[iv*12+month]*yts; 529 } 530 } 531 532 /*Recover delta18o and Delta18oSurface at present day, lgm and at time t*/ 533 this->parameters->FindParam(&Delta18oPresent,SmbDelta18oEnum,finaltime); 534 this->parameters->FindParam(&Delta18oLgm,SmbDelta18oEnum,(finaltime-(21000*yts))); 535 this->parameters->FindParam(&Delta18oTime,SmbDelta18oEnum,time); 536 this->parameters->FindParam(&Delta18oSurfacePresent,SmbDelta18oSurfaceEnum,finaltime); 537 this->parameters->FindParam(&Delta18oSurfaceLgm,SmbDelta18oSurfaceEnum,(finaltime-(21000*yts))); 538 this->parameters->FindParam(&Delta18oSurfaceTime,SmbDelta18oSurfaceEnum,time); 539 540 /*Compute the temperature and precipitation*/ 541 for(int iv=0;iv<numvertices;iv++){ 542 ComputeDelta18oTemperaturePrecipitation(Delta18oSurfacePresent, Delta18oSurfaceLgm, Delta18oSurfaceTime, 543 Delta18oPresent, Delta18oLgm, Delta18oTime, 544 &PrecipitationsPresentday[iv*12], 545 &TemperaturesLgm[iv*12], &TemperaturesPresentday[iv*12], 546 &monthlytemperatures[iv*12], &monthlyprec[iv*12]); 547 } 548 549 /*Update inputs*/ 550 TransientInput* NewTemperatureInput = new TransientInput(SmbMonthlytemperaturesEnum); 551 TransientInput* NewPrecipitationInput = new TransientInput(SmbPrecipitationEnum); 552 for (int imonth=0;imonth<12;imonth++) { 553 for(i=0;i<numvertices;i++) tmp[i]=monthlytemperatures[i*12+imonth]; 554 switch(this->ObjectEnum()){ 555 case TriaEnum: NewTemperatureInput->AddTimeInput(new TriaInput(SmbMonthlytemperaturesEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break; 556 case PentaEnum: NewTemperatureInput->AddTimeInput(new PentaInput(SmbMonthlytemperaturesEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break; 557 case TetraEnum: NewTemperatureInput->AddTimeInput(new TetraInput(SmbMonthlytemperaturesEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break; 558 default: _error_("Not implemented yet"); 559 } 560 for(i=0;i<numvertices;i++) tmp[i]=monthlyprec[i*12+imonth]/yts; 561 switch(this->ObjectEnum()){ 562 case TriaEnum: NewPrecipitationInput->AddTimeInput(new TriaInput(SmbPrecipitationEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break; 563 case PentaEnum: NewPrecipitationInput->AddTimeInput(new PentaInput(SmbPrecipitationEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break; 564 case TetraEnum: NewPrecipitationInput->AddTimeInput(new TetraInput(SmbPrecipitationEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break; 565 default: _error_("Not implemented yet"); 566 } 567 } 568 NewTemperatureInput->Configure(this->parameters); 569 NewPrecipitationInput->Configure(this->parameters); 570 571 this->inputs->AddInput(NewTemperatureInput); 572 this->inputs->AddInput(NewPrecipitationInput); 573 574 switch(this->ObjectEnum()){ 575 case TriaEnum: break; 576 case PentaEnum: 577 case TetraEnum: 578 this->InputExtrude(SmbMonthlytemperaturesEnum,-1); 579 this->InputExtrude(SmbPrecipitationEnum,-1); 580 break; 581 default: _error_("Not implemented yet"); 582 } 583 584 /*clean-up*/ 585 delete gauss; 586 xDelete<IssmDouble>(monthlytemperatures); 587 xDelete<IssmDouble>(monthlyprec); 588 xDelete<IssmDouble>(TemperaturesPresentday); 589 xDelete<IssmDouble>(TemperaturesLgm); 590 xDelete<IssmDouble>(PrecipitationsPresentday); 591 xDelete<IssmDouble>(tmp); 592 593 } 594 /*}}}*/ 595 void Element::MungsmtpParameterization(void){/*{{{*/ 596 /*Are we on the base? If not, return*/ 597 if(!IsOnBase()) return; 598 599 int numvertices = this->GetNumberOfVertices(); 600 601 int i; 602 IssmDouble* monthlytemperatures=xNew<IssmDouble>(12*numvertices); 603 IssmDouble* monthlyprec=xNew<IssmDouble>(12*numvertices); 604 IssmDouble* TemperaturesPresentday=xNew<IssmDouble>(12*numvertices); 605 IssmDouble* TemperaturesLgm=xNew<IssmDouble>(12*numvertices); 606 IssmDouble* PrecipitationsPresentday=xNew<IssmDouble>(12*numvertices); 607 IssmDouble* PrecipitationsLgm=xNew<IssmDouble>(12*numvertices); 608 IssmDouble* tmp=xNew<IssmDouble>(numvertices); 609 IssmDouble TdiffTime,PfacTime; 610 IssmDouble time,yts,time_yr; 611 this->parameters->FindParam(&time,TimeEnum); 612 this->parameters->FindParam(&yts,ConstantsYtsEnum); 613 time_yr=floor(time/yts)*yts; 614 615 /*Recover present day temperature and precipitation*/ 616 Input* input=this->inputs->GetInput(SmbTemperaturesPresentdayEnum); _assert_(input); 617 Input* input2=this->inputs->GetInput(SmbTemperaturesLgmEnum); _assert_(input2); 618 Input* input3=this->inputs->GetInput(SmbPrecipitationsPresentdayEnum); _assert_(input3); 619 Input* input4=this->inputs->GetInput(SmbPrecipitationsLgmEnum); _assert_(input4); 620 /*loop over vertices: */ 621 Gauss* gauss=this->NewGauss(); 622 for(int month=0;month<12;month++) { 623 for(int iv=0;iv<numvertices;iv++) { 624 gauss->GaussVertex(iv); 625 input->GetInputValue(&TemperaturesPresentday[iv*12+month],gauss,month/12.*yts); 626 input2->GetInputValue(&TemperaturesLgm[iv*12+month],gauss,month/12.*yts); 627 input3->GetInputValue(&PrecipitationsPresentday[iv*12+month],gauss,month/12.*yts); 628 input4->GetInputValue(&PrecipitationsLgm[iv*12+month],gauss,month/12.*yts); 629 630 PrecipitationsPresentday[iv*12+month]=PrecipitationsPresentday[iv*12+month]*yts; 631 PrecipitationsLgm[iv*12+month]=PrecipitationsLgm[iv*12+month]*yts; 632 } 633 } 634 635 /*Recover interpolation parameters at time t*/ 636 this->parameters->FindParam(&TdiffTime,SmbTdiffEnum,time); 637 this->parameters->FindParam(&PfacTime,SmbPfacEnum,time); 638 639 /*Compute the temperature and precipitation*/ 640 for(int iv=0;iv<numvertices;iv++){ 641 ComputeMungsmTemperaturePrecipitation(TdiffTime,PfacTime, 642 &PrecipitationsLgm[iv*12],&PrecipitationsPresentday[iv*12], 643 &TemperaturesLgm[iv*12], &TemperaturesPresentday[iv*12], 644 &monthlytemperatures[iv*12], &monthlyprec[iv*12]); 645 } 646 647 /*Update inputs*/ 648 TransientInput* NewTemperatureInput = new TransientInput(SmbMonthlytemperaturesEnum); 649 TransientInput* NewPrecipitationInput = new TransientInput(SmbPrecipitationEnum); 650 for (int imonth=0;imonth<12;imonth++) { 651 for(i=0;i<numvertices;i++) tmp[i]=monthlytemperatures[i*12+imonth]; 652 switch(this->ObjectEnum()){ 653 case TriaEnum: NewTemperatureInput->AddTimeInput(new TriaInput(SmbMonthlytemperaturesEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break; 654 case PentaEnum: NewTemperatureInput->AddTimeInput(new PentaInput(SmbMonthlytemperaturesEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break; 655 case TetraEnum: NewTemperatureInput->AddTimeInput(new TetraInput(SmbMonthlytemperaturesEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break; 656 default: _error_("Not implemented yet"); 657 } 658 for(i=0;i<numvertices;i++) tmp[i]=monthlyprec[i*12+imonth]/yts; 659 switch(this->ObjectEnum()){ 660 case TriaEnum: NewPrecipitationInput->AddTimeInput(new TriaInput(SmbPrecipitationEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break; 661 case PentaEnum: NewPrecipitationInput->AddTimeInput(new PentaInput(SmbPrecipitationEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break; 662 case TetraEnum: NewPrecipitationInput->AddTimeInput(new TetraInput(SmbPrecipitationEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break; 663 default: _error_("Not implemented yet"); 664 } 665 } 666 NewTemperatureInput->Configure(this->parameters); 667 NewPrecipitationInput->Configure(this->parameters); 668 669 this->inputs->AddInput(NewTemperatureInput); 670 this->inputs->AddInput(NewPrecipitationInput); 671 672 switch(this->ObjectEnum()){ 673 case TriaEnum: break; 674 case PentaEnum: 675 case TetraEnum: 676 this->InputExtrude(SmbMonthlytemperaturesEnum,-1); 677 this->InputExtrude(SmbPrecipitationEnum,-1); 678 break; 679 default: _error_("Not implemented yet"); 680 } 681 682 /*clean-up*/ 683 delete gauss; 684 xDelete<IssmDouble>(monthlytemperatures); 685 xDelete<IssmDouble>(monthlyprec); 686 xDelete<IssmDouble>(TemperaturesPresentday); 687 xDelete<IssmDouble>(TemperaturesLgm); 688 xDelete<IssmDouble>(PrecipitationsPresentday); 689 xDelete<IssmDouble>(PrecipitationsLgm); 690 xDelete<IssmDouble>(tmp); 691 692 } 693 /*}}}*/ 694 void Element::Delta18opdParameterization(void){/*{{{*/ 695 /*Are we on the base? If not, return*/ 696 if(!IsOnBase()) return; 697 698 int numvertices = this->GetNumberOfVertices(); 699 700 int i; 701 IssmDouble* monthlytemperatures=xNew<IssmDouble>(12*numvertices); 702 IssmDouble* monthlyprec=xNew<IssmDouble>(12*numvertices); 703 IssmDouble* TemperaturesPresentday=xNew<IssmDouble>(12*numvertices); 704 IssmDouble* PrecipitationsPresentday=xNew<IssmDouble>(12*numvertices); 705 IssmDouble* tmp=xNew<IssmDouble>(numvertices); 706 IssmDouble Delta18oTime; 707 IssmDouble dpermil; 708 IssmDouble time,yts,time_yr,month; 709 this->parameters->FindParam(&time,TimeEnum); 710 this->parameters->FindParam(&yts,ConstantsYtsEnum); 711 time_yr=floor(time/yts)*yts; 712 713 /*Get some pdd parameters*/ 714 dpermil=this->matpar->GetMaterialParameter(SmbDpermilEnum); 715 716 /*Recover present day temperature and precipitation*/ 717 Input* input=this->inputs->GetInput(SmbTemperaturesPresentdayEnum); _assert_(input); 718 Input* input2=this->inputs->GetInput(SmbPrecipitationsPresentdayEnum); _assert_(input2); 719 /*loop over vertices: */ 720 Gauss* gauss=this->NewGauss(); 721 for(int month=0;month<12;month++) { 722 for(int iv=0;iv<numvertices;iv++) { 723 gauss->GaussVertex(iv); 724 input->GetInputValue(&TemperaturesPresentday[iv*12+month],gauss,month/12.*yts); 725 input2->GetInputValue(&PrecipitationsPresentday[iv*12+month],gauss,month/12.*yts); 726 727 PrecipitationsPresentday[iv*12+month]=PrecipitationsPresentday[iv*12+month]*yts; 728 } 729 } 730 731 /*Recover interpolation parameters at time t*/ 732 this->parameters->FindParam(&Delta18oTime,SmbDelta18oEnum,time); 733 734 /*Compute the temperature and precipitation*/ 735 for(int iv=0;iv<numvertices;iv++){ 736 ComputeD18OTemperaturePrecipitationFromPD(Delta18oTime,dpermil, 737 &PrecipitationsPresentday[iv*12], &TemperaturesPresentday[iv*12], 738 &monthlytemperatures[iv*12], &monthlyprec[iv*12]); 739 } 740 741 /*Update inputs*/ 742 TransientInput* NewTemperatureInput = new TransientInput(SmbMonthlytemperaturesEnum); 743 TransientInput* NewPrecipitationInput = new TransientInput(SmbPrecipitationEnum); 744 for (int imonth=0;imonth<12;imonth++) { 745 for(i=0;i<numvertices;i++) tmp[i]=monthlytemperatures[i*12+imonth]; 746 switch(this->ObjectEnum()){ 747 case TriaEnum: NewTemperatureInput->AddTimeInput(new TriaInput(SmbMonthlytemperaturesEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break; 748 case PentaEnum: NewTemperatureInput->AddTimeInput(new PentaInput(SmbMonthlytemperaturesEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break; 749 case TetraEnum: NewTemperatureInput->AddTimeInput(new TetraInput(SmbMonthlytemperaturesEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break; 750 default: _error_("Not implemented yet"); 751 } 752 for(i=0;i<numvertices;i++) tmp[i]=monthlyprec[i*12+imonth]/yts; 753 switch(this->ObjectEnum()){ 754 case TriaEnum: NewPrecipitationInput->AddTimeInput(new TriaInput(SmbPrecipitationEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break; 755 case PentaEnum: NewPrecipitationInput->AddTimeInput(new PentaInput(SmbPrecipitationEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break; 756 case TetraEnum: NewPrecipitationInput->AddTimeInput(new TetraInput(SmbPrecipitationEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break; 757 default: _error_("Not implemented yet"); 758 } 759 } 760 NewTemperatureInput->Configure(this->parameters); 761 NewPrecipitationInput->Configure(this->parameters); 762 763 this->inputs->AddInput(NewTemperatureInput); 764 this->inputs->AddInput(NewPrecipitationInput); 765 766 switch(this->ObjectEnum()){ 767 case TriaEnum: break; 768 case PentaEnum: 769 case TetraEnum: 770 this->InputExtrude(SmbMonthlytemperaturesEnum,-1); 771 this->InputExtrude(SmbPrecipitationEnum,-1); 772 break; 773 default: _error_("Not implemented yet"); 774 } 775 776 /*clean-up*/ 777 delete gauss; 778 xDelete<IssmDouble>(monthlytemperatures); 779 xDelete<IssmDouble>(monthlyprec); 780 xDelete<IssmDouble>(TemperaturesPresentday); 781 xDelete<IssmDouble>(PrecipitationsPresentday); 782 xDelete<IssmDouble>(tmp); 783 505 784 } 506 785 /*}}}*/ … … 1090 1369 } 1091 1370 else if(vector_type==2){ //element vector 1371 1372 IssmDouble value; 1373 1092 1374 /*Are we in transient or static? */ 1093 1375 if(M==iomodel->numberofelements){ … … 1103 1385 else _error_("could not recognize nature of vector from code " << code); 1104 1386 } 1105 else { 1106 _error_("transient element inputs not supported yet!"); 1107 } 1108 } 1109 else{ 1110 _error_("Cannot add input for vector type " << vector_type << " (not supported)"); 1111 } 1112 }/*}}}*/ 1387 else if(M==iomodel->numberofelements+1){ 1388 /*create transient input: */ 1389 IssmDouble* times = xNew<IssmDouble>(N); 1390 for(t=0;t<N;t++) times[t] = vector[(M-1)*N+t]; 1391 TransientInput* transientinput=new TransientInput(vector_enum,times,N); 1392 TriaInput* bof=NULL; 1393 for(t=0;t<N;t++){ 1394 value=vector[N*this->Sid()+t]; 1395 switch(this->ObjectEnum()){ 1396 case TriaEnum: transientinput->AddTimeInput(new TriaInput( vector_enum,&value,P0Enum)); break; 1397 case PentaEnum: transientinput->AddTimeInput(new PentaInput(vector_enum,&value,P0Enum)); break; 1398 case TetraEnum: transientinput->AddTimeInput(new TetraInput(vector_enum,&value,P0Enum)); break; 1399 default: _error_("Not implemented yet"); 1400 } 1401 } 1402 this->inputs->AddInput(transientinput); 1403 xDelete<IssmDouble>(times); 1404 } 1405 else _error_("element vector is either numberofelements or numberofelements+1 long. Field provided (" << EnumToStringx(vector_enum) << ") is " << M << " long"); 1406 } 1407 else _error_("Cannot add input for vector type " << vector_type << " (not supported)"); 1408 } 1409 /*}}}*/ 1113 1410 void Element::InputDuplicate(int original_enum,int new_enum){/*{{{*/ 1114 1411 … … 1169 1466 bool Element::IsIceInElement(){/*{{{*/ 1170 1467 return (this->inputs->Min(MaskIceLevelsetEnum)<0.); 1468 } 1469 /*}}}*/ 1470 bool Element::IsWaterInElement(){/*{{{*/ 1471 return (this->inputs->Max(MaskOceanLevelsetEnum)>0.); 1472 } 1473 /*}}}*/ 1474 bool Element::IsLandInElement(){/*{{{*/ 1475 return (this->inputs->Max(MaskLandLevelsetEnum)>0.); 1171 1476 } 1172 1477 /*}}}*/ … … 1182 1487 name==SurfaceSlopeXEnum || 1183 1488 name==SurfaceSlopeYEnum || 1184 name==SurfaceforcingsMassBalanceEnum || 1489 name==SmbMassBalanceEnum || 1490 name==SmbAccumulationEnum || 1491 name==SmbRunoffEnum || 1492 name==SmbMeltEnum || 1493 name==SmbRefreezeEnum || 1494 name==SmbEvaporationEnum || 1185 1495 name==BasalforcingsGroundediceMeltingRateEnum || 1186 1496 name==BasalforcingsFloatingiceMeltingRateEnum || … … 1200 1510 name==InversionVzObsEnum || 1201 1511 name==TemperatureEnum || 1512 name==TemperaturePDDEnum || 1202 1513 name==EnthalpyEnum || 1203 1514 name==EnthalpyPicardEnum || … … 1206 1517 name==FrictionCoefficientEnum || 1207 1518 name==FrictionAsEnum || 1519 name==FrictionEffectivePressureEnum || 1208 1520 name==MaskGroundediceLevelsetEnum || 1209 1521 name==MaskIceLevelsetEnum || … … 1219 1531 name==MaterialsRheologyBbarEnum || 1220 1532 name==MaterialsRheologyNEnum || 1533 name==SealevelEnum || 1534 name==SealevelEustaticEnum || 1535 name==SealevelriseDeltathicknessEnum || 1221 1536 name==GiaWEnum || 1222 1537 name==GiadWdtEnum || … … 1263 1578 1264 1579 }/*}}}*/ 1580 void Element::MantlePlumeGeothermalFlux(){/*{{{*/ 1581 1582 int numvertices = this->GetNumberOfVertices(); 1583 IssmDouble mantleconductivity,nusselt,dtbg,plumeradius,topplumedepth,bottomplumedepth,plumex,plumey; 1584 IssmDouble crustthickness,uppercrustthickness,uppercrustheat,lowercrustheat; 1585 IssmDouble crustheat,plumeheat,dt,middleplumedepth,a,e,eprime,A0,lambda,Alambda,dAlambda; 1586 IssmDouble x,y,z,c; 1587 IssmDouble* values = xNew<IssmDouble>(numvertices); 1588 IssmDouble *xyz_list = NULL; 1589 1590 parameters->FindParam(&mantleconductivity,BasalforcingsMantleconductivityEnum); 1591 parameters->FindParam(&nusselt,BasalforcingsNusseltEnum); 1592 parameters->FindParam(&dtbg,BasalforcingsDtbgEnum); 1593 parameters->FindParam(&plumeradius,BasalforcingsPlumeradiusEnum); 1594 parameters->FindParam(&topplumedepth,BasalforcingsTopplumedepthEnum); 1595 parameters->FindParam(&bottomplumedepth,BasalforcingsBottomplumedepthEnum); 1596 parameters->FindParam(&plumex,BasalforcingsPlumexEnum); 1597 parameters->FindParam(&plumey,BasalforcingsPlumeyEnum); 1598 parameters->FindParam(&crustthickness,BasalforcingsCrustthicknessEnum); 1599 parameters->FindParam(&uppercrustthickness,BasalforcingsUppercrustthicknessEnum); 1600 parameters->FindParam(&uppercrustheat,BasalforcingsUppercrustheatEnum); 1601 parameters->FindParam(&lowercrustheat,BasalforcingsLowercrustheatEnum); 1602 1603 this->GetVerticesCoordinates(&xyz_list); 1604 c=plumeradius; 1605 a=(bottomplumedepth-topplumedepth)/2.; 1606 e=pow(a*a-c*c,1./2.)/a; 1607 A0=(1-pow(e,2.))/pow(e,3.)*(1./2.*log((1+e)/(1-e))-e); 1608 for(int i=0;i<numvertices;i++){ 1609 y=xyz_list[i*3+0]-plumex; 1610 z=xyz_list[i*3+1]-plumey; 1611 x=-(a+topplumedepth+crustthickness); 1612 lambda=(-pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2)+sqrt(pow(-pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2),2)+4*(pow(c,2)*pow(x,2)+pow(a,2)*(-pow(c,2)+pow(y,2)+pow(z,2)))))/2; 1613 dAlambda=(-8*a*pow(c,2)*x*(-2*pow(a,2)+2*pow(c,2)+sqrt(2)*sqrt((a-c)*(a+c))*sqrt(pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2)+sqrt(pow(-pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2),2)+4*(pow(c,2)*pow(x,2)+pow(a,2)*(-pow(c,2)+pow(y,2)+pow(z,2))))))*(pow(a,4)*(pow(y,2)+pow(z,2))+pow(c,4)*(pow(y,2)+pow(z,2))+pow(pow(x,2)+pow(y,2)+pow(z,2),2)*(pow(x,2)+pow(y,2)+pow(z,2)+sqrt(pow(-pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2),2)+4*(pow(c,2)*pow(x,2)+pow(a,2)*(-pow(c,2)+pow(y,2)+pow(z,2)))))+pow(c,2)*(pow(x,4)-pow(x,2)*(pow(y,2)+pow(z,2))-(pow(y,2)+pow(z,2))*(2*pow(y,2)+2*pow(z,2)+sqrt(pow(-pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2),2)+4*(pow(c,2)*pow(x,2)+pow(a,2)*(-pow(c,2)+pow(y,2)+pow(z,2))))))+pow(a,2)*(-pow(x,4)+pow(x,2)*(pow(y,2)+pow(z,2))+(pow(y,2)+pow(z,2))*(-2*pow(c,2)+2*(pow(y,2)+pow(z,2))+sqrt(pow(-pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2),2)+4*(pow(c,2)*pow(x,2)+pow(a,2)*(-pow(c,2)+pow(y,2)+pow(z,2))))))))/(sqrt((a-c)*(a+c))*sqrt(pow(-pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2),2)+4*(pow(c,2)*pow(x,2)+pow(a,2)*(-pow(c,2)+pow(y,2)+pow(z,2))))*pow(pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2)+sqrt(pow(-pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2),2)+4*(pow(c,2)*pow(x,2)+pow(a,2)*(-pow(c,2)+pow(y,2)+pow(z,2)))),3.5)*pow(-(sqrt(2)*sqrt((a-c)*(a+c)))+sqrt(pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2)+sqrt(pow(-pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2),2)+4*(pow(c,2)*pow(x,2)+pow(a,2)*(-pow(c,2)+pow(y,2)+pow(z,2))))),2)*(sqrt(2)*sqrt((a-c)*(a+c))+sqrt(pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2)+sqrt(pow(-pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2),2)+4*(pow(c,2)*pow(x,2)+pow(a,2)*(-pow(c,2)+pow(y,2)+pow(z,2))))))); 1614 eprime=pow((a*a-plumeradius*plumeradius)/(a*a+lambda),1./2.); 1615 Alambda=(1.-e*e)/(e*e*e)*(1./2.*log((1.+eprime)/(1.-eprime))-eprime); 1616 dt=dtbg-(nusselt-1.)/(1.+A0*(nusselt-1.))*(Alambda*dtbg+x*dtbg*dAlambda); 1617 plumeheat=mantleconductivity*dt; 1618 crustheat=uppercrustheat*uppercrustthickness+lowercrustheat*(crustthickness-uppercrustthickness); 1619 values[i]=crustheat+plumeheat; 1620 } 1621 1622 this->AddInput(BasalforcingsGeothermalfluxEnum,values,P1Enum); 1623 xDelete<IssmDouble>(xyz_list); 1624 xDelete<IssmDouble>(values); 1625 1626 }/*}}}*/ 1627 void Element::MarshallElement(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction,int numanalyses){/*{{{*/ 1628 1629 _assert_(this); 1630 if(marshall_direction==MARSHALLING_BACKWARD){ 1631 inputs=new Inputs(); 1632 nodes = NULL; 1633 } 1634 1635 MARSHALLING_ENUM(ElementEnum); 1636 1637 MARSHALLING(id); 1638 MARSHALLING(sid); 1639 MARSHALLING(element_type); 1640 MARSHALLING_DYNAMIC(element_type_list,int,numanalyses); 1641 inputs->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction); 1642 1643 } 1644 /*}}}*/ 1265 1645 void Element::MigrateGroundingLine(IssmDouble* phi_ungrounding){/*{{{*/ 1266 1646 … … 1275 1655 IssmDouble* b = xNew<IssmDouble>(numvertices); 1276 1656 IssmDouble* r = xNew<IssmDouble>(numvertices); 1657 IssmDouble* sl = xNew<IssmDouble>(numvertices); 1277 1658 1278 1659 /*Recover info at the vertices: */ … … 1283 1664 GetInputListOnVertices(&b[0],BaseEnum); 1284 1665 GetInputListOnVertices(&r[0],BedEnum); 1666 GetInputListOnVertices(&sl[0],SealevelEnum); 1285 1667 GetInputListOnVertices(&phi[0],MaskGroundediceLevelsetEnum); 1286 1668 rho_water = matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum); … … 1306 1688 /*Change only if AggressiveMigration or if the ice sheet is in contact with the ocean*/ 1307 1689 else{ // phi>0 1308 bed_hydro=-density*h[i] ;1690 bed_hydro=-density*h[i]+sl[i]; 1309 1691 if (bed_hydro>r[i]){ 1310 1692 /*Unground only if the element is connected to the ice shelf*/ 1311 1693 if(migration_style==AggressiveMigrationEnum || migration_style==SubelementMigrationEnum || migration_style==SubelementMigration2Enum){ 1312 s[i] = (1-density)*h[i] ;1313 b[i] = -density*h[i] ;1694 s[i] = (1-density)*h[i]+sl[i]; 1695 b[i] = -density*h[i]+sl[i]; 1314 1696 } 1315 1697 else if(migration_style==SoftMigrationEnum && phi_ungrounding[vertices[i]->Pid()]<0.){ 1316 s[i] = (1-density)*h[i] ;1317 b[i] = -density*h[i] ;1698 s[i] = (1-density)*h[i]+sl[i]; 1699 b[i] = -density*h[i]+sl[i]; 1318 1700 } 1319 1701 else{ … … 1327 1709 for(i=0;i<numvertices;i++){ 1328 1710 if(migration_style==SoftMigrationEnum){ 1329 bed_hydro=-density*h[i] ;1711 bed_hydro=-density*h[i]+sl[i]; 1330 1712 if(phi[i]<0. || bed_hydro<=r[i] || phi_ungrounding[vertices[i]->Pid()]<0.){ 1331 phi[i]=h[i]+ r[i]/density;1713 phi[i]=h[i]+(r[i]-sl[i])/density; 1332 1714 } 1333 1715 } 1334 else if(migration_style!=ContactEnum) phi[i]=h[i]+ r[i]/density;1716 else if(migration_style!=ContactEnum) phi[i]=h[i]+(r[i]-sl[i])/density; 1335 1717 else{ 1336 1718 /*do nothing*/ … … 1349 1731 xDelete<IssmDouble>(b); 1350 1732 xDelete<IssmDouble>(s); 1733 xDelete<IssmDouble>(sl); 1351 1734 xDelete<IssmDouble>(h); 1352 1735 1353 1736 } 1354 1737 /*}}}*/ 1738 void Element::MismipFloatingiceMeltingRate(){/*{{{*/ 1739 1740 int numvertices = this->GetNumberOfVertices(); 1741 IssmDouble meltratefactor,thresholdthickness,upperdepthmelt; 1742 IssmDouble* base = xNew<IssmDouble>(numvertices); 1743 IssmDouble* bed = xNew<IssmDouble>(numvertices); 1744 IssmDouble* values = xNew<IssmDouble>(numvertices); 1745 1746 parameters->FindParam(&meltratefactor,BasalforcingsMeltrateFactorEnum); 1747 parameters->FindParam(&thresholdthickness,BasalforcingsThresholdThicknessEnum); 1748 parameters->FindParam(&upperdepthmelt,BasalforcingsUpperdepthMeltEnum); 1749 1750 this->GetInputListOnVertices(base,BaseEnum); 1751 this->GetInputListOnVertices(bed,BedEnum); 1752 for(int i=0;i<numvertices;i++){ 1753 if(base[i]>upperdepthmelt){ 1754 values[i]=0; 1755 } 1756 else{ 1757 values[i]=meltratefactor*tanh((base[i]-bed[i])/thresholdthickness)*(upperdepthmelt-base[i]); 1758 } 1759 } 1760 1761 this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,values,P1Enum); 1762 xDelete<IssmDouble>(base); 1763 xDelete<IssmDouble>(bed); 1764 xDelete<IssmDouble>(values); 1765 1766 }/*}}}*/ 1355 1767 ElementMatrix* Element::NewElementMatrix(int approximation_enum){/*{{{*/ 1356 1768 return new ElementMatrix(nodes,this->GetNumberOfNodes(),this->parameters,approximation_enum); … … 1365 1777 } 1366 1778 /*}}}*/ 1779 void Element::PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm){/*{{{*/ 1780 1781 int numvertices = this->GetNumberOfVertices(); 1782 1783 int i; 1784 IssmDouble* agd=xNew<IssmDouble>(numvertices); // surface mass balance 1785 IssmDouble* melt=xNew<IssmDouble>(numvertices); // surface mass balance 1786 IssmDouble* accu=xNew<IssmDouble>(numvertices); // surface mass balance 1787 IssmDouble* monthlytemperatures=xNew<IssmDouble>(12*numvertices); 1788 IssmDouble* monthlyprec=xNew<IssmDouble>(12*numvertices); 1789 IssmDouble* yearlytemperatures=xNew<IssmDouble>(numvertices); memset(yearlytemperatures, 0., numvertices*sizeof(IssmDouble)); 1790 IssmDouble* tmp=xNew<IssmDouble>(numvertices); 1791 IssmDouble* h=xNew<IssmDouble>(numvertices); 1792 IssmDouble* s=xNew<IssmDouble>(numvertices); 1793 IssmDouble* s0p=xNew<IssmDouble>(numvertices); 1794 IssmDouble* s0t=xNew<IssmDouble>(numvertices); 1795 IssmDouble rho_water,rho_ice,desfac,rlaps,rlapslgm; 1796 IssmDouble PfacTime,TdiffTime,sealevTime; 1797 IssmDouble mavg=1./12.; //factor for monthly average 1798 1799 /*Get material parameters :*/ 1800 rho_water=this->matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum); 1801 rho_ice=this->matpar->GetMaterialParameter(MaterialsRhoIceEnum); 1802 1803 /*Get some pdd parameters*/ 1804 desfac=this->matpar->GetMaterialParameter(SmbDesfacEnum); 1805 rlaps=this->matpar->GetMaterialParameter(SmbRlapsEnum); 1806 rlapslgm=this->matpar->GetMaterialParameter(SmbRlapslgmEnum); 1807 1808 /*Recover monthly temperatures and precipitation and compute the yearly mean temperatures*/ 1809 Input* input=this->inputs->GetInput(SmbMonthlytemperaturesEnum); _assert_(input); 1810 Input* input2=this->inputs->GetInput(SmbPrecipitationEnum); _assert_(input2); 1811 IssmDouble time,yts,time_yr; 1812 this->parameters->FindParam(&time,TimeEnum); 1813 this->parameters->FindParam(&yts,ConstantsYtsEnum); 1814 time_yr=floor(time/yts)*yts; 1815 1816 /*loop over vertices: */ 1817 Gauss* gauss=this->NewGauss(); 1818 for(int month=0;month<12;month++) { 1819 for(int iv=0;iv<numvertices;iv++) { 1820 gauss->GaussVertex(iv); 1821 input->GetInputValue(&monthlytemperatures[iv*12+month],gauss,time_yr+month/12.*yts); 1822 // yearlytemperatures[iv]=yearlytemperatures[iv]+monthlytemperatures[iv*12+month]*mavg; // Has to be in Kelvin 1823 monthlytemperatures[iv*12+month]=monthlytemperatures[iv*12+month]-273.15; // conversion from Kelvin to celcius for PDD module 1824 input2->GetInputValue(&monthlyprec[iv*12+month],gauss,time_yr+month/12.*yts); 1825 monthlyprec[iv*12+month]=monthlyprec[iv*12+month]*yts; 1826 } 1827 } 1828 1829 /*Recover Pfac, Tdiff and sealev at time t: 1830 * This parameters are used to interpolate the temperature 1831 * and precipitaton between PD and LGM when ismungsm==1 */ 1832 if (ismungsm==1){ 1833 this->parameters->FindParam(&TdiffTime,SmbTdiffEnum,time); 1834 this->parameters->FindParam(&sealevTime,SmbSealevEnum,time); 1835 } 1836 else { 1837 TdiffTime=0; 1838 sealevTime=0; 1839 } 1840 1841 /*Recover info at the vertices: */ 1842 GetInputListOnVertices(&h[0],ThicknessEnum); 1843 GetInputListOnVertices(&s[0],SurfaceEnum); 1844 GetInputListOnVertices(&s0p[0],SmbS0pEnum); 1845 GetInputListOnVertices(&s0t[0],SmbS0tEnum); 1846 1847 /*measure the surface mass balance*/ 1848 for (int iv = 0; iv<numvertices; iv++){ 1849 agd[iv]=PddSurfaceMassBalance(&monthlytemperatures[iv*12], &monthlyprec[iv*12], 1850 pdds, pds, &melt[iv], &accu[iv], signorm, yts, h[iv], s[iv], 1851 desfac, s0t[iv], s0p[iv],rlaps,rlapslgm,TdiffTime,sealevTime, 1852 rho_water,rho_ice); 1853 /*Get yearlytemperatures */ 1854 for(int month=0;month<12;month++) { 1855 yearlytemperatures[iv]=yearlytemperatures[iv]+(monthlytemperatures[iv*12+month]+273.15)*mavg; // Has to be in Kelvin 1856 } 1857 } 1858 1859 /*Update inputs*/ 1860 // TransientInput* NewTemperatureInput = new TransientInput(SmbMonthlytemperaturesEnum); 1861 // TransientInput* NewPrecipitationInput = new TransientInput(SmbPrecipitationEnum); 1862 // for (int imonth=0;imonth<12;imonth++) { 1863 // for(i=0;i<numvertices;i++) tmp[i]=monthlytemperatures[i*12+imonth]; 1864 // TriaInput* newmonthinput1 = new TriaInput(SmbMonthlytemperaturesEnum,&tmp[0],P1Enum); 1865 // NewTemperatureInput->AddTimeInput(newmonthinput1,time+imonth/12.*yts); 1866 // 1867 // for(i=0;i<numvertices;i++) tmp[i]=monthlyprec[i*12+imonth]/yts; 1868 // TriaInput* newmonthinput2 = new TriaInput(SmbPrecipitationEnum,&tmp[0],P1Enum); 1869 // NewPrecipitationInput->AddTimeInput(newmonthinput2,time+imonth/12.*yts); 1870 // } 1871 // NewTemperatureInput->Configure(this->parameters); 1872 // NewPrecipitationInput->Configure(this->parameters); 1873 1874 switch(this->ObjectEnum()){ 1875 case TriaEnum: 1876 // this->inputs->AddInput(new TriaInput(TemperatureEnum,&yearlytemperatures[0],P1Enum)); 1877 this->inputs->AddInput(new TriaInput(TemperaturePDDEnum,&yearlytemperatures[0],P1Enum)); 1878 this->inputs->AddInput(new TriaInput(SmbMassBalanceEnum,&agd[0],P1Enum)); 1879 this->inputs->AddInput(new TriaInput(SmbAccumulationEnum,&accu[0],P1Enum)); 1880 this->inputs->AddInput(new TriaInput(SmbMeltEnum,&melt[0],P1Enum)); 1881 break; 1882 case PentaEnum: 1883 if(IsOnSurface()){ 1884 GetInputListOnVertices(&s[0],TemperatureEnum); 1885 yearlytemperatures[0] = s[0]; 1886 yearlytemperatures[1] = s[1]; 1887 yearlytemperatures[2] = s[2]; 1888 this->inputs->AddInput(new PentaInput(TemperatureEnum,&yearlytemperatures[0],P1Enum)); 1889 } 1890 this->inputs->AddInput(new PentaInput(SmbMassBalanceEnum,&agd[0],P1Enum)); 1891 this->inputs->AddInput(new PentaInput(TemperaturePDDEnum,&yearlytemperatures[0],P1Enum)); 1892 this->InputExtrude(TemperaturePDDEnum,-1); 1893 this->InputExtrude(SmbMassBalanceEnum,-1); 1894 break; 1895 case TetraEnum: 1896 if(IsOnSurface()){ 1897 GetInputListOnVertices(&s[0],TemperatureEnum); 1898 yearlytemperatures[0] = s[0]; 1899 yearlytemperatures[1] = s[1]; 1900 yearlytemperatures[2] = s[2]; 1901 this->inputs->AddInput(new TetraInput(TemperatureEnum,&yearlytemperatures[0],P1Enum)); 1902 } 1903 this->inputs->AddInput(new TetraInput(SmbMassBalanceEnum,&agd[0],P1Enum)); 1904 this->inputs->AddInput(new TetraInput(TemperaturePDDEnum,&yearlytemperatures[0],P1Enum)); 1905 this->InputExtrude(TemperaturePDDEnum,-1); 1906 this->InputExtrude(SmbMassBalanceEnum,-1); 1907 break; 1908 default: _error_("Not implemented yet"); 1909 } 1910 // this->inputs->AddInput(NewTemperatureInput); 1911 // this->inputs->AddInput(NewPrecipitationInput); 1912 // this->inputs->AddInput(new TriaVertexInput(ThermalSpcTemperatureEnum,&Tsurf[0])); 1913 1914 //this->InputExtrude(SmbMassBalanceEnum,-1); 1915 // this->InputExtrude(SmbMonthlytemperaturesEnum,-1); 1916 // this->InputExtrude(SmbPrecipitationEnum,-1); 1917 1918 /*clean-up*/ 1919 delete gauss; 1920 xDelete<IssmDouble>(monthlytemperatures); 1921 xDelete<IssmDouble>(monthlyprec); 1922 xDelete<IssmDouble>(agd); 1923 xDelete<IssmDouble>(melt); 1924 xDelete<IssmDouble>(accu); 1925 xDelete<IssmDouble>(yearlytemperatures); 1926 xDelete<IssmDouble>(h); 1927 xDelete<IssmDouble>(s); 1928 xDelete<IssmDouble>(s0t); 1929 xDelete<IssmDouble>(s0p); 1930 xDelete<IssmDouble>(tmp); 1931 1932 } 1933 /*}}}*/ 1367 1934 IssmDouble Element::PureIceEnthalpy(IssmDouble pressure){/*{{{*/ 1368 1935 return this->matpar->PureIceEnthalpy(pressure); 1369 1936 }/*}}}*/ 1370 void Element::ResultInterpolation(int* pinterpolation,int* pnodesperelement,int output_enum){/*{{{*/ 1371 1937 void Element::ResultInterpolation(int* pinterpolation,int* pnodesperelement,int* parray_size, int output_enum){/*{{{*/ 1938 1939 /*Some intputs need to be computed, even if they are already in inputs, they might not be up to date!*/ 1940 switch(output_enum){ 1941 case ViscousHeatingEnum: this->ViscousHeatingCreateInput(); break; 1942 case StressMaxPrincipalEnum: this->StressMaxPrincipalCreateInput(); break; 1943 case StressTensorxxEnum: 1944 case StressTensorxyEnum: 1945 case StressTensorxzEnum: 1946 case StressTensoryyEnum: 1947 case StressTensoryzEnum: 1948 case StressTensorzzEnum: this->ComputeStressTensor(); break; 1949 case StrainRatexxEnum: 1950 case StrainRatexyEnum: 1951 case StrainRatexzEnum: 1952 case StrainRateyyEnum: 1953 case StrainRateyzEnum: 1954 case StrainRatezzEnum: 1955 case StrainRateeffectiveEnum: this->ComputeStrainRate(); break; 1956 case DeviatoricStressxxEnum: 1957 case DeviatoricStressxyEnum: 1958 case DeviatoricStressxzEnum: 1959 case DeviatoricStressyyEnum: 1960 case DeviatoricStressyzEnum: 1961 case DeviatoricStresszzEnum: 1962 case DeviatoricStresseffectiveEnum: this->ComputeDeviatoricStressTensor(); break; 1963 case SigmaNNEnum: this->ComputeSigmaNN(); break; 1964 case NewDamageEnum: this->ComputeNewDamage(); break; 1965 case StressIntensityFactorEnum: this->StressIntensityFactor(); break; 1966 case CalvingratexEnum: 1967 case CalvingrateyEnum: 1968 case CalvingCalvingrateEnum: 1969 this->StrainRateparallel(); 1970 this->StrainRateperpendicular(); 1971 int calvinglaw; 1972 this->FindParam(&calvinglaw,CalvingLawEnum); 1973 switch(calvinglaw){ 1974 case DefaultCalvingEnum: 1975 //do nothing 1976 break; 1977 case CalvingLevermannEnum: 1978 this->CalvingRateLevermann(); 1979 break; 1980 case CalvingDevEnum: 1981 this->CalvingRateDev(); 1982 break; 1983 default: 1984 _error_("Calving law "<<EnumToStringx(calvinglaw)<<" not supported yet"); 1985 } 1986 break; 1987 case StrainRateparallelEnum: this->StrainRateparallel(); break; 1988 case StrainRateperpendicularEnum: this->StrainRateperpendicular(); break; 1989 } 1990 1991 /*Find input*/ 1372 1992 Input* input=this->inputs->GetInput(output_enum); 1373 1993 1374 1994 /*If this input is not already in Inputs, maybe it needs to be computed?*/ 1375 if(!input){ 1376 switch(output_enum){ 1377 case ViscousHeatingEnum: 1378 this->ViscousHeatingCreateInput(); 1379 input=this->inputs->GetInput(output_enum); 1380 break; 1381 case StressMaxPrincipalEnum: 1382 this->StressMaxPrincipalCreateInput(); 1383 input=this->inputs->GetInput(output_enum); 1384 break; 1385 case StressTensorxxEnum: 1386 case StressTensorxyEnum: 1387 case StressTensorxzEnum: 1388 case StressTensoryyEnum: 1389 case StressTensoryzEnum: 1390 case StressTensorzzEnum: 1391 this->ComputeStressTensor(); 1392 input=this->inputs->GetInput(output_enum); 1393 break; 1394 case StrainRatexxEnum: 1395 case StrainRatexyEnum: 1396 case StrainRatexzEnum: 1397 case StrainRateyyEnum: 1398 case StrainRateyzEnum: 1399 case StrainRatezzEnum: 1400 case StrainRateeffectiveEnum: 1401 this->ComputeStrainRate(); 1402 input=this->inputs->GetInput(output_enum); 1403 break; 1404 case DeviatoricStressxxEnum: 1405 case DeviatoricStressxyEnum: 1406 case DeviatoricStressxzEnum: 1407 case DeviatoricStressyyEnum: 1408 case DeviatoricStressyzEnum: 1409 case DeviatoricStresszzEnum: 1410 this->ComputeDeviatoricStressTensor(); 1411 input=this->inputs->GetInput(output_enum); 1412 break; 1413 case SigmaNNEnum: 1414 this->ComputeSigmaNN(); 1415 input=this->inputs->GetInput(output_enum); 1416 break; 1417 case NewDamageEnum: 1418 this->ComputeNewDamage(); 1419 input=this->inputs->GetInput(output_enum); 1420 break; 1421 case StressIntensityFactorEnum: 1422 this->StressIntensityFactor(); 1423 input=this->inputs->GetInput(output_enum); 1424 break; 1425 case CalvingratexEnum: 1426 case CalvingrateyEnum: 1427 case CalvingCalvingrateEnum: 1428 this->StrainRateparallel(); 1429 this->StrainRateperpendicular(); 1430 int calvinglaw; 1431 this->FindParam(&calvinglaw,CalvingLawEnum); 1432 switch(calvinglaw){ 1433 case DefaultCalvingEnum: 1434 //do nothing 1435 break; 1436 case CalvingLevermannEnum: 1437 this->CalvingRateLevermann(); 1438 break; 1439 case CalvingPiEnum: 1440 this->CalvingRatePi(); 1441 break; 1442 default: 1443 _error_("Calving law "<<EnumToStringx(calvinglaw)<<" not supported yet"); 1444 } 1445 input=this->inputs->GetInput(output_enum); 1446 break; 1447 case StrainRateparallelEnum: 1448 this->StrainRateparallel(); 1449 input=this->inputs->GetInput(output_enum); 1450 break; 1451 case StrainRateperpendicularEnum: 1452 this->StrainRateperpendicular(); 1453 input=this->inputs->GetInput(output_enum); 1454 break; 1455 default: 1456 _error_("input "<<EnumToStringx(output_enum)<<" not found in element"); 1457 } 1458 } 1995 if(!input) _error_("input "<<EnumToStringx(output_enum)<<" not found in element"); 1459 1996 1460 1997 /*Assign output pointer*/ 1461 _assert_(input);1462 1998 *pinterpolation = input->GetResultInterpolation(); 1463 1999 *pnodesperelement = input->GetResultNumberOfNodes(); 2000 *parray_size = input->GetResultArraySize(); 1464 2001 }/*}}}*/ 1465 2002 void Element::ResultToPatch(IssmDouble* values,int nodesperelement,int output_enum){/*{{{*/ … … 1469 2006 1470 2007 input->ResultToPatch(values,nodesperelement,this->Sid()); 2008 2009 } /*}}}*/ 2010 void Element::ResultToMatrix(IssmDouble* values,int ncols,int output_enum){/*{{{*/ 2011 2012 Input* input=this->inputs->GetInput(output_enum); 2013 if(!input) _error_("input "<<EnumToStringx(output_enum)<<" not found in element"); 2014 2015 input->ResultToMatrix(values,ncols,this->Sid()); 1471 2016 1472 2017 } /*}}}*/ … … 1595 2140 } 1596 2141 /*}}}*/ 2142 void Element::SmbGemb(){/*{{{*/ 2143 2144 /*Intermediary variables: {{{*/ 2145 bool isinitialized=false; 2146 IssmDouble zTop,dzTop,zMax,zMin,zY,dzMin; 2147 IssmDouble Tmean; 2148 IssmDouble C; 2149 IssmDouble Tz,Vz; 2150 IssmDouble rho_ice, rho_water,aSnow,aIce; 2151 IssmDouble time,dt; 2152 IssmDouble t,smb_dt; 2153 IssmDouble yts; 2154 IssmDouble Ta,V,dlw,dsw,P,eAir,pAir; 2155 int aIdx=0; 2156 int denIdx=0; 2157 int swIdx=0; 2158 IssmDouble cldFrac,t0wet, t0dry, K; 2159 IssmDouble ulw; 2160 IssmDouble netSW; 2161 IssmDouble netLW; 2162 IssmDouble lhf, shf, dayEC; 2163 IssmDouble initMass; 2164 IssmDouble sumR, sumM, sumEC, sumP, sumW,sumMassAdd; 2165 IssmDouble sumMass,dMass; 2166 bool isgraingrowth,isalbedo,isshortwave,isthermal,isaccumulation,ismelt,isdensification,isturbulentflux; 2167 IssmDouble init_scaling=1.0; 2168 /*}}}*/ 2169 /*Output variables:{{{ */ 2170 IssmDouble* dz=NULL; 2171 IssmDouble* d = NULL; 2172 IssmDouble* re = NULL; 2173 IssmDouble* gdn = NULL; 2174 IssmDouble* gsp = NULL; 2175 IssmDouble EC = 0; 2176 IssmDouble* W = NULL; 2177 IssmDouble* a = NULL; 2178 IssmDouble* swf=NULL; 2179 IssmDouble* T = NULL; 2180 IssmDouble T_bottom; 2181 IssmDouble M; 2182 IssmDouble R; 2183 IssmDouble mAdd; 2184 int m; 2185 int count=0; 2186 /*}}}*/ 2187 2188 /*only compute SMB at the surface: */ 2189 if (!IsOnSurface()) return; 2190 2191 2192 /*Retrieve material properties and parameters:{{{ */ 2193 rho_ice = matpar->GetMaterialParameter(MaterialsRhoIceEnum); 2194 rho_water = matpar->GetMaterialParameter(MaterialsRhoFreshwaterEnum); 2195 parameters->FindParam(&aSnow,SmbASnowEnum); 2196 parameters->FindParam(&aIce,SmbAIceEnum); 2197 parameters->FindParam(&time,TimeEnum); /*transient core time at which we run the smb core*/ 2198 parameters->FindParam(&dt,TimesteppingTimeStepEnum); /*transient core time step*/ 2199 parameters->FindParam(&smb_dt,SmbDtEnum); /*time period for the smb solution, usually smaller than the glaciological dt*/ 2200 parameters->FindParam(&aIdx,SmbAIdxEnum); 2201 parameters->FindParam(&denIdx,SmbDenIdxEnum); 2202 parameters->FindParam(&swIdx,SmbSwIdxEnum); 2203 parameters->FindParam(&cldFrac,SmbCldFracEnum); 2204 parameters->FindParam(&t0wet,SmbT0wetEnum); 2205 parameters->FindParam(&t0dry,SmbT0dryEnum); 2206 parameters->FindParam(&K,SmbKEnum); 2207 parameters->FindParam(&isgraingrowth,SmbIsgraingrowthEnum); 2208 parameters->FindParam(&isalbedo,SmbIsalbedoEnum); 2209 parameters->FindParam(&isshortwave,SmbIsshortwaveEnum); 2210 parameters->FindParam(&isthermal,SmbIsthermalEnum); 2211 parameters->FindParam(&isaccumulation,SmbIsaccumulationEnum); 2212 parameters->FindParam(&ismelt,SmbIsmeltEnum); 2213 parameters->FindParam(&isdensification,SmbIsdensificationEnum); 2214 parameters->FindParam(&isturbulentflux,SmbIsturbulentfluxEnum); 2215 parameters->FindParam(&init_scaling,SmbInitDensityScalingEnum); 2216 /*}}}*/ 2217 /*Retrieve inputs: {{{*/ 2218 Input* zTop_input=this->GetInput(SmbZTopEnum); _assert_(zTop_input); 2219 Input* dzTop_input=this->GetInput(SmbDzTopEnum); _assert_(dzTop_input); 2220 Input* dzMin_input=this->GetInput(SmbDzMinEnum); _assert_(dzMin_input); 2221 Input* zMax_input=this->GetInput(SmbZMaxEnum); _assert_(zMax_input); 2222 Input* zMin_input=this->GetInput(SmbZMinEnum); _assert_(zMin_input); 2223 Input* zY_input=this->GetInput(SmbZYEnum); _assert_(zY_input); 2224 Input* Tmean_input=this->GetInput(SmbTmeanEnum); _assert_(Tmean_input); 2225 Input* C_input=this->GetInput(SmbCEnum); _assert_(C_input); 2226 Input* Tz_input=this->GetInput(SmbTzEnum); _assert_(Tz_input); 2227 Input* Vz_input=this->GetInput(SmbVzEnum); _assert_(Vz_input); 2228 Input* Ta_input=this->GetInput(SmbTaEnum); _assert_(Ta_input); 2229 Input* V_input=this->GetInput(SmbVEnum); _assert_(V_input); 2230 Input* Dlwr_input=this->GetInput(SmbDlwrfEnum); _assert_(Dlwr_input); 2231 Input* Dswr_input=this->GetInput(SmbDswrfEnum); _assert_(Dswr_input); 2232 Input* P_input=this->GetInput(SmbPEnum); _assert_(P_input); 2233 Input* eAir_input=this->GetInput(SmbEAirEnum); _assert_(eAir_input); 2234 Input* pAir_input=this->GetInput(SmbPAirEnum); _assert_(pAir_input); 2235 Input* isinitialized_input=this->inputs->GetInput(SmbIsInitializedEnum); 2236 2237 /*Retrieve input values:*/ 2238 Gauss* gauss=this->NewGauss(1); gauss->GaussPoint(0); 2239 2240 zTop_input->GetInputValue(&zTop,gauss); 2241 dzTop_input->GetInputValue(&dzTop,gauss); 2242 dzMin_input->GetInputValue(&dzMin,gauss); 2243 zMax_input->GetInputValue(&zMax,gauss); 2244 zMin_input->GetInputValue(&zMin,gauss); 2245 zY_input->GetInputValue(&zY,gauss); 2246 Tmean_input->GetInputValue(&Tmean,gauss); 2247 C_input->GetInputValue(&C,gauss); 2248 Tz_input->GetInputValue(&Tz,gauss); 2249 Vz_input->GetInputValue(&Vz,gauss); 2250 /*}}}*/ 2251 2252 /*First, check that the initial structures have been setup in GEMB. If not, initialize profile variables: layer thickness dz, * density d, temperature T, etc. {{{*/ 2253 if(!isinitialized_input){ 2254 2255 if(VerboseSmb() && this->Sid()==0)_printf0_("smb core: Initializing grid\n"); 2256 GembgridInitialize(&dz, &m, zTop, dzTop, zMax, zY); 2257 //if(this->Sid()==1) for(int i=0;i<m;i++)_printf_("z[" << i << "]=" << 2258 //dz[i] << "\n"); 2259 2260 /*initialize profile variables:*/ 2261 d = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)d[i]=rho_ice*init_scaling; //ice density scaled by a factor 2262 re = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)re[i]=2.5; //set grain size to old snow [mm] 2263 gdn = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)gdn[i]=0; //set grain dentricity to old snow 2264 gsp = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)gsp[i]=0; //set grain sphericity to old snow 2265 EC = 0; //surface evaporation (-) condensation (+) [kg m-2] 2266 W = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)W[i]=0; //set water content to zero [kg m-2] 2267 a = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)a[i]=aSnow; //set albedo equal to fresh snow [fraction] 2268 T = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)T[i]=Tmean; //set initial grid cell temperature to the annual mean temperature [K] 2269 2270 //fixed lower temperatuer bounday condition - T is fixed 2271 T_bottom=T[m-1]; 2272 2273 /*Flag the initialization:*/ 2274 this->AddInput(new BoolInput(SmbIsInitializedEnum,true)); 2275 } 2276 else{ 2277 /*Recover inputs: */ 2278 DoubleArrayInput* dz_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbDzEnum)); _assert_(dz_input); 2279 DoubleArrayInput* d_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbDEnum));_assert_(d_input); 2280 DoubleArrayInput* re_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbReEnum));_assert_(re_input); 2281 DoubleArrayInput* gdn_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbGdnEnum));_assert_(gdn_input); 2282 DoubleArrayInput* gsp_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbGspEnum));_assert_(gsp_input); 2283 DoubleInput* EC_input= dynamic_cast<DoubleInput*>(this->GetInput(SmbECEnum));_assert_(EC_input); 2284 DoubleArrayInput* W_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbWEnum));_assert_(W_input); 2285 DoubleArrayInput* a_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbAEnum));_assert_(a_input); 2286 DoubleArrayInput* T_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbTEnum));_assert_(T_input); 2287 2288 /*Recover arrays: */ 2289 dz_input->GetValues(&dz,&m); 2290 d_input->GetValues(&d,&m); 2291 re_input->GetValues(&re,&m); 2292 gdn_input->GetValues(&gdn,&m); 2293 gsp_input->GetValues(&gsp,&m); 2294 EC_input->GetInputValue(&EC); 2295 W_input->GetValues(&W,&m); 2296 a_input->GetValues(&a,&m); 2297 T_input->GetValues(&T,&m); 2298 2299 //fixed lower temperatuer bounday condition - T is fixed 2300 T_bottom=T[m-1]; 2301 2302 } /*}}}*/ 2303 2304 // determine initial mass [kg] 2305 initMass=0; for(int i=0;i<m;i++) initMass += dz[i]*d[i] + W[i]; 2306 2307 // initialize cumulative variables 2308 sumR = 0; sumM = 0; sumEC = 0; sumP = 0; sumMassAdd = 0; 2309 2310 //before starting loop, realize that the transient core runs this smb_core at time = time +deltaT. 2311 //go back to time - deltaT: 2312 time-=dt; 2313 2314 /*Start loop: */ 2315 count=1; 2316 for (t=time;t<=time+dt;t=t+smb_dt){ 2317 2318 if(VerboseSmb() && this->Sid()==0 && IssmComm::GetRank()==0)_printf0_("Time: t=" << setprecision(8) << t/365.0/24.0/3600.0 << " yr/" << (time+dt)/365.0/24.0/3600.0 << " yr" << setprecision(3) << " Step: " << count << "\n"); 2319 2320 /*extract daily data:{{{*/ 2321 Ta_input->GetInputValue(&Ta,gauss,t);//screen level air temperature [K] 2322 V_input->GetInputValue(&V,gauss,t); //wind speed [m s-1] 2323 Dlwr_input->GetInputValue(&dlw,gauss,t); //downward longwave radiation flux [W m-2] 2324 Dswr_input->GetInputValue(&dsw,gauss,t); //downward shortwave radiation flux [W m-2] 2325 P_input->GetInputValue(&P,gauss,t); //precipitation [kg m-2] 2326 eAir_input->GetInputValue(&eAir,gauss,t); //screen level vapor pressure [Pa] 2327 pAir_input->GetInputValue(&pAir,gauss,t); // screen level air pressure [Pa] 2328 //_printf_("Time: " << t << " Ta: " << Ta << " V: " << V << " dlw: " << dlw << " dsw: " << dsw << " P: " << P << " eAir: " << eAir << " pAir: " << pAir << "\n"); 2329 /*}}}*/ 2330 2331 /*Snow grain metamorphism:*/ 2332 if(isgraingrowth)grainGrowth(re, gdn, gsp, T, dz, d, W, smb_dt, m, aIdx,this->Sid()); 2333 2334 /*Snow, firn and ice albedo:*/ 2335 if(isalbedo)albedo(a,aIdx,re,d,cldFrac,aIce, aSnow,T,W,P,EC,t0wet,t0dry,K,smb_dt,m,this->Sid()); 2336 2337 2338 /*Distribution of absorbed short wave radation with depth:*/ 2339 if(isshortwave)shortwave(&swf, swIdx, aIdx, dsw, a[0], d, dz, re,m,this->Sid()); 2340 2341 /*Calculate net shortwave [W m-2]*/ 2342 netSW = cellsum(swf,m); 2343 2344 /*Thermal profile computation:*/ 2345 if(isthermal)thermo(&EC, T, dz, d, swf, dlw, Ta, V, eAir, pAir, W[0], smb_dt, m, Vz, Tz,this->Sid()); 2346 2347 /*Change in thickness of top cell due to evaporation/condensation assuming same density as top cell. 2348 * need to fix this in case all or more of cell evaporates */ 2349 dz[0] = dz[0] + EC / d[0]; 2350 2351 /*Add snow/rain to top grid cell adjusting cell depth, temperature and density*/ 2352 if(isaccumulation)accumulation(&T, &dz, &d, &W, &a, &re, &gdn, &gsp, &m, Ta, P, dzMin, aSnow,this->Sid()); 2353 2354 /*Calculate water production, M [kg m-2] resulting from snow/ice temperature exceeding 273.15 deg K 2355 * (> 0 deg C), runoff R [kg m-2] and resulting changes in density and determine wet compaction [m]*/ 2356 if(ismelt)melt(&M, &R, &mAdd, &T, &d, &dz, &W, &a, &re, &gdn, &gsp, &m, dzMin, zMax, zMin,this->Sid()); 2357 2358 /*Allow non-melt densification and determine compaction [m]*/ 2359 if(isdensification)densification(d,dz, T, re, denIdx, C, smb_dt, Tmean,rho_ice,m,this->Sid()); 2360 2361 /*Calculate upward longwave radiation flux [W m-2] not used in energy balance. Calculated for every 2362 * sub-time step in thermo equations*/ 2363 ulw = 5.67E-8 * pow(T[0],4.0); 2364 2365 /*Calculate net longwave [W m-2]*/ 2366 netLW = dlw - ulw; 2367 2368 /*Calculate turbulent heat fluxes [W m-2]*/ 2369 if(isturbulentflux)turbulentFlux(&shf, &lhf, &dayEC, Ta, T[0], V, eAir, pAir, d[0], W[0], Vz, Tz,this->Sid()); 2370 2371 /*Verbose some resuls in debug mode: {{{*/ 2372 if(VerboseSmb() && 0){ 2373 _printf_("smb log: count[" << count << "] m[" << m << "] " 2374 << setprecision(16) << "T[" << cellsum(T,m) << "] " 2375 << "d[" << cellsum(d,m) << "] " 2376 << "dz[" << cellsum(dz,m) << "] " 2377 << "a[" << cellsum(a,m) << "] " 2378 << "W[" << cellsum(W,m) << "] " 2379 << "re[" << cellsum(re,m) << "] " 2380 << "gdn[" << cellsum(gdn,m) << "] " 2381 << "gsp[" << cellsum(gsp,m) << "] " 2382 << "swf[" << netSW << "] " 2383 << "\n"); 2384 } /*}}}*/ 2385 2386 /*Sum component mass changes [kg m-2]*/ 2387 sumMassAdd = mAdd + sumMassAdd; 2388 sumM = M + sumM; 2389 sumR = R + sumR; 2390 sumW = cellsum(W,m); 2391 sumP = P + sumP; 2392 sumEC = sumEC + EC; // evap (-)/cond(+) 2393 2394 /*Calculate total system mass:*/ 2395 sumMass=0; for(int i=0;i<m;i++) sumMass += dz[i]*d[i]; 2396 2397 #ifndef _HAVE_ADOLC_ //we want to avoid the round operation at all cost. Not differentiable. 2398 dMass = sumMass + sumR + sumW - sumP - sumEC - initMass - sumMassAdd; 2399 dMass = round(dMass * 100.0)/100.0; 2400 2401 /*Check mass conservation:*/ 2402 if (dMass != 0.0) _printf_("total system mass not conserved in MB function"); 2403 #endif 2404 2405 /*Check bottom grid cell T is unchanged:*/ 2406 if (T[m-1]!=T_bottom) _printf_("T(end)~=T_bottom" << "\n"); 2407 2408 /*Free ressources: */ 2409 xDelete<IssmDouble>(swf); 2410 2411 /*increase counter:*/ 2412 count++; 2413 2414 } //for (t=time;t<=time+dt;t=t+smb_dt) 2415 2416 2417 /*Save generated inputs: */ 2418 this->AddInput(new DoubleArrayInput(SmbDzEnum,dz,m)); 2419 this->AddInput(new DoubleArrayInput(SmbDEnum,d,m)); 2420 this->AddInput(new DoubleArrayInput(SmbReEnum,re,m)); 2421 this->AddInput(new DoubleArrayInput(SmbGdnEnum,gdn,m)); 2422 this->AddInput(new DoubleArrayInput(SmbGspEnum,gsp,m)); 2423 this->AddInput(new DoubleArrayInput(SmbTEnum,T,m)); 2424 this->AddInput(new DoubleInput(SmbECEnum,EC)); 2425 this->AddInput(new DoubleArrayInput(SmbWEnum,W,m)); 2426 this->AddInput(new DoubleArrayInput(SmbAEnum,a,m)); 2427 this->AddInput(new DoubleInput(SmbMassBalanceEnum,(sumP + sumEC -sumR)/rho_water/dt)); 2428 this->AddInput(new DoubleInput(SmbRunoffEnum,sumR/rho_water/dt)); 2429 this->AddInput(new DoubleInput(SmbPrecipitationEnum,sumP/rho_water/dt)); 2430 this->AddInput(new DoubleInput(SmbCondensationEnum,sumEC/rho_water/dt)); 2431 2432 2433 2434 /*Free allocations:{{{*/ 2435 xDelete<IssmDouble>(dz); 2436 xDelete<IssmDouble>(d); 2437 xDelete<IssmDouble>(re); 2438 xDelete<IssmDouble>(gdn); 2439 xDelete<IssmDouble>(gsp); 2440 xDelete<IssmDouble>(W); 2441 xDelete<IssmDouble>(a); 2442 xDelete<IssmDouble>(T); 2443 delete gauss; 2444 /*}}}*/ 2445 } 2446 /*}}}*/ 1597 2447 void Element::StrainRateFS(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/ 1598 2448 /*Compute the 3d Strain Rate (6 components): … … 1794 2644 /*Clean up and return*/ 1795 2645 xDelete<IssmDouble>(maxprincipal); 2646 xDelete<IssmDouble>(xyz_list); 1796 2647 delete gauss; 1797 2648 } … … 2223 3074 /*Clean up and return*/ 2224 3075 xDelete<IssmDouble>(viscousheating); 3076 xDelete<IssmDouble>(xyz_list); 2225 3077 delete gauss; 2226 3078 }
Note:
See TracChangeset
for help on using the changeset viewer.
