[issm-svn] r12293 - issm/trunk/src
lemorzad at issm.ess.uci.edu
lemorzad at issm.ess.uci.edu
Tue May 29 09:20:06 PDT 2012
Author: lemorzad
Date: 2012-05-29 09:20:06 -0700 (Tue, 29 May 2012)
New Revision: 12293
Added:
issm/trunk/src/m/classes/clusters/acenet.m
issm/trunk/src/m/enum/EnumToModelField.m
Modified:
issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h
issm/trunk/src/c/modules/EnumToStringx/EnumToStringx.cpp
issm/trunk/src/c/modules/ModelProcessorx/CreateParameters.cpp
issm/trunk/src/c/modules/ModelProcessorx/Prognostic/UpdateElementsPrognostic.cpp
issm/trunk/src/c/modules/PositiveDegreeDayx/PositiveDegreeDayx.cpp
issm/trunk/src/c/modules/StringToEnumx/StringToEnumx.cpp
issm/trunk/src/c/modules/modules.h
issm/trunk/src/c/objects/Elements/Penta.cpp
issm/trunk/src/c/objects/Elements/Tria.cpp
issm/trunk/src/c/solutions/prognostic_core.cpp
issm/trunk/src/m/classes/surfaceforcings.m
issm/trunk/src/m/utils/BC/SetMarineIceSheetBC.m
Log:
Modified: issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h
===================================================================
--- issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h 2012-05-26 00:14:19 UTC (rev 12292)
+++ issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h 2012-05-29 16:20:06 UTC (rev 12293)
@@ -483,6 +483,7 @@
/*Rheology law (move too Material) {{{1*/
PatersonEnum,
ArrheniusEnum,
+ SurfaceforcingsIspddEnum,
/*}}}*/
MaximumNumberOfEnums
};
Modified: issm/trunk/src/c/modules/EnumToStringx/EnumToStringx.cpp
===================================================================
--- issm/trunk/src/c/modules/EnumToStringx/EnumToStringx.cpp 2012-05-26 00:14:19 UTC (rev 12292)
+++ issm/trunk/src/c/modules/EnumToStringx/EnumToStringx.cpp 2012-05-29 16:20:06 UTC (rev 12293)
@@ -451,6 +451,7 @@
case OptionLogicalEnum : return "OptionLogical";
case PatersonEnum : return "Paterson";
case ArrheniusEnum : return "Arrhenius";
+ case SurfaceforcingsIspddEnum : return "SurfaceforcingsIspdd";
default : return "unknown";
}
Modified: issm/trunk/src/c/modules/ModelProcessorx/CreateParameters.cpp
===================================================================
--- issm/trunk/src/c/modules/ModelProcessorx/CreateParameters.cpp 2012-05-26 00:14:19 UTC (rev 12292)
+++ issm/trunk/src/c/modules/ModelProcessorx/CreateParameters.cpp 2012-05-29 16:20:06 UTC (rev 12293)
@@ -88,6 +88,7 @@
parameters->AddObject(iomodel->CopyConstantObject(QmuIsdakotaEnum));
parameters->AddObject(iomodel->CopyConstantObject(InversionIscontrolEnum));
parameters->AddObject(iomodel->CopyConstantObject(InversionTaoEnum));
+ parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsIspddEnum));
/*some parameters that did not come with the iomodel: */
parameters->AddObject(new IntParam(SolutionTypeEnum,solution_type));
Modified: issm/trunk/src/c/modules/ModelProcessorx/Prognostic/UpdateElementsPrognostic.cpp
===================================================================
--- issm/trunk/src/c/modules/ModelProcessorx/Prognostic/UpdateElementsPrognostic.cpp 2012-05-26 00:14:19 UTC (rev 12292)
+++ issm/trunk/src/c/modules/ModelProcessorx/Prognostic/UpdateElementsPrognostic.cpp 2012-05-29 16:20:06 UTC (rev 12293)
@@ -19,6 +19,7 @@
int numberofelements;
int stabilization;
bool dakota_analysis;
+ bool ispdd;
/*Fetch data needed: */
iomodel->Constant(&dim,MeshDimensionEnum);
@@ -26,6 +27,7 @@
iomodel->Constant(&stabilization,PrognosticStabilizationEnum);
iomodel->Constant(&dakota_analysis,QmuIsdakotaEnum);
iomodel->FetchData(1,MeshElementsEnum);
+ iomodel->Constant(&ispdd,SurfaceforcingsIspddEnum);
/*Update elements: */
int counter=0;
@@ -43,6 +45,7 @@
iomodel->FetchDataToInput(elements,BathymetryEnum);
iomodel->FetchDataToInput(elements,MaskElementonfloatingiceEnum);
iomodel->FetchDataToInput(elements,MaskElementonwaterEnum);
+ iomodel->FetchDataToInput(elements,SurfaceforcingsPrecipitationEnum);
iomodel->FetchDataToInput(elements,SurfaceforcingsMassBalanceEnum);
iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum);
iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateCorrectionEnum);
@@ -65,7 +68,11 @@
iomodel->FetchDataToInput(elements,PressureEnum);
iomodel->FetchDataToInput(elements,TemperatureEnum);
}
-
+ if(ispdd){
+ iomodel->FetchDataToInput(elements,VyEnum);
+ iomodel->FetchDataToInput(elements,ThermalSpctemperatureEnum);
+ iomodel->FetchDataToInput(elements,SurfaceforcingsPrecipitationEnum);
+ }
/*Free data: */
iomodel->DeleteData(1,MeshElementsEnum);
}
Modified: issm/trunk/src/c/modules/PositiveDegreeDayx/PositiveDegreeDayx.cpp
===================================================================
--- issm/trunk/src/c/modules/PositiveDegreeDayx/PositiveDegreeDayx.cpp 2012-05-26 00:14:19 UTC (rev 12292)
+++ issm/trunk/src/c/modules/PositiveDegreeDayx/PositiveDegreeDayx.cpp 2012-05-29 16:20:06 UTC (rev 12293)
@@ -19,14 +19,14 @@
// ,NTIME)
// monthly mean precip rate (m/yr water equivalent): vPrec(NA,NTIME)
// OUTPUT: mass-balance (m/yr ice): agd(NA)
-// mean annual surface temperature (degrees C): Tsurf(NA)
+// mean annual sssurface temperature (degrees C): Tsurf(NA)
int i, it, jj, itm;
double DT = 0.02, sigfac, snormfac;
double signorm = 5.5; // signorm : sigma of the temperature distribution for a normal day
double siglim; // sigma limit for the integration which is equal to 2.5 sigmanorm
- double signormc; // sigma of the temperature distribution for cloudy day
- double siglimc=0, siglim0, siglim0c;
+ double signormc = signorm - 0.5; // sigma of the temperature distribution for cloudy day
+ double siglimc, siglim0, siglim0c;
double tstep, tsint, tint, tstepc;
int NPDMAX = 1504, NPDCMAX = 1454;
//double pdds[NPDMAX]={0};
@@ -42,16 +42,17 @@
pdds=(double*)xmalloc(NPDMAX*sizeof(double)+1);
pds=(double*)xmalloc(NPDCMAX*sizeof(double)+1);
-
- // PDD constant
- siglim = 2.5*signorm;
-
// initialize PDD (creation of a lookup table)
tstep = 0.1;
tsint = tstep*0.5;
sigfac = -1.0/(2.0*pow(signorm,2));
snormfac = 1.0/(signorm*sqrt(2.0*acos(-1.0)));
siglim = 2.5*signorm;
+ siglimc = 2.5*signormc
+ siglim0 = siglim/DT + 0.5
+ siglim0c = siglimc/DT + 0.5
+ PDup = siglimc+PDCUT
+
itm = (int)(2*siglim/DT + 1.5);
if (itm >= NPDMAX){
Modified: issm/trunk/src/c/modules/StringToEnumx/StringToEnumx.cpp
===================================================================
--- issm/trunk/src/c/modules/StringToEnumx/StringToEnumx.cpp 2012-05-26 00:14:19 UTC (rev 12292)
+++ issm/trunk/src/c/modules/StringToEnumx/StringToEnumx.cpp 2012-05-29 16:20:06 UTC (rev 12293)
@@ -461,6 +461,7 @@
else if (strcmp(name,"OptionLogical")==0) return OptionLogicalEnum;
else if (strcmp(name,"Paterson")==0) return PatersonEnum;
else if (strcmp(name,"Arrhenius")==0) return ArrheniusEnum;
+ else if (strcmp(name,"SurfaceforcingsIspdd")==0) return SurfaceforcingsIspddEnum;
else stage=5;
}
/*If we reach this point, the string provided has not been found*/
Modified: issm/trunk/src/c/modules/modules.h
===================================================================
--- issm/trunk/src/c/modules/modules.h 2012-05-26 00:14:19 UTC (rev 12292)
+++ issm/trunk/src/c/modules/modules.h 2012-05-29 16:20:06 UTC (rev 12293)
@@ -92,6 +92,7 @@
#include "./OutputResultsx/OutputResultsx.h"
#include "./ConstraintsStatex/ConstraintsStatex.h"
#include "./PointCloudFindNeighborsx/PointCloudFindNeighborsx.h"
+#include "./PositiveDegreeDayx/PositiveDegreeDayx.h"
#include "./PropagateFlagsFromConnectivityx/PropagateFlagsFromConnectivityx.h"
#include "./Dakotax/Dakotax.h"
#include "./Reduceloadx/Reduceloadx.h"
@@ -121,4 +122,5 @@
#include "./VerticesDofx/VerticesDofx.h"
#include "./VecMergex/VecMergex.h"
+
#endif
Modified: issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- issm/trunk/src/c/objects/Elements/Penta.cpp 2012-05-26 00:14:19 UTC (rev 12292)
+++ issm/trunk/src/c/objects/Elements/Penta.cpp 2012-05-29 16:20:06 UTC (rev 12293)
@@ -2433,7 +2433,8 @@
// PDD and PD constants and variables
double siglim; // sigma limit for the integration which is equal to 2.5 sigmanorm
- double siglimc=0, siglim0, siglim0c;
+ double signormc = signorm - 0.5; // sigma of the temperature distribution for cloudy day
+ double siglimc, siglim0, siglim0c;
double PDup, pddsig, PDCUT = 2.0; // PDcut: rain/snow cutoff temperature (C)
double DT = 0.02;
double pddt, pd; // pd: snow/precip fraction, precipitation falling as snow
@@ -2479,13 +2480,13 @@
/*Get material parameters :*/
rho_ice=matpar->GetRhoIce();
- rho_water=matpar->GetRhoWater();
- density=rho_ice/rho_water;
+ //rho_freshwater=matpar->GetRhoWater();
- sconv=(rho_water/rho_ice)/12.; //rhow_rain/rhoi / 12 months
+ sconv=(1000/rho_ice)/12.; //rhow_rain/rhoi / 12 months
/*PDD constant*/
siglim = 2.5*signorm;
+ siglimc = 2.5*signormc
siglim0 = siglim/DT + 0.5;
siglim0c = siglimc/DT + 0.5;
PDup = siglimc+PDCUT;
Modified: issm/trunk/src/c/objects/Elements/Tria.cpp
===================================================================
--- issm/trunk/src/c/objects/Elements/Tria.cpp 2012-05-26 00:14:19 UTC (rev 12292)
+++ issm/trunk/src/c/objects/Elements/Tria.cpp 2012-05-29 16:20:06 UTC (rev 12293)
@@ -2237,7 +2237,7 @@
void Tria::PositiveDegreeDay(double* pdds,double* pds,double signorm){
int i,iqj,imonth;
- double agd[NUMVERTICES]; // surface and basal
+ double agd[NUMVERTICES]; // surface mass balance
double saccu[NUMVERTICES] = {0}; // yearly surface accumulation
double smelt[NUMVERTICES] = {0}; // yearly melt
double precrunoff[NUMVERTICES]; // yearly runoff
@@ -2246,18 +2246,19 @@
double runoff; //meltwater only, does not include rain
double sconv; //rhow_rain/rhoi / 12 months
- double rho_water,rho_ice,density;
- double lapser=6.5/1000, sealev=0; // lapse rate. degrees per meter.
- double desfac = 0.5; //desert elevation factor
- double s0p[NUMVERTICES]={0}; //should be set to elevation from precip source
- double s0t[NUMVERTICES]={0}; //should be set to elevation from temperature source
+ double rho_water,rho_ice,density;
+ double lapser=6.5/1000, sealev=0; // lapse rate. degrees per meter. 7.5 lev's 99 paper, 9 Marshall 99 paper
+ double desfac = 0.5; // desert elevation factor
+ double s0p[NUMVERTICES]={0}; // should be set to elevation from precip source
+ double s0t[NUMVERTICES]={0}; // should be set to elevation from temperature source
double st; // elevation between altitude of the temp record and current altitude
double sp; // elevation between altitude of the prec record and current altitude
// PDD and PD constants and variables
double siglim; // sigma limit for the integration which is equal to 2.5 sigmanorm
- double siglimc=0, siglim0, siglim0c;
+ double signormc = signorm - 0.5; // sigma of the temperature distribution for cloudy day
+ double siglimc, siglim0, siglim0c;
double PDup, pddsig, PDCUT = 2.0; // PDcut: rain/snow cutoff temperature (C)
double DT = 0.02;
double pddt, pd; // pd: snow/precip fraction, precipitation falling as snow
@@ -2274,7 +2275,7 @@
double Tsurf[NUMVERTICES] = {0}; // average annual temperature
double h[NUMVERTICES],s[NUMVERTICES],ttmp[NUMVERTICES],prectmp[NUMVERTICES]; // ,b[NUMVERTICES]
- double t[NUMVERTICES][12],prec[NUMVERTICES][12];
+ double t[NUMVERTICES+1][12],prec[NUMVERTICES+1][12];
double deltm=1/12;
int ismon[12]={12,1,2,3,4,5,6,7,8,9,10,11};
@@ -2303,13 +2304,13 @@
/*Get material parameters :*/
rho_ice=matpar->GetRhoIce();
- rho_water=matpar->GetRhoWater();
- density=rho_ice/rho_water;
+ //rho_freshwater=matpar->GetRhoWater();
- sconv=(rho_water/rho_ice)/12.; //rhow_rain/rhoi / 12 months
+ sconv=(1000/rho_ice)/12.; //rhow_rain/rhoi / 12 months
/*PDD constant*/
siglim = 2.5*signorm;
+ siglimc = 2.5*signormc
siglim0 = siglim/DT + 0.5;
siglim0c = siglimc/DT + 0.5;
PDup = siglimc+PDCUT;
@@ -2334,7 +2335,7 @@
if (sp>0.0){q = exp(-desfac*sp);}
else {q = 1.0;}
- qmt[i]= qmt[i] + prec[i][imonth]*sconv; //*sconv to convert in m of ice equivalent
+ qmt[i]= qmt[i] + prec[i][imonth]*sconv; //*sconv to convert in m of ice equivalent per month
qmpt= q*prec[i][imonth]*sconv;
qmp[i]= qmp[i] + qmpt;
qm[i]= qm[i] + qmpt*pd;
Modified: issm/trunk/src/c/solutions/prognostic_core.cpp
===================================================================
--- issm/trunk/src/c/solutions/prognostic_core.cpp 2012-05-26 00:14:19 UTC (rev 12292)
+++ issm/trunk/src/c/solutions/prognostic_core.cpp 2012-05-29 16:20:06 UTC (rev 12293)
@@ -15,16 +15,23 @@
/*parameters: */
bool save_results;
+ bool ispdd;
/*activate formulation: */
femmodel->SetCurrentConfiguration(PrognosticAnalysisEnum);
/*recover parameters: */
femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
+ femmodel->parameters->FindParam(&ispdd,SurfaceforcingsIspddEnum);
+ if(ispdd){
+ _printf_(VerboseSolution()," call positive degree day module\n");
+ PositiveDegreeDayx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+ }
+
_printf_(VerboseSolution()," call computational core\n");
solver_linear(femmodel);
-
+
if(save_results){
_printf_(VerboseSolution()," saving results\n");
InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum);
Modified: issm/trunk/src/m/classes/surfaceforcings.m
===================================================================
--- issm/trunk/src/m/classes/surfaceforcings.m 2012-05-26 00:14:19 UTC (rev 12292)
+++ issm/trunk/src/m/classes/surfaceforcings.m 2012-05-29 16:20:06 UTC (rev 12293)
@@ -7,6 +7,7 @@
properties (SetAccess=public)
precipitation = NaN;
mass_balance = NaN;
+ ispdd = NaN;
end
methods
function obj = surfaceforcings(varargin) % {{{
@@ -18,6 +19,9 @@
end
end % }}}
function obj = setdefaultparameters(obj) % {{{
+
+ %pdd methode not use in default mode
+ obj.ispdd=0;
end % }}}
function checkconsistency(obj,md,solution,analyses) % {{{
@@ -28,17 +32,23 @@
if ismember(BalancethicknessAnalysisEnum,analyses),
checkfield(md,'surfaceforcings.mass_balance','size',[md.mesh.numberofvertices 1],'NaN',1);
end
+ if (ismember(PrognosticAnalysisEnum,analyses) & md.surfaceforcings.ispdd),
+ checkfield(md,'surfaceforcings.ispdd','numel',1,'values',[0 1]);
+ end
end % }}}
function disp(obj) % {{{
disp(sprintf(' surface forcings parameters:'));
fielddisplay(obj,'precipitation','surface precipitation [m/yr water eq]');
fielddisplay(obj,'mass_balance','surface mass balance [m/yr ice eq]');
+ fielddisplay(obj,'ispdd','is pdd activated (0 or 1, default is 0)');
end % }}}
function marshall(obj,fid) % {{{
WriteData(fid,'object',obj,'fieldname','precipitation','format','DoubleMat','mattype',1);
WriteData(fid,'object',obj,'fieldname','mass_balance','format','DoubleMat','mattype',1);
+ WriteData(fid,'object',obj,'fieldname','ispdd','format','Boolean');
+
end % }}}
end
end
Modified: issm/trunk/src/m/utils/BC/SetMarineIceSheetBC.m
===================================================================
--- issm/trunk/src/m/utils/BC/SetMarineIceSheetBC.m 2012-05-26 00:14:19 UTC (rev 12292)
+++ issm/trunk/src/m/utils/BC/SetMarineIceSheetBC.m 2012-05-29 16:20:06 UTC (rev 12293)
@@ -57,6 +57,7 @@
%segment on Neumann (Ice Front)
pos=find(vertexonicefront(md.mesh.segments(:,1)) | vertexonicefront(md.mesh.segments(:,2)));
+
if (md.mesh.dimension==2)
pressureload=md.mesh.segments(pos,:);
elseif md.mesh.dimension==3
@@ -105,3 +106,4 @@
else
disp(' no thermal boundary conditions created: no observed temperature found');
end
+
More information about the issm-svn
mailing list