Changeset 23906


Ignore:
Timestamp:
05/03/19 11:10:58 (6 years ago)
Author:
tpelle
Message:

NEW: Implemented local version of ISMIP6 melt rate parameterization

Location:
issm/trunk-jpl
Files:
8 edited

Legend:

Unmodified
Added
Removed

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

    r23904 r23906  
    27262726        IssmDouble  tf,gamma0,base,delta_t_basin,mean_tf_basin,absval;
    27272727        IssmDouble  basalmeltrate[NUMVERTICES];
     2728        bool        islocal;
    27282729        IssmDouble* delta_t = NULL;
    27292730        IssmDouble* mean_tf = NULL;
    … …  
    27442745        this->parameters->FindParam(&gamma0,BasalforcingsIsmp6Gamma0Enum);
    27452746        this->parameters->FindParam(&delta_t,&M,BasalforcingsIsmp6DeltaTEnum);    _assert_(M==num_basins);
    2746         this->parameters->FindParam(&mean_tf,&N,BasalforcingsIsmp6AverageTfEnum); _assert_(N==num_basins);
     2747        this->parameters->FindParam(&islocal,BasalforcingsIsmp6IsLocalEnum);
     2748        if(!islocal) {
     2749                this->parameters->FindParam(&mean_tf,&N,BasalforcingsIsmp6AverageTfEnum); _assert_(N==num_basins);
     2750        }
    27472751        Input* tf_input=this->GetInput(BasalforcingsIsmp6TfShelfEnum);            _assert_(tf_input);
    2748         delta_t_basin = delta_t[basinid]; mean_tf_basin = mean_tf[basinid];
     2752        delta_t_basin = delta_t[basinid];
     2753        if(!islocal) mean_tf_basin = mean_tf[basinid];
    27492754
    27502755        /* Compute melt rate */
    … …  
    27532758                gauss->GaussVertex(i);
    27542759                tf_input->GetInputValue(&tf,gauss);
    2755                 absval = mean_tf_basin+delta_t_basin;
    2756                 if (absval<0) {absval = -1.*absval;}
    2757                 basalmeltrate[i] = gamma0*pow((rhow*cp)/(rhoi*lf),2)*(tf+delta_t_basin)*absval;
     2760                if(!islocal) {
     2761                        absval = mean_tf_basin+delta_t_basin;
     2762                        if (absval<0) {absval = -1.*absval;}
     2763                        basalmeltrate[i] = gamma0*pow((rhow*cp)/(rhoi*lf),2)*(tf+delta_t_basin)*absval;}
     2764                else {basalmeltrate[i] = gamma0*pow((rhow*cp)/(rhoi*lf),2)*pow(max(tf+delta_t_basin,0.),2);}
    27582765        }
    27592766
    … …  
    34893496
    34903497                /*Compute ice shelf base slope (radians)*/
    3491                 alpha = atan(sqrt(slopex*slopex + slopey*slopey));
     3498           alpha = atan(sqrt(slopex*slopex + slopey*slopey));
    34923499                if(alpha>=M_PI) alpha = M_PI - 0.001;               //ensure sin(alpha) > 0 for meltrate calculations
    34933500       

  • issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.cpp

    r23799 r23906  
    7575        int         num_basins, basinid,num_depths;
    7676        IssmDouble  area, tf, base, time;
     77        bool        islocal;
    7778        IssmDouble* tf_depths = NULL;
    7879
    7980        femmodel->parameters->FindParam(&num_basins,BasalforcingsIsmp6NumBasinsEnum);
    8081        femmodel->parameters->FindParam(&tf_depths,&num_depths,BasalforcingsIsmp6TfDepthsEnum); _assert_(tf_depths);
     82        femmodel->parameters->FindParam(&islocal,BasalforcingsIsmp6IsLocalEnum);
    8183
    8284        /*Binary search works for vectors that are sorted in increasing order only, make depths positive*/
    … …  
    148150        }
    149151
    150         /*Compute sums of tf*area and shelf-area per cpu*/
    151         for(int i=0;i<femmodel->elements->Size();i++){
    152                 Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
    153                 if(!element->IsIceInElement() || !element->IsFloating()) continue;
    154                 Input* tf_input=element->GetInput(BasalforcingsIsmp6TfShelfEnum); _assert_(tf_input);
    155                 element->inputs->GetInputValue(&basinid,BasalforcingsIsmp6BasinIdEnum);
    156                 Gauss* gauss=element->NewGauss(1); gauss->GaussPoint(0);
    157                 tf_input->GetInputValue(&tf,gauss);
    158                 delete gauss;
    159                 area=element->GetHorizontalSurfaceArea();
    160                 tf_weighted_avg[basinid]+=tf*area;
    161                 areas_summed[basinid]   +=area;
     152        if(!islocal) {
     153                /*Compute sums of tf*area and shelf-area per cpu*/
     154                for(int i=0;i<femmodel->elements->Size();i++){
     155                        Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
     156                        if(!element->IsIceInElement() || !element->IsFloating()) continue;
     157                        Input* tf_input=element->GetInput(BasalforcingsIsmp6TfShelfEnum); _assert_(tf_input);
     158                        element->inputs->GetInputValue(&basinid,BasalforcingsIsmp6BasinIdEnum);
     159                        Gauss* gauss=element->NewGauss(1); gauss->GaussPoint(0);
     160                        tf_input->GetInputValue(&tf,gauss);
     161                        delete gauss;
     162                        area=element->GetHorizontalSurfaceArea();
     163                        tf_weighted_avg[basinid]+=tf*area;
     164                        areas_summed[basinid]   +=area;
     165                }
     166
     167                /*Syncronize across cpus*/
     168                ISSM_MPI_Allreduce(tf_weighted_avg,tf_weighted_avg_cpu,num_basins,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
     169                ISSM_MPI_Allreduce(areas_summed,areas_summed_cpu,num_basins,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
     170
     171                /*Compute weighted means and save*/
     172                for(int k=0;k<num_basins;k++){tf_weighted_avg_cpu[k] = tf_weighted_avg_cpu[k]/areas_summed_cpu[k];}
     173                femmodel->parameters->AddObject(new DoubleVecParam(BasalforcingsIsmp6AverageTfEnum,tf_weighted_avg_cpu,num_basins));
    162174        }
    163 
    164         /*Syncronize across cpus*/
    165         ISSM_MPI_Allreduce(tf_weighted_avg,tf_weighted_avg_cpu,num_basins,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
    166         ISSM_MPI_Allreduce(areas_summed,areas_summed_cpu,num_basins,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
    167 
    168         /*Compute weighted means and save*/
    169         for(int k=0;k<num_basins;k++){tf_weighted_avg_cpu[k] = tf_weighted_avg_cpu[k]/areas_summed_cpu[k];}
    170         femmodel->parameters->AddObject(new DoubleVecParam(BasalforcingsIsmp6AverageTfEnum,tf_weighted_avg_cpu,num_basins));
    171175
    172176        /*Cleanup and return */

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp

    r23843 r23906  
    238238                        parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.num_basins",BasalforcingsIsmp6NumBasinsEnum));
    239239                        parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.gamma_0",BasalforcingsIsmp6Gamma0Enum));
     240                   parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.islocal",BasalforcingsIsmp6IsLocalEnum));       
    240241                        iomodel->FetchData(&transparam,&M,&N,"md.basalforcings.delta_t");
    241242                        parameters->AddObject(new DoubleVecParam(BasalforcingsIsmp6DeltaTEnum,transparam,N));

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

    r23900 r23906  
    6767   BasalforcingsIsmp6DeltaTEnum,
    6868   BasalforcingsIsmp6Gamma0Enum,
     69        BasalforcingsIsmp6IsLocalEnum,
    6970   BasalforcingsIsmp6NumBasinsEnum,
    7071   BasalforcingsIsmp6TfDepthsEnum,

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

    r23900 r23906  
    7575                case BasalforcingsIsmp6DeltaTEnum : return "BasalforcingsIsmp6DeltaT";
    7676                case BasalforcingsIsmp6Gamma0Enum : return "BasalforcingsIsmp6Gamma0";
     77                case BasalforcingsIsmp6IsLocalEnum : return "BasalforcingsIsmp6IsLocal";
    7778                case BasalforcingsIsmp6NumBasinsEnum : return "BasalforcingsIsmp6NumBasins";
    7879                case BasalforcingsIsmp6TfDepthsEnum : return "BasalforcingsIsmp6TfDepths";

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

    r23900 r23906  
    7575              else if (strcmp(name,"BasalforcingsIsmp6DeltaT")==0) return BasalforcingsIsmp6DeltaTEnum;
    7676              else if (strcmp(name,"BasalforcingsIsmp6Gamma0")==0) return BasalforcingsIsmp6Gamma0Enum;
     77              else if (strcmp(name,"BasalforcingsIsmp6IsLocal")==0) return BasalforcingsIsmp6IsLocalEnum;
    7778              else if (strcmp(name,"BasalforcingsIsmp6NumBasins")==0) return BasalforcingsIsmp6NumBasinsEnum;
    7879              else if (strcmp(name,"BasalforcingsIsmp6TfDepths")==0) return BasalforcingsIsmp6TfDepthsEnum;
    … …  
    136137              else if (strcmp(name,"FemModelComm")==0) return FemModelCommEnum;
    137138              else if (strcmp(name,"FlowequationFeFS")==0) return FlowequationFeFSEnum;
    138               else if (strcmp(name,"FlowequationIsFS")==0) return FlowequationIsFSEnum;
    139139         else stage=2;
    140140   }
    141141   if(stage==2){
    142               if (strcmp(name,"FlowequationIsHO")==0) return FlowequationIsHOEnum;
     142              if (strcmp(name,"FlowequationIsFS")==0) return FlowequationIsFSEnum;
     143              else if (strcmp(name,"FlowequationIsHO")==0) return FlowequationIsHOEnum;
    143144              else if (strcmp(name,"FlowequationIsL1L2")==0) return FlowequationIsL1L2Enum;
    144145              else if (strcmp(name,"FlowequationIsSIA")==0) return FlowequationIsSIAEnum;
    … …  
    259260              else if (strcmp(name,"MaterialsTemperateiceconductivity")==0) return MaterialsTemperateiceconductivityEnum;
    260261              else if (strcmp(name,"MaterialsEffectiveconductivityAveraging")==0) return MaterialsEffectiveconductivityAveragingEnum;
    261               else if (strcmp(name,"MaterialsThermalconductivity")==0) return MaterialsThermalconductivityEnum;
    262262         else stage=3;
    263263   }
    264264   if(stage==3){
    265               if (strcmp(name,"MaterialsThermalExchangeVelocity")==0) return MaterialsThermalExchangeVelocityEnum;
     265              if (strcmp(name,"MaterialsThermalconductivity")==0) return MaterialsThermalconductivityEnum;
     266              else if (strcmp(name,"MaterialsThermalExchangeVelocity")==0) return MaterialsThermalExchangeVelocityEnum;
    266267              else if (strcmp(name,"MeltingOffset")==0) return MeltingOffsetEnum;
    267268              else if (strcmp(name,"MeshAverageVertexConnectivity")==0) return MeshAverageVertexConnectivityEnum;
    … …  
    382383              else if (strcmp(name,"SteadystateNumRequestedOutputs")==0) return SteadystateNumRequestedOutputsEnum;
    383384              else if (strcmp(name,"SteadystateReltol")==0) return SteadystateReltolEnum;
    384               else if (strcmp(name,"SteadystateRequestedOutputs")==0) return SteadystateRequestedOutputsEnum;
    385385         else stage=4;
    386386   }
    387387   if(stage==4){
    388               if (strcmp(name,"Step")==0) return StepEnum;
     388              if (strcmp(name,"SteadystateRequestedOutputs")==0) return SteadystateRequestedOutputsEnum;
     389              else if (strcmp(name,"Step")==0) return StepEnum;
    389390              else if (strcmp(name,"StressbalanceAbstol")==0) return StressbalanceAbstolEnum;
    390391              else if (strcmp(name,"StressbalanceFSreconditioning")==0) return StressbalanceFSreconditioningEnum;
    … …  
    505506              else if (strcmp(name,"EffectivePressureHydrostep")==0) return EffectivePressureHydrostepEnum;
    506507              else if (strcmp(name,"EffectivePressureStacked")==0) return EffectivePressureStackedEnum;
    507               else if (strcmp(name,"EffectivePressure")==0) return EffectivePressureEnum;
    508508         else stage=5;
    509509   }
    510510   if(stage==5){
    511               if (strcmp(name,"Enthalpy")==0) return EnthalpyEnum;
     511              if (strcmp(name,"EffectivePressure")==0) return EffectivePressureEnum;
     512              else if (strcmp(name,"Enthalpy")==0) return EnthalpyEnum;
    512513              else if (strcmp(name,"EnthalpyPicard")==0) return EnthalpyPicardEnum;
    513514              else if (strcmp(name,"EplHeadStacked")==0) return EplHeadStackedEnum;
    … …  
    628629              else if (strcmp(name,"SealevelRSLRate")==0) return SealevelRSLRateEnum;
    629630              else if (strcmp(name,"SealevelUEsa")==0) return SealevelUEsaEnum;
    630               else if (strcmp(name,"SealevelRSLEustaticRate")==0) return SealevelRSLEustaticRateEnum;
    631631         else stage=6;
    632632   }
    633633   if(stage==6){
    634               if (strcmp(name,"SealevelriseSpcthickness")==0) return SealevelriseSpcthicknessEnum;
     634              if (strcmp(name,"SealevelRSLEustaticRate")==0) return SealevelRSLEustaticRateEnum;
     635              else if (strcmp(name,"SealevelriseSpcthickness")==0) return SealevelriseSpcthicknessEnum;
    635636              else if (strcmp(name,"SealevelriseStericRate")==0) return SealevelriseStericRateEnum;
    636637              else if (strcmp(name,"SealevelNEsa")==0) return SealevelNEsaEnum;
    … …  
    751752              else if (strcmp(name,"StressTensorxy")==0) return StressTensorxyEnum;
    752753              else if (strcmp(name,"StressTensorxz")==0) return StressTensorxzEnum;
    753               else if (strcmp(name,"StressTensoryy")==0) return StressTensoryyEnum;
    754754         else stage=7;
    755755   }
    756756   if(stage==7){
    757               if (strcmp(name,"StressTensoryz")==0) return StressTensoryzEnum;
     757              if (strcmp(name,"StressTensoryy")==0) return StressTensoryyEnum;
     758              else if (strcmp(name,"StressTensoryz")==0) return StressTensoryzEnum;
    758759              else if (strcmp(name,"StressTensorzz")==0) return StressTensorzzEnum;
    759760              else if (strcmp(name,"SurfaceAbsMisfit")==0) return SurfaceAbsMisfitEnum;
    … …  
    874875              else if (strcmp(name,"Outputdefinition75")==0) return Outputdefinition75Enum;
    875876              else if (strcmp(name,"Outputdefinition76")==0) return Outputdefinition76Enum;
    876               else if (strcmp(name,"Outputdefinition77")==0) return Outputdefinition77Enum;
    877877         else stage=8;
    878878   }
    879879   if(stage==8){
    880               if (strcmp(name,"Outputdefinition78")==0) return Outputdefinition78Enum;
     880              if (strcmp(name,"Outputdefinition77")==0) return Outputdefinition77Enum;
     881              else if (strcmp(name,"Outputdefinition78")==0) return Outputdefinition78Enum;
    881882              else if (strcmp(name,"Outputdefinition79")==0) return Outputdefinition79Enum;
    882883              else if (strcmp(name,"Outputdefinition7")==0) return Outputdefinition7Enum;
    … …  
    997998              else if (strcmp(name,"FreeSurfaceTopAnalysis")==0) return FreeSurfaceTopAnalysisEnum;
    998999              else if (strcmp(name,"FrontalForcingsDefault")==0) return FrontalForcingsDefaultEnum;
    999               else if (strcmp(name,"FrontalForcingsRignot")==0) return FrontalForcingsRignotEnum;
    10001000         else stage=9;
    10011001   }
    10021002   if(stage==9){
    1003               if (strcmp(name,"FSApproximation")==0) return FSApproximationEnum;
     1003              if (strcmp(name,"FrontalForcingsRignot")==0) return FrontalForcingsRignotEnum;
     1004              else if (strcmp(name,"FSApproximation")==0) return FSApproximationEnum;
    10041005              else if (strcmp(name,"Fset")==0) return FsetEnum;
    10051006              else if (strcmp(name,"FSpressure")==0) return FSpressureEnum;
    … …  
    11201121              else if (strcmp(name,"Mumps")==0) return MumpsEnum;
    11211122              else if (strcmp(name,"Nodal")==0) return NodalEnum;
    1122               else if (strcmp(name,"Nodalvalue")==0) return NodalvalueEnum;
    11231123         else stage=10;
    11241124   }
    11251125   if(stage==10){
    1126               if (strcmp(name,"NodeSId")==0) return NodeSIdEnum;
     1126              if (strcmp(name,"Nodalvalue")==0) return NodalvalueEnum;
     1127              else if (strcmp(name,"NodeSId")==0) return NodeSIdEnum;
    11271128              else if (strcmp(name,"NoFrictionOnPartiallyFloating")==0) return NoFrictionOnPartiallyFloatingEnum;
    11281129              else if (strcmp(name,"NoMeltOnPartiallyFloating")==0) return NoMeltOnPartiallyFloatingEnum;
    … …  
    12431244              else if (strcmp(name,"VertexLId")==0) return VertexLIdEnum;
    12441245              else if (strcmp(name,"Vertices")==0) return VerticesEnum;
    1245               else if (strcmp(name,"ViscousHeating")==0) return ViscousHeatingEnum;
    12461246         else stage=11;
    12471247   }
    12481248   if(stage==11){
    1249               if (strcmp(name,"Water")==0) return WaterEnum;
     1249              if (strcmp(name,"ViscousHeating")==0) return ViscousHeatingEnum;
     1250              else if (strcmp(name,"Water")==0) return WaterEnum;
    12501251              else if (strcmp(name,"XTaylorHood")==0) return XTaylorHoodEnum;
    12511252              else if (strcmp(name,"XY")==0) return XYEnum;

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

    r23799 r23906  
    1212                tf_depths                 = NaN;
    1313                delta_t                   = NaN;
     14                islocal                   = 0;
    1415                geothermalflux            = NaN;
    1516                groundedice_melting_rate  = NaN;
    … …  
    4748                function self = setdefaultparameters(self) % {{{
    4849                        self.gamma_0 = 14477; %m/yr
     50                        self.islocal = false;
    4951                end % }}}
    5052                function md = checkconsistency(self,md,solution,analyses) % {{{
    … …  
    5557                        md = checkfield(md,'fieldname','basalforcings.tf_depths','NaN',1,'Inf',1);
    5658                        md = checkfield(md,'fieldname','basalforcings.delta_t','NaN',1,'Inf',1,'numel',md.basalforcings.num_basins);
     59                        md = checkfield(md,'fieldname','basalforcings.islocal','values',[0 1]);
    5760                        md = checkfield(md,'fieldname','basalforcings.geothermalflux','NaN',1,'Inf',1,'>=',0,'timeseries',1);
    5861                        md = checkfield(md,'fieldname','basalforcings.groundedice_melting_rate','NaN',1,'Inf',1,'timeseries',1);
    … …  
    7275                        fielddisplay(self,'tf','thermal forcing (ocean temperature minus freezing point) (degrees C)');
    7376                        fielddisplay(self,'delta_t','Ocean temperature correction per basin (degrees C)');
     77                        fielddisplay(self,'islocal','boolean to use the local version of the ISMIP6 melt rate parameterization (default false)');
    7478                        fielddisplay(self,'geothermalflux','geothermal heat flux (W/m^2)');
    7579                        fielddisplay(self,'groundedice_melting_rate','basal melting rate (positive if melting) (m/yr)');
    … …  
    8791                        WriteData(fid,prefix,'object',self,'fieldname','tf','format','MatArray','name','md.basalforcings.tf','timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
    8892                        WriteData(fid,prefix,'object',self,'fieldname','delta_t','format','DoubleMat','name','md.basalforcings.delta_t','timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
     93                        WriteData(fid,prefix,'object',self,'fieldname','islocal','format','Boolean');
    8994                        WriteData(fid,prefix,'object',self,'fieldname','geothermalflux','format','DoubleMat','name','md.basalforcings.geothermalflux','mattype',1,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts);
    9095                        WriteData(fid,prefix,'object',self,'fieldname','groundedice_melting_rate','format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);

  • issm/trunk-jpl/test/NightlyRun/test472.m

    r23799 r23906  
    2424md.basalforcings.tf_depths = [0 -1000 -2000];
    2525md.basalforcings.gamma_0 = 14477;
     26md.basalforcings.islocal = 0;
    2627
    2728%Build an artificial tf field (for times 0 and 1, 3 depth layers)
Note: See TracChangeset for help on using the changeset viewer.