Changeset 21469


Ignore:
Timestamp:
01/09/17 10:17:06 (8 years ago)
Author:
hakesson
Message:

NEW: added SMB gradients ela method

Location:
issm/trunk-jpl
Files:
6 added
6 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp

    r21232 r21469  
    105105                        iomodel->FetchDataToInput(elements,"md.smb.b_neg",SmbBNegEnum);
    106106                        break;
     107                case SMBgradientselaEnum:
     108                        iomodel->FetchDataToInput(elements,"md.smb.ela",SmbElaEnum);
     109                        iomodel->FetchDataToInput(elements,"md.smb.b_pos",SmbBPosEnum);
     110                        iomodel->FetchDataToInput(elements,"md.smb.b_neg",SmbBNegEnum);
     111                        iomodel->FetchDataToInput(elements,"md.smb.b_max",SmbBMaxEnum);
     112                        iomodel->FetchDataToInput(elements,"md.smb.b_min",SmbBMinEnum);
     113                        break;
    107114                case SMBhenningEnum:
    108115                        iomodel->FetchDataToInput(elements,"md.smb.smbref",SmbSmbrefEnum,0.);
    … …  
    217224                        /*Nothing to add to parameters*/
    218225                        break;
     226                case SMBgradientselaEnum:
     227                        /*Nothing to add to parameters*/
     228                        break;
    219229                case SMBhenningEnum:
    220230                        /*Nothing to add to parameters*/
    … …  
    282292                        SmbGradientsx(femmodel);
    283293                        break;
     294                case SMBgradientselaEnum:
     295                        if(VerboseSolution())_printf0_("        call smb gradients ela module\n");
     296                        SmbGradientsElax(femmodel);
     297                        break;
    284298                case SMBhenningEnum:
    285299                        if(VerboseSolution())_printf0_("  call smb Henning module\n");

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

    r20690 r21469  
    1717        IssmDouble rho_water;                   // density of fresh water
    1818        IssmDouble rho_ice;                     // density of ice
     19        IssmDouble yts;                                                         // conversion factor year to second
    1920
    2021        /*Loop over all the elements of this partition*/
    … …  
    4445                rho_water=element->matpar->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
    4546
     47                /* Get constants */
     48                femmodel->parameters->FindParam(&yts,ConstantsYtsEnum);
     49
    4650                // loop over all vertices
    4751                for(v=0;v<numvertices;v++){
    … …  
    5256                                smb[v]=Smbref[v]+b_neg[v]*(s[v]-Href[v]);
    5357                        }
     58
    5459                        smb[v]=smb[v]/1000*rho_water/rho_ice;      // SMB in m/y ice
    5560                }  //end of the loop over the vertices
    … …  
    6368                xDelete<IssmDouble>(s);
    6469                xDelete<IssmDouble>(smb);
     70        }
     71
     72}/*}}}*/
     73void SmbGradientsElax(FemModel* femmodel){/*{{{*/
     74
     75        // void SurfaceMassBalancex(hd,agd,ni){
     76        //    INPUT parameters: ni: working size of arrays
     77        //    INPUT: surface elevation (m): hd(NA)
     78        //    OUTPUT: surface mass-balance (m/yr ice): agd(NA)
     79        int v;
     80
     81        /*Loop over all the elements of this partition*/
     82        for(int i=0;i<femmodel->elements->Size();i++){
     83                Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
     84
     85                /*Allocate all arrays*/
     86                int         numvertices = element->GetNumberOfVertices();
     87                IssmDouble* ela       = xNew<IssmDouble>(numvertices); // Equilibrium Line Altitude (m a.s.l) to which deviations are used to calculate the SMB
     88                IssmDouble* b_pos       = xNew<IssmDouble>(numvertices); // SMB gradient above ELA (m ice eq. per m elevation change)
     89                IssmDouble* b_neg       = xNew<IssmDouble>(numvertices); // SMB gradient below ELA (m ice eq. per m elevation change)
     90                IssmDouble* b_max       = xNew<IssmDouble>(numvertices); // Upper cap on SMB rate (m/y ice eq.)
     91                IssmDouble* b_min       = xNew<IssmDouble>(numvertices); // Lower cap on SMB rate (m/y ice eq.)
     92                IssmDouble* s           = xNew<IssmDouble>(numvertices); // Surface elevation (m a.s.l.)
     93                IssmDouble* smb         = xNew<IssmDouble>(numvertices); // SMB (m/y ice eq.)
     94
     95                /*Recover ELA, SMB gradients, and caps*/
     96                element->GetInputListOnVertices(ela,SmbElaEnum);
     97                element->GetInputListOnVertices(b_pos,SmbBPosEnum);
     98                element->GetInputListOnVertices(b_neg,SmbBNegEnum);
     99                element->GetInputListOnVertices(b_max,SmbBMaxEnum);
     100                element->GetInputListOnVertices(b_min,SmbBMinEnum);
     101
     102                /*Recover surface elevation at vertices: */
     103                element->GetInputListOnVertices(s,SurfaceEnum);
     104
     105                /*Loop over all vertices, calculate SMB*/
     106                for(v=0;v<numvertices;v++){
     107                        // if surface is above the ELA
     108                        if(s[v]>ela[v]){               
     109                                smb[v]=b_pos[v]*(s[v]-ela[v]);
     110                        }
     111                        // if surface is below or equal to the ELA
     112                        else{
     113                                smb[v]=b_neg[v]*(s[v]-ela[v]);
     114                        }
     115
     116                        // if SMB is larger than upper cap, set SMB to upper cap
     117                        if(smb[v]>b_max[v]){
     118                                smb[v]=b_max[v];
     119                        }
     120                        // if SMB is smaller than lower cap, set SMB to lower cap
     121                        if(smb[v]<b_min[v]){
     122                                smb[v]=b_min[v];
     123                        }
     124                }  //end of the loop over the vertices
     125
     126                /*Add input to element and Free memory*/
     127                element->AddInput(SmbMassBalanceEnum,smb,P1Enum);
     128                xDelete<IssmDouble>(ela);
     129                xDelete<IssmDouble>(b_pos);
     130                xDelete<IssmDouble>(b_neg);
     131                xDelete<IssmDouble>(b_max);
     132                xDelete<IssmDouble>(b_min);
     133                xDelete<IssmDouble>(s);
     134                xDelete<IssmDouble>(smb);
     135
    65136        }
    66137

  • issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h

    r21216 r21469  
    1111void SurfaceMassBalancex(FemModel* femmodel);
    1212void SmbGradientsx(FemModel* femmodel);
     13void SmbGradientsElax(FemModel* femmodel);
    1314void Delta18oParameterizationx(FemModel* femmodel);
    1415void MungsmtpParameterizationx(FemModel* femmodel);

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

    r21463 r21469  
    464464        SmbRefreezeEnum,
    465465        SMBgcmEnum,
     466        /*SMBgradientsela*/
     467        SMBgradientselaEnum,
     468        SmbElaEnum,
     469        SmbBMaxEnum,
     470        SmbBMinEnum,
    466471        /*}}}*/
    467472        /*Inputs {{{*/

  • issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

    r21464 r21469  
    465465                case SmbRefreezeEnum : return "SmbRefreeze";
    466466                case SMBgcmEnum : return "SMBgcm";
     467                case SMBgradientselaEnum : return "SMBgradientsela";
     468                case SmbElaEnum : return "SmbEla";
     469                case SmbBMaxEnum : return "SmbBMax";
     470                case SmbBMinEnum : return "SmbBMin";
    467471                case AdjointpEnum : return "Adjointp";
    468472                case AdjointxEnum : return "Adjointx";

  • issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

    r21464 r21469  
    474474              else if (strcmp(name,"SmbRefreeze")==0) return SmbRefreezeEnum;
    475475              else if (strcmp(name,"SMBgcm")==0) return SMBgcmEnum;
     476              else if (strcmp(name,"SMBgradientsela")==0) return SMBgradientselaEnum;
     477              else if (strcmp(name,"SmbEla")==0) return SmbElaEnum;
     478              else if (strcmp(name,"SmbBMax")==0) return SmbBMaxEnum;
     479              else if (strcmp(name,"SmbBMin")==0) return SmbBMinEnum;
    476480              else if (strcmp(name,"Adjointp")==0) return AdjointpEnum;
    477481              else if (strcmp(name,"Adjointx")==0) return AdjointxEnum;
    … …  
    502506              else if (strcmp(name,"TemperaturePicard")==0) return TemperaturePicardEnum;
    503507              else if (strcmp(name,"TemperaturePDD")==0) return TemperaturePDDEnum;
    504               else if (strcmp(name,"ThicknessAbsMisfit")==0) return ThicknessAbsMisfitEnum;
     508         else stage=5;
     509   }
     510   if(stage==5){
     511              if (strcmp(name,"ThicknessAbsMisfit")==0) return ThicknessAbsMisfitEnum;
    505512              else if (strcmp(name,"SurfaceAbsMisfit")==0) return SurfaceAbsMisfitEnum;
    506513              else if (strcmp(name,"Vel")==0) return VelEnum;
    507514              else if (strcmp(name,"Velocity")==0) return VelocityEnum;
    508          else stage=5;
    509    }
    510    if(stage==5){
    511               if (strcmp(name,"VxAverage")==0) return VxAverageEnum;
     515              else if (strcmp(name,"VxAverage")==0) return VxAverageEnum;
    512516              else if (strcmp(name,"Vx")==0) return VxEnum;
    513517              else if (strcmp(name,"VxPicard")==0) return VxPicardEnum;
    … …  
    625629              else if (strcmp(name,"Outputdefinition40")==0) return Outputdefinition40Enum;
    626630              else if (strcmp(name,"Outputdefinition41")==0) return Outputdefinition41Enum;
    627               else if (strcmp(name,"Outputdefinition42")==0) return Outputdefinition42Enum;
     631         else stage=6;
     632   }
     633   if(stage==6){
     634              if (strcmp(name,"Outputdefinition42")==0) return Outputdefinition42Enum;
    628635              else if (strcmp(name,"Outputdefinition43")==0) return Outputdefinition43Enum;
    629636              else if (strcmp(name,"Outputdefinition44")==0) return Outputdefinition44Enum;
    630637              else if (strcmp(name,"Outputdefinition45")==0) return Outputdefinition45Enum;
    631          else stage=6;
    632    }
    633    if(stage==6){
    634               if (strcmp(name,"Outputdefinition46")==0) return Outputdefinition46Enum;
     638              else if (strcmp(name,"Outputdefinition46")==0) return Outputdefinition46Enum;
    635639              else if (strcmp(name,"Outputdefinition47")==0) return Outputdefinition47Enum;
    636640              else if (strcmp(name,"Outputdefinition48")==0) return Outputdefinition48Enum;
    … …  
    748752              else if (strcmp(name,"Gsl")==0) return GslEnum;
    749753              else if (strcmp(name,"Cuffey")==0) return CuffeyEnum;
    750               else if (strcmp(name,"BuddJacka")==0) return BuddJackaEnum;
     754         else stage=7;
     755   }
     756   if(stage==7){
     757              if (strcmp(name,"BuddJacka")==0) return BuddJackaEnum;
    751758              else if (strcmp(name,"CuffeyTemperate")==0) return CuffeyTemperateEnum;
    752759              else if (strcmp(name,"Paterson")==0) return PatersonEnum;
    753760              else if (strcmp(name,"Arrhenius")==0) return ArrheniusEnum;
    754          else stage=7;
    755    }
    756    if(stage==7){
    757               if (strcmp(name,"LliboutryDuval")==0) return LliboutryDuvalEnum;
     761              else if (strcmp(name,"LliboutryDuval")==0) return LliboutryDuvalEnum;
    758762              else if (strcmp(name,"ExtrapolationVariable")==0) return ExtrapolationVariableEnum;
    759763              else if (strcmp(name,"IceMaskNodeActivation")==0) return IceMaskNodeActivationEnum;
    … …  
    871875              else if (strcmp(name,"BalancethicknessSolution")==0) return BalancethicknessSolutionEnum;
    872876              else if (strcmp(name,"Balancethickness2Analysis")==0) return Balancethickness2AnalysisEnum;
    873               else if (strcmp(name,"Balancethickness2Solution")==0) return Balancethickness2SolutionEnum;
     877         else stage=8;
     878   }
     879   if(stage==8){
     880              if (strcmp(name,"Balancethickness2Solution")==0) return Balancethickness2SolutionEnum;
    874881              else if (strcmp(name,"BalancethicknessSoftAnalysis")==0) return BalancethicknessSoftAnalysisEnum;
    875882              else if (strcmp(name,"BalancethicknessSoftSolution")==0) return BalancethicknessSoftSolutionEnum;
    876883              else if (strcmp(name,"BalancevelocityAnalysis")==0) return BalancevelocityAnalysisEnum;
    877          else stage=8;
    878    }
    879    if(stage==8){
    880               if (strcmp(name,"BalancevelocitySolution")==0) return BalancevelocitySolutionEnum;
     884              else if (strcmp(name,"BalancevelocitySolution")==0) return BalancevelocitySolutionEnum;
    881885              else if (strcmp(name,"L2ProjectionEPLAnalysis")==0) return L2ProjectionEPLAnalysisEnum;
    882886              else if (strcmp(name,"L2ProjectionBaseAnalysis")==0) return L2ProjectionBaseAnalysisEnum;
    … …  
    994998              else if (strcmp(name,"Parameters")==0) return ParametersEnum;
    995999              else if (strcmp(name,"Vertices")==0) return VerticesEnum;
    996               else if (strcmp(name,"Results")==0) return ResultsEnum;
     1000         else stage=9;
     1001   }
     1002   if(stage==9){
     1003              if (strcmp(name,"Results")==0) return ResultsEnum;
    9971004              else if (strcmp(name,"MaximumNumberOfDefinitions")==0) return MaximumNumberOfDefinitionsEnum;
    998          else stage=9;
     1005         else stage=10;
    9991006   }
    10001007        /*If we reach this point, the string provided has not been found*/
Note: See TracChangeset for help on using the changeset viewer.