source: issm/oecreview/Archive/22819-23185/ISSM-22851-22852.diff

Last change on this file was 23186, checked in by Mathieu Morlighem, 7 years ago

CHG: added Archive/22819-23185
File size: 27.0 KB

  • ../trunk-jpl/src/m/classes/SMBd18opdd.m

     
    1717                ismungsm                  = 0;
    1818                isd18opd                  = 0;
    1919                issetpddfac               = 0;
    20                 istemperaturescaled       = 0;
     20                istemperaturescaled       = 1;
     21                isprecipscaled            = 1;
    2122                delta18o                  = NaN;
    2223                delta18o_surface          = NaN;
    2324                temperatures_presentday   = NaN;
    2425                precipitations_presentday = NaN;
     26                temperatures_reconstructed = NaN;
     27                precipitations_reconstructed = NaN;
    2528                pddfac_snow               = NaN;
    2629                pddfac_ice                = NaN;
    2730                requested_outputs      = {};
    … …  
    3841                function self = extrude(self,md) % {{{
    3942                        if(self.isd18opd),self.temperatures_presentday=project3d(md,'vector',self.temperatures_presentday,'type','node');end
    4043                        if(self.isd18opd),self.precipitations_presentday=project3d(md,'vector',self.precipitations_presentday,'type','node');end
     44                        if(self.istemperaturescaled==0),self.temperatures_reconstructed=project3d(md,'vector',self.temperatures_reconstructed,'type','node');end
     45                        if(self.isprecipscaled),self.precipitations_reconstructed=project3d(md,'vector',self.precipitations_reconstructed,'type','node');end
    4146                        if(self.issetpddfac), self.pddfac_snow=project3d(md,'vector',self.pddfac_snow,'type','node');end
    4247                        if(self.issetpddfac), self.pddfac_ice=project3d(md,'vector',self.pddfac_ice,'type','node');end
    4348                        self.s0p=project3d(md,'vector',self.s0p,'type','node');
    … …  
    6671                  self.ismungsm   = 0;
    6772                  self.isd18opd   = 1;
    6873                  self.istemperaturescaled = 1;
     74                  self.isprecipscaled = 1;
    6975                  self.desfac     = 0.5;
    7076                  self.rlaps      = 6.5;
    7177                  self.rlapslgm   = 6.5;
    … …  
    9298                                        md = checkfield(md,'fieldname','smb.delta18o','NaN',1,'Inf',1,'size',[2,NaN],'singletimeseries',1);
    9399                                        md = checkfield(md,'fieldname','smb.dpermil','>=',0,'numel',1);
    94100                                   md = checkfield(md,'fieldname','smb.f','>=',0,'numel',1);
     101                                        if(self.istemperaturescaled==0)
     102                                                lent=size(self.temperatures_reconstructed,2);
     103                                                multt=ceil(lent/12)*12;
     104                                                md = checkfield(md,'fieldname','smb.temperatures_reconstructed','size',[md.mesh.numberofvertices+1 multt],'NaN',1,'Inf',1,'timeseries',1);
     105                                        end
     106                                        if(self.isprecipscaled==0)
     107                                                lenp=size(self.temperatures_reconstructed,2);
     108                                                multp=ceil(lent/12)*12;
     109                                                md = checkfield(md,'fieldname','smb.precipitations_reconstructed','size',[md.mesh.numberofvertices+1 multt],'NaN',1,'Inf',1,'timeseries',1);
     110                                        end
    95111                                end
    96112                                if(self.issetpddfac==1)
    97113                                        md = checkfield(md,'fieldname','smb.pddfac_snow','>=',0,'NaN',1,'Inf',1);
    … …  
    113129                        if(self.isd18opd==1)
    114130                                fielddisplay(self,'temperatures_presentday','monthly present day surface temperatures [K], required if delta18o/mungsm/d18opd is activated');
    115131                                fielddisplay(self,'precipitations_presentday','monthly surface precipitation [m/yr water eq], required if delta18o/mungsm/d18opd is activated');
    116                                 fielddisplay(self,'istemperaturescaled','if delta18o parametrisation from present day temperature and precipitation is activated, is temperature scaled to delta18o value? (0 or 1, default is 0)');
     132                                fielddisplay(self,'istemperaturescaled','if delta18o parametrisation from present day temperature and precipitation is activated, is temperature scaled to delta18o value? (0 or 1, default is 1)');
     133                                fielddisplay(self,'isprecipscaled','if delta18o parametrisation from present day temperature and precipitation is activated, is precip scaled to delta18o value? (0 or 1, default is 1)');
     134                                if(self.istemperaturescaled==0)
     135                                        fielddisplay(self,'temperatures_reconstructed','monthly historical surface temperatures [K], required if delta18o/mungsm/d18opd is activated and istemperaturescaled is not activated');
     136                                end
     137                                if(self.isprecipscaled==0)
     138                                        fielddisplay(self,'precipitations_reconstructed','monthly historical precipitation [m/yr water eq], required if delta18o/mungsm/d18opd is activated and isprecipscaled is not activated');
     139                                end
    117140                                fielddisplay(self,'delta18o','delta18o [per mil], required if pdd is activated and d18opd activated'); 
    118141                                fielddisplay(self,'dpermil','degree per mil, required if d18opd is activated');                           
    119142                           fielddisplay(self,'f','precip/temperature scaling factor, required if d18opd is activated');
    … …  
    148171                                WriteData(fid,prefix,'object',self,'class','smb','fieldname','temperatures_presentday','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
    149172                                WriteData(fid,prefix,'object',self,'class','smb','fieldname','precipitations_presentday','format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
    150173                                WriteData(fid,prefix,'object',self,'class','smb','fieldname','istemperaturescaled','format','Boolean');
     174                                WriteData(fid,prefix,'object',self,'class','smb','fieldname','isprecipscaled','format','Boolean');
    151175                                WriteData(fid,prefix,'object',self,'class','smb','fieldname','delta18o','format','DoubleMat','mattype',1,'timeserieslength',2,'yts',md.constants.yts);
    152176                                WriteData(fid,prefix,'object',self,'class','smb','fieldname','dpermil','format','Double');
    153177                           WriteData(fid,prefix,'object',self,'class','smb','fieldname','f','format','Double');
     178                                if self.istemperaturescaled==0
     179                                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','temperatures_reconstructed','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
     180                                end
     181                                if self.isprecipscaled==0
     182                                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','precipitations_reconstructed','format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
     183                                end
    154184                        end
    155185                        if self.issetpddfac==1
    156186                                WriteData(fid,prefix,'object',self,'class','smb','fieldname','pddfac_snow','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);

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

     
    278278        SmbIssetpddfacEnum,
    279279        SmbIsshortwaveEnum,
    280280        SmbIstemperaturescaledEnum,
     281        SmbIsprecipscaledEnum,
    281282        SmbIsthermalEnum,
    282283        SmbIsturbulentfluxEnum,
    283284        SmbKEnum,
    … …  
    532533        SmbPrecipitationEnum,
    533534        SmbPrecipitationsLgmEnum,
    534535        SmbPrecipitationsPresentdayEnum,
     536        SmbPrecipitationsReconstructedEnum,
    535537        SmbReEnum,
    536538        SmbRefreezeEnum,
    537539        SmbReiniEnum,
    … …  
    543545        SmbTaEnum,
    544546        SmbTemperaturesLgmEnum,
    545547        SmbTemperaturesPresentdayEnum,
     548        SmbTemperaturesReconstructedEnum,
    546549        SmbTEnum,
    547550        SmbTeValueEnum,
    548551        SmbTiniEnum,

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

     
    286286                case SmbIssetpddfacEnum : return "SmbIssetpddfac";
    287287                case SmbIsshortwaveEnum : return "SmbIsshortwave";
    288288                case SmbIstemperaturescaledEnum : return "SmbIstemperaturescaled";
     289                case SmbIsprecipscaledEnum : return "SmbIsprecipscaled";
    289290                case SmbIsthermalEnum : return "SmbIsthermal";
    290291                case SmbIsturbulentfluxEnum : return "SmbIsturbulentflux";
    291292                case SmbKEnum : return "SmbK";
    … …  
    538539                case SmbPrecipitationEnum : return "SmbPrecipitation";
    539540                case SmbPrecipitationsLgmEnum : return "SmbPrecipitationsLgm";
    540541                case SmbPrecipitationsPresentdayEnum : return "SmbPrecipitationsPresentday";
     542                case SmbPrecipitationsReconstructedEnum : return "SmbPrecipitationsReconstructed";
    541543                case SmbReEnum : return "SmbRe";
    542544                case SmbRefreezeEnum : return "SmbRefreeze";
    543545                case SmbReiniEnum : return "SmbReini";
    … …  
    549551                case SmbTaEnum : return "SmbTa";
    550552                case SmbTemperaturesLgmEnum : return "SmbTemperaturesLgm";
    551553                case SmbTemperaturesPresentdayEnum : return "SmbTemperaturesPresentday";
     554                case SmbTemperaturesReconstructedEnum : return "SmbTemperaturesReconstructed";
    552555                case SmbTEnum : return "SmbT";
    553556                case SmbTeValueEnum : return "SmbTeValue";
    554557                case SmbTiniEnum : return "SmbTini";

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

     
    292292              else if (strcmp(name,"SmbIssetpddfac")==0) return SmbIssetpddfacEnum;
    293293              else if (strcmp(name,"SmbIsshortwave")==0) return SmbIsshortwaveEnum;
    294294              else if (strcmp(name,"SmbIstemperaturescaled")==0) return SmbIstemperaturescaledEnum;
     295              else if (strcmp(name,"SmbIsprecipscaled")==0) return SmbIsprecipscaledEnum;
    295296              else if (strcmp(name,"SmbIsthermal")==0) return SmbIsthermalEnum;
    296297              else if (strcmp(name,"SmbIsturbulentflux")==0) return SmbIsturbulentfluxEnum;
    297298              else if (strcmp(name,"SmbK")==0) return SmbKEnum;
    … …  
    381382              else if (strcmp(name,"BasalforcingsGeothermalflux")==0) return BasalforcingsGeothermalfluxEnum;
    382383              else if (strcmp(name,"BasalforcingsGroundediceMeltingRate")==0) return BasalforcingsGroundediceMeltingRateEnum;
    383384              else if (strcmp(name,"BasalforcingsPicoBasinId")==0) return BasalforcingsPicoBasinIdEnum;
    384               else if (strcmp(name,"BasalforcingsPicoBoxId")==0) return BasalforcingsPicoBoxIdEnum;
    385385         else stage=4;
    386386   }
    387387   if(stage==4){
    388               if (strcmp(name,"BasalforcingsPicoSubShelfOceanOverturning")==0) return BasalforcingsPicoSubShelfOceanOverturningEnum;
     388              if (strcmp(name,"BasalforcingsPicoBoxId")==0) return BasalforcingsPicoBoxIdEnum;
     389              else if (strcmp(name,"BasalforcingsPicoSubShelfOceanOverturning")==0) return BasalforcingsPicoSubShelfOceanOverturningEnum;
    389390              else if (strcmp(name,"BasalforcingsPicoSubShelfOceanSalinity")==0) return BasalforcingsPicoSubShelfOceanSalinityEnum;
    390391              else if (strcmp(name,"BasalforcingsPicoSubShelfOceanTemp")==0) return BasalforcingsPicoSubShelfOceanTempEnum;
    391392              else if (strcmp(name,"Base")==0) return BaseEnum;
    … …  
    504505              else if (strcmp(name,"RheologyBAbsGradient")==0) return RheologyBAbsGradientEnum;
    505506              else if (strcmp(name,"RheologyBbarAbsGradient")==0) return RheologyBbarAbsGradientEnum;
    506507              else if (strcmp(name,"Sealevel")==0) return SealevelEnum;
    507               else if (strcmp(name,"SealevelriseDeltathickness")==0) return SealevelriseDeltathicknessEnum;
    508508         else stage=5;
    509509   }
    510510   if(stage==5){
    511               if (strcmp(name,"SedimentHead")==0) return SedimentHeadEnum;
     511              if (strcmp(name,"SealevelriseDeltathickness")==0) return SealevelriseDeltathicknessEnum;
     512              else if (strcmp(name,"SedimentHead")==0) return SedimentHeadEnum;
    512513              else if (strcmp(name,"SedimentHeadOld")==0) return SedimentHeadOldEnum;
    513514              else if (strcmp(name,"SedimentHeadResidual")==0) return SedimentHeadResidualEnum;
    514515              else if (strcmp(name,"SigmaNN")==0) return SigmaNNEnum;
    … …  
    550551              else if (strcmp(name,"SmbPrecipitation")==0) return SmbPrecipitationEnum;
    551552              else if (strcmp(name,"SmbPrecipitationsLgm")==0) return SmbPrecipitationsLgmEnum;
    552553              else if (strcmp(name,"SmbPrecipitationsPresentday")==0) return SmbPrecipitationsPresentdayEnum;
     554              else if (strcmp(name,"SmbPrecipitationsReconstructed")==0) return SmbPrecipitationsReconstructedEnum;
    553555              else if (strcmp(name,"SmbRe")==0) return SmbReEnum;
    554556              else if (strcmp(name,"SmbRefreeze")==0) return SmbRefreezeEnum;
    555557              else if (strcmp(name,"SmbReini")==0) return SmbReiniEnum;
    … …  
    561563              else if (strcmp(name,"SmbTa")==0) return SmbTaEnum;
    562564              else if (strcmp(name,"SmbTemperaturesLgm")==0) return SmbTemperaturesLgmEnum;
    563565              else if (strcmp(name,"SmbTemperaturesPresentday")==0) return SmbTemperaturesPresentdayEnum;
     566              else if (strcmp(name,"SmbTemperaturesReconstructed")==0) return SmbTemperaturesReconstructedEnum;
    564567              else if (strcmp(name,"SmbT")==0) return SmbTEnum;
    565568              else if (strcmp(name,"SmbTeValue")==0) return SmbTeValueEnum;
    566569              else if (strcmp(name,"SmbTini")==0) return SmbTiniEnum;
    … …  
    625628              else if (strcmp(name,"VzSSA")==0) return VzSSAEnum;
    626629              else if (strcmp(name,"Watercolumn")==0) return WatercolumnEnum;
    627630              else if (strcmp(name,"WaterColumnOld")==0) return WaterColumnOldEnum;
    628               else if (strcmp(name,"WaterfractionDrainage")==0) return WaterfractionDrainageEnum;
    629               else if (strcmp(name,"WaterfractionDrainageIntegrated")==0) return WaterfractionDrainageIntegratedEnum;
    630               else if (strcmp(name,"Waterfraction")==0) return WaterfractionEnum;
    631631         else stage=6;
    632632   }
    633633   if(stage==6){
    634               if (strcmp(name,"Waterheight")==0) return WaterheightEnum;
     634              if (strcmp(name,"WaterfractionDrainage")==0) return WaterfractionDrainageEnum;
     635              else if (strcmp(name,"WaterfractionDrainageIntegrated")==0) return WaterfractionDrainageIntegratedEnum;
     636              else if (strcmp(name,"Waterfraction")==0) return WaterfractionEnum;
     637              else if (strcmp(name,"Waterheight")==0) return WaterheightEnum;
    635638              else if (strcmp(name,"WeightsSurfaceObservation")==0) return WeightsSurfaceObservationEnum;
    636639              else if (strcmp(name,"InputsEND")==0) return InputsENDEnum;
    637640              else if (strcmp(name,"Absolute")==0) return AbsoluteEnum;
    … …  
    748751              else if (strcmp(name,"FSSolver")==0) return FSSolverEnum;
    749752              else if (strcmp(name,"FSvelocity")==0) return FSvelocityEnum;
    750753              else if (strcmp(name,"GaussPenta")==0) return GaussPentaEnum;
    751               else if (strcmp(name,"GaussSeg")==0) return GaussSegEnum;
    752               else if (strcmp(name,"GaussTetra")==0) return GaussTetraEnum;
    753               else if (strcmp(name,"GaussTria")==0) return GaussTriaEnum;
    754754         else stage=7;
    755755   }
    756756   if(stage==7){
    757               if (strcmp(name,"GenericOption")==0) return GenericOptionEnum;
     757              if (strcmp(name,"GaussSeg")==0) return GaussSegEnum;
     758              else if (strcmp(name,"GaussTetra")==0) return GaussTetraEnum;
     759              else if (strcmp(name,"GaussTria")==0) return GaussTriaEnum;
     760              else if (strcmp(name,"GenericOption")==0) return GenericOptionEnum;
    758761              else if (strcmp(name,"GenericParam")==0) return GenericParamEnum;
    759762              else if (strcmp(name,"GiadWdt")==0) return GiadWdtEnum;
    760763              else if (strcmp(name,"GiaIvinsAnalysis")==0) return GiaIvinsAnalysisEnum;
    … …  
    871874              else if (strcmp(name,"MaxDivergence")==0) return MaxDivergenceEnum;
    872875              else if (strcmp(name,"MaxVel")==0) return MaxVelEnum;
    873876              else if (strcmp(name,"MaxVx")==0) return MaxVxEnum;
    874               else if (strcmp(name,"MaxVy")==0) return MaxVyEnum;
    875               else if (strcmp(name,"MaxVz")==0) return MaxVzEnum;
    876               else if (strcmp(name,"Melange")==0) return MelangeEnum;
    877877         else stage=8;
    878878   }
    879879   if(stage==8){
    880               if (strcmp(name,"MeltingAnalysis")==0) return MeltingAnalysisEnum;
     880              if (strcmp(name,"MaxVy")==0) return MaxVyEnum;
     881              else if (strcmp(name,"MaxVz")==0) return MaxVzEnum;
     882              else if (strcmp(name,"Melange")==0) return MelangeEnum;
     883              else if (strcmp(name,"MeltingAnalysis")==0) return MeltingAnalysisEnum;
    881884              else if (strcmp(name,"MeshElements")==0) return MeshElementsEnum;
    882885              else if (strcmp(name,"MeshX")==0) return MeshXEnum;
    883886              else if (strcmp(name,"MeshY")==0) return MeshYEnum;
    … …  
    994997              else if (strcmp(name,"Outputdefinition89")==0) return Outputdefinition89Enum;
    995998              else if (strcmp(name,"Outputdefinition8")==0) return Outputdefinition8Enum;
    996999              else if (strcmp(name,"Outputdefinition90")==0) return Outputdefinition90Enum;
    997               else if (strcmp(name,"Outputdefinition91")==0) return Outputdefinition91Enum;
    998               else if (strcmp(name,"Outputdefinition92")==0) return Outputdefinition92Enum;
    999               else if (strcmp(name,"Outputdefinition93")==0) return Outputdefinition93Enum;
    10001000         else stage=9;
    10011001   }
    10021002   if(stage==9){
    1003               if (strcmp(name,"Outputdefinition94")==0) return Outputdefinition94Enum;
     1003              if (strcmp(name,"Outputdefinition91")==0) return Outputdefinition91Enum;
     1004              else if (strcmp(name,"Outputdefinition92")==0) return Outputdefinition92Enum;
     1005              else if (strcmp(name,"Outputdefinition93")==0) return Outputdefinition93Enum;
     1006              else if (strcmp(name,"Outputdefinition94")==0) return Outputdefinition94Enum;
    10041007              else if (strcmp(name,"Outputdefinition95")==0) return Outputdefinition95Enum;
    10051008              else if (strcmp(name,"Outputdefinition96")==0) return Outputdefinition96Enum;
    10061009              else if (strcmp(name,"Outputdefinition97")==0) return Outputdefinition97Enum;
    … …  
    11171120              else if (strcmp(name,"TransientArrayParam")==0) return TransientArrayParamEnum;
    11181121              else if (strcmp(name,"TransientInput")==0) return TransientInputEnum;
    11191122              else if (strcmp(name,"TransientParam")==0) return TransientParamEnum;
    1120               else if (strcmp(name,"TransientSolution")==0) return TransientSolutionEnum;
    1121               else if (strcmp(name,"Tria")==0) return TriaEnum;
    1122               else if (strcmp(name,"TriaInput")==0) return TriaInputEnum;
    11231123         else stage=10;
    11241124   }
    11251125   if(stage==10){
    1126               if (strcmp(name,"UzawaPressureAnalysis")==0) return UzawaPressureAnalysisEnum;
     1126              if (strcmp(name,"TransientSolution")==0) return TransientSolutionEnum;
     1127              else if (strcmp(name,"Tria")==0) return TriaEnum;
     1128              else if (strcmp(name,"TriaInput")==0) return TriaInputEnum;
     1129              else if (strcmp(name,"UzawaPressureAnalysis")==0) return UzawaPressureAnalysisEnum;
    11271130              else if (strcmp(name,"VectorParam")==0) return VectorParamEnum;
    11281131              else if (strcmp(name,"Vertex")==0) return VertexEnum;
    11291132              else if (strcmp(name,"VertexPId")==0) return VertexPIdEnum;

  • ../trunk-jpl/src/c/shared/Elements/elements.h

     
    3333                                           IssmDouble* TemperaturesLgm, IssmDouble* TemperaturesPresentday,
    3434                                           IssmDouble* monthlytemperaturesout, IssmDouble* monthlyprecout);
    3535void ComputeD18OTemperaturePrecipitationFromPD(IssmDouble d018,IssmDouble dpermil,bool isTemperatureScaled,
    36                                      IssmDouble f, IssmDouble* PrecipitationPresentday,IssmDouble* TemperaturePresentday,
    37                                                IssmDouble* monthlytemperaturesout, IssmDouble* monthlyprecout); 
     36                        bool isPrecipScaled, IssmDouble f, IssmDouble* PrecipitationPresentday,IssmDouble* TemperaturePresentday,
     37                        IssmDouble* PrecipitationReconstructed,IssmDouble* TemperatureReconstructed, IssmDouble* monthlytemperaturesout,
     38                        IssmDouble* monthlyprecout);
    3839IssmDouble DrainageFunctionWaterfraction(IssmDouble waterfraction, IssmDouble dt=0.);
    3940IssmDouble StressIntensityIntegralWeight(IssmDouble depth, IssmDouble water_depth, IssmDouble thickness);
    4041

  • ../trunk-jpl/src/c/shared/Elements/ComputeD18OTemperaturePrecipitationFromPD.cpp

     
    77#include "../Numerics/numerics.h"
    88
    99void ComputeD18OTemperaturePrecipitationFromPD(IssmDouble d018,IssmDouble dpermil,bool isTemperatureScaled,
    10                                      IssmDouble f, IssmDouble* PrecipitationPresentday,IssmDouble* TemperaturePresentday,
    11                                                IssmDouble* monthlytemperaturesout, IssmDouble* monthlyprecout){
     10                        bool isPrecipScaled, IssmDouble f, IssmDouble* PrecipitationPresentday,IssmDouble* TemperaturePresentday,
     11                        IssmDouble* PrecipitationReconstructed,IssmDouble* TemperatureReconstructed, IssmDouble* monthlytemperaturesout,
     12                        IssmDouble* monthlyprecout){
    1213 
    1314  IssmDouble monthlytemperaturestmp[12],monthlyprectmp[12];
    1415  IssmDouble deltaTemp;
    … …  
    2223  for (int imonth = 0; imonth<12; imonth++){
    2324   
    2425         if (isTemperatureScaled)monthlytemperaturestmp[imonth] = TemperaturePresentday[imonth] + deltaTemp;
    25          else monthlytemperaturestmp[imonth] = TemperaturePresentday[imonth];
     26         else{
     27                 monthlytemperaturestmp[imonth] = TemperatureReconstructed[imonth];
     28                 deltaTemp=TemperatureReconstructed[imonth]-TemperaturePresentday[imonth];
     29         }
    2630
    27     monthlyprectmp[imonth] = PrecipitationPresentday[imonth]*exp((f/dpermil)*deltaTemp);
    28    
     31         if (isPrecipScaled)monthlyprectmp[imonth] = PrecipitationPresentday[imonth]*exp((f/dpermil)*deltaTemp);
     32         else monthlyprectmp[imonth] = PrecipitationReconstructed[imonth];
     33         
    2934    /*Assign output pointer*/
    3035    *(monthlytemperaturesout+imonth) = monthlytemperaturestmp[imonth];
    3136    *(monthlyprecout+imonth) = monthlyprectmp[imonth];

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

     
    534534        IssmDouble* monthlyprec=xNew<IssmDouble>(12*numvertices);
    535535        IssmDouble* TemperaturesPresentday=xNew<IssmDouble>(12*numvertices);
    536536        IssmDouble* PrecipitationsPresentday=xNew<IssmDouble>(12*numvertices);
     537        IssmDouble* TemperaturesReconstructed=xNew<IssmDouble>(12*numvertices);
     538        IssmDouble* PrecipitationsReconstructed=xNew<IssmDouble>(12*numvertices);
    537539        IssmDouble* tmp=xNew<IssmDouble>(numvertices);
    538540        IssmDouble Delta18oTime;
    539541        IssmDouble dpermil,f;
    540542        IssmDouble time,yts,time_yr,month,time_clim;
    541543        bool isTemperatureScaled=true;
     544        bool isPrecipScaled=true;
    542545        this->parameters->FindParam(&time,TimeEnum);
    543546        this->parameters->FindParam(&yts,ConstantsYtsEnum);
    544547        this->parameters->FindParam(&f,SmbFEnum);
    … …  
    548551        /*Get some pdd parameters*/
    549552        dpermil=this->matpar->GetMaterialParameter(SmbDpermilEnum);
    550553
     554        this->parameters->FindParam(&isTemperatureScaled,SmbIstemperaturescaledEnum);
     555        this->parameters->FindParam(&isPrecipScaled,SmbIsprecipscaledEnum);
     556
    551557        /*Recover present day temperature and precipitation*/
    552558        Input*     input=this->inputs->GetInput(SmbTemperaturesPresentdayEnum);    _assert_(input);
    553559        Input*     input2=this->inputs->GetInput(SmbPrecipitationsPresentdayEnum); _assert_(input2);
     560        Input*     input3=NULL;
     561        if(!isTemperatureScaled){
     562                input3=this->inputs->GetInput(SmbTemperaturesPresentdayEnum);           _assert_(input3);
     563        }
     564        Input*     input4=NULL;
     565        if(!isPrecipScaled){
     566                input4=this->inputs->GetInput(SmbPrecipitationsPresentdayEnum);         _assert_(input4);
     567        }
    554568        int offset;
    555569
    556570        offset=dynamic_cast<TransientInput*>(input)->GetTimeInputOffset(time_yr);
    … …  
    565579                        gauss->GaussVertex(iv);
    566580                        input->GetInputValue(&TemperaturesPresentday[iv*12+month],gauss,time_clim+month/12.*yts);
    567581                        input2->GetInputValue(&PrecipitationsPresentday[iv*12+month],gauss,time_clim+month/12.*yts);
     582                        PrecipitationsPresentday[iv*12+month]=PrecipitationsPresentday[iv*12+month]*yts;
    568583
    569                         PrecipitationsPresentday[iv*12+month]=PrecipitationsPresentday[iv*12+month]*yts;
     584                        if(!isTemperatureScaled){
     585                                input3->GetInputValue(&TemperaturesReconstructed[iv*12+month],gauss,time_clim+month/12.*yts);
     586                        }
     587                        if(!isPrecipScaled){
     588                                input4->GetInputValue(&PrecipitationsReconstructed[iv*12+month],gauss,time_clim+month/12.*yts);
     589                                PrecipitationsReconstructed[iv*12+month]=PrecipitationsReconstructed[iv*12+month]*yts;
     590                        }
     591
    570592                }
    571593        }
    572594
    573595        /*Recover interpolation parameters at time t*/
    574596        this->parameters->FindParam(&Delta18oTime,SmbDelta18oEnum,time);
    575         this->parameters->FindParam(&isTemperatureScaled,SmbIstemperaturescaledEnum);
    576597
    577598        /*Compute the temperature and precipitation*/
    578599        for(int iv=0;iv<numvertices;iv++){
    579                 ComputeD18OTemperaturePrecipitationFromPD(Delta18oTime,dpermil,isTemperatureScaled,f,
    580                                         &PrecipitationsPresentday[iv*12], &TemperaturesPresentday[iv*12],
     600                ComputeD18OTemperaturePrecipitationFromPD(Delta18oTime,dpermil,isTemperatureScaled,isPrecipScaled,
     601                                        f,&PrecipitationsPresentday[iv*12], &TemperaturesPresentday[iv*12],
     602                                        &PrecipitationsReconstructed[iv*12], &TemperaturesReconstructed[iv*12],
    581603                                        &monthlytemperatures[iv*12], &monthlyprec[iv*12]);
    582604        }
    583605
    … …  
    622644        xDelete<IssmDouble>(monthlyprec);
    623645        xDelete<IssmDouble>(TemperaturesPresentday);
    624646        xDelete<IssmDouble>(PrecipitationsPresentday);
     647        xDelete<IssmDouble>(TemperaturesReconstructed);
     648        xDelete<IssmDouble>(PrecipitationsReconstructed);
    625649        xDelete<IssmDouble>(tmp);
    626650       
    627651}

  • ../trunk-jpl/src/c/analyses/SmbAnalysis.cpp

     
    2020void SmbAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
    2121       
    2222        int    smb_model;
    23         bool   isdelta18o,ismungsm,isd18opd,issetpddfac;
     23        bool   isdelta18o,ismungsm,isd18opd,issetpddfac,isprecipscaled,istemperaturescaled;
    2424       
    2525        /*Update elements: */
    2626        int counter=0;
    … …  
    8282                                iomodel->FetchDataToInput(elements,"md.smb.temperatures_presentday",SmbTemperaturesPresentdayEnum);
    8383                                iomodel->FetchDataToInput(elements,"md.smb.precipitations_presentday",SmbPrecipitationsPresentdayEnum);
    8484                                iomodel->FetchDataToInput(elements,"md.smb.precipitations_lgm",SmbPrecipitationsLgmEnum);
    85                         }
    86                         else{
     85                        }else{
    8786                                iomodel->FetchDataToInput(elements,"md.smb.precipitation",SmbPrecipitationEnum);
    8887                                iomodel->FetchDataToInput(elements,"md.smb.monthlytemperatures",SmbMonthlytemperaturesEnum);
    8988                        }
    9089                        break;
    9190                case SMBd18opddEnum:
     91                        iomodel->FindConstant(&istemperaturescaled,"md.smb.istemperaturescaled");
     92                        iomodel->FindConstant(&isprecipscaled,"md.smb.isprecipscaled");
    9293                        iomodel->FindConstant(&ismungsm,"md.smb.ismungsm");
    9394                        iomodel->FindConstant(&isd18opd,"md.smb.isd18opd");
    9495                        iomodel->FindConstant(&issetpddfac,"md.smb.issetpddfac");
    … …  
    9899                        if(isd18opd){
    99100                                iomodel->FetchDataToInput(elements,"md.smb.temperatures_presentday",SmbTemperaturesPresentdayEnum);
    100101                                iomodel->FetchDataToInput(elements,"md.smb.precipitations_presentday",SmbPrecipitationsPresentdayEnum);
     102                                if(!istemperaturescaled){
     103                                        iomodel->FetchDataToInput(elements,"md.smb.temperatures_reconstructed",SmbTemperaturesReconstructedEnum);
     104                                }
     105                                if(!isprecipscaled){
     106                                        iomodel->FetchDataToInput(elements,"md.smb.precipitations_reconstructed",SmbPrecipitationsReconstructedEnum);
     107                                }
    101108                        }
    102109                        if(issetpddfac){
    103110                                iomodel->FetchDataToInput(elements,"md.smb.pddfac_snow",SmbPddfacSnowEnum,-1.);
    … …  
    219226                        if(isd18opd){
    220227                                parameters->AddObject(iomodel->CopyConstantObject("md.smb.f",SmbFEnum));
    221228                                parameters->AddObject(iomodel->CopyConstantObject("md.smb.istemperaturescaled",SmbIstemperaturescaledEnum));
     229                                parameters->AddObject(iomodel->CopyConstantObject("md.smb.isprecipscaled",SmbIsprecipscaledEnum));
    222230                                iomodel->FetchData(&temp,&N,&M,"md.smb.delta18o"); _assert_(N==2);
    223231                                parameters->AddObject(new TransientParam(SmbDelta18oEnum,&temp[0],&temp[M],interp,M));
    224232                                iomodel->DeleteData(temp,"md.smb.delta18o");
Note: See TracBrowser for help on using the repository browser.