Changeset 22539


Ignore:
Timestamp:
03/14/18 18:46:50 (7 years ago)
Author:
schlegel
Message:

CHG: edits to the GEMB thermo module, add FAC calculation, fix enums to remove underscore

Location:
issm/trunk-jpl
Files:
11 edited

Legend:

Unmodified
Added
Removed

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

    r22515 r22539  
    28162816        IssmDouble sumMassAdd=0.0;
    28172817        IssmDouble sumdz_add=0.0;
     2818        IssmDouble fac=0.0;
    28182819        IssmDouble sumMass=0.0;
    28192820        IssmDouble dMass=0.0;
    … …  
    31173118                sumP = P +  sumP;
    31183119                sumEC = sumEC + EC;  // evap (-)/cond(+)
    3119                 sumdz_add=dz_add+sumdz_add;
    31203120
    31213121                /*Calculate total system mass:*/
    3122                 sumMass=0; for(int i=0;i<m;i++) sumMass += dz[i]*d[i];
     3122                sumMass=0;
     3123                fac=0;
     3124                for(int i=0;i<m;i++){
     3125                        sumMass += dz[i]*d[i];
     3126                        fac += dz[i]*(rho_ice - fmin(d[i],rho_ice));
     3127                }
    31233128
    31243129                #ifndef _HAVE_ADOLC_ //we want to avoid the round operation at all cost. Not differentiable.
    … …  
    31543159        this->AddInput(new DoubleInput(SmbRunoffEnum,sumR/dt/rho_ice));
    31553160        this->AddInput(new DoubleInput(SmbPrecipitationEnum,sumP/dt/rho_ice));
    3156         this->AddInput(new DoubleInput(SmbDz_addEnum,sumdz_add/yts));
    3157         this->AddInput(new DoubleInput(SmbM_addEnum,sumMassAdd/dt));
     3161        this->AddInput(new DoubleInput(SmbDzAddEnum,sumdz_add));
     3162        this->AddInput(new DoubleInput(SmbMAddEnum,sumMassAdd/dt));
     3163        this->AddInput(new DoubleInput(SmbFACEnum,fac/1000)); // output in meters
    31583164
    31593165        /*Free allocations:{{{*/

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

    r22497 r22539  
    583583        /*outputs:*/
    584584        IssmDouble EC=0.0;
    585         IssmDouble* T=NULL;
     585        IssmDouble* T=*pT;
    586586
    587587        if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("   thermal module\n");
    588 
    589         /*Recover pointers: */
    590         T=*pT;
    591588
    592589        ds = d[0];      // density of top grid cell
    … …  
    686683        }
    687684        dt=max_fdt;
    688        
     685
    689686        // determine mean (harmonic mean) of K/dz for u, d, & p
    690687        Au = xNew<IssmDouble>(m);
    … …  
    783780               
    784781                // calculate heat/wind 'coef_H' stability factor
    785                 if (Ri < -0.03-Ttol) coefH = 1.3 * coefM;
     782                if (Ri <= -0.03+Ttol) coefH = coefM/1.3;
    786783                else coefH = coefM;
    787784               
    … …  
    790787                shf = C * CA * (Ta - Ts);
    791788
    792                 // adjust using Monin–Obukhov stability theory
     789                // adjust using Monin-Obukhov stability theory
    793790                shf = shf / (coefM * coefH);
    794791
    … …  
    796793                // determine if snow pack is melting & calcualte surface vapour pressure over ice or liquid water
    797794                if (Ts >= CtoK-Ttol){
    798                         L = LV;
    799 
     795                        L = LV; //for liquid water at 273.15 k to vapor
     796
     797                        //for liquid surface (assume liquid on surface when Ts == 0 deg C)
     798                        // Wright (1997), US Meteorological Handbook from Murphy and Koop, 2005 Appendix A
     799                        eS = 611.21 * exp(17.502 * (Ts - CtoK) / (240.97 + Ts - CtoK));
     800                }
     801                else{
     802                        L = LS; // latent heat of sublimation
     803                       
    800804                        // for an ice surface Murphy and Koop, 2005 [Equation 7]
    801805                        eS = exp(9.550426 - 5723.265/Ts + 3.53068 * log(Ts) - 0.00728332 * Ts);
    802806                }
    803                 else{
    804                         L = LS; // latent heat of sublimation for liquid surface (assume liquid on surface when Ts == 0 deg C)
    805                         // Wright (1997), US Meteorological Handbook from Murphy and Koop, 2005 Appendix A
    806                         eS = 611.21 * exp(17.502 * (Ts - CtoK) / (240.97 + Ts - CtoK));
    807                 }
    808                
     807
    809808                // Latent heat flux [W m-2]
    810809                lhf = C * L * (eAir - eS) * 0.622 / pAir;
    811                
     810
    812811                // adjust using Monin-Obukhov stability theory (if lhf '+' then there is energy and mass gained at the surface,
    813812                // if '-' then there is mass and energy loss at the surface.
    … …  
    832831               
    833832                // temperature diffusion
    834                 for(int j=0;j<m;j++)T0[1+j]=T[j];
     833                for(int j=0;j<m;j++) T0[1+j]=T[j];
    835834                for(int j=0;j<m;j++) Tu[j] = T0[j];
    836835                for(int j=0;j<m;j++) Td[j] = T0[2+j];
    … …  
    839838                // calculate cumulative evaporation (+)/condensation(-)
    840839                EC = EC + (EC_day/dts)*dt;
    841    
     840
    842841                /* CHECK FOR ENERGY (E) CONSERVATION [UNITS: J]
    843842                //energy flux across lower boundary (energy supplied by underling ice)
    … …  
    854853                */
    855854        }
    856        
     855
    857856        /*Free ressources:*/
    858857        xDelete<IssmDouble>(K);
    … …  
    16801679        // WORKS FINE BUT HAS BEEN DISABLED FOR CONVIENCE OF MODEL OUTPUT
    16811680        // INTERPRETATION
    1682         // // calculate total model depth
     1681        // calculate total model depth
    16831682        Ztot = cellsum(dz,n);
    16841683   
    … …  
    19931992
    19941993        // if V = 0 goes to infinity therfore if V = 0 change
    1995         if(V < 0.01) V=0.01;
     1994        if(V<0.01-Dtol)V=0.01;
    19961995
    19971996        // calculate the Bulk Richardson Number (Ri)
    … …  
    19991998
    20001999        // calculate Monin-Obukhov stability factors 'coef_M' and 'coef_H'
    2001 
     2000       
    20022001        // do not allow Ri to exceed 0.19
    2003         if(Ri>0.19-Ttol)Ri= 0.19;
    2004 
    2005         // calculate momentum 'coef_M' stability factor
    2006         if (Ri > 0.0) coef_M = pow(1.0-5.2*Ri,-1.0); // if stable
    2007         else coef_M = pow(1.0-18*Ri,-0.25);
     2002        Ri = fmin(Ri, 0.19);
     2003
     2004        // calculate momentum 'coefM' stability factor
     2005        if (Ri > 0.0+Ttol){
     2006                // if stable
     2007                coef_M = 1.0/(1.0-5.2*Ri);
     2008        }
     2009        else {
     2010                coef_M =pow (1.0-18*Ri,-0.25);
     2011        }
    20082012
    20092013        // calculate heat/wind 'coef_H' stability factor
    2010         if (Ri < -0.03) coef_H = 1.3 * coef_M;
     2014        if (Ri <= -0.03+Ttol) coef_H = coef_M/1.3;
    20112015        else coef_H = coef_M;
    2012                
     2016
    20132017        //// Bulk-transfer coefficient
    20142018        An =  pow(0.4,2) / pow(log(Tz/z0),2);     // Bulk-transfer coefficient
    … …  
    20222026        shf = shf / (coef_M * coef_H);
    20232027
    2024         //// Latent Heat
    2025         // determine if snow pack is melting & calcualte surface vapour pressure
    2026         // over ice or liquid water
     2028        // Latent Heat
     2029        // determine if snow pack is melting & calcualte surface vapour pressure over ice or liquid water
    20272030        if (Ts >= CtoK-Ttol){
    2028                 L = LV;
    2029 
    2030                 // for an ice surface Murphy and Koop, 2005 [Equation 7]
    2031                 eS = exp(9.550426 - 5723.265/Ts + 3.53068 * log(Ts)- 0.00728332 * Ts);
     2031                L = LV; //for liquid water at 273.15 k to vapor
     2032
     2033                //for liquid surface (assume liquid on surface when Ts == 0 deg C)
     2034                // Wright (1997), US Meteorological Handbook from Murphy and Koop, 2005 Appendix A
     2035                eS = 611.21 * exp(17.502 * (Ts - CtoK) / (240.97 + Ts - CtoK));
    20322036        }
    20332037        else{
    20342038                L = LS; // latent heat of sublimation
    2035                 // for liquid surface (assume liquid on surface when Ts == 0 deg C)
    2036                 // Wright (1997), US Meteorological Handbook from Murphy and Koop,
    2037                 // 2005 Apendix A
    2038                 eS = 611.21 * exp(17.502 * (Ts - CtoK) / (240.97 + Ts - CtoK));
    2039         }
     2039
     2040                // for an ice surface Murphy and Koop, 2005 [Equation 7]
     2041                eS = exp(9.550426 - 5723.265/Ts + 3.53068 * log(Ts) - 0.00728332 * Ts);
     2042        }
     2043
    20402044
    20412045        // Latent heat flux [W m-2]

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

    r22507 r22539  
    484484        SmbIsdensificationEnum,
    485485        SmbIsturbulentfluxEnum,
    486         SmbDz_addEnum,
    487         SmbM_addEnum,
     486        SmbDzAddEnum,
     487        SmbMAddEnum,
     488        SmbFACEnum,
    488489        /*SMBpdd*/
    489490        SMBpddEnum,     

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

    r22507 r22539  
    486486                case SmbIsdensificationEnum : return "SmbIsdensification";
    487487                case SmbIsturbulentfluxEnum : return "SmbIsturbulentflux";
    488                 case SmbDz_addEnum : return "SmbDz_add";
    489                 case SmbM_addEnum : return "SmbM_add";
     488                case SmbDzAddEnum : return "SmbDzAdd";
     489                case SmbMAddEnum : return "SmbMAdd";
     490                case SmbFACEnum : return "SmbFAC";
    490491                case SMBpddEnum : return "SMBpdd";
    491492                case SmbDelta18oEnum : return "SmbDelta18o";

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

    r22507 r22539  
    495495              else if (strcmp(name,"SmbIsdensification")==0) return SmbIsdensificationEnum;
    496496              else if (strcmp(name,"SmbIsturbulentflux")==0) return SmbIsturbulentfluxEnum;
    497               else if (strcmp(name,"SmbDz_add")==0) return SmbDz_addEnum;
    498               else if (strcmp(name,"SmbM_add")==0) return SmbM_addEnum;
     497              else if (strcmp(name,"SmbDzAdd")==0) return SmbDzAddEnum;
     498              else if (strcmp(name,"SmbMAdd")==0) return SmbMAddEnum;
     499              else if (strcmp(name,"SmbFAC")==0) return SmbFACEnum;
    499500              else if (strcmp(name,"SMBpdd")==0) return SMBpddEnum;
    500501              else if (strcmp(name,"SmbDelta18o")==0) return SmbDelta18oEnum;
    … …  
    505506              else if (strcmp(name,"SmbIsd18opd")==0) return SmbIsd18opdEnum;
    506507              else if (strcmp(name,"SmbIstemperaturescaled")==0) return SmbIstemperaturescaledEnum;
    507               else if (strcmp(name,"SmbPrecipitationsPresentday")==0) return SmbPrecipitationsPresentdayEnum;
    508508         else stage=5;
    509509   }
    510510   if(stage==5){
    511               if (strcmp(name,"SmbPrecipitationsLgm")==0) return SmbPrecipitationsLgmEnum;
     511              if (strcmp(name,"SmbPrecipitationsPresentday")==0) return SmbPrecipitationsPresentdayEnum;
     512              else if (strcmp(name,"SmbPrecipitationsLgm")==0) return SmbPrecipitationsLgmEnum;
    512513              else if (strcmp(name,"SmbTemperaturesPresentday")==0) return SmbTemperaturesPresentdayEnum;
    513514              else if (strcmp(name,"SmbTemperaturesLgm")==0) return SmbTemperaturesLgmEnum;
    … …  
    628629              else if (strcmp(name,"DeviatoricStressxz")==0) return DeviatoricStressxzEnum;
    629630              else if (strcmp(name,"DeviatoricStressyy")==0) return DeviatoricStressyyEnum;
    630               else if (strcmp(name,"DeviatoricStressyz")==0) return DeviatoricStressyzEnum;
    631631         else stage=6;
    632632   }
    633633   if(stage==6){
    634               if (strcmp(name,"DeviatoricStresszz")==0) return DeviatoricStresszzEnum;
     634              if (strcmp(name,"DeviatoricStressyz")==0) return DeviatoricStressyzEnum;
     635              else if (strcmp(name,"DeviatoricStresszz")==0) return DeviatoricStresszzEnum;
    635636              else if (strcmp(name,"DeviatoricStresseffective")==0) return DeviatoricStresseffectiveEnum;
    636637              else if (strcmp(name,"LambdaS")==0) return LambdaSEnum;
    … …  
    751752              else if (strcmp(name,"Outputdefinition89")==0) return Outputdefinition89Enum;
    752753              else if (strcmp(name,"Outputdefinition90")==0) return Outputdefinition90Enum;
    753               else if (strcmp(name,"Outputdefinition91")==0) return Outputdefinition91Enum;
    754754         else stage=7;
    755755   }
    756756   if(stage==7){
    757               if (strcmp(name,"Outputdefinition92")==0) return Outputdefinition92Enum;
     757              if (strcmp(name,"Outputdefinition91")==0) return Outputdefinition91Enum;
     758              else if (strcmp(name,"Outputdefinition92")==0) return Outputdefinition92Enum;
    758759              else if (strcmp(name,"Outputdefinition93")==0) return Outputdefinition93Enum;
    759760              else if (strcmp(name,"Outputdefinition94")==0) return Outputdefinition94Enum;
    … …  
    874875              else if (strcmp(name,"SealevelriseRequestedOutputs")==0) return SealevelriseRequestedOutputsEnum;
    875876              else if (strcmp(name,"SealevelriseNumRequestedOutputs")==0) return SealevelriseNumRequestedOutputsEnum;
    876               else if (strcmp(name,"LoveNfreq")==0) return LoveNfreqEnum;
    877877         else stage=8;
    878878   }
    879879   if(stage==8){
    880               if (strcmp(name,"LoveFrequencies")==0) return LoveFrequenciesEnum;
     880              if (strcmp(name,"LoveNfreq")==0) return LoveNfreqEnum;
     881              else if (strcmp(name,"LoveFrequencies")==0) return LoveFrequenciesEnum;
    881882              else if (strcmp(name,"LoveShNmax")==0) return LoveShNmaxEnum;
    882883              else if (strcmp(name,"LoveShNmin")==0) return LoveShNminEnum;
    … …  
    997998              else if (strcmp(name,"Tria")==0) return TriaEnum;
    998999              else if (strcmp(name,"TriaInput")==0) return TriaInputEnum;
    999               else if (strcmp(name,"Tetra")==0) return TetraEnum;
    10001000         else stage=9;
    10011001   }
    10021002   if(stage==9){
    1003               if (strcmp(name,"TetraInput")==0) return TetraInputEnum;
     1003              if (strcmp(name,"Tetra")==0) return TetraEnum;
     1004              else if (strcmp(name,"TetraInput")==0) return TetraInputEnum;
    10041005              else if (strcmp(name,"Penta")==0) return PentaEnum;
    10051006              else if (strcmp(name,"PentaInput")==0) return PentaInputEnum;
    … …  
    11201121              else if (strcmp(name,"P1xP3")==0) return P1xP3Enum;
    11211122              else if (strcmp(name,"P1xP4")==0) return P1xP4Enum;
    1122               else if (strcmp(name,"P2xP4")==0) return P2xP4Enum;
    11231123         else stage=10;
    11241124   }
    11251125   if(stage==10){
    1126               if (strcmp(name,"P1P1")==0) return P1P1Enum;
     1126              if (strcmp(name,"P2xP4")==0) return P2xP4Enum;
     1127              else if (strcmp(name,"P1P1")==0) return P1P1Enum;
    11271128              else if (strcmp(name,"P1P1GLS")==0) return P1P1GLSEnum;
    11281129              else if (strcmp(name,"MINI")==0) return MINIEnum;

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

    r22383 r22539  
    199199        elseif strcmp(fieldname,'SmbRunoff'),
    200200                field = field*yts;
     201        elseif strcmp(fieldname,'SmbEvaporation'),
     202                field = field*yts;
     203        elseif strcmp(fieldname,'SmbRefreeze'),
     204                field = field*yts;
    201205        elseif strcmp(fieldname,'SmbEC'),
    202206                field = field*yts;
    … …  
    205209        elseif strcmp(fieldname,'SmbMelt'),
    206210                field = field*yts;
    207     elseif strcmp(fieldname,'SmbDz_add'),
    208         field = field*yts;
    209     elseif strcmp(fieldname,'SmbM_add'),
    210         field = field*yts;
     211        elseif strcmp(fieldname,'SmbMAdd'),
     212                field = field*yts;
    211213        elseif strcmp(fieldname,'CalvingCalvingrate'),
    212214                field = field*yts;

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

    r22399 r22539  
    209209                elif fieldname=='SmbRunoff':
    210210                        field = field*yts
     211                elif fieldname=='SmbEvaporation':
     212                        field = field*yts;
     213                elif fieldname=='SmbRefreeze':
     214                        field = field*yts;
    211215                elif fieldname=='SmbEC':
    212216                        field = field*yts
    … …  
    215219                elif fieldname=='SmbMelt':
    216220                        field = field*yts
    217                 elif fieldname=='SmbDz_add':
    218                         field = field*yts
    219                 elif fieldname=='SmbM_add':
     221                elif fieldname=='SmbMAdd':
    220222                        field = field*yts
    221223                elif fieldname=='CalvingCalvingrate':

  • issm/trunk-jpl/test/NightlyRun/test243.m

    r22432 r22539  
    2727
    2828%smb settings
    29 md.smb.requested_outputs={'SmbDz','SmbT','SmbD','SmbRe','SmbGdn','SmbGsp','SmbEC','SmbA','SmbMassBalance'};
     29md.smb.requested_outputs={'SmbDz','SmbT','SmbD','SmbRe','SmbGdn','SmbGsp','SmbEC','SmbA','SmbMassBalance','SmbMAdd','SmbDzAdd','SmbFAC'};
    3030
    3131%only run smb core:
    … …  
    4444
    4545%Fields and tolerances to track changes
    46 field_names      ={'SmbDz','SmbT' ,'SmbD' ,'SmbRe','SmbGdn','SmbGsp','SmbA' ,'SmbEC','SmbMassBalance'};
    47 field_tolerances ={1e-11,1e-12,1e-11,1e-11,1e-11,1e-11,1e-12,1e-12,1e-12};
     46field_names      ={'SmbDz','SmbT' ,'SmbD' ,'SmbRe','SmbGdn','SmbGsp','SmbA' ,'SmbEC','SmbMassBalance','SmbMAdd','SmbDzAdd','SmbFAC'};
     47field_tolerances ={1e-11,1e-12,1e-11,1e-11,1e-11,1e-11,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12};
    4848
    4949field_values={...
    … …  
    5757        (md.results.TransientSolution(end).SmbEC(1)),...
    5858        (md.results.TransientSolution(end).SmbMassBalance(1)),...
     59        (md.results.TransientSolution(end).SmbMAdd(1)),...
     60        (md.results.TransientSolution(end).SmbDzAdd(1)),...
     61        (md.results.TransientSolution(end).SmbFAC(1))
    5962        };

  • issm/trunk-jpl/test/NightlyRun/test243.py

    r22432 r22539  
    4040
    4141#smb settings
    42 md.smb.requested_outputs = ['SmbDz','SmbT','SmbD','SmbRe','SmbGdn','SmbGsp','SmbEC','SmbA','SmbMassBalance']
     42md.smb.requested_outputs = ['SmbDz','SmbT','SmbD','SmbRe','SmbGdn','SmbGsp','SmbEC','SmbA','SmbMassBalance','SmbMAdd','SmbDzAdd','SmbFAC']
    4343
    4444#only run smb core:
    … …  
    5757
    5858#Fields and tolerances to track changes
    59 field_names      = ['SmbDz','SmbT' ,'SmbD' ,'SmbRe','SmbGdn','SmbGsp','SmbA' ,'SmbEC','SmbMassBalance']
    60 field_tolerances = [1e-11,1e-12,1e-11,1e-11,1e-11,1e-11,1e-12,1e-12,1e-12]
     59field_names      = ['SmbDz','SmbT' ,'SmbD' ,'SmbRe','SmbGdn','SmbGsp','SmbA' ,'SmbEC','SmbMassBalance','SmbMAdd','SmbDzAdd','SmbFAC']
     60field_tolerances = [1e-11,1e-12,1e-11,1e-11,1e-11,1e-11,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12]
    6161#shape is different in python solution (fixed using reshape) which can cause test failure:
    6262field_values = [
    … …  
    7070        md.results.TransientSolution[-1].SmbEC[0],
    7171        md.results.TransientSolution[-1].SmbMassBalance[0],
     72        md.results.TransientSolution[-1].SmbMAdd[0],
     73        md.results.TransientSolution[-1].SmbDzAdd[0],
     74        md.results.TransientSolution[-1].SmbFAC[0],
    7275        ]
Note: See TracChangeset for help on using the changeset viewer.