Changeset 26947
- Timestamp:
- 04/09/22 06:47:21 (3 years ago)
- Location:
- issm/trunk-jpl
- Files:
-
- 14 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/classes/Elements/Element.cpp
r26850 r26947 2296 2296 2297 2297 }/*}}}*/ 2298 void Element::LapseRateBasinSMB( IssmDouble* lapseratepos, IssmDouble* lapserateneg,IssmDouble* refelevation){/*{{{*/2298 void Element::LapseRateBasinSMB(int numelevbins, IssmDouble* lapserates, IssmDouble* elevbins,IssmDouble* refelevation){/*{{{*/ 2299 2299 2300 2300 /*Variable declaration*/ 2301 2301 const int numvertices = this->GetNumberOfVertices(); 2302 bool isadjustsmb = true; 2303 int basinid; 2304 IssmDouble ela,refelevation_b,lapseratepos_b,lapserateneg_b; 2305 IssmDouble* surfacelist = xNew<IssmDouble>(numvertices); 2306 IssmDouble* smblist = xNew<IssmDouble>(numvertices); 2302 bool isadjustsmb = false; 2303 int basinid,bb1,bb2; 2304 IssmDouble ela,refelevation_b; 2305 IssmDouble* surfacelist = xNew<IssmDouble>(numvertices); 2306 IssmDouble* smblist = xNew<IssmDouble>(numvertices); 2307 /* numelevbins values of lapse rates */ 2308 IssmDouble* lapserates_b = xNew<IssmDouble>(numelevbins); 2309 /* (numelevbins-1) limits between elevation bins (be cautious with indexing) */ 2310 IssmDouble* elevbins_b = xNew<IssmDouble>(numelevbins-1); 2307 2311 2308 2312 /*Retrieve SMB values non-adjusted for SMB lapse rate*/ 2309 2313 this->GetInputListOnVertices(smblist,SmbMassBalanceEnum); 2314 /*Get surface elevation on vertices*/ 2315 this->GetInputListOnVertices(surfacelist,SurfaceEnum); 2310 2316 /*Get basin-specific parameters of the element*/ 2311 2317 this->GetInputValue(&basinid,SmbBasinsIdEnum); 2312 2318 refelevation_b = refelevation[basinid]; 2313 lapseratepos_b = lapseratepos[basinid]; 2314 lapserateneg_b = lapserateneg[basinid]; 2315 /*Get surface elevation on vertices*/ 2316 this->GetInputListOnVertices(surfacelist,SurfaceEnum); 2317 2318 for(int v=0;v<numvertices;v++){ 2319 /*Find ELA*/ 2320 if(smblist[v]>=0 && lapseratepos_b!=0) ela = refelevation_b-smblist[v]/lapseratepos_b; 2321 else if(smblist[v]<0 && lapserateneg_b!=0) ela = refelevation_b-smblist[v]/lapserateneg_b; 2322 /*Lapse rate is zero: SMB not adjusted*/ 2323 else isadjustsmb = false; 2324 /*Adjust SMB value*/ 2325 if(isadjustsmb) smblist[v] = lapseratepos_b*max(surfacelist[v]-ela,0.0)+lapserateneg_b*min(surfacelist[v]-ela,0.0); 2319 for(int ii=0;ii<(numelevbins-1);ii++) elevbins_b[ii] = elevbins[basinid*(numelevbins-1)+ii]; 2320 for(int ii=0;ii<numelevbins;ii++){ 2321 lapserates_b[ii] = lapserates[basinid*numelevbins+ii]; 2322 if(lapserates_b[ii]!=0) isadjustsmb=true; 2323 } 2324 /*Adjust SMB if any lapse rate value is non-zero*/ 2325 if(isadjustsmb){ 2326 2327 for(int v=0;v<numvertices;v++){ 2328 /*Find elevation bin of Reference elevation and of Vertex*/ 2329 for(int ii=0;ii<(numelevbins-1);ii++){ 2330 if(surfacelist[v]<=elevbins_b[ii]) bb1 = ii; 2331 if(refelevation_b<=elevbins_b[ii]) bb2 = ii; 2332 } 2333 /*Check for elevations above highest bin limit */ 2334 if(surfacelist[v]>elevbins_b[numelevbins-1-1]) bb1 = numelevbins-1; 2335 if(refelevation_b>elevbins_b[numelevbins-1-1]) bb2 = numelevbins-1; 2336 /*Vertex and Reference elevation in same elevation bin*/ 2337 if(bb1==bb2){ 2338 smblist[v] = smblist[v]+(surfacelist[v]-refelevation_b)*lapserates_b[bb2]; 2339 } 2340 /*Vertex in lower elevation bin than Reference elevation*/ 2341 if(bb1<bb2){ 2342 smblist[v] = smblist[v]+(elevbins_b[bb2-1]-refelevation_b)*lapserates_b[bb2]; 2343 for(int ii=bb2-1;ii>bb1;ii--) smblist[v] = smblist[v]+(elevbins_b[ii-1]-elevbins_b[ii])*lapserates_b[ii]; 2344 smblist[v] = smblist[v]+(surfacelist[v]-elevbins_b[bb1])*lapserates_b[bb1]; 2345 } 2346 /*Vertex in higher elevation bin than Reference elevation*/ 2347 if(bb1>bb2){ 2348 smblist[v] = smblist[v]+(elevbins_b[bb2]-refelevation_b)*lapserates_b[bb2]; 2349 for(int ii=bb2+1;ii<bb1;ii++) smblist[v] = smblist[v]+(elevbins_b[ii]-elevbins_b[ii-1])*lapserates_b[ii]; 2350 smblist[v] = smblist[v]+(surfacelist[v]-elevbins_b[bb1-1])*lapserates_b[bb1]; 2351 } 2352 } 2326 2353 } 2327 2354 … … 2330 2357 2331 2358 /*Cleanup*/ 2359 xDelete<IssmDouble>(lapserates_b); 2360 xDelete<IssmDouble>(elevbins_b); 2332 2361 xDelete<IssmDouble>(surfacelist); 2333 2362 xDelete<IssmDouble>(smblist); -
issm/trunk-jpl/src/c/classes/Elements/Element.h
r26850 r26947 154 154 bool IsLandInElement(); 155 155 void Ismip6FloatingiceMeltingRate(); 156 void LapseRateBasinSMB( IssmDouble* lapseratepos, IssmDouble* lapserateneg,IssmDouble* refelevation);156 void LapseRateBasinSMB(int numelevbins, IssmDouble* lapserates, IssmDouble* elevbins,IssmDouble* refelevation); 157 157 void LinearFloatingiceMeltingRate(); 158 158 void SpatialLinearFloatingiceMeltingRate(); … … 229 229 virtual void AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part)=0; 230 230 virtual void BasalNodeIndices(int* pnumindices,int** pindices,int finiteelement){_error_("not implemented yet");}; 231 virtual void CalvingProjectionXY(void){_error_("not implemented yet");}; 231 232 virtual void CalvingRateParameterization(void){_error_("not implemented yet");}; 232 233 virtual void CalvingRateVonmises(void){_error_("not implemented yet");}; -
issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp
r26859 r26947 428 428 parameters->AddObject(iomodel->CopyConstantObject("md.smb.ar_initialtime",SmbAutoregressionInitialTimeEnum)); 429 429 parameters->AddObject(iomodel->CopyConstantObject("md.smb.ar_timestep",SmbAutoregressionTimestepEnum)); 430 parameters->AddObject(iomodel->CopyConstantObject("md.smb.num_bins",SmbNumElevationBinsEnum)); 430 431 iomodel->FetchData(&transparam,&M,&N,"md.smb.beta0"); 431 432 parameters->AddObject(new DoubleVecParam(SmbBeta0Enum,transparam,N)); … … 437 438 parameters->AddObject(new DoubleMatParam(SmbPhiEnum,transparam,M,N)); 438 439 xDelete<IssmDouble>(transparam); 439 iomodel->FetchData(&transparam,&M,&N,"md.smb.lapserate _pos");440 parameters->AddObject(new Double VecParam(SmbLapseRatePosEnum,transparam,N));441 xDelete<IssmDouble>(transparam); 442 iomodel->FetchData(&transparam,&M,&N,"md.smb. lapserate_neg");443 parameters->AddObject(new Double VecParam(SmbLapseRateNegEnum,transparam,N));440 iomodel->FetchData(&transparam,&M,&N,"md.smb.lapserates"); 441 parameters->AddObject(new DoubleMatParam(SmbLapseRatesEnum,transparam,M,N)); 442 xDelete<IssmDouble>(transparam); 443 iomodel->FetchData(&transparam,&M,&N,"md.smb.elevationbins"); 444 parameters->AddObject(new DoubleMatParam(SmbElevationBinsEnum,transparam,M,N)); 444 445 xDelete<IssmDouble>(transparam); 445 446 iomodel->FetchData(&transparam,&M,&N,"md.smb.refelevation"); -
issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp
r26810 r26947 198 198 bool isstochastic; 199 199 bool issmbstochastic = false; 200 int M,N,Nphi,arorder,numbasins, my_rank;200 int M,N,Nphi,arorder,numbasins,numelevbins,my_rank; 201 201 femmodel->parameters->FindParam(&numbasins,SmbNumBasinsEnum); 202 202 femmodel->parameters->FindParam(&arorder,SmbAutoregressiveOrderEnum); 203 femmodel->parameters->FindParam(&numelevbins,SmbNumElevationBinsEnum); 203 204 IssmDouble tinit_ar; 204 IssmDouble* beta0 = NULL;205 IssmDouble* beta1 = NULL;206 IssmDouble* phi = NULL;207 IssmDouble* lapserate pos= NULL;208 IssmDouble* lapserateneg= NULL;209 IssmDouble* refelevation = NULL;205 IssmDouble* beta0 = NULL; 206 IssmDouble* beta1 = NULL; 207 IssmDouble* phi = NULL; 208 IssmDouble* lapserates = NULL; 209 IssmDouble* elevbins = NULL; 210 IssmDouble* refelevation = NULL; 210 211 211 212 femmodel->parameters->FindParam(&tinit_ar,SmbAutoregressionInitialTimeEnum); 212 femmodel->parameters->FindParam(&beta0,&M,SmbBeta0Enum); _assert_(M==numbasins);213 femmodel->parameters->FindParam(&beta1,&M,SmbBeta1Enum); _assert_(M==numbasins);214 femmodel->parameters->FindParam(&phi,&M,&Nphi,SmbPhiEnum); _assert_(M==numbasins); _assert_(Nphi==arorder);215 femmodel->parameters->FindParam(&lapserate pos,&M,SmbLapseRatePosEnum); _assert_(M==numbasins);216 femmodel->parameters->FindParam(& lapserateneg,&M,SmbLapseRateNegEnum); _assert_(M==numbasins);217 femmodel->parameters->FindParam(&refelevation,&M,SmbRefElevationEnum); _assert_(M==numbasins);213 femmodel->parameters->FindParam(&beta0,&M,SmbBeta0Enum); _assert_(M==numbasins); 214 femmodel->parameters->FindParam(&beta1,&M,SmbBeta1Enum); _assert_(M==numbasins); 215 femmodel->parameters->FindParam(&phi,&M,&Nphi,SmbPhiEnum); _assert_(M==numbasins); _assert_(Nphi==arorder); 216 femmodel->parameters->FindParam(&lapserates,&M,&N,SmbLapseRatesEnum); _assert_(M==numbasins); _assert_(N==numelevbins); 217 femmodel->parameters->FindParam(&elevbins,&M,&N,SmbElevationBinsEnum); _assert_(M==numbasins); _assert_(N==numelevbins-1); 218 femmodel->parameters->FindParam(&refelevation,&M,SmbRefElevationEnum); _assert_(M==numbasins); 218 219 219 220 femmodel->parameters->FindParam(&isstochastic,StochasticForcingIsStochasticForcingEnum); … … 237 238 element->Autoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,issmbstochastic,SMBautoregressionEnum); 238 239 /*Compute lapse rate adjustment*/ 239 element->LapseRateBasinSMB( lapseratepos,lapserateneg,refelevation);240 element->LapseRateBasinSMB(numelevbins,lapserates,elevbins,refelevation); 240 241 } 241 242 … … 244 245 xDelete<IssmDouble>(beta1); 245 246 xDelete<IssmDouble>(phi); 246 xDelete<IssmDouble>(lapserate pos);247 xDelete<IssmDouble>( lapserateneg);247 xDelete<IssmDouble>(lapserates); 248 xDelete<IssmDouble>(elevbins); 248 249 xDelete<IssmDouble>(refelevation); 249 250 }/*}}}*/ -
issm/trunk-jpl/src/c/shared/Enum/Enum.vim
r26874 r26947 465 465 syn keyword cConstant SmbDpermilEnum 466 466 syn keyword cConstant SmbDsnowIdxEnum 467 syn keyword cConstant SmbElevationBinsEnum 467 468 syn keyword cConstant SmbCldFracEnum 468 469 syn keyword cConstant SmbDelta18oEnum … … 492 493 syn keyword cConstant SmbIsturbulentfluxEnum 493 494 syn keyword cConstant SmbKEnum 494 syn keyword cConstant SmbLapseRateNegEnum 495 syn keyword cConstant SmbLapseRatePosEnum 495 syn keyword cConstant SmbLapseRatesEnum 496 496 syn keyword cConstant SmbNumBasinsEnum 497 syn keyword cConstant SmbNumElevationBinsEnum 497 498 syn keyword cConstant SmbNumRequestedOutputsEnum 498 499 syn keyword cConstant SmbPfacEnum … … 1086 1087 syn keyword cConstant TransientAccumulatedDeltaIceThicknessEnum 1087 1088 syn keyword cConstant VelEnum 1089 syn keyword cConstant VVmismipglposinputEnum 1088 1090 syn keyword cConstant VxAverageEnum 1089 1091 syn keyword cConstant VxBaseEnum … … 1623 1625 syn keyword cType Cfsurfacesquare 1624 1626 syn keyword cType Channel 1627 syn keyword cType classes 1625 1628 syn keyword cType Constraint 1626 1629 syn keyword cType Constraints … … 1629 1632 syn keyword cType ControlInput 1630 1633 syn keyword cType Covertree 1634 syn keyword cType DatasetInput 1631 1635 syn keyword cType DataSetParam 1632 syn keyword cType DatasetInput1633 1636 syn keyword cType Definition 1634 1637 syn keyword cType DependentObject … … 1643 1646 syn keyword cType ElementInput 1644 1647 syn keyword cType ElementMatrix 1648 syn keyword cType Elements 1645 1649 syn keyword cType ElementVector 1646 syn keyword cType Elements1647 1650 syn keyword cType ExponentialVariogram 1648 1651 syn keyword cType ExternalResult … … 1651 1654 syn keyword cType Friction 1652 1655 syn keyword cType Gauss 1656 syn keyword cType GaussianVariogram 1657 syn keyword cType gaussobjects 1653 1658 syn keyword cType GaussPenta 1654 1659 syn keyword cType GaussSeg 1655 1660 syn keyword cType GaussTetra 1656 1661 syn keyword cType GaussTria 1657 syn keyword cType GaussianVariogram1658 1662 syn keyword cType GenericExternalResult 1659 1663 syn keyword cType GenericOption … … 1671 1675 syn keyword cType IssmDirectApplicInterface 1672 1676 syn keyword cType IssmParallelDirectApplicInterface 1677 syn keyword cType krigingobjects 1673 1678 syn keyword cType Load 1674 1679 syn keyword cType Loads … … 1681 1686 syn keyword cType Matice 1682 1687 syn keyword cType Matlitho 1688 syn keyword cType matrixobjects 1683 1689 syn keyword cType MatrixParam 1684 1690 syn keyword cType Misfit … … 1693 1699 syn keyword cType Observations 1694 1700 syn keyword cType Option 1701 syn keyword cType Options 1695 1702 syn keyword cType OptionUtilities 1696 syn keyword cType Options1697 1703 syn keyword cType Param 1698 1704 syn keyword cType Parameters … … 1708 1714 syn keyword cType Regionaloutput 1709 1715 syn keyword cType Results 1716 syn keyword cType Riftfront 1710 1717 syn keyword cType RiftStruct 1711 syn keyword cType Riftfront1712 1718 syn keyword cType SealevelGeometry 1713 1719 syn keyword cType Seg 1714 1720 syn keyword cType SegInput 1721 syn keyword cType Segment 1715 1722 syn keyword cType SegRef 1716 syn keyword cType Segment1717 1723 syn keyword cType SpcDynamic 1718 1724 syn keyword cType SpcStatic … … 1733 1739 syn keyword cType Vertex 1734 1740 syn keyword cType Vertices 1735 syn keyword cType classes1736 syn keyword cType gaussobjects1737 syn keyword cType krigingobjects1738 syn keyword cType matrixobjects1739 1741 syn keyword cType AdjointBalancethickness2Analysis 1740 1742 syn keyword cType AdjointBalancethicknessAnalysis -
issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
r26859 r26947 459 459 SmbDpermilEnum, 460 460 SmbDsnowIdxEnum, 461 SmbElevationBinsEnum, 461 462 SmbCldFracEnum, 462 463 SmbDelta18oEnum, … … 486 487 SmbIsturbulentfluxEnum, 487 488 SmbKEnum, 488 SmbLapseRateNegEnum, 489 SmbLapseRatePosEnum, 489 SmbLapseRatesEnum, 490 490 SmbNumBasinsEnum, 491 SmbNumElevationBinsEnum, 491 492 SmbNumRequestedOutputsEnum, 492 493 SmbPfacEnum, -
issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
r26859 r26947 467 467 case SmbDpermilEnum : return "SmbDpermil"; 468 468 case SmbDsnowIdxEnum : return "SmbDsnowIdx"; 469 case SmbElevationBinsEnum : return "SmbElevationBins"; 469 470 case SmbCldFracEnum : return "SmbCldFrac"; 470 471 case SmbDelta18oEnum : return "SmbDelta18o"; … … 494 495 case SmbIsturbulentfluxEnum : return "SmbIsturbulentflux"; 495 496 case SmbKEnum : return "SmbK"; 496 case SmbLapseRateNegEnum : return "SmbLapseRateNeg"; 497 case SmbLapseRatePosEnum : return "SmbLapseRatePos"; 497 case SmbLapseRatesEnum : return "SmbLapseRates"; 498 498 case SmbNumBasinsEnum : return "SmbNumBasins"; 499 case SmbNumElevationBinsEnum : return "SmbNumElevationBins"; 499 500 case SmbNumRequestedOutputsEnum : return "SmbNumRequestedOutputs"; 500 501 case SmbPfacEnum : return "SmbPfac"; … … 1088 1089 case TransientAccumulatedDeltaIceThicknessEnum : return "TransientAccumulatedDeltaIceThickness"; 1089 1090 case VelEnum : return "Vel"; 1091 case VVmismipglposinputEnum : return "VVmismipglposinput"; 1090 1092 case VxAverageEnum : return "VxAverage"; 1091 1093 case VxBaseEnum : return "VxBase"; -
issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim
r26859 r26947 458 458 syn keyword juliaConstC SmbDpermilEnum 459 459 syn keyword juliaConstC SmbDsnowIdxEnum 460 syn keyword juliaConstC SmbElevationBinsEnum 460 461 syn keyword juliaConstC SmbCldFracEnum 461 462 syn keyword juliaConstC SmbDelta18oEnum … … 485 486 syn keyword juliaConstC SmbIsturbulentfluxEnum 486 487 syn keyword juliaConstC SmbKEnum 487 syn keyword juliaConstC SmbLapseRateNegEnum 488 syn keyword juliaConstC SmbLapseRatePosEnum 488 syn keyword juliaConstC SmbLapseRatesEnum 489 489 syn keyword juliaConstC SmbNumBasinsEnum 490 syn keyword juliaConstC SmbNumElevationBinsEnum 490 491 syn keyword juliaConstC SmbNumRequestedOutputsEnum 491 492 syn keyword juliaConstC SmbPfacEnum … … 1079 1080 syn keyword juliaConstC TransientAccumulatedDeltaIceThicknessEnum 1080 1081 syn keyword juliaConstC VelEnum 1082 syn keyword juliaConstC VVmismipglposinputEnum 1081 1083 syn keyword juliaConstC VxAverageEnum 1082 1084 syn keyword juliaConstC VxBaseEnum -
issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
r26859 r26947 476 476 else if (strcmp(name,"SmbDpermil")==0) return SmbDpermilEnum; 477 477 else if (strcmp(name,"SmbDsnowIdx")==0) return SmbDsnowIdxEnum; 478 else if (strcmp(name,"SmbElevationBins")==0) return SmbElevationBinsEnum; 478 479 else if (strcmp(name,"SmbCldFrac")==0) return SmbCldFracEnum; 479 480 else if (strcmp(name,"SmbDelta18o")==0) return SmbDelta18oEnum; … … 503 504 else if (strcmp(name,"SmbIsturbulentflux")==0) return SmbIsturbulentfluxEnum; 504 505 else if (strcmp(name,"SmbK")==0) return SmbKEnum; 505 else if (strcmp(name,"SmbLapseRateNeg")==0) return SmbLapseRateNegEnum; 506 else if (strcmp(name,"SmbLapseRatePos")==0) return SmbLapseRatePosEnum; 506 else if (strcmp(name,"SmbLapseRates")==0) return SmbLapseRatesEnum; 507 507 else if (strcmp(name,"SmbNumBasins")==0) return SmbNumBasinsEnum; 508 508 else stage=5; 509 509 } 510 510 if(stage==5){ 511 if (strcmp(name,"SmbNumRequestedOutputs")==0) return SmbNumRequestedOutputsEnum; 511 if (strcmp(name,"SmbNumElevationBins")==0) return SmbNumElevationBinsEnum; 512 else if (strcmp(name,"SmbNumRequestedOutputs")==0) return SmbNumRequestedOutputsEnum; 512 513 else if (strcmp(name,"SmbPfac")==0) return SmbPfacEnum; 513 514 else if (strcmp(name,"SmbPhi")==0) return SmbPhiEnum; … … 628 629 else if (strcmp(name,"BasalforcingsGeothermalflux")==0) return BasalforcingsGeothermalfluxEnum; 629 630 else if (strcmp(name,"BasalforcingsGroundediceMeltingRate")==0) return BasalforcingsGroundediceMeltingRateEnum; 630 else if (strcmp(name,"BasalforcingsLinearBasinId")==0) return BasalforcingsLinearBasinIdEnum;631 631 else stage=6; 632 632 } 633 633 if(stage==6){ 634 if (strcmp(name,"BasalforcingsPerturbationMeltingRate")==0) return BasalforcingsPerturbationMeltingRateEnum; 634 if (strcmp(name,"BasalforcingsLinearBasinId")==0) return BasalforcingsLinearBasinIdEnum; 635 else if (strcmp(name,"BasalforcingsPerturbationMeltingRate")==0) return BasalforcingsPerturbationMeltingRateEnum; 635 636 else if (strcmp(name,"BasalforcingsSpatialDeepwaterElevation")==0) return BasalforcingsSpatialDeepwaterElevationEnum; 636 637 else if (strcmp(name,"BasalforcingsSpatialDeepwaterMeltingRate")==0) return BasalforcingsSpatialDeepwaterMeltingRateEnum; … … 751 752 else if (strcmp(name,"FrictionP")==0) return FrictionPEnum; 752 753 else if (strcmp(name,"FrictionPressureAdjustedTemperature")==0) return FrictionPressureAdjustedTemperatureEnum; 753 else if (strcmp(name,"FrictionQ")==0) return FrictionQEnum;754 754 else stage=7; 755 755 } 756 756 if(stage==7){ 757 if (strcmp(name,"FrictionSedimentCompressibilityCoefficient")==0) return FrictionSedimentCompressibilityCoefficientEnum; 757 if (strcmp(name,"FrictionQ")==0) return FrictionQEnum; 758 else if (strcmp(name,"FrictionSedimentCompressibilityCoefficient")==0) return FrictionSedimentCompressibilityCoefficientEnum; 758 759 else if (strcmp(name,"FrictionTillFrictionAngle")==0) return FrictionTillFrictionAngleEnum; 759 760 else if (strcmp(name,"FrictionWaterLayer")==0) return FrictionWaterLayerEnum; … … 874 875 else if (strcmp(name,"SealevelBarystaticIceWeights")==0) return SealevelBarystaticIceWeightsEnum; 875 876 else if (strcmp(name,"SealevelBarystaticIceArea")==0) return SealevelBarystaticIceAreaEnum; 876 else if (strcmp(name,"SealevelBarystaticIceLatbar")==0) return SealevelBarystaticIceLatbarEnum;877 877 else stage=8; 878 878 } 879 879 if(stage==8){ 880 if (strcmp(name,"SealevelBarystaticIceLongbar")==0) return SealevelBarystaticIceLongbarEnum; 880 if (strcmp(name,"SealevelBarystaticIceLatbar")==0) return SealevelBarystaticIceLatbarEnum; 881 else if (strcmp(name,"SealevelBarystaticIceLongbar")==0) return SealevelBarystaticIceLongbarEnum; 881 882 else if (strcmp(name,"SealevelBarystaticIceLoad")==0) return SealevelBarystaticIceLoadEnum; 882 883 else if (strcmp(name,"SealevelBarystaticHydroMask")==0) return SealevelBarystaticHydroMaskEnum; … … 997 998 else if (strcmp(name,"SmbGdnini")==0) return SmbGdniniEnum; 998 999 else if (strcmp(name,"SmbGsp")==0) return SmbGspEnum; 999 else if (strcmp(name,"SmbGspini")==0) return SmbGspiniEnum;1000 1000 else stage=9; 1001 1001 } 1002 1002 if(stage==9){ 1003 if (strcmp(name,"SmbHref")==0) return SmbHrefEnum; 1003 if (strcmp(name,"SmbGspini")==0) return SmbGspiniEnum; 1004 else if (strcmp(name,"SmbHref")==0) return SmbHrefEnum; 1004 1005 else if (strcmp(name,"SmbIsInitialized")==0) return SmbIsInitializedEnum; 1005 1006 else if (strcmp(name,"SmbMAdd")==0) return SmbMAddEnum; … … 1112 1113 else if (strcmp(name,"TransientAccumulatedDeltaIceThickness")==0) return TransientAccumulatedDeltaIceThicknessEnum; 1113 1114 else if (strcmp(name,"Vel")==0) return VelEnum; 1115 else if (strcmp(name,"VVmismipglposinput")==0) return VVmismipglposinputEnum; 1114 1116 else if (strcmp(name,"VxAverage")==0) return VxAverageEnum; 1115 1117 else if (strcmp(name,"VxBase")==0) return VxBaseEnum; … … 1119 1121 else if (strcmp(name,"VxShear")==0) return VxShearEnum; 1120 1122 else if (strcmp(name,"VxSurface")==0) return VxSurfaceEnum; 1121 else if (strcmp(name,"VyAverage")==0) return VyAverageEnum;1122 else if (strcmp(name,"VyBase")==0) return VyBaseEnum;1123 1123 else stage=10; 1124 1124 } 1125 1125 if(stage==10){ 1126 if (strcmp(name,"Vy")==0) return VyEnum; 1126 if (strcmp(name,"VyAverage")==0) return VyAverageEnum; 1127 else if (strcmp(name,"VyBase")==0) return VyBaseEnum; 1128 else if (strcmp(name,"Vy")==0) return VyEnum; 1127 1129 else if (strcmp(name,"VyMesh")==0) return VyMeshEnum; 1128 1130 else if (strcmp(name,"VyObs")==0) return VyObsEnum; … … 1242 1244 else if (strcmp(name,"Outputdefinition95")==0) return Outputdefinition95Enum; 1243 1245 else if (strcmp(name,"Outputdefinition96")==0) return Outputdefinition96Enum; 1244 else if (strcmp(name,"Outputdefinition97")==0) return Outputdefinition97Enum;1245 else if (strcmp(name,"Outputdefinition98")==0) return Outputdefinition98Enum;1246 1246 else stage=11; 1247 1247 } 1248 1248 if(stage==11){ 1249 if (strcmp(name,"Outputdefinition99")==0) return Outputdefinition99Enum; 1249 if (strcmp(name,"Outputdefinition97")==0) return Outputdefinition97Enum; 1250 else if (strcmp(name,"Outputdefinition98")==0) return Outputdefinition98Enum; 1251 else if (strcmp(name,"Outputdefinition99")==0) return Outputdefinition99Enum; 1250 1252 else if (strcmp(name,"Outputdefinition9")==0) return Outputdefinition9Enum; 1251 1253 else if (strcmp(name,"Outputdefinition100")==0) return Outputdefinition100Enum; … … 1365 1367 else if (strcmp(name,"GaussSeg")==0) return GaussSegEnum; 1366 1368 else if (strcmp(name,"GaussTetra")==0) return GaussTetraEnum; 1367 else if (strcmp(name,"GaussTria")==0) return GaussTriaEnum;1368 else if (strcmp(name,"GenericOption")==0) return GenericOptionEnum;1369 1369 else stage=12; 1370 1370 } 1371 1371 if(stage==12){ 1372 if (strcmp(name,"GenericParam")==0) return GenericParamEnum; 1372 if (strcmp(name,"GaussTria")==0) return GaussTriaEnum; 1373 else if (strcmp(name,"GenericOption")==0) return GenericOptionEnum; 1374 else if (strcmp(name,"GenericParam")==0) return GenericParamEnum; 1373 1375 else if (strcmp(name,"GenericExternalResult")==0) return GenericExternalResultEnum; 1374 1376 else if (strcmp(name,"Gradient1")==0) return Gradient1Enum; … … 1488 1490 else if (strcmp(name,"NoMeltOnPartiallyFloating")==0) return NoMeltOnPartiallyFloatingEnum; 1489 1491 else if (strcmp(name,"Nodal")==0) return NodalEnum; 1490 else if (strcmp(name,"Nodalvalue")==0) return NodalvalueEnum;1491 else if (strcmp(name,"NodeSId")==0) return NodeSIdEnum;1492 1492 else stage=13; 1493 1493 } 1494 1494 if(stage==13){ 1495 if (strcmp(name,"NoneApproximation")==0) return NoneApproximationEnum; 1495 if (strcmp(name,"Nodalvalue")==0) return NodalvalueEnum; 1496 else if (strcmp(name,"NodeSId")==0) return NodeSIdEnum; 1497 else if (strcmp(name,"NoneApproximation")==0) return NoneApproximationEnum; 1496 1498 else if (strcmp(name,"None")==0) return NoneEnum; 1497 1499 else if (strcmp(name,"Numberedcostfunction")==0) return NumberedcostfunctionEnum; … … 1611 1613 else if (strcmp(name,"TotalSmb")==0) return TotalSmbEnum; 1612 1614 else if (strcmp(name,"TotalSmbScaled")==0) return TotalSmbScaledEnum; 1613 else if (strcmp(name,"TransientArrayParam")==0) return TransientArrayParamEnum;1614 else if (strcmp(name,"TransientInput")==0) return TransientInputEnum;1615 1615 else stage=14; 1616 1616 } 1617 1617 if(stage==14){ 1618 if (strcmp(name,"TransientParam")==0) return TransientParamEnum; 1618 if (strcmp(name,"TransientArrayParam")==0) return TransientArrayParamEnum; 1619 else if (strcmp(name,"TransientInput")==0) return TransientInputEnum; 1620 else if (strcmp(name,"TransientParam")==0) return TransientParamEnum; 1619 1621 else if (strcmp(name,"TransientSolution")==0) return TransientSolutionEnum; 1620 1622 else if (strcmp(name,"Tria")==0) return TriaEnum; -
issm/trunk-jpl/src/m/classes/SMBautoregression.m
r26832 r26947 14 14 phi = NaN; 15 15 basin_id = NaN; 16 lapserate _pos= NaN;17 lapserate_neg= NaN;16 lapserates = NaN; 17 elevationbins = NaN; 18 18 refelevation = NaN; 19 19 steps_per_step = 1; … … 62 62 self.ar_order = 0.0; %autoregression model of order 0 63 63 end % }}} 64 function md = checkconsistency(self,md,solution,analyses) 64 function md = checkconsistency(self,md,solution,analyses) % {{{ 65 65 66 66 if ismember('MasstransportAnalysis',analyses), … … 73 73 md = checkfield(md,'fieldname','smb.ar_timestep','numel',1,'NaN',1,'Inf',1,'>=',md.timestepping.time_step); %autoregression time step cannot be finer than ISSM timestep 74 74 md = checkfield(md,'fieldname','smb.phi','NaN',1,'Inf',1,'size',[md.smb.num_basins,md.smb.ar_order]); 75 76 if (any(isnan(md.smb. lapserate_pos)==0) || numel(md.smb.lapserate_pos)>1)77 md = checkfield(md,'fieldname','smb. lapserate_pos','NaN',1,'Inf',1,'size',[1,md.smb.num_basins],'numel',md.smb.num_basins);75 76 if (any(isnan(md.smb.refelevation)==0) || numel(md.smb.refelevation)>1) 77 md = checkfield(md,'fieldname','smb.refelevation','NaN',1,'Inf',1,'>=',0,'size',[1,md.smb.num_basins],'numel',md.smb.num_basins); 78 78 end 79 if (any(isnan(md.smb.lapserate_neg)==0) || numel(md.smb.lapserate_neg)>1) 80 md = checkfield(md,'fieldname','smb.lapserate_neg','NaN',1,'Inf',1,'size',[1,md.smb.num_basins],'numel',md.smb.num_basins); 81 end 82 if (any(isnan(md.smb.refelevation)==0) || numel(md.smb.refelevation)>1) 83 md = checkfield(md,'fieldname','smb.refelevation','NaN',1,'Inf',1,'>=',0,'size',[1,md.smb.num_basins],'numel',md.smb.num_basins); 79 [nbas,nbins] = size(md.smb.lapserates); 80 if (any(isnan(reshape(md.smb.lapserates,[1,nbas*nbins]))==0) || numel(md.smb.lapserates)>1) 81 md = checkfield(md,'fieldname','smb.lapserates','NaN',1,'Inf',1,'size',[md.smb.num_basins,nbins],'numel',md.smb.num_basins*nbins); 82 md = checkfield(md,'fieldname','smb.elevationbins','NaN',1,'Inf',1,'size',[md.smb.num_basins,nbins-1],'numel',md.smb.num_basins*(nbins-1)); 83 if(issorted(md.smb.elevationbins,2)==0) 84 error('md.smb.elevationbins should have rows in order of increasing elevation'); 85 end 86 elseif (isnan(md.smb.elevationbins(1,1))==0 || numel(md.smb.elevationbins)>1) 87 %elevationbins specified but not lapserates: this will inevitably lead to inconsistencies 88 [nbas,nbins] = size(md.smb.elevationbins); 89 nbins = nbins+1; 90 md = checkfield(md,'fieldname','smb.lapserates','NaN',1,'Inf',1,'size',[md.smb.num_basins,nbins],'numel',md.smb.num_basins*nbins); 91 md = checkfield(md,'fieldname','smb.elevationbins','NaN',1,'Inf',1,'size',[md.smb.num_basins,nbins-1],'numel',md.smb.num_basins*(nbins-1)); 84 92 end 85 86 93 end 87 94 md = checkfield(md,'fieldname','smb.steps_per_step','>=',1,'numel',[1]); 88 95 md = checkfield(md,'fieldname','smb.averaging','numel',[1],'values',[0 1 2]); 89 96 md = checkfield(md,'fieldname','smb.requested_outputs','stringrow',1); 90 end 97 end % }}} 91 98 function disp(self) % {{{ 92 99 disp(sprintf(' surface forcings parameters:')); … … 99 106 fielddisplay(self,'ar_timestep','time resolution of the autoregressive model [yr]'); 100 107 fielddisplay(self,'phi','basin-specific vectors of lag coefficients [unitless]'); 101 fielddisplay(self,'lapserate _pos','basin-specific SMB lapse rates applied in range of SMB>=0[m ice eq yr^-1 m^-1] (default: no lapse rate)');102 fielddisplay(self,' lapserate_neg','basin-specific SMB lapse rates applied in range of SMB<0 [m ice eq yr^-1 m^-1] (default: no lapse rate');103 fielddisplay(self,'refelevation','basin-specific reference elevations on which SMB lapse rates are applied(default: basin mean elevation) [m]');108 fielddisplay(self,'lapserates','basin-specific SMB lapse rates applied in each elevation bin, 1 row per basin, 1 column per bin [m ice eq yr^-1 m^-1] (default: no lapse rate)'); 109 fielddisplay(self,'elevationbins','basin-specific separations between elevation bins, 1 row per basin, 1 column per limit between bins [m] (default: no basin separation)'); 110 fielddisplay(self,'refelevation','basin-specific reference elevations at which SMB is calculated, and from which SMB is downscaled using lapserates (default: basin mean elevation) [m]'); 104 111 fielddisplay(self, 'steps_per_step', 'number of smb steps per time step'); 105 112 fielddisplay(self, 'averaging', 'averaging methods from short to long steps'); … … 114 121 yts=md.constants.yts; 115 122 116 templapserate _pos = md.smb.lapserate_pos;117 temp lapserate_neg = md.smb.lapserate_neg;123 templapserates = md.smb.lapserates; 124 tempelevationbins = md.smb.elevationbins; 118 125 temprefelevation = md.smb.refelevation; 119 if(any(isnan(md.smb.lapserate_pos))) 120 templapserate_pos = zeros(1,md.smb.num_basins); 121 disp(' smb.lapserate_pos not specified: set to 0'); 126 [nbas,nbins] = size(md.smb.lapserates); 127 if(any(isnan(reshape(md.smb.lapserates,[1,nbas*nbins])))) 128 templapserates = zeros(md.smb.num_basins,2); 129 disp(' smb.lapserates not specified: set to 0'); 130 tempelevationbins = zeros(md.smb.num_basins,1); %dummy elevation bins 122 131 end 123 if(any(isnan(md.smb.lapserate_neg)))124 templapserate_neg = zeros(1,md.smb.num_basins);125 disp(' smb.lapserate_neg not specified: set to 0');126 end127 132 if(any(isnan(md.smb.refelevation))) 128 133 temprefelevation = zeros(1,md.smb.num_basins); … … 136 141 temprefelevation(ii) = sum(areas(indices).*elemsh)/sum(areas(indices)); 137 142 end 138 if(any( [templapserate_pos,templapserate_neg])~=0)143 if(any(reshape(md.smb.lapserates,[1,nbas*nbins])~=0)) 139 144 disp(' smb.refelevation not specified: Reference elevations set to mean surface elevation of basins'); 140 145 end 141 146 end 147 [nbas,nbins] = size(templapserates); 142 148 143 149 WriteData(fid,prefix,'name','md.smb.model','data',13,'format','Integer'); … … 150 156 WriteData(fid,prefix,'object',self,'class','smb','fieldname','beta1','format','DoubleMat','name','md.smb.beta1','scale',1./(yts^2),'yts',yts); 151 157 WriteData(fid,prefix,'object',self,'class','smb','fieldname','phi','format','DoubleMat','name','md.smb.phi','yts',yts); 152 WriteData(fid,prefix,'data',templapserate _pos,'format','DoubleMat','name','md.smb.lapserate_pos','scale',1./yts,'yts',yts);153 WriteData(fid,prefix,'data',temp lapserate_neg,'format','DoubleMat','name','md.smb.lapserate_neg','scale',1./yts,'yts',yts);158 WriteData(fid,prefix,'data',templapserates,'format','DoubleMat','name','md.smb.lapserates','scale',1./yts,'yts',yts); 159 WriteData(fid,prefix,'data',tempelevationbins,'format','DoubleMat','name','md.smb.elevationbins'); 154 160 WriteData(fid,prefix,'data',temprefelevation,'format','DoubleMat','name','md.smb.refelevation'); 161 WriteData(fid,prefix,'data',nbins,'format','Integer','name','md.smb.num_bins'); 155 162 WriteData(fid,prefix,'object',self,'fieldname','steps_per_step','format','Integer'); 156 163 WriteData(fid,prefix,'object',self,'fieldname','averaging','format','Integer'); -
issm/trunk-jpl/src/m/classes/SMBautoregression.py
r26905 r26947 23 23 self.phi = np.nan 24 24 self.basin_id = np.nan 25 self.lapserate _pos = np.nan26 self. lapserate_neg= np.nan25 self.lapserates = np.nan 26 self.elevationbins = np.nan 27 27 self.refelevation = np.nan 28 28 self.steps_per_step = 1 … … 47 47 s += '{}\n'.format(fielddisplay(self, 'ar_timestep', 'time resolution of the autoregressive model [yr]')) 48 48 s += '{}\n'.format(fielddisplay(self, 'phi', 'basin-specific vectors of lag coefficients [unitless]')) 49 s += '{}\n'.format(fielddisplay(self, 'lapserate _pos', 'basin-specific SMB lapse rates applied in range of SMB>=0[m ice eq yr^-1 m^-1] (default: no lapse rate)'))50 s += '{}\n'.format(fielddisplay(self, ' lapserate_neg', 'basin-specific SMB lapse rates applied in range of SMB<0 [m ice eq yr^-1 m^-1] (default: no lapse rate)'))51 s += '{}\n'.format(fielddisplay(self, 'refelevation', 'basin-specific reference elevations on which SMB lapse rates are applied(default: basin mean elevation) [m]'))49 s += '{}\n'.format(fielddisplay(self, 'lapserates', 'basin-specific SMB lapse rates applied in each elevation bin, 1 row per basin, 1 column per bin [m ice eq yr^-1 m^-1] (default: no lapse rate)')) 50 s += '{}\n'.format(fielddisplay(self, 'elevationbins', 'basin-specific SMB lapse rates applied in range of SMB<0 [m ice eq yr^-1 m^-1] (default: no lapse rate)')) 51 s += '{}\n'.format(fielddisplay(self, 'refelevation', 'basin-specific reference elevations at which SMB is calculated, and from which SMB is downscaled using lapserates (default: basin mean elevation) [m]')) 52 52 s += '{}\n'.format(fielddisplay(self, 'steps_per_step', 'number of smb steps per time step')) 53 53 s += '{}\n'.format(fielddisplay(self, 'averaging', 'averaging methods from short to long steps')) … … 104 104 md = checkfield(md, 'fieldname', 'smb.ar_timestep', 'numel', 1, 'NaN', 1, 'Inf', 1, '>=', md.timestepping.time_step) # Autoregression time step cannot be finer than ISSM timestep 105 105 md = checkfield(md, 'fieldname', 'smb.phi', 'NaN', 1, 'Inf', 1, 'size', [md.smb.num_basins, md.smb.ar_order]) 106 107 if(np.any(np.isnan(self.lapserate_pos) is False) or np.size(self.lapserate_pos) > 1): 108 if len(np.shape(self.lapserate_pos)) == 1: 109 self.lapserate_pos = np.array([self.lapserate_pos]) 110 md = checkfield(md, 'fieldname', 'smb.lapserate_pos', 'NaN', 1, 'Inf', 1, 'size', [1, md.smb.num_basins], 'numel', md.smb.num_basins) 111 112 if(np.any(np.isnan(self.lapserate_neg) is False) or np.size(self.lapserate_neg) > 1): 113 if len(np.shape(self.lapserate_neg)) == 1: 114 self.lapserate_neg = np.array([self.lapserate_neg]) 115 md = checkfield(md, 'fieldname', 'smb.lapserate_neg', 'NaN', 1, 'Inf', 1, 'size', [1, md.smb.num_basins], 'numel', md.smb.num_basins) 116 106 117 107 if(np.any(np.isnan(self.refelevation) is False) or np.size(self.refelevation) > 1): 118 108 if len(np.shape(self.refelevation)) == 1: 119 109 self.refelevation = np.array([self.refelevation]) 120 110 md = checkfield(md, 'fieldname', 'smb.refelevation', 'NaN', 1, 'Inf', 1, '>=', 0, 'size', [1, md.smb.num_basins], 'numel', md.smb.num_basins) 111 112 nbins = np.shape(self.lapserates)[1] 113 if(np.any(np.isnan(self.lapserates) is False) or np.size(self.lapserates) > 1): 114 if len(np.shape(self.lapserates)) == 1: 115 self.lapserates = np.array([self.lapserates]) 116 if len(np.shape(self.elevationbins)) == 1: 117 self.elevationbins = np.array([self.elevationbins]) 118 md = checkfield(md, 'fieldname', 'smb.lapserates', 'NaN', 1, 'Inf', 1, 'size', [md.smb.num_basins, nbins], 'numel', md.smb.num_basins*nbins) 119 md = checkfield(md, 'fieldname', 'smb.elevationbins', 'NaN', 1, 'Inf', 1, 'size', [md.smb.num_basins, nbins-1], 'numel', md.smb.num_basins*(nbins-1)) 120 for rr in range(md.smb.num_basins): 121 if(np.all(self.elevationbins[rr,0:-1]<=self.elevationbins[rr,1:])==False): 122 raise TypeError('md.smb.elevationbins should have rows in order of increasing elevation') 123 elif(np.any(np.isnan(self.elevationbins) is False) or np.size(self.elevationbins) > 1): 124 #elevationbins specified but not lapserates: this will inevitably lead to inconsistencies 125 nbins = np.shape(self.elevationbins)[1] 126 nbins += 1 127 md = checkfield(md, 'fieldname', 'smb.lapserates', 'NaN', 1, 'Inf', 1, 'size', [md.smb.num_basins, nbins], 'numel', md.smb.num_basins*nbins) 128 md = checkfield(md, 'fieldname', 'smb.elevationbins', 'NaN', 1, 'Inf', 1, 'size', [md.smb.num_basins, nbins-1], 'numel', md.smb.num_basins*(nbins-1)) 121 129 122 130 md = checkfield(md, 'fieldname', 'smb.steps_per_step', '>=', 1, 'numel', [1]) … … 129 137 yts = md.constants.yts 130 138 131 templapserate _pos = np.copy(md.smb.lapserate_pos)132 temp lapserate_neg = np.copy(md.smb.lapserate_neg)139 templapserates = np.copy(md.smb.lapserates) 140 tempelevationbins = np.copy(md.smb.elevationbins) 133 141 temprefelevation = np.copy(md.smb.refelevation) 134 if(np.any(np.isnan(md.smb.lapserate_pos))): 135 templapserate_pos = np.zeros((md.smb.num_basins)).reshape(1,md.smb.num_basins) 136 print(' smb.lapserate_pos not specified: set to 0') 137 if(np.any(np.isnan(md.smb.lapserate_neg))): 138 templapserate_neg = np.zeros((md.smb.num_basins)).reshape(1,md.smb.num_basins) 139 print(' smb.lapserate_neg not specified: set to 0') 142 if(np.any(np.isnan(md.smb.lapserates))): 143 templapserates = np.zeros((md.smb.num_basins,2)) 144 print(' smb.lapserates not specified: set to 0') 145 tempelevationbins = np.zeros((md.smb.num_basins,1)) #dummy elevation bins 140 146 if(np.any(np.isnan(md.smb.refelevation))): 141 147 temprefelevation = np.zeros((md.smb.num_basins)).reshape(1,md.smb.num_basins) … … 147 153 elemsh[jj] = np.mean(md.geometry.surface[md.mesh.elements[indices[jj], :] - 1]) 148 154 temprefelevation[0, ii] = np.sum(areas[indices] * elemsh) / np.sum(areas[indices]) 149 if(np.any(templapserate _pos != 0) or np.any(templapserate_neg!= 0)):155 if(np.any(templapserates != 0)): 150 156 print(' smb.refelevation not specified: Reference elevations set to mean surface elevation of basins') 157 nbins = np.shape(templapserates)[1] 151 158 152 159 WriteData(fid, prefix, 'name', 'md.smb.model', 'data', 13, 'format', 'Integer') … … 159 166 WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'beta1', 'format', 'DoubleMat', 'name', 'md.smb.beta1', 'scale', 1 / (yts ** 2), 'yts', yts) 160 167 WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'phi', 'format', 'DoubleMat', 'name', 'md.smb.phi', 'yts', yts) 161 WriteData(fid, prefix, 'data', templapserate _pos, 'name', 'md.smb.lapserate_pos', 'format', 'DoubleMat', 'scale', 1 / yts, 'yts', yts)162 WriteData(fid, prefix, 'data', temp lapserate_neg, 'name', 'md.smb.lapserate_neg', 'format', 'DoubleMat', 'scale', 1 / yts, 'yts', yts)168 WriteData(fid, prefix, 'data', templapserates, 'name', 'md.smb.lapserates', 'format', 'DoubleMat', 'scale', 1 / yts, 'yts', yts) 169 WriteData(fid, prefix, 'data', tempelevationbins, 'name', 'md.smb.elevationbins', 'format', 'DoubleMat') 163 170 WriteData(fid, prefix, 'data', temprefelevation, 'name', 'md.smb.refelevation', 'format', 'DoubleMat') 171 WriteData(fid, prefix, 'data', nbins, 'name', 'md.smb.num_bins', 'format', 'Integer') 164 172 WriteData(fid, prefix, 'object', self, 'fieldname', 'steps_per_step', 'format', 'Integer') 165 173 WriteData(fid, prefix, 'object', self, 'fieldname', 'averaging', 'format', 'Integer') -
issm/trunk-jpl/test/NightlyRun/test257.m
r26810 r26947 45 45 md.smb.ar_timestep = 2.0; %timestep of the autoregressive model [yr] 46 46 md.smb.phi = [[0.2,0.1,0.05,0.01];[0.4,0.2,-0.2,0.1];[0.4,-0.4,0.1,-0.1]]; 47 md.smb.lapserate _pos = [0.01,0.0,0.0];48 md.smb. lapserate_neg = [0.01,0.0,0.0];47 md.smb.lapserates = [0.01,0.0;0.01,-0.01;0.0,-0.01]; 48 md.smb.elevationbins = [100;150;100]; 49 49 50 50 %Stochastic forcing … … 53 53 md.stochasticforcing.covariance = [[0.15 0.08 -0.02];[0.08 0.12 -0.05];[-0.02 -0.05 0.1]]; %global covariance among- and between-fields 54 54 md.stochasticforcing.randomflag = 0; %fixed random seeds 55 md.stochasticforcing.stochastictimestep = 1.0; 55 56 56 57 md=solve(md,'Transient'); -
issm/trunk-jpl/test/NightlyRun/test257.py
r26907 r26947 50 50 md.smb.ar_timestep = 2.0 #timestep of the autoregressive model [yr] 51 51 md.smb.phi = np.array([[0.2, 0.1, 0.05, 0.01], [0.4, 0.2, -0.2, 0.1], [0.4, -0.4, 0.1, -0.1]]) 52 md.smb.lapserate _pos = np.array([[0.01,0,0]])53 md.smb. lapserate_neg = np.array([[0.01,0,0]])52 md.smb.lapserates = np.array([[0.01,0.0],[0.01,-0.01],[0.0,-0.01]]) 53 md.smb.elevationbins = np.array([100,150,100]).reshape(md.smb.num_basins,1) 54 54 55 55 # Stochastic forcing … … 58 58 md.stochasticforcing.covariance = np.array([[0.15, 0.08, -0.02], [0.08, 0.12, -0.05], [-0.02, -0.05, 0.1]]) # global covariance among- and between-fields 59 59 md.stochasticforcing.randomflag = 0 # fixed random seeds 60 60 md.stochasticforcing.stochastictimestep = 1.0 61 61 62 62 md = solve(md, 'Transient')
Note:
See TracChangeset
for help on using the changeset viewer.