Index: /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp (revision 25996) +++ /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp (revision 25997) @@ -46,4 +46,5 @@ iomodel->FetchDataToInput(inputs,elements,"md.smb.V",SmbVEnum); iomodel->FetchDataToInput(inputs,elements,"md.smb.dswrf",SmbDswrfEnum); + iomodel->FetchDataToInput(inputs,elements,"md.smb.dswdiffrf",SmbDswdiffrfEnum); iomodel->FetchDataToInput(inputs,elements,"md.smb.dlwrf",SmbDlwrfEnum); iomodel->FetchDataToInput(inputs,elements,"md.smb.P",SmbPEnum); @@ -70,8 +71,13 @@ iomodel->FetchDataToInput(inputs,elements,"md.smb.Wini",SmbWiniEnum); iomodel->FetchDataToInput(inputs,elements,"md.smb.Aini",SmbAiniEnum); + iomodel->FetchDataToInput(inputs,elements,"md.smb.Adiffini",SmbAdiffiniEnum); iomodel->FetchDataToInput(inputs,elements,"md.smb.Tini",SmbTiniEnum); iomodel->FetchDataToInput(inputs,elements,"md.smb.Sizeini",SmbSizeiniEnum); iomodel->FetchDataToInput(inputs,elements,"md.smb.aValue",SmbAValueEnum); iomodel->FetchDataToInput(inputs,elements,"md.smb.teValue",SmbTeValueEnum); + iomodel->FetchDataToInput(inputs,elements,"md.smb.szaValue",SmbSzaValueEnum); + iomodel->FetchDataToInput(inputs,elements,"md.smb.cotValue",SmbCotValueEnum); + iomodel->FetchDataToInput(inputs,elements,"md.smb.ccsnowValue",SmbCcsnowValueEnum); + iomodel->FetchDataToInput(inputs,elements,"md.smb.cciceValue",SmbCciceValueEnum); break; case SMBpddEnum: Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 25996) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 25997) @@ -3546,4 +3546,5 @@ IssmDouble dlw=0.0; IssmDouble dsw=0.0; + IssmDouble dswdiff=0.0; IssmDouble P=0.0; IssmDouble eAir=0.0; @@ -3602,4 +3603,5 @@ IssmDouble* W = NULL; IssmDouble* a = NULL; + IssmDouble* adiff = NULL; IssmDouble* swf=NULL; IssmDouble* T = NULL; @@ -3617,4 +3619,5 @@ IssmDouble* Wini = NULL; IssmDouble* aini = NULL; + IssmDouble* adiffini = NULL; IssmDouble* Tini = NULL; int m=0; @@ -3695,4 +3698,5 @@ this->inputs->GetArray(SmbWiniEnum,this->lid,&Wini,&m); this->inputs->GetArray(SmbAiniEnum,this->lid,&aini,&m); + this->inputs->GetArray(SmbAdiffiniEnum,this->lid,&adiffini,&m); this->inputs->GetArray(SmbTiniEnum,this->lid,&Tini,&m); EC_input = this->GetInput(SmbECiniEnum); _assert_(EC_input); @@ -3714,4 +3718,5 @@ W = xNew(m); for(int i=0;i(m); for(int i=0;i(m); for(int i=0;i(m); for(int i=0;i(m);for(int i=0;i(m);for(int i=0;i(m);for(int i=0;i(m);for(int i=0;iinputs->GetArray(SmbWEnum,this->lid,&W,&m); this->inputs->GetArray(SmbAEnum,this->lid,&a,&m); + this->inputs->GetArray(SmbAdiffEnum,this->lid,&adiff,&m); this->inputs->GetArray(SmbTEnum,this->lid,&T,&m); EC_input = this->GetInput(SmbECDtEnum); _assert_(EC_input); @@ -3821,4 +3828,5 @@ Input *Dlwr_input= this->GetInput(SmbDlwrfEnum,timeinputs); _assert_(Dlwr_input); Input *Dswr_input= this->GetInput(SmbDswrfEnum,timeinputs); _assert_(Dswr_input); + Input *Dswrdiff_input= this->GetInput(SmbDswdiffrfEnum,timeinputs); _assert_(Dswrdiff_input); Input *P_input = this->GetInput(SmbPEnum,timeinputs); _assert_(P_input); Input *eAir_input= this->GetInput(SmbEAirEnum,timeinputs); _assert_(eAir_input); @@ -3826,4 +3834,8 @@ Input *teValue_input= this->GetInput(SmbTeValueEnum,timeinputs); _assert_(teValue_input); Input *aValue_input= this->GetInput(SmbAValueEnum,timeinputs); _assert_(aValue_input); + Input *szaValue_input= this->GetInput(SmbSzaValueEnum,timeinputs); _assert_(szaValue_input); + Input *cotValue_input= this->GetInput(SmbCotValueEnum,timeinputs); _assert_(cotValue_input); + Input *ccsnowValue_input= this->GetInput(SmbCcsnowValueEnum,timeinputs); _assert_(ccsnowValue_input); + Input *cciceValue_input= this->GetInput(SmbCciceValueEnum,timeinputs); _assert_(cciceValue_input); /*extract daily data:{{{*/ @@ -3832,4 +3844,5 @@ Dlwr_input->GetInputValue(&dlw,gauss); //downward longwave radiation flux [W m-2] Dswr_input->GetInputValue(&dsw,gauss); //downward shortwave radiation flux [W m-2] + Dswrdiff_input->GetInputValue(&dswdiff,gauss); //downward shortwave diffuse radiation flux [W m-2] P_input->GetInputValue(&P,gauss); //precipitation [kg m-2] eAir_input->GetInputValue(&eAir,gauss); //screen level vapor pressure [Pa] @@ -3837,9 +3850,8 @@ teValue_input->GetInputValue(&teValue,gauss); // Emissivity [0-1] aValue_input->GetInputValue(&aValue,gauss); // Albedo [0 1] - //Not implemented as input - //szaValue; // Solar Zenith Angle [degree] - //cotValue; // Cloud Optical Thickness - //ccsnowValue;//concentration of light absorbing carbon for snow [ppm1] - //cciceValue;//concentration of light absorbing carbon for ice [ppm1] + szaValue_input->GetInputValue(&szaValue,gauss); // Solar Zenith Angle [degree] + cotValue_input->GetInputValue(&cotValue,gauss); // Cloud Optical Thickness + ccsnowValue_input->GetInputValue(&ccsnowValue,gauss); //concentration of light absorbing carbon for snow [ppm1] + cciceValue_input->GetInputValue(&cciceValue,gauss); //concentration of light absorbing carbon for ice [ppm1] //_printf_("Time: " << t << " Ta: " << Ta << " V: " << V << " dlw: " << dlw << " dsw: " << dsw << " P: " << P << " eAir: " << eAir << " pAir: " << pAir << "\n"); /*}}}*/ @@ -3849,8 +3861,8 @@ /*Snow, firn and ice albedo:*/ - 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()); + 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()); /*Distribution of absorbed short wave radation with depth:*/ - if(isshortwave)shortwave(&swf, swIdx, aIdx, dsw, a[0], d, dz, re,rho_ice,m,this->Sid()); + if(isshortwave)shortwave(&swf, swIdx, aIdx, dsw, dswdiff, a[0], adiff[0], d, dz, re,rho_ice,m,this->Sid()); /*Calculate net shortwave [W m-2]*/ @@ -3869,9 +3881,9 @@ /*Add snow/rain to top grid cell adjusting cell depth, temperature and density*/ - 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()); + 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()); /*Calculate water production, M [kg m-2] resulting from snow/ice temperature exceeding 273.15 deg K * (> 0 deg C), runoff R [kg m-2] and resulting changes in density and determine wet compaction [m]*/ - 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()); + 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()); /*Allow non-melt densification and determine compaction [m]*/ @@ -3903,4 +3915,5 @@ << "lwf[" << netLW << "] " << "a[" << a << "] " + << "adiff[" << adiff << "] " << "te[" << teValue << "] " << "\n"); @@ -3975,4 +3988,5 @@ this->inputs->SetArrayInput(SmbWEnum,this->lid,W,m); this->inputs->SetArrayInput(SmbAEnum,this->lid,a,m); + this->inputs->SetArrayInput(SmbAdiffEnum,this->lid,adiff,m); this->SetElementInput(SmbECEnum,sumEC/dt/rho_ice); this->SetElementInput(SmbMassBalanceEnum,(sumP + sumEC -sumR)/dt/rho_ice); @@ -4001,4 +4015,5 @@ if(W) xDelete(W); if(a) xDelete(a); + if(adiff) xDelete(adiff); if(T) xDelete(T); if(dzini) xDelete(dzini); @@ -4009,4 +4024,5 @@ if(Wini) xDelete(Wini); if(aini) xDelete(aini); + if(adiffini) xDelete(adiffini); if(Tini) xDelete(Tini); if(swf) xDelete(swf); Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp (revision 25996) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp (revision 25997) @@ -197,4 +197,80 @@ } /*}}}*/ +IssmDouble gardnerAlb(IssmDouble* re, IssmDouble* dz, IssmDouble* d, IssmDouble clabSnow, IssmDouble clabIce, IssmDouble SZA, IssmDouble COT, IssmDouble dPHC, int m){ /*{{{*/ + //gardnerAlb(S1, c1, SZA, t, z1, S2, c2) + //This is an implementation of the snow and ice broadband albedo + // parameterization developed by Alex Gardner. + //Created By: Alex S. Gardner, Jet Propulsion Laboratory [alex.s.gardner@jpl.nasa.gov] + // Last Modified: June, 2014 + //Full Reference: Gardner, A. S., and Sharp, M. J.: A review of snow and + // ice albedo and the development of a new physically based broadband albedo + // parameterization, J. Geophys. Res., 115, F01009, 10.1029/2009jf001444, + // 2010. + + //INPUTS + // ONE LAYER + // - S1 : specific surface area of the snow or ice [cm^2 g-1] + // - c1 : concentration of light absorbing carbon [ppm1] + // - SZA : solar zenith angle of the incident radiation [deg] + // - t : cloud optical thickness + // TWO LAYER + // - z1 : depth of snow suface layer [mm w.e.] + // - S2 : specific surface area of bottom ice layer [cm^2 g-1] + // - c2 : concentration of light absorbing carbon of bottom ice + // layer [ppm1] + IssmDouble c1=clabSnow; + IssmDouble c2=clabIce; + IssmDouble t=COT; + IssmDouble a=0.0; + + //Single layer albedo parameterization + //convert effective radius to specific surface area [cm2 g-1] + IssmDouble S1 = 3.0 / (0.091 * re[0]); + + //effective solar zenith angle + IssmDouble x = min(pow(t/(3*cos(Pi*SZA/180.0)),0.5), 1.0); + IssmDouble u = 0.64*x + (1-x)*cos(Pi*SZA/180.0); + + // pure snow albedo + IssmDouble as = 1.48 - pow(S1,-0.07); + + //change in pure snow albedo due to soot loading + 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))); + + //Two layer albedo parameterization + // do two layer calculation if there is more than 1 layer + IssmDouble z1=0.0; + int lice=0; + for(int l=0;(l0 & lice Dtol){ + // determine albedo values for bottom layer + IssmDouble S2 = 3.0 / (0.091 * re[lice]); + + // pure snow albedo + IssmDouble as2 = 1.48 - pow(S2,-0.07); + + // change in pure snow albedo due to soot loading + 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))); + + // determine the effective change due to finite depth and soot loading + IssmDouble A = min(1.0, (2.1 * pow(z1,1.35*(1-as) - 0.1*c1 - 0.13))); + + dac = (as2 + dac2 - as) + A*((as + dac) - (as2 + dac2)); + } + + // change in albedo due to solar zenith angle + IssmDouble dasz = 0.53*as*(1 - (as + dac))*pow(1 - u,1.2); + + // change in albedo due to cloud (apart from change in diffuse fraction) + IssmDouble dat = (0.1*t*pow(as + dac,1.3)) / (pow(1 + 1.5*t,as)); + + // Broadband albedo + a = as + dac + dasz + dat; + + return a; +} /*}}}*/ void grainGrowth(IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, IssmDouble* T,IssmDouble* dz,IssmDouble* d, IssmDouble* W,IssmDouble smb_dt,int m,int aIdx,int sid){ /*{{{*/ @@ -410,5 +486,5 @@ } /*}}}*/ -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) { /*{{{*/ +void 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) { /*{{{*/ //// Calculates Snow, firn and ice albedo as a function of: @@ -471,4 +547,5 @@ /*output: */ IssmDouble* a=NULL; + IssmDouble* adiff=NULL; if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_(" albedo module\n"); @@ -476,4 +553,5 @@ /*Recover pointers: */ a=*pa; + adiff=*padiff; //some constants: @@ -490,76 +568,7 @@ if(aIdx==1){ //function of effective grain radius - - //gardnerAlb(S1, c1, SZA, t, z1, S2, c2) - //This is an implementation of the snow and ice broadband albedo - // parameterization developed by Alex Gardner. - //Created By: Alex S. Gardner, Jet Propulsion Laboratory [alex.s.gardner@jpl.nasa.gov] - // Last Modified: June, 2014 - //Full Reference: Gardner, A. S., and Sharp, M. J.: A review of snow and - // ice albedo and the development of a new physically based broadband albedo - // parameterization, J. Geophys. Res., 115, F01009, 10.1029/2009jf001444, - // 2010. - - //INPUTS - // ONE LAYER - // - S1 : specific surface area of the snow or ice [cm^2 g-1] - // - c1 : concentration of light absorbing carbon [ppm1] - // - SZA : solar zenith angle of the incident radiation [deg] - // - t : cloud optical thickness - // TWO LAYER - // - z1 : depth of snow suface layer [mm w.e.] - // - S2 : specific surface area of bottom ice layer [cm^2 g-1] - // - c2 : concentration of light absorbing carbon of bottom ice - // layer [ppm1] - IssmDouble c1=clabSnow; - IssmDouble c2=clabIce; - IssmDouble t=COT; - - //Single layer albedo parameterization - //convert effective radius to specific surface area [cm2 g-1] - IssmDouble S1 = 3.0 / (0.091 * re[0]); - - //effective solar zenith angle - IssmDouble x = min(pow(t/(3*cos(Pi*SZA/180.0)),0.5), 1.0); - IssmDouble u = 0.64*x + (1-x)*cos(Pi*SZA/180.0); - - // pure snow albedo - IssmDouble as = 1.48 - pow(S1,-0.07); - - //change in pure snow albedo due to soot loading - 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))); - - //Two layer albedo parameterization - // do two layer calculation if there is more than 1 layer - IssmDouble z1=0.0; - int lice=0; - for(int l=0;(l0 & lice Dtol){ - // determine albedo values for bottom layer - IssmDouble S2 = 3.0 / (0.091 * re[lice]); - - // pure snow albedo - IssmDouble as2 = 1.48 - pow(S2,-0.07); - - // change in pure snow albedo due to soot loading - 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))); - - // determine the effective change due to finite depth and soot loading - IssmDouble A = min(1.0, (2.1 * pow(z1,1.35*(1-as) - 0.1*c1 - 0.13))); - - dac = (as2 + dac2 - as) + A*((as + dac) - (as2 + dac2)); - } - - // change in albedo due to solar zenith angle - IssmDouble dasz = 0.53*as*(1 - (as + dac))*pow(1 - u,1.2); - - // change in albedo due to cloud (apart from change in diffuse fraction) - IssmDouble dat = (0.1*t*pow(as + dac,1.3)) / (pow(1 + 1.5*t,as)); - - // Broadband albedo - a[0] = as + dac + dasz + dat; + // clabSnow, IssmDouble clabIce, IssmDouble SZA, IssmDouble COT, int m + a[0]=gardnerAlb(re, dz, d, clabSnow, clabIce, SZA, COT, dPHC, m); + adiff[0]=gardnerAlb(re, dz, d, clabSnow, clabIce, 50, COT, dPHC, m); } else if(aIdx==2){ @@ -698,4 +707,5 @@ /*Assign output pointers:*/ *pa=a; + *padiff=adiff; } /*}}}*/ @@ -1078,5 +1088,5 @@ } /*}}}*/ -void shortwave(IssmDouble** pswf, int swIdx, int aIdx, IssmDouble dsw, IssmDouble as, IssmDouble* d, IssmDouble* dz, IssmDouble* re, IssmDouble dIce, int m, int sid){ /*{{{*/ +void 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){ /*{{{*/ // DISTRIBUTES ABSORBED SHORTWAVE RADIATION WITHIN SNOW/ICE @@ -1092,5 +1102,7 @@ // aIdx = method for calculating albedo (1-4) // dsw = downward shortwave radiative flux [w m-2] + // dswdir = downward shortwave direct radiative flux [w m-2] // as = surface albedo + // asdiff = surface albedo for diffuse radiation // d = grid cell density [kg m-3] // dz = grid cell depth [m] @@ -1113,5 +1125,10 @@ // calculate surface shortwave radiation fluxes [W m-2] - swf[0] = (1.0 - as) * dsw; + if (aIdx == 1){ + swf[0] = (1.0 - as) * dsw + (1.0 - asdiff) * dswdiff; + } + else{ + swf[0] = (1.0 - as) * dsw; + } } else{ // sw radation is absorbed at depth within the glacier @@ -1277,5 +1294,5 @@ } /*}}}*/ -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){ /*{{{*/ +void 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){ /*{{{*/ // Adds precipitation and deposition to the model grid @@ -1318,4 +1335,5 @@ IssmDouble* W=NULL; IssmDouble* a=NULL; + IssmDouble* adiff=NULL; IssmDouble* re=NULL; IssmDouble* gdn=NULL; @@ -1331,4 +1349,5 @@ W=*pW; a=*pa; + adiff=*padiff; re=*pre; gdn=*pgdn; @@ -1389,4 +1408,5 @@ newcell(&W,0.0,top,m); //new cell W newcell(&a,aSnow,top,m); //new cell a + newcell(&adiff,aSnow,top,m); //new cell a newcell(&re,refall,top,m); //new cell grain size newcell(&gdn,dfall,top,m); //new cell grain dendricity @@ -1457,4 +1477,5 @@ *pW=W; *pa=a; + *padiff=adiff; *pre=re; *pgdn=gdn; @@ -1462,5 +1483,5 @@ *pm=m; } /*}}}*/ -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){ /*{{{*/ +void 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){ /*{{{*/ //// MELT ROUTINE @@ -1504,4 +1525,5 @@ IssmDouble W_bot=0.0; IssmDouble a_bot=0.0; + IssmDouble adiff_bot=0.0; IssmDouble re_bot=0.0; IssmDouble gdn_bot=0.0; @@ -1533,4 +1555,5 @@ IssmDouble* W=*pW; IssmDouble* a=*pa; + IssmDouble* adiff=*padiff; IssmDouble* re=*pre; IssmDouble* gdn=*pgdn; @@ -1783,4 +1806,5 @@ celldelete(&T,n,D,D_size); celldelete(&a,n,D,D_size); + celldelete(&adiff,n,D,D_size); celldelete(&re,n,D,D_size); celldelete(&gdn,n,D,D_size); @@ -1845,4 +1869,5 @@ T[X1] = (T[X2]*m[X2] + T[X1]*m[X1]) / m_new; a[X1] = (a[X2]*m[X2] + a[X1]*m[X1]) / m_new; + adiff[X1] = (adiff[X2]*m[X2] + adiff[X1]*m[X1]) / m_new; re[X1] = (re[X2]*m[X2] + re[X1]*m[X1]) / m_new; gdn[X1] = (gdn[X2]*m[X2] + gdn[X1]*m[X1]) / m_new; @@ -1871,4 +1896,5 @@ celldelete(&T,n,D,D_size); celldelete(&a,n,D,D_size); + celldelete(&adiff,n,D,D_size); celldelete(&re,n,D,D_size); celldelete(&gdn,n,D,D_size); @@ -1908,4 +1934,5 @@ cellsplit(&d, n, j,1.0); cellsplit(&a, n, j,1.0); + cellsplit(&adiff, n, j,1.0); cellsplit(&EI, n, j,.5); cellsplit(&EW, n, j,.5); @@ -1936,4 +1963,5 @@ d_bot=d[n-1]; a_bot=a[n-1]; + adiff_bot=adiff[n-1]; re_bot=re[n-1]; gdn_bot=gdn[n-1]; @@ -1950,4 +1978,5 @@ newcell(&d,d_bot,top,n); newcell(&a,a_bot,top,n); + newcell(&adiff,adiff_bot,top,n); newcell(&re,re_bot,top,n); newcell(&gdn,gdn_bot,top,n); @@ -1975,4 +2004,5 @@ celldelete(&d,n,D,D_size); celldelete(&a,n,D,D_size); + celldelete(&adiff,n,D,D_size); celldelete(&re,n,D,D_size); celldelete(&gdn,n,D,D_size); @@ -2034,4 +2064,5 @@ *pW=W; *pa=a; + *padiff=adiff; *pre=re; *pgdn=gdn; Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h (revision 25996) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h (revision 25997) @@ -28,10 +28,11 @@ void GembgridInitialize(IssmDouble** pdz, int* psize, IssmDouble zTop, IssmDouble dzTop, IssmDouble zMax, IssmDouble zY); IssmDouble Marbouty(IssmDouble T, IssmDouble d, IssmDouble dT); +IssmDouble gardnerAlb(IssmDouble* re, IssmDouble* dz, IssmDouble* d, IssmDouble clabSnow, IssmDouble clabIce, IssmDouble SZA, IssmDouble COT, IssmDouble dPHC, int m); void grainGrowth(IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, IssmDouble* T,IssmDouble* dz,IssmDouble* d, IssmDouble* W,IssmDouble smb_dt,int m,int aIdx, int sid); -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); -void shortwave(IssmDouble** pswf, int swIdx, int aIdx, IssmDouble dsw, IssmDouble as, IssmDouble* d, IssmDouble* dz, IssmDouble* re, IssmDouble dIce, int m, int sid); +void 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); +void 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); void 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); -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); -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); +void 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); +void 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); void 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); void 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); Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 25996) +++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 25997) @@ -757,6 +757,8 @@ syn keyword cConstant SmbAccumulatedRunoffEnum syn keyword cConstant SmbAEnum +syn keyword cConstant SmbAdiffEnum syn keyword cConstant SmbAValueEnum syn keyword cConstant SmbAccumulationEnum +syn keyword cConstant SmbAdiffiniEnum syn keyword cConstant SmbAiniEnum syn keyword cConstant SmbBMaxEnum @@ -765,4 +767,7 @@ syn keyword cConstant SmbBPosEnum syn keyword cConstant SmbCEnum +syn keyword cConstant SmbCcsnowValueEnum +syn keyword cConstant SmbCciceValueEnum +syn keyword cConstant SmbCotValueEnum syn keyword cConstant SmbDEnum syn keyword cConstant SmbDailyairdensityEnum @@ -778,4 +783,5 @@ syn keyword cConstant SmbDlwrfEnum syn keyword cConstant SmbDswrfEnum +syn keyword cConstant SmbDswdiffrfEnum syn keyword cConstant SmbDzAddEnum syn keyword cConstant SmbDzEnum @@ -830,4 +836,5 @@ syn keyword cConstant SmbSmbCorrEnum syn keyword cConstant SmbSmbrefEnum +syn keyword cConstant SmbSzaValueEnum syn keyword cConstant SmbTEnum syn keyword cConstant SmbTaEnum Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 25996) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 25997) @@ -753,6 +753,8 @@ SmbAccumulatedRunoffEnum, SmbAEnum, + SmbAdiffEnum, SmbAValueEnum, SmbAccumulationEnum, + SmbAdiffiniEnum, SmbAiniEnum, SmbBMaxEnum, @@ -761,4 +763,7 @@ SmbBPosEnum, SmbCEnum, + SmbCcsnowValueEnum, + SmbCciceValueEnum, + SmbCotValueEnum, SmbDEnum, SmbDailyairdensityEnum, @@ -774,4 +779,5 @@ SmbDlwrfEnum, SmbDswrfEnum, + SmbDswdiffrfEnum, SmbDzAddEnum, SmbDzEnum, @@ -827,4 +833,5 @@ SmbSmbCorrEnum, SmbSmbrefEnum, + SmbSzaValueEnum, SmbTEnum, SmbTaEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 25996) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 25997) @@ -759,6 +759,8 @@ case SmbAccumulatedRunoffEnum : return "SmbAccumulatedRunoff"; case SmbAEnum : return "SmbA"; + case SmbAdiffEnum : return "SmbAdiff"; case SmbAValueEnum : return "SmbAValue"; case SmbAccumulationEnum : return "SmbAccumulation"; + case SmbAdiffiniEnum : return "SmbAdiffini"; case SmbAiniEnum : return "SmbAini"; case SmbBMaxEnum : return "SmbBMax"; @@ -767,4 +769,7 @@ case SmbBPosEnum : return "SmbBPos"; case SmbCEnum : return "SmbC"; + case SmbCcsnowValueEnum : return "SmbCcsnowValue"; + case SmbCciceValueEnum : return "SmbCciceValue"; + case SmbCotValueEnum : return "SmbCotValue"; case SmbDEnum : return "SmbD"; case SmbDailyairdensityEnum : return "SmbDailyairdensity"; @@ -780,4 +785,5 @@ case SmbDlwrfEnum : return "SmbDlwrf"; case SmbDswrfEnum : return "SmbDswrf"; + case SmbDswdiffrfEnum : return "SmbDswdiffrf"; case SmbDzAddEnum : return "SmbDzAdd"; case SmbDzEnum : return "SmbDz"; @@ -832,4 +838,5 @@ case SmbSmbCorrEnum : return "SmbSmbCorr"; case SmbSmbrefEnum : return "SmbSmbref"; + case SmbSzaValueEnum : return "SmbSzaValue"; case SmbTEnum : return "SmbT"; case SmbTaEnum : return "SmbTa"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 25996) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 25997) @@ -777,6 +777,8 @@ else if (strcmp(name,"SmbAccumulatedRunoff")==0) return SmbAccumulatedRunoffEnum; else if (strcmp(name,"SmbA")==0) return SmbAEnum; + else if (strcmp(name,"SmbAdiff")==0) return SmbAdiffEnum; else if (strcmp(name,"SmbAValue")==0) return SmbAValueEnum; else if (strcmp(name,"SmbAccumulation")==0) return SmbAccumulationEnum; + else if (strcmp(name,"SmbAdiffini")==0) return SmbAdiffiniEnum; else if (strcmp(name,"SmbAini")==0) return SmbAiniEnum; else if (strcmp(name,"SmbBMax")==0) return SmbBMaxEnum; @@ -785,4 +787,7 @@ else if (strcmp(name,"SmbBPos")==0) return SmbBPosEnum; else if (strcmp(name,"SmbC")==0) return SmbCEnum; + else if (strcmp(name,"SmbCcsnowValue")==0) return SmbCcsnowValueEnum; + else if (strcmp(name,"SmbCciceValue")==0) return SmbCciceValueEnum; + else if (strcmp(name,"SmbCotValue")==0) return SmbCotValueEnum; else if (strcmp(name,"SmbD")==0) return SmbDEnum; else if (strcmp(name,"SmbDailyairdensity")==0) return SmbDailyairdensityEnum; @@ -798,4 +803,5 @@ else if (strcmp(name,"SmbDlwrf")==0) return SmbDlwrfEnum; else if (strcmp(name,"SmbDswrf")==0) return SmbDswrfEnum; + else if (strcmp(name,"SmbDswdiffrf")==0) return SmbDswdiffrfEnum; else if (strcmp(name,"SmbDzAdd")==0) return SmbDzAddEnum; else if (strcmp(name,"SmbDz")==0) return SmbDzEnum; @@ -850,4 +856,5 @@ else if (strcmp(name,"SmbSmbCorr")==0) return SmbSmbCorrEnum; else if (strcmp(name,"SmbSmbref")==0) return SmbSmbrefEnum; + else if (strcmp(name,"SmbSzaValue")==0) return SmbSzaValueEnum; else if (strcmp(name,"SmbT")==0) return SmbTEnum; else if (strcmp(name,"SmbTa")==0) return SmbTaEnum; @@ -868,5 +875,8 @@ else if (strcmp(name,"SmbZMax")==0) return SmbZMaxEnum; else if (strcmp(name,"SmbZMin")==0) return SmbZMinEnum; - else if (strcmp(name,"SmbZTop")==0) return SmbZTopEnum; + else stage=8; + } + if(stage==8){ + if (strcmp(name,"SmbZTop")==0) return SmbZTopEnum; else if (strcmp(name,"SmbZY")==0) return SmbZYEnum; else if (strcmp(name,"SolidearthExternalDisplacementEastRate")==0) return SolidearthExternalDisplacementEastRateEnum; @@ -875,8 +885,5 @@ else if (strcmp(name,"SolidearthExternalGeoidRate")==0) return SolidearthExternalGeoidRateEnum; else if (strcmp(name,"SolidearthExternalBarystaticSeaLevelRate")==0) return SolidearthExternalBarystaticSeaLevelRateEnum; - else stage=8; - } - if(stage==8){ - if (strcmp(name,"StrainRateeffective")==0) return StrainRateeffectiveEnum; + else if (strcmp(name,"StrainRateeffective")==0) return StrainRateeffectiveEnum; else if (strcmp(name,"StrainRateparallel")==0) return StrainRateparallelEnum; else if (strcmp(name,"StrainRateperpendicular")==0) return StrainRateperpendicularEnum; @@ -991,5 +998,8 @@ else if (strcmp(name,"Outputdefinition52")==0) return Outputdefinition52Enum; else if (strcmp(name,"Outputdefinition53")==0) return Outputdefinition53Enum; - else if (strcmp(name,"Outputdefinition54")==0) return Outputdefinition54Enum; + else stage=9; + } + if(stage==9){ + if (strcmp(name,"Outputdefinition54")==0) return Outputdefinition54Enum; else if (strcmp(name,"Outputdefinition55")==0) return Outputdefinition55Enum; else if (strcmp(name,"Outputdefinition56")==0) return Outputdefinition56Enum; @@ -998,8 +1008,5 @@ else if (strcmp(name,"Outputdefinition59")==0) return Outputdefinition59Enum; else if (strcmp(name,"Outputdefinition5")==0) return Outputdefinition5Enum; - else stage=9; - } - if(stage==9){ - if (strcmp(name,"Outputdefinition60")==0) return Outputdefinition60Enum; + else if (strcmp(name,"Outputdefinition60")==0) return Outputdefinition60Enum; else if (strcmp(name,"Outputdefinition61")==0) return Outputdefinition61Enum; else if (strcmp(name,"Outputdefinition62")==0) return Outputdefinition62Enum; @@ -1114,5 +1121,8 @@ else if (strcmp(name,"DepthAverageAnalysis")==0) return DepthAverageAnalysisEnum; else if (strcmp(name,"DeviatoricStressErrorEstimator")==0) return DeviatoricStressErrorEstimatorEnum; - else if (strcmp(name,"Divergence")==0) return DivergenceEnum; + else stage=10; + } + if(stage==10){ + if (strcmp(name,"Divergence")==0) return DivergenceEnum; else if (strcmp(name,"Domain3Dsurface")==0) return Domain3DsurfaceEnum; else if (strcmp(name,"DoubleArrayInput")==0) return DoubleArrayInputEnum; @@ -1121,8 +1131,5 @@ else if (strcmp(name,"DoubleMatArrayParam")==0) return DoubleMatArrayParamEnum; else if (strcmp(name,"DoubleMatExternalResult")==0) return DoubleMatExternalResultEnum; - else stage=10; - } - if(stage==10){ - if (strcmp(name,"DoubleMatParam")==0) return DoubleMatParamEnum; + else if (strcmp(name,"DoubleMatParam")==0) return DoubleMatParamEnum; else if (strcmp(name,"DoubleParam")==0) return DoubleParamEnum; else if (strcmp(name,"DoubleVecParam")==0) return DoubleVecParamEnum; @@ -1237,5 +1244,8 @@ else if (strcmp(name,"MantlePlumeGeothermalFlux")==0) return MantlePlumeGeothermalFluxEnum; else if (strcmp(name,"MassFlux")==0) return MassFluxEnum; - else if (strcmp(name,"Masscon")==0) return MassconEnum; + else stage=11; + } + if(stage==11){ + if (strcmp(name,"Masscon")==0) return MassconEnum; else if (strcmp(name,"Massconaxpby")==0) return MassconaxpbyEnum; else if (strcmp(name,"Massfluxatgate")==0) return MassfluxatgateEnum; @@ -1244,8 +1254,5 @@ else if (strcmp(name,"Matdamageice")==0) return MatdamageiceEnum; else if (strcmp(name,"Matenhancedice")==0) return MatenhancediceEnum; - else stage=11; - } - if(stage==11){ - if (strcmp(name,"Materials")==0) return MaterialsEnum; + else if (strcmp(name,"Materials")==0) return MaterialsEnum; else if (strcmp(name,"Matestar")==0) return MatestarEnum; else if (strcmp(name,"Matice")==0) return MaticeEnum; @@ -1360,5 +1367,8 @@ else if (strcmp(name,"SpatialLinearFloatingMeltRate")==0) return SpatialLinearFloatingMeltRateEnum; else if (strcmp(name,"SpcDynamic")==0) return SpcDynamicEnum; - else if (strcmp(name,"SpcStatic")==0) return SpcStaticEnum; + else stage=12; + } + if(stage==12){ + if (strcmp(name,"SpcStatic")==0) return SpcStaticEnum; else if (strcmp(name,"SpcTransient")==0) return SpcTransientEnum; else if (strcmp(name,"Sset")==0) return SsetEnum; @@ -1367,8 +1377,5 @@ else if (strcmp(name,"StressIntensityFactor")==0) return StressIntensityFactorEnum; else if (strcmp(name,"StressbalanceAnalysis")==0) return StressbalanceAnalysisEnum; - else stage=12; - } - if(stage==12){ - if (strcmp(name,"StressbalanceConvergenceNumSteps")==0) return StressbalanceConvergenceNumStepsEnum; + else if (strcmp(name,"StressbalanceConvergenceNumSteps")==0) return StressbalanceConvergenceNumStepsEnum; else if (strcmp(name,"StressbalanceSIAAnalysis")==0) return StressbalanceSIAAnalysisEnum; else if (strcmp(name,"StressbalanceSolution")==0) return StressbalanceSolutionEnum; Index: /issm/trunk-jpl/src/m/classes/SMBgemb.m =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBgemb.m (revision 25996) +++ /issm/trunk-jpl/src/m/classes/SMBgemb.m (revision 25997) @@ -41,5 +41,5 @@ %optional inputs: - aValue = NaN; %Albedo forcing at every element. Used only if aIdx == 0. + aValue = NaN; %Albedo forcing at every element. Used only if aIdx == 0, or density exceeds adThresh teValue = NaN; %Outward longwave radiation thermal emissivity forcing at every element (default in code is 1) @@ -53,4 +53,5 @@ Wini = NaN; %Water content (kg m-2) Aini = NaN; %albedo (0-1) + Adiffini = NaN; %albedo, diffusive radiation (0-1) Tini = NaN; %snow temperature (K) Sizeini = NaN; %Number of layers @@ -58,27 +59,27 @@ %settings: aIdx = NaN; %method for calculating albedo and subsurface absorption (default is 1) - % 0: direct input from aValue parameter - % 1: effective grain radius [Gardner & Sharp, 2009] - % 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) - % 3: density and cloud amount [Greuell & Konzelmann, 1994] - % 4: exponential time decay & wetness [Bougamont & Bamber, 2005] + % 0: direct input from aValue parameter + % 1: effective grain radius [Gardner & Sharp, 2009] + % 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) + % 3: density and cloud amount [Greuell & Konzelmann, 1994] + % 4: exponential time decay & wetness [Bougamont & Bamber, 2005] 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)) denIdx = NaN; %densification model to use (default is 2): - % 1 = emperical model of Herron and Langway (1980) - % 2 = semi-emperical model of Anthern et al. (2010) - % 3 = DO NOT USE: physical model from Appendix B of Anthern et al. (2010) - % 4 = DO NOT USE: emperical model of Li and Zwally (2004) - % 5 = DO NOT USE: modified emperical model (4) by Helsen et al. (2008) - % 6 = Antarctica semi-emperical model of Ligtenberg et al. (2011) - % 7 = Greenland semi-emperical model of Kuipers Munneke et al. (2015) + % 1 = emperical model of Herron and Langway (1980) + % 2 = semi-emperical model of Anthern et al. (2010) + % 3 = DO NOT USE: physical model from Appendix B of Anthern et al. (2010) + % 4 = DO NOT USE: emperical model of Li and Zwally (2004) + % 5 = DO NOT USE: modified emperical model (4) by Helsen et al. (2008) + % 6 = Antarctica semi-emperical model of Ligtenberg et al. (2011) + % 7 = Greenland semi-emperical model of Kuipers Munneke et al. (2015) dsnowIdx = NaN; %model for fresh snow accumulation density (default is 1): - % 0 = Original GEMB value, 150 kg/m^3 - % 1 = Antarctica value of fresh snow density, 350 kg/m^3 - % 2 = Greenland value of fresh snow density, 315 kg/m^3, Fausto et al. (2008) - % 3 = Antarctica model of Kaspers et al. (2004) - % 4 = Greenland model of Kuipers Munneke et al. (2015) + % 0 = Original GEMB value, 150 kg/m^3 + % 1 = Antarctica value of fresh snow density, 350 kg/m^3 + % 2 = Greenland value of fresh snow density, 315 kg/m^3, Fausto et al. (2008) + % 3 = Antarctica model of Kaspers et al. (2004) + % 4 = Greenland model of Kuipers Munneke et al. (2015) zTop = NaN; % depth over which grid length is constant at the top of the snopack (default 10) [m] @@ -92,16 +93,22 @@ %specific albedo parameters: + %Method 1 + dswdiffrf = NaN; %downward diffusive shortwave radiation flux [W/m^2] + szaValue = NaN; %Solar Zenith Angle [degree] + cotValue = NaN; %Cloud Optical Thickness + ccsnowValue = NaN; %concentration of light absorbing carbon for snow [ppm1] + cciceValue = NaN; %concentration of light absorbing carbon for ice [ppm1] %Method 1 and 2: aSnow = NaN; % new snow albedo (0.64 - 0.89) aIce = NaN; % range 0.27-0.58 for old snow - %Method 3: Radiation Correction Factors -> only used for met station data and Greuell & Konzelmann, 1994 albedo + %Method 3: Radiation Correction Factors -> only used for met station data and Greuell & Konzelmann, 1994 albedo cldFrac = NaN; % average cloud amount - %Method 4: additonal tuning parameters albedo as a funtion of age and water content (Bougamont et al., 2005) + %Method 4: additonal tuning parameters albedo as a funtion of age and water content (Bougamont et al., 2005) t0wet = NaN; % time scale for wet snow (15-21.9) t0dry = NaN; % warm snow timescale (30) K = NaN; % time scale temperature coef. (7) adThresh = NaN; %Apply aIdx method to all areas with densities below this value, - %or else apply direct input value from aValue, allowing albedo to be altered. - %Default value is rho water (1023 kg m-3). + %or else apply direct input value from aValue, allowing albedo to be altered. + %Default value is rho water (1023 kg m-3). %densities: @@ -134,17 +141,63 @@ function self = extrude(self,md) % {{{ - self.Ta=project3d(md,'vector',self.Ta,'type','node'); - self.V=project3d(md,'vector',self.V,'type','node'); - self.dswrf=project3d(md,'vector',self.dswrf,'type','node'); - self.dswrf=project3d(md,'vector',self.dswrf,'type','node'); - self.P=project3d(md,'vector',self.P,'type','node'); - self.eAir=project3d(md,'vector',self.eAir,'type','node'); - self.pAir=project3d(md,'vector',self.pAir,'type','node'); - + self.Ta=project3d(md,'vector',self.Ta,'type','element'); + self.V=project3d(md,'vector',self.V,'type','element'); + self.dswrf=project3d(md,'vector',self.dswrf,'type','element'); + self.dslrf=project3d(md,'vector',self.dslrf,'type','element'); + self.P=project3d(md,'vector',self.P,'type','element'); + self.eAir=project3d(md,'vector',self.eAir,'type','element'); + self.pAir=project3d(md,'vector',self.pAir,'type','element'); + + if ~isnan(self.Dzini) + self.Dzini=project3d(md,'vector',self.Dzini,'type','element'); + end + if ~isnan(self.Dini) + self.Dini=project3d(md,'vector',self.Dini,'type','element'); + end + if ~isnan(self.Reini) + self.Reini=project3d(md,'vector',self.Reini,'type','element'); + end + if ~isnan(self.Gdnini) + self.Gdnini=project3d(md,'vector',self.Gdnini,'type','element'); + end + if ~isnan(self.Gspini) + self.Gspini=project3d(md,'vector',self.Gspini,'type','element'); + end + if ~isnan(self.ECini) + self.ECini=project3d(md,'vector',self.ECini,'type','element'); + end + if ~isnan(self.Wini) + self.Wini=project3d(md,'vector',self.Wini,'type','element'); + end + if ~isnan(self.Aini) + self.Aini=project3d(md,'vector',self.Aini,'type','element'); + end + if ~isnan(self.Adiffini) + self.Adiffini=project3d(md,'vector',self.Adiffini,'type','element'); + end + if ~isnan(self.Tini) + self.Tini=project3d(md,'vector',self.Tini,'type','element'); + end + + if ~isnan(self.dswdiffrf) + self.dswdiffrf=project3d(md,'vector',self.dswdiffrf,'type','element'); + end + if ~isnan(self.szaValue) + self.szaValue=project3d(md,'vector',self.szaValue,'type','element'); + end + if ~isnan(self.cotValue) + self.cotValue=project3d(md,'vector',self.cotValue,'type','element'); + end + if ~isnan(self.ccsnowValue) + self.ccsnowValue=project3d(md,'vector',self.ccsnowValue,'type','element'); + end + if ~isnan(self.cciceValue) + self.cciceValue=project3d(md,'vector',self.cciceValue,'type','element'); + end if (aIdx == 0) & ~isnan(self.aValue) - self.aValue=project3d(md,'vector',self.aValue,'type','node'); + self.aValue=project3d(md,'vector',self.aValue,'type','element'); end if ~isnan(self.teValue) - self.teValue=project3d(md,'vector',self.teValue,'type','node'); + self.teValue=project3d(md,'vector',self.teValue,'type','element'); end @@ -152,5 +205,5 @@ end % }}} function list = defaultoutputs(self,md) % {{{ - list = {'SmbMassBalance'}; + list = {'SmbMassBalance','SmbAccumulatedMassBalance'}; end % }}} function self = setdefaultparameters(self,mesh,geometry) % {{{ @@ -196,4 +249,10 @@ self.aValue = self.aSnow*ones(mesh.numberofelements,1); + self.dswdiffrf=0.0*ones(mesh.numberofelements,1); + self.szaValue=0.0*ones(mesh.numberofelements,1); + self.cotValue=0.0*ones(mesh.numberofelements,1); + self.ccsnowValue=0.0*ones(mesh.numberofelements,1); + self.cciceValue=0.0*ones(mesh.numberofelements,1); + self.Dzini=0.05*ones(mesh.numberofelements,2); self.Dini=910.0*ones(mesh.numberofelements,2); @@ -204,4 +263,5 @@ self.Wini=0.0*ones(mesh.numberofelements,2); self.Aini=self.aSnow*ones(mesh.numberofelements,2); + self.Adiffini=ones(mesh.numberofelements,2); self.Tini=273.15*ones(mesh.numberofelements,2); % /!\ Default value of Tini must be equal to Tmean but don't know Tmean yet (computed when atmospheric forcings are interpolated on mesh). @@ -226,4 +286,5 @@ md = checkfield(md,'fieldname','smb.V','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<',45,'size',size(self.Ta)); %max 500 km/h md = checkfield(md,'fieldname','smb.dswrf','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<=',1400,'size',size(self.Ta)); + md = checkfield(md,'fieldname','smb.dswdiffrf','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<=',1400); md = checkfield(md,'fieldname','smb.dlwrf','timeseries',1,'NaN',1,'Inf',1,'>=',0,'size',size(self.Ta)); md = checkfield(md,'fieldname','smb.P','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<=',200,'size',size(self.Ta)); @@ -232,15 +293,15 @@ if (self.isclimatology) md = checkfield(md,'fieldname', 'smb.Ta', 'size',[md.mesh.numberofelements+1],... - 'message',['Ta must have md.mesh.numberofelements+1 rows in order to force a climatology']); + 'message',['Ta must have md.mesh.numberofelements+1 rows in order to force a climatology']); md = checkfield(md,'fieldname', 'smb.V', 'size',[md.mesh.numberofelements+1],... - 'message',['V must have md.mesh.numberofelements+1 rows in order to force a climatology']); + 'message',['V must have md.mesh.numberofelements+1 rows in order to force a climatology']); md = checkfield(md,'fieldname', 'smb.dswrf', 'size',[md.mesh.numberofelements+1],... - 'message',['dswrf must have md.mesh.numberofelements+1 rows in order to force a climatology']); + 'message',['dswrf must have md.mesh.numberofelements+1 rows in order to force a climatology']); md = checkfield(md,'fieldname', 'smb.dlwrf', 'size',[md.mesh.numberofelements+1],... - 'message',['dlwrf must have md.mesh.numberofelements+1 rows in order to force a climatology']); + 'message',['dlwrf must have md.mesh.numberofelements+1 rows in order to force a climatology']); md = checkfield(md,'fieldname', 'smb.P', 'size',[md.mesh.numberofelements+1],... - 'message',['P must have md.mesh.numberofelements+1 rows in order to force a climatology']); + 'message',['P must have md.mesh.numberofelements+1 rows in order to force a climatology']); md = checkfield(md,'fieldname', 'smb.eAir', 'size',[md.mesh.numberofelements+1],... - 'message',['eAir must have md.mesh.numberofelements+1 rows in order to force a climatology']); + 'message',['eAir must have md.mesh.numberofelements+1 rows in order to force a climatology']); end @@ -273,4 +334,10 @@ md = checkfield(md,'fieldname','smb.aSnow','NaN',1,'Inf',1,'>=',.64,'<=',.89); md = checkfield(md,'fieldname','smb.aIce','NaN',1,'Inf',1,'>=',.27,'<=',.58); + if self.aIdx==1 + md = checkfield(md,'fieldname','smb.szaValue','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<=',90); + md = checkfield(md,'fieldname','smb.cotValue','timeseries',1,'NaN',1,'Inf',1,'>=',0); + md = checkfield(md,'fieldname','smb.ccsnowValue','timeseries',1,'NaN',1,'Inf',1,'>=',0); + md = checkfield(md,'fieldname','smb.cciceValue','timeseries',1,'NaN',1,'Inf',1,'>=',0); + end case 3 md = checkfield(md,'fieldname','smb.cldFrac','NaN',1,'Inf',1,'>=',0,'<=',1); @@ -307,6 +374,7 @@ fielddisplay(self,'Ta','2 m air temperature, in Kelvin'); fielddisplay(self,'V','wind speed (m s-1)'); - fielddisplay(self,'dlwrf','downward shortwave radiation flux [W/m^2]'); - fielddisplay(self,'dswrf','downward longwave radiation flux [W/m^2]'); + fielddisplay(self,'dswrf','downward shortwave radiation flux [W/m^2]'); + fielddisplay(self,'dswdiffrf','downward diffusive portion of shortwave radiation flux (default to 0) [W/m^2]'); + fielddisplay(self,'dlwrf','downward longwave radiation flux [W/m^2]'); fielddisplay(self,'P','precipitation [mm w.e. / m^2]'); fielddisplay(self,'eAir','screen level vapor pressure [Pa]'); @@ -328,9 +396,9 @@ 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.'}); fielddisplay(self,'aIdx',{'method for calculating albedo and subsurface absorption (default is 1)',... - '0: direct input from aValue parameter',... - '1: effective grain radius [Gardner & Sharp, 2009]',... - '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)',... - '3: density and cloud amount [Greuell & Konzelmann, 1994]',... - '4: exponential time decay & wetness [Bougamont & Bamber, 2005]'}) + '0: direct input from aValue parameter',... + '1: effective grain radius [Gardner & Sharp, 2009]',... + '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)',... + '3: density and cloud amount [Greuell & Konzelmann, 1994]',... + '4: exponential time decay & wetness [Bougamont & Bamber, 2005]'}) fielddisplay(self,'teValue','Outward longwave radiation thermal emissivity forcing at every element (default in code is 1)'); @@ -345,13 +413,20 @@ fielddisplay(self,'Wini','Initial snow water content when restart [kg m-2]'); fielddisplay(self,'Aini','Initial albedo when restart [-]'); + fielddisplay(self,'Adiffini','Initial diffusive radiation albedo when restart (default to 1) [-]'); fielddisplay(self,'Tini','Initial snow temperature when restart [K]'); fielddisplay(self,'Sizeini','Initial number of layers when restart [-]'); %additional albedo parameters: - fielddisplay(self,'aValue','Albedo forcing at every element.'); + fielddisplay(self,'aValue','Albedo forcing at every element'); switch self.aIdx case {1 2} fielddisplay(self,'aSnow','new snow albedo (0.64 - 0.89)'); fielddisplay(self,'aIce','albedo of ice (0.27-0.58)'); + if self.aIdx==1 + fielddisplay(self,'szaValue','Solar Zenith Angle [degree]'); + fielddisplay(self,'cotValue','Cloud Optical Thickness'); + fielddisplay(self,'ccsnowValue','concentration of light absorbing carbon for snow [ppm1]'); + fielddisplay(self,'cciceValue','concentration of light absorbing carbon for snow [ppm1]'); + end case 3 fielddisplay(self,'cldFrac','average cloud amount'); @@ -364,17 +439,17 @@ 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)]'); fielddisplay(self,'denIdx',{'densification model to use (default is 2):',... - '1 = emperical model of Herron and Langway (1980)',... - '2 = semi-emperical model of Anthern et al. (2010)',... - '3 = DO NOT USE: physical model from Appendix B of Anthern et al. (2010)',... - '4 = DO NOT USE: emperical model of Li and Zwally (2004)',... - '5 = DO NOT USE: modified emperical model (4) by Helsen et al. (2008)',... - '6 = Antarctica semi-emperical model of Ligtenberg et al. (2011)',... - '7 = Greenland semi-emperical model of Kuipers Munneke et al. (2015)'}); + '1 = emperical model of Herron and Langway (1980)',... + '2 = semi-emperical model of Anthern et al. (2010)',... + '3 = DO NOT USE: physical model from Appendix B of Anthern et al. (2010)',... + '4 = DO NOT USE: emperical model of Li and Zwally (2004)',... + '5 = DO NOT USE: modified emperical model (4) by Helsen et al. (2008)',... + '6 = Antarctica semi-emperical model of Ligtenberg et al. (2011)',... + '7 = Greenland semi-emperical model of Kuipers Munneke et al. (2015)'}); fielddisplay(self,'dsnowIdx',{'model for fresh snow accumulation density (default is 1):',... - '0 = Original GEMB value, 150 kg/m^3',... - '1 = Antarctica value of fresh snow density, 350 kg/m^3',... - '2 = Greenland value of fresh snow density, 315 kg/m^3, Fausto et al. (2008)',... - '3 = Antarctica model of Kaspers et al. (2004), Make sure to set Vmean accurately',... - '4 = Greenland model of Kuipers Munneke et al. (2015)'}); + '0 = Original GEMB value, 150 kg/m^3',... + '1 = Antarctica value of fresh snow density, 350 kg/m^3',... + '2 = Greenland value of fresh snow density, 315 kg/m^3, Fausto et al. (2008)',... + '3 = Antarctica model of Kaspers et al. (2004), Make sure to set Vmean accurately',... + '4 = Greenland model of Kuipers Munneke et al. (2015)'}); fielddisplay(self, 'steps_per_step', 'number of smb steps per time step'); @@ -407,4 +482,5 @@ WriteData(fid,prefix,'object',self,'class','smb','fieldname','V','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts); WriteData(fid,prefix,'object',self,'class','smb','fieldname','dswrf','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','dswdiffrf','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts); WriteData(fid,prefix,'object',self,'class','smb','fieldname','dlwrf','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts); WriteData(fid,prefix,'object',self,'class','smb','fieldname','P','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts); @@ -441,4 +517,8 @@ WriteData(fid,prefix,'object',self,'class','smb','fieldname','aValue','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts); WriteData(fid,prefix,'object',self,'class','smb','fieldname','teValue','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','szaValue','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','cotValue','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','ccsnowValue','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','cciceValue','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts); %snow properties init @@ -451,4 +531,5 @@ WriteData(fid,prefix,'object',self,'class','smb','fieldname','Wini','format','DoubleMat','mattype',3); WriteData(fid,prefix,'object',self,'class','smb','fieldname','Aini','format','DoubleMat','mattype',3); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','Adiffini','format','DoubleMat','mattype',3); WriteData(fid,prefix,'object',self,'class','smb','fieldname','Tini','format','DoubleMat','mattype',3); WriteData(fid,prefix,'object',self,'class','smb','fieldname','Sizeini','format','IntMat','mattype',2); @@ -464,4 +545,25 @@ error('All GEMB forcings (Ta, P, V, dswrf, dlwrf, eAir, pAir) must have the same time steps in the final row!'); end + if size(md.smb.teValue,2)>1 & any(md.smb.teValue(end,:) - md.smb.Ta(end,:) ~= 0) + error('If GEMB forcing teValue is transient, it must have the same time steps as input Ta in the final row!'); + end + if size(md.smb.dswdiffrf,2)>1 & any(md.smb.dswdiffrf(end,:) - md.smb.Ta(end,:) ~= 0) + error('If GEMB forcing dswdiffrf is transient, it must have the same time steps as input Ta in the final row!'); + end + if size(md.smb.aValue,2)>1 & any(md.smb.aValue(end,:) - md.smb.Ta(end,:) ~= 0) + error('If GEMB forcing aValue is transient, it must have the same time steps as input Ta in the final row!'); + end + if size(md.smb.szaValue,2)>1 & any(md.smb.szaValue(end,:) - md.smb.Ta(end,:) ~= 0) + error('If GEMB forcing szaValue is transient, it must have the same time steps as input Ta in the final row!'); + end + if size(md.smb.cotValue,2)>1 & any(md.smb.cotValue(end,:) - md.smb.Ta(end,:) ~= 0) + error('If GEMB forcing cotValue is transient, it must have the same time steps as input Ta in the final row!'); + end + if size(md.smb.ccsnowValue,2)>1 & any(md.smb.ccsnowValue(end,:) - md.smb.Ta(end,:) ~= 0) + error('If GEMB forcing ccsnowValue is transient, it must have the same time steps as input Ta in the final row!'); + end + if size(md.smb.cciceValue,2)>1 & any(md.smb.cciceValue(end,:) - md.smb.Ta(end,:) ~= 0) + error('If GEMB forcing cciceValue is transient, it must have the same time steps as input Ta in the final row!'); + end time=self.Ta(end,:); %assume all forcings are on the same time step dtime=diff(time,1); Index: /issm/trunk-jpl/src/m/classes/SMBgemb.py =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBgemb.py (revision 25996) +++ /issm/trunk-jpl/src/m/classes/SMBgemb.py (revision 25997) @@ -47,5 +47,5 @@ #optional inputs: - self.aValue = np.nan # Albedo forcing at every element. Used only if aIdx == 0. + self.aValue = np.nan # Albedo forcing at every element. Used only if aIdx == 0, or density exceeds adThresh self.teValue = np.nan # Outward longwave radiation thermal emissivity forcing at every element (default in code is 1) @@ -59,4 +59,5 @@ self.Wini = np.nan # Water content (kg m-2) self.Aini = np.nan # albedo (0-1) + self.Adiffini = np.nan # albedo, diffusive radiation (0-1) self.Tini = np.nan # snow temperature (K) self.Sizeini = np.nan # Number of layers @@ -96,4 +97,10 @@ #specific albedo parameters: + #Method 1 + dswdiffrf = np.nan #downward diffusive shortwave radiation flux [W/m^2] + szaValue = np.nan #Solar Zenith Angle [degree] + cotValue = np.nan #Cloud Optical Thickness + ccsnowValue = np.nan #concentration of light absorbing carbon for snow [ppm1] + cciceValue = np.nan #concentration of light absorbing carbon for ice [ppm1] #Method 1 and 2: self.aSnow = np.nan # new snow albedo (0.64 - 0.89) @@ -152,6 +159,7 @@ string = "%s\n%s" % (string, fielddisplay(self, 'Ta', '2 m air temperature, in Kelvin')) string = "%s\n%s" % (string, fielddisplay(self, 'V', 'wind speed (m s-1)')) - string = "%s\n%s" % (string, fielddisplay(self, 'dlwrf', 'downward shortwave radiation flux [W/m^2]')) - string = "%s\n%s" % (string, fielddisplay(self, 'dswrf', 'downward longwave radiation flux [W/m^2]')) + string = "%s\n%s" % (string, fielddisplay(self, 'dswrf', 'downward shortwave radiation flux [W/m^2]')) + string = "%s\n%s" % (string, fielddisplay(self, 'dswdiffrf', 'downward diffusive portion of shortwave radiation flux (default to 0) [W/m^2]')) + string = "%s\n%s" % (string, fielddisplay(self, 'dlwrf', 'downward longwave radiation flux [W/m^2]')) string = "%s\n%s" % (string, fielddisplay(self, 'P', 'precipitation [mm w.e. / m^2]')) string = "%s\n%s" % (string, fielddisplay(self, 'eAir', 'screen level vapor pressure [Pa]')) @@ -188,12 +196,19 @@ string = "%s\n%s" % (string, fielddisplay(self, 'Wini', 'Initial snow water content when restart [kg m-2]')) string = "%s\n%s" % (string, fielddisplay(self, 'Aini', 'Initial albedo when restart [-]')) + string = "%s\n%s" % (string, fielddisplay(self, 'Adiffini', 'Initial diffusive radiation albedo when restart (default to 1) [-]')) string = "%s\n%s" % (string, fielddisplay(self, 'Tini', 'Initial snow temperature when restart [K]')) string = "%s\n%s" % (string, fielddisplay(self, 'Sizeini', 'Initial number of layers when restart [-]')) #additional albedo parameters: - string = "%s\n%s" % (string, fielddisplay(self, 'aValue', 'Albedo forcing at every element.')) + string = "%s\n%s" % (string, fielddisplay(self, 'aValue', 'Albedo forcing at every element')) 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)): string = "%s\n%s" % (string, fielddisplay(self, 'aSnow', 'new snow albedo (0.64 - 0.89)')) string = "%s\n%s" % (string, fielddisplay(self, 'aIce', 'albedo of ice (0.27-0.58)')) + if self.aIdx == 1: + string = "%s\n%s" % (string, fielddisplay(self,'szaValue','Solar Zenith Angle [degree]')) + string = "%s\n%s" % (string, fielddisplay(self,'cotValue','Cloud Optical Thickness')) + string = "%s\n%s" % (string, fielddisplay(self,'ccsnowValue','concentration of light absorbing carbon for snow [ppm1]')) + string = "%s\n%s" % (string, fielddisplay(self,'cciceValue','concentration of light absorbing carbon for snow [ppm1]')) + elif self.aIdx == 3: string = "%s\n%s" % (string, fielddisplay(self, 'cldFrac', 'average cloud amount')) @@ -228,16 +243,47 @@ def extrude(self, md): # {{{ - self.Ta = project3d(md, 'vector', self.Ta, 'type', 'node') - self.V = project3d(md, 'vector', self.V, 'type', 'node') - self.dswrf = project3d(md, 'vector', self.dswrf, 'type', 'node') - self.dswrf = project3d(md, 'vector', self.dswrf, 'type', 'node') - self.P = project3d(md, 'vector', self.P, 'type', 'node') - self.eAir = project3d(md, 'vector', self.eAir, 'type', 'node') - self.pAir = project3d(md, 'vector', self.pAir, 'type', 'node') - - if (self.aIdx == 0) and np.isnan(self.aValue): - self.aValue = project3d(md, 'vector', self.aValue, 'type', 'node') - if np.isnan(self.teValue): - self.teValue = project3d(md, 'vector', self.teValue, 'type', 'node') + self.Ta = project3d(md, 'vector', self.Ta, 'type', 'element') + self.V = project3d(md, 'vector', self.V, 'type', 'element') + self.dswrf = project3d(md, 'vector', self.dswrf, 'type', 'element') + self.dswrf = project3d(md, 'vector', self.dswrf, 'type', 'element') + self.P = project3d(md, 'vector', self.P, 'type', 'element') + self.eAir = project3d(md, 'vector', self.eAir, 'type', 'element') + self.pAir = project3d(md, 'vector', self.pAir, 'type', 'element') + + if not np.isnan(self.Dzini): + self.self.Dzini=project3d(md,'vector',self.self.Dzini,'type','element'); + if not np.isnan(self.Dini): + self.self.Dini=project3d(md,'vector',self.Dini,'type','element'); + if not np.isnan(self.Reini): + self.self.Reini=project3d(md,'vector',self.Reini,'type','element'); + if not np.isnan(self.Gdnini): + self.Gdnini=project3d(md,'vector',self.Gdnini,'type','element'); + if not np.isnan(self.Gspini): + self.Gspini=project3d(md,'vector',self.Gspini,'type','element'); + if not np.isnan(self.ECini): + self.ECini=project3d(md,'vector',self.ECini,'type','element'); + if not np.isnan(self.Wini): + self.Wini=project3d(md,'vector',self.Wini,'type','element'); + if not np.isnan(self.Aini): + self.Aini=project3d(md,'vector',self.Aini,'type','element'); + if not np.isnan(self.Adiffini): + self.Adiffini=project3d(md,'vector',self.Adiffini,'type','element'); + if not np.isnan(self.Tini): + self.Tini=project3d(md,'vector',self.Tini,'type','element'); + if not np.isnan(self.dswdiffrf): + self.dswdiffrf=project3d(md,'vector',self.dswdiffrf,'type','element'); + if not np.isnan(self.szaValue): + self.szaValue=project3d(md,'vector',self.szaValue,'type','element'); + if not np.isnan(self.cotValue): + self.cotValue=project3d(md,'vector',self.cotValue,'type','element'); + if not np.isnan(self.ccsnowValue): + self.ccsnowValue=project3d(md,'vector',self.ccsnowValue,'type','element'); + if not np.isnan(self.cciceValue): + self.cciceValue=project3d(md,'vector',self.cciceValue,'type','element'); + + if (self.aIdx == 0) and (not np.isnan(self.aValue)): + self.aValue = project3d(md, 'vector', self.aValue, 'type', 'element') + if not np.isnan(self.teValue): + self.teValue = project3d(md, 'vector', self.teValue, 'type', 'element') return self @@ -245,5 +291,5 @@ def defaultoutputs(self, md): # {{{ - return ['SmbMassBalance'] + return ['SmbMassBalance','SmbAccumulatedMassBalance'] #}}} @@ -290,4 +336,10 @@ self.aValue = self.aSnow * np.ones(mesh.numberofelements,) + self.dswdiffrf = 0.0 * np.ones(mesh.numberofelements,) + self.szaValue = 0.0 * np.ones(mesh.numberofelements,) + self.cotValue = 0.0 * np.ones(mesh.numberofelements,) + self.ccsnowValue = 0.0 * np.ones(mesh.numberofelements,) + self.cciceValue = 0.0 * np.ones(mesh.numberofelements,) + self.Dzini = 0.05 * np.ones((mesh.numberofelements, 2)) self.Dini = 910.0 * np.ones((mesh.numberofelements, 2)) @@ -298,4 +350,5 @@ self.Wini = 0.0 * np.ones((mesh.numberofelements, 2)) self.Aini = self.aSnow * np.ones((mesh.numberofelements, 2)) + self.Adiffini = np.ones((mesh.numberofelements, 2)) self.Tini = 273.15 * np.ones((mesh.numberofelements, 2)) # /!\ Default value of Tini must be equal to Tmean but don't know Tmean yet (computed when atmospheric forcings are interpolated on mesh). @@ -321,4 +374,5 @@ md = checkfield(md, 'fieldname', 'smb.V', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '> = ', 0, '<', 45, 'size', np.shape(self.Ta)) #max 500 km/h md = checkfield(md, 'fieldname', 'smb.dswrf', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '> = ', 0, '< = ', 1400, 'size', np.shape(self.Ta)) + md = checkfield(md, 'fieldname', 'smb.dswdiffrf', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '> = ', 0, '< = ', 1400) md = checkfield(md, 'fieldname', 'smb.dlwrf', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '> = ', 0, 'size', np.shape(self.Ta)) md = checkfield(md, 'fieldname', 'smb.P', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '> = ', 0, '< = ', 200, 'size', np.shape(self.Ta)) @@ -358,4 +412,10 @@ md = checkfield(md, 'fieldname', 'smb.aSnow', 'NaN', 1, 'Inf', 1, '> = ', .64, '< = ', .89) md = checkfield(md, 'fieldname', 'smb.aIce', 'NaN', 1, 'Inf', 1, '> = ', .27, '< = ', .58) + if self.aIdx == 1: + md = checkfield(md,'fieldname','smb.szaValue','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<=',90) + md = checkfield(md,'fieldname','smb.cotValue','timeseries',1,'NaN',1,'Inf',1,'>=',0) + md = checkfield(md,'fieldname','smb.ccsnowValue','timeseries',1,'NaN',1,'Inf',1,'>=',0) + md = checkfield(md,'fieldname','smb.cciceValue','timeseries',1,'NaN',1,'Inf',1,'>=',0) + elif self.aIdx == 0: md = checkfield(md, 'fieldname', 'smb.aValue', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '>=', 0, '<=', 1) @@ -396,4 +456,5 @@ WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'V', 'format', 'DoubleMat', 'mattype', 2, 'timeserieslength', md.mesh.numberofelements + 1, 'yts', yts) WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'dswrf', 'format', 'DoubleMat', 'mattype', 2, 'timeserieslength', md.mesh.numberofelements + 1, 'yts', yts) + WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'dswdiffrf', 'format', 'DoubleMat', 'mattype', 2, 'timeserieslength', md.mesh.numberofelements + 1, 'yts', yts) WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'dlwrf', 'format', 'DoubleMat', 'mattype', 2, 'timeserieslength', md.mesh.numberofelements + 1, 'yts', yts) WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'P', 'format', 'DoubleMat', 'mattype', 2, 'timeserieslength', md.mesh.numberofelements + 1, 'yts', yts) @@ -430,4 +491,8 @@ WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'aValue', 'format', 'DoubleMat', 'mattype', 2, 'timeserieslength', md.mesh.numberofelements + 1, 'yts', yts) WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'teValue', 'format', 'DoubleMat', 'mattype', 2, 'timeserieslength', md.mesh.numberofelements + 1, 'yts', yts) + WriteData(fid,prefix,'object',self,'class','smb','fieldname','szaValue','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts) + WriteData(fid,prefix,'object',self,'class','smb','fieldname','cotValue','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts) + WriteData(fid,prefix,'object',self,'class','smb','fieldname','ccsnowValue','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts) + WriteData(fid,prefix,'object',self,'class','smb','fieldname','cciceValue','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts) #snow properties init @@ -440,4 +505,5 @@ WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'Wini', 'format', 'DoubleMat', 'mattype', 3) WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'Aini', 'format', 'DoubleMat', 'mattype', 3) + WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'Adiffini', 'format', 'DoubleMat', 'mattype', 3) WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'Tini', 'format', 'DoubleMat', 'mattype', 3) WriteData(fid, prefix, 'object', self, 'class', 'smb', 'fieldname', 'Sizeini', 'format', 'IntMat', 'mattype', 2) @@ -448,4 +514,19 @@ raise IOError('All GEMB forcings (Ta, P, V, dswrf, dlwrf, eAir, pAir) must have the same time steps in the final row!') + if ((np.ndim(self.teValue)>1) & np.any(self.teValue[-1] - self.Ta[-1] != 0)): + raise IOError('If GEMB forcing teValue is transient, it must have the same time steps as input Ta in the final row!') + if ((np.ndim(self.dswdiffrf)>1) & np.any(self.dswdiffrf[-1] - self.Ta[-1] != 0)): + raise IOError('If GEMB forcing dswdiffrf is transient, it must have the same time steps as input Ta in the final row!') + if ((np.ndim(self.aValue)>1) & np.any(self.aValue[-1] - self.Ta[-1] != 0)): + raise IOError('If GEMB forcing aValue is transient, it must have the same time steps as input Ta in the final row!') + if ((np.ndim(self.szaValue)>1) & np.any(self.szaValue[-1] - self.Ta[-1] != 0)): + raise IOError('If GEMB forcing szaValue is transient, it must have the same time steps as input Ta in the final row!') + if ((np.ndim(self.cotValue)>1) & np.any(self.cotValue[-1] - self.Ta[-1] != 0)): + raise IOError('If GEMB forcing cotValue is transient, it must have the same time steps as input Ta in the final row!') + if ((np.ndim(self.ccsnowValue)>1) & np.any(self.ccsnowValue[-1] - self.Ta[-1] != 0)): + raise IOError('If GEMB forcing ccsnowValue is transient, it must have the same time steps as input Ta in the final row!') + if ((np.ndim(self.cciceValue)>1) & np.any(self.cciceValue[-1] - self.Ta[-1] != 0)): + raise IOError('If GEMB forcing cciceValue is transient, it must have the same time steps as input Ta in the final row!') + time = self.Ta[-1] #assume all forcings are on the same time step dtime = np.diff(time, n=1, axis=0)