Changeset 22330


Ignore:
Timestamp:
01/04/18 15:14:01 (7 years ago)
Author:
seroussi
Message:

NEW: finished implementing scaled correction functions

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

Legend:

Unmodified
Added
Removed

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

    r22326 r22330  
    690690        /*Intermediaries*/
    691691        int         domaintype;
    692         IssmDouble  phi,base_area;
     692        IssmDouble  phi,base_area,scalefactor,floatingarea;
    693693        IssmDouble  xyz_list[NUMVERTICES][3];
    694694
    … …  
    703703        base_area= 1./2.*fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1]));
    704704
     705        floatingarea=(1-phi)*base_area;
     706
     707        if(scaled==true){
     708                Input* scalefactor_input = inputs->GetInput(MeshScaleFactorEnum); _assert_(scalefactor_input);
     709                scalefactor_input->GetInputAverage(&scalefactor);
     710                floatingarea=floatingarea*scalefactor;
     711        }
     712
    705713        /*Clean up and return*/
    706         return (1-phi)*base_area;
     714        return floatingarea;
    707715}
    708716/*}}}*/
    … …  
    11641172        /*Intermediaries*/
    11651173        int         domaintype;
    1166         IssmDouble  phi,base_area;
     1174        IssmDouble  phi,base_area,scalefactor,groundedarea;
    11671175        IssmDouble  xyz_list[NUMVERTICES][3];
    11681176
    … …  
    11771185        base_area= 1./2.*fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1]));
    11781186
     1187        groundedarea=phi*base_area;
     1188
     1189        if(scaled==true){
     1190                Input* scalefactor_input = inputs->GetInput(MeshScaleFactorEnum); _assert_(scalefactor_input);
     1191                scalefactor_input->GetInputAverage(&scalefactor);
     1192                groundedarea=groundedarea*scalefactor;
     1193        }
    11791194        /*Clean up and return*/
    1180         return phi*base_area;
     1195        return groundedarea;
    11811196}
    11821197/*}}}*/
    … …  
    11841199
    11851200        /*The volume of a troncated prism is base * 1/3 sum(length of edges)*/
    1186         IssmDouble base,height;
     1201        IssmDouble base,height,scalefactor;
    11871202        IssmDouble xyz_list[NUMVERTICES][3];
    11881203
    … …  
    11961211        base = 1./2.*fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1]));
    11971212
     1213        if(scaled==true){ //scale for area projection correction
     1214                Input* scalefactor_input = inputs->GetInput(MeshScaleFactorEnum); _assert_(scalefactor_input);
     1215                scalefactor_input->GetInputAverage(&scalefactor);
     1216                base=base*scalefactor;
     1217        }
     1218
    11981219        /*Now get the average height*/
    11991220        height = 1./3.*((xyz_list[3][2]-xyz_list[0][2])+(xyz_list[4][2]-xyz_list[1][2])+(xyz_list[5][2]-xyz_list[2][2]));
    … …  
    12071228        /*Volume above floatation: H + rho_water/rho_ice*bathymetry for nodes on the bed*/
    12081229        IssmDouble rho_ice,rho_water;
    1209         IssmDouble base,bed,surface,bathymetry;
     1230        IssmDouble base,bed,surface,bathymetry,scalefactor;
    12101231        IssmDouble xyz_list[NUMVERTICES][3];
    12111232
    … …  
    12201241         * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/
    12211242        base = 1./2.*fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1]));
     1243        if(scaled==true){
     1244                Input* scalefactor_input = inputs->GetInput(MeshScaleFactorEnum); _assert_(scalefactor_input);
     1245                scalefactor_input->GetInputAverage(&scalefactor);
     1246                base=base*scalefactor;
     1247        }
    12221248
    12231249        /*Now get the average height above floatation*/
    … …  
    26392665        bool       mainlyfloating;
    26402666        IssmDouble fbmb=0;
    2641         IssmDouble rho_ice,fraction1,fraction2,floatingmelt,Jdet;
     2667        IssmDouble rho_ice,fraction1,fraction2,floatingmelt,Jdet,scalefactor;
    26422668        IssmDouble Total_Fbmb=0;
    26432669        IssmDouble xyz_list[NUMVERTICES][3];
    … …  
    26502676        Input* floatingmelt_input = this->GetInput(BasalforcingsFloatingiceMeltingRateEnum); _assert_(floatingmelt_input);
    26512677        Input* gllevelset_input = this->GetInput(MaskGroundediceLevelsetEnum); _assert_(gllevelset_input);
     2678        Input* scalefactor_input = NULL;
     2679        if(scaled==true){
     2680                scalefactor_input = this->GetInput(MeshScaleFactorEnum); _assert_(scalefactor_input);
     2681        }
    26522682        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    26532683
    … …  
    26602690                this->JacobianDeterminantBase(&Jdet,&xyz_list[0][0],gauss);
    26612691                floatingmelt_input->GetInputValue(&floatingmelt,gauss);
    2662                 fbmb+=floatingmelt*Jdet*gauss->weight;
     2692                if(scaled==true){
     2693                        scalefactor_input->GetInputValue(&scalefactor,gauss);
     2694                }
     2695                else scalefactor=1;
     2696                fbmb+=floatingmelt*Jdet*gauss->weight*scalefactor;
    26632697        }
    26642698
    … …  
    26762710        bool       mainlyfloating;
    26772711        IssmDouble gbmb=0;
    2678         IssmDouble rho_ice,fraction1,fraction2,groundedmelt,Jdet;
     2712        IssmDouble rho_ice,fraction1,fraction2,groundedmelt,Jdet,scalefactor;
    26792713        IssmDouble Total_Gbmb=0;
    26802714        IssmDouble xyz_list[NUMVERTICES][3];
    … …  
    26872721        Input* groundedmelt_input = this->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(groundedmelt_input);
    26882722        Input* gllevelset_input = this->GetInput(MaskGroundediceLevelsetEnum); _assert_(gllevelset_input);
     2723        Input* scalefactor_input = NULL;
     2724        if(scaled==true){
     2725                scalefactor_input = this->GetInput(MeshScaleFactorEnum); _assert_(scalefactor_input);
     2726        }
    26892727        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    26902728
    … …  
    26972735                this->JacobianDeterminantBase(&Jdet,&xyz_list[0][0],gauss);
    26982736                groundedmelt_input->GetInputValue(&groundedmelt,gauss);
    2699                 gbmb+=groundedmelt*Jdet*gauss->weight;
     2737                if(scaled==true){
     2738                        scalefactor_input->GetInputValue(&scalefactor,gauss);
     2739                }
     2740                else scalefactor=1;
     2741                gbmb+=groundedmelt*Jdet*gauss->weight*scalefactor;
    27002742        }
    27012743
    … …  
    27102752
    27112753        /*The smb[Gt yr-1] of one element is area[m2] * smb [ m ice yr^-1] * rho_ice [kg m-3] / 1e+10^12 */
    2712         IssmDouble base,smb,rho_ice;
     2754        IssmDouble base,smb,rho_ice,scalefactor;
    27132755        IssmDouble Total_Smb=0;
    27142756        IssmDouble xyz_list[NUMVERTICES][3];
    … …  
    27302772
    27312773        smb_input->GetInputAverage(&smb);
    2732         Total_Smb=rho_ice*base*smb;// smb on element in kg s-1
     2774        if(scaled==true){
     2775                Input* scalefactor_input = inputs->GetInput(MeshScaleFactorEnum); _assert_(scalefactor_input);
     2776                scalefactor_input->GetInputAverage(&scalefactor);// average scalefactor on element
     2777        }
     2778        else{
     2779                scalefactor=1.;
     2780        }
     2781        Total_Smb=rho_ice*base*smb*scalefactor;// smb on element in kg s-1
    27332782
    27342783        /*Return: */

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

    r22326 r22330  
    894894        /*Intermediaries*/
    895895        int         domaintype;
    896         IssmDouble  phi;
     896        IssmDouble  phi,scalefactor,floatingarea;
    897897        IssmDouble *xyz_list  = NULL;
    898898
    … …  
    905905        this->GetVerticesCoordinates(&xyz_list);
    906906        phi=this->GetGroundedPortion(xyz_list);
     907        floatingarea=(1-phi)*this->GetArea();
     908        if(scaled==true){
     909                Input* scalefactor_input = inputs->GetInput(MeshScaleFactorEnum); _assert_(scalefactor_input);
     910                scalefactor_input->GetInputAverage(&scalefactor);
     911                floatingarea=floatingarea*scalefactor;
     912        }
    907913
    908914        /*Clean up and return*/
    909915        xDelete<IssmDouble>(xyz_list);
    910         return (1-phi)*this->GetArea();
     916        return floatingarea;
    911917}
    912918/*}}}*/
    … …  
    16141620        /*Intermediaries*/
    16151621        int         domaintype;
    1616         IssmDouble  phi;
     1622        IssmDouble  phi,scalefactor,groundedarea;
    16171623        IssmDouble *xyz_list  = NULL;
    16181624
    … …  
    16251631        this->GetVerticesCoordinates(&xyz_list);
    16261632        phi=this->GetGroundedPortion(xyz_list);
     1633        groundedarea=phi*this->GetArea();
     1634        if(scaled==true){
     1635                Input* scalefactor_input = inputs->GetInput(MeshScaleFactorEnum); _assert_(scalefactor_input);
     1636                scalefactor_input->GetInputAverage(&scalefactor);
     1637                groundedarea=groundedarea*scalefactor;
     1638        }
    16271639
    16281640        /*Clean up and return*/
    16291641        xDelete<IssmDouble>(xyz_list);
    1630         return phi*this->GetArea();
     1642        return groundedarea;
    16311643}
    16321644/*}}}*/
    … …  
    16791691        /*Intermediaries*/
    16801692        int i, numiceverts;
    1681         IssmDouble area_base,surface,base,Haverage;
     1693        IssmDouble area_base,surface,base,Haverage,scalefactor;
    16821694        IssmDouble Haux[NUMVERTICES], surfaces[NUMVERTICES], bases[NUMVERTICES];
     1695        IssmDouble SFaux[NUMVERTICES], scalefactors[NUMVERTICES];
    16831696        IssmDouble s[2]; // s:fraction of intersected triangle edges, that lies inside ice
    16841697        int* indices=NULL;
    16851698        IssmDouble* H=NULL;
     1699        IssmDouble* SF=NULL;
    16861700
    16871701        if(!IsIceInElement())return 0.;
    … …  
    16911705
    16921706        if(false && IsIcefront()){
    1693                 area_base=this->GetAreaIce();
    16941707                //Assumption: linear ice thickness profile on element.
    16951708                //Hence ice thickness at intersection of levelset function with triangle edge is linear interpolation of ice thickness at vertices.
    16961709                this->GetLevelsetIntersection(&indices, &numiceverts, s, MaskIceLevelsetEnum, 0.);
     1710                int numthk=numiceverts+2;
     1711                H=xNew<IssmDouble>(numthk);
     1712                //Correct area distortion caused by projection if requestion
     1713                area_base=this->GetAreaIce();
     1714                if(scaled==true){
     1715                        GetInputListOnVertices(&scalefactors[0],MeshScaleFactorEnum);
     1716                        for(i=0;i<NUMVERTICES;i++) SFaux[i]= scalefactors[indices[i]]; //sort thicknesses in ice/noice
     1717                        switch(numiceverts){
     1718                                case 1: // average over triangle
     1719                                        SF[0]=SFaux[0];
     1720                                        SF[1]=SFaux[0]+s[0]*(SFaux[1]-SFaux[0]);
     1721                                        SF[2]=SFaux[0]+s[1]*(SFaux[2]-SFaux[0]);
     1722                                        break;
     1723                                case 2: // average over quadrangle
     1724                                        SF[0]=SFaux[0];
     1725                                        SF[1]=SFaux[1];
     1726                                        SF[2]=SFaux[0]+s[0]*(SFaux[2]-SFaux[0]);
     1727                                        SF[3]=SFaux[1]+s[1]*(SFaux[2]-SFaux[1]);
     1728                                        break;
     1729                                default:
     1730                                        _error_("Number of ice covered vertices wrong in Tria::IceVolume()");
     1731                                        break;
     1732                        }
     1733                        scalefactor=0.;
     1734                        for(i=0;i<numthk;i++)   scalefactor+=SF[i];
     1735                        scalefactor/=IssmDouble(numthk);
     1736                        area_base=area_base*scalefactor;
     1737                }
    16971738                GetInputListOnVertices(&surfaces[0],SurfaceEnum);
    16981739                GetInputListOnVertices(&bases[0],BaseEnum);
    16991740                for(i=0;i<NUMVERTICES;i++) Haux[i]= surfaces[indices[i]]-bases[indices[i]]; //sort thicknesses in ice/noice
    1700                 int numthk=numiceverts+2;
    1701                 H=xNew<IssmDouble>(numthk);
    17021741                switch(numiceverts){
    17031742                        case 1: // average over triangle
    … …  
    17231762                /*First get back the area of the base*/
    17241763                area_base=this->GetArea();
     1764                if(scaled==true){
     1765                        Input* scalefactor_input = inputs->GetInput(MeshScaleFactorEnum); _assert_(scalefactor_input);
     1766                        scalefactor_input->GetInputAverage(&scalefactor);
     1767                        area_base=area_base*scalefactor;
     1768                }
    17251769
    17261770                /*Now get the average height*/
    … …  
    17351779        xDelete<int>(indices);
    17361780        xDelete<IssmDouble>(H);
     1781        xDelete<IssmDouble>(SF);
    17371782
    17381783        if(domaintype==Domain2DverticalEnum){
    … …  
    17481793        /*The volume above floatation: H + rho_water/rho_ice * bathymetry */
    17491794        IssmDouble rho_ice,rho_water;
    1750         IssmDouble base,surface,bed,bathymetry;
     1795        IssmDouble base,surface,bed,bathymetry,scalefactor;
    17511796        IssmDouble xyz_list[NUMVERTICES][3];
    17521797
    … …  
    17611806         * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/
    17621807        base = 1./2. * fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1]));
     1808        if(scaled==true){
     1809                Input* scalefactor_input = inputs->GetInput(MeshScaleFactorEnum); _assert_(scalefactor_input);
     1810                scalefactor_input->GetInputAverage(&scalefactor);
     1811                base=base*scalefactor;
     1812        }
    17631813
    17641814        /*Now get the average height and bathymetry*/
    … …  
    31213171        bool       mainlyfloating;
    31223172        IssmDouble fbmb=0;
    3123         IssmDouble rho_ice,fraction1,fraction2,floatingmelt,Jdet;
     3173        IssmDouble rho_ice,fraction1,fraction2,floatingmelt,Jdet,scalefactor;
    31243174        IssmDouble Total_Fbmb=0;
    31253175        IssmDouble xyz_list[NUMVERTICES][3];
    … …  
    31323182        Input* floatingmelt_input = this->GetInput(BasalforcingsFloatingiceMeltingRateEnum); _assert_(floatingmelt_input);
    31333183        Input* gllevelset_input = this->GetInput(MaskGroundediceLevelsetEnum); _assert_(gllevelset_input);
     3184        Input* scalefactor_input = NULL;
     3185        if(scaled==true){
     3186                scalefactor_input = this->GetInput(MeshScaleFactorEnum); _assert_(scalefactor_input);
     3187        }
    31343188        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    31353189
    … …  
    31423196                this->JacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
    31433197                floatingmelt_input->GetInputValue(&floatingmelt,gauss);
    3144                 fbmb+=floatingmelt*Jdet*gauss->weight;
     3198                if(scaled==true){
     3199                        scalefactor_input->GetInputValue(&scalefactor,gauss);
     3200                }
     3201                else scalefactor=1;
     3202                fbmb+=floatingmelt*Jdet*gauss->weight*scalefactor;
    31453203        }
    31463204
    … …  
    31583216        bool       mainlyfloating;
    31593217        IssmDouble gbmb=0;
    3160         IssmDouble rho_ice,fraction1,fraction2,groundedmelt,Jdet;
     3218        IssmDouble rho_ice,fraction1,fraction2,groundedmelt,Jdet,scalefactor;
    31613219        IssmDouble Total_Gbmb=0;
    31623220        IssmDouble xyz_list[NUMVERTICES][3];
    … …  
    31693227        Input* groundedmelt_input = this->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(groundedmelt_input);
    31703228        Input* gllevelset_input = this->GetInput(MaskGroundediceLevelsetEnum); _assert_(gllevelset_input);
     3229        Input* scalefactor_input = NULL;
     3230        if(scaled==true){
     3231                scalefactor_input = this->GetInput(MeshScaleFactorEnum); _assert_(scalefactor_input);
     3232        }
    31713233        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    31723234
    … …  
    31793241                this->JacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
    31803242                groundedmelt_input->GetInputValue(&groundedmelt,gauss);
    3181                 gbmb+=groundedmelt*Jdet*gauss->weight;
     3243                if(scaled==true){
     3244                        scalefactor_input->GetInputValue(&scalefactor,gauss);
     3245                }
     3246                else scalefactor=1;
     3247                gbmb+=groundedmelt*Jdet*gauss->weight*scalefactor;
    31823248        }
    31833249
    … …  
    31923258
    31933259        /*The smb[kg yr-1] of one element is area[m2] * smb [kg m^-2 yr^-1]*/
    3194         IssmDouble base,smb,rho_ice;
     3260        IssmDouble base,smb,rho_ice,scalefactor;
    31953261        IssmDouble Total_Smb=0;
    31963262        IssmDouble xyz_list[NUMVERTICES][3];
    … …  
    32103276        /*Now get the average SMB over the element*/
    32113277        Input* smb_input = inputs->GetInput(SmbMassBalanceEnum); _assert_(smb_input);
    3212         smb_input->GetInputAverage(&smb);                                                                                                                                                                                               // average smb on element in m ice s-1
    3213    Total_Smb=rho_ice*base*smb;                                                                                                                                                                                                                  // smb on element in kg s-1
     3278        smb_input->GetInputAverage(&smb);       // average smb on element in m ice s-1
     3279        if(scaled==true){
     3280                Input* scalefactor_input = inputs->GetInput(MeshScaleFactorEnum); _assert_(scalefactor_input);
     3281                scalefactor_input->GetInputAverage(&scalefactor);// average scalefactor on element
     3282        }
     3283        else{
     3284                scalefactor=1.;
     3285        }
     3286   Total_Smb=rho_ice*base*smb*scalefactor;      // smb on element in kg s-1
    32143287
    32153288        /*Return: */

  • issm/trunk-jpl/src/m/solve/parseresultsfrombuffer.js

    r21674 r22330  
    9999                else if (fieldname == 'BasalforcingsFloatingiceMeltingRate') for (var i=0;i<field.length;i++)field[i]= field[i]*yts;
    100100                else if (fieldname == 'TotalSmb') for (var i=0;i<field.length;i++)field[i]= field[i]/Math.pow(10,12)*yts; //(GigaTon/year)
     101                else if (fieldname == 'TotalSmbScaled') for (var i=0;i<field.length;i++)field[i]= field[i]/Math.pow(10,12)*yts; //(GigaTon/year)
    101102                else if (fieldname == 'SmbMassBalance') for (var i=0;i<field.length;i++)field[i]= field[i]*yts;
    102103                else if (fieldname == 'SmbPrecipitation') for (var i=0;i<field.length;i++)field[i]= field[i]*yts;

  • issm/trunk-jpl/src/m/solve/parseresultsfromdisk.m

    r21674 r22330  
    183183        elseif strcmp(fieldname,'TotalFloatingBmb'),
    184184                field = field/10.^12*yts; %(GigaTon/year)
     185        elseif strcmp(fieldname,'TotalFloatingBmbScaled'),
     186                field = field/10.^12*yts; %(GigaTon/year)
    185187        elseif strcmp(fieldname,'TotalGroundedBmb'),
    186188                field = field/10.^12*yts; %(GigaTon/year)
     189        elseif strcmp(fieldname,'TotalGroundedBmbScaled'),
     190                field = field/10.^12*yts; %(GigaTon/year)
    187191        elseif strcmp(fieldname,'TotalSmb'),
     192                field = field/10.^12*yts; %(GigaTon/year)
     193        elseif strcmp(fieldname,'TotalSmbScaled'),
    188194                field = field/10.^12*yts; %(GigaTon/year)
    189195        elseif strcmp(fieldname,'SmbMassBalance'),

  • issm/trunk-jpl/src/m/solve/parseresultsfromdisk.py

    r22267 r22330  
    193193                elif fieldname=='TotalFloatingBmb':
    194194                        field = field/10.**12*yts #(GigaTon/year)
     195                elif fieldname=='TotalFloatingBmbScaled':
     196                        field = field/10.**12*yts #(GigaTon/year)
    195197                elif fieldname=='TotalGroundedBmb':
    196198                        field = field/10.**12*yts #(GigaTon/year)
     199                elif fieldname=='TotalGroundedBmbScaled':
     200                        field = field/10.**12*yts #(GigaTon/year)
    197201                elif fieldname=='TotalSmb':
     202                        field = field/10.**12*yts #(GigaTon/year)
     203                elif fieldname=='TotalSmbScaled':
    198204                        field = field/10.**12*yts #(GigaTon/year)
    199205                elif fieldname=='SmbMassBalance':
Note: See TracChangeset for help on using the changeset viewer.