Changeset 25997


Ignore:
Timestamp:
02/17/21 18:10:09 (4 years ago)
Author:
schlegel
Message:

CHG: implement inputs for gardner albedo scheme

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

Legend:

Unmodified
Added
Removed

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

    r25539 r25997  
    4646                        iomodel->FetchDataToInput(inputs,elements,"md.smb.V",SmbVEnum);
    4747                        iomodel->FetchDataToInput(inputs,elements,"md.smb.dswrf",SmbDswrfEnum);
     48                        iomodel->FetchDataToInput(inputs,elements,"md.smb.dswdiffrf",SmbDswdiffrfEnum);
    4849                        iomodel->FetchDataToInput(inputs,elements,"md.smb.dlwrf",SmbDlwrfEnum);
    4950                        iomodel->FetchDataToInput(inputs,elements,"md.smb.P",SmbPEnum);
    … …  
    7071                        iomodel->FetchDataToInput(inputs,elements,"md.smb.Wini",SmbWiniEnum);
    7172                        iomodel->FetchDataToInput(inputs,elements,"md.smb.Aini",SmbAiniEnum);
     73                        iomodel->FetchDataToInput(inputs,elements,"md.smb.Adiffini",SmbAdiffiniEnum);
    7274                        iomodel->FetchDataToInput(inputs,elements,"md.smb.Tini",SmbTiniEnum);
    7375                        iomodel->FetchDataToInput(inputs,elements,"md.smb.Sizeini",SmbSizeiniEnum);
    7476                        iomodel->FetchDataToInput(inputs,elements,"md.smb.aValue",SmbAValueEnum);
    7577                        iomodel->FetchDataToInput(inputs,elements,"md.smb.teValue",SmbTeValueEnum);
     78                        iomodel->FetchDataToInput(inputs,elements,"md.smb.szaValue",SmbSzaValueEnum);
     79                        iomodel->FetchDataToInput(inputs,elements,"md.smb.cotValue",SmbCotValueEnum);
     80                        iomodel->FetchDataToInput(inputs,elements,"md.smb.ccsnowValue",SmbCcsnowValueEnum);
     81                        iomodel->FetchDataToInput(inputs,elements,"md.smb.cciceValue",SmbCciceValueEnum);
    7682                        break;
    7783                case SMBpddEnum:

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

    r25993 r25997  
    35463546        IssmDouble dlw=0.0;
    35473547        IssmDouble dsw=0.0;
     3548        IssmDouble dswdiff=0.0;
    35483549        IssmDouble P=0.0;
    35493550        IssmDouble eAir=0.0;
    … …  
    36023603        IssmDouble* W = NULL;
    36033604        IssmDouble* a = NULL;
     3605        IssmDouble* adiff = NULL;
    36043606        IssmDouble* swf=NULL;
    36053607        IssmDouble* T = NULL;
    … …  
    36173619        IssmDouble* Wini = NULL;
    36183620        IssmDouble* aini = NULL;
     3621        IssmDouble* adiffini = NULL;
    36193622        IssmDouble* Tini = NULL;
    36203623        int         m=0;
    … …  
    36953698                this->inputs->GetArray(SmbWiniEnum,this->lid,&Wini,&m);
    36963699                this->inputs->GetArray(SmbAiniEnum,this->lid,&aini,&m);
     3700                this->inputs->GetArray(SmbAdiffiniEnum,this->lid,&adiffini,&m);
    36973701                this->inputs->GetArray(SmbTiniEnum,this->lid,&Tini,&m);
    36983702                EC_input = this->GetInput(SmbECiniEnum);  _assert_(EC_input);
    … …  
    37143718                        W = xNew<IssmDouble>(m); for(int i=0;i<m;i++)W[i]=Wini[0];             //set water content to zero [kg m-2]
    37153719                        a = xNew<IssmDouble>(m); for(int i=0;i<m;i++)a[i]=aini[0];         //set albedo equal to fresh snow [fraction]
     3720                        adiff = xNew<IssmDouble>(m); for(int i=0;i<m;i++)adiff[i]=adiffini[0];         //set diffusive albedo equal to 1 [fraction]
    37163721                        T = xNew<IssmDouble>(m); for(int i=0;i<m;i++)T[i]=Tmean;         //set initial grid cell temperature to the annual mean temperature [K]
    37173722                        /*/!\ Default value of T can not be retrived from SMBgemb.m (like other snow properties)
    … …  
    37323737                        W = xNew<IssmDouble>(m);for(int i=0;i<m;i++)W[i]=Wini[i];
    37333738                        a = xNew<IssmDouble>(m);for(int i=0;i<m;i++)a[i]=aini[i];
     3739                        adiff = xNew<IssmDouble>(m);for(int i=0;i<m;i++)adiff[i]=adiffini[i];
    37343740                        T = xNew<IssmDouble>(m);for(int i=0;i<m;i++)T[i]=Tini[i];
    37353741
    … …  
    37573763                this->inputs->GetArray(SmbWEnum,this->lid,&W,&m);
    37583764                this->inputs->GetArray(SmbAEnum,this->lid,&a,&m);
     3765                this->inputs->GetArray(SmbAdiffEnum,this->lid,&adiff,&m);
    37593766                this->inputs->GetArray(SmbTEnum,this->lid,&T,&m);
    37603767                EC_input = this->GetInput(SmbECDtEnum);  _assert_(EC_input);
    … …  
    38213828        Input *Dlwr_input= this->GetInput(SmbDlwrfEnum,timeinputs); _assert_(Dlwr_input);
    38223829        Input *Dswr_input= this->GetInput(SmbDswrfEnum,timeinputs); _assert_(Dswr_input);
     3830        Input *Dswrdiff_input= this->GetInput(SmbDswdiffrfEnum,timeinputs); _assert_(Dswrdiff_input);
    38233831        Input *P_input   = this->GetInput(SmbPEnum,timeinputs);     _assert_(P_input);
    38243832        Input *eAir_input= this->GetInput(SmbEAirEnum,timeinputs);  _assert_(eAir_input);
    … …  
    38263834        Input *teValue_input= this->GetInput(SmbTeValueEnum,timeinputs); _assert_(teValue_input);
    38273835        Input *aValue_input= this->GetInput(SmbAValueEnum,timeinputs); _assert_(aValue_input);
     3836        Input *szaValue_input= this->GetInput(SmbSzaValueEnum,timeinputs); _assert_(szaValue_input);
     3837        Input *cotValue_input= this->GetInput(SmbCotValueEnum,timeinputs); _assert_(cotValue_input);
     3838        Input *ccsnowValue_input= this->GetInput(SmbCcsnowValueEnum,timeinputs); _assert_(ccsnowValue_input);
     3839        Input *cciceValue_input= this->GetInput(SmbCciceValueEnum,timeinputs); _assert_(cciceValue_input);
    38283840
    38293841        /*extract daily data:{{{*/
    … …  
    38323844        Dlwr_input->GetInputValue(&dlw,gauss);   //downward longwave radiation flux [W m-2]
    38333845        Dswr_input->GetInputValue(&dsw,gauss);   //downward shortwave radiation flux [W m-2]
     3846        Dswrdiff_input->GetInputValue(&dswdiff,gauss);   //downward shortwave diffuse radiation flux [W m-2]
    38343847        P_input->GetInputValue(&P,gauss);        //precipitation [kg m-2]
    38353848        eAir_input->GetInputValue(&eAir,gauss);  //screen level vapor pressure [Pa]
    … …  
    38373850        teValue_input->GetInputValue(&teValue,gauss);  // Emissivity [0-1]
    38383851        aValue_input->GetInputValue(&aValue,gauss);  // Albedo [0 1]
    3839         //Not implemented as input
    3840         //szaValue;  // Solar Zenith Angle [degree]
    3841         //cotValue;  // Cloud Optical Thickness
    3842         //ccsnowValue;//concentration of light absorbing carbon for snow [ppm1]
    3843         //cciceValue;//concentration of light absorbing carbon for ice [ppm1]
     3852        szaValue_input->GetInputValue(&szaValue,gauss);  // Solar Zenith Angle [degree]
     3853        cotValue_input->GetInputValue(&cotValue,gauss);  // Cloud Optical Thickness
     3854        ccsnowValue_input->GetInputValue(&ccsnowValue,gauss); //concentration of light absorbing carbon for snow [ppm1]
     3855        cciceValue_input->GetInputValue(&cciceValue,gauss); //concentration of light absorbing carbon for ice [ppm1]
    38443856        //_printf_("Time: " << t << " Ta: " << Ta << " V: " << V << " dlw: " << dlw << " dsw: " << dsw << " P: " << P << " eAir: " << eAir << " pAir: " << pAir << "\n");
    38453857        /*}}}*/
    … …  
    38493861
    38503862        /*Snow, firn and ice albedo:*/
    3851         if(isalbedo)albedo(&a,aIdx,re,dz,d,cldFrac,aIce,aSnow,aValue,adThresh,T,W,P,EC,Msurf,ccsnowValue,cciceValue,szaValue,cotValue,t0wet,t0dry,K,smb_dt,rho_ice,m,this->Sid());
     3863        if(isalbedo)albedo(&a,&adiff,aIdx,re,dz,d,cldFrac,aIce,aSnow,aValue,adThresh,T,W,P,EC,Msurf,ccsnowValue,cciceValue,szaValue,cotValue,t0wet,t0dry,K,smb_dt,rho_ice,m,this->Sid());
    38523864
    38533865        /*Distribution of absorbed short wave radation with depth:*/
    3854         if(isshortwave)shortwave(&swf, swIdx, aIdx, dsw, a[0], d, dz, re,rho_ice,m,this->Sid());
     3866        if(isshortwave)shortwave(&swf, swIdx, aIdx, dsw, dswdiff, a[0], adiff[0], d, dz, re,rho_ice,m,this->Sid());
    38553867
    38563868        /*Calculate net shortwave [W m-2]*/
    … …  
    38693881
    38703882        /*Add snow/rain to top grid cell adjusting cell depth, temperature and density*/
    3871         if(isaccumulation)accumulation(&T, &dz, &d, &W, &a, &re, &gdn, &gsp, &m, aIdx, dsnowIdx, Tmean, Ta, P, dzMin, aSnow, C, V, Vmean, rho_ice,this->Sid());
     3883        if(isaccumulation)accumulation(&T, &dz, &d, &W, &a, &adiff, &re, &gdn, &gsp, &m, aIdx, dsnowIdx, Tmean, Ta, P, dzMin, aSnow, C, V, Vmean, rho_ice,this->Sid());
    38723884
    38733885        /*Calculate water production, M [kg m-2] resulting from snow/ice temperature exceeding 273.15 deg K
    38743886         * (> 0 deg C), runoff R [kg m-2] and resulting changes in density and determine wet compaction [m]*/
    3875         if(ismelt)melt(&M, &Msurf, &R, &mAdd, &dz_add, &T, &d, &dz, &W, &a, &re, &gdn, &gsp, &m, dzMin, zMax, zMin, zTop, zY, rho_ice,this->Sid());
     3887        if(ismelt)melt(&M, &Msurf, &R, &mAdd, &dz_add, &T, &d, &dz, &W, &a, &adiff, &re, &gdn, &gsp, &m, dzMin, zMax, zMin, zTop, zY, rho_ice,this->Sid());
    38763888
    38773889        /*Allow non-melt densification and determine compaction [m]*/
    … …  
    39033915                                        << "lwf[" << netLW << "] "
    39043916                                        << "a[" << a << "] "
     3917                                        << "adiff[" << adiff << "] "
    39053918                                        << "te[" << teValue << "] "
    39063919                                        << "\n");
    … …  
    39753988        this->inputs->SetArrayInput(SmbWEnum,this->lid,W,m);
    39763989        this->inputs->SetArrayInput(SmbAEnum,this->lid,a,m);
     3990        this->inputs->SetArrayInput(SmbAdiffEnum,this->lid,adiff,m);
    39773991        this->SetElementInput(SmbECEnum,sumEC/dt/rho_ice);
    39783992        this->SetElementInput(SmbMassBalanceEnum,(sumP + sumEC -sumR)/dt/rho_ice);
    … …  
    40014015        if(W) xDelete<IssmDouble>(W);
    40024016        if(a) xDelete<IssmDouble>(a);
     4017        if(adiff) xDelete<IssmDouble>(adiff);
    40034018        if(T) xDelete<IssmDouble>(T);
    40044019        if(dzini) xDelete<IssmDouble>(dzini);
    … …  
    40094024        if(Wini) xDelete<IssmDouble>(Wini);
    40104025        if(aini) xDelete<IssmDouble>(aini);
     4026        if(adiffini) xDelete<IssmDouble>(adiffini);
    40114027        if(Tini) xDelete<IssmDouble>(Tini);
    40124028        if(swf) xDelete<IssmDouble>(swf);

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

    r25993 r25997  
    197197
    198198} /*}}}*/
     199IssmDouble gardnerAlb(IssmDouble* re, IssmDouble* dz, IssmDouble* d, IssmDouble clabSnow, IssmDouble clabIce, IssmDouble SZA, IssmDouble COT, IssmDouble dPHC, int m){ /*{{{*/
     200        //gardnerAlb(S1, c1, SZA, t, z1, S2, c2)
     201        //This is an implementation of the snow and ice broadband albedo
     202        //  parameterization developed by Alex Gardner.
     203        //Created By: Alex S. Gardner, Jet Propulsion Laboratory [alex.s.gardner@jpl.nasa.gov]
     204        //  Last Modified: June, 2014
     205        //Full Reference: Gardner, A. S., and Sharp, M. J.: A review of snow and
     206        //  ice albedo and the development of a new physically based broadband albedo
     207        //  parameterization, J. Geophys. Res., 115, F01009, 10.1029/2009jf001444,
     208        //  2010.
     209
     210        //INPUTS
     211        // ONE LAYER
     212        //  - S1    : specific surface area of the snow or ice [cm^2 g-1]
     213        //  - c1    : concentration of light absorbing carbon  [ppm1]
     214        //  - SZA   : solar zenith angle of the incident radiation [deg]
     215        //  - t     : cloud optical thickness
     216        // TWO LAYER
     217        //  - z1    : depth of snow suface layer [mm w.e.]
     218        //  - S2    : specific surface area of bottom ice layer [cm^2 g-1]
     219        //  - c2    : concentration of light absorbing carbon of bottom ice
     220        //             layer [ppm1]
     221        IssmDouble c1=clabSnow;
     222        IssmDouble c2=clabIce;
     223        IssmDouble t=COT;
     224        IssmDouble a=0.0;
     225
     226        //Single layer albedo parameterization
     227        //convert effective radius to specific surface area [cm2 g-1]
     228        IssmDouble S1 = 3.0 / (0.091 * re[0]);
     229
     230        //effective solar zenith angle
     231        IssmDouble x = min(pow(t/(3*cos(Pi*SZA/180.0)),0.5), 1.0);
     232        IssmDouble u = 0.64*x + (1-x)*cos(Pi*SZA/180.0);
     233
     234        // pure snow albedo
     235        IssmDouble as = 1.48 - pow(S1,-0.07);
     236
     237        //change in pure snow albedo due to soot loading
     238        IssmDouble dac = max(0.04 - as, pow(-c1,0.55)/(0.16 + 0.6*pow(S1,0.5) + (1.8*pow(c1,0.6))*pow(S1,-0.25)));
     239
     240        //Two layer albedo parameterization
     241        //  do two layer calculation if there is more than 1 layer
     242        IssmDouble z1=0.0;
     243        int lice=0;
     244        for(int l=0;(l<m & d[l]<dPHC-Dtol);l++){
     245                z1=z1+dz[l]*d[l]; //mm
     246                lice=l+1;
     247        }
     248        if (m>0 & lice<m & z1 > Dtol){
     249                // determine albedo values for bottom layer
     250                IssmDouble S2 = 3.0 / (0.091 * re[lice]);
     251
     252                // pure snow albedo
     253                IssmDouble as2 = 1.48 - pow(S2,-0.07);
     254
     255                // change in pure snow albedo due to soot loading
     256                IssmDouble dac2 = max(0.04 - as2, pow(-c2,0.55)/(0.16 + 0.6*pow(S1,0.5) + (1.8*pow(c2,0.6))*pow(S2,-0.25)));
     257
     258                // determine the effective change due to finite depth and soot loading
     259                IssmDouble A = min(1.0, (2.1 * pow(z1,1.35*(1-as) - 0.1*c1 - 0.13)));
     260
     261                dac =  (as2 + dac2 - as) + A*((as + dac) - (as2 + dac2));
     262        }
     263
     264        // change in albedo due to solar zenith angle
     265        IssmDouble dasz = 0.53*as*(1 - (as + dac))*pow(1 - u,1.2);
     266
     267        // change in albedo due to cloud (apart from change in diffuse fraction)
     268        IssmDouble dat = (0.1*t*pow(as + dac,1.3)) / (pow(1 + 1.5*t,as));
     269
     270        // Broadband albedo
     271        a = as + dac + dasz + dat;
     272
     273        return a;
     274}   /*}}}*/
    199275void grainGrowth(IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, IssmDouble* T,IssmDouble* dz,IssmDouble* d, IssmDouble* W,IssmDouble smb_dt,int m,int aIdx,int sid){ /*{{{*/
    200276
    … …  
    410486
    411487}  /*}}}*/
    412 void albedo(IssmDouble** pa, int aIdx, IssmDouble* re, IssmDouble* dz, IssmDouble* d, IssmDouble cldFrac, IssmDouble aIce, IssmDouble aSnow, IssmDouble aValue, IssmDouble adThresh, IssmDouble* TK, IssmDouble* W, IssmDouble P, IssmDouble EC, IssmDouble Msurf, IssmDouble clabSnow, IssmDouble clabIce, IssmDouble SZA, IssmDouble COT, IssmDouble t0wet, IssmDouble t0dry, IssmDouble K, IssmDouble dt, IssmDouble dIce, int m,int sid) { /*{{{*/
     488void albedo(IssmDouble** pa, IssmDouble** padiff, int aIdx, IssmDouble* re, IssmDouble* dz, IssmDouble* d, IssmDouble cldFrac, IssmDouble aIce, IssmDouble aSnow, IssmDouble aValue, IssmDouble adThresh, IssmDouble* TK, IssmDouble* W, IssmDouble P, IssmDouble EC, IssmDouble Msurf, IssmDouble clabSnow, IssmDouble clabIce, IssmDouble SZA, IssmDouble COT, IssmDouble t0wet, IssmDouble t0dry, IssmDouble K, IssmDouble dt, IssmDouble dIce, int m,int sid) { /*{{{*/
    413489
    414490        //// Calculates Snow, firn and ice albedo as a function of:
    … …  
    471547        /*output: */
    472548        IssmDouble* a=NULL;
     549        IssmDouble* adiff=NULL;
    473550
    474551        if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("   albedo module\n");
    … …  
    476553        /*Recover pointers: */
    477554        a=*pa;
     555        adiff=*padiff;
    478556
    479557        //some constants:
    … …  
    490568                if(aIdx==1){
    491569                        //function of effective grain radius
    492 
    493                         //gardnerAlb(S1, c1, SZA, t, z1, S2, c2)
    494                         //This is an implementation of the snow and ice broadband albedo
    495                         //  parameterization developed by Alex Gardner.
    496                         //Created By: Alex S. Gardner, Jet Propulsion Laboratory [alex.s.gardner@jpl.nasa.gov]
    497                         //  Last Modified: June, 2014
    498                         //Full Reference: Gardner, A. S., and Sharp, M. J.: A review of snow and
    499                         //  ice albedo and the development of a new physically based broadband albedo
    500                         //  parameterization, J. Geophys. Res., 115, F01009, 10.1029/2009jf001444,
    501                         //  2010.
    502 
    503                         //INPUTS
    504                         // ONE LAYER
    505                         //  - S1    : specific surface area of the snow or ice [cm^2 g-1]
    506                         //  - c1    : concentration of light absorbing carbon  [ppm1]
    507                         //  - SZA   : solar zenith angle of the incident radiation [deg]
    508                         //  - t     : cloud optical thickness
    509                         // TWO LAYER
    510                         //  - z1    : depth of snow suface layer [mm w.e.]
    511                         //  - S2    : specific surface area of bottom ice layer [cm^2 g-1]
    512                         //  - c2    : concentration of light absorbing carbon of bottom ice
    513                         //             layer [ppm1]
    514                         IssmDouble c1=clabSnow;
    515                         IssmDouble c2=clabIce;
    516                         IssmDouble t=COT;
    517 
    518                         //Single layer albedo parameterization
    519                         //convert effective radius to specific surface area [cm2 g-1]
    520                         IssmDouble S1 = 3.0 / (0.091 * re[0]);
    521 
    522                         //effective solar zenith angle
    523                         IssmDouble x = min(pow(t/(3*cos(Pi*SZA/180.0)),0.5), 1.0);
    524                         IssmDouble u = 0.64*x + (1-x)*cos(Pi*SZA/180.0);
    525 
    526                         // pure snow albedo
    527                         IssmDouble as = 1.48 - pow(S1,-0.07);
    528 
    529                         //change in pure snow albedo due to soot loading
    530                         IssmDouble dac = max(0.04 - as, pow(-c1,0.55)/(0.16 + 0.6*pow(S1,0.5) + (1.8*pow(c1,0.6))*pow(S1,-0.25)));
    531 
    532                         //Two layer albedo parameterization
    533                         //  do two layer calculation if there is more than 1 layer
    534                         IssmDouble z1=0.0;
    535                         int lice=0;
    536                         for(int l=0;(l<m & d[l]<dPHC-Dtol);l++){
    537                                 z1=z1+dz[l]*d[l]; //mm
    538                                 lice=l+1;
    539                         }
    540                         if (m>0 & lice<m & z1 > Dtol){
    541                                 // determine albedo values for bottom layer
    542                                 IssmDouble S2 = 3.0 / (0.091 * re[lice]);
    543 
    544                                 // pure snow albedo
    545                                 IssmDouble as2 = 1.48 - pow(S2,-0.07);
    546 
    547                                 // change in pure snow albedo due to soot loading
    548                                 IssmDouble dac2 = max(0.04 - as2, pow(-c2,0.55)/(0.16 + 0.6*pow(S1,0.5) + (1.8*pow(c2,0.6))*pow(S2,-0.25)));
    549 
    550                                 // determine the effective change due to finite depth and soot loading
    551                                 IssmDouble A = min(1.0, (2.1 * pow(z1,1.35*(1-as) - 0.1*c1 - 0.13)));
    552 
    553                                 dac =  (as2 + dac2 - as) + A*((as + dac) - (as2 + dac2));
    554                         }
    555 
    556                         // change in albedo due to solar zenith angle
    557                         IssmDouble dasz = 0.53*as*(1 - (as + dac))*pow(1 - u,1.2);
    558 
    559                         // change in albedo due to cloud (apart from change in diffuse fraction)
    560                         IssmDouble dat = (0.1*t*pow(as + dac,1.3)) / (pow(1 + 1.5*t,as));
    561 
    562                         // Broadband albedo
    563                         a[0] = as + dac + dasz + dat;
     570                        // clabSnow, IssmDouble clabIce, IssmDouble SZA, IssmDouble COT, int m
     571                        a[0]=gardnerAlb(re, dz, d, clabSnow, clabIce, SZA, COT, dPHC, m);
     572                        adiff[0]=gardnerAlb(re, dz, d, clabSnow, clabIce, 50, COT, dPHC, m);
    564573                }
    565574                else if(aIdx==2){
    … …  
    698707        /*Assign output pointers:*/
    699708        *pa=a;
     709        *padiff=adiff;
    700710
    701711}  /*}}}*/
    … …  
    10781088
    10791089}  /*}}}*/
    1080 void shortwave(IssmDouble** pswf, int swIdx, int aIdx, IssmDouble dsw, IssmDouble as, IssmDouble* d, IssmDouble* dz, IssmDouble* re, IssmDouble dIce, int m, int sid){ /*{{{*/
     1090void shortwave(IssmDouble** pswf, int swIdx, int aIdx, IssmDouble dsw, IssmDouble dswdiff, IssmDouble as, IssmDouble asdiff, IssmDouble* d, IssmDouble* dz, IssmDouble* re, IssmDouble dIce, int m, int sid){ /*{{{*/
    10811091
    10821092        // DISTRIBUTES ABSORBED SHORTWAVE RADIATION WITHIN SNOW/ICE
    … …  
    10921102        //   aIdx    = method for calculating albedo (1-4)
    10931103        //   dsw     = downward shortwave radiative flux [w m-2]
     1104        //   dswdir  = downward shortwave direct radiative flux [w m-2]
    10941105        //   as      = surface albedo
     1106        //   asdiff  = surface albedo for diffuse radiation
    10951107        //   d       = grid cell density [kg m-3]
    10961108        //   dz      = grid cell depth [m]
    … …  
    11131125
    11141126                // calculate surface shortwave radiation fluxes [W m-2]
    1115                 swf[0] = (1.0 - as) * dsw;
     1127                if (aIdx == 1){
     1128                        swf[0] = (1.0 - as) * dsw +  (1.0 - asdiff) * dswdiff;
     1129                }
     1130                else{
     1131                        swf[0] = (1.0 - as) * dsw;
     1132                }
    11161133        }
    11171134        else{ // sw radation is absorbed at depth within the glacier
    … …  
    12771294
    12781295} /*}}}*/
    1279 void accumulation(IssmDouble** pT, IssmDouble** pdz, IssmDouble** pd, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pm, int aIdx, int dsnowIdx, IssmDouble Tmean, IssmDouble T_air, IssmDouble P, IssmDouble dzMin, IssmDouble aSnow, IssmDouble C, IssmDouble V, IssmDouble Vmean, IssmDouble dIce, int sid){ /*{{{*/
     1296void accumulation(IssmDouble** pT, IssmDouble** pdz, IssmDouble** pd, IssmDouble** pW, IssmDouble** pa, IssmDouble** padiff, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pm, int aIdx, int dsnowIdx, IssmDouble Tmean, IssmDouble T_air, IssmDouble P, IssmDouble dzMin, IssmDouble aSnow, IssmDouble C, IssmDouble V, IssmDouble Vmean, IssmDouble dIce, int sid){ /*{{{*/
    12801297
    12811298        // Adds precipitation and deposition to the model grid
    … …  
    13181335        IssmDouble* W=NULL;
    13191336        IssmDouble* a=NULL;
     1337        IssmDouble* adiff=NULL;
    13201338        IssmDouble* re=NULL;
    13211339        IssmDouble* gdn=NULL;
    … …  
    13311349        W=*pW;
    13321350        a=*pa;
     1351        adiff=*padiff;
    13331352        re=*pre;
    13341353        gdn=*pgdn;
    … …  
    13891408                                newcell(&W,0.0,top,m); //new cell W
    13901409                                newcell(&a,aSnow,top,m); //new cell a
     1410                                newcell(&adiff,aSnow,top,m); //new cell a
    13911411                                newcell(&re,refall,top,m); //new cell grain size
    13921412                                newcell(&gdn,dfall,top,m); //new cell grain dendricity
    … …  
    14571477        *pW=W;
    14581478        *pa=a;
     1479        *padiff=adiff;
    14591480        *pre=re;
    14601481        *pgdn=gdn;
    … …  
    14621483        *pm=m;
    14631484} /*}}}*/
    1464 void melt(IssmDouble* pM, IssmDouble* pMs, IssmDouble* pR, IssmDouble* pmAdd, IssmDouble* pdz_add, IssmDouble** pT, IssmDouble** pd, IssmDouble** pdz, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pn, IssmDouble dzMin, IssmDouble zMax, IssmDouble zMin, IssmDouble zTop, IssmDouble zY, IssmDouble dIce, int sid){ /*{{{*/
     1485void melt(IssmDouble* pM, IssmDouble* pMs, IssmDouble* pR, IssmDouble* pmAdd, IssmDouble* pdz_add, IssmDouble** pT, IssmDouble** pd, IssmDouble** pdz, IssmDouble** pW, IssmDouble** pa, IssmDouble** padiff, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pn, IssmDouble dzMin, IssmDouble zMax, IssmDouble zMin, IssmDouble zTop, IssmDouble zY, IssmDouble dIce, int sid){ /*{{{*/
    14651486
    14661487        //// MELT ROUTINE
    … …  
    15041525        IssmDouble W_bot=0.0;
    15051526        IssmDouble a_bot=0.0;
     1527        IssmDouble adiff_bot=0.0;
    15061528        IssmDouble re_bot=0.0;
    15071529        IssmDouble gdn_bot=0.0;
    … …  
    15331555        IssmDouble* W=*pW;
    15341556        IssmDouble* a=*pa;
     1557        IssmDouble* adiff=*padiff;
    15351558        IssmDouble* re=*pre;
    15361559        IssmDouble* gdn=*pgdn;
    … …  
    17831806                celldelete(&T,n,D,D_size);
    17841807                celldelete(&a,n,D,D_size);
     1808                celldelete(&adiff,n,D,D_size);
    17851809                celldelete(&re,n,D,D_size);
    17861810                celldelete(&gdn,n,D,D_size);
    … …  
    18451869                        T[X1] = (T[X2]*m[X2] + T[X1]*m[X1]) / m_new;
    18461870                        a[X1] = (a[X2]*m[X2] + a[X1]*m[X1]) / m_new;
     1871                        adiff[X1] = (adiff[X2]*m[X2] + adiff[X1]*m[X1]) / m_new;
    18471872                        re[X1] = (re[X2]*m[X2] + re[X1]*m[X1]) / m_new;
    18481873                        gdn[X1] = (gdn[X2]*m[X2] + gdn[X1]*m[X1]) / m_new;
    … …  
    18711896        celldelete(&T,n,D,D_size);
    18721897        celldelete(&a,n,D,D_size);
     1898        celldelete(&adiff,n,D,D_size);
    18731899        celldelete(&re,n,D,D_size);
    18741900        celldelete(&gdn,n,D,D_size);
    … …  
    19081934                        cellsplit(&d, n, j,1.0);
    19091935                        cellsplit(&a, n, j,1.0);
     1936                        cellsplit(&adiff, n, j,1.0);
    19101937                        cellsplit(&EI, n, j,.5);
    19111938                        cellsplit(&EW, n, j,.5);
    … …  
    19361963                d_bot=d[n-1];
    19371964                a_bot=a[n-1];
     1965                adiff_bot=adiff[n-1];
    19381966                re_bot=re[n-1];
    19391967                gdn_bot=gdn[n-1];
    … …  
    19501978                newcell(&d,d_bot,top,n);
    19511979                newcell(&a,a_bot,top,n);
     1980                newcell(&adiff,adiff_bot,top,n);
    19521981                newcell(&re,re_bot,top,n);
    19531982                newcell(&gdn,gdn_bot,top,n);
    … …  
    19752004                celldelete(&d,n,D,D_size);
    19762005                celldelete(&a,n,D,D_size);
     2006                celldelete(&adiff,n,D,D_size);
    19772007                celldelete(&re,n,D,D_size);
    19782008                celldelete(&gdn,n,D,D_size);
    … …  
    20342064        *pW=W;
    20352065        *pa=a;
     2066        *padiff=adiff;
    20362067        *pre=re;
    20372068        *pgdn=gdn;

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

    r25993 r25997  
    2828void       GembgridInitialize(IssmDouble** pdz, int* psize, IssmDouble zTop, IssmDouble dzTop, IssmDouble zMax, IssmDouble zY);
    2929IssmDouble Marbouty(IssmDouble T, IssmDouble d, IssmDouble dT);
     30IssmDouble gardnerAlb(IssmDouble* re, IssmDouble* dz, IssmDouble* d, IssmDouble clabSnow, IssmDouble clabIce, IssmDouble SZA, IssmDouble COT, IssmDouble dPHC, int m);
    3031void grainGrowth(IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, IssmDouble* T,IssmDouble* dz,IssmDouble* d, IssmDouble* W,IssmDouble smb_dt,int m,int aIdx, int sid);
    31 void albedo(IssmDouble** a,int aIdx, IssmDouble* re, IssmDouble* dz, IssmDouble* d, IssmDouble cldFrac, IssmDouble aIce, IssmDouble aSnow, IssmDouble aValue, IssmDouble adThresh, IssmDouble* T, IssmDouble* W, IssmDouble P, IssmDouble EC, IssmDouble Msurf, IssmDouble clabSnow, IssmDouble clabIce, IssmDouble SZA, IssmDouble COT, IssmDouble t0wet, IssmDouble t0dry, IssmDouble K, IssmDouble dt, IssmDouble dIce, int m, int sid);
    32 void shortwave(IssmDouble** pswf, int swIdx, int aIdx, IssmDouble dsw, IssmDouble as, IssmDouble* d, IssmDouble* dz, IssmDouble* re, IssmDouble dIce, int m, int sid);
     32void albedo(IssmDouble** a, IssmDouble** adiff, int aIdx, IssmDouble* re, IssmDouble* dz, IssmDouble* d, IssmDouble cldFrac, IssmDouble aIce, IssmDouble aSnow, IssmDouble aValue, IssmDouble adThresh, IssmDouble* T, IssmDouble* W, IssmDouble P, IssmDouble EC, IssmDouble Msurf, IssmDouble clabSnow, IssmDouble clabIce, IssmDouble SZA, IssmDouble COT, IssmDouble t0wet, IssmDouble t0dry, IssmDouble K, IssmDouble dt, IssmDouble dIce, int m, int sid);
     33void shortwave(IssmDouble** pswf, int swIdx, int aIdx, IssmDouble dsw, IssmDouble dswdiff, IssmDouble as, IssmDouble asdiff, IssmDouble* d, IssmDouble* dz, IssmDouble* re, IssmDouble dIce, int m, int sid);
    3334void thermo(IssmDouble* pEC, IssmDouble** T, IssmDouble* pulwrf, IssmDouble* dz, IssmDouble* d, IssmDouble* swf, IssmDouble dlw, IssmDouble Ta, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble teValue, IssmDouble Ws, IssmDouble dt0, int m, IssmDouble Vz, IssmDouble Tz, IssmDouble thermo_scaling, IssmDouble dIce, int sid, bool isconstrainsurfaceT);
    34 void accumulation(IssmDouble** pT, IssmDouble** pdz, IssmDouble** pd, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pm, int aIdx, int dsnowIdx, IssmDouble Tmean, IssmDouble Ta, IssmDouble P, IssmDouble dzMin, IssmDouble aSnow, IssmDouble C, IssmDouble V, IssmDouble Vmean, IssmDouble dIce, int sid);
    35 void melt(IssmDouble* pM, IssmDouble* pMs, IssmDouble* pR, IssmDouble* pmAdd, IssmDouble* pdz_add, IssmDouble** pT, IssmDouble** pd, IssmDouble** pdz, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pn, IssmDouble dzMin, IssmDouble zMax, IssmDouble zMin, IssmDouble zTop, IssmDouble zY, IssmDouble dIce, int sid);
     35void accumulation(IssmDouble** pT, IssmDouble** pdz, IssmDouble** pd, IssmDouble** pW, IssmDouble** pa, IssmDouble** padiff, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pm, int aIdx, int dsnowIdx, IssmDouble Tmean, IssmDouble Ta, IssmDouble P, IssmDouble dzMin, IssmDouble aSnow, IssmDouble C, IssmDouble V, IssmDouble Vmean, IssmDouble dIce, int sid);
     36void melt(IssmDouble* pM, IssmDouble* pMs, IssmDouble* pR, IssmDouble* pmAdd, IssmDouble* pdz_add, IssmDouble** pT, IssmDouble** pd, IssmDouble** pdz, IssmDouble** pW, IssmDouble** pa, IssmDouble** padiff, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pn, IssmDouble dzMin, IssmDouble zMax, IssmDouble zMin, IssmDouble zTop, IssmDouble zY, IssmDouble dIce, int sid);
    3637void densification(IssmDouble** pd,IssmDouble** pdz, IssmDouble* T, IssmDouble* re, int denIdx, int aIdx, int swIdx, IssmDouble adThresh, IssmDouble C, IssmDouble dt, IssmDouble Tmean, IssmDouble dIce, int m, int sid);
    3738void turbulentFlux(IssmDouble* pshf, IssmDouble* plhf, IssmDouble* pEC, IssmDouble Ta, IssmDouble Ts, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble ds, IssmDouble Ws, IssmDouble Vz, IssmDouble Tz, IssmDouble dIce, int sid);

  • issm/trunk-jpl/src/c/shared/Enum/Enum.vim

    r25958 r25997  
    757757syn keyword cConstant SmbAccumulatedRunoffEnum
    758758syn keyword cConstant SmbAEnum
     759syn keyword cConstant SmbAdiffEnum
    759760syn keyword cConstant SmbAValueEnum
    760761syn keyword cConstant SmbAccumulationEnum
     762syn keyword cConstant SmbAdiffiniEnum
    761763syn keyword cConstant SmbAiniEnum
    762764syn keyword cConstant SmbBMaxEnum
    … …  
    765767syn keyword cConstant SmbBPosEnum
    766768syn keyword cConstant SmbCEnum
     769syn keyword cConstant SmbCcsnowValueEnum
     770syn keyword cConstant SmbCciceValueEnum
     771syn keyword cConstant SmbCotValueEnum
    767772syn keyword cConstant SmbDEnum
    768773syn keyword cConstant SmbDailyairdensityEnum
    … …  
    778783syn keyword cConstant SmbDlwrfEnum
    779784syn keyword cConstant SmbDswrfEnum
     785syn keyword cConstant SmbDswdiffrfEnum
    780786syn keyword cConstant SmbDzAddEnum
    781787syn keyword cConstant SmbDzEnum
    … …  
    830836syn keyword cConstant SmbSmbCorrEnum
    831837syn keyword cConstant SmbSmbrefEnum
     838syn keyword cConstant SmbSzaValueEnum
    832839syn keyword cConstant SmbTEnum
    833840syn keyword cConstant SmbTaEnum

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

    r25958 r25997  
    753753        SmbAccumulatedRunoffEnum,
    754754        SmbAEnum,
     755        SmbAdiffEnum,
    755756        SmbAValueEnum,
    756757        SmbAccumulationEnum,
     758        SmbAdiffiniEnum,
    757759        SmbAiniEnum,
    758760        SmbBMaxEnum,
    … …  
    761763        SmbBPosEnum,
    762764        SmbCEnum,
     765        SmbCcsnowValueEnum,
     766        SmbCciceValueEnum,
     767        SmbCotValueEnum,
    763768        SmbDEnum,
    764769        SmbDailyairdensityEnum,
    … …  
    774779        SmbDlwrfEnum,
    775780        SmbDswrfEnum,
     781        SmbDswdiffrfEnum,
    776782        SmbDzAddEnum,
    777783        SmbDzEnum,
    … …  
    827833        SmbSmbCorrEnum,
    828834        SmbSmbrefEnum,
     835        SmbSzaValueEnum,
    829836        SmbTEnum,
    830837        SmbTaEnum,

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

    r25958 r25997  
    759759                case SmbAccumulatedRunoffEnum : return "SmbAccumulatedRunoff";
    760760                case SmbAEnum : return "SmbA";
     761                case SmbAdiffEnum : return "SmbAdiff";
    761762                case SmbAValueEnum : return "SmbAValue";
    762763                case SmbAccumulationEnum : return "SmbAccumulation";
     764                case SmbAdiffiniEnum : return "SmbAdiffini";
    763765                case SmbAiniEnum : return "SmbAini";
    764766                case SmbBMaxEnum : return "SmbBMax";
    … …  
    767769                case SmbBPosEnum : return "SmbBPos";
    768770                case SmbCEnum : return "SmbC";
     771                case SmbCcsnowValueEnum : return "SmbCcsnowValue";
     772                case SmbCciceValueEnum : return "SmbCciceValue";
     773                case SmbCotValueEnum : return "SmbCotValue";
    769774                case SmbDEnum : return "SmbD";
    770775                case SmbDailyairdensityEnum : return "SmbDailyairdensity";
    … …  
    780785                case SmbDlwrfEnum : return "SmbDlwrf";
    781786                case SmbDswrfEnum : return "SmbDswrf";
     787                case SmbDswdiffrfEnum : return "SmbDswdiffrf";
    782788                case SmbDzAddEnum : return "SmbDzAdd";
    783789                case SmbDzEnum : return "SmbDz";
    … …  
    832838                case SmbSmbCorrEnum : return "SmbSmbCorr";
    833839                case SmbSmbrefEnum : return "SmbSmbref";
     840                case SmbSzaValueEnum : return "SmbSzaValue";
    834841                case SmbTEnum : return "SmbT";
    835842                case SmbTaEnum : return "SmbTa";

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

    r25958 r25997  
    777777              else if (strcmp(name,"SmbAccumulatedRunoff")==0) return SmbAccumulatedRunoffEnum;
    778778              else if (strcmp(name,"SmbA")==0) return SmbAEnum;
     779              else if (strcmp(name,"SmbAdiff")==0) return SmbAdiffEnum;
    779780              else if (strcmp(name,"SmbAValue")==0) return SmbAValueEnum;
    780781              else if (strcmp(name,"SmbAccumulation")==0) return SmbAccumulationEnum;
     782              else if (strcmp(name,"SmbAdiffini")==0) return SmbAdiffiniEnum;
    781783              else if (strcmp(name,"SmbAini")==0) return SmbAiniEnum;
    782784              else if (strcmp(name,"SmbBMax")==0) return SmbBMaxEnum;
    … …  
    785787              else if (strcmp(name,"SmbBPos")==0) return SmbBPosEnum;
    786788              else if (strcmp(name,"SmbC")==0) return SmbCEnum;
     789              else if (strcmp(name,"SmbCcsnowValue")==0) return SmbCcsnowValueEnum;
     790              else if (strcmp(name,"SmbCciceValue")==0) return SmbCciceValueEnum;
     791              else if (strcmp(name,"SmbCotValue")==0) return SmbCotValueEnum;
    787792              else if (strcmp(name,"SmbD")==0) return SmbDEnum;
    788793              else if (strcmp(name,"SmbDailyairdensity")==0) return SmbDailyairdensityEnum;
    … …  
    798803              else if (strcmp(name,"SmbDlwrf")==0) return SmbDlwrfEnum;
    799804              else if (strcmp(name,"SmbDswrf")==0) return SmbDswrfEnum;
     805              else if (strcmp(name,"SmbDswdiffrf")==0) return SmbDswdiffrfEnum;
    800806              else if (strcmp(name,"SmbDzAdd")==0) return SmbDzAddEnum;
    801807              else if (strcmp(name,"SmbDz")==0) return SmbDzEnum;
    … …  
    850856              else if (strcmp(name,"SmbSmbCorr")==0) return SmbSmbCorrEnum;
    851857              else if (strcmp(name,"SmbSmbref")==0) return SmbSmbrefEnum;
     858              else if (strcmp(name,"SmbSzaValue")==0) return SmbSzaValueEnum;
    852859              else if (strcmp(name,"SmbT")==0) return SmbTEnum;
    853860              else if (strcmp(name,"SmbTa")==0) return SmbTaEnum;
    … …  
    868875              else if (strcmp(name,"SmbZMax")==0) return SmbZMaxEnum;
    869876              else if (strcmp(name,"SmbZMin")==0) return SmbZMinEnum;
    870               else if (strcmp(name,"SmbZTop")==0) return SmbZTopEnum;
     877         else stage=8;
     878   }
     879   if(stage==8){
     880              if (strcmp(name,"SmbZTop")==0) return SmbZTopEnum;
    871881              else if (strcmp(name,"SmbZY")==0) return SmbZYEnum;
    872882              else if (strcmp(name,"SolidearthExternalDisplacementEastRate")==0) return SolidearthExternalDisplacementEastRateEnum;
    … …  
    875885              else if (strcmp(name,"SolidearthExternalGeoidRate")==0) return SolidearthExternalGeoidRateEnum;
    876886              else if (strcmp(name,"SolidearthExternalBarystaticSeaLevelRate")==0) return SolidearthExternalBarystaticSeaLevelRateEnum;
    877          else stage=8;
    878    }
    879    if(stage==8){
    880               if (strcmp(name,"StrainRateeffective")==0) return StrainRateeffectiveEnum;
     887              else if (strcmp(name,"StrainRateeffective")==0) return StrainRateeffectiveEnum;
    881888              else if (strcmp(name,"StrainRateparallel")==0) return StrainRateparallelEnum;
    882889              else if (strcmp(name,"StrainRateperpendicular")==0) return StrainRateperpendicularEnum;
    … …  
    991998              else if (strcmp(name,"Outputdefinition52")==0) return Outputdefinition52Enum;
    992999              else if (strcmp(name,"Outputdefinition53")==0) return Outputdefinition53Enum;
    993               else if (strcmp(name,"Outputdefinition54")==0) return Outputdefinition54Enum;
     1000         else stage=9;
     1001   }
     1002   if(stage==9){
     1003              if (strcmp(name,"Outputdefinition54")==0) return Outputdefinition54Enum;
    9941004              else if (strcmp(name,"Outputdefinition55")==0) return Outputdefinition55Enum;
    9951005              else if (strcmp(name,"Outputdefinition56")==0) return Outputdefinition56Enum;
    … …  
    9981008              else if (strcmp(name,"Outputdefinition59")==0) return Outputdefinition59Enum;
    9991009              else if (strcmp(name,"Outputdefinition5")==0) return Outputdefinition5Enum;
    1000          else stage=9;
    1001    }
    1002    if(stage==9){
    1003               if (strcmp(name,"Outputdefinition60")==0) return Outputdefinition60Enum;
     1010              else if (strcmp(name,"Outputdefinition60")==0) return Outputdefinition60Enum;
    10041011              else if (strcmp(name,"Outputdefinition61")==0) return Outputdefinition61Enum;
    10051012              else if (strcmp(name,"Outputdefinition62")==0) return Outputdefinition62Enum;
    … …  
    11141121              else if (strcmp(name,"DepthAverageAnalysis")==0) return DepthAverageAnalysisEnum;
    11151122              else if (strcmp(name,"DeviatoricStressErrorEstimator")==0) return DeviatoricStressErrorEstimatorEnum;
    1116               else if (strcmp(name,"Divergence")==0) return DivergenceEnum;
     1123         else stage=10;
     1124   }
     1125   if(stage==10){
     1126              if (strcmp(name,"Divergence")==0) return DivergenceEnum;
    11171127              else if (strcmp(name,"Domain3Dsurface")==0) return Domain3DsurfaceEnum;
    11181128              else if (strcmp(name,"DoubleArrayInput")==0) return DoubleArrayInputEnum;
    … …  
    11211131              else if (strcmp(name,"DoubleMatArrayParam")==0) return DoubleMatArrayParamEnum;
    11221132              else if (strcmp(name,"DoubleMatExternalResult")==0) return DoubleMatExternalResultEnum;
    1123          else stage=10;
    1124    }
    1125    if(stage==10){
    1126               if (strcmp(name,"DoubleMatParam")==0) return DoubleMatParamEnum;
     1133              else if (strcmp(name,"DoubleMatParam")==0) return DoubleMatParamEnum;
    11271134              else if (strcmp(name,"DoubleParam")==0) return DoubleParamEnum;
    11281135              else if (strcmp(name,"DoubleVecParam")==0) return DoubleVecParamEnum;
    … …  
    12371244              else if (strcmp(name,"MantlePlumeGeothermalFlux")==0) return MantlePlumeGeothermalFluxEnum;
    12381245              else if (strcmp(name,"MassFlux")==0) return MassFluxEnum;
    1239               else if (strcmp(name,"Masscon")==0) return MassconEnum;
     1246         else stage=11;
     1247   }
     1248   if(stage==11){
     1249              if (strcmp(name,"Masscon")==0) return MassconEnum;
    12401250              else if (strcmp(name,"Massconaxpby")==0) return MassconaxpbyEnum;
    12411251              else if (strcmp(name,"Massfluxatgate")==0) return MassfluxatgateEnum;
    … …  
    12441254              else if (strcmp(name,"Matdamageice")==0) return MatdamageiceEnum;
    12451255              else if (strcmp(name,"Matenhancedice")==0) return MatenhancediceEnum;
    1246          else stage=11;
    1247    }
    1248    if(stage==11){
    1249               if (strcmp(name,"Materials")==0) return MaterialsEnum;
     1256              else if (strcmp(name,"Materials")==0) return MaterialsEnum;
    12501257              else if (strcmp(name,"Matestar")==0) return MatestarEnum;
    12511258              else if (strcmp(name,"Matice")==0) return MaticeEnum;
    … …  
    13601367              else if (strcmp(name,"SpatialLinearFloatingMeltRate")==0) return SpatialLinearFloatingMeltRateEnum;
    13611368              else if (strcmp(name,"SpcDynamic")==0) return SpcDynamicEnum;
    1362               else if (strcmp(name,"SpcStatic")==0) return SpcStaticEnum;
     1369         else stage=12;
     1370   }
     1371   if(stage==12){
     1372              if (strcmp(name,"SpcStatic")==0) return SpcStaticEnum;
    13631373              else if (strcmp(name,"SpcTransient")==0) return SpcTransientEnum;
    13641374              else if (strcmp(name,"Sset")==0) return SsetEnum;
    … …  
    13671377              else if (strcmp(name,"StressIntensityFactor")==0) return StressIntensityFactorEnum;
    13681378              else if (strcmp(name,"StressbalanceAnalysis")==0) return StressbalanceAnalysisEnum;
    1369          else stage=12;
    1370    }
    1371    if(stage==12){
    1372               if (strcmp(name,"StressbalanceConvergenceNumSteps")==0) return StressbalanceConvergenceNumStepsEnum;
     1379              else if (strcmp(name,"StressbalanceConvergenceNumSteps")==0) return StressbalanceConvergenceNumStepsEnum;
    13731380              else if (strcmp(name,"StressbalanceSIAAnalysis")==0) return StressbalanceSIAAnalysisEnum;
    13741381              else if (strcmp(name,"StressbalanceSolution")==0) return StressbalanceSolutionEnum;

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

    r25933 r25997  
    4141
    4242                %optional inputs:
    43                 aValue = NaN; %Albedo forcing at every element.  Used only if aIdx == 0.
     43                aValue = NaN; %Albedo forcing at every element.  Used only if aIdx == 0, or density exceeds adThresh
    4444                teValue = NaN; %Outward longwave radiation thermal emissivity forcing at every element (default in code is 1)
    4545
    … …  
    5353                Wini = NaN; %Water content (kg m-2)
    5454                Aini = NaN; %albedo (0-1)
     55                Adiffini = NaN; %albedo, diffusive radiation (0-1)
    5556                Tini = NaN; %snow temperature (K)
    5657                Sizeini = NaN; %Number of layers
    … …  
    5859                %settings:
    5960                aIdx   = NaN; %method for calculating albedo and subsurface absorption (default is 1)
    60                               % 0: direct input from aValue parameter
    61                               % 1: effective grain radius [Gardner & Sharp, 2009]
    62                               % 2: effective grain radius [Brun et al., 1992; LeFebre et al., 2003]], with swIdx=1, SW penetration follows grain size in 3 spectral bands (Brun et al., 1992)
    63                               % 3: density and cloud amount [Greuell & Konzelmann, 1994]
    64                               % 4: exponential time decay & wetness [Bougamont & Bamber, 2005]
     61                % 0: direct input from aValue parameter
     62                % 1: effective grain radius [Gardner & Sharp, 2009]
     63                % 2: effective grain radius [Brun et al., 1992; LeFebre et al., 2003]], with swIdx=1, SW penetration follows grain size in 3 spectral bands (Brun et al., 1992)
     64                % 3: density and cloud amount [Greuell & Konzelmann, 1994]
     65                % 4: exponential time decay & wetness [Bougamont & Bamber, 2005]
    6566
    6667                swIdx  = NaN; % apply all SW to top grid cell (0) or allow SW to penetrate surface (1) (default 1, with snow density (taken from Bassford, 2004))
    6768
    6869                denIdx = NaN; %densification model to use (default is 2):
    69                               % 1 = emperical model of Herron and Langway (1980)
    70                               % 2 = semi-emperical model of Anthern et al. (2010)
    71                               % 3 = DO NOT USE: physical model from Appendix B of Anthern et al. (2010)
    72                               % 4 = DO NOT USE: emperical model of Li and Zwally (2004)
    73                               % 5 = DO NOT USE: modified emperical model (4) by Helsen et al. (2008)
    74                               % 6 = Antarctica semi-emperical model of Ligtenberg et al. (2011)
    75                               % 7 = Greenland semi-emperical model of Kuipers Munneke et al. (2015)
     70                % 1 = emperical model of Herron and Langway (1980)
     71                % 2 = semi-emperical model of Anthern et al. (2010)
     72                % 3 = DO NOT USE: physical model from Appendix B of Anthern et al. (2010)
     73                % 4 = DO NOT USE: emperical model of Li and Zwally (2004)
     74                % 5 = DO NOT USE: modified emperical model (4) by Helsen et al. (2008)
     75                % 6 = Antarctica semi-emperical model of Ligtenberg et al. (2011)
     76                % 7 = Greenland semi-emperical model of Kuipers Munneke et al. (2015)
    7677
    7778                dsnowIdx = NaN; %model for fresh snow accumulation density (default is 1):
    78                                 % 0 = Original GEMB value, 150 kg/m^3
    79                                 % 1 = Antarctica value of fresh snow density, 350 kg/m^3
    80                                 % 2 = Greenland value of fresh snow density, 315 kg/m^3, Fausto et al. (2008)
    81                                 % 3 = Antarctica model of Kaspers et al. (2004)
    82                                 % 4 = Greenland model of Kuipers Munneke et al. (2015)
     79                % 0 = Original GEMB value, 150 kg/m^3
     80                % 1 = Antarctica value of fresh snow density, 350 kg/m^3
     81                % 2 = Greenland value of fresh snow density, 315 kg/m^3, Fausto et al. (2008)
     82                % 3 = Antarctica model of Kaspers et al. (2004)
     83                % 4 = Greenland model of Kuipers Munneke et al. (2015)
    8384
    8485                zTop  = NaN; % depth over which grid length is constant at the top of the snopack (default 10) [m]
    … …  
    9293
    9394                %specific albedo parameters:
     95                %Method 1
     96                dswdiffrf = NaN; %downward diffusive shortwave radiation flux [W/m^2]
     97                szaValue = NaN; %Solar Zenith Angle [degree]
     98                cotValue = NaN; %Cloud Optical Thickness
     99                ccsnowValue = NaN; %concentration of light absorbing carbon for snow [ppm1]
     100                cciceValue = NaN; %concentration of light absorbing carbon for ice [ppm1]
    94101                %Method 1 and 2:
    95102                aSnow = NaN; % new snow albedo (0.64 - 0.89)
    96103                aIce  = NaN; % range 0.27-0.58 for old snow
    97                              %Method 3: Radiation Correction Factors -> only used for met station data and Greuell & Konzelmann, 1994 albedo
     104                %Method 3: Radiation Correction Factors -> only used for met station data and Greuell & Konzelmann, 1994 albedo
    98105                cldFrac = NaN; % average cloud amount
    99                                %Method 4: additonal tuning parameters albedo as a funtion of age and water content (Bougamont et al., 2005)
     106                %Method 4: additonal tuning parameters albedo as a funtion of age and water content (Bougamont et al., 2005)
    100107                t0wet = NaN; % time scale for wet snow (15-21.9)
    101108                t0dry = NaN; % warm snow timescale (30)
    102109                K     = NaN; % time scale temperature coef. (7)
    103110                adThresh = NaN; %Apply aIdx method to all areas with densities below this value,
    104                                 %or else apply direct input value from aValue, allowing albedo to be altered.
    105                                 %Default value is rho water (1023 kg m-3).
     111                %or else apply direct input value from aValue, allowing albedo to be altered.
     112                %Default value is rho water (1023 kg m-3).
    106113
    107114                %densities:
    … …  
    134141                function self = extrude(self,md) % {{{
    135142
    136                         self.Ta=project3d(md,'vector',self.Ta,'type','node');
    137                         self.V=project3d(md,'vector',self.V,'type','node');
    138                         self.dswrf=project3d(md,'vector',self.dswrf,'type','node');
    139                         self.dswrf=project3d(md,'vector',self.dswrf,'type','node');
    140                         self.P=project3d(md,'vector',self.P,'type','node');
    141                         self.eAir=project3d(md,'vector',self.eAir,'type','node');
    142                         self.pAir=project3d(md,'vector',self.pAir,'type','node');
    143 
     143                        self.Ta=project3d(md,'vector',self.Ta,'type','element');
     144                        self.V=project3d(md,'vector',self.V,'type','element');
     145                        self.dswrf=project3d(md,'vector',self.dswrf,'type','element');
     146                        self.dslrf=project3d(md,'vector',self.dslrf,'type','element');
     147                        self.P=project3d(md,'vector',self.P,'type','element');
     148                        self.eAir=project3d(md,'vector',self.eAir,'type','element');
     149                        self.pAir=project3d(md,'vector',self.pAir,'type','element');
     150
     151                        if ~isnan(self.Dzini)
     152                                self.Dzini=project3d(md,'vector',self.Dzini,'type','element');
     153                        end
     154                        if ~isnan(self.Dini)
     155                                self.Dini=project3d(md,'vector',self.Dini,'type','element');
     156                        end
     157                        if ~isnan(self.Reini)
     158                                self.Reini=project3d(md,'vector',self.Reini,'type','element');
     159                        end
     160                        if ~isnan(self.Gdnini)
     161                                self.Gdnini=project3d(md,'vector',self.Gdnini,'type','element');
     162                        end
     163                        if ~isnan(self.Gspini)
     164                                self.Gspini=project3d(md,'vector',self.Gspini,'type','element');
     165                        end
     166                        if ~isnan(self.ECini)
     167                                self.ECini=project3d(md,'vector',self.ECini,'type','element');
     168                        end
     169                        if ~isnan(self.Wini)
     170                                self.Wini=project3d(md,'vector',self.Wini,'type','element');
     171                        end
     172                        if ~isnan(self.Aini)
     173                                self.Aini=project3d(md,'vector',self.Aini,'type','element');
     174                        end
     175                        if ~isnan(self.Adiffini)
     176                                self.Adiffini=project3d(md,'vector',self.Adiffini,'type','element');
     177                        end
     178                        if ~isnan(self.Tini)
     179                                self.Tini=project3d(md,'vector',self.Tini,'type','element');
     180                        end
     181
     182                        if ~isnan(self.dswdiffrf)
     183                                self.dswdiffrf=project3d(md,'vector',self.dswdiffrf,'type','element');
     184                        end
     185                        if ~isnan(self.szaValue)
     186                                self.szaValue=project3d(md,'vector',self.szaValue,'type','element');
     187                        end
     188                        if ~isnan(self.cotValue)
     189                                self.cotValue=project3d(md,'vector',self.cotValue,'type','element');
     190                        end
     191                        if ~isnan(self.ccsnowValue)
     192                                self.ccsnowValue=project3d(md,'vector',self.ccsnowValue,'type','element');
     193                        end
     194                        if ~isnan(self.cciceValue)
     195                                self.cciceValue=project3d(md,'vector',self.cciceValue,'type','element');
     196                        end
    144197                        if (aIdx == 0) & ~isnan(self.aValue)
    145                                 self.aValue=project3d(md,'vector',self.aValue,'type','node');
     198                                self.aValue=project3d(md,'vector',self.aValue,'type','element');
    146199                        end
    147200                        if ~isnan(self.teValue)
    148                                 self.teValue=project3d(md,'vector',self.teValue,'type','node');
     201                                self.teValue=project3d(md,'vector',self.teValue,'type','element');
    149202                        end
    150203
    … …  
    152205                end % }}}
    153206                function list = defaultoutputs(self,md) % {{{
    154                         list = {'SmbMassBalance'};
     207                        list = {'SmbMassBalance','SmbAccumulatedMassBalance'};
    155208                end % }}}
    156209                function self = setdefaultparameters(self,mesh,geometry) % {{{
    … …  
    196249                        self.aValue = self.aSnow*ones(mesh.numberofelements,1);
    197250
     251                        self.dswdiffrf=0.0*ones(mesh.numberofelements,1);
     252                        self.szaValue=0.0*ones(mesh.numberofelements,1);
     253                        self.cotValue=0.0*ones(mesh.numberofelements,1);
     254                        self.ccsnowValue=0.0*ones(mesh.numberofelements,1);
     255                        self.cciceValue=0.0*ones(mesh.numberofelements,1);
     256
    198257                        self.Dzini=0.05*ones(mesh.numberofelements,2);
    199258                        self.Dini=910.0*ones(mesh.numberofelements,2);
    … …  
    204263                        self.Wini=0.0*ones(mesh.numberofelements,2);
    205264                        self.Aini=self.aSnow*ones(mesh.numberofelements,2);
     265                        self.Adiffini=ones(mesh.numberofelements,2);
    206266                        self.Tini=273.15*ones(mesh.numberofelements,2);
    207267                        %               /!\ Default value of Tini must be equal to Tmean but don't know Tmean yet (computed when atmospheric forcings are interpolated on mesh).
    … …  
    226286                        md = checkfield(md,'fieldname','smb.V','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<',45,'size',size(self.Ta)); %max 500 km/h
    227287                        md = checkfield(md,'fieldname','smb.dswrf','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<=',1400,'size',size(self.Ta));
     288                        md = checkfield(md,'fieldname','smb.dswdiffrf','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<=',1400);
    228289                        md = checkfield(md,'fieldname','smb.dlwrf','timeseries',1,'NaN',1,'Inf',1,'>=',0,'size',size(self.Ta));
    229290                        md = checkfield(md,'fieldname','smb.P','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<=',200,'size',size(self.Ta));
    … …  
    232293                        if (self.isclimatology)
    233294                                md = checkfield(md,'fieldname', 'smb.Ta', 'size',[md.mesh.numberofelements+1],...
    234                                                 'message',['Ta must have md.mesh.numberofelements+1 rows in order to force a climatology']);
     295                                        'message',['Ta must have md.mesh.numberofelements+1 rows in order to force a climatology']);
    235296                                md = checkfield(md,'fieldname', 'smb.V', 'size',[md.mesh.numberofelements+1],...
    236                                                 'message',['V must have md.mesh.numberofelements+1 rows in order to force a climatology']);
     297                                        'message',['V must have md.mesh.numberofelements+1 rows in order to force a climatology']);
    237298                                md = checkfield(md,'fieldname', 'smb.dswrf', 'size',[md.mesh.numberofelements+1],...
    238                                                 'message',['dswrf must have md.mesh.numberofelements+1 rows in order to force a climatology']);
     299                                        'message',['dswrf must have md.mesh.numberofelements+1 rows in order to force a climatology']);
    239300                                md = checkfield(md,'fieldname', 'smb.dlwrf', 'size',[md.mesh.numberofelements+1],...
    240                                                 'message',['dlwrf must have md.mesh.numberofelements+1 rows in order to force a climatology']);
     301                                        'message',['dlwrf must have md.mesh.numberofelements+1 rows in order to force a climatology']);
    241302                                md = checkfield(md,'fieldname', 'smb.P', 'size',[md.mesh.numberofelements+1],...
    242                                                 'message',['P must have md.mesh.numberofelements+1 rows in order to force a climatology']);
     303                                        'message',['P must have md.mesh.numberofelements+1 rows in order to force a climatology']);
    243304                                md = checkfield(md,'fieldname', 'smb.eAir', 'size',[md.mesh.numberofelements+1],...
    244                                                 'message',['eAir must have md.mesh.numberofelements+1 rows in order to force a climatology']);
     305                                        'message',['eAir must have md.mesh.numberofelements+1 rows in order to force a climatology']);
    245306                        end
    246307
    … …  
    273334                                        md = checkfield(md,'fieldname','smb.aSnow','NaN',1,'Inf',1,'>=',.64,'<=',.89);
    274335                                        md = checkfield(md,'fieldname','smb.aIce','NaN',1,'Inf',1,'>=',.27,'<=',.58);
     336                                        if self.aIdx==1
     337                                                md = checkfield(md,'fieldname','smb.szaValue','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<=',90);
     338                                                md = checkfield(md,'fieldname','smb.cotValue','timeseries',1,'NaN',1,'Inf',1,'>=',0);
     339                                                md = checkfield(md,'fieldname','smb.ccsnowValue','timeseries',1,'NaN',1,'Inf',1,'>=',0);
     340                                                md = checkfield(md,'fieldname','smb.cciceValue','timeseries',1,'NaN',1,'Inf',1,'>=',0);
     341                                        end
    275342                                case 3
    276343                                        md = checkfield(md,'fieldname','smb.cldFrac','NaN',1,'Inf',1,'>=',0,'<=',1);
    … …  
    307374                        fielddisplay(self,'Ta','2 m air temperature, in Kelvin');
    308375                        fielddisplay(self,'V','wind speed (m s-1)');
    309                         fielddisplay(self,'dlwrf','downward shortwave radiation flux [W/m^2]');
    310                         fielddisplay(self,'dswrf','downward longwave radiation flux [W/m^2]');
     376                        fielddisplay(self,'dswrf','downward shortwave radiation flux [W/m^2]');
     377                        fielddisplay(self,'dswdiffrf','downward diffusive portion of shortwave radiation flux (default to 0) [W/m^2]');
     378                        fielddisplay(self,'dlwrf','downward longwave radiation flux [W/m^2]');
    311379                        fielddisplay(self,'P','precipitation [mm w.e. / m^2]');
    312380                        fielddisplay(self,'eAir','screen level vapor pressure [Pa]');
    … …  
    328396                        fielddisplay(self,'adThresh',{'Apply aIdx method to all areas with densities below this value,','or else apply direct input value from aValue, allowing albedo to be altered.'});
    329397                        fielddisplay(self,'aIdx',{'method for calculating albedo and subsurface absorption (default is 1)',...
    330                                             '0: direct input from aValue parameter',...
    331                                             '1: effective grain radius [Gardner & Sharp, 2009]',...
    332                                             '2: effective grain radius [Brun et al., 1992; LeFebre et al., 2003], with swIdx=1, SW penetration follows grain size in 3 spectral bands (Brun et al., 1992)',...
    333                                             '3: density and cloud amount [Greuell & Konzelmann, 1994]',...
    334                                             '4: exponential time decay & wetness [Bougamont & Bamber, 2005]'})
     398                                '0: direct input from aValue parameter',...
     399                                '1: effective grain radius [Gardner & Sharp, 2009]',...
     400                                '2: effective grain radius [Brun et al., 1992; LeFebre et al., 2003], with swIdx=1, SW penetration follows grain size in 3 spectral bands (Brun et al., 1992)',...
     401                                '3: density and cloud amount [Greuell & Konzelmann, 1994]',...
     402                                '4: exponential time decay & wetness [Bougamont & Bamber, 2005]'})
    335403
    336404                        fielddisplay(self,'teValue','Outward longwave radiation thermal emissivity forcing at every element (default in code is 1)');
    … …  
    345413                        fielddisplay(self,'Wini','Initial snow water content when restart [kg m-2]');
    346414                        fielddisplay(self,'Aini','Initial albedo when restart [-]');
     415                        fielddisplay(self,'Adiffini','Initial diffusive radiation albedo when restart (default to 1) [-]');
    347416                        fielddisplay(self,'Tini','Initial snow temperature when restart [K]');
    348417                        fielddisplay(self,'Sizeini','Initial number of layers when restart [-]');
    349418
    350419                        %additional albedo parameters:
    351                         fielddisplay(self,'aValue','Albedo forcing at every element.');
     420                        fielddisplay(self,'aValue','Albedo forcing at every element');
    352421                        switch self.aIdx
    353422                                case {1 2}
    354423                                        fielddisplay(self,'aSnow','new snow albedo (0.64 - 0.89)');
    355424                                        fielddisplay(self,'aIce','albedo of ice (0.27-0.58)');
     425                                        if self.aIdx==1
     426                                                fielddisplay(self,'szaValue','Solar Zenith Angle [degree]');
     427                                                fielddisplay(self,'cotValue','Cloud Optical Thickness');
     428                                                fielddisplay(self,'ccsnowValue','concentration of light absorbing carbon for snow [ppm1]');
     429                                                fielddisplay(self,'cciceValue','concentration of light absorbing carbon for snow [ppm1]');
     430                                        end
    356431                                case 3
    357432                                        fielddisplay(self,'cldFrac','average cloud amount');
    … …  
    364439                        fielddisplay(self,'swIdx','apply all SW to top grid cell (0) or allow SW to penetrate surface (1) [default 1, with snow density (taken from Bassford, 2004)]');
    365440                        fielddisplay(self,'denIdx',{'densification model to use (default is 2):',...
    366                                             '1 = emperical model of Herron and Langway (1980)',...
    367                                             '2 = semi-emperical model of Anthern et al. (2010)',...
    368                                             '3 = DO NOT USE: physical model from Appendix B of Anthern et al. (2010)',...
    369                                             '4 = DO NOT USE: emperical model of Li and Zwally (2004)',...
    370                                             '5 = DO NOT USE: modified emperical model (4) by Helsen et al. (2008)',...
    371                                             '6 = Antarctica semi-emperical model of Ligtenberg et al. (2011)',...
    372                                             '7 = Greenland semi-emperical model of Kuipers Munneke et al. (2015)'});
     441                                '1 = emperical model of Herron and Langway (1980)',...
     442                                '2 = semi-emperical model of Anthern et al. (2010)',...
     443                                '3 = DO NOT USE: physical model from Appendix B of Anthern et al. (2010)',...
     444                                '4 = DO NOT USE: emperical model of Li and Zwally (2004)',...
     445                                '5 = DO NOT USE: modified emperical model (4) by Helsen et al. (2008)',...
     446                                '6 = Antarctica semi-emperical model of Ligtenberg et al. (2011)',...
     447                                '7 = Greenland semi-emperical model of Kuipers Munneke et al. (2015)'});
    373448                        fielddisplay(self,'dsnowIdx',{'model for fresh snow accumulation density (default is 1):',...
    374                                             '0 = Original GEMB value, 150 kg/m^3',...
    375                                             '1 = Antarctica value of fresh snow density, 350 kg/m^3',...
    376                                             '2 = Greenland value of fresh snow density, 315 kg/m^3, Fausto et al. (2008)',...
    377                                             '3 = Antarctica model of Kaspers et al. (2004), Make sure to set Vmean accurately',...
    378                                             '4 = Greenland model of Kuipers Munneke et al. (2015)'});
     449                                '0 = Original GEMB value, 150 kg/m^3',...
     450                                '1 = Antarctica value of fresh snow density, 350 kg/m^3',...
     451                                '2 = Greenland value of fresh snow density, 315 kg/m^3, Fausto et al. (2008)',...
     452                                '3 = Antarctica model of Kaspers et al. (2004), Make sure to set Vmean accurately',...
     453                                '4 = Greenland model of Kuipers Munneke et al. (2015)'});
    379454
    380455                        fielddisplay(self, 'steps_per_step', 'number of smb steps per time step');
    … …  
    407482                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','V','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts);
    408483                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','dswrf','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts);
     484                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','dswdiffrf','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts);
    409485                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','dlwrf','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts);
    410486                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','P','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts);
    … …  
    441517                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','aValue','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts);
    442518                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','teValue','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts);
     519                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','szaValue','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts);
     520                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','cotValue','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts);
     521                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','ccsnowValue','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts);
     522                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','cciceValue','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts);
    443523
    444524                        %snow properties init
    … …  
    451531                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','Wini','format','DoubleMat','mattype',3);
    452532                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','Aini','format','DoubleMat','mattype',3);
     533                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','Adiffini','format','DoubleMat','mattype',3);
    453534                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','Tini','format','DoubleMat','mattype',3);
    454535                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','Sizeini','format','IntMat','mattype',2);
    … …  
    464545                                error('All GEMB forcings (Ta, P, V, dswrf, dlwrf, eAir, pAir) must have the same time steps in the final row!');
    465546                        end
     547                        if size(md.smb.teValue,2)>1 & any(md.smb.teValue(end,:) - md.smb.Ta(end,:) ~= 0)
     548                                error('If GEMB forcing teValue is transient, it must have the same time steps as input Ta in the final row!');
     549                        end
     550                        if size(md.smb.dswdiffrf,2)>1 & any(md.smb.dswdiffrf(end,:) - md.smb.Ta(end,:) ~= 0)
     551                                error('If GEMB forcing dswdiffrf is transient, it must have the same time steps as input Ta in the final row!');
     552                        end
     553                        if size(md.smb.aValue,2)>1 & any(md.smb.aValue(end,:) - md.smb.Ta(end,:) ~= 0)
     554                                error('If GEMB forcing aValue is transient, it must have the same time steps as input Ta in the final row!');
     555                        end
     556                        if size(md.smb.szaValue,2)>1 & any(md.smb.szaValue(end,:) - md.smb.Ta(end,:) ~= 0)
     557                                error('If GEMB forcing szaValue is transient, it must have the same time steps as input Ta in the final row!');
     558                        end
     559                        if size(md.smb.cotValue,2)>1 & any(md.smb.cotValue(end,:) - md.smb.Ta(end,:) ~= 0)
     560                                error('If GEMB forcing cotValue is transient, it must have the same time steps as input Ta in the final row!');
     561                        end
     562                        if size(md.smb.ccsnowValue,2)>1 & any(md.smb.ccsnowValue(end,:) - md.smb.Ta(end,:) ~= 0)
     563                                error('If GEMB forcing ccsnowValue is transient, it must have the same time steps as input Ta in the final row!');
     564                        end
     565                        if size(md.smb.cciceValue,2)>1 & any(md.smb.cciceValue(end,:) - md.smb.Ta(end,:) ~= 0)
     566                                error('If GEMB forcing cciceValue is transient, it must have the same time steps as input Ta in the final row!');
     567                        end
    466568                        time=self.Ta(end,:); %assume all forcings are on the same time step
    467569                        dtime=diff(time,1);

  • issm/trunk-jpl/src/m/classes/SMBgemb.py

    r25933 r25997  
    4747
    4848        #optional inputs:
    49         self.aValue                 = np.nan    # Albedo forcing at every element.  Used only if aIdx == 0.
     49        self.aValue                 = np.nan    # Albedo forcing at every element.  Used only if aIdx == 0, or density exceeds adThresh
    5050        self.teValue                = np.nan    # Outward longwave radiation thermal emissivity forcing at every element (default in code is 1)
    5151
    … …  
    5959        self.Wini                   = np.nan    # Water content (kg m-2)
    6060        self.Aini                   = np.nan    # albedo (0-1)
     61        self.Adiffini               = np.nan    # albedo, diffusive radiation (0-1)
    6162        self.Tini                   = np.nan    # snow temperature (K)
    6263        self.Sizeini                = np.nan    # Number of layers
    … …  
    9697
    9798        #specific albedo parameters:
     99        #Method 1
     100        dswdiffrf                   = np.nan    #downward diffusive shortwave radiation flux [W/m^2]
     101        szaValue                    = np.nan    #Solar Zenith Angle [degree]
     102        cotValue                    = np.nan    #Cloud Optical Thickness
     103        ccsnowValue                 = np.nan    #concentration of light absorbing carbon for snow [ppm1]
     104        cciceValue                  = np.nan    #concentration of light absorbing carbon for ice [ppm1]
    98105        #Method 1 and 2:
    99106        self.aSnow                  = np.nan    # new snow albedo (0.64 - 0.89)
    … …  
    152159        string = "%s\n%s" % (string, fielddisplay(self, 'Ta', '2 m air temperature, in Kelvin'))
    153160        string = "%s\n%s" % (string, fielddisplay(self, 'V', 'wind speed (m s-1)'))
    154         string = "%s\n%s" % (string, fielddisplay(self, 'dlwrf', 'downward shortwave radiation flux [W/m^2]'))
    155         string = "%s\n%s" % (string, fielddisplay(self, 'dswrf', 'downward longwave radiation flux [W/m^2]'))
     161        string = "%s\n%s" % (string, fielddisplay(self, 'dswrf', 'downward shortwave radiation flux [W/m^2]'))
     162        string = "%s\n%s" % (string, fielddisplay(self, 'dswdiffrf', 'downward diffusive portion of shortwave radiation flux (default to 0) [W/m^2]'))
     163        string = "%s\n%s" % (string, fielddisplay(self, 'dlwrf', 'downward longwave radiation flux [W/m^2]'))
    156164        string = "%s\n%s" % (string, fielddisplay(self, 'P', 'precipitation [mm w.e. / m^2]'))
    157165        string = "%s\n%s" % (string, fielddisplay(self, 'eAir', 'screen level vapor pressure [Pa]'))
    … …  
    188196        string = "%s\n%s" % (string, fielddisplay(self, 'Wini', 'Initial snow water content when restart [kg m-2]'))
    189197        string = "%s\n%s" % (string, fielddisplay(self, 'Aini', 'Initial albedo when restart [-]'))
     198        string = "%s\n%s" % (string, fielddisplay(self, 'Adiffini', 'Initial diffusive radiation albedo when restart (default to 1) [-]'))
    190199        string = "%s\n%s" % (string, fielddisplay(self, 'Tini', 'Initial snow temperature when restart [K]'))
    191200        string = "%s\n%s" % (string, fielddisplay(self, 'Sizeini', 'Initial number of layers when restart [-]'))
    192201
    193202        #additional albedo parameters:
    194         string = "%s\n%s" % (string, fielddisplay(self, 'aValue', 'Albedo forcing at every element.'))
     203        string = "%s\n%s" % (string, fielddisplay(self, 'aValue', 'Albedo forcing at every element'))
    195204        if isinstance(self.aIdx, (list, type(np.array([1, 2])))) and (self.aIdx == [1, 2] or (1 in self.aIdx and 2 in self.aIdx)):
    196205            string = "%s\n%s" % (string, fielddisplay(self, 'aSnow', 'new snow albedo (0.64 - 0.89)'))
    197206            string = "%s\n%s" % (string, fielddisplay(self, 'aIce', 'albedo of ice (0.27-0.58)'))
     207            if self.aIdx == 1:
     208                string = "%s\n%s" % (string, fielddisplay(self,'szaValue','Solar Zenith Angle [degree]'))
     209                string = "%s\n%s" % (string, fielddisplay(self,'cotValue','Cloud Optical Thickness'))
     210                string = "%s\n%s" % (string, fielddisplay(self,'ccsnowValue','concentration of light absorbing carbon for snow [ppm1]'))
     211                string = "%s\n%s" % (string, fielddisplay(self,'cciceValue','concentration of light absorbing carbon for snow [ppm1]'))
     212
    198213        elif self.aIdx == 3:
    199214            string = "%s\n%s" % (string, fielddisplay(self, 'cldFrac', 'average cloud amount'))
    … …  
    228243
    229244    def extrude(self, md):  # {{{
    230         self.Ta = project3d(md, 'vector', self.Ta, 'type', 'node')
    231         self.V = project3d(md, 'vector', self.V, 'type', 'node')
    232         self.dswrf = project3d(md, 'vector', self.dswrf, 'type', 'node')
    233         self.dswrf = project3d(md, 'vector', self.dswrf, 'type', 'node')
    234         self.P = project3d(md, 'vector', self.P, 'type', 'node')
    235         self.eAir = project3d(md, 'vector', self.eAir, 'type', 'node')
    236         self.pAir = project3d(md, 'vector', self.pAir, 'type', 'node')
    237 
    238         if (self.aIdx == 0) and np.isnan(self.aValue):
    239             self.aValue = project3d(md, 'vector', self.aValue, 'type', 'node')
    240         if np.isnan(self.teValue):
    241             self.teValue = project3d(md, 'vector', self.teValue, 'type', 'node')
     245        self.Ta = project3d(md, 'vector', self.Ta, 'type', 'element')
     246        self.V = project3d(md, 'vector', self.V, 'type', 'element')
     247        self.dswrf = project3d(md, 'vector', self.dswrf, 'type', 'element')
     248        self.dswrf = project3d(md, 'vector', self.dswrf, 'type', 'element')
     249        self.P = project3d(md, 'vector', self.P, 'type', 'element')
     250        self.eAir = project3d(md, 'vector', self.eAir, 'type', 'element')
     251        self.pAir = project3d(md, 'vector', self.pAir, 'type', 'element')
     252
     253        if not np.isnan(self.Dzini):
     254            self.self.Dzini=project3d(md,'vector',self.self.Dzini,'type','element');
     255        if not np.isnan(self.Dini):
     256            self.self.Dini=project3d(md,'vector',self.Dini,'type','element');
     257        if not np.isnan(self.Reini):
     258            self.self.Reini=project3d(md,'vector',self.Reini,'type','element');
     259        if not np.isnan(self.Gdnini):
     260            self.Gdnini=project3d(md,'vector',self.Gdnini,'type','element');
     261        if not np.isnan(self.Gspini):
     262            self.Gspini=project3d(md,'vector',self.Gspini,'type','element');
     263        if not np.isnan(self.ECini):
     264            self.ECini=project3d(md,'vector',self.ECini,'type','element');
     265        if not np.isnan(self.Wini):
     266            self.Wini=project3d(md,'vector',self.Wini,'type','element');
     267        if not np.isnan(self.Aini):
     268            self.Aini=project3d(md,'vector',self.Aini,'type','element');
     269        if not np.isnan(self.Adiffini):
     270            self.Adiffini=project3d(md,'vector',self.Adiffini,'type','element');
     271        if not np.isnan(self.Tini):
     272            self.Tini=project3d(md,'vector',self.Tini,'type','element');
     273        if not np.isnan(self.dswdiffrf):
     274            self.dswdiffrf=project3d(md,'vector',self.dswdiffrf,'type','element');
     275        if not np.isnan(self.szaValue):
     276            self.szaValue=project3d(md,'vector',self.szaValue,'type','element');
     277        if not np.isnan(self.cotValue):
     278            self.cotValue=project3d(md,'vector',self.cotValue,'type','element');
     279        if not np.isnan(self.ccsnowValue):
     280            self.ccsnowValue=project3d(md,'vector',self.ccsnowValue,'type','element');
     281        if not np.isnan(self.cciceValue):
     282            self.cciceValue=project3d(md,'vector',self.cciceValue,'type','element');
     283
     284        if (self.aIdx == 0) and (not np.isnan(self.aValue)):
     285            self.aValue = project3d(md, 'vector', self.aValue, 'type', 'element')
     286        if not np.isnan(self.teValue):
     287            self.teValue = project3d(md, 'vector', self.teValue, 'type', 'element')
    242288
    243289        return self
    … …  
    245291
    246292    def defaultoutputs(self, md):  # {{{
    247         return ['SmbMassBalance']
     293        return ['SmbMassBalance','SmbAccumulatedMassBalance']
    248294    #}}}
    249295
    … …  
    290336        self.aValue = self.aSnow * np.ones(mesh.numberofelements,)
    291337
     338        self.dswdiffrf = 0.0 * np.ones(mesh.numberofelements,)
     339        self.szaValue = 0.0 * np.ones(mesh.numberofelements,)
     340        self.cotValue = 0.0 * np.ones(mesh.numberofelements,)
     341        self.ccsnowValue = 0.0 * np.ones(mesh.numberofelements,)
     342        self.cciceValue = 0.0 * np.ones(mesh.numberofelements,)
     343
    292344        self.Dzini = 0.05 * np.ones((mesh.numberofelements, 2))
    293345        self.Dini = 910.0 * np.ones((mesh.numberofelements, 2))
    … …  
    298350        self.Wini = 0.0 * np.ones((mesh.numberofelements, 2))
    299351        self.Aini = self.aSnow * np.ones((mesh.numberofelements, 2))
     352        self.Adiffini = np.ones((mesh.numberofelements, 2))
    300353        self.Tini = 273.15 * np.ones((mesh.numberofelements, 2))
    301354#       /!\ Default value of Tini must be equal to Tmean but don't know Tmean yet (computed when atmospheric forcings are interpolated on mesh).
    … …  
    321374        md = checkfield(md, 'fieldname', 'smb.V', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '> = ', 0, '<', 45, 'size', np.shape(self.Ta))  #max 500 km/h
    322375        md = checkfield(md, 'fieldname', 'smb.dswrf', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '> = ', 0, '< = ', 1400, 'size', np.shape(self.Ta))
     376        md = checkfield(md, 'fieldname', 'smb.dswdiffrf', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '> = ', 0, '< = ', 1400)
    323377        md = checkfield(md, 'fieldname', 'smb.dlwrf', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '> = ', 0, 'size', np.shape(self.Ta))
    324378        md = checkfield(md, 'fieldname', 'smb.P', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '> = ', 0, '< = ', 200, 'size', np.shape(self.Ta))
    … …  
    358412            md = checkfield(md, 'fieldname', 'smb.aSnow', 'NaN', 1, 'Inf', 1, '> = ', .64, '< = ', .89)
    359413            md = checkfield(md, 'fieldname', 'smb.aIce', 'NaN', 1, 'Inf', 1, '> = ', .27, '< = ', .58)
     414            if self.aIdx == 1:
     415                md = checkfield(md,'fieldname','smb.szaValue','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<=',90)
     416                md = checkfield(md,'fieldname','smb.cotValue','timeseries',1,'NaN',1,'Inf',1,'>=',0)
     417                md = checkfield(md,'fieldname','smb.ccsnowValue','timeseries',1,'NaN',1,'Inf',1,'>=',0)
     418                md = checkfield(md,'fieldname','smb.cciceValue','timeseries',1,'NaN',1,'Inf',1,'>=',0)
     419
    360420        elif self.aIdx == 0:
    361421            md = checkfield(md, 'fieldname', 'smb.aValue', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '>=', 0, '<=', 1)
    … …  
    396456        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'V', 'format', 'DoubleMat', 'mattype', 2, 'timeserieslength', md.mesh.numberofelements + 1, 'yts', yts)
    397457        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'dswrf', 'format', 'DoubleMat', 'mattype', 2, 'timeserieslength', md.mesh.numberofelements + 1, 'yts', yts)
     458        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'dswdiffrf', 'format', 'DoubleMat', 'mattype', 2, 'timeserieslength', md.mesh.numberofelements + 1, 'yts', yts)
    398459        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'dlwrf', 'format', 'DoubleMat', 'mattype', 2, 'timeserieslength', md.mesh.numberofelements + 1, 'yts', yts)
    399460        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'P', 'format', 'DoubleMat', 'mattype', 2, 'timeserieslength', md.mesh.numberofelements + 1, 'yts', yts)
    … …  
    430491        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'aValue', 'format', 'DoubleMat', 'mattype', 2, 'timeserieslength', md.mesh.numberofelements + 1, 'yts', yts)
    431492        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'teValue', 'format', 'DoubleMat', 'mattype', 2, 'timeserieslength', md.mesh.numberofelements + 1, 'yts', yts)
     493        WriteData(fid,prefix,'object',self,'class','smb','fieldname','szaValue','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts)
     494        WriteData(fid,prefix,'object',self,'class','smb','fieldname','cotValue','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts)
     495        WriteData(fid,prefix,'object',self,'class','smb','fieldname','ccsnowValue','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts)
     496        WriteData(fid,prefix,'object',self,'class','smb','fieldname','cciceValue','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts)
    432497
    433498        #snow properties init
    … …  
    440505        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'Wini', 'format', 'DoubleMat', 'mattype', 3)
    441506        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'Aini', 'format', 'DoubleMat', 'mattype', 3)
     507        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'Adiffini', 'format', 'DoubleMat', 'mattype', 3)
    442508        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'Tini', 'format', 'DoubleMat', 'mattype', 3)
    443509        WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'Sizeini', 'format', 'IntMat', 'mattype', 2)
    … …  
    448514            raise IOError('All GEMB forcings (Ta, P, V, dswrf, dlwrf, eAir, pAir) must have the same time steps in the final row!')
    449515
     516        if ((np.ndim(self.teValue)>1) & np.any(self.teValue[-1] - self.Ta[-1] != 0)):
     517            raise IOError('If GEMB forcing teValue is transient, it must have the same time steps as input Ta in the final row!')
     518        if ((np.ndim(self.dswdiffrf)>1) & np.any(self.dswdiffrf[-1] - self.Ta[-1] != 0)):
     519            raise IOError('If GEMB forcing dswdiffrf is transient, it must have the same time steps as input Ta in the final row!')
     520        if ((np.ndim(self.aValue)>1) & np.any(self.aValue[-1] - self.Ta[-1] != 0)):
     521            raise IOError('If GEMB forcing aValue is transient, it must have the same time steps as input Ta in the final row!')
     522        if ((np.ndim(self.szaValue)>1) & np.any(self.szaValue[-1] - self.Ta[-1] != 0)):
     523            raise IOError('If GEMB forcing szaValue is transient, it must have the same time steps as input Ta in the final row!')
     524        if ((np.ndim(self.cotValue)>1) & np.any(self.cotValue[-1] - self.Ta[-1] != 0)):
     525            raise IOError('If GEMB forcing cotValue is transient, it must have the same time steps as input Ta in the final row!')
     526        if ((np.ndim(self.ccsnowValue)>1) & np.any(self.ccsnowValue[-1] - self.Ta[-1] != 0)):
     527            raise IOError('If GEMB forcing ccsnowValue is transient, it must have the same time steps as input Ta in the final row!')
     528        if ((np.ndim(self.cciceValue)>1) & np.any(self.cciceValue[-1] - self.Ta[-1] != 0)):
     529            raise IOError('If GEMB forcing cciceValue is transient, it must have the same time steps as input Ta in the final row!')
     530
    450531        time = self.Ta[-1]  #assume all forcings are on the same time step
    451532        dtime = np.diff(time, n=1, axis=0)
Note: See TracChangeset for help on using the changeset viewer.