Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 24682) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 24683) @@ -562,5 +562,5 @@ this->parameters->FindParam(&isPrecipScaled,SmbIsprecipscaledEnum); - /*Recover present day temperature and precipitation*/ + /*Recover present day temperature and precipitation*/ DatasetInput2 *dinput3 = NULL; DatasetInput2 *dinput4 = NULL; @@ -3588,5 +3588,5 @@ } /*}}}*/ -void Element::SmbGemb(){/*{{{*/ +void Element::SmbGemb(IssmDouble timeinputs, int count){/*{{{*/ /*only compute SMB at the surface: */ @@ -3605,7 +3605,4 @@ IssmDouble C=0.0; IssmDouble Tz,Vz=0.0; - IssmDouble time,dt,starttime,finaltime; - IssmDouble timeclim=0.0; - IssmDouble t,smb_dt; IssmDouble yts; IssmDouble Ta=0.0; @@ -3618,4 +3615,5 @@ IssmDouble teValue=1.0; IssmDouble aValue=0.0; + IssmDouble dt,time,smb_dt; int aIdx=0; int denIdx=0; @@ -3627,11 +3625,10 @@ IssmDouble dayEC=0.0; IssmDouble initMass=0.0; - IssmDouble sumR=0.0; - IssmDouble sumM=0.0; - IssmDouble sumEC=0.0; - IssmDouble sumP=0.0; - IssmDouble sumW=0.0; - IssmDouble sumMassAdd=0.0; - IssmDouble sumdz_add=0.0; + IssmDouble sumR=0.0; + IssmDouble sumM=0.0; + IssmDouble sumEC=0.0; + IssmDouble sumP=0.0; + IssmDouble sumW=0.0; + IssmDouble sumMassAdd=0.0; IssmDouble fac=0.0; IssmDouble sumMass=0.0; @@ -3642,6 +3639,4 @@ IssmDouble thermo_scaling=1.0; IssmDouble adThresh=1023.0; - int offsetend=-1; - IssmDouble time0, timeend, delta; /*}}}*/ /*Output variables:{{{ */ @@ -3676,5 +3671,4 @@ IssmDouble* Tini = NULL; int m=0; - int count=0; /*}}}*/ @@ -3687,6 +3681,4 @@ parameters->FindParam(&dt,TimesteppingTimeStepEnum); /*transient core time step*/ parameters->FindParam(&yts,ConstantsYtsEnum); - parameters->FindParam(&finaltime,TimesteppingFinalTimeEnum); - parameters->FindParam(&starttime,TimesteppingStartTimeEnum); parameters->FindParam(&smb_dt,SmbDtEnum); /*time period for the smb solution, usually smaller than the glaciological dt*/ parameters->FindParam(&aIdx,SmbAIdxEnum); @@ -3698,5 +3690,4 @@ parameters->FindParam(&t0dry,SmbT0dryEnum); parameters->FindParam(&K,SmbKEnum); - parameters->FindParam(&isclimatology,SmbIsclimatologyEnum); parameters->FindParam(&isgraingrowth,SmbIsgraingrowthEnum); parameters->FindParam(&isalbedo,SmbIsalbedoEnum); @@ -3723,6 +3714,4 @@ Input2 *Tz_input = this->GetInput2(SmbTzEnum); _assert_(Tz_input); Input2 *Vz_input = this->GetInput2(SmbVzEnum); _assert_(Vz_input); - Input2 *teValue_input = this->GetInput2(SmbTeValueEnum); _assert_(teValue_input); - Input2 *aValue_input = this->GetInput2(SmbAValueEnum); _assert_(aValue_input); Input2 *EC_input = NULL; @@ -3742,6 +3731,4 @@ Tz_input->GetInputValue(&Tz,gauss); Vz_input->GetInputValue(&Vz,gauss); - teValue_input->GetInputValue(&teValue,gauss); - aValue_input->GetInputValue(&aValue,gauss); /*}}}*/ @@ -3763,5 +3750,5 @@ EC_input->GetInputAverage(&EC); - /*Retrive the correct value of m (without the zeroes at the end)*/ + /*Retrieve the correct value of m (without the zeroes at the end)*/ this->GetInput2Value(&m,SmbSizeiniEnum); @@ -3816,7 +3803,6 @@ this->inputs2->GetArray(SmbAEnum,this->lid,&a,&m); this->inputs2->GetArray(SmbTEnum,this->lid,&T,&m); - EC_input = this->GetInput2(SmbECEnum); _assert_(EC_input); + EC_input = this->GetInput2(SmbECDtEnum); _assert_(EC_input); EC_input->GetInputAverage(&EC); - EC=EC*dt*rho_ice; //fixed lower temperature bounday condition - T is fixed @@ -3830,149 +3816,164 @@ // initialize cumulative variables sumR = 0; sumM = 0; sumEC = 0; sumP = 0; sumMassAdd = 0; - sumdz_add=0; //before starting loop, realize that the transient core runs this smb_core at time = time +deltaT. - //go back to time - deltaT: - time-=dt; - - timeclim=time; - if (isclimatology){ - //If this is a climatology, we need to repeat the forcing after the final time - TransientInput2* Ta_input_tr = this->inputs2->GetTransientInput(SmbTaEnum); _assert_(Ta_input_tr); - offsetend = Ta_input_tr->GetTimeInputOffset(finaltime); - time0 = Ta_input_tr->GetTimeByOffset(-1); - timeend = Ta_input_tr->GetTimeByOffset(offsetend); - if (time>time0 & timeend>time0){ - delta=(time-time0) - (timeend-time0)*(reCast((time-time0)/(timeend-time0))); - timeclim=time0+delta; - } - } - - /*Start loop: */ - count=1; - for (t=time;tSid()==0 && IssmComm::GetRank()==0)_printf0_("Time: t=" << setprecision(8) << t/365.0/24.0/3600.0 << " yr/" << (time+dt)/365.0/24.0/3600.0 << " yr" << setprecision(3) << " Step: " << count << "\n"); - - Input2 *Ta_input = this->GetInput2(SmbTaEnum,t-time+timeclim); _assert_(Ta_input); - Input2 *V_input = this->GetInput2(SmbVEnum,t-time+timeclim); _assert_(V_input); - Input2 *Dlwr_input= this->GetInput2(SmbDlwrfEnum,t-time+timeclim); _assert_(Dlwr_input); - Input2 *Dswr_input= this->GetInput2(SmbDswrfEnum,t-time+timeclim); _assert_(Dswr_input); - Input2 *P_input = this->GetInput2(SmbPEnum,t-time+timeclim); _assert_(P_input); - Input2 *eAir_input= this->GetInput2(SmbEAirEnum,t-time+timeclim); _assert_(eAir_input); - Input2 *pAir_input= this->GetInput2(SmbPAirEnum,t-time+timeclim); _assert_(pAir_input); - - /*extract daily data:{{{*/ - Ta_input->GetInputValue(&Ta,gauss);//screen level air temperature [K] - V_input->GetInputValue(&V,gauss); //wind speed [m s-1] - Dlwr_input->GetInputValue(&dlw,gauss); //downward longwave radiation flux [W m-2] - Dswr_input->GetInputValue(&dsw,gauss); //downward shortwave radiation flux [W m-2] - P_input->GetInputValue(&P,gauss); //precipitation [kg m-2] - eAir_input->GetInputValue(&eAir,gauss); //screen level vapor pressure [Pa] - pAir_input->GetInputValue(&pAir,gauss); // screen level air pressure [Pa] - teValue_input->GetInputValue(&teValue,gauss); // Emissivity [0-1] - aValue_input->GetInputValue(&aValue,gauss); // Albedo [0 1] - //_printf_("Time: " << t << " Ta: " << Ta << " V: " << V << " dlw: " << dlw << " dsw: " << dsw << " P: " << P << " eAir: " << eAir << " pAir: " << pAir << "\n"); - /*}}}*/ - - /*Snow grain metamorphism:*/ - if(isgraingrowth)grainGrowth(&re, &gdn, &gsp, T, dz, d, W, smb_dt, m, aIdx,this->Sid()); - - /*Snow, firn and ice albedo:*/ - if(isalbedo)albedo(&a,aIdx,re,d,cldFrac,aIce,aSnow,aValue,adThresh,T,W,P,EC,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()); - - /*Calculate net shortwave [W m-2]*/ - netSW = netSW + cellsum(swf,m)*smb_dt/dt; - - /*Thermal profile computation:*/ - if(isthermal)thermo(&EC, &T, &ulw, dz, d, swf, dlw, Ta, V, eAir, pAir, teValue, W[0], smb_dt, m, Vz, Tz, thermo_scaling,rho_ice,this->Sid()); - - /*Change in thickness of top cell due to evaporation/condensation assuming same density as top cell. - * need to fix this in case all or more of cell evaporates */ - dz[0] = dz[0] + EC / d[0]; - - /*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()); - - /*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, &R, &mAdd, &dz_add, &T, &d, &dz, &W, &a, &re, &gdn, &gsp, &m, dzMin, zMax, zMin, zTop,rho_ice,this->Sid()); - - /*Allow non-melt densification and determine compaction [m]*/ - if(isdensification)densification(&d,&dz, T, re, denIdx, C, smb_dt, Tmean,rho_ice,m,this->Sid()); - - /*Calculate upward longwave radiation flux [W m-2] not used in energy balance. Calculated for every - * sub-time step in thermo equations*/ - //ulw = 5.67E-8 * pow(T[0],4.0) * teValue; // + deltatest here - - /*Calculate net longwave [W m-2]*/ - meanULW = meanULW + ulw*smb_dt/dt; - netLW = netLW + (dlw - ulw)*smb_dt/dt; - - /*Calculate turbulent heat fluxes [W m-2]*/ - if(isturbulentflux)turbulentFlux(&shf, &lhf, &dayEC, Ta, T[0], V, eAir, pAir, d[0], W[0], Vz, Tz,rho_ice,this->Sid()); - - /*Verbose some results in debug mode: {{{*/ - if(VerboseSmb() && 0){ - _printf_("smb log: count[" << count << "] m[" << m << "] " - << setprecision(16) << "T[" << cellsum(T,m) << "] " - << "d[" << cellsum(d,m) << "] " - << "dz[" << cellsum(dz,m) << "] " - << "a[" << cellsum(a,m) << "] " - << "W[" << cellsum(W,m) << "] " - << "re[" << cellsum(re,m) << "] " - << "gdn[" << cellsum(gdn,m) << "] " - << "gsp[" << cellsum(gsp,m) << "] " - << "swf[" << netSW << "] " - << "lwf[" << netLW << "] " - << "a[" << a << "] " - << "te[" << teValue << "] " - << "\n"); - } /*}}}*/ - - meanLHF = meanLHF + lhf*smb_dt/dt; - meanSHF = meanSHF + shf*smb_dt/dt; - - /*Sum component mass changes [kg m-2]*/ - sumMassAdd = mAdd + sumMassAdd; - sumM = M + sumM; - sumR = R + sumR; - sumW = cellsum(W,m); - sumP = P + sumP; - sumEC = sumEC + EC; // evap (-)/cond(+) - - /*Calculate total system mass:*/ - sumMass=0; - fac=0; - for(int i=0;iSid()==0 && IssmComm::GetRank()==0){ - if (T[m-1]!=T_bottom) _printf_("T(end)~=T_bottom" << "\n"); - } - - /*Free ressources: */ - xDelete(swf); - - /*increase counter:*/ - count++; - } //for (t=time;tSid()==0 && IssmComm::GetRank()==0)_printf0_("Time: t=" << setprecision(8) << timeinputs/365.0/24.0/3600.0 << " yr/" << (time+dt)/365.0/24.0/3600.0 << " yr" << setprecision(3) << " Step: " << count << "\n"); + + /*Get daily accumulated inputs {{{*/ + if (count>1){ + Input2 *sumEC_input = this->GetInput2(SmbECEnum); _assert_(sumEC_input); + Input2 *sumM_input = this->GetInput2(SmbMeltEnum); _assert_(sumM_input); + Input2 *sumR_input = this->GetInput2(SmbRunoffEnum); _assert_(sumR_input); + Input2 *sumP_input = this->GetInput2(SmbPrecipitationEnum); _assert_(sumP_input); + Input2 *ULW_input = this->GetInput2(SmbMeanULWEnum); _assert_(ULW_input); + Input2 *LW_input = this->GetInput2(SmbNetLWEnum); _assert_(LW_input); + Input2 *SW_input = this->GetInput2(SmbNetSWEnum); _assert_(SW_input); + Input2 *LHF_input = this->GetInput2(SmbMeanLHFEnum); _assert_(LHF_input); + Input2 *SHF_input = this->GetInput2(SmbMeanSHFEnum); _assert_(SHF_input); + Input2 *DzAdd_input = this->GetInput2(SmbDzAddEnum); _assert_(DzAdd_input); + Input2 *MassAdd_input = this->GetInput2(SmbMAddEnum); _assert_(MassAdd_input); + Input2 *InitMass_input = this->GetInput2(SmbMInitnum); _assert_(InitMass_input); + + ULW_input->GetInputAverage(&meanULW); + LW_input->GetInputAverage(&netLW); + SW_input->GetInputAverage(&netSW); + LHF_input->GetInputAverage(&meanLHF); + SHF_input->GetInputAverage(&meanSHF); + DzAdd_input->GetInputAverage(&dz_add); + MassAdd_input->GetInputAverage(&sumMassAdd); + sumMassAdd=sumMassAdd*dt; + InitMass_input->GetInputAverage(&initMass); + sumEC_input->GetInputAverage(&sumEC); + sumEC=sumEC*dt*rho_ice; + sumM_input->GetInputAverage(&sumM); + sumM=sumM*dt*rho_ice; + sumR_input->GetInputAverage(&sumR); + sumR=sumR*dt*rho_ice; + sumP_input->GetInputAverage(&sumP); + sumP=sumP*dt*rho_ice; + } + /*}}}*/ + + // Get time forcing inputs + Input2 *Ta_input = this->GetInput2(SmbTaEnum,timeinputs); _assert_(Ta_input); + Input2 *V_input = this->GetInput2(SmbVEnum,timeinputs); _assert_(V_input); + Input2 *Dlwr_input= this->GetInput2(SmbDlwrfEnum,timeinputs); _assert_(Dlwr_input); + Input2 *Dswr_input= this->GetInput2(SmbDswrfEnum,timeinputs); _assert_(Dswr_input); + Input2 *P_input = this->GetInput2(SmbPEnum,timeinputs); _assert_(P_input); + Input2 *eAir_input= this->GetInput2(SmbEAirEnum,timeinputs); _assert_(eAir_input); + Input2 *pAir_input= this->GetInput2(SmbPAirEnum,timeinputs); _assert_(pAir_input); + Input2 *teValue_input= this->GetInput2(SmbTeValueEnum,timeinputs); _assert_(teValue_input); + Input2 *aValue_input= this->GetInput2(SmbAValueEnum,timeinputs); _assert_(aValue_input); + + /*extract daily data:{{{*/ + Ta_input->GetInputValue(&Ta,gauss);//screen level air temperature [K] + V_input->GetInputValue(&V,gauss); //wind speed [m s-1] + Dlwr_input->GetInputValue(&dlw,gauss); //downward longwave radiation flux [W m-2] + Dswr_input->GetInputValue(&dsw,gauss); //downward shortwave radiation flux [W m-2] + P_input->GetInputValue(&P,gauss); //precipitation [kg m-2] + eAir_input->GetInputValue(&eAir,gauss); //screen level vapor pressure [Pa] + pAir_input->GetInputValue(&pAir,gauss); // screen level air pressure [Pa] + teValue_input->GetInputValue(&teValue,gauss); // Emissivity [0-1] + aValue_input->GetInputValue(&aValue,gauss); // Albedo [0 1] + //_printf_("Time: " << t << " Ta: " << Ta << " V: " << V << " dlw: " << dlw << " dsw: " << dsw << " P: " << P << " eAir: " << eAir << " pAir: " << pAir << "\n"); + /*}}}*/ + + /*Snow grain metamorphism:*/ + if(isgraingrowth)grainGrowth(&re, &gdn, &gsp, T, dz, d, W, smb_dt, m, aIdx,this->Sid()); + + /*Snow, firn and ice albedo:*/ + if(isalbedo)albedo(&a,aIdx,re,d,cldFrac,aIce,aSnow,aValue,adThresh,T,W,P,EC,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()); + + /*Calculate net shortwave [W m-2]*/ + netSW = netSW + cellsum(swf,m)*smb_dt/dt; + + /*Thermal profile computation:*/ + if(isthermal)thermo(&EC, &T, &ulw, dz, d, swf, dlw, Ta, V, eAir, pAir, teValue, W[0], smb_dt, m, Vz, Tz, thermo_scaling,rho_ice,this->Sid()); + + /*Change in thickness of top cell due to evaporation/condensation assuming same density as top cell. + * need to fix this in case all or more of cell evaporates */ + dz[0] = dz[0] + EC / d[0]; + + /*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()); + + /*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, &R, &mAdd, &dz_add, &T, &d, &dz, &W, &a, &re, &gdn, &gsp, &m, dzMin, zMax, zMin, zTop,rho_ice,this->Sid()); + + /*Allow non-melt densification and determine compaction [m]*/ + if(isdensification)densification(&d,&dz, T, re, denIdx, C, smb_dt, Tmean,rho_ice,m,this->Sid()); + + /*Calculate upward longwave radiation flux [W m-2] not used in energy balance. Calculated for every + * sub-time step in thermo equations*/ + //ulw = 5.67E-8 * pow(T[0],4.0) * teValue; // + deltatest here + + /*Calculate net longwave [W m-2]*/ + meanULW = meanULW + ulw*smb_dt/dt; + netLW = netLW + (dlw - ulw)*smb_dt/dt; + + /*Calculate turbulent heat fluxes [W m-2]*/ + if(isturbulentflux)turbulentFlux(&shf, &lhf, &dayEC, Ta, T[0], V, eAir, pAir, d[0], W[0], Vz, Tz,rho_ice,this->Sid()); + + /*Verbose some results in debug mode: {{{*/ + if(VerboseSmb() && 0){ + _printf_("smb log: count[" << count << "] m[" << m << "] " + << setprecision(16) << "T[" << cellsum(T,m) << "] " + << "d[" << cellsum(d,m) << "] " + << "dz[" << cellsum(dz,m) << "] " + << "a[" << cellsum(a,m) << "] " + << "W[" << cellsum(W,m) << "] " + << "re[" << cellsum(re,m) << "] " + << "gdn[" << cellsum(gdn,m) << "] " + << "gsp[" << cellsum(gsp,m) << "] " + << "swf[" << netSW << "] " + << "lwf[" << netLW << "] " + << "a[" << a << "] " + << "te[" << teValue << "] " + << "\n"); + } /*}}}*/ + + meanLHF = meanLHF + lhf*smb_dt/dt; + meanSHF = meanSHF + shf*smb_dt/dt; + + /*Sum component mass changes [kg m-2]*/ + sumMassAdd = mAdd + sumMassAdd; + sumM = M + sumM; + sumR = R + sumR; + sumW = cellsum(W,m); + sumP = P + sumP; + sumEC = sumEC + EC; // evap (-)/cond(+) + + /*Calculate total system mass:*/ + sumMass=0; + fac=0; + for(int i=0;iSid()==0 && IssmComm::GetRank()==0){ + if (T[m-1]!=T_bottom) _printf_("T(end)~=T_bottom" << "\n"); + } /*Save generated inputs: */ @@ -3995,7 +3996,10 @@ this->SetElementInput(SmbMeanLHFEnum,meanLHF); this->SetElementInput(SmbMeanSHFEnum,meanSHF); - this->SetElementInput(SmbDzAddEnum,sumdz_add); + this->SetElementInput(SmbDzAddEnum,dz_add); + this->SetElementInput(SmbMInitnum,initMass); this->SetElementInput(SmbMAddEnum,sumMassAdd/dt); + this->SetElementInput(SmbWAddEnum,sumW/dt); this->SetElementInput(SmbFACEnum,fac/1000.); // output in meters + this->SetElementInput(SmbECDtEnum,EC); /*Free allocations:{{{*/ @@ -4016,4 +4020,5 @@ if(aini) xDelete(aini); if(Tini) xDelete(Tini); + if(swf) xDelete(swf); delete gauss; Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 24682) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 24683) @@ -172,5 +172,5 @@ void SmbSemic(); int Sid(); - void SmbGemb(); + void SmbGemb(IssmDouble timeinputs, int count); void StrainRateESA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input2* vx_input,Input2* vy_input); void StrainRateFS(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input2* vx_input,Input2* vy_input,Input2* vz_input); Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp (revision 24682) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp (revision 24683) @@ -6,4 +6,6 @@ #include "../../shared/shared.h" #include "../../toolkits/toolkits.h" +#include "../modules.h" +#include "../../classes/Inputs2/TransientInput2.h" const double Pi = 3.141592653589793; @@ -32,7 +34,49 @@ void Gembx(FemModel* femmodel){ /*{{{*/ - for(int i=0;ielements->Size();i++){ - Element* element=xDynamicCast(femmodel->elements->GetObjectByOffset(i)); - element->SmbGemb(); + int count=0; + IssmDouble time,dt,finaltime,starttime; + IssmDouble timeclim=0.0; + IssmDouble t,smb_dt; + IssmDouble delta; + bool isclimatology=false; + + femmodel->parameters->FindParam(&time,TimeEnum); /*transient core time at which we run the smb core*/ + femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum); /*transient core time step*/ + femmodel->parameters->FindParam(&finaltime,TimesteppingFinalTimeEnum); + femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum); + femmodel->parameters->FindParam(&smb_dt,SmbDtEnum); /*time period for the smb solution, usually smaller than the glaciological dt*/ + femmodel->parameters->FindParam(&isclimatology,SmbIsclimatologyEnum); + + //before starting loop, realize that the transient core runs this smb_core at time = time +deltaT. + //go back to time - deltaT: + time-=dt; + + IssmDouble timeinputs = time; + + /*Start loop: */ + count=1; + for (t=time;telements->Size();i++){ + Element* element=xDynamicCast(femmodel->elements->GetObjectByOffset(i)); + + timeclim=time; + if (isclimatology){ + //If this is a climatology, we need to repeat the forcing after the final time + TransientInput2* Ta_input_tr = element->inputs2->GetTransientInput(SmbTaEnum); _assert_(Ta_input_tr); + + /*Get temperature climatology value*/ + int offsetend = Ta_input_tr->GetTimeInputOffset(finaltime); + IssmDouble time0 = Ta_input_tr->GetTimeByOffset(-1); + IssmDouble timeend = Ta_input_tr->GetTimeByOffset(offsetend); + if (time>time0 & timeend>time0){ + delta=(time-time0) - (timeend-time0)*(reCast((time-time0)/(timeend-time0))); + timeclim=time0+delta; + } + } + timeinputs = t-time+timeclim; + element->SmbGemb(timeinputs,count); + } + count=count+1; } @@ -1045,5 +1089,5 @@ // SWs and SWss coefficients need to be better constranted. Greuell // and Konzelmann 1994 used SWs = 0.36 and SWss = 0.64 as this the - // the // of SW radiation with wavelengths > and < 800 nm + // the % of SW radiation with wavelengths > and < 800 nm // respectively. This, however, may not account for the fact that // the albedo of wavelengths > 800 nm has a much lower albedo. @@ -2107,7 +2151,7 @@ // calculate the Bulk Richardson Number (Ri) - Ri = (2.0*9.81* (Vz - z0) * (Ta - Ts)) / ((Ta + Ts)* pow(V,2)); - - // calculate Monin-Obukhov stability factors 'coef_M' and 'coef_H' + Ri = (2.0*9.81* (Vz - z0) * (Ta - Ts)) / ((Ta + Ts)* pow(V,2.0)); + + // calculate Monin-Obukhov stability factors 'coefM' and 'coefH' // do not allow Ri to exceed 0.19 Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 24682) +++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 24683) @@ -149,4 +149,5 @@ syn keyword cConstant FrictionCouplingEnum syn keyword cConstant FrictionDeltaEnum +syn keyword cConstant FrictionEffectivePressureLimitEnum syn keyword cConstant FrictionFEnum syn keyword cConstant FrictionGammaEnum @@ -717,4 +718,5 @@ syn keyword cConstant SmbEAirEnum syn keyword cConstant SmbECEnum +syn keyword cConstant SmbECDtEnum syn keyword cConstant SmbECiniEnum syn keyword cConstant SmbElaEnum @@ -774,4 +776,5 @@ syn keyword cConstant SmbVzEnum syn keyword cConstant SmbWEnum +syn keyword cConstant SmbWAddEnum syn keyword cConstant SmbWiniEnum syn keyword cConstant SmbZMaxEnum @@ -1335,5 +1338,4 @@ syn keyword cType Cfsurfacesquare syn keyword cType Channel -syn keyword cType classes syn keyword cType Constraint syn keyword cType Constraints @@ -1342,6 +1344,6 @@ syn keyword cType ControlInput2 syn keyword cType Covertree +syn keyword cType DataSetParam syn keyword cType DatasetInput2 -syn keyword cType DataSetParam syn keyword cType Definition syn keyword cType DependentObject @@ -1355,6 +1357,6 @@ syn keyword cType ElementInput2 syn keyword cType ElementMatrix +syn keyword cType ElementVector syn keyword cType Elements -syn keyword cType ElementVector syn keyword cType ExponentialVariogram syn keyword cType ExternalResult @@ -1363,10 +1365,9 @@ syn keyword cType Friction syn keyword cType Gauss -syn keyword cType GaussianVariogram -syn keyword cType gaussobjects syn keyword cType GaussPenta syn keyword cType GaussSeg syn keyword cType GaussTetra syn keyword cType GaussTria +syn keyword cType GaussianVariogram syn keyword cType GenericExternalResult syn keyword cType GenericOption @@ -1383,5 +1384,4 @@ syn keyword cType IssmDirectApplicInterface syn keyword cType IssmParallelDirectApplicInterface -syn keyword cType krigingobjects syn keyword cType Load syn keyword cType Loads @@ -1394,5 +1394,4 @@ syn keyword cType Matice syn keyword cType Matlitho -syn keyword cType matrixobjects syn keyword cType MatrixParam syn keyword cType Misfit @@ -1407,6 +1406,6 @@ syn keyword cType Observations syn keyword cType Option +syn keyword cType OptionUtilities syn keyword cType Options -syn keyword cType OptionUtilities syn keyword cType Param syn keyword cType Parameters @@ -1422,10 +1421,10 @@ syn keyword cType Regionaloutput syn keyword cType Results +syn keyword cType RiftStruct syn keyword cType Riftfront -syn keyword cType RiftStruct syn keyword cType Seg syn keyword cType SegInput2 +syn keyword cType SegRef syn keyword cType Segment -syn keyword cType SegRef syn keyword cType SpcDynamic syn keyword cType SpcStatic @@ -1446,4 +1445,8 @@ syn keyword cType Vertex syn keyword cType Vertices +syn keyword cType classes +syn keyword cType gaussobjects +syn keyword cType krigingobjects +syn keyword cType matrixobjects syn keyword cType AdjointBalancethickness2Analysis syn keyword cType AdjointBalancethicknessAnalysis @@ -1464,6 +1467,6 @@ syn keyword cType FreeSurfaceBaseAnalysis syn keyword cType FreeSurfaceTopAnalysis +syn keyword cType GLheightadvectionAnalysis syn keyword cType GiaIvinsAnalysis -syn keyword cType GLheightadvectionAnalysis syn keyword cType HydrologyDCEfficientAnalysis syn keyword cType HydrologyDCInefficientAnalysis Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 24682) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 24683) @@ -715,4 +715,5 @@ SmbEAirEnum, SmbECEnum, + SmbECDtEnum, SmbECiniEnum, SmbElaEnum, @@ -734,4 +735,5 @@ SmbMeanULWEnum, SmbMeltEnum, + SmbMInitnum, SmbMonthlytemperaturesEnum, SmbNetLWEnum, @@ -772,4 +774,5 @@ SmbVzEnum, SmbWEnum, + SmbWAddEnum, SmbWiniEnum, SmbZMaxEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 24682) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 24683) @@ -720,4 +720,5 @@ case SmbEAirEnum : return "SmbEAir"; case SmbECEnum : return "SmbEC"; + case SmbECDtEnum : return "SmbECDt"; case SmbECiniEnum : return "SmbECini"; case SmbElaEnum : return "SmbEla"; @@ -777,4 +778,5 @@ case SmbVzEnum : return "SmbVz"; case SmbWEnum : return "SmbW"; + case SmbWAddEnum : return "SmbWAdd"; case SmbWiniEnum : return "SmbWini"; case SmbZMaxEnum : return "SmbZMax"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 24682) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 24683) @@ -735,4 +735,5 @@ else if (strcmp(name,"SmbEAir")==0) return SmbEAirEnum; else if (strcmp(name,"SmbEC")==0) return SmbECEnum; + else if (strcmp(name,"SmbECDt")==0) return SmbECDtEnum; else if (strcmp(name,"SmbECini")==0) return SmbECiniEnum; else if (strcmp(name,"SmbEla")==0) return SmbElaEnum; @@ -751,9 +752,9 @@ else if (strcmp(name,"SmbMassBalanceTransient")==0) return SmbMassBalanceTransientEnum; else if (strcmp(name,"SmbMeanLHF")==0) return SmbMeanLHFEnum; - else if (strcmp(name,"SmbMeanSHF")==0) return SmbMeanSHFEnum; else stage=7; } if(stage==7){ - if (strcmp(name,"SmbMeanULW")==0) return SmbMeanULWEnum; + if (strcmp(name,"SmbMeanSHF")==0) return SmbMeanSHFEnum; + else if (strcmp(name,"SmbMeanULW")==0) return SmbMeanULWEnum; else if (strcmp(name,"SmbMelt")==0) return SmbMeltEnum; else if (strcmp(name,"SmbMonthlytemperatures")==0) return SmbMonthlytemperaturesEnum; @@ -795,4 +796,5 @@ else if (strcmp(name,"SmbVz")==0) return SmbVzEnum; else if (strcmp(name,"SmbW")==0) return SmbWEnum; + else if (strcmp(name,"SmbWAdd")==0) return SmbWAddEnum; else if (strcmp(name,"SmbWini")==0) return SmbWiniEnum; else if (strcmp(name,"SmbZMax")==0) return SmbZMaxEnum; @@ -873,10 +875,10 @@ else if (strcmp(name,"Outputdefinition17")==0) return Outputdefinition17Enum; else if (strcmp(name,"Outputdefinition18")==0) return Outputdefinition18Enum; - else if (strcmp(name,"Outputdefinition19")==0) return Outputdefinition19Enum; - else if (strcmp(name,"Outputdefinition20")==0) return Outputdefinition20Enum; else stage=8; } if(stage==8){ - if (strcmp(name,"Outputdefinition21")==0) return Outputdefinition21Enum; + if (strcmp(name,"Outputdefinition19")==0) return Outputdefinition19Enum; + else if (strcmp(name,"Outputdefinition20")==0) return Outputdefinition20Enum; + else if (strcmp(name,"Outputdefinition21")==0) return Outputdefinition21Enum; else if (strcmp(name,"Outputdefinition22")==0) return Outputdefinition22Enum; else if (strcmp(name,"Outputdefinition23")==0) return Outputdefinition23Enum; @@ -996,10 +998,10 @@ else if (strcmp(name,"BoolInput2")==0) return BoolInput2Enum; else if (strcmp(name,"IntInput2")==0) return IntInput2Enum; - else if (strcmp(name,"BoolParam")==0) return BoolParamEnum; - else if (strcmp(name,"Boundary")==0) return BoundaryEnum; else stage=9; } if(stage==9){ - if (strcmp(name,"BuddJacka")==0) return BuddJackaEnum; + if (strcmp(name,"BoolParam")==0) return BoolParamEnum; + else if (strcmp(name,"Boundary")==0) return BoundaryEnum; + else if (strcmp(name,"BuddJacka")==0) return BuddJackaEnum; else if (strcmp(name,"CalvingDev2")==0) return CalvingDev2Enum; else if (strcmp(name,"CalvingHab")==0) return CalvingHabEnum; @@ -1119,10 +1121,10 @@ else if (strcmp(name,"IcefrontMassFlux")==0) return IcefrontMassFluxEnum; else if (strcmp(name,"IcefrontMassFluxLevelset")==0) return IcefrontMassFluxLevelsetEnum; - else if (strcmp(name,"Incremental")==0) return IncrementalEnum; - else if (strcmp(name,"Indexed")==0) return IndexedEnum; else stage=10; } if(stage==10){ - if (strcmp(name,"IntExternalResult")==0) return IntExternalResultEnum; + if (strcmp(name,"Incremental")==0) return IncrementalEnum; + else if (strcmp(name,"Indexed")==0) return IndexedEnum; + else if (strcmp(name,"IntExternalResult")==0) return IntExternalResultEnum; else if (strcmp(name,"IntInput")==0) return IntInputEnum; else if (strcmp(name,"ElementInput2")==0) return ElementInput2Enum; @@ -1242,10 +1244,10 @@ else if (strcmp(name,"Regionaloutput")==0) return RegionaloutputEnum; else if (strcmp(name,"Regular")==0) return RegularEnum; - else if (strcmp(name,"RecoveryAnalysis")==0) return RecoveryAnalysisEnum; - else if (strcmp(name,"Riftfront")==0) return RiftfrontEnum; else stage=11; } if(stage==11){ - if (strcmp(name,"SIAApproximation")==0) return SIAApproximationEnum; + if (strcmp(name,"RecoveryAnalysis")==0) return RecoveryAnalysisEnum; + else if (strcmp(name,"Riftfront")==0) return RiftfrontEnum; + else if (strcmp(name,"SIAApproximation")==0) return SIAApproximationEnum; else if (strcmp(name,"SMBcomponents")==0) return SMBcomponentsEnum; else if (strcmp(name,"SMBd18opdd")==0) return SMBd18opddEnum;