Changeset 13521


Ignore:
Timestamp:
10/03/12 11:51:45 (12 years ago)
Author:
helsen
Message:

update of the SMB gradients code

Location:
issm/trunk-jpl/src
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/EnumDefinitions/EnumDefinitions.h

    r13516 r13521  
    197197        SurfaceforcingsMonthlytemperaturesEnum,
    198198        SurfaceforcingsHcEnum,
     199        SurfaceforcingsHrefEnum,
     200        SurfaceforcingsSmbrefEnum,
    199201        SurfaceforcingsSmbPosMaxEnum,
    200202        SurfaceforcingsSmbPosMinEnum,

  • issm/trunk-jpl/src/c/classes/objects/Elements/Penta.cpp

    r13414 r13521  
    26572657        IssmDouble b_neg[NUMVERTICES];                          // Hs-SMB relation paremeter
    26582658        IssmDouble Hc[NUMVERTICES];                                     // elevation of transition between accumulation regime and ablation regime
     2659        IssmDouble Href[NUMVERTICES];                                   // reference elevation from which deviations are used to calculate the SMB adjustment
     2660        IssmDouble Smbref[NUMVERTICES];                         // reference SMB to which deviations are added
    26592661        IssmDouble smb_pos_max[NUMVERTICES];            // maximum SMB value in the accumulation regime
    26602662        IssmDouble smb_pos_min[NUMVERTICES];            // minimum SMB value in the accumulation regime
    … …  
    26672669        /*Recover SmbGradients*/
    26682670        GetInputListOnVertices(&Hc[0],SurfaceforcingsHcEnum);
     2671        GetInputListOnVertices(&Href[0],SurfaceforcingsHrefEnum);
     2672        GetInputListOnVertices(&Smbref[0],SurfaceforcingsSmbrefEnum);
    26692673        GetInputListOnVertices(&smb_pos_max[0],SurfaceforcingsSmbPosMaxEnum);
    26702674        GetInputListOnVertices(&smb_pos_min[0],SurfaceforcingsSmbPosMinEnum);
    … …  
    26832687                       
    26842688   // loop over all vertices
    2685    for(i=0;i<NUMVERTICES;i++){
     2689 for(i=0;i<NUMVERTICES;i++){
    26862690     if(s[i]>Hc[i]){
    2687             smb[i]=a_pos[i]+b_pos[i]*s[i];
    2688                  if(smb[i]>smb_pos_max[i]){smb[i]=smb_pos_max[i];}
    2689                  if(smb[i]<smb_pos_min[i]){smb[i]=smb_pos_min[i];}
     2691                  if(Href[i]>Hc[i]){smb[i]=Smbref[i]+b_pos[i]*(s[i]-Href[i]);}
     2692                  if(Href[i]<=Hc[i]){smb[i]=a_pos[i]+b_pos[i]*s[i];}
     2693                  if(smb[i]>smb_pos_max[i]){smb[i]=smb_pos_max[i];}
     2694                  if(smb[i]<smb_pos_min[i]){smb[i]=smb_pos_min[i];}
    26902695          }
    26912696          else{
    2692             smb[i]=a_neg[i]+b_neg[i]*s[i];
     2697                  if(Href[i]>Hc[i]){smb[i]=a_neg[i]+b_neg[i]*s[i];}
     2698                  if(Href[i]<=Hc[i]){smb[i]=Smbref[i]+b_neg[i]*(s[i]-Href[i]);}
    26932699          }
    26942700          smb[i]=smb[i]/rho_ice;      // SMB in m/y ice         
    2695          
    26962701        }  //end of the loop over the vertices
    26972702          /*Update inputs*/
    … …  
    36543659                        vel=sqrt(pow(vx,2.)+pow(vy,2.)+pow(vz,2.))+1.e-14;
    36553660                        h=sqrt( pow(hx*vx/vel,2.) + pow(hy*vy/vel,2.) + pow(hz*vz/vel,2.));
    3656 
    3657                         K[0][0]=h/(2*vel)*fabs(vx*vx);  K[0][1]=h/(2*vel)*fabs(vx*vy); K[0][2]=h/(2*vel)*fabs(vx*vz);
    3658                         K[1][0]=h/(2*vel)*fabs(vy*vx);  K[1][1]=h/(2*vel)*fabs(vy*vy); K[1][2]=h/(2*vel)*fabs(vy*vz);
    3659                         K[2][0]=h/(2*vel)*fabs(vz*vx);  K[2][1]=h/(2*vel)*fabs(vz*vy); K[2][2]=h/(2*vel)*fabs(vz*vz);
     3661                        vz=0;
     3662                        K[0][0]=h/(2*vel)*(vx*vx);  K[0][1]=h/(2*vel)*(vx*vy); K[0][2]=h/(2*vel)*(vx*vz);
     3663                        K[1][0]=h/(2*vel)*(vy*vx);  K[1][1]=h/(2*vel)*(vy*vy); K[1][2]=h/(2*vel)*(vy*vz);
     3664                        K[2][0]=h/(2*vel)*(vz*vx);  K[2][1]=h/(2*vel)*(vz*vy); K[2][2]=h/(2*vel)*(vz*vz);
    36603665
    36613666                        D_scalar_stab=gauss->weight*Jdet;
    … …  
    43144319        IssmDouble B_average,s_average;
    43154320        int*   doflist=NULL;
    4316         //IssmDouble pressure[numdof];
     4321        IssmDouble pressure[numdof];
    43174322
    43184323        /*Get dof list: */
    … …  
    43224327        for(i=0;i<numdof;i++){
    43234328                values[i]=solution[doflist[i]];
    4324                 //GetInputListOnVertices(&pressure[0],PressureEnum);
    4325                 //if(values[i]>matpar->TMeltingPoint(pressure[i])) values[i]=matpar->TMeltingPoint(pressure[i]);
    4326                 //if(values[i]<matpar->TMeltingPoint(pressure[i])-50) values[i]=matpar->TMeltingPoint(pressure[i])-50;
     4329                GetInputListOnVertices(&pressure[0],PressureEnum);
     4330                if(values[i]>matpar->TMeltingPoint(pressure[i])) values[i]=matpar->TMeltingPoint(pressure[i]);
     4331                if(values[i]<matpar->TMeltingPoint(pressure[i])-50) values[i]=matpar->TMeltingPoint(pressure[i])-50;
    43274332
    43284333                /*Check solution*/

  • issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp

    r13415 r13521  
    22632263        IssmDouble b_neg[NUMVERTICES];                          // Hs-SMB relation paremeter
    22642264        IssmDouble Hc[NUMVERTICES];                                     // elevation of transition between accumulation regime and ablation regime
     2265        IssmDouble Href[NUMVERTICES];                                   // reference elevation from which deviations are used to calculate the SMB adjustment
     2266        IssmDouble Smbref[NUMVERTICES];                         // reference SMB to which deviations are added
    22652267        IssmDouble smb_pos_max[NUMVERTICES];            // maximum SMB value in the accumulation regime
    22662268        IssmDouble smb_pos_min[NUMVERTICES];            // minimum SMB value in the accumulation regime
    … …  
    22732275        /*Recover SmbGradients*/
    22742276        GetInputListOnVertices(&Hc[0],SurfaceforcingsHcEnum);
     2277        GetInputListOnVertices(&Href[0],SurfaceforcingsHrefEnum);
     2278        GetInputListOnVertices(&Smbref[0],SurfaceforcingsSmbrefEnum);
    22752279        GetInputListOnVertices(&smb_pos_max[0],SurfaceforcingsSmbPosMaxEnum);
    22762280        GetInputListOnVertices(&smb_pos_min[0],SurfaceforcingsSmbPosMinEnum);
    … …  
    22912295   for(i=0;i<NUMVERTICES;i++){
    22922296     if(s[i]>Hc[i]){
    2293             smb[i]=a_pos[i]+b_pos[i]*s[i];
    2294                  if(smb[i]>smb_pos_max[i]){smb[i]=smb_pos_max[i];}
    2295                  if(smb[i]<smb_pos_min[i]){smb[i]=smb_pos_min[i];}
     2297                  if(Href[i]>Hc[i]){smb[i]=Smbref[i]+b_pos[i]*(s[i]-Href[i]);}
     2298                  if(Href[i]<=Hc[i]){smb[i]=a_pos[i]+b_pos[i]*s[i];}
     2299                  if(smb[i]>smb_pos_max[i]){smb[i]=smb_pos_max[i];}
     2300                  if(smb[i]<smb_pos_min[i]){smb[i]=smb_pos_min[i];}
    22962301          }
    22972302          else{
    2298             smb[i]=a_neg[i]+b_neg[i]*s[i];
     2303                  if(Href[i]>Hc[i]){smb[i]=a_neg[i]+b_neg[i]*s[i];}
     2304                  if(Href[i]<=Hc[i]){smb[i]=Smbref[i]+b_neg[i]*(s[i]-Href[i]);}
    22992305          }
    23002306          smb[i]=smb[i]/rho_ice;      // SMB in m/y ice         
    2301         }  //end of the loop over the vertices
     2307  /*   printf("s %e \n",s[i]);
     2308     printf("Hsref %e \n",Href[i]);
     2309     printf("Hc %e \n",Hc[i]);
     2310     printf("Smbref %e \n",Smbref[i]);
     2311     printf("b_neg %e \n",b_neg[i]);
     2312     printf("smb %e \n",smb[i]);
     2313          _error_("stop-in-code"); */
     2314                }  //end of the loop over the vertices
    23022315          /*Update inputs*/
    23032316          this->inputs->AddInput(new TriaP1Input(SurfaceforcingsMassBalanceEnum,&smb[0]));

  • issm/trunk-jpl/src/c/modules/EnumToStringx/EnumToStringx.cpp

    r13516 r13521  
    202202                case SurfaceforcingsMonthlytemperaturesEnum : return "SurfaceforcingsMonthlytemperatures";
    203203                case SurfaceforcingsHcEnum : return "SurfaceforcingsHc";
     204                case SurfaceforcingsHrefEnum : return "SurfaceforcingsHref";
     205                case SurfaceforcingsSmbrefEnum : return "SurfaceforcingsSmbref";
    204206                case SurfaceforcingsSmbPosMaxEnum : return "SurfaceforcingsSmbPosMax";
    205207                case SurfaceforcingsSmbPosMinEnum : return "SurfaceforcingsSmbPosMin";

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/Prognostic/UpdateElementsPrognostic.cpp

    r13073 r13521  
    7373        if(issmbgradients){
    7474                iomodel->FetchDataToInput(elements,SurfaceforcingsHcEnum);
     75                iomodel->FetchDataToInput(elements,SurfaceforcingsHrefEnum);
     76                iomodel->FetchDataToInput(elements,SurfaceforcingsSmbrefEnum);
    7577                iomodel->FetchDataToInput(elements,SurfaceforcingsSmbPosMaxEnum);
    7678                iomodel->FetchDataToInput(elements,SurfaceforcingsSmbPosMinEnum);

  • issm/trunk-jpl/src/c/modules/StringToEnumx/StringToEnumx.cpp

    r13516 r13521  
    206206              else if (strcmp(name,"SurfaceforcingsMonthlytemperatures")==0) return SurfaceforcingsMonthlytemperaturesEnum;
    207207              else if (strcmp(name,"SurfaceforcingsHc")==0) return SurfaceforcingsHcEnum;
     208              else if (strcmp(name,"SurfaceforcingsHref")==0) return SurfaceforcingsHrefEnum;
     209              else if (strcmp(name,"SurfaceforcingsSmbref")==0) return SurfaceforcingsSmbrefEnum;
    208210              else if (strcmp(name,"SurfaceforcingsSmbPosMax")==0) return SurfaceforcingsSmbPosMaxEnum;
    209211              else if (strcmp(name,"SurfaceforcingsSmbPosMin")==0) return SurfaceforcingsSmbPosMinEnum;
    … …  
    259261              else if (strcmp(name,"PrognosticSolution")==0) return PrognosticSolutionEnum;
    260262              else if (strcmp(name,"SteadystateSolution")==0) return SteadystateSolutionEnum;
    261               else if (strcmp(name,"SurfaceSlopeAnalysis")==0) return SurfaceSlopeAnalysisEnum;
    262               else if (strcmp(name,"SurfaceSlopeSolution")==0) return SurfaceSlopeSolutionEnum;
    263263         else stage=3;
    264264   }
    265265   if(stage==3){
    266               if (strcmp(name,"SurfaceSlopeXAnalysis")==0) return SurfaceSlopeXAnalysisEnum;
     266              if (strcmp(name,"SurfaceSlopeAnalysis")==0) return SurfaceSlopeAnalysisEnum;
     267              else if (strcmp(name,"SurfaceSlopeSolution")==0) return SurfaceSlopeSolutionEnum;
     268              else if (strcmp(name,"SurfaceSlopeXAnalysis")==0) return SurfaceSlopeXAnalysisEnum;
    267269              else if (strcmp(name,"SurfaceSlopeYAnalysis")==0) return SurfaceSlopeYAnalysisEnum;
    268270              else if (strcmp(name,"ThermalAnalysis")==0) return ThermalAnalysisEnum;
    … …  
    382384              else if (strcmp(name,"SurfaceArea")==0) return SurfaceAreaEnum;
    383385              else if (strcmp(name,"SurfaceAverageVelMisfit")==0) return SurfaceAverageVelMisfitEnum;
    384               else if (strcmp(name,"SurfaceLogVelMisfit")==0) return SurfaceLogVelMisfitEnum;
    385               else if (strcmp(name,"SurfaceLogVxVyMisfit")==0) return SurfaceLogVxVyMisfitEnum;
    386386         else stage=4;
    387387   }
    388388   if(stage==4){
    389               if (strcmp(name,"SurfaceRelVelMisfit")==0) return SurfaceRelVelMisfitEnum;
     389              if (strcmp(name,"SurfaceLogVelMisfit")==0) return SurfaceLogVelMisfitEnum;
     390              else if (strcmp(name,"SurfaceLogVxVyMisfit")==0) return SurfaceLogVxVyMisfitEnum;
     391              else if (strcmp(name,"SurfaceRelVelMisfit")==0) return SurfaceRelVelMisfitEnum;
    390392              else if (strcmp(name,"SurfaceSlopeX")==0) return SurfaceSlopeXEnum;
    391393              else if (strcmp(name,"SurfaceSlopeY")==0) return SurfaceSlopeYEnum;
    … …  
    505507              else if (strcmp(name,"Separate")==0) return SeparateEnum;
    506508              else if (strcmp(name,"Sset")==0) return SsetEnum;
    507               else if (strcmp(name,"Verbose")==0) return VerboseEnum;
    508               else if (strcmp(name,"TriangleInterp")==0) return TriangleInterpEnum;
    509509         else stage=5;
    510510   }
    511511   if(stage==5){
    512               if (strcmp(name,"BilinearInterp")==0) return BilinearInterpEnum;
     512              if (strcmp(name,"Verbose")==0) return VerboseEnum;
     513              else if (strcmp(name,"TriangleInterp")==0) return TriangleInterpEnum;
     514              else if (strcmp(name,"BilinearInterp")==0) return BilinearInterpEnum;
    513515              else if (strcmp(name,"NearestInterp")==0) return NearestInterpEnum;
    514516              else if (strcmp(name,"XY")==0) return XYEnum;

  • issm/trunk-jpl/src/c/shared/Numerics/UnitConversion.cpp

    r13056 r13521  
    7272                case SurfaceforcingsANegEnum:                                           scale=yts;break;     //m/yr
    7373                case SurfaceforcingsBNegEnum:                                           scale=yts;break;     //m/yr
     74                case SurfaceforcingsSmbrefEnum:                                 scale=yts;break;     //m/yr
    7475                case MisfitEnum:                             scale=pow(yts,2);break; //(m/yr)^2
    7576                case MassFluxEnum:                           scale=pow((IssmDouble)10,-12)*yts;break; // (GigaTon/year)

  • issm/trunk-jpl/src/m/classes/surfaceforcings.m

    r13040 r13521  
    1212                isdelta18o = 0;
    1313                hc = NaN;
     14                href = NaN;
     15                smbref = NaN;
    1416                smb_pos_max = NaN;
    1517                smb_pos_min = NaN;
    … …  
    6062                                elseif(obj.issmbgradients)
    6163                                        md = checkfield(md,'surfaceforcings.hc','forcing',1,'NaN',1);
     64                                        md = checkfield(md,'surfaceforcings.href','forcing',1,'NaN',1);
     65                                        md = checkfield(md,'surfaceforcings.smbref','forcing',1,'NaN',1);
    6266                                        md = checkfield(md,'surfaceforcings.smb_pos_max','forcing',1,'NaN',1);
    6367                                        md = checkfield(md,'surfaceforcings.smb_pos_min','forcing',1,'NaN',1);
    … …  
    8993                        fielddisplay(obj,'issmbgradients','is smb gradients method activated (0 or 1, default is 0)');
    9094                        fielddisplay(obj,'hc',' elevation of intersection between accumulation and ablation regime required if smb gradients is activated');
     95                        fielddisplay(obj,'href',' reference elevation from which deviation is used to calculate SMB adjustment in smb gradients method');
     96                        fielddisplay(obj,'smbref',' reference smb from which deviation is calculated in smb gradients method');
    9197                        fielddisplay(obj,'smb_pos_max',' maximum value of positive smb required if smb gradients is activated');
    9298                        fielddisplay(obj,'smb_pos_min',' minimum value of positive smb required if smb gradients is activated');
    … …  
    117123                        if obj.issmbgradients,
    118124                                WriteData(fid,'object',obj,'fieldname','hc','format','DoubleMat','mattype',1);
     125                                WriteData(fid,'object',obj,'fieldname','href','format','DoubleMat','mattype',1);
     126                                WriteData(fid,'object',obj,'fieldname','smbref','format','DoubleMat','mattype',1);
    119127                                WriteData(fid,'object',obj,'fieldname','smb_pos_max','format','DoubleMat','mattype',1);
    120128                                WriteData(fid,'object',obj,'fieldname','smb_pos_min','format','DoubleMat','mattype',1);
Note: See TracChangeset for help on using the changeset viewer.