Index: ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 24682) +++ ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 24683) @@ -714,6 +714,7 @@ SmbDziniEnum, SmbEAirEnum, SmbECEnum, + SmbECDtEnum, SmbECiniEnum, SmbElaEnum, SmbEvaporationEnum, @@ -733,6 +734,7 @@ SmbMeanSHFEnum, SmbMeanULWEnum, SmbMeltEnum, + SmbMInitnum, SmbMonthlytemperaturesEnum, SmbNetLWEnum, SmbNetSWEnum, @@ -771,6 +773,7 @@ SmbVmeanEnum, SmbVzEnum, SmbWEnum, + SmbWAddEnum, SmbWiniEnum, SmbZMaxEnum, SmbZMinEnum, Index: ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 24682) +++ ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 24683) @@ -719,6 +719,7 @@ case SmbDziniEnum : return "SmbDzini"; case SmbEAirEnum : return "SmbEAir"; case SmbECEnum : return "SmbEC"; + case SmbECDtEnum : return "SmbECDt"; case SmbECiniEnum : return "SmbECini"; case SmbElaEnum : return "SmbEla"; case SmbEvaporationEnum : return "SmbEvaporation"; @@ -776,6 +777,7 @@ case SmbVmeanEnum : return "SmbVmean"; case SmbVzEnum : return "SmbVz"; case SmbWEnum : return "SmbW"; + case SmbWAddEnum : return "SmbWAdd"; case SmbWiniEnum : return "SmbWini"; case SmbZMaxEnum : return "SmbZMax"; case SmbZMinEnum : return "SmbZMin"; Index: ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 24682) +++ ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 24683) @@ -734,6 +734,7 @@ else if (strcmp(name,"SmbDzini")==0) return SmbDziniEnum; 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; else if (strcmp(name,"SmbEvaporation")==0) return SmbEvaporationEnum; @@ -750,11 +751,11 @@ else if (strcmp(name,"SmbMassBalanceSubstep")==0) return SmbMassBalanceSubstepEnum; 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; else if (strcmp(name,"SmbNetLW")==0) return SmbNetLWEnum; @@ -794,6 +795,7 @@ else if (strcmp(name,"SmbVmean")==0) return SmbVmeanEnum; 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; else if (strcmp(name,"SmbZMin")==0) return SmbZMinEnum; @@ -872,12 +874,12 @@ else if (strcmp(name,"Outputdefinition16")==0) return Outputdefinition16Enum; 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; else if (strcmp(name,"Outputdefinition24")==0) return Outputdefinition24Enum; @@ -995,12 +997,12 @@ else if (strcmp(name,"BoolInput")==0) return BoolInputEnum; 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; else if (strcmp(name,"CalvingLevermann")==0) return CalvingLevermannEnum; @@ -1118,12 +1120,12 @@ else if (strcmp(name,"IceVolumeScaled")==0) return IceVolumeScaledEnum; 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; else if (strcmp(name,"SegInput2")==0) return SegInput2Enum; @@ -1241,12 +1243,12 @@ else if (strcmp(name,"ProfilingSolutionTime")==0) return ProfilingSolutionTimeEnum; 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; else if (strcmp(name,"SMBforcing")==0) return SMBforcingEnum; Index: ../trunk-jpl/src/c/shared/Enum/Enum.vim =================================================================== --- ../trunk-jpl/src/c/shared/Enum/Enum.vim (revision 24682) +++ ../trunk-jpl/src/c/shared/Enum/Enum.vim (revision 24683) @@ -148,6 +148,7 @@ syn keyword cConstant FlowequationIsSSAEnum syn keyword cConstant FrictionCouplingEnum syn keyword cConstant FrictionDeltaEnum +syn keyword cConstant FrictionEffectivePressureLimitEnum syn keyword cConstant FrictionFEnum syn keyword cConstant FrictionGammaEnum syn keyword cConstant FrictionLawEnum @@ -716,6 +717,7 @@ syn keyword cConstant SmbDziniEnum syn keyword cConstant SmbEAirEnum syn keyword cConstant SmbECEnum +syn keyword cConstant SmbECDtEnum syn keyword cConstant SmbECiniEnum syn keyword cConstant SmbElaEnum syn keyword cConstant SmbEvaporationEnum @@ -773,6 +775,7 @@ syn keyword cConstant SmbVmeanEnum syn keyword cConstant SmbVzEnum syn keyword cConstant SmbWEnum +syn keyword cConstant SmbWAddEnum syn keyword cConstant SmbWiniEnum syn keyword cConstant SmbZMaxEnum syn keyword cConstant SmbZMinEnum @@ -1334,7 +1337,6 @@ syn keyword cType Cfsurfacelogvel syn keyword cType Cfsurfacesquare syn keyword cType Channel -syn keyword cType classes syn keyword cType Constraint syn keyword cType Constraints syn keyword cType Contour @@ -1341,8 +1343,8 @@ syn keyword cType Contours 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 syn keyword cType DoubleMatArrayParam @@ -1354,8 +1356,8 @@ syn keyword cType ElementHook 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 syn keyword cType FemModel @@ -1362,12 +1364,11 @@ syn keyword cType FileParam 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 syn keyword cType GenericParam @@ -1382,7 +1383,6 @@ syn keyword cType IoModel syn keyword cType IssmDirectApplicInterface syn keyword cType IssmParallelDirectApplicInterface -syn keyword cType krigingobjects syn keyword cType Load syn keyword cType Loads syn keyword cType Masscon @@ -1393,7 +1393,6 @@ syn keyword cType Matestar syn keyword cType Matice syn keyword cType Matlitho -syn keyword cType matrixobjects syn keyword cType MatrixParam syn keyword cType Misfit syn keyword cType Moulin @@ -1406,8 +1405,8 @@ syn keyword cType Observation 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 syn keyword cType Pengrid @@ -1421,12 +1420,12 @@ syn keyword cType Radar 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 syn keyword cType SpcTransient @@ -1445,6 +1444,10 @@ syn keyword cType VectorParam 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 syn keyword cType AdjointHorizAnalysis @@ -1463,8 +1466,8 @@ syn keyword cType ExtrudeFromTopAnalysis 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 syn keyword cType HydrologyGlaDSAnalysis Index: ../trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp =================================================================== --- ../trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp (revision 24682) +++ ../trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp (revision 24683) @@ -5,6 +5,8 @@ #include "./SurfaceMassBalancex.h" #include "../../shared/shared.h" #include "../../toolkits/toolkits.h" +#include "../modules.h" +#include "../../classes/Inputs2/TransientInput2.h" const double Pi = 3.141592653589793; const double CtoK = 273.15; // Kelvin to Celcius conversion/ice melt. point T in K @@ -31,9 +33,51 @@ 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; } } /*}}}*/ @@ -1044,7 +1088,7 @@ // 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. @@ -2106,9 +2150,9 @@ if(V< 0.01-Dtol)V=0.01; // calculate the Bulk Richardson Number (Ri) - Ri = (2.0*9.81* (Vz - z0) * (Ta - Ts)) / ((Ta + Ts)* pow(V,2)); + Ri = (2.0*9.81* (Vz - z0) * (Ta - Ts)) / ((Ta + Ts)* pow(V,2.0)); - // calculate Monin-Obukhov stability factors 'coef_M' and 'coef_H' + // calculate Monin-Obukhov stability factors 'coefM' and 'coefH' // do not allow Ri to exceed 0.19 Ri = min(Ri, 0.19); Index: ../trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Element.h (revision 24682) +++ ../trunk-jpl/src/c/classes/Elements/Element.h (revision 24683) @@ -171,7 +171,7 @@ void SetIntInput(Inputs2* inputs2,int enum_in,int value); 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); void StrainRateHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input2* vx_input,Input2* vy_input); Index: ../trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 24682) +++ ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 24683) @@ -561,7 +561,7 @@ this->parameters->FindParam(&isTemperatureScaled,SmbIstemperaturescaledEnum); this->parameters->FindParam(&isPrecipScaled,SmbIsprecipscaledEnum); - /*Recover present day temperature and precipitation*/ + /*Recover present day temperature and precipitation*/ DatasetInput2 *dinput3 = NULL; DatasetInput2 *dinput4 = NULL; int offset_t,offset_p,N; @@ -3587,7 +3587,7 @@ } /*}}}*/ -void Element::SmbGemb(){/*{{{*/ +void Element::SmbGemb(IssmDouble timeinputs, int count){/*{{{*/ /*only compute SMB at the surface: */ if (!IsOnSurface()) return; @@ -3604,9 +3604,6 @@ IssmDouble Vmean=0.0; 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; IssmDouble V=0.0; @@ -3617,6 +3614,7 @@ IssmDouble pAir=0.0; IssmDouble teValue=1.0; IssmDouble aValue=0.0; + IssmDouble dt,time,smb_dt; int aIdx=0; int denIdx=0; int dsnowIdx=0; @@ -3626,13 +3624,12 @@ IssmDouble shf=0.0; 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; IssmDouble dMass=0.0; @@ -3641,8 +3638,6 @@ IssmDouble init_scaling=0.0; IssmDouble thermo_scaling=1.0; IssmDouble adThresh=1023.0; - int offsetend=-1; - IssmDouble time0, timeend, delta; /*}}}*/ /*Output variables:{{{ */ IssmDouble* dz=NULL; @@ -3675,7 +3670,6 @@ IssmDouble* aini = NULL; IssmDouble* Tini = NULL; int m=0; - int count=0; /*}}}*/ /*Retrieve material properties and parameters:{{{ */ @@ -3686,8 +3680,6 @@ parameters->FindParam(&time,TimeEnum); /*transient core time at which we run the smb core*/ 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); parameters->FindParam(&denIdx,SmbDenIdxEnum); @@ -3697,7 +3689,6 @@ parameters->FindParam(&t0wet,SmbT0wetEnum); parameters->FindParam(&t0dry,SmbT0dryEnum); parameters->FindParam(&K,SmbKEnum); - parameters->FindParam(&isclimatology,SmbIsclimatologyEnum); parameters->FindParam(&isgraingrowth,SmbIsgraingrowthEnum); parameters->FindParam(&isalbedo,SmbIsalbedoEnum); parameters->FindParam(&isshortwave,SmbIsshortwaveEnum); @@ -3722,8 +3713,6 @@ Input2 *C_input = this->GetInput2(SmbCEnum); _assert_(C_input); 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; /*Retrieve input values:*/ @@ -3741,8 +3730,6 @@ C_input->GetInputValue(&C,gauss); Tz_input->GetInputValue(&Tz,gauss); Vz_input->GetInputValue(&Vz,gauss); - teValue_input->GetInputValue(&teValue,gauss); - aValue_input->GetInputValue(&aValue,gauss); /*}}}*/ /*First, check that the initial structures have been setup in GEMB. If not, initialize profile variables: layer thickness dz, * density d, temperature T, etc. {{{*/ @@ -3762,7 +3749,7 @@ EC_input = this->GetInput2(SmbECiniEnum); _assert_(EC_input); 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); if(m==2){ //Snow properties are initialized with default values. Vertical grid has to be initialized too @@ -3815,9 +3802,8 @@ this->inputs2->GetArray(SmbWEnum,this->lid,&W,&m); 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 _assert_(m>0); @@ -3829,152 +3815,167 @@ // 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; + //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; - } - } + if(VerboseSmb() && this->Sid()==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"); - /*Start loop: */ - count=1; - for (t=time;t1){ + 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); - if(VerboseSmb() && this->Sid()==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"); + 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; + } + /*}}}*/ - 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); + // 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"); - /*}}}*/ + /*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 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()); + /*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()); + /*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; + /*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()); + /*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]; + /*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()); + /*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()); + /*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()); + /*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 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 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()); + /*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"); - } /*}}}*/ + /*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; + 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(+) + /*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"); - } + /*Check bottom grid cell T is unchanged:*/ + if(VerboseSmb() && this->Sid()==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;tinputs2->SetArrayInput(SmbDzEnum,this->lid,dz,m); this->inputs2->SetArrayInput(SmbDEnum,this->lid,d,m); @@ -3994,9 +3995,12 @@ this->SetElementInput(SmbNetSWEnum,netSW); 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:{{{*/ if(dz) xDelete(dz); @@ -4015,6 +4019,7 @@ if(Wini) xDelete(Wini); if(aini) xDelete(aini); if(Tini) xDelete(Tini); + if(swf) xDelete(swf); delete gauss; /*}}}*/