Changeset 26789


Ignore:
Timestamp:
01/17/22 10:32:37 (3 years ago)
Author:
vverjans
Message:

BUG: fixed spatiallinearbasalforcings() class and added parameters for UpperWaterMeltRate and PerturbationMeltRate

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

Legend:

Unmodified
Added
Removed

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

    r26640 r26789  
    108108                        break;
    109109                case SpatialLinearFloatingMeltRateEnum:
    110                         iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.deepwater_melting_rate",BasalforcingsDeepwaterMeltingRateEnum);
    111                         iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.deepwater_elevation",BasalforcingsDeepwaterElevationEnum);
    112                         iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.upperwater_elevation",BasalforcingsUpperwaterElevationEnum);
     110                        iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.deepwater_melting_rate",BasalforcingsSpatialDeepwaterMeltingRateEnum);
     111                        iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.deepwater_elevation",BasalforcingsSpatialDeepwaterElevationEnum);
     112                        iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.upperwater_melting_rate",BasalforcingsSpatialUpperwaterMeltingRateEnum);
     113                        iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.upperwater_elevation",BasalforcingsSpatialUpperwaterElevationEnum);
    113114                        break;
    114115                case BasalforcingsPicoEnum:

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

    r26640 r26789  
    182182                        break;
    183183                case SpatialLinearFloatingMeltRateEnum:
    184                         iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.deepwater_melting_rate",BasalforcingsDeepwaterMeltingRateEnum);
    185                         iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.deepwater_elevation",BasalforcingsDeepwaterElevationEnum);
    186                         iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.upperwater_elevation",BasalforcingsUpperwaterElevationEnum);
     184                        iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.deepwater_melting_rate",BasalforcingsSpatialDeepwaterMeltingRateEnum);
     185                        iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.deepwater_elevation",BasalforcingsSpatialDeepwaterElevationEnum);
     186                        iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.upperwater_melting_rate",BasalforcingsSpatialUpperwaterMeltingRateEnum);
     187                        iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.upperwater_elevation",BasalforcingsSpatialUpperwaterElevationEnum);
     188                        iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.perturbation_melting_rate",BasalforcingsPerturbationMeltingRateEnum,0.);
    187189                        break;
    188190                case BasalforcingsPicoEnum:

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

    r26640 r26789  
    836836                                break;
    837837                        case SpatialLinearFloatingMeltRateEnum:
    838                                 iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.deepwater_melting_rate",BasalforcingsDeepwaterMeltingRateEnum);
    839                                 iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.deepwater_elevation",BasalforcingsDeepwaterElevationEnum);
    840                                 iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.upperwater_elevation",BasalforcingsUpperwaterElevationEnum);
     838                                iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.deepwater_melting_rate",BasalforcingsSpatialDeepwaterMeltingRateEnum);
     839                                iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.deepwater_elevation",BasalforcingsSpatialDeepwaterElevationEnum);
     840                                iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.upperwater_melting_rate",BasalforcingsSpatialUpperwaterMeltingRateEnum);
     841                                iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.upperwater_elevation",BasalforcingsSpatialUpperwaterElevationEnum);
    841842                                break;
    842843                        case BasalforcingsPicoEnum:

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

    r26640 r26789  
    139139                        break;
    140140                case SpatialLinearFloatingMeltRateEnum:
    141                         iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.deepwater_melting_rate",BasalforcingsDeepwaterMeltingRateEnum);
    142                         iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.deepwater_elevation",BasalforcingsDeepwaterElevationEnum);
    143                         iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.upperwater_elevation",BasalforcingsUpperwaterElevationEnum);
     141                        iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.deepwater_melting_rate",BasalforcingsSpatialDeepwaterMeltingRateEnum);
     142                        iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.deepwater_elevation",BasalforcingsSpatialDeepwaterElevationEnum);
     143                        iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.upperwater_melting_rate",BasalforcingsSpatialUpperwaterMeltingRateEnum);
     144                        iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.upperwater_elevation",BasalforcingsSpatialUpperwaterElevationEnum);
    144145                        break;
    145146                case BasalforcingsPicoEnum:

  • issm/trunk-jpl/src/c/classes/Elements/Element.cpp

    r26615 r26789  
    22982298        const int NUM_VERTICES = this->GetNumberOfVertices();
    22992299
     2300        IssmDouble alpha;
    23002301        IssmDouble deepwatermelt[MAXVERTICES];
    2301         IssmDouble deepwaterel[MAXVERTICES];;
     2302        IssmDouble deepwaterel[MAXVERTICES];
     2303        IssmDouble upperwatermelt[MAXVERTICES];
    23022304        IssmDouble upperwaterel[MAXVERTICES];
     2305        IssmDouble perturbation[MAXVERTICES];
    23032306        IssmDouble base[MAXVERTICES];
    23042307        IssmDouble values[MAXVERTICES];
    23052308
    23062309        this->GetInputListOnVertices(&base[0],BaseEnum);
    2307         this->GetInputListOnVertices(&deepwatermelt[0],BasalforcingsDeepwaterMeltingRateEnum);
    2308         this->GetInputListOnVertices(&deepwaterel[0],BasalforcingsDeepwaterElevationEnum);
    2309         this->GetInputListOnVertices(&upperwaterel[0],BasalforcingsUpperwaterElevationEnum);
     2310        this->GetInputListOnVertices(&deepwatermelt[0],BasalforcingsSpatialDeepwaterMeltingRateEnum);
     2311        this->GetInputListOnVertices(&deepwaterel[0],BasalforcingsSpatialDeepwaterElevationEnum);
     2312        this->GetInputListOnVertices(&upperwatermelt[0],BasalforcingsSpatialUpperwaterMeltingRateEnum);
     2313        this->GetInputListOnVertices(&upperwaterel[0],BasalforcingsSpatialUpperwaterElevationEnum);
     2314   this->GetInputListOnVertices(&perturbation[0],BasalforcingsPerturbationMeltingRateEnum);
    23102315
    23112316        for(int i=0;i<NUM_VERTICES;i++){
    2312                 if(base[i]>upperwaterel[i])      values[i]=0;
    2313                 else if (base[i]<deepwaterel[i]) values[i]=deepwatermelt[i];
    2314                 else values[i]=deepwatermelt[i]*(base[i]-upperwaterel[i])/(deepwaterel[i]-upperwaterel[i]);
     2317                if(base[i]>=upperwaterel[i]){
     2318                        values[i]=upperwatermelt[i];
     2319                }
     2320                else if (base[i]<deepwaterel[i]){
     2321                        values[i]=deepwatermelt[i];
     2322                }
     2323                else{
     2324                        alpha = (base[i]-upperwaterel[i])/(deepwaterel[i]-upperwaterel[i]);
     2325                        values[i]=deepwatermelt[i]*alpha+(1.-alpha)*upperwatermelt[i];
     2326                }
     2327                values[i] += perturbation[i];
    23152328        }
    23162329

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

    r26784 r26789  
    582582syn keyword cConstant BasalforcingsGroundediceMeltingRateEnum
    583583syn keyword cConstant BasalforcingsPerturbationMeltingRateEnum
     584syn keyword cConstant BasalforcingsSpatialDeepwaterElevationEnum
     585syn keyword cConstant BasalforcingsSpatialDeepwaterMeltingRateEnum
     586syn keyword cConstant BasalforcingsSpatialUpperwaterElevationEnum
     587syn keyword cConstant BasalforcingsSpatialUpperwaterMeltingRateEnum
    584588syn keyword cConstant BasalforcingsIsmip6BasinIdEnum
    585589syn keyword cConstant BasalforcingsIsmip6TfEnum

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

    r26784 r26789  
    578578        BasalforcingsGroundediceMeltingRateEnum,
    579579        BasalforcingsPerturbationMeltingRateEnum,
     580        BasalforcingsSpatialDeepwaterElevationEnum,
     581        BasalforcingsSpatialDeepwaterMeltingRateEnum,
     582        BasalforcingsSpatialUpperwaterElevationEnum,
     583        BasalforcingsSpatialUpperwaterMeltingRateEnum,
    580584        BasalforcingsIsmip6BasinIdEnum,
    581585        BasalforcingsIsmip6TfEnum,

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

    r26784 r26789  
    575575syn keyword juliaConstC BasalforcingsGroundediceMeltingRateEnum
    576576syn keyword juliaConstC BasalforcingsPerturbationMeltingRateEnum
     577syn keyword juliaConstC BasalforcingsSpatialDeepwaterElevationEnum
     578syn keyword juliaConstC BasalforcingsSpatialDeepwaterMeltingRateEnum
     579syn keyword juliaConstC BasalforcingsSpatialUpperwaterElevationEnum
     580syn keyword juliaConstC BasalforcingsSpatialUpperwaterMeltingRateEnum
    577581syn keyword juliaConstC BasalforcingsIsmip6BasinIdEnum
    578582syn keyword juliaConstC BasalforcingsIsmip6TfEnum

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

    r26750 r26789  
    595595              else if (strcmp(name,"BasalforcingsGroundediceMeltingRate")==0) return BasalforcingsGroundediceMeltingRateEnum;
    596596              else if (strcmp(name,"BasalforcingsPerturbationMeltingRate")==0) return BasalforcingsPerturbationMeltingRateEnum;
     597              else if (strcmp(name,"BasalforcingsSpatialDeepwaterElevation")==0) return BasalforcingsSpatialDeepwaterElevationEnum;
     598              else if (strcmp(name,"BasalforcingsSpatialDeepwaterMeltingRate")==0) return BasalforcingsSpatialDeepwaterMeltingRateEnum;
     599              else if (strcmp(name,"BasalforcingsSpatialUpperwaterElevation")==0) return BasalforcingsSpatialUpperwaterElevationEnum;
     600              else if (strcmp(name,"BasalforcingsSpatialUpperwaterMeltingRate")==0) return BasalforcingsSpatialUpperwaterMeltingRateEnum;
    597601              else if (strcmp(name,"BasalforcingsIsmip6BasinId")==0) return BasalforcingsIsmip6BasinIdEnum;
    598602              else if (strcmp(name,"BasalforcingsIsmip6Tf")==0) return BasalforcingsIsmip6TfEnum;
    … …  
    625629              else if (strcmp(name,"BedNorth")==0) return BedNorthEnum;
    626630              else if (strcmp(name,"BedNorthGRD")==0) return BedNorthGRDEnum;
    627               else if (strcmp(name,"BedSlopeX")==0) return BedSlopeXEnum;
     631         else stage=6;
     632   }
     633   if(stage==6){
     634              if (strcmp(name,"BedSlopeX")==0) return BedSlopeXEnum;
    628635              else if (strcmp(name,"BedSlopeY")==0) return BedSlopeYEnum;
    629636              else if (strcmp(name,"BottomPressure")==0) return BottomPressureEnum;
    630637              else if (strcmp(name,"BottomPressureOld")==0) return BottomPressureOldEnum;
    631          else stage=6;
    632    }
    633    if(stage==6){
    634               if (strcmp(name,"CalvingCalvingrate")==0) return CalvingCalvingrateEnum;
     638              else if (strcmp(name,"CalvingCalvingrate")==0) return CalvingCalvingrateEnum;
    635639              else if (strcmp(name,"CalvingHabFraction")==0) return CalvingHabFractionEnum;
    636640              else if (strcmp(name,"CalvingMeltingrate")==0) return CalvingMeltingrateEnum;
    … …  
    748752              else if (strcmp(name,"HydrologyGapHeightY")==0) return HydrologyGapHeightYEnum;
    749753              else if (strcmp(name,"HydrologyGapHeightYY")==0) return HydrologyGapHeightYYEnum;
    750               else if (strcmp(name,"HydrologyHead")==0) return HydrologyHeadEnum;
     754         else stage=7;
     755   }
     756   if(stage==7){
     757              if (strcmp(name,"HydrologyHead")==0) return HydrologyHeadEnum;
    751758              else if (strcmp(name,"HydrologyHeadOld")==0) return HydrologyHeadOldEnum;
    752759              else if (strcmp(name,"HydrologyMoulinInput")==0) return HydrologyMoulinInputEnum;
    753760              else if (strcmp(name,"HydrologyNeumannflux")==0) return HydrologyNeumannfluxEnum;
    754          else stage=7;
    755    }
    756    if(stage==7){
    757               if (strcmp(name,"HydrologyReynolds")==0) return HydrologyReynoldsEnum;
     761              else if (strcmp(name,"HydrologyReynolds")==0) return HydrologyReynoldsEnum;
    758762              else if (strcmp(name,"HydrologySheetConductivity")==0) return HydrologySheetConductivityEnum;
    759763              else if (strcmp(name,"HydrologySheetThickness")==0) return HydrologySheetThicknessEnum;
    … …  
    871875              else if (strcmp(name,"SealevelGNrotm2")==0) return SealevelGNrotm2Enum;
    872876              else if (strcmp(name,"SealevelGNrotm3")==0) return SealevelGNrotm3Enum;
    873               else if (strcmp(name,"SealevelGErotm1")==0) return SealevelGErotm1Enum;
     877         else stage=8;
     878   }
     879   if(stage==8){
     880              if (strcmp(name,"SealevelGErotm1")==0) return SealevelGErotm1Enum;
    874881              else if (strcmp(name,"SealevelGErotm2")==0) return SealevelGErotm2Enum;
    875882              else if (strcmp(name,"SealevelGErotm3")==0) return SealevelGErotm3Enum;
    876883              else if (strcmp(name,"SealevelRSLBarystatic")==0) return SealevelRSLBarystaticEnum;
    877          else stage=8;
    878    }
    879    if(stage==8){
    880               if (strcmp(name,"SealevelRSLRate")==0) return SealevelRSLRateEnum;
     884              else if (strcmp(name,"SealevelRSLRate")==0) return SealevelRSLRateEnum;
    881885              else if (strcmp(name,"SealevelUGrd")==0) return SealevelUGrdEnum;
    882886              else if (strcmp(name,"SealevelNGrd")==0) return SealevelNGrdEnum;
    … …  
    994998              else if (strcmp(name,"SmbReini")==0) return SmbReiniEnum;
    995999              else if (strcmp(name,"SmbRunoff")==0) return SmbRunoffEnum;
    996               else if (strcmp(name,"SmbRunoffSubstep")==0) return SmbRunoffSubstepEnum;
     1000         else stage=9;
     1001   }
     1002   if(stage==9){
     1003              if (strcmp(name,"SmbRunoffSubstep")==0) return SmbRunoffSubstepEnum;
    9971004              else if (strcmp(name,"SmbRunoffTransient")==0) return SmbRunoffTransientEnum;
    9981005              else if (strcmp(name,"SmbS0gcm")==0) return SmbS0gcmEnum;
    9991006              else if (strcmp(name,"SmbS0p")==0) return SmbS0pEnum;
    1000          else stage=9;
    1001    }
    1002    if(stage==9){
    1003               if (strcmp(name,"SmbS0t")==0) return SmbS0tEnum;
     1007              else if (strcmp(name,"SmbS0t")==0) return SmbS0tEnum;
    10041008              else if (strcmp(name,"SmbSizeini")==0) return SmbSizeiniEnum;
    10051009              else if (strcmp(name,"SmbSmbCorr")==0) return SmbSmbCorrEnum;
    … …  
    11171121              else if (strcmp(name,"Outputdefinition12")==0) return Outputdefinition12Enum;
    11181122              else if (strcmp(name,"Outputdefinition13")==0) return Outputdefinition13Enum;
    1119               else if (strcmp(name,"Outputdefinition14")==0) return Outputdefinition14Enum;
     1123         else stage=10;
     1124   }
     1125   if(stage==10){
     1126              if (strcmp(name,"Outputdefinition14")==0) return Outputdefinition14Enum;
    11201127              else if (strcmp(name,"Outputdefinition15")==0) return Outputdefinition15Enum;
    11211128              else if (strcmp(name,"Outputdefinition16")==0) return Outputdefinition16Enum;
    11221129              else if (strcmp(name,"Outputdefinition17")==0) return Outputdefinition17Enum;
    1123          else stage=10;
    1124    }
    1125    if(stage==10){
    1126               if (strcmp(name,"Outputdefinition18")==0) return Outputdefinition18Enum;
     1130              else if (strcmp(name,"Outputdefinition18")==0) return Outputdefinition18Enum;
    11271131              else if (strcmp(name,"Outputdefinition19")==0) return Outputdefinition19Enum;
    11281132              else if (strcmp(name,"Outputdefinition20")==0) return Outputdefinition20Enum;
    … …  
    12401244              else if (strcmp(name,"BeckmannGoosseFloatingMeltRate")==0) return BeckmannGoosseFloatingMeltRateEnum;
    12411245              else if (strcmp(name,"BedSlopeSolution")==0) return BedSlopeSolutionEnum;
    1242               else if (strcmp(name,"BoolExternalResult")==0) return BoolExternalResultEnum;
     1246         else stage=11;
     1247   }
     1248   if(stage==11){
     1249              if (strcmp(name,"BoolExternalResult")==0) return BoolExternalResultEnum;
    12431250              else if (strcmp(name,"BoolInput")==0) return BoolInputEnum;
    12441251              else if (strcmp(name,"IntInput")==0) return IntInputEnum;
    12451252              else if (strcmp(name,"DoubleInput")==0) return DoubleInputEnum;
    1246          else stage=11;
    1247    }
    1248    if(stage==11){
    1249               if (strcmp(name,"BoolParam")==0) return BoolParamEnum;
     1253              else if (strcmp(name,"BoolParam")==0) return BoolParamEnum;
    12501254              else if (strcmp(name,"Boundary")==0) return BoundaryEnum;
    12511255              else if (strcmp(name,"BuddJacka")==0) return BuddJackaEnum;
    … …  
    13631367              else if (strcmp(name,"IceVolumeAboveFloatation")==0) return IceVolumeAboveFloatationEnum;
    13641368              else if (strcmp(name,"IceVolumeAboveFloatationScaled")==0) return IceVolumeAboveFloatationScaledEnum;
    1365               else if (strcmp(name,"IceVolume")==0) return IceVolumeEnum;
     1369         else stage=12;
     1370   }
     1371   if(stage==12){
     1372              if (strcmp(name,"IceVolume")==0) return IceVolumeEnum;
    13661373              else if (strcmp(name,"IceVolumeScaled")==0) return IceVolumeScaledEnum;
    13671374              else if (strcmp(name,"IcefrontMassFlux")==0) return IcefrontMassFluxEnum;
    13681375              else if (strcmp(name,"IcefrontMassFluxLevelset")==0) return IcefrontMassFluxLevelsetEnum;
    1369          else stage=12;
    1370    }
    1371    if(stage==12){
    1372               if (strcmp(name,"Incremental")==0) return IncrementalEnum;
     1376              else if (strcmp(name,"Incremental")==0) return IncrementalEnum;
    13731377              else if (strcmp(name,"Indexed")==0) return IndexedEnum;
    13741378              else if (strcmp(name,"IntExternalResult")==0) return IntExternalResultEnum;
    … …  
    14861490              else if (strcmp(name,"Profiler")==0) return ProfilerEnum;
    14871491              else if (strcmp(name,"ProfilingCurrentFlops")==0) return ProfilingCurrentFlopsEnum;
    1488               else if (strcmp(name,"ProfilingCurrentMem")==0) return ProfilingCurrentMemEnum;
     1492         else stage=13;
     1493   }
     1494   if(stage==13){
     1495              if (strcmp(name,"ProfilingCurrentMem")==0) return ProfilingCurrentMemEnum;
    14891496              else if (strcmp(name,"ProfilingSolutionTime")==0) return ProfilingSolutionTimeEnum;
    14901497              else if (strcmp(name,"Regionaloutput")==0) return RegionaloutputEnum;
    14911498              else if (strcmp(name,"Regular")==0) return RegularEnum;
    1492          else stage=13;
    1493    }
    1494    if(stage==13){
    1495               if (strcmp(name,"RecoveryAnalysis")==0) return RecoveryAnalysisEnum;
     1499              else if (strcmp(name,"RecoveryAnalysis")==0) return RecoveryAnalysisEnum;
    14961500              else if (strcmp(name,"Riftfront")==0) return RiftfrontEnum;
    14971501              else if (strcmp(name,"SamplingAnalysis")==0) return SamplingAnalysisEnum;

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

    r25688 r26789  
    99                deepwater_melting_rate    = NaN;
    1010                deepwater_elevation       = NaN;
     11                upperwater_melting_rate   = NaN;
    1112                upperwater_elevation      = NaN;
     13                perturbation_melting_rate = NaN;
    1214                geothermalflux            = NaN;
    1315        end
    … …  
    2527                                                self.deepwater_elevation=lb.deepwater_elevation*ones(nvertices,1);
    2628                                                self.deepwater_melting_rate=lb.deepwater_melting_rate*ones(nvertices,1);
     29                                                self.upperwater_melting_rate=lb.upperwater_melting_rate*ones(nvertices,1);
    2730                                                self.upperwater_elevation=lb.upperwater_elevation*ones(nvertices,1);
     31                                                self.perturbation_melting_rate=lb.perturbation_melting_rate*ones(nvertices,1);
    2832                                        else
    2933                                                self=structtoobj(spatiallinearbasalforcings(),varargin{1});
    … …  
    3741                        self.deepwater_melting_rate=project3d(md,'vector',self.deepwater_melting_rate,'type','node','layer',1);
    3842                        self.deepwater_elevation=project3d(md,'vector',self.deepwater_elevation,'type','node','layer',1);
     43                        self.upperwater_melting_rate=project3d(md,'vector',self.upperwater_melting_rate,'type','node','layer',1);
    3944                        self.upperwater_elevation=project3d(md,'vector',self.upperwater_elevation,'type','node','layer',1);
    4045                        self.geothermalflux=project3d(md,'vector',self.geothermalflux,'type','node','layer',1); %bedrock only gets geothermal flux
     46                        self.perturbation_melting_rate=project3d(md,'vector',self.perturbation_melting_rate,'type','node','layer',1);
    4147                end % }}}
    4248                function self = initialize(self,md) % {{{
    … …  
    5864                function md = checkconsistency(self,md,solution,analyses) % {{{
    5965
     66                        if numel(md.basalforcings.perturbation_melting_rate)>1
     67            md = checkfield(md,'fieldname','basalforcings.perturbation_melting_rate','NaN',1,'Inf',1,'timeseries',1);
     68         end
     69
    6070                        if ismember('MasstransportAnalysis',analyses) & ~(strcmp(solution,'TransientSolution') & md.transient.ismasstransport==0),
    6171                                md = checkfield(md,'fieldname','basalforcings.groundedice_melting_rate','NaN',1,'Inf',1,'timeseries',1);
    6272                                md = checkfield(md,'fieldname','basalforcings.deepwater_melting_rate','NaN',1,'Inf',1,'timeseries',1,'>=',0);
    6373                                md = checkfield(md,'fieldname','basalforcings.deepwater_elevation','NaN',1,'Inf',1,'timeseries',1);
     74                                md = checkfield(md,'fieldname','basalforcings.upperwater_melting_rate','NaN',1,'Inf',1,'timeseries',1,'>=',0);
    6475                                md = checkfield(md,'fieldname','basalforcings.upperwater_elevation','NaN',1,'Inf',1,'timeseries',1,'<',0);
    6576                        end
    … …  
    6980                                md = checkfield(md,'fieldname','basalforcings.deepwater_melting_rate','>=',0,'numel',1);
    7081                                md = checkfield(md,'fieldname','basalforcings.deepwater_elevation','<','basalforcings.upperwater_elevation','numel',1);
     82                                md = checkfield(md,'fieldname','basalforcings.upperwater_melting_rate','>=',0,'numel',1);
    7183                                md = checkfield(md,'fieldname','basalforcings.upperwater_elevation','<=',0,'numel',1);
    7284                        end
    … …  
    7688                                md = checkfield(md,'fieldname','basalforcings.deepwater_melting_rate','>=',0,'numel',1);
    7789                                md = checkfield(md,'fieldname','basalforcings.deepwater_elevation','<','basalforcings.upperwater_elevation','numel',1);
     90                                md = checkfield(md,'fieldname','basalforcings.upperwater_melting_rate','>=',0,'numel',1);
    7891                                md = checkfield(md,'fieldname','basalforcings.upperwater_elevation','<',0,'numel',1);
    7992                                md = checkfield(md,'fieldname','basalforcings.geothermalflux','NaN',1,'Inf',1,'timeseries',1,'>=',0);
    … …  
    8699                        fielddisplay(self,'deepwater_melting_rate','basal melting rate (positive if melting applied for floating ice whith base < deepwater_elevation) [m/yr]');
    87100                        fielddisplay(self,'deepwater_elevation','elevation of ocean deepwater [m]');
     101                        fielddisplay(self,'upperwater_melting_rate','basal melting rate (positive if melting applied for floating ice whith base >= upperwater_elevation) [m/yr]');
    88102                        fielddisplay(self,'upperwater_elevation','elevation of ocean upperwater [m]');
     103                        fielddisplay(self,'perturbation_melting_rate','basal melting rate perturbation added to computed melting rate (positive if melting) [m/yr]');
    89104                        fielddisplay(self,'geothermalflux','geothermal heat flux [W/m^2]');
    90105
    … …  
    94109                        yts=md.constants.yts;
    95110
    96                         floatingice_melting_rate=zeros(md.mesh.numberofvertices,1);
    97                         pos=find(md.geometry.base<=md.basalforcings.deepwater_elevation);
    98                         floatingice_melting_rate(pos)=md.basalforcings.deepwater_melting_rate(pos);
    99                         pos=find(md.geometry.base>md.basalforcings.deepwater_elevation & md.geometry.base<md.basalforcings.upperwater_elevation);
    100                         floatingice_melting_rate(pos)=md.basalforcings.deepwater_melting_rate(pos).*(md.geometry.base(pos)-md.basalforcings.upperwater_elevation(pos))./(md.basalforcings.deepwater_elevation(pos)-md.basalforcings.upperwater_elevation(pos));
    101111                        WriteData(fid,prefix,'name','md.basalforcings.model','data',6,'format','Integer');
    102                         WriteData(fid,prefix,'data',floatingice_melting_rate,'format','DoubleMat','name','md.basalforcings.floatingice_melting_rate','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
    103112                        WriteData(fid,prefix,'object',self,'fieldname','groundedice_melting_rate','format','DoubleMat','name','md.basalforcings.groundedice_melting_rate','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
    104113                        WriteData(fid,prefix,'object',self,'fieldname','geothermalflux','name','md.basalforcings.geothermalflux','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
    105114                        WriteData(fid,prefix,'object',self,'fieldname','deepwater_melting_rate','format','DoubleMat','name','md.basalforcings.deepwater_melting_rate','scale',1./yts,'mattype',1);
    106115                        WriteData(fid,prefix,'object',self,'fieldname','deepwater_elevation','format','DoubleMat','name','md.basalforcings.deepwater_elevation','mattype',1);
     116                        WriteData(fid,prefix,'object',self,'fieldname','upperwater_melting_rate','format','DoubleMat','name','md.basalforcings.upperwater_melting_rate','scale',1./yts,'mattype',1);
    107117                        WriteData(fid,prefix,'object',self,'fieldname','upperwater_elevation','format','DoubleMat','name','md.basalforcings.upperwater_elevation','mattype',1);
     118                        WriteData(fid,prefix,'object',self,'fieldname','perturbation_melting_rate','format','DoubleMat','name','md.basalforcings.perturbation_melting_rate','scale',1./yts,'mattype',1);
    108119                end % }}}
    109120        end

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

    r25762 r26789  
    2121            self.deepwater_melting_rate     = np.nan
    2222            self.deepwater_elevation        = np.nan
     23            self.upperwater_melting_rate    = np.nan
    2324            self.upperwater_elevation       = np.nan
    2425            self.geothermalflux             = np.nan
     26            self.perturbation_melting_rate  = np.nan
    2527
    2628            self.setdefaultparameters()
    … …  
    3032                nvertices = len(lb.groundedice_melting_rate)
    3133                self.groundedice_melting_rate   = lb.groundedice_melting_rate
    32                 self.deepwater_melting_rate     = lb.geothermalflux
     34                self.geothermalflux             = lb.geothermalflux
    3335                self.deepwater_elevation        = lb.deepwater_elevation * np.ones((nvertices, ))
    34                 self.upperwater_elevation       = lb.deepwater_melting_rate * np.ones((nvertices, ))
    35                 self.geothermalflux             = lb.upperwater_elevation * np.ones((nvertices, ))
     36                self.deepwater_melting_rate     = lb.deepwater_melting_rate * np.ones((nvertices, ))
     37                self.upperwater_melting_rate    = lb.upperwater_melting_rate * np.ones((nvertices, ))
     38                self.upperwater_elevation       = lb.upperwater_elevation * np.ones((nvertices, ))
     39                self.perturbation_melting_rate  = lb.perturbation_melting_rate * np.ones((nvertices, ))
    3640            else:
    3741                # TODO: This has not been tested
    … …  
    4650        s += '{}\n'.format(fielddisplay(self, 'deepwater_melting_rate', 'basal melting rate (positive if melting applied for floating ice whith base < deepwater_elevation) [m/yr]'))
    4751        s += '{}\n'.format(fielddisplay(self, 'deepwater_elevation', 'elevation of ocean deepwater [m]'))
     52        s += '{}\n'.format(fielddisplay(self, 'upperwater_melting_rate', 'basal melting rate (positive if melting applied for floating ice whith base >= upperwater_elevation) [m/yr]'))
    4853        s += '{}\n'.format(fielddisplay(self, 'upperwater_elevation', 'elevation of ocean upperwater [m]'))
     54        s += '{}\n'.format(fielddisplay(self, 'perturbation_melting_rate', 'basal melting rate perturbation added to computed melting rate (positive if melting) [m/yr]'))
    4955        s += '{}\n'.format(fielddisplay(self, 'geothermalflux', 'geothermal heat flux [W/m^2]'))
    5056        return s
    … …  
    5561        self.deepwater_melting_rate = project3d(md, 'vector', self.deepwater_melting_rate, 'type', 'node', 'layer', 1)
    5662        self.deepwater_elevation = project3d(md, 'vector', self.deepwater_elevation, 'type', 'node', 'layer', 1)
     63        self.upperwater_melting_rate = project3d(md, 'vector', self.upperwater_melting_rate, 'type', 'node', 'layer', 1)
    5764        self.upperwater_elevation = project3d(md, 'vector', self.upperwater_elevation, 'type', 'node', 'layer', 1)
    5865        self.geothermalflux = project3d(md, 'vector', self.geothermalflux, 'type', 'node', 'layer', 1) # Bedrock only gets geothermal flux
     66        self.perturbation_melting_rate = project3d(md, 'vector', self.upperwater_melting_rate, 'type', 'node', 'layer', 1)
    5967        return self
    6068    #}}}
    … …  
    7280
    7381    def checkconsistency(self, md, solution, analyses): #{{{
     82        if not np.all(np.isnan(self.perturbation_melting_rate)):
     83            md = checkfield(md, 'fieldname', 'basalforcings.perturbation_melting_rate', 'NaN', 1, 'Inf', 1, 'timeseries', 1)
    7484        if 'MasstransportAnalysis' in analyses and not solution == 'TransientSolution' and not md.transient.ismasstransport:
    7585            md = checkfield(md, 'fieldname', 'basalforcings.groundedice_melting_rate', 'NaN', 1, 'Inf', 1, 'timeseries', 1)
    7686            md = checkfield(md, 'fieldname', 'basalforcings.deepwater_melting_rate', 'NaN', 1, 'Inf', 1, 'timeseries', 1, '>=', 0)
    7787            md = checkfield(md, 'fieldname', 'basalforcings.deepwater_elevation', 'NaN', 1, 'Inf', 1, 'timeseries', 1)
     88            md = checkfield(md, 'fieldname', 'basalforcings.upperwater_melting_rate', 'NaN', 1, 'Inf', 1, 'timeseries', 1, '>=', 0)
    7889            md = checkfield(md, 'fieldname', 'basalforcings.upperwater_elevation', 'NaN', 1, 'Inf', 1, 'timeseries', 1, '<', 0)
    7990        if 'BalancethicknessAnalysis' in analyses:
    … …  
    8293            md = checkfield(md, 'fieldname', 'basalforcings.deepwater_melting_rate', '>=', 0, 'numel', 1)
    8394            md = checkfield(md, 'fieldname', 'basalforcings.deepwater_elevation', '<', 'basalforcings.upperwater_elevation', 'numel', 1)
     95            md = checkfield(md, 'fieldname', 'basalforcings.upperwater_melting_rate', '>=', 0, 'numel', 1)
    8496            md = checkfield(md, 'fieldname', 'basalforcings.upperwater_elevation', '<=', 0, 'numel', 1)
    8597        if 'ThermalAnalysis' in analyses and not solution == 'TransientSolution' and not md.transient.isthermal:
    … …  
    89101            md = checkfield(md, 'fieldname', 'basalforcings.deepwater_elevation', '<', 'basalforcings.upperwater_elevation', 'numel', 1)
    90102            md = checkfield(md, 'fieldname', 'basalforcings.upperwater_elevation', '<', 0, 'numel', 1)
     103            md = checkfield(md, 'fieldname', 'basalforcings.deepwater_melting_rate', '>=', 0, 'numel', 1)
    91104            md = checkfield(md, 'fieldname', 'basalforcings.geothermalflux', 'NaN', 1, 'Inf', 1, 'timeseries', 1, '>=', 0)
    92105        return md
    … …  
    96109        yts = md.constants.yts
    97110
    98         floatingice_melting_rate = np.zeros((md.mesh.numberofvertices, ))
    99         pos = np.where(md.geometry.base <= md.basalforcings.deepwater_elevation)[0]
    100         floatingice_melting_rate[pos] = md.basalforcings.deepwater_melting_rate[pos]
    101         pos = np.where(np.logical_and.reduce(md.geometry.base > md.basalforcings.deepwater_elevation, md.geometry.base < md.basalforcings.upperwater_elevation))
    102         floatingice_melting_rate[pos] = md.basalforcings.deepwater_melting_rate[pos] * (md.geometry.base[pos] - md.basalforcings.upperwater_elevation[pos]) / (md.basalforcings.deepwater_elevation[pos] - md.basalforcings.upperwater_elevation[pos])
    103111        WriteData(fid, prefix, 'name', 'md.basalforcings.model', 'data', 6, 'format', 'Integer')
    104         WriteData(fid, prefix, 'data', floatingice_melting_rate, 'format', 'DoubleMat', 'name', 'md.basalforcings.floatingice_melting_rate', 'mattype', 1, 'scale', 1. / yts, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
    105112        WriteData(fid, prefix, 'object', self, 'fieldname', 'groundedice_melting_rate', 'format', 'DoubleMat', 'name', 'md.basalforcings.groundedice_melting_rate', 'mattype', 1, 'scale', 1. / yts, 'timeserieslength', md.mesh.numberofvertices + 1,'yts', md.constants.yts)
    106113        WriteData(fid, prefix,'object', self, 'fieldname', 'geothermalflux', 'name', 'md.basalforcings.geothermalflux', 'format', 'DoubleMat', 'mattype', 1,'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
    107114        WriteData(fid, prefix, 'object', self, 'fieldname', 'deepwater_melting_rate', 'format', 'DoubleMat', 'name', 'md.basalforcings.deepwater_melting_rate', 'scale', 1. / yts, 'mattype', 1)
    108115        WriteData(fid, prefix, 'object', self, 'fieldname', 'deepwater_elevation', 'format', 'DoubleMat', 'name', 'md.basalforcings.deepwater_elevation', 'mattype', 1)
     116        WriteData(fid, prefix, 'object', self, 'fieldname', 'upperwater_melting_rate', 'format', 'DoubleMat', 'name', 'md.basalforcings.upperwater_melting_rate', 'scale', 1. / yts, 'mattype', 1)
    109117        WriteData(fid, prefix, 'object', self, 'fieldname', 'upperwater_elevation', 'format', 'DoubleMat', 'name', 'md.basalforcings.upperwater_elevation', 'mattype', 1)
     118        WriteData(fid, prefix, 'object', self, 'fieldname', 'perturbation_melting_rate', 'format', 'DoubleMat', 'name', 'md.basalforcings.perturbation_melting_rate', 'scale', 1. / yts, 'mattype', 1)
    110119    #}}}
Note: See TracChangeset for help on using the changeset viewer.