Changeset 12293

Timestamp:

05/29/12 09:20:06 (13 years ago)

Author:

lemorzad

Message:

 

Location:

issm/trunk/src

Files:

• c/EnumDefinitions/EnumDefinitions.h (modified) (1 diff)
• c/modules/EnumToStringx/EnumToStringx.cpp (modified) (1 diff)
• c/modules/ModelProcessorx/CreateParameters.cpp (modified) (1 diff)
• c/modules/ModelProcessorx/Prognostic/UpdateElementsPrognostic.cpp (modified) (4 diffs)
• c/modules/PositiveDegreeDayx/PositiveDegreeDayx.cpp (modified) (4 diffs)
• c/modules/StringToEnumx/StringToEnumx.cpp (modified) (1 diff)
• c/modules/modules.h (modified) (2 diffs)
• c/objects/Elements/Penta.cpp (modified) (2 diffs)
• c/objects/Elements/Tria.cpp (modified) (6 diffs)
• c/solutions/prognostic_core.cpp (modified) (2 diffs)
• m/classes/clusters/acenet.m (added)
• m/classes/surfaceforcings.m (modified) (5 diffs)
• m/enum/EnumToModelField.m (copied) (copied from issm/trunk/src/m/enum/EnumToModelField.m ) (1 diff)
• m/utils/BC/SetMarineIceSheetBC.m (modified) (2 diffs)

issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h ¶

@@ -484 +484 @@
         PatersonEnum,
         ArrheniusEnum,
+        SurfaceforcingsIspddEnum,
         /*}}}*/
         MaximumNumberOfEnums

issm/trunk/src/c/modules/EnumToStringx/EnumToStringx.cpp ¶

@@ -452 +452 @@
                 case PatersonEnum : return "Paterson";
                 case ArrheniusEnum : return "Arrhenius";
+                case SurfaceforcingsIspddEnum : return "SurfaceforcingsIspdd";
                 default : return "unknown";
 

issm/trunk/src/c/modules/ModelProcessorx/CreateParameters.cpp ¶

@@ -89 +89 @@
         parameters->AddObject(iomodel->CopyConstantObject(InversionIscontrolEnum));
         parameters->AddObject(iomodel->CopyConstantObject(InversionTaoEnum));
+        parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsIspddEnum));
 
         /*some parameters that did not come with the iomodel: */

issm/trunk/src/c/modules/ModelProcessorx/Prognostic/UpdateElementsPrognostic.cpp ¶

@@ -20 +20 @@
         int    stabilization;
         bool   dakota_analysis;
+        bool   ispdd;
 
         /*Fetch data needed: */
@@ -27 +28 @@
         iomodel->Constant(&dakota_analysis,QmuIsdakotaEnum);
         iomodel->FetchData(1,MeshElementsEnum);
+        iomodel->Constant(&ispdd,SurfaceforcingsIspddEnum);
 
         /*Update elements: */
@@ -44 +46 @@
         iomodel->FetchDataToInput(elements,MaskElementonfloatingiceEnum);
         iomodel->FetchDataToInput(elements,MaskElementonwaterEnum);
+        iomodel->FetchDataToInput(elements,SurfaceforcingsPrecipitationEnum);
         iomodel->FetchDataToInput(elements,SurfaceforcingsMassBalanceEnum);
         iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum);
@@ -66 +69 @@
                 iomodel->FetchDataToInput(elements,TemperatureEnum);
         }
-
+        if(ispdd){
+          iomodel->FetchDataToInput(elements,VyEnum); 
+          iomodel->FetchDataToInput(elements,ThermalSpctemperatureEnum);
+          iomodel->FetchDataToInput(elements,SurfaceforcingsPrecipitationEnum);
+        }
         /*Free data: */
         iomodel->DeleteData(1,MeshElementsEnum);

issm/trunk/src/c/modules/PositiveDegreeDayx/PositiveDegreeDayx.cpp ¶

@@ -20 +20 @@
 //    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;
@@ -26 +26 @@
   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;
@@ -43 +43 @@
   pds=(double*)xmalloc(NPDCMAX*sizeof(double)+1); 
   
-  
-  // PDD constant
-  siglim = 2.5*signorm; 
-  
   // initialize PDD (creation of a lookup table)
   tstep = 0.1;
@@ -53 +49 @@
   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);
   

issm/trunk/src/c/modules/StringToEnumx/StringToEnumx.cpp ¶

@@ -462 +462 @@
               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;
    }

issm/trunk/src/c/modules/modules.h ¶

@@ -93 +93 @@
 #include "./ConstraintsStatex/ConstraintsStatex.h"
 #include "./PointCloudFindNeighborsx/PointCloudFindNeighborsx.h"
+#include "./PositiveDegreeDayx/PositiveDegreeDayx.h"
 #include "./PropagateFlagsFromConnectivityx/PropagateFlagsFromConnectivityx.h"
 #include "./Dakotax/Dakotax.h"
@@ -122 +123 @@
 #include "./VecMergex/VecMergex.h"
 
+
 #endif

issm/trunk/src/c/objects/Elements/Penta.cpp ¶

@@ -2434 +2434 @@
    // 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;
@@ -2480 +2481 @@
    /*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;

issm/trunk/src/c/objects/Elements/Tria.cpp ¶

@@ -2238 +2238 @@
 
    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
@@ -2247 +2247 @@
    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
@@ -2258 +2258 @@
    // 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;
@@ -2275 +2276 @@
    
    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};
@@ -2304 +2305 @@
    /*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;
@@ -2335 +2336 @@
        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;

issm/trunk/src/c/solutions/prognostic_core.cpp ¶

@@ -16 +16 @@
         /*parameters: */
         bool save_results;
+        bool ispdd;
 
         /*activate formulation: */
@@ -22 +23 @@
         /*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");

issm/trunk/src/m/classes/surfaceforcings.m ¶

@@ -8 +8 @@
                 precipitation = NaN;
                 mass_balance  = NaN;
+                ispdd = NaN;
         end
         methods
@@ -19 +20 @@
                 end % }}}
                 function obj = setdefaultparameters(obj) % {{{
+                  
+                  %pdd methode not use in default mode
+                  obj.ispdd=0;
 
                 end % }}}
@@ -29 +33 @@
                                 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) % {{{
@@ -35 +42 @@
                         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 % }}}
@@ -40 +48 @@
                         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

issm/trunk/src/m/enum/EnumToModelField.m ¶

@@ -21 +21 @@
                 case InversionVyObsEnum(), string='vy_obs'; return
                 case BasalforcingsMeltingRateEnum(), string='basal_melting_rate'; return
-      case SurfaceforcingsAccumulationRateEnum(), string='surface_accumulation_rate'; return
-                case SurfaceforcingsAblationRateEnum(), string='surface_ablation_rate'; return
                 case SurfaceforcingsMassBalanceEnum(), string='surface_mass_balance'; return
                 otherwise, error(['Enum ' num2str(enum)  ' not found associated to any model field']);

issm/trunk/src/m/utils/BC/SetMarineIceSheetBC.m ¶

@@ -58 +58 @@
 %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,:);
@@ -106 +107 @@
         disp('      no thermal boundary conditions created: no observed temperature found');
 end
+