Index: /issm/trunk-jpl/src/c/shared/Elements/ComputeD18OTemperaturePrecipitationFromPD.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Elements/ComputeD18OTemperaturePrecipitationFromPD.cpp (revision 19173) +++ /issm/trunk-jpl/src/c/shared/Elements/ComputeD18OTemperaturePrecipitationFromPD.cpp (revision 19173) @@ -0,0 +1,31 @@ +/* file: ComputeTemperaturePrecipitationfrom018.cpp + Scale present day monthly precipitation and temperature fields + along the NGRIP oxygen isotope record. + */ + +#include "./elements.h" +#include "../Numerics/numerics.h" + +void ComputeD18OTemperaturePrecipitationFromPD(IssmDouble d018,IssmDouble dpermil, + IssmDouble* PrecipitationPresentday,IssmDouble* TemperaturePresentday, + IssmDouble* monthlytemperaturesout, IssmDouble* monthlyprecout){ + + IssmDouble monthlytemperaturestmp[12],monthlyprectmp[12]; + IssmDouble deltaTemp; + + /* Constants */ + // dpermil = 2.4;/*degrees C per mil*/ + + /*Create Delta Temp to be applied to monthly temps and used in precip scaling*/ + deltaTemp = dpermil * (d018+34.83); + + for (int imonth = 0; imonth<12; imonth++){ + + monthlytemperaturesout[imonth] = TemperaturePresentday[imonth] + deltaTemp; + monthlyprecout[imonth] = PrecipitationPresentday[imonth]*exp((0.169/dpermil)*deltaTemp); + + /*Assign output pointer*/ + *(monthlytemperaturesout+imonth) = monthlytemperaturestmp[imonth]; + *(monthlyprecout+imonth) = monthlyprectmp[imonth]; + } +} Index: /issm/trunk-jpl/src/m/enum/SurfaceforcingsDpermilEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/SurfaceforcingsDpermilEnum.m (revision 19173) +++ /issm/trunk-jpl/src/m/enum/SurfaceforcingsDpermilEnum.m (revision 19173) @@ -0,0 +1,11 @@ +function macro=SurfaceforcingsDpermilEnum() +%SURFACEFORCINGSDPERMILENUM - Enum of SurfaceforcingsDpermil +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=SurfaceforcingsDpermilEnum() + +macro=StringToEnum('SurfaceforcingsDpermil'); Index: /issm/trunk-jpl/src/m/enum/SurfaceforcingsIsd18opdEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/SurfaceforcingsIsd18opdEnum.m (revision 19173) +++ /issm/trunk-jpl/src/m/enum/SurfaceforcingsIsd18opdEnum.m (revision 19173) @@ -0,0 +1,11 @@ +function macro=SurfaceforcingsIsd18opdEnum() +%SURFACEFORCINGSISD18OPDENUM - Enum of SurfaceforcingsIsd18opd +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=SurfaceforcingsIsd18opdEnum() + +macro=StringToEnum('SurfaceforcingsIsd18opd'); Index: /issm/trunk-jpl/test/NightlyRun/test2361.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2361.m (revision 19173) +++ /issm/trunk-jpl/test/NightlyRun/test2361.m (revision 19173) @@ -0,0 +1,88 @@ +%Test Name: SquareShelfTranIspddIsdeltaSSA2d +md=triangle(model(),'../Exp/Square.exp',150000.); +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.isd18opd=1; +%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 +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 temperatures +% Same temperature over the all region: +tmonth(1:12)=238.15+20.; +for imonth=0:11 + md.surfaceforcings.temperatures_presentday(1:md.mesh.numberofvertices,imonth+1)=tmonth(imonth+1); + % Time for the last line: + md.surfaceforcings.temperatures_presentday(md.mesh.numberofvertices+1,imonth+1)=((imonth+1)/12); +end + +% creating initialization and spc temperatures initialization and +% spc +md.thermal.spctemperature=mean(md.surfaceforcings.temperatures_presentday(1:md.mesh.numberofvertices,1:12),2)-10; %-10*ones(md.mesh.numberofvertices,1); + +md.initialization.temperature=md.thermal.spctemperature; %md.surfaceforcings.temperatures_presentday(1:md.mesh.numberofvertices,1); + +% 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; + % Time for the last line: + md.surfaceforcings.precipitations_presentday(md.mesh.numberofvertices+1,imonth+1)=((imonth+1)/12); +end + +% time steps and resolution +md.timestepping.time_step=20; +md.settings.output_frequency=1; +md.timestepping.final_time=60; + +% +md.transient.requested_outputs={'default','SurfaceforcingsMonthlytemperatures'}; +md=setflowequation(md,'SSA','all'); +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','Vel1','Pressure1','Bed1','Surface1','Thickness1','Vx2','Vy2','Vel2','Pressure2','Bed2','Surface2','Thickness2','Vx3','Vy3','Vel3','Pressure3','Bed3','Surface3','Thickness3'}; +field_tolerances={1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,... + 1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,... + 1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13}; +field_values={... + (md.results.TransientSolution(1).Vx),... + (md.results.TransientSolution(1).Vy),... + (md.results.TransientSolution(1).Vel),... + (md.results.TransientSolution(1).Pressure),... + (md.results.TransientSolution(1).Base),... + (md.results.TransientSolution(1).Surface),... + (md.results.TransientSolution(1).Thickness),... + (md.results.TransientSolution(1).SurfaceforcingsMonthlytemperatures),... + (md.results.TransientSolution(1).SurfaceforcingsMassBalance),... + (md.results.TransientSolution(2).Vx),... + (md.results.TransientSolution(2).Vy),... + (md.results.TransientSolution(2).Vel),... + (md.results.TransientSolution(2).Pressure),... + (md.results.TransientSolution(2).Base),... + (md.results.TransientSolution(2).Surface),... + (md.results.TransientSolution(2).Thickness),... + (md.results.TransientSolution(2).SurfaceforcingsMonthlytemperatures),... + (md.results.TransientSolution(2).SurfaceforcingsMassBalance),... + (md.results.TransientSolution(3).Vx),... + (md.results.TransientSolution(3).Vy),... + (md.results.TransientSolution(3).Vel),... + (md.results.TransientSolution(3).Pressure),... + (md.results.TransientSolution(3).Base),... + (md.results.TransientSolution(3).Surface),... + (md.results.TransientSolution(3).Thickness),... + (md.results.TransientSolution(3).SurfaceforcingsMonthlytemperatures),... + (md.results.TransientSolution(3).SurfaceforcingsMassBalance),... + };