Index: ../trunk-jpl/test/NightlyRun/test236.py =================================================================== --- ../trunk-jpl/test/NightlyRun/test236.py (revision 18967) +++ ../trunk-jpl/test/NightlyRun/test236.py (revision 18968) @@ -15,6 +15,7 @@ # Use of ispdd and isdelta18o methods md.surfaceforcings = SMBpdd(); md.surfaceforcings.isdelta18o=1 +md.surfaceforcings.ismungsm=0 # Add temperature, precipitation and delta18o needed to measure the surface mass balance # creating delta18o @@ -49,14 +50,20 @@ for imonth in xrange(0,12): md.surfaceforcings.precipitations_presentday[0:md.mesh.numberofvertices,imonth]=-0.4*10**(-6)*md.mesh.y+0.5 md.surfaceforcings.precipitations_presentday[md.mesh.numberofvertices,imonth]=((float(imonth)+1.)/12.) + md.surfaceforcings.precipitations_lgm[0:md.mesh.numberofvertices,imonth]=-0.4*10**(-6)*md.mesh.y+0.5 + md.surfaceforcings.precipitations_lgm[md.mesh.numberofvertices,imonth]=((float(imonth)+1.)/12.) +# Interpolation factors +md.surfaceforcings.Tdiff[1,1:md.timestepping.final_time]=0.5; +md.surfaceforcings.sealev[1,1:md.timestepping.final_time]=0.5; +# Year of each data point +md.surfaceforcings.Tdiff[2,1:md.timestepping.final_time]=1:1:md.timestepping.final_time; +md.surfaceforcings.sealev[2,1:md.timestepping.final_time]=1:1:md.timestepping.final_time; + # time steps and resolution md.timestepping.time_step=20. md.timestepping.final_time=60. -md.surfaceforcings.Pfac=[1,1] -md.surfaceforcings.Tdiff=[1,1] -md.surfaceforcings.sealev=[1,1] # md.transient.requested_outputs=['default','SurfaceforcingsMonthlytemperatures'] Index: ../trunk-jpl/test/NightlyRun/test237.py =================================================================== --- ../trunk-jpl/test/NightlyRun/test237.py (revision 18967) +++ ../trunk-jpl/test/NightlyRun/test237.py (revision 18968) @@ -11,20 +11,21 @@ md=triangle(model(),'../Exp/Square.exp',600000) #180000 md=setmask(md,'all','') +md=parameterize(md,'../Par/SquareShelf.py') # Use of ispdd and isdelta18o methods md.surfaceforcings = SMBpdd(); -md.surfaceforcings.isdelta18o=1 +md.surfaceforcings.isdelta18o=0 +md.surfaceforcings.ismungsm=1 -md=parameterize(md,'../Par/SquareShelf.py') -# Add temperature, precipitation and delta18o needed to measure the surface mass balance -# creating delta18o -delta18o=numpy.loadtxt('../Data/delta18o.data') -md.surfaceforcings.delta18o=delta18o -# creating delta18oSurface -md.surfaceforcings.delta18o_surface = numpy.zeros((2,numpy.size(delta18o,axis=1))) -md.surfaceforcings.delta18o_surface[1,:] = delta18o[1,:] +# # Add temperature, precipitation and delta18o needed to measure the surface mass balance +# # creating delta18o +# delta18o=numpy.loadtxt('../Data/delta18o.data') +# md.surfaceforcings.delta18o=delta18o +# # creating delta18oSurface +# md.surfaceforcings.delta18o_surface = numpy.zeros((2,numpy.size(delta18o,axis=1))) +# md.surfaceforcings.delta18o_surface[1,:] = delta18o[1,:] # creating Present day and lgm temperatures # Same temperature over the all region: @@ -51,16 +52,23 @@ for imonth in xrange(0,12): md.surfaceforcings.precipitations_presentday[0:md.mesh.numberofvertices,imonth]=-0.4*10**(-6)*md.mesh.y+0.5 md.surfaceforcings.precipitations_presentday[md.mesh.numberofvertices,imonth]=((float(imonth)+1.)/12.) + md.surfaceforcings.precipitations_lgm[0:md.mesh.numberofvertices,imonth]=-0.4*10**(-6)*md.mesh.y+0.5 + md.surfaceforcings.precipitations_lgm[md.mesh.numberofvertices,imonth]=((float(imonth)+1.)/12.) +# Interpolation factors +md.surfaceforcings.Pfac[1,1:md.timestepping.final_time]=0.5; +md.surfaceforcings.Tdiff[1,1:md.timestepping.final_time]=0.5; +md.surfaceforcings.sealev[1,1:md.timestepping.final_time]=0.5; +# Year of each data point +md.surfaceforcings.Pfac[2,1:md.timestepping.final_time]=1:1:md.timestepping.final_time; +md.surfaceforcings.Tdiff[2,1:md.timestepping.final_time]=1:1:md.timestepping.final_time; +md.surfaceforcings.sealev[2,1:md.timestepping.final_time]=1:1:md.timestepping.final_time; + # time steps and resolution md.timestepping.time_step=20. md.settings.output_frequency=1 md.timestepping.final_time=60. -md.surfaceforcings.Pfac=[1,1] -md.surfaceforcings.Tdiff=[1,1] -md.surfaceforcings.sealev=[1,1] - # md.transient.requested_outputs=['default','SurfaceforcingsMonthlytemperatures'] md.extrude(3,1.) Index: ../trunk-jpl/test/NightlyRun/test236.m =================================================================== --- ../trunk-jpl/test/NightlyRun/test236.m (revision 18967) +++ ../trunk-jpl/test/NightlyRun/test236.m (revision 18968) @@ -2,12 +2,18 @@ md=setmask(md,'all',''); md=parameterize(md,'../Par/SquareShelf.par'); +%md.verbose=verbose('all'); + % Use of ispdd and isdelta18o methods md.surfaceforcings = SMBpdd(); md.surfaceforcings.isdelta18o=1; +md.surfaceforcings.ismungsm=0; +%md.surfaceforcings.precipitation(1:md.mesh.numberofvertices,1:12)=0; +%md.surfaceforcings.monthlytemperatures(1:md.mesh.numberofvertices,1:12)=273; + % Add temperature, precipitation and delta18o needed to measure the surface mass balance -% creating delta18o +% creating delta18o load '../Data/delta18o.data' md.surfaceforcings.delta18o=delta18o; % creating delta18oSurface @@ -25,7 +31,8 @@ md.surfaceforcings.temperatures_lgm(md.mesh.numberofvertices+1,imonth+1)=((imonth+1)/12); end -% creating initialization and spc temperatures initialization and spc +% creating initialization and spc temperatures initialization and +% spc md.thermal.spctemperature=mean(md.surfaceforcings.temperatures_lgm(1:md.mesh.numberofvertices,1:12),2); %-10*ones(md.mesh.numberofvertices,1); md.thermal.spctemperature=repmat(md.thermal.spctemperature,1,md.timestepping.final_time/md.timestepping.time_step); itemp=0:md.timestepping.time_step:md.timestepping.final_time-md.timestepping.time_step; @@ -36,21 +43,28 @@ % creating precipitation for imonth=0:11 md.surfaceforcings.precipitations_presentday(1:md.mesh.numberofvertices,imonth+1)=-0.4*10^(-6)*md.mesh.y+0.5; + md.surfaceforcings.precipitations_lgm(1:md.mesh.numberofvertices,imonth+1)=-0.4*10^(-6)*md.mesh.y+0.5; + % Time for the last line: md.surfaceforcings.precipitations_presentday(md.mesh.numberofvertices+1,imonth+1)=((imonth+1)/12); + md.surfaceforcings.precipitations_lgm(md.mesh.numberofvertices+1,imonth+1)=((imonth+1)/12); end +% Interpolation factors +md.surfaceforcings.Tdiff(1,1:md.timestepping.final_time)=0.5; +md.surfaceforcings.sealev(1,1:md.timestepping.final_time)=0.5; +% Year of each data point +md.surfaceforcings.Tdiff(2,1:md.timestepping.final_time)=1:1:md.timestepping.final_time; +md.surfaceforcings.sealev(2,1:md.timestepping.final_time)=1:1:md.timestepping.final_time; + % time steps and resolution md.timestepping.time_step=20; +md.settings.output_frequency=1; md.timestepping.final_time=60; -md.surfaceforcings.Pfac=[1;1]; -md.surfaceforcings.Tdiff=[1;1]; -md.surfaceforcings.sealev=[1;1]; - % md.transient.requested_outputs={'default','SurfaceforcingsMonthlytemperatures'}; md=setflowequation(md,'SSA','all'); -md.cluster=generic('name',oshostname(),'np',3); +md.cluster=generic('name',oshostname(),'np',1); % 3 for the cluster md=solve(md,TransientSolutionEnum()); %Fields and tolerances to track changes Index: ../trunk-jpl/test/NightlyRun/test237.m =================================================================== --- ../trunk-jpl/test/NightlyRun/test237.m (revision 18967) +++ ../trunk-jpl/test/NightlyRun/test237.m (revision 18968) @@ -1,19 +1,27 @@ md=triangle(model(),'../Exp/Square.exp',600000.);%180000 md=setmask(md,'all',''); +md=parameterize(md,'../Par/SquareShelf.par'); -% Use of ispdd and isdelta18o methods +%md.verbose=verbose('all'); + +% Use of ispdd methods md.surfaceforcings = SMBpdd(); -md.surfaceforcings.isdelta18o=1; +md.surfaceforcings.isdelta18o=0; +md.surfaceforcings.ismungsm=1; -md=parameterize(md,'../Par/SquareShelf.par'); +if md.surfaceforcings.isdelta18o==0 & md.surfaceforcings.ismungsm==0 + md.surfaceforcings.precipitation=zeros(md.mesh.numberofvertices,1); + md.surfaceforcings.monthlytemperatures=273*ones(md.mesh.numberofvertices,1); +end + % Add temperature, precipitation and delta18o needed to measure the surface mass balance -% creating delta18o -load '../Data/delta18o.data' -md.surfaceforcings.delta18o=delta18o; -% creating delta18oSurface -md.surfaceforcings.delta18o_surface(1,1:(length(delta18o))) = 0; -md.surfaceforcings.delta18o_surface(2,:) = delta18o(2,:); +% % creating delta18o +% load '../Data/delta18o.data' +% md.surfaceforcings.delta18o=delta18o; +% % creating delta18oSurface +% md.surfaceforcings.delta18o_surface(1,1:(length(delta18o))) = 0; +% md.surfaceforcings.delta18o_surface(2,:) = delta18o(2,:); % creating Present day and lgm temperatures % Same temperature over the all region: @@ -37,12 +45,19 @@ % creating precipitation for imonth=0:11 md.surfaceforcings.precipitations_presentday(1:md.mesh.numberofvertices,imonth+1)=-0.4*10^(-6)*md.mesh.y+0.5; + md.surfaceforcings.precipitations_lgm(1:md.mesh.numberofvertices,imonth+1)=-0.4*10^(-6)*md.mesh.y+0.5; + % Time for the last line: md.surfaceforcings.precipitations_presentday(md.mesh.numberofvertices+1,imonth+1)=((imonth+1)/12); + md.surfaceforcings.precipitations_lgm(md.mesh.numberofvertices+1,imonth+1)=((imonth+1)/12); end -md.surfaceforcings.Pfac=[1;1]; -md.surfaceforcings.Tdiff=[1;1]; -md.surfaceforcings.sealev=[1;1]; +md.surfaceforcings.Pfac(1,1:md.timestepping.final_time)=0.5; +md.surfaceforcings.Tdiff(1,1:md.timestepping.final_time)=0.5; +md.surfaceforcings.sealev(1,1:md.timestepping.final_time)=0.5; +% Year of each data point +md.surfaceforcings.Pfac(2,1:md.timestepping.final_time)=1:1:md.timestepping.final_time; +md.surfaceforcings.Tdiff(2,1:md.timestepping.final_time)=1:1:md.timestepping.final_time; +md.surfaceforcings.sealev(2,1:md.timestepping.final_time)=1:1:md.timestepping.final_time; % time steps and resolution md.timestepping.time_step=20; @@ -51,10 +66,11 @@ % md.transient.requested_outputs={'default','SurfaceforcingsMonthlytemperatures'}; -md=extrude(md,3,1.); +md=extrude(md,3,1); + md=setflowequation(md,'SSA','all'); -md.cluster=generic('name',oshostname(),'np',1); -md=solve(md,TransientSolutionEnum()); +md.cluster=generic('name',oshostname(),'np',1); % 3 for the cluster +md=solve(md,TransientSolutionEnum); %Fields and tolerances to track changes field_names ={'Vx1','Vy1','Vz1','Vel1','Pressure1','Bed1','Surface1','Thickness1','Temperature1','BasalforcingsGroundediceMeltingRate1', ... Index: ../trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 18967) +++ ../trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 18968) @@ -119,7 +119,7 @@ int stabilization,finiteelement,smb_model; bool dakota_analysis; - bool isdelta18o; + bool isdelta18o,ismungsm; bool isgroundingline; bool islevelset; @@ -178,17 +178,19 @@ break; case SMBpddEnum: iomodel->Constant(&isdelta18o,SurfaceforcingsIsdelta18oEnum); + iomodel->Constant(&ismungsm,SurfaceforcingsIsmungsmEnum); iomodel->FetchDataToInput(elements,ThermalSpctemperatureEnum); - if(isdelta18o){ + if(isdelta18o || ismungsm){ iomodel->FetchDataToInput(elements,SurfaceforcingsTemperaturesLgmEnum); iomodel->FetchDataToInput(elements,SurfaceforcingsTemperaturesPresentdayEnum); iomodel->FetchDataToInput(elements,SurfaceforcingsPrecipitationsPresentdayEnum); iomodel->FetchDataToInput(elements,SurfaceforcingsPrecipitationsLgmEnum); } else{ - iomodel->FetchDataToInput(elements,SurfaceforcingsPrecipitationEnum); + iomodel->FetchDataToInput(elements,SurfaceforcingsPrecipitationEnum); iomodel->FetchDataToInput(elements,SurfaceforcingsMonthlytemperaturesEnum); } + break; case SMBgradientsEnum: iomodel->FetchDataToInput(elements,SurfaceforcingsHrefEnum); Index: ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 18967) +++ ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 18968) @@ -343,6 +343,7 @@ SurfaceforcingsDelta18oEnum, SurfaceforcingsDelta18oSurfaceEnum, SurfaceforcingsIsdelta18oEnum, + SurfaceforcingsIsmungsmEnum, SurfaceforcingsPrecipitationsPresentdayEnum, SurfaceforcingsPrecipitationsLgmEnum, SurfaceforcingsTemperaturesPresentdayEnum, Index: ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 18967) +++ ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 18968) @@ -349,6 +349,7 @@ case SurfaceforcingsDelta18oEnum : return "SurfaceforcingsDelta18o"; case SurfaceforcingsDelta18oSurfaceEnum : return "SurfaceforcingsDelta18oSurface"; case SurfaceforcingsIsdelta18oEnum : return "SurfaceforcingsIsdelta18o"; + case SurfaceforcingsIsmungsmEnum : return "SurfaceforcingsIsmungsm"; case SurfaceforcingsPrecipitationsPresentdayEnum : return "SurfaceforcingsPrecipitationsPresentday"; case SurfaceforcingsPrecipitationsLgmEnum : return "SurfaceforcingsPrecipitationsLgm"; case SurfaceforcingsTemperaturesPresentdayEnum : return "SurfaceforcingsTemperaturesPresentday"; Index: ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 18967) +++ ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 18968) @@ -355,6 +355,7 @@ else if (strcmp(name,"SurfaceforcingsDelta18o")==0) return SurfaceforcingsDelta18oEnum; else if (strcmp(name,"SurfaceforcingsDelta18oSurface")==0) return SurfaceforcingsDelta18oSurfaceEnum; else if (strcmp(name,"SurfaceforcingsIsdelta18o")==0) return SurfaceforcingsIsdelta18oEnum; + else if (strcmp(name,"SurfaceforcingsIsmungsm")==0) return SurfaceforcingsIsmungsmEnum; else if (strcmp(name,"SurfaceforcingsPrecipitationsPresentday")==0) return SurfaceforcingsPrecipitationsPresentdayEnum; else if (strcmp(name,"SurfaceforcingsPrecipitationsLgm")==0) return SurfaceforcingsPrecipitationsLgmEnum; else if (strcmp(name,"SurfaceforcingsTemperaturesPresentday")==0) return SurfaceforcingsTemperaturesPresentdayEnum; @@ -381,11 +382,11 @@ else if (strcmp(name,"SurfaceforcingsRunoff")==0) return SurfaceforcingsRunoffEnum; else if (strcmp(name,"SMBmeltcomponents")==0) return SMBmeltcomponentsEnum; else if (strcmp(name,"SurfaceforcingsMelt")==0) return SurfaceforcingsMeltEnum; - else if (strcmp(name,"SurfaceforcingsRefreeze")==0) return SurfaceforcingsRefreezeEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"SurfaceforcingsIspdd")==0) return SurfaceforcingsIspddEnum; + if (strcmp(name,"SurfaceforcingsRefreeze")==0) return SurfaceforcingsRefreezeEnum; + else if (strcmp(name,"SurfaceforcingsIspdd")==0) return SurfaceforcingsIspddEnum; else if (strcmp(name,"SurfaceforcingsIssmbgradients")==0) return SurfaceforcingsIssmbgradientsEnum; else if (strcmp(name,"SolutionType")==0) return SolutionTypeEnum; else if (strcmp(name,"AnalysisType")==0) return AnalysisTypeEnum; @@ -504,11 +505,11 @@ else if (strcmp(name,"MatrixParam")==0) return MatrixParamEnum; else if (strcmp(name,"Masscon")==0) return MassconEnum; else if (strcmp(name,"MassconName")==0) return MassconNameEnum; - else if (strcmp(name,"MassconDefinitionenum")==0) return MassconDefinitionenumEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"MassconLevelset")==0) return MassconLevelsetEnum; + if (strcmp(name,"MassconDefinitionenum")==0) return MassconDefinitionenumEnum; + else if (strcmp(name,"MassconLevelset")==0) return MassconLevelsetEnum; else if (strcmp(name,"Massconaxpby")==0) return MassconaxpbyEnum; else if (strcmp(name,"MassconaxpbyName")==0) return MassconaxpbyNameEnum; else if (strcmp(name,"MassconaxpbyDefinitionenum")==0) return MassconaxpbyDefinitionenumEnum; @@ -627,11 +628,11 @@ else if (strcmp(name,"DeviatoricStressxz")==0) return DeviatoricStressxzEnum; else if (strcmp(name,"DeviatoricStressyy")==0) return DeviatoricStressyyEnum; else if (strcmp(name,"DeviatoricStressyz")==0) return DeviatoricStressyzEnum; - else if (strcmp(name,"DeviatoricStresszz")==0) return DeviatoricStresszzEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"StrainRate")==0) return StrainRateEnum; + if (strcmp(name,"DeviatoricStresszz")==0) return DeviatoricStresszzEnum; + else if (strcmp(name,"StrainRate")==0) return StrainRateEnum; else if (strcmp(name,"StrainRatexx")==0) return StrainRatexxEnum; else if (strcmp(name,"StrainRatexy")==0) return StrainRatexyEnum; else if (strcmp(name,"StrainRatexz")==0) return StrainRatexzEnum; @@ -750,11 +751,11 @@ else if (strcmp(name,"GroundinglineMigration")==0) return GroundinglineMigrationEnum; else if (strcmp(name,"Gset")==0) return GsetEnum; else if (strcmp(name,"Index")==0) return IndexEnum; - else if (strcmp(name,"Indexed")==0) return IndexedEnum; else stage=7; } if(stage==7){ - if (strcmp(name,"Intersect")==0) return IntersectEnum; + if (strcmp(name,"Indexed")==0) return IndexedEnum; + else if (strcmp(name,"Intersect")==0) return IntersectEnum; else if (strcmp(name,"Nodal")==0) return NodalEnum; else if (strcmp(name,"OldGradient")==0) return OldGradientEnum; else if (strcmp(name,"OutputFilePointer")==0) return OutputFilePointerEnum; Index: ../trunk-jpl/src/c/shared/Elements/elements.h =================================================================== --- ../trunk-jpl/src/c/shared/Elements/elements.h (revision 18967) +++ ../trunk-jpl/src/c/shared/Elements/elements.h (revision 18968) @@ -12,18 +12,21 @@ IssmDouble Arrhenius(IssmDouble temperature,IssmDouble depth,IssmDouble n); IssmDouble LliboutryDuval(IssmDouble enthalpy, IssmDouble pressure, IssmDouble n, IssmDouble betaCC, IssmDouble referencetemperature, IssmDouble heatcapacity, IssmDouble latentheat); // IssmDouble LliboutryDuval(IssmDouble temperature, IssmDouble waterfraction, IssmDouble depth,IssmDouble n); -IssmDouble PddSurfaceMassBlance(IssmDouble* monthlytemperatures, IssmDouble* monthlyprec, - IssmDouble* TemperaturesLgm, IssmDouble* TemperaturesPresentday, - IssmDouble* PrecipitationsLgm, IssmDouble* PrecipitationsPresentday, - IssmDouble* pdds, IssmDouble* pds,IssmDouble signorm, IssmDouble yts, - IssmDouble h, IssmDouble s, IssmDouble desfac, - IssmDouble s0t,IssmDouble s0p, IssmDouble rlaps, IssmDouble rlapslgm, - IssmDouble PfacTime,IssmDouble TdiffTime,IssmDouble sealevTime); +IssmDouble PddSurfaceMassBalance(IssmDouble* monthlytemperatures, IssmDouble* monthlyprec, + IssmDouble* pdds, IssmDouble* pds,IssmDouble signorm, IssmDouble yts, + IssmDouble h, IssmDouble s, IssmDouble desfac,IssmDouble s0t, + IssmDouble s0p, IssmDouble rlaps, IssmDouble rlapslgm, + IssmDouble TdiffTime,IssmDouble sealevTime, + IssmDouble rho_water, IssmDouble rho_ice); void ComputeDelta18oTemperaturePrecipitation(IssmDouble Delta18oSurfacePresent, IssmDouble Delta18oSurfaceLgm, IssmDouble Delta18oSurfaceTime, - IssmDouble Delta18oPresent, IssmDouble Delta18oLgm, IssmDouble Delta18oTime, - IssmDouble* PrecipitationsPresentday, - IssmDouble* TemperaturesLgm, IssmDouble* TemperaturesPresentday, - IssmDouble* monthlytemperaturesout, IssmDouble* monthlyprecout); + IssmDouble Delta18oPresent, IssmDouble Delta18oLgm, IssmDouble Delta18oTime, + IssmDouble* PrecipitationsPresentday, + IssmDouble* TemperaturesLgm, IssmDouble* TemperaturesPresentday, + IssmDouble* monthlytemperaturesout, IssmDouble* monthlyprecout); +void ComputeMungsmTemperaturePrecipitation(IssmDouble TdiffTime, IssmDouble PfacTime, + IssmDouble* PrecipitationsLgm,IssmDouble* PrecipitationsPresentday, + IssmDouble* TemperaturesLgm, IssmDouble* TemperaturesPresentday, + IssmDouble* monthlytemperaturesout, IssmDouble* monthlyprecout); IssmDouble DrainageFunctionWaterfraction(IssmDouble waterfraction, IssmDouble dt=0.); IssmDouble StressIntensityIntegralWeight(IssmDouble depth, IssmDouble water_depth, IssmDouble thickness); Index: ../trunk-jpl/src/c/shared/Elements/PddSurfaceMassBalance.cpp =================================================================== --- ../trunk-jpl/src/c/shared/Elements/PddSurfaceMassBalance.cpp (revision 18967) +++ ../trunk-jpl/src/c/shared/Elements/PddSurfaceMassBalance.cpp (revision 18968) @@ -1,22 +1,18 @@ /* file: PddSurfaceMassBlance.cpp Calculating the surface mass balance using the positive degree day method. + Updating the precipitation and temperature to the new elevation */ #include "./elements.h" #include "../Numerics/numerics.h" -IssmDouble PddSurfaceMassBlance(IssmDouble* monthlytemperatures, IssmDouble* monthlyprec, - IssmDouble* TemperaturesLgm, IssmDouble* TemperaturesPresentday, - IssmDouble* PrecipitationsLgm, IssmDouble* PrecipitationsPresentday, - IssmDouble* pdds, IssmDouble* pds, IssmDouble signorm, - IssmDouble yts, IssmDouble h, IssmDouble s, IssmDouble desfac, - IssmDouble s0t, - IssmDouble s0p, IssmDouble rlaps, IssmDouble rlapslgm, - IssmDouble PfacTime,IssmDouble TdiffTime,IssmDouble sealevTime){ +IssmDouble PddSurfaceMassBalance(IssmDouble* monthlytemperatures, IssmDouble* monthlyprec, + IssmDouble* pdds, IssmDouble* pds, IssmDouble signorm, + IssmDouble yts, IssmDouble h, IssmDouble s, IssmDouble desfac, + IssmDouble s0t,IssmDouble s0p, IssmDouble rlaps,IssmDouble rlapslgm, + IssmDouble TdiffTime,IssmDouble sealevTime, + IssmDouble rho_water,IssmDouble rho_ice){ - // CURENTLY monthlytemperatures and monthlyprec ARE NOT USED HERE. - // THEY ARE REPLACE BY tdiffh No usage of deltO18 anymore - // output: IssmDouble B; // surface mass balance, melt+accumulation int iqj,imonth; @@ -27,15 +23,15 @@ IssmDouble prect; // total precipitation during 1 year taking into account des. ef. IssmDouble water; //water=rain + snowmelt IssmDouble runoff; //meltwater only, does not include rain + IssmDouble sconv; //rhow_rain/rhoi / 12 months //IssmDouble sealev=0.; // degrees per meter. 7.5 lev's 99 paper, 9 Marshall 99 paper //IssmDouble Pfac=0.5,Tdiff=0.5; - IssmDouble rtlaps; + IssmDouble rtlaps; // IssmDouble lapser=6.5 // lapse rate // IssmDouble desfac = 0.3; // desert elevation factor // IssmDouble s0p=0.; // should be set to elevation from precip source // IssmDouble s0t=0.; // should be set to elevation from temperature source - IssmDouble tdiffh; IssmDouble st; // elevation between altitude of the temp record and current altitude IssmDouble sp; // elevation between altitude of the prec record and current altitude IssmDouble deselcut=1.0; @@ -48,7 +44,6 @@ IssmDouble DT = 0.02; IssmDouble pddt, pd; // pd: snow/precip fraction, precipitation falling as snow - IssmDouble premt; // monthly precipitation (m of ice equivalent) IssmDouble q, qmpt; // q is desert/elev. fact, hnpfac is huybrect fact, and pd is normal dist. IssmDouble qm = 0.; // snow part of the precipitation IssmDouble qmt = 0.; // precipitation without desertification effect adjustment @@ -72,6 +67,7 @@ IssmDouble fsupT=0.5, fsupndd=0.5; // Tsurf mode factors for supice IssmDouble pddtj, hmx2; + sconv=(rho_water/rho_ice)/12.; //rhow_rain/rhoi / 12 months /*PDD constant*/ siglim = 2.5*signorm; @@ -89,8 +85,8 @@ rtlaps=TdiffTime*rlapslgm + (1.-TdiffTime)*rlaps; // lapse rate /*********compute Surface temperature *******/ - tdiffh = TdiffTime*( TemperaturesLgm[imonth] - TemperaturesPresentday[imonth] ); - tstar = tdiffh + TemperaturesPresentday[imonth] - rtlaps *max(st,sealevTime*0.001); + monthlytemperatures[imonth]=monthlytemperatures[imonth] - rtlaps *max(st,sealevTime*0.001); + tstar = monthlytemperatures[imonth]; Tsurf = tstar*deltm+Tsurf; /*********compute PD ****************/ @@ -105,10 +101,8 @@ if (sp>0.0){q = exp(-desfac*sp);} else {q = 1.0;} - premt =min(1.5, PrecipitationsPresentday[imonth] * pow(PrecipitationsLgm[imonth],PfacTime)); // [m/month] - - qmt= qmt + premt; - qmpt= q*premt; + qmt= qmt + monthlyprec[imonth]*sconv; //*sconv to convert in m of ice equivalent per month + qmpt= q*monthlyprec[imonth]*sconv; qmp= qmp + qmpt; qm= qm + qmpt*pd; @@ -123,6 +117,10 @@ pdd = pdd + pddsig*deltm; frzndd = frzndd - (tstar-pddsig)*deltm;} else{frzndd = frzndd - tstar*deltm; } + + /*Assign output pointer*/ + *(monthlytemperatures+imonth) = monthlytemperatures[imonth]; + *(monthlyprec+imonth) = monthlyprec[imonth]; } // end of seasonal loop //****************************************************************** Index: ../trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- ../trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 18967) +++ ../trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 18968) @@ -21,7 +21,7 @@ int numoutputs,materialtype,smb_model,basalforcing_model; char** requestedoutputs = NULL; IssmDouble time; - bool isdelta18o; + bool isdelta18o,ismungsm; /*parameters for mass flux:*/ int mass_flux_num_profiles = 0; @@ -107,27 +107,31 @@ break; case SMBpddEnum: parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsIsdelta18oEnum)); + parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsIsmungsmEnum)); parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsDesfacEnum)); parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsS0pEnum)); parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsS0tEnum)); parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsRlapsEnum)); parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsRlapslgmEnum)); iomodel->Constant(&isdelta18o,SurfaceforcingsIsdelta18oEnum); + iomodel->Constant(&ismungsm,SurfaceforcingsIsmungsmEnum); - iomodel->FetchData(&temp,&N,&M,SurfaceforcingsPfacEnum); _assert_(N==2); - for(i=0;iAddObject(new TransientParam(SurfaceforcingsPfacEnum,&temp[0],&temp[M],M)); - iomodel->DeleteData(temp,SurfaceforcingsPfacEnum); + if(ismungsm){ + iomodel->FetchData(&temp,&N,&M,SurfaceforcingsPfacEnum); _assert_(N==2); + for(i=0;iAddObject(new TransientParam(SurfaceforcingsPfacEnum,&temp[0],&temp[M],M)); + iomodel->DeleteData(temp,SurfaceforcingsPfacEnum); - iomodel->FetchData(&temp,&N,&M,SurfaceforcingsTdiffEnum); _assert_(N==2); - for(i=0;iAddObject(new TransientParam(SurfaceforcingsTdiffEnum,&temp[0],&temp[M],M)); - iomodel->DeleteData(temp,SurfaceforcingsTdiffEnum); + iomodel->FetchData(&temp,&N,&M,SurfaceforcingsTdiffEnum); _assert_(N==2); + for(i=0;iAddObject(new TransientParam(SurfaceforcingsTdiffEnum,&temp[0],&temp[M],M)); + iomodel->DeleteData(temp,SurfaceforcingsTdiffEnum); - iomodel->FetchData(&temp,&N,&M,SurfaceforcingsSealevEnum); _assert_(N==2); - for(i=0;iAddObject(new TransientParam(SurfaceforcingsSealevEnum,&temp[0],&temp[M],M)); - iomodel->DeleteData(temp,SurfaceforcingsSealevEnum); + iomodel->FetchData(&temp,&N,&M,SurfaceforcingsSealevEnum); _assert_(N==2); + for(i=0;iAddObject(new TransientParam(SurfaceforcingsSealevEnum,&temp[0],&temp[M],M)); + iomodel->DeleteData(temp,SurfaceforcingsSealevEnum); + } if(isdelta18o){ iomodel->Constant(&yts,ConstantsYtsEnum); Index: ../trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h =================================================================== --- ../trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h (revision 18967) +++ ../trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h (revision 18968) @@ -11,6 +11,7 @@ void SurfaceMassBalancex(FemModel* femmodel); void SmbGradientsx(FemModel* femmodel); void Delta18oParameterizationx(FemModel* femmodel); +void MungsmtpParameterizationx(FemModel* femmodel); void PositiveDegreeDayx(FemModel* femmodel); void SmbHenningx(FemModel* femmodel); void SmbComponentsx(FemModel* femmodel); Index: ../trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp =================================================================== --- ../trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 18967) +++ ../trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 18968) @@ -10,7 +10,7 @@ /*Intermediaties*/ int smb_model; - bool isdelta18o; + bool isdelta18o,ismungsm; /*First, get SMB model from parameters*/ femmodel->parameters->FindParam(&smb_model,SurfaceforcingsEnum); @@ -22,10 +22,15 @@ break; case SMBpddEnum: femmodel->parameters->FindParam(&isdelta18o,SurfaceforcingsIsdelta18oEnum); + femmodel->parameters->FindParam(&ismungsm,SurfaceforcingsIsmungsmEnum); if(isdelta18o){ - if(VerboseSolution()) _printf0_(" call Delta18oParametrization module\n"); + if(VerboseSolution()) _printf0_(" call Delta18oParameterization module\n"); Delta18oParameterizationx(femmodel); } + if(ismungsm){ + if(VerboseSolution()) _printf0_(" call MungsmtpParameterization module\n"); + MungsmtpParameterizationx(femmodel); + } if(VerboseSolution()) _printf0_(" call positive degree day module\n"); PositiveDegreeDayx(femmodel); break; @@ -117,6 +122,14 @@ } }/*}}}*/ +void MungsmtpParameterizationx(FemModel* femmodel){/*{{{*/ + + for(int i=0;ielements->Size();i++){ + Element* element=xDynamicCast(femmodel->elements->GetObjectByOffset(i)); + element->MungsmtpParameterization(); + } + +}/*}}}*/ void PositiveDegreeDayx(FemModel* femmodel){/*{{{*/ // void PositiveDegreeDayx(hd,vTempsea,vPrec,agd,Tsurf,ni){ @@ -143,6 +156,8 @@ IssmDouble PDup, PDCUT = 2.0; // PDcut: rain/snow cutoff temperature (C) IssmDouble tstar; // monthly mean surface temp + bool ismungsm; + IssmDouble *pdds = NULL; IssmDouble *pds = NULL; Element *element = NULL; @@ -150,6 +165,9 @@ pdds=xNew(NPDMAX+1); pds=xNew(NPDCMAX+1); + // Get ismungsm parameter + femmodel->parameters->FindParam(&ismungsm,SurfaceforcingsIsmungsmEnum); + /* initialize PDD (creation of a lookup table)*/ tstep = 0.1; tsint = tstep*0.5; @@ -204,9 +222,8 @@ for(i=0;ielements->Size();i++){ element=xDynamicCast(femmodel->elements->GetObjectByOffset(i)); - element->PositiveDegreeDay(pdds,pds,signorm); + element->PositiveDegreeDay(pdds,pds,signorm,ismungsm); } - /*free ressouces: */ xDelete(pdds); xDelete(pds); Index: ../trunk-jpl/src/c/Makefile.am =================================================================== --- ../trunk-jpl/src/c/Makefile.am (revision 18967) +++ ../trunk-jpl/src/c/Makefile.am (revision 18968) @@ -219,6 +219,7 @@ ./shared/Elements/PrintArrays.cpp\ ./shared/Elements/PddSurfaceMassBalance.cpp\ ./shared/Elements/ComputeDelta18oTemperaturePrecipitation.cpp\ + ./shared/Elements/ComputeMungsmTemperaturePrecipitation.cpp\ ./shared/Elements/DrainageFunctionWaterfraction.cpp\ ./shared/String/sharedstring.h\ ./shared/String/DescriptorIndex.cpp\ Index: ../trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Element.h (revision 18967) +++ ../trunk-jpl/src/c/classes/Elements/Element.h (revision 18968) @@ -179,6 +179,7 @@ virtual void ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index)=0; virtual void ControlToVectors(Vector* vector_control, Vector* vector_gradient,int control_enum)=0; virtual void Delta18oParameterization(void)=0; + virtual void MungsmtpParameterization(void)=0; virtual void ElementResponse(IssmDouble* presponse,int response_enum)=0; virtual void ElementSizes(IssmDouble* phx,IssmDouble* phy,IssmDouble* phz)=0; virtual int FiniteElement(void)=0; @@ -250,7 +251,7 @@ virtual void NormalTop(IssmDouble* normal,IssmDouble* xyz_list)=0; virtual int NumberofNodesPressure(void)=0; virtual int NumberofNodesVelocity(void)=0; - virtual void PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm)=0; + virtual void PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm)=0; virtual void PotentialUngrounding(Vector* potential_sheet_ungrounding)=0; virtual int PressureInterpolation()=0; virtual void ReduceMatrices(ElementMatrix* Ke,ElementVector* pe)=0; Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 18967) +++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 18968) @@ -608,7 +608,6 @@ input->GetInputValue(&TemperaturesPresentday[iv][month],gauss,month/12.*yts); input2->GetInputValue(&TemperaturesLgm[iv][month],gauss,month/12.*yts); input3->GetInputValue(&PrecipitationsPresentday[iv][month],gauss,month/12.*yts); - PrecipitationsPresentday[iv][month]=PrecipitationsPresentday[iv][month]/yts; // converion in m/sec } } @@ -650,6 +649,69 @@ /*clean-up*/ delete gauss; } +void Tria::MungsmtpParameterization(void){/*{{{*/ + /*Are we on the base? If not, return*/ + if(!IsOnBase()) return; + + int i; + IssmDouble monthlytemperatures[NUMVERTICES][12],monthlyprec[NUMVERTICES][12]; + IssmDouble TemperaturesPresentday[NUMVERTICES][12],TemperaturesLgm[NUMVERTICES][12]; + IssmDouble PrecipitationsPresentday[NUMVERTICES][12],PrecipitationsLgm[NUMVERTICES][12]; + IssmDouble tmp[NUMVERTICES]; + IssmDouble TdiffTime,PfacTime; + IssmDouble time,yts; + this->parameters->FindParam(&time,TimeEnum); + this->parameters->FindParam(&yts,ConstantsYtsEnum); + + /*Recover present day temperature and precipitation*/ + Input* input=inputs->GetInput(SurfaceforcingsTemperaturesPresentdayEnum); _assert_(input); + Input* input2=inputs->GetInput(SurfaceforcingsTemperaturesLgmEnum); _assert_(input2); + Input* input3=inputs->GetInput(SurfaceforcingsPrecipitationsPresentdayEnum); _assert_(input3); + Input* input4=inputs->GetInput(SurfaceforcingsPrecipitationsLgmEnum); _assert_(input4); + GaussTria* gauss=new GaussTria(); + for(int month=0;month<12;month++) { + for(int iv=0;ivGaussVertex(iv); + input->GetInputValue(&TemperaturesPresentday[iv][month],gauss,month/12.*yts); + input2->GetInputValue(&TemperaturesLgm[iv][month],gauss,month/12.*yts); + input3->GetInputValue(&PrecipitationsPresentday[iv][month],gauss,month/12.*yts); + input4->GetInputValue(&PrecipitationsLgm[iv][month],gauss,month/12.*yts); + } + } + + /*Recover interpolation parameters at time t*/ + this->parameters->FindParam(&TdiffTime,SurfaceforcingsTdiffEnum,time); + this->parameters->FindParam(&PfacTime,SurfaceforcingsPfacEnum,time); + + /*Compute the temperature and precipitation*/ + for(int iv=0;ivAddTimeInput(newmonthinput1,time+imonth/12.*yts); + + for(i=0;iAddTimeInput(newmonthinput2,time+imonth/12.*yts); + } + NewTemperatureInput->Configure(this->parameters); + NewPrecipitationInput->Configure(this->parameters); + + this->inputs->AddInput(NewTemperatureInput); + this->inputs->AddInput(NewPrecipitationInput); + + /*clean-up*/ + delete gauss; +} /*}}}*/ int Tria::EdgeOnBaseIndex(void){/*{{{*/ @@ -2441,83 +2503,94 @@ return TriaRef::NumberofNodes(this->VelocityInterpolation()); } /*}}}*/ -void Tria::PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm){/*{{{*/ - +void Tria::PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm){/*{{{*/ + + int i; IssmDouble agd[NUMVERTICES]; // surface mass balance IssmDouble monthlytemperatures[NUMVERTICES][12],monthlyprec[NUMVERTICES][12]; - IssmDouble TemperaturesPresentday[NUMVERTICES][12],TemperaturesLgm[NUMVERTICES][12]; - IssmDouble PrecipitationsPresentday[NUMVERTICES][12],PrecipitationsLgm[NUMVERTICES][12]; + IssmDouble tmp[NUMVERTICES]; IssmDouble h[NUMVERTICES],s[NUMVERTICES]; IssmDouble rho_water,rho_ice,desfac,s0p,s0t,rlaps,rlapslgm; IssmDouble PfacTime,TdiffTime,sealevTime; - IssmDouble sconv; //rhow_rain/rhoi / 12 months /*Get material parameters :*/ rho_water=matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum); rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum); - //rho_ice=matpar->GetRhoIce(); - //rho_water=matpar->GetRhoFreshwater(); - sconv=(rho_water/rho_ice)/12.; //to convert monbthly prec in m of ice equivalent per month + /*Get some pdd parameters*/ + desfac=matpar->GetMaterialParameter(SurfaceforcingsDesfacEnum); + s0p=matpar->GetMaterialParameter(SurfaceforcingsS0pEnum); + s0t=matpar->GetMaterialParameter(SurfaceforcingsS0tEnum); + rlaps=matpar->GetMaterialParameter(SurfaceforcingsRlapsEnum); + rlapslgm=matpar->GetMaterialParameter(SurfaceforcingsRlapslgmEnum); /*Recover monthly temperatures and precipitation and Present day and LGm ones*/ Input* input=inputs->GetInput(SurfaceforcingsMonthlytemperaturesEnum); _assert_(input); Input* input2=inputs->GetInput(SurfaceforcingsPrecipitationEnum); _assert_(input2); - Input* input3=inputs->GetInput(SurfaceforcingsTemperaturesPresentdayEnum); _assert_(input3); - Input* input4=inputs->GetInput(SurfaceforcingsTemperaturesLgmEnum); _assert_(input4); - Input* input5=inputs->GetInput(SurfaceforcingsPrecipitationsPresentdayEnum); _assert_(input5); - Input* input6=inputs->GetInput(SurfaceforcingsPrecipitationsLgmEnum); _assert_(input6); GaussTria* gauss=new GaussTria(); IssmDouble time,yts; this->parameters->FindParam(&time,TimeEnum); this->parameters->FindParam(&yts,ConstantsYtsEnum); + + for(int month=0;month<12;month++) { for(int iv=0;ivGaussVertex(iv); input->GetInputValue(&monthlytemperatures[iv][month],gauss,time+month/12.*yts); monthlytemperatures[iv][month]=monthlytemperatures[iv][month]-273.15; // conversion from Kelvin to celcius input2->GetInputValue(&monthlyprec[iv][month],gauss,time+month/12.*yts); - monthlyprec[iv][month]=monthlyprec[iv][month]*sconv*yts; // convertion in m/y - input3->GetInputValue(&TemperaturesPresentday[iv][month],gauss,month/12.*yts); - TemperaturesPresentday[iv][month]=TemperaturesPresentday[iv][month]-273.15; // conversion from Kelvin to celcius - input4->GetInputValue(&TemperaturesLgm[iv][month],gauss,month/12.*yts); - TemperaturesLgm[iv][month]=TemperaturesLgm[iv][month]-273.15; // conversion from Kelvin to celcius - input5->GetInputValue(&PrecipitationsPresentday[iv][month],gauss,month/12.*yts); - PrecipitationsPresentday[iv][month]=PrecipitationsPresentday[iv][month]; - input6->GetInputValue(&PrecipitationsLgm[iv][month],gauss,month/12.*yts); - PrecipitationsLgm[iv][month]=PrecipitationsLgm[iv][month]; } } - /*Recover Pfac, Tdiff and sealev at time t:*/ - this->parameters->FindParam(&PfacTime,SurfaceforcingsPfacEnum,time); - this->parameters->FindParam(&TdiffTime,SurfaceforcingsTdiffEnum,time); - this->parameters->FindParam(&sealevTime,SurfaceforcingsSealevEnum,time); + /*Recover Pfac, Tdiff and sealev at time t: + This parameters are used to interpolate the temperature + and precipitaton between PD and LGM when ismungsm==1 */ + if (ismungsm==1){ + this->parameters->FindParam(&TdiffTime,SurfaceforcingsTdiffEnum,time); + this->parameters->FindParam(&sealevTime,SurfaceforcingsSealevEnum,time); + } + else { + TdiffTime=0; + sealevTime=0; + } /*Recover info at the vertices: */ GetInputListOnVertices(&h[0],ThicknessEnum); GetInputListOnVertices(&s[0],SurfaceEnum); - - /*Get other pdd parameters*/ - desfac=matpar->GetMaterialParameter(SurfaceforcingsDesfacEnum); - s0p=matpar->GetMaterialParameter(SurfaceforcingsS0pEnum); - s0t=matpar->GetMaterialParameter(SurfaceforcingsS0tEnum); - rlaps=matpar->GetMaterialParameter(SurfaceforcingsRlapsEnum); - rlapslgm=matpar->GetMaterialParameter(SurfaceforcingsRlapslgmEnum); /*measure the surface mass balance*/ - for (int iv = 0; ivAddTimeInput(newmonthinput1,time+imonth/12.*yts); + // + // for(i=0;iAddTimeInput(newmonthinput2,time+imonth/12.*yts); + // } + // NewTemperatureInput->Configure(this->parameters); + // NewPrecipitationInput->Configure(this->parameters); + + this->inputs->AddInput(new TriaInput(SurfaceforcingsMassBalanceEnum,&agd[0],P1Enum)); + // this->inputs->AddInput(NewTemperatureInput); + // this->inputs->AddInput(NewPrecipitationInput); // this->inputs->AddInput(new TriaVertexInput(ThermalSpcTemperatureEnum,&Tsurf[0])); + //this->InputExtrude(SurfaceforcingsMassBalanceEnum,-1); + // this->InputExtrude(SurfaceforcingsMonthlytemperaturesEnum,-1); + // this->InputExtrude(SurfaceforcingsPrecipitationEnum,-1); + /*clean-up*/ delete gauss; } Index: ../trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 18967) +++ ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 18968) @@ -62,6 +62,7 @@ void ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index); void ControlToVectors(Vector* vector_control, Vector* vector_gradient,int control_enum); void Delta18oParameterization(void); + void MungsmtpParameterization(void); int EdgeOnBaseIndex(); void EdgeOnBaseIndices(int* pindex1,int* pindex); int EdgeOnSurfaceIndex(); @@ -108,7 +109,7 @@ int NodalValue(IssmDouble* pvalue, int index, int natureofdataenum); int NumberofNodesPressure(void); int NumberofNodesVelocity(void); - void PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm); + void PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm); void PotentialUngrounding(Vector* potential_sheet_ungrounding); int PressureInterpolation(); void ReduceMatrices(ElementMatrix* Ke,ElementVector* pe); Index: ../trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 18967) +++ ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 18968) @@ -607,7 +607,6 @@ input->GetInputValue(&TemperaturesPresentday[iv][month],gauss,month/12.*yts); input2->GetInputValue(&TemperaturesLgm[iv][month],gauss,month/12.*yts); input3->GetInputValue(&PrecipitationsPresentday[iv][month],gauss,month/12.*yts); - PrecipitationsPresentday[iv][month]=PrecipitationsPresentday[iv][month]/yts; // converion in m/sec } } @@ -652,6 +651,72 @@ /*clean-up*/ delete gauss; } +void Penta::MungsmtpParameterization(void){/*{{{*/ + /*Are we on the base? If not, return*/ + if(!IsOnBase()) return; + + int i; + IssmDouble monthlytemperatures[NUMVERTICES][12],monthlyprec[NUMVERTICES][12]; + IssmDouble TemperaturesPresentday[NUMVERTICES][12],TemperaturesLgm[NUMVERTICES][12]; + IssmDouble PrecipitationsPresentday[NUMVERTICES][12],PrecipitationsLgm[NUMVERTICES][12]; + IssmDouble tmp[NUMVERTICES]; + IssmDouble TdiffTime,PfacTime; + IssmDouble time,yts; + this->parameters->FindParam(&time,TimeEnum); + this->parameters->FindParam(&yts,ConstantsYtsEnum); + + /*Recover present day temperature and precipitation*/ + Input* input=inputs->GetInput(SurfaceforcingsTemperaturesPresentdayEnum); _assert_(input); + Input* input2=inputs->GetInput(SurfaceforcingsTemperaturesLgmEnum); _assert_(input2); + Input* input3=inputs->GetInput(SurfaceforcingsPrecipitationsPresentdayEnum); _assert_(input3); + Input* input4=inputs->GetInput(SurfaceforcingsPrecipitationsLgmEnum); _assert_(input4); + GaussPenta* gauss=new GaussPenta(); + for(int month=0;month<12;month++) { + for(int iv=0;ivGaussVertex(iv); + input->GetInputValue(&TemperaturesPresentday[iv][month],gauss,month/12.*yts); + input2->GetInputValue(&TemperaturesLgm[iv][month],gauss,month/12.*yts); + input3->GetInputValue(&PrecipitationsPresentday[iv][month],gauss,month/12.*yts); + input4->GetInputValue(&PrecipitationsLgm[iv][month],gauss,month/12.*yts); + } + } + + /*Recover interpolation parameters at time t*/ + this->parameters->FindParam(&TdiffTime,SurfaceforcingsTdiffEnum,time); + this->parameters->FindParam(&PfacTime,SurfaceforcingsPfacEnum,time); + + /*Compute the temperature and precipitation*/ + for(int iv=0;ivAddTimeInput(newmonthinput1,time+imonth/12.*yts); + + for(i=0;iAddTimeInput(newmonthinput2,time+imonth/12.*yts); + } + NewTemperatureInput->Configure(this->parameters); + NewPrecipitationInput->Configure(this->parameters); + + this->inputs->AddInput(NewTemperatureInput); + this->inputs->AddInput(NewPrecipitationInput); + + this->InputExtrude(SurfaceforcingsMonthlytemperaturesEnum,-1); + this->InputExtrude(SurfaceforcingsPrecipitationEnum,-1); + + /*clean-up*/ + delete gauss; +} /*}}}*/ void Penta::ElementResponse(IssmDouble* presponse,int response_enum){/*{{{*/ @@ -2081,32 +2146,30 @@ } /*}}}*/ -void Penta::PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm){/*{{{*/ +void Penta::PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm){/*{{{*/ + int i; IssmDouble agd[NUMVERTICES]; // surface mass balance IssmDouble monthlytemperatures[NUMVERTICES][12],monthlyprec[NUMVERTICES][12]; - IssmDouble TemperaturesPresentday[NUMVERTICES][12],TemperaturesLgm[NUMVERTICES][12]; - IssmDouble PrecipitationsPresentday[NUMVERTICES][12],PrecipitationsLgm[NUMVERTICES][12]; + IssmDouble tmp[NUMVERTICES]; IssmDouble h[NUMVERTICES],s[NUMVERTICES]; IssmDouble rho_water,rho_ice,desfac,s0p,s0t,rlaps,rlapslgm; IssmDouble PfacTime,TdiffTime,sealevTime; - IssmDouble sconv; //rhow_rain/rhoi / 12 months /*Get material parameters :*/ rho_water=matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum); rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum); - //rho_ice=matpar->GetRhoIce(); - //rho_water=matpar->GetRhoFreshwater(); - sconv=(rho_water/rho_ice)/12.; //to convert monbthly prec in m of ice equivalent per month + /*Get some pdd parameters*/ + desfac=matpar->GetMaterialParameter(SurfaceforcingsDesfacEnum); + s0p=matpar->GetMaterialParameter(SurfaceforcingsS0pEnum); + s0t=matpar->GetMaterialParameter(SurfaceforcingsS0tEnum); + rlaps=matpar->GetMaterialParameter(SurfaceforcingsRlapsEnum); + rlapslgm=matpar->GetMaterialParameter(SurfaceforcingsRlapslgmEnum); /*Recover monthly temperatures and precipitation*/ Input* input=inputs->GetInput(SurfaceforcingsMonthlytemperaturesEnum); _assert_(input); Input* input2=inputs->GetInput(SurfaceforcingsPrecipitationEnum); _assert_(input2); - Input* input3=inputs->GetInput(SurfaceforcingsTemperaturesPresentdayEnum); _assert_(input3); - Input* input4=inputs->GetInput(SurfaceforcingsTemperaturesLgmEnum); _assert_(input4); - Input* input5=inputs->GetInput(SurfaceforcingsPrecipitationsPresentdayEnum); _assert_(input5); - Input* input6=inputs->GetInput(SurfaceforcingsPrecipitationsLgmEnum); _assert_(input6); GaussPenta* gauss=new GaussPenta(); IssmDouble time,yts; this->parameters->FindParam(&time,TimeEnum); @@ -2118,50 +2181,61 @@ input->GetInputValue(&monthlytemperatures[iv][month],gauss,time+month/12.*yts); monthlytemperatures[iv][month]=monthlytemperatures[iv][month]-273.15; // conversion from Kelvin to celcius input2->GetInputValue(&monthlyprec[iv][month],gauss,time+month/12.*yts); - monthlyprec[iv][month]=monthlyprec[iv][month]*sconv; // convertion to m/yr - input3->GetInputValue(&TemperaturesPresentday[iv][month],gauss,month/12.*yts); - TemperaturesPresentday[iv][month]=TemperaturesPresentday[iv][month]-273.15; // conversion from Kelvin to celcius - input4->GetInputValue(&TemperaturesLgm[iv][month],gauss,month/12.*yts); - TemperaturesLgm[iv][month]=TemperaturesLgm[iv][month]-273.15; // conversion from Kelvin to celcius - input5->GetInputValue(&PrecipitationsPresentday[iv][month],gauss,month/12.*yts); - PrecipitationsPresentday[iv][month]=PrecipitationsPresentday[iv][month]; - input6->GetInputValue(&PrecipitationsLgm[iv][month],gauss,month/12.*yts); - PrecipitationsLgm[iv][month]=PrecipitationsLgm[iv][month]; } } - /*Recover Pfac, Tdiff and sealev at time t:*/ - this->parameters->FindParam(&PfacTime,SurfaceforcingsPfacEnum,time); - this->parameters->FindParam(&TdiffTime,SurfaceforcingsTdiffEnum,time); - this->parameters->FindParam(&sealevTime,SurfaceforcingsSealevEnum,time); + /*Recover Pfac, Tdiff and sealev at time t: + This parameters are used to interpolate the temperature + and precipitaton between PD and LGM when ismungsm==1 */ + if (ismungsm==1){ + this->parameters->FindParam(&TdiffTime,SurfaceforcingsTdiffEnum,time); + this->parameters->FindParam(&sealevTime,SurfaceforcingsSealevEnum,time); + } + else { + TdiffTime=0; + sealevTime=0; + } /*Recover info at the vertices: */ GetInputListOnVertices(&h[0],ThicknessEnum); GetInputListOnVertices(&s[0],SurfaceEnum); - /*Get other pdd parameters*/ - desfac=matpar->GetMaterialParameter(SurfaceforcingsDesfacEnum); - s0p=matpar->GetMaterialParameter(SurfaceforcingsS0pEnum); - s0t=matpar->GetMaterialParameter(SurfaceforcingsS0tEnum); - rlaps=matpar->GetMaterialParameter(SurfaceforcingsRlapsEnum); - rlapslgm=matpar->GetMaterialParameter(SurfaceforcingsRlapslgmEnum); /*measure the surface mass balance*/ for (int iv = 0; iv < NUMVERTICES; iv++){ - agd[iv]=PddSurfaceMassBlance(&monthlytemperatures[iv][0], &monthlyprec[iv][0], - &TemperaturesLgm[iv][0], &TemperaturesPresentday[iv][0], - &PrecipitationsLgm[iv][0], &PrecipitationsPresentday[iv][0], + agd[iv]=PddSurfaceMassBalance(&monthlytemperatures[iv][0], &monthlyprec[iv][0], pdds,pds, signorm, yts, h[iv], s[iv], - desfac, s0t, s0p,rlaps,rlapslgm,PfacTime,TdiffTime,sealevTime); + desfac, s0t, s0p,rlaps,rlapslgm,TdiffTime,sealevTime, + rho_water,rho_ice); } /*Update inputs*/ + // TransientInput* NewTemperatureInput = new TransientInput(SurfaceforcingsMonthlytemperaturesEnum); + // TransientInput* NewPrecipitationInput = new TransientInput(SurfaceforcingsPrecipitationEnum); + // for (int imonth=0;imonth<12;imonth++) { + // for(i=0;iAddTimeInput(newmonthinput1,time+imonth/12.*yts); + // + // for(i=0;iAddTimeInput(newmonthinput2,time+imonth/12.*yts); + // } + // NewTemperatureInput->Configure(this->parameters); + // NewPrecipitationInput->Configure(this->parameters); + + + this->inputs->AddInput(new PentaInput(SurfaceforcingsMassBalanceEnum,&agd[0],P1Enum)); - //this->inputs->AddInput(new PentaVertexInput(ThermalSpcTemperatureEnum,&Tsurf[0])); - this->InputExtrude(SurfaceforcingsMassBalanceEnum,-1); + // this->inputs->AddInput(NewTemperatureInput); + // this->inputs->AddInput(NewPrecipitationInput); + // //this->inputs->AddInput(new PentaVertexInput(ThermalSpcTemperatureEnum,&Tsurf[0])); + this->InputExtrude(SurfaceforcingsMassBalanceEnum,-1); + // this->InputExtrude(SurfaceforcingsMonthlytemperaturesEnum,-1); + // this->InputExtrude(SurfaceforcingsPrecipitationEnum,-1); - /*clean-up*/ - delete gauss; + /*clean-up*/ + delete gauss; } /*}}}*/ void Penta::PotentialUngrounding(Vector* potential_ungrounding){/*{{{*/ Index: ../trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Penta.h (revision 18967) +++ ../trunk-jpl/src/c/classes/Elements/Penta.h (revision 18968) @@ -58,6 +58,7 @@ void ControlToVectors(Vector* vector_control, Vector* vector_gradient,int control_enum); ElementMatrix* CreateBasalMassMatrix(void); void Delta18oParameterization(void); + void MungsmtpParameterization(void); void ElementResponse(IssmDouble* presponse,int response_enum); void ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz); int FiniteElement(void); @@ -136,7 +137,7 @@ int NodalValue(IssmDouble* pvalue, int index, int natureofdataenum); int NumberofNodesPressure(void); int NumberofNodesVelocity(void); - void PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm); + void PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm); void PotentialUngrounding(Vector* potential_sheet_ungrounding); int PressureInterpolation(); void ReduceMatrices(ElementMatrix* Ke,ElementVector* pe); Index: ../trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Seg.h (revision 18967) +++ ../trunk-jpl/src/c/classes/Elements/Seg.h (revision 18968) @@ -53,6 +53,7 @@ void ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index){_error_("not implemented yet");}; void ControlToVectors(Vector* vector_control, Vector* vector_gradient,int control_enum){_error_("not implemented yet");}; void Delta18oParameterization(void){_error_("not implemented yet");}; + void MungsmtpParameterization(void){_error_("not implemented yet");}; void ElementResponse(IssmDouble* presponse,int response_enum){_error_("not implemented yet");}; void ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){_error_("not implemented yet");}; void FSContactMigration(Vector* vertexgrounded,Vector* vertexfloating){_error_("not implemented yet");}; @@ -128,7 +129,7 @@ void NormalTop(IssmDouble* normal,IssmDouble* xyz_list){_error_("not implemented yet");}; int NumberofNodesPressure(void){_error_("not implemented yet");}; int NumberofNodesVelocity(void){_error_("not implemented yet");}; - void PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm){_error_("not implemented yet");}; + void PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm){_error_("not implemented yet");}; void PotentialUngrounding(Vector* potential_sheet_ungrounding){_error_("not implemented yet");}; int PressureInterpolation(void){_error_("not implemented yet");}; void ReduceMatrices(ElementMatrix* Ke,ElementVector* pe){_error_("not implemented yet");}; Index: ../trunk-jpl/src/c/classes/Elements/Tetra.h =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Tetra.h (revision 18967) +++ ../trunk-jpl/src/c/classes/Elements/Tetra.h (revision 18968) @@ -53,6 +53,7 @@ void ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index){_error_("not implemented yet");}; void ControlToVectors(Vector* vector_control, Vector* vector_gradient,int control_enum){_error_("not implemented yet");}; void Delta18oParameterization(void){_error_("not implemented yet");}; + void MungsmtpParameterization(void){_error_("not implemented yet");}; IssmDouble DragCoefficientAbsGradient(void){_error_("not implemented yet");}; void ElementResponse(IssmDouble* presponse,int response_enum){_error_("not implemented yet");}; void ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz); @@ -134,7 +135,7 @@ void NormalTop(IssmDouble* normal,IssmDouble* xyz_list); int NumberofNodesPressure(void); int NumberofNodesVelocity(void); - void PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm){_error_("not implemented yet");}; + void PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm){_error_("not implemented yet");}; void PotentialUngrounding(Vector* potential_sheet_ungrounding){_error_("not implemented yet");}; int PressureInterpolation(void); void ResetFSBasalBoundaryCondition(void); Index: ../trunk-jpl/src/m/classes/SMBpdd.py =================================================================== --- ../trunk-jpl/src/m/classes/SMBpdd.py (revision 18967) +++ ../trunk-jpl/src/m/classes/SMBpdd.py (revision 18968) @@ -25,11 +25,13 @@ self.Tdiff = float('NaN') self.sealev = float('NaN') self.isdelta18o = 0 + self.ismungsm = 0 self.delta18o = float('NaN') self.delta18o_surface = float('NaN') self.temperatures_presentday = float('NaN') self.temperatures_lgm = float('NaN') self.precipitations_presentday = float('NaN') + self.precipitations_lgm = float('NaN') #set defaults self.setdefaultparameters() @@ -37,46 +39,63 @@ def __repr__(self): # {{{ string=" surface forcings parameters:" - string="%s\n%s"%(string,fielddisplay(self,'precipitation','surface precipitation [m/yr water eq]')) + string="%s\n%s"%(string,fielddisplay(self,'isdelta18o','is temperature and precipitation delta18o parametrisation activated (0 or 1, default is 0)')) + string="%s\n%s"%(string,fielddisplay(self,'ismungsm','is temperature and precipitation mungsm parametrisation activated (0 or 1, default is 0)')) string="%s\n%s"%(string,fielddisplay(self,'desfac','desertification elevation factor (between 0 and 1, default is 0.5) [m]')) - string="%s\n%s"%(string,fielddisplay(self,'isdelta18o','is temperature and precipitation delta18o parametrisation activated (0 or 1, default is 0)')) string="%s\n%s"%(string,fielddisplay(self,'s0p','should be set to elevation from precip source (between 0 and a few 1000s m, default is 0) [m]')) string="%s\n%s"%(string,fielddisplay(self,'s0t','should be set to elevation from temperature source (between 0 and a few 1000s m, default is 0) [m]')) string="%s\n%s"%(string,fielddisplay(self,'rlaps','present day lapse rate [degree/km]')) string="%s\n%s"%(string,fielddisplay(self,'rlapslgm','LGM lapse rate [degree/km]')) - string="%s\n%s"%(string,fielddisplay(self,'Pfac','time interpolation parameter for precipitation, 1D (year)')) - string="%s\n%s"%(string,fielddisplay(self,'Tdiff','time interpolation parameter for temperature, 1D (year)')) - string="%s\n%s"%(string,fielddisplay(self,'sealev','sea level [m], 1D(year)')) - string="%s\n%s"%(string,fielddisplay(self,'monthlytemperatures','monthly surface temperatures [K], required if pdd is activated and delta18o not activated')) - string="%s\n%s"%(string,fielddisplay(self,'temperatures_presentday','monthly present day surface temperatures [K], required if pdd is activated and delta18o activated')) - string="%s\n%s"%(string,fielddisplay(self,'temperatures_lgm','monthly LGM surface temperatures [K], required if pdd is activated and delta18o activated')) - string="%s\n%s"%(string,fielddisplay(self,'precipitations_presentday','monthly surface precipitation [m/yr water eq], required if pdd is activated and delta18o activated')) - string="%s\n%s"%(string,fielddisplay(self,'delta18o','delta18o, required if pdd is activated and delta18o activated')) - string="%s\n%s"%(string,fielddisplay(self,'delta18o_surface','surface elevation of the delta18o site, required if pdd is activated and delta18o activated')) - + if not (self.isdelta18o and self.ismungsm): + string="%s\n%s"%(string,fielddisplay(self,'monthlytemperatures',['monthly surface temperatures [K], required if pdd is activated and delta18o not activated'])) + string="%s\n%s"%(string,fielddisplay(self,'precipitation',['monthly surface precipitation [m/yr water eq], required if pdd is activated and delta18o or mungsm not activated'])) + if self.isdelta18o: + string="%s\n%s"%(string,fielddisplay(self,'delta18o','delta18o, required if pdd is activated and delta18o activated')) + string="%s\n%s"%(string,fielddisplay(self,'delta18o_surface','surface elevation of the delta18o site, required if pdd is activated and delta18o activated')) + string="%s\n%s"%(string,fielddisplay(self,'temperatures_presentday','monthly present day surface temperatures [K], required if delta18o/mungsm is activated')) + string="%s\n%s"%(string,fielddisplay(self,'temperatures_lgm','monthly LGM surface temperatures [K], required if delta18o or mungsm is activated')) + string="%s\n%s"%(string,fielddisplay(self,'precipitations_presentday','monthly surface precipitation [m/yr water eq], required if delta18o or mungsm is activated')) + string="%s\n%s"%(string,fielddisplay(self,'precipitations_lgm','monthly surface precipitation [m/yr water eq], required if delta18o or mungsm is activated')) + string="%s\n%s"%(string,fielddisplay(self,'Tdiff','time interpolation parameter for temperature, 1D(year), required if mungsm is activated')) + string="%s\n%s"%(string,fielddisplay(self,'sealev','sea level [m], 1D(year), required if mungsm is activated')) + if self.ismungsm: + string="%s\n%s"%(string,fielddisplay(self,'temperatures_presentday','monthly present day surface temperatures [K], required if delta18o/mungsm is activated')) + string="%s\n%s"%(string,fielddisplay(self,'temperatures_lgm','monthly LGM surface temperatures [K], required if delta18o or mungsm is activated')) + string="%s\n%s"%(string,fielddisplay(self,'precipitations_presentday','monthly surface precipitation [m/yr water eq], required if delta18o or mungsm is activated')) + string="%s\n%s"%(string,fielddisplay(self,'precipitations_lgm','monthly surface precipitation [m/yr water eq], required if delta18o or mungsm is activated')) + string="%s\n%s"%(string,fielddisplay(self,'Pfac','time interpolation parameter for precipitation, 1D(year), required if mungsm is activated')) + string="%s\n%s"%(string,fielddisplay(self,'Tdiff','time interpolation parameter for temperature, 1D(year), required if mungsm is activated')) + string="%s\n%s"%(string,fielddisplay(self,'sealev','sea level [m], 1D(year), required if mungsm is activated')) return string #}}} def extrude(self,md): # {{{ - self.precipitation=project3d(md,'vector',self.precipitation,'type','node'); - self.monthlytemperatures=project3d(md,'vector',self.monthlytemperatures,'type','node'); + if not (self.isdelta18o and self.ismungsm): + self.precipitation=project3d(md,'vector',self.precipitation,'type','node') + self.monthlytemperatures=project3d(md,'vector',self.monthlytemperatures,'type','node') if self.isdelta18o: self.temperatures_lgm=project3d(md,'vector',self.temperatures_lgm,'type','node') if self.isdelta18o: self.temperatures_presentday=project3d(md,'vector',self.temperatures_presentday,'type','node') if self.isdelta18o: self.precipitations_presentday=project3d(md,'vector',self.precipitations_presentday,'type','node') + if self.isdelta18o: self.precipitations_lgm=project3d(md,'vector',self.precipitations_lgm,'type','node') + if self.ismungsm: self.temperatures_lgm=project3d(md,'vector',self.temperatures_lgm,'type','node') + if self.ismungsm: self.temperatures_presentday=project3d(md,'vector',self.temperatures_presentday,'type','node') + if self.ismungsm: self.precipitations_presentday=project3d(md,'vector',self.precipitations_presentday,'type','node') + if self.ismungsm: self.precipitations_lgm=project3d(md,'vector',self.precipitations_lgm,'type','node') return self #}}} def initialize(self,md): # {{{ - if numpy.all(numpy.isnan(self.precipitation)): - self.precipitation=numpy.zeros((md.mesh.numberofvertices,1)) - print " no SMBpdd.precipitation specified: values set as zero" - - return self + # if numpy.all(numpy.isnan(self.precipitation)): + # self.precipitation=numpy.zeros((md.mesh.numberofvertices,1)) + # print " no SMBpdd.precipitation specified: values set as zero" + # + # return self #}}} def setdefaultparameters(self): # {{{ #pdd method not used in default mode self.isdelta18o = 0 + self.ismungsm = 0 self.desfac = 0.5 self.s0p = 0. self.s0t = 0. @@ -88,23 +107,32 @@ def checkconsistency(self,md,solution,analyses): # {{{ if MasstransportAnalysisEnum() in analyses: - md = checkfield(md,'fieldname','surfaceforcings.desfac','<=',1,'numel',[1]); - md = checkfield(md,'fieldname','surfaceforcings.s0p','>=',0,'numel',[1]); - md = checkfield(md,'fieldname','surfaceforcings.s0t','>=',0,'numel',[1]); - md = checkfield(md,'fieldname','surfaceforcings.rlaps','>=',0,'numel',[1]); - md = checkfield(md,'fieldname','surfaceforcings.rlapslgm','>=',0,'numel',[1]); - md = checkfield(md,'fieldname','surfaceforcings.Pfac','NaN',1,'size',[2,numpy.nan]); - md = checkfield(md,'fieldname','surfaceforcings.Tdiff','NaN',1,'size',[2,numpy.nan]); - md = checkfield(md,'fieldname','surfaceforcings.sealev','NaN',1,'size',[2,numpy.nan]); - if not self.isdelta18o: - md = checkfield(md,'fieldname','surfaceforcings.monthlytemperatures','forcing',1,'NaN',1) - md = checkfield(md,'fieldname','surfaceforcings.precipitation','forcing',1,'NaN',1) - else: + md = checkfield(md,'fieldname','surfaceforcings.desfac','<=',1,'numel',[1]) + md = checkfield(md,'fieldname','surfaceforcings.s0p','>=',0,'numel',[1]) + md = checkfield(md,'fieldname','surfaceforcings.s0t','>=',0,'numel',[1]) + md = checkfield(md,'fieldname','surfaceforcings.rlaps','>=',0,'numel',[1]) + md = checkfield(md,'fieldname','surfaceforcings.rlapslgm','>=',0,'numel',[1]) + + if not (self.isdelta18o and self.ismungsm): + md = checkfield(md,'fieldname','surfaceforcings.monthlytemperatures','NaN',1) + md = checkfield(md,'fieldname','surfaceforcings.precipitation','NaN',1) + elif self.isdelta18o: md = checkfield(md,'fieldname','surfaceforcings.delta18o','NaN',1) md = checkfield(md,'fieldname','surfaceforcings.delta18o_surface','NaN',1) md = checkfield(md,'fieldname','surfaceforcings.temperatures_presentday','size',[md.mesh.numberofvertices+1,12],'NaN',1) md = checkfield(md,'fieldname','surfaceforcings.temperatures_lgm','size',[md.mesh.numberofvertices+1,12],'NaN',1) md = checkfield(md,'fieldname','surfaceforcings.precipitations_presentday','size',[md.mesh.numberofvertices+1,12],'NaN',1) + md = checkfield(md,'fieldname','surfaceforcings.precipitations_lgm','size',[md.mesh.numberofvertices+1 12],'NaN',1) + md = checkfield(md,'fieldname','surfaceforcings.Tdiff','NaN',1,'size',[2,numpy.nan]) + md = checkfield(md,'fieldname','surfaceforcings.sealev','NaN',1,'size',[2,numpy.nan]) + elif self.ismungsm: + md = checkfield(md,'fieldname','surfaceforcings.temperatures_presentday','size',[md.mesh.numberofvertices+1 12],'NaN',1) + md = checkfield(md,'fieldname','surfaceforcings.temperatures_lgm','size',[md.mesh.numberofvertices+1 12],'NaN',1) + md = checkfield(md,'fieldname','surfaceforcings.precipitations_presentday','size',[md.mesh.numberofvertices+1 12],'NaN',1) + md = checkfield(md,'fieldname','surfaceforcings.precipitations_lgm','size',[md.mesh.numberofvertices+1 12],'NaN',1) + md = checkfield(md,'fieldname','surfaceforcings.Pfac','NaN',1,'size',[2,numpy.nan]) + md = checkfield(md,'fieldname','surfaceforcings.Tdiff','NaN',1,'size',[2,numpy.nan]) + md = checkfield(md,'fieldname','surfaceforcings.sealev','NaN',1,'size',[2,numpy.nan]) return md # }}} @@ -112,24 +140,34 @@ yts=365.0*24.0*3600.0 - WriteData(fid,'enum',SurfaceforcingsEnum(),'data',SMBpddEnum(),'format','Integer'); + WriteData(fid,'enum',SurfaceforcingsEnum(),'data',SMBpddEnum(),'format','Integer') - WriteData(fid,'object',self,'class','surfaceforcings','fieldname','precipitation','format','DoubleMat','mattype',1,'scale',1./yts,'forcinglength',md.mesh.numberofvertices+1) WriteData(fid,'object',self,'class','surfaceforcings','fieldname','isdelta18o','format','Boolean') - WriteData(fid,'object',self,'class','surfaceforcings','fieldname','desfac','format','Double'); - WriteData(fid,'object',self,'class','surfaceforcings','fieldname','s0p','format','Double'); - WriteData(fid,'object',self,'class','surfaceforcings','fieldname','s0t','format','Double'); - WriteData(fid,'object',self,'class','surfaceforcings','fieldname','rlaps','format','Double'); - WriteData(fid,'object',self,'class','surfaceforcings','fieldname','rlapslgm','format','Double'); - WriteData(fid,'object',self,'class','surfaceforcings','fieldname','Pfac','format','DoubleMat','mattype',1); - WriteData(fid,'object',self,'class','surfaceforcings','fieldname','Tdiff','format','DoubleMat','mattype',1); - WriteData(fid,'object',self,'class','surfaceforcings','fieldname','sealev','format','DoubleMat','mattype',1); - if self.isdelta18o: + WriteData(fid,'object',self,'class','surfaceforcings','fieldname','ismungsm','format','Boolean') + WriteData(fid,'object',self,'class','surfaceforcings','fieldname','desfac','format','Double') + WriteData(fid,'object',self,'class','surfaceforcings','fieldname','s0p','format','Double') + WriteData(fid,'object',self,'class','surfaceforcings','fieldname','s0t','format','Double') + WriteData(fid,'object',self,'class','surfaceforcings','fieldname','rlaps','format','Double') + WriteData(fid,'object',self,'class','surfaceforcings','fieldname','rlapslgm','format','Double') + + if not (self.isdelta18o and self.ismungsm): + WriteData(fid,'object',self,'class','surfaceforcings','fieldname','monthlytemperatures','format','DoubleMat','mattype',1) + WriteData(fid,'object',self,'class','surfaceforcings','fieldname','precipitation','format','DoubleMat','mattype',1,'scale',1./yts,'forcinglength',md.mesh.numberofvertices+1) + elif self.isdelta18o: WriteData(fid,'object',self,'class','surfaceforcings','fieldname','temperatures_presentday','format','DoubleMat','mattype',1) WriteData(fid,'object',self,'class','surfaceforcings','fieldname','temperatures_lgm','format','DoubleMat','mattype',1) WriteData(fid,'object',self,'class','surfaceforcings','fieldname','precipitations_presentday','format','DoubleMat','mattype',1) + WriteData(fid,'object',self,'class','surfaceforcings','fieldname','precipitations_lgm','format','DoubleMat','mattype',1) WriteData(fid,'object',self,'class','surfaceforcings','fieldname','delta18o_surface','format','DoubleMat','mattype',1) WriteData(fid,'object',self,'class','surfaceforcings','fieldname','delta18o','format','DoubleMat','mattype',1) - else: - WriteData(fid,'object',self,'class','surfaceforcings','fieldname','monthlytemperatures','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1) + WriteData(fid,'object',self,'class','surfaceforcings','fieldname','Tdiff','format','DoubleMat','mattype',1) + WriteData(fid,'object',self,'class','surfaceforcings','fieldname','sealev','format','DoubleMat','mattype',1) + elif self.ismungsm: + WriteData(fid,'object',self,'class','surfaceforcings','fieldname','temperatures_presentday','format','DoubleMat','mattype',1) + WriteData(fid,'object',self,'class','surfaceforcings','fieldname','temperatures_lgm','format','DoubleMat','mattype',1) + WriteData(fid,'object',self,'class','surfaceforcings','fieldname','precipitations_presentday','format','DoubleMat','mattype',1) + WriteData(fid,'object',self,'class','surfaceforcings','fieldname','precipitations_lgm','format','DoubleMat','mattype',1) + WriteData(fid,'object',self,'class','surfaceforcings','fieldname','Pfac','format','DoubleMat','mattype',1) + WriteData(fid,'object',self,'class','surfaceforcings','fieldname','Tdiff','format','DoubleMat','mattype',1) + WriteData(fid,'object',self,'class','surfaceforcings','fieldname','sealev','format','DoubleMat','mattype',1) # }}} Index: ../trunk-jpl/src/m/classes/SMBpdd.m =================================================================== --- ../trunk-jpl/src/m/classes/SMBpdd.m (revision 18967) +++ ../trunk-jpl/src/m/classes/SMBpdd.m (revision 18968) @@ -16,6 +16,7 @@ Tdiff = NaN; sealev = NaN; isdelta18o = 0; + ismungsm = 0; delta18o = NaN; delta18o_surface = NaN; temperatures_presentday = NaN; @@ -33,27 +34,30 @@ end end % }}} function self = extrude(self,md) % {{{ - - self.precipitation=project3d(md,'vector',self.precipitation,'type','node'); - self.monthlytemperatures=project3d(md,'vector',self.monthlytemperatures,'type','node'); + if(self.isdelta18o==0 & self.ismungsm==0),self.precipitation=project3d(md,'vector',self.precipitation,'type','node');end + if(self.isdelta18o==0 & self.ismungsm==0),self.monthlytemperatures=project3d(md,'vector',self.monthlytemperatures,'type','node');end if(self.isdelta18o),self.temperatures_lgm=project3d(md,'vector',self.temperatures_lgm,'type','node');end if(self.isdelta18o),self.temperatures_presentday=project3d(md,'vector',self.temperatures_presentday,'type','node');end if(self.isdelta18o),self.precipitations_presentday=project3d(md,'vector',self.precipitations_presentday,'type','node');end if(self.isdelta18o),self.precipitations_lgm=project3d(md,'vector',self.precipitations_lgm,'type','node');end + if(self.ismungsm),self.temperatures_lgm=project3d(md,'vector',self.temperatures_lgm,'type','node');end + if(self.ismungsm),self.temperatures_presentday=project3d(md,'vector',self.temperatures_presentday,'type','node');end + if(self.ismungsm),self.precipitations_presentday=project3d(md,'vector',self.precipitations_presentday,'type','node');end + if(self.ismungsm),self.precipitations_lgm=project3d(md,'vector',self.precipitations_lgm,'type','node');end - end % }}} function self = initialize(self,md) % {{{ + + % if isnan(self.precipitation), + % self.precipitation=zeros(md.mesh.numberofvertices,1); + % disp(' no SMBpdd.precipitation specified: values set as zero'); + % end - if isnan(self.precipitation), - self.precipitation=zeros(md.mesh.numberofvertices,1); - disp(' no SMBpdd.precipitation specified: values set as zero'); - end - end % }}} function obj = setdefaultparameters(obj) % {{{ obj.isdelta18o = 0; + obj.ismungsm = 0; obj.desfac = 0.5; obj.s0p = 0; obj.s0t = 0; @@ -69,44 +73,61 @@ md = checkfield(md,'fieldname','surfaceforcings.s0t','>=',0,'numel',1); md = checkfield(md,'fieldname','surfaceforcings.rlaps','>=',0,'numel',1); md = checkfield(md,'fieldname','surfaceforcings.rlapslgm','>=',0,'numel',1); - md = checkfield(md,'fieldname','surfaceforcings.Pfac','NaN',1,'size',[2,NaN]); - md = checkfield(md,'fieldname','surfaceforcings.Tdiff','NaN',1,'size',[2,NaN]); - md = checkfield(md,'fieldname','surfaceforcings.sealev','NaN',1,'size',[2,NaN]); - if(obj.isdelta18o==0) - md = checkfield(md,'fieldname','surfaceforcings.monthlytemperatures','forcing',1,'NaN',1); - md = checkfield(md,'fieldname','surfaceforcings.precipitation','forcing',1,'NaN',1); - else - md = checkfield(md,'fieldname','surfaceforcings.delta18o','NaN',1); - md = checkfield(md,'fieldname','surfaceforcings.delta18o_surface','NaN',1); - md = checkfield(md,'fieldname','surfaceforcings.temperatures_presentday','size',[md.mesh.numberofvertices+1 12],'NaN',1); - md = checkfield(md,'fieldname','surfaceforcings.temperatures_lgm','size',[md.mesh.numberofvertices+1 12],'NaN',1); - md = checkfield(md,'fieldname','surfaceforcings.precipitations_presentday','size',[md.mesh.numberofvertices+1 12],'NaN',1); - md = checkfield(md,'fieldname','surfaceforcings.precipitations_lgm','size',[md.mesh.numberofvertices+1 12],'NaN',1); + if(obj.isdelta18o==0 & obj.ismungsm==0) + md = checkfield(md,'fieldname','surfaceforcings.monthlytemperatures','forcing',1,'NaN',1); + md = checkfield(md,'fieldname','surfaceforcings.precipitation','forcing',1,'NaN',1); + elseif(obj.isdelta18o==1) + md = checkfield(md,'fieldname','surfaceforcings.delta18o','NaN',1); + md = checkfield(md,'fieldname','surfaceforcings.delta18o_surface','NaN',1); + md = checkfield(md,'fieldname','surfaceforcings.temperatures_presentday','size',[md.mesh.numberofvertices+1 12],'NaN',1); + md = checkfield(md,'fieldname','surfaceforcings.temperatures_lgm','size',[md.mesh.numberofvertices+1 12],'NaN',1); + md = checkfield(md,'fieldname','surfaceforcings.precipitations_presentday','size',[md.mesh.numberofvertices+1 12],'NaN',1); + md = checkfield(md,'fieldname','surfaceforcings.precipitations_lgm','size',[md.mesh.numberofvertices+1 12],'NaN',1); + md = checkfield(md,'fieldname','surfaceforcings.Tdiff','NaN',1); + md = checkfield(md,'fieldname','surfaceforcings.sealev','NaN',1); + elseif(obj.ismungsm==1) + md = checkfield(md,'fieldname','surfaceforcings.temperatures_presentday','size',[md.mesh.numberofvertices+1 12],'NaN',1); + md = checkfield(md,'fieldname','surfaceforcings.temperatures_lgm','size',[md.mesh.numberofvertices+1 12],'NaN',1); + md = checkfield(md,'fieldname','surfaceforcings.precipitations_presentday','size',[md.mesh.numberofvertices+1 12],'NaN',1); + md = checkfield(md,'fieldname','surfaceforcings.precipitations_lgm','size',[md.mesh.numberofvertices+1 12],'NaN',1); + md = checkfield(md,'fieldname','surfaceforcings.Pfac','NaN',1,'size',[2,NaN]); + md = checkfield(md,'fieldname','surfaceforcings.Tdiff','NaN',1); + md = checkfield(md,'fieldname','surfaceforcings.sealev','NaN',1); end - end + end end % }}} function disp(obj) % {{{ disp(sprintf(' surface forcings parameters:')); disp(sprintf('\n PDD and deltaO18 parameters:')); fielddisplay(obj,'isdelta18o','is temperature and precipitation delta18o parametrisation activated (0 or 1, default is 0)'); + fielddisplay(obj,'ismungsm','is temperature and precipitation mungsm parametrisation activated (0 or 1, default is 0)'); fielddisplay(obj,'desfac','desertification elevation factor (between 0 and 1, default is 0.5) [m]'); fielddisplay(obj,'s0p','should be set to elevation from precip source (between 0 and a few 1000s m, default is 0) [m]'); fielddisplay(obj,'s0t','should be set to elevation from temperature source (between 0 and a few 1000s m, default is 0) [m]'); fielddisplay(obj,'rlaps','present day lapse rate [degree/km]'); fielddisplay(obj,'rlapslgm','LGM lapse rate [degree/km]'); - fielddisplay(obj,'Pfac','time interpolation parameter for precipitation, 1D(year)'); - fielddisplay(obj,'Tdiff','time interpolation parameter for temperature, 1D(year)'); - fielddisplay(obj,'sealev','sea level [m], 1D(year)'); - fielddisplay(obj,'monthlytemperatures',['CURRENTLY NOT USED monthly surface temperatures [K], required if pdd is activated and delta18o not activated']); - fielddisplay(obj,'precipitation',['CURRENTLY NOT USED monthly surface precipitation [m/yr water eq]']); - fielddisplay(obj,'temperatures_presentday','monthly present day surface temperatures [K], required if pdd is activated and delta18o activated'); - fielddisplay(obj,'temperatures_lgm','monthly LGM surface temperatures [K], required if pdd is activated and delta18o activated'); - fielddisplay(obj,'precipitations_presentday','monthly surface precipitation [m/yr water eq], required if pdd is activated and delta18o activated'); - fielddisplay(obj,'precipitations_lgm','monthly surface precipitation [m/yr water eq], required if pdd is activated and delta18o activated'); - fielddisplay(obj,'delta18o','delta18o, required if pdd is activated and delta18o activated'); - fielddisplay(obj,'delta18o_surface','surface elevation of the delta18o site, required if pdd is activated and delta18o activated'); - + if(obj.isdelta18o==0 & obj.ismungsm==0) + fielddisplay(obj,'monthlytemperatures',['monthly surface temperatures [K], required if pdd is activated and delta18o not activated']); + fielddisplay(obj,'precipitation',['monthly surface precipitation [m/yr water eq], required if pdd is activated and delta18o or mungsm not activated']); + elseif(obj.isdelta18o==1) + fielddisplay(obj,'delta18o','delta18o, required if pdd is activated and delta18o activated'); + fielddisplay(obj,'delta18o_surface','surface elevation of the delta18o site, required if pdd is activated and delta18o activated'); + fielddisplay(obj,'temperatures_presentday','monthly present day surface temperatures [K], required if delta18o/mungsm is activated'); + fielddisplay(obj,'temperatures_lgm','monthly LGM surface temperatures [K], required if delta18o or mungsm is activated'); + fielddisplay(obj,'precipitations_presentday','monthly surface precipitation [m/yr water eq], required if delta18o or mungsm is activated'); + fielddisplay(obj,'precipitations_lgm','monthly surface precipitation [m/yr water eq], required if delta18o or mungsm is activated'); + fielddisplay(obj,'Tdiff','time interpolation parameter for temperature, 1D(year), required if mungsm is activated'); + fielddisplay(obj,'sealev','sea level [m], 1D(year), required if mungsm is activated'); + elseif(obj.ismungsm==1) + fielddisplay(obj,'temperatures_presentday','monthly present day surface temperatures [K], required if delta18o/mungsm is activated'); + fielddisplay(obj,'temperatures_lgm','monthly LGM surface temperatures [K], required if delta18o or mungsm is activated'); + fielddisplay(obj,'precipitations_presentday','monthly surface precipitation [m/yr water eq], required if delta18o or mungsm is activated'); + fielddisplay(obj,'precipitations_lgm','monthly surface precipitation [m/yr water eq], required if delta18o or mungsm is activated'); + fielddisplay(obj,'Pfac','time interpolation parameter for precipitation, 1D(year), required if mungsm is activated'); + fielddisplay(obj,'Tdiff','time interpolation parameter for temperature, 1D(year), required if mungsm is activated'); + fielddisplay(obj,'sealev','sea level [m], 1D(year), required if mungsm is activated'); + end end % }}} function marshall(obj,md,fid) % {{{ @@ -114,25 +135,35 @@ WriteData(fid,'enum',SurfaceforcingsEnum(),'data',SMBpddEnum(),'format','Integer'); - WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','precipitation','format','DoubleMat','mattype',1,'scale',1./yts,'forcinglength',md.mesh.numberofvertices+1); + WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','isdelta18o','format','Boolean'); + WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','ismungsm','format','Boolean'); WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','desfac','format','Double'); WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','s0p','format','Double'); WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','s0t','format','Double'); WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','rlaps','format','Double'); WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','rlapslgm','format','Double'); - WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','Pfac','format','DoubleMat','mattype',1); - WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','Tdiff','format','DoubleMat','mattype',1); - WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','sealev','format','DoubleMat','mattype',1); - WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','isdelta18o','format','Boolean'); - if obj.isdelta18o + + if(obj.isdelta18o==0 & obj.ismungsm==0) + %WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','monthlytemperatures','format','DoubleMat','mattype',1); + WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','monthlytemperatures','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1); + WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','precipitation','format','DoubleMat','mattype',1,'scale',1./yts,'forcinglength',md.mesh.numberofvertices+1); + elseif obj.isdelta18o WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','temperatures_presentday','format','DoubleMat','mattype',1); WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','temperatures_lgm','format','DoubleMat','mattype',1); WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','precipitations_presentday','format','DoubleMat','mattype',1); WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','precipitations_lgm','format','DoubleMat','mattype',1); WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','delta18o_surface','format','DoubleMat','mattype',1); WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','delta18o','format','DoubleMat','mattype',1); - else - WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','monthlytemperatures','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1); + WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','Tdiff','format','DoubleMat','mattype',1); + WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','sealev','format','DoubleMat','mattype',1); + elseif obj.ismungsm + WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','temperatures_presentday','format','DoubleMat','mattype',1); + WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','temperatures_lgm','format','DoubleMat','mattype',1); + WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','precipitations_presentday','format','DoubleMat','mattype',1); + WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','precipitations_lgm','format','DoubleMat','mattype',1); + WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','Pfac','format','DoubleMat','mattype',1); + WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','Tdiff','format','DoubleMat','mattype',1); + WriteData(fid,'object',obj,'class','surfaceforcings','fieldname','sealev','format','DoubleMat','mattype',1); end end % }}} end Index: ../trunk-jpl/src/m/enum/EnumDefinitions.py =================================================================== --- ../trunk-jpl/src/m/enum/EnumDefinitions.py (revision 18967) +++ ../trunk-jpl/src/m/enum/EnumDefinitions.py (revision 18968) @@ -341,6 +341,7 @@ def SurfaceforcingsDelta18oEnum(): return StringToEnum("SurfaceforcingsDelta18o")[0] def SurfaceforcingsDelta18oSurfaceEnum(): return StringToEnum("SurfaceforcingsDelta18oSurface")[0] def SurfaceforcingsIsdelta18oEnum(): return StringToEnum("SurfaceforcingsIsdelta18o")[0] +def SurfaceforcingsIsmungsmEnum(): return StringToEnum("SurfaceforcingsIsmungsm")[0] def SurfaceforcingsPrecipitationsPresentdayEnum(): return StringToEnum("SurfaceforcingsPrecipitationsPresentday")[0] def SurfaceforcingsPrecipitationsLgmEnum(): return StringToEnum("SurfaceforcingsPrecipitationsLgm")[0] def SurfaceforcingsTemperaturesPresentdayEnum(): return StringToEnum("SurfaceforcingsTemperaturesPresentday")[0]