Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 19553) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 19554) @@ -58,4 +58,5 @@ ./classes/Inputs/IntInput.cpp\ ./classes/Inputs/DoubleInput.cpp\ + ./classes/Inputs/DoubleArrayInput.cpp\ ./classes/Inputs/DatasetInput.cpp\ ./classes/Materials/Materials.cpp\ @@ -171,4 +172,5 @@ ./modules/SpcNodesx/SpcNodesx.cpp\ ./modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp\ + ./modules/SurfaceMassBalancex/Gembx.cpp\ ./modules/Reducevectorgtofx/Reducevectorgtofx.cpp\ ./modules/Reduceloadx/Reduceloadx.cpp\ Index: /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp (revision 19553) +++ /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp (revision 19554) @@ -23,7 +23,4 @@ bool isdelta18o,ismungsm,isd18opd; - /*Figure out smb model: */ - iomodel->Constant(&smb_model,SmbEnum); - /*Update elements: */ int counter=0; @@ -36,7 +33,32 @@ } + /*Figure out smb model: */ + iomodel->Constant(&smb_model,SmbEnum); + + + Input* bof=NULL; + Element* element=NULL; switch(smb_model){ case SMBforcingEnum: iomodel->FetchDataToInput(elements,SmbMassBalanceEnum,0.); + break; + case SMBgembEnum: + iomodel->FetchDataToInput(elements,SmbTaEnum); + iomodel->FetchDataToInput(elements,SmbVEnum); + iomodel->FetchDataToInput(elements,SmbDswrfEnum); + iomodel->FetchDataToInput(elements,SmbDlwrfEnum); + iomodel->FetchDataToInput(elements,SmbPEnum); + iomodel->FetchDataToInput(elements,SmbEAirEnum); + iomodel->FetchDataToInput(elements,SmbPAirEnum); + iomodel->FetchDataToInput(elements,SmbZTopEnum); + iomodel->FetchDataToInput(elements,SmbDzTopEnum); + iomodel->FetchDataToInput(elements,SmbDzMinEnum); + iomodel->FetchDataToInput(elements,SmbZYEnum); + iomodel->FetchDataToInput(elements,SmbZMaxEnum); + iomodel->FetchDataToInput(elements,SmbZMinEnum); + iomodel->FetchDataToInput(elements,SmbTmeanEnum); + iomodel->FetchDataToInput(elements,SmbCEnum); + iomodel->FetchDataToInput(elements,SmbTzEnum); + iomodel->FetchDataToInput(elements,SmbVzEnum); break; case SMBpddEnum: @@ -92,4 +114,6 @@ _error_("Surface mass balance model "<AddObject(iomodel->CopyConstantObject(SmbAIdxEnum)); + parameters->AddObject(iomodel->CopyConstantObject(SmbSwIdxEnum)); + parameters->AddObject(iomodel->CopyConstantObject(SmbDenIdxEnum)); + parameters->AddObject(iomodel->CopyConstantObject(SmbOutputFreqEnum)); + parameters->AddObject(iomodel->CopyConstantObject(SmbCldFracEnum)); + parameters->AddObject(iomodel->CopyConstantObject(SmbT0wetEnum)); + parameters->AddObject(iomodel->CopyConstantObject(SmbT0dryEnum)); + parameters->AddObject(iomodel->CopyConstantObject(SmbKEnum)); + parameters->AddObject(iomodel->CopyConstantObject(SmbASnowEnum)); + parameters->AddObject(iomodel->CopyConstantObject(SmbAIceEnum)); + parameters->AddObject(iomodel->CopyConstantObject(SmbDtEnum)); + parameters->AddObject(iomodel->CopyConstantObject(SmbIsgraingrowthEnum)); + parameters->AddObject(iomodel->CopyConstantObject(SmbIsalbedoEnum)); + parameters->AddObject(iomodel->CopyConstantObject(SmbIsshortwaveEnum)); + parameters->AddObject(iomodel->CopyConstantObject(SmbIsthermalEnum)); + parameters->AddObject(iomodel->CopyConstantObject(SmbIsaccumulationEnum)); + parameters->AddObject(iomodel->CopyConstantObject(SmbIsmeltEnum)); + parameters->AddObject(iomodel->CopyConstantObject(SmbIsdensificationEnum)); + parameters->AddObject(iomodel->CopyConstantObject(SmbIsturbulentfluxEnum)); break; case SMBpddEnum: @@ -196,4 +241,7 @@ /*Nothing to be done*/ break; + case SMBgembEnum: + Gembx(femmodel); + break; case SMBpddEnum: bool isdelta18o,ismungsm; Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 19553) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 19554) @@ -14,4 +14,5 @@ #include "../classes.h" #include "../../shared/shared.h" +#include "../../modules/SurfaceMassBalancex/SurfaceMassBalancex.h" /*}}}*/ @@ -1366,4 +1367,7 @@ } else if(vector_type==2){ //element vector + + IssmDouble value; + /*Are we in transient or static? */ if(M==iomodel->numberofelements){ @@ -1379,12 +1383,27 @@ else _error_("could not recognize nature of vector from code " << code); } - else { - _error_("transient element inputs not supported yet!"); - } - } - else{ - _error_("Cannot add input for vector type " << vector_type << " (not supported)"); - } -}/*}}}*/ + else if(M==iomodel->numberofelements+1){ + /*create transient input: */ + IssmDouble* times = xNew(N); + for(t=0;tSid()+t]; + switch(this->ObjectEnum()){ + case TriaEnum: transientinput->AddTimeInput(new TriaInput( vector_enum,&value,P0Enum)); break; + case PentaEnum: transientinput->AddTimeInput(new PentaInput(vector_enum,&value,P0Enum)); break; + case TetraEnum: transientinput->AddTimeInput(new TetraInput(vector_enum,&value,P0Enum)); break; + default: _error_("Not implemented yet"); + } + } + this->inputs->AddInput(transientinput); + xDelete(times); + } + else _error_("element vector is either numberofelements or numberofelements+1 long. Field provided (" << EnumToStringx(vector_enum) << ") is " << M << " long"); + } + else _error_("Cannot add input for vector type " << vector_type << " (not supported)"); +} +/*}}}*/ void Element::InputDuplicate(int original_enum,int new_enum){/*{{{*/ @@ -1866,5 +1885,5 @@ }/*}}}*/ -void Element::ResultInterpolation(int* pinterpolation,int* pnodesperelement,int output_enum){/*{{{*/ +void Element::ResultInterpolation(int* pinterpolation,int* pnodesperelement,int* parray_size, int output_enum){/*{{{*/ Input* input=this->inputs->GetInput(output_enum); @@ -1960,4 +1979,5 @@ *pinterpolation = input->GetResultInterpolation(); *pnodesperelement = input->GetResultNumberOfNodes(); + *parray_size = input->GetResultArraySize(); }/*}}}*/ void Element::ResultToPatch(IssmDouble* values,int nodesperelement,int output_enum){/*{{{*/ @@ -1967,4 +1987,12 @@ input->ResultToPatch(values,nodesperelement,this->Sid()); + +} /*}}}*/ +void Element::ResultToMatrix(IssmDouble* values,int ncols,int output_enum){/*{{{*/ + + Input* input=this->inputs->GetInput(output_enum); + if(!input) _error_("input "<ResultToMatrix(values,ncols,this->Sid()); } /*}}}*/ @@ -2093,4 +2121,261 @@ } /*}}}*/ +void Element::SmbGemb(){/*{{{*/ + + /*Intermediary variables: {{{*/ + bool isinitialized=false; + IssmDouble zTop,dzTop,zMax,zMin,zY,dzMin; + IssmDouble Tmean; + IssmDouble C; + IssmDouble Tz,Vz; + IssmDouble rho_ice, aSnow,aIce; + IssmDouble time,dt; + IssmDouble t,smb_dt; + IssmDouble Ta,V,dlw,dsw,P,eAir,pAir; + int aIdx=0; + int denIdx=0; + int swIdx=0; + IssmDouble cldFrac,t0wet, t0dry, K; + IssmDouble comp1,comp2; + IssmDouble ulw; + IssmDouble netSW; + IssmDouble netLW; + IssmDouble lhf, shf, dayEC; + IssmDouble initMass; + IssmDouble sumR, sumM, sumEC, sumP, sumW,sumMassAdd; + IssmDouble sumMass,dMass; + bool isgraingrowth,isalbedo,isshortwave,isthermal,isaccumulation,ismelt,isdensification,isturbulentflux; + /*}}}*/ + /*Output variables:{{{ */ + IssmDouble* dz=NULL; + IssmDouble* d = NULL; + IssmDouble* re = NULL; + IssmDouble* gdn = NULL; + IssmDouble* gsp = NULL; + IssmDouble EC = 0; + IssmDouble* W = NULL; + IssmDouble* a = NULL; + IssmDouble* swf=NULL; + IssmDouble* T = NULL; + IssmDouble T_bottom; + IssmDouble M; + IssmDouble R; + IssmDouble mAdd; + int m; + /*}}}*/ + + /*only compute SMB at the surface: */ + if (!IsOnSurface()) return; + + /*Retrieve material properties and parameters:{{{ */ + rho_ice = matpar->GetMaterialParameter(MaterialsRhoIceEnum); + parameters->FindParam(&aSnow,SmbASnowEnum); + parameters->FindParam(&aIce,SmbAIceEnum); + parameters->FindParam(&time,TimeEnum); /*transient core time at which we run the smb core*/ + parameters->FindParam(&dt,TimesteppingTimeStepEnum); /*transient core time step*/ + 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); + parameters->FindParam(&swIdx,SmbSwIdxEnum); + parameters->FindParam(&cldFrac,SmbCldFracEnum); + parameters->FindParam(&t0wet,SmbT0wetEnum); + parameters->FindParam(&t0dry,SmbT0dryEnum); + parameters->FindParam(&K,SmbKEnum); + parameters->FindParam(&isgraingrowth,SmbIsgraingrowthEnum); + parameters->FindParam(&isalbedo,SmbIsalbedoEnum); + parameters->FindParam(&isshortwave,SmbIsshortwaveEnum); + parameters->FindParam(&isthermal,SmbIsthermalEnum); + parameters->FindParam(&isaccumulation,SmbIsaccumulationEnum); + parameters->FindParam(&ismelt,SmbIsmeltEnum); + parameters->FindParam(&isdensification,SmbIsdensificationEnum); + parameters->FindParam(&isturbulentflux,SmbIsturbulentfluxEnum); + /*}}}*/ + /*Retrieve inputs: {{{*/ + Input* zTop_input=this->GetInput(SmbZTopEnum); _assert_(zTop_input); + Input* dzTop_input=this->GetInput(SmbDzTopEnum); _assert_(dzTop_input); + Input* dzMin_input=this->GetInput(SmbDzMinEnum); _assert_(dzMin_input); + Input* zMax_input=this->GetInput(SmbZMaxEnum); _assert_(zMax_input); + Input* zMin_input=this->GetInput(SmbZMinEnum); _assert_(zMin_input); + Input* zY_input=this->GetInput(SmbZYEnum); _assert_(zY_input); + Input* Tmean_input=this->GetInput(SmbTmeanEnum); _assert_(Tmean_input); + Input* C_input=this->GetInput(SmbCEnum); _assert_(C_input); + Input* Tz_input=this->GetInput(SmbTzEnum); _assert_(Tz_input); + Input* Vz_input=this->GetInput(SmbVzEnum); _assert_(Vz_input); + Input* Ta_input=this->GetInput(SmbTaEnum); _assert_(Ta_input); + Input* V_input=this->GetInput(SmbVEnum); _assert_(V_input); + Input* Dlwr_input=this->GetInput(SmbDlwrfEnum); _assert_(Dlwr_input); + Input* Dswr_input=this->GetInput(SmbDswrfEnum); _assert_(Dswr_input); + Input* P_input=this->GetInput(SmbPEnum); _assert_(P_input); + Input* eAir_input=this->GetInput(SmbEAirEnum); _assert_(eAir_input); + Input* pAir_input=this->GetInput(SmbPAirEnum); _assert_(pAir_input); + Input* isinitialized_input=this->inputs->GetInput(SmbIsInitializedEnum); + + /*Retrieve input values:*/ + Gauss* gauss=this->NewGauss(1); gauss->GaussPoint(0); + + zTop_input->GetInputValue(&zTop,gauss); + dzTop_input->GetInputValue(&dzTop,gauss); + dzMin_input->GetInputValue(&dzMin,gauss); + zMax_input->GetInputValue(&zMax,gauss); + zMin_input->GetInputValue(&zMin,gauss); + zY_input->GetInputValue(&zY,gauss); + Tmean_input->GetInputValue(&Tmean,gauss); + C_input->GetInputValue(&C,gauss); + Tz_input->GetInputValue(&Tz,gauss); + Vz_input->GetInputValue(&Vz,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. {{{*/ + if(!isinitialized_input){ + + /*Initialize grid:*/ + GembgridInitialize(&dz, &m, zTop, dzTop, zMax, zY); + //if(this->Sid()==1) 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;iAddInput(new BoolInput(SmbIsInitializedEnum,true)); + } + else{ + /*Recover inputs: */ + DoubleArrayInput* dz_input= dynamic_cast(this->GetInput(SmbDzEnum)); + DoubleArrayInput* d_input= dynamic_cast(this->GetInput(SmbDEnum)); + DoubleArrayInput* re_input= dynamic_cast(this->GetInput(SmbReEnum)); + DoubleArrayInput* gdn_input= dynamic_cast(this->GetInput(SmbGdnEnum)); + DoubleArrayInput* gsp_input= dynamic_cast(this->GetInput(SmbGspEnum)); + DoubleInput* EC_input= dynamic_cast(this->GetInput(SmbECEnum)); + DoubleArrayInput* W_input= dynamic_cast(this->GetInput(SmbWEnum)); + DoubleArrayInput* a_input= dynamic_cast(this->GetInput(SmbAEnum)); + DoubleArrayInput* T_input= dynamic_cast(this->GetInput(SmbTEnum)); + + /*Recover arrays: */ + dz_input->GetValues(&dz,&m); + d_input->GetValues(&d,&m); + re_input->GetValues(&re,&m); + gdn_input->GetValues(&gdn,&m); + gsp_input->GetValues(&gsp,&m); + EC_input->GetInputValue(&EC); + W_input->GetValues(&W,&m); + a_input->GetValues(&a,&m); + T_input->GetValues(&T,&m); + + } /*}}}*/ + + // determine initial mass [kg] + initMass=0; for(int i=0;iGetInputValue(&Ta,gauss,t);//screen level air temperature [K] + V_input->GetInputValue(&V,gauss,t); //wind speed [m s-1] + Dlwr_input->GetInputValue(&dlw,gauss,t); //downward longwave radiation flux [W m-2] + Dswr_input->GetInputValue(&dsw,gauss,t); //downward shortwave radiation flux [W m-2] + P_input->GetInputValue(&P,gauss,t); //precipitation [kg m-2] + eAir_input->GetInputValue(&eAir,gauss,t); //screen level vapor pressure [Pa] + pAir_input->GetInputValue(&pAir,gauss,t); // screen level air pressure [Pa] + //_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); + + /*Snow, firn and ice albedo:*/ + if(isalbedo)albedo(a,aIdx,re,d,cldFrac,aIce, aSnow,T,W,P,EC,t0wet,t0dry,K,smb_dt,m); + + /*Distribution of absorbed short wave radation with depth:*/ + if(isshortwave)shortwave(&swf, swIdx, aIdx, dsw, a[0], d, dz, re,m); + + /*Thermal profile computation:*/ + if(isthermal)thermo(&EC, T, dz, d, swf, dlw, Ta, V, eAir, pAir, W[0], smb_dt, Vz, Tz,m); + + /*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, Ta, P, dzMin, aSnow); + + /*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]*/ + comp2 = cellsum(dz,m); + if(ismelt)melt(&M, &R, &mAdd, &T, &d, &dz, &W, &a, &re, &gdn, &gsp, &m, dzMin, zMax, zMin); + comp2 = (comp2 - cellsum(dz,m)); + + /*Allow non-melt densification and determine compaction [m]*/ + comp1 = cellsum(dz,m); + if(isdensification)densification(d,dz, T, re, denIdx, C, smb_dt, Tmean,m); + comp1 = (comp1 - cellsum(dz,m)); + + /*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); + + /*Calculate net shortwave and longwave [W m-2]*/ + netSW = cellsum(swf,m); + netLW = dlw - ulw; + + /*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); + + /*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; for(int i=0;iAddInput(new DoubleArrayInput(SmbDzEnum,dz,m)); + this->AddInput(new DoubleArrayInput(SmbDEnum,d,m)); + this->AddInput(new DoubleArrayInput(SmbReEnum,re,m)); + this->AddInput(new DoubleArrayInput(SmbGdnEnum,gdn,m)); + this->AddInput(new DoubleArrayInput(SmbGspEnum,gsp,m)); + this->AddInput(new DoubleInput(SmbECEnum,EC)); + this->AddInput(new DoubleArrayInput(SmbWEnum,W,m)); + this->AddInput(new DoubleArrayInput(SmbAEnum,a,m)); + this->AddInput(new DoubleArrayInput(SmbTEnum,T,m)); + + /*Free allocations:{{{*/ + xDelete(dz); + xDelete(d); + xDelete(re); + xDelete(gdn); + xDelete(gsp); + xDelete(W); + xDelete(a); + xDelete(T); + /*}}}*/ +} +/*}}}*/ void Element::StrainRateFS(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/ /*Compute the 3d Strain Rate (6 components): Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 19553) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 19554) @@ -130,9 +130,11 @@ IssmDouble PureIceEnthalpy(IssmDouble pressure); void ResetIceLevelset(void); - void ResultInterpolation(int* pinterpolation,int*nodesperelement,int output_enum); + void ResultInterpolation(int* pinterpolation,int*nodesperelement,int* parray_size, int output_enum); void ResultToPatch(IssmDouble* values,int nodesperelement,int output_enum); + void ResultToMatrix(IssmDouble* values,int ncols,int output_enum); void ResultToVector(Vector* vector,int output_enum); void SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum); int Sid(); + void SmbGemb(); void StrainRateFS(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input); void StrainRateHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 19553) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 19554) @@ -1040,19 +1040,24 @@ /*Vector layout*/ - int interpolation,nodesperelement,size; - int rank_interpolation=-1,rank_nodesperelement=-1; + int interpolation,nodesperelement,size,nlines,ncols,array_size; + int rank_interpolation=-1,rank_nodesperelement=-1,rank_arraysize=-1,max_rank_arraysize=0; + bool isarray=false; /*Get interpolation (and compute input if necessary)*/ for(int j=0;jSize();j++){ Element* element=xDynamicCast(this->elements->GetObjectByOffset(j)); - element->ResultInterpolation(&rank_interpolation,&rank_nodesperelement,output_enum); + element->ResultInterpolation(&rank_interpolation,&rank_nodesperelement,&rank_arraysize,output_enum); + if (rank_arraysize>max_rank_arraysize)max_rank_arraysize=rank_arraysize; } + rank_arraysize=max_rank_arraysize; /*Broadcast for cpus that do not have any elements*/ ISSM_MPI_Reduce(&rank_interpolation,&interpolation,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm()); ISSM_MPI_Reduce(&rank_nodesperelement,&nodesperelement,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm()); + ISSM_MPI_Reduce(&rank_arraysize,&array_size,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm()); ISSM_MPI_Bcast(&interpolation,1,ISSM_MPI_INT,0,IssmComm::GetComm()); ISSM_MPI_Bcast(&nodesperelement,1,ISSM_MPI_INT,0,IssmComm::GetComm()); - + ISSM_MPI_Bcast(&array_size,1,ISSM_MPI_INT,0,IssmComm::GetComm()); + if(results_on_nodes){ @@ -1080,19 +1085,38 @@ /*Allocate vector depending on interpolation*/ switch(interpolation){ - case P0Enum: size = this->elements->NumberOfElements(); break; - case P1Enum: size = this->vertices->NumberOfVertices(); break; + case P0Enum: isarray = false; size = this->elements->NumberOfElements(); break; + case P1Enum: isarray = false; size = this->vertices->NumberOfVertices(); break; + case P0ArrayEnum: isarray = true; nlines = this->elements->NumberOfElements(); ncols= array_size; break; default: _error_("Interpolation "< *vector_result = new Vector(size); - - /*Fill in vector*/ - for(int j=0;jSize();j++){ - Element* element=(Element*)elements->GetObjectByOffset(j); - element->ResultToVector(vector_result,output_enum); + if (!isarray){ + Vector *vector_result = new Vector(size); + + /*Fill in vector*/ + for(int j=0;jSize();j++){ + Element* element=(Element*)elements->GetObjectByOffset(j); + element->ResultToVector(vector_result,output_enum); + } + vector_result->Assemble(); + + if(save_results)results->AddResult(new GenericExternalResult*>(results->Size()+1,output_enum,vector_result,step,time)); } - vector_result->Assemble(); - - if(save_results)results->AddResult(new GenericExternalResult*>(results->Size()+1,output_enum,vector_result,step,time)); + else{ + IssmDouble* values = xNewZeroInit(nlines*ncols); + IssmDouble* allvalues = xNew(nlines*ncols); + + /*Fill-in matrix*/ + for(int j=0;jSize();j++){ + Element* element=xDynamicCast(this->elements->GetObjectByOffset(j)); + element->ResultToMatrix(values,ncols, output_enum); + } + /*Gather from all cpus*/ + ISSM_MPI_Allreduce((void*)values,(void*)allvalues,ncols*nlines,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm()); + xDelete(values); + + if(save_results)results->AddResult(new GenericExternalResult(results->Size()+1,output_enum,allvalues,nlines,ncols,step,time)); + xDelete(allvalues); + } } isvec = true; Index: /issm/trunk-jpl/src/c/classes/Inputs/BoolInput.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/BoolInput.h (revision 19553) +++ /issm/trunk-jpl/src/c/classes/Inputs/BoolInput.h (revision 19554) @@ -38,4 +38,5 @@ int GetResultInterpolation(void){return P0Enum;}; int GetResultNumberOfNodes(void){return 1;}; + int GetResultArraySize(void){return 1;}; void ResultToPatch(IssmDouble* values,int nodesperelement,int sid){_error_("not supported yet");}; void Configure(Parameters* parameters); Index: /issm/trunk-jpl/src/c/classes/Inputs/ControlInput.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/ControlInput.h (revision 19553) +++ /issm/trunk-jpl/src/c/classes/Inputs/ControlInput.h (revision 19554) @@ -79,4 +79,5 @@ int GetResultInterpolation(void); int GetResultNumberOfNodes(void); + int GetResultArraySize(void){return 1;}; void ResultToPatch(IssmDouble* values,int nodesperelement,int sid){_error_("not supported yet");}; void GetGradient(Vector* gradient_vec,int* doflist); Index: /issm/trunk-jpl/src/c/classes/Inputs/DatasetInput.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/DatasetInput.h (revision 19553) +++ /issm/trunk-jpl/src/c/classes/Inputs/DatasetInput.h (revision 19554) @@ -74,4 +74,5 @@ int GetResultInterpolation(void){_error_("not implemented yet");}; int GetResultNumberOfNodes(void){_error_("not implemented yet");}; + int GetResultArraySize(void){_error_("not implemented yet");}; void ResultToPatch(IssmDouble* values,int nodesperelement,int sid){_error_("not supported yet");}; void GetGradient(Vector* gradient_vec,int* doflist){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Inputs/DoubleArrayInput.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/DoubleArrayInput.cpp (revision 19554) +++ /issm/trunk-jpl/src/c/classes/Inputs/DoubleArrayInput.cpp (revision 19554) @@ -0,0 +1,111 @@ +/*!\file DoubleArrayInput.c + * \brief: implementation of the DoubleArrayInput object + */ + +#ifdef HAVE_CONFIG_H + #include +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "../classes.h" +#include "../../shared/shared.h" + +/*DoubleArrayInput constructors and destructor*/ +DoubleArrayInput::DoubleArrayInput(){/*{{{*/ + return; +} +/*}}}*/ +DoubleArrayInput::DoubleArrayInput(int in_enum_type,IssmDouble* in_values, int in_m){/*{{{*/ + + enum_type=in_enum_type; + m=in_m; + values=xNew(m); + xMemCpy(values,in_values,m); + +} +/*}}}*/ +DoubleArrayInput::~DoubleArrayInput(){/*{{{*/ + return; +} +/*}}}*/ + +/*Object virtual functions definitions:*/ +void DoubleArrayInput::Echo(void){/*{{{*/ + this->DeepEcho(); +} +/*}}}*/ +void DoubleArrayInput::DeepEcho(void){/*{{{*/ + + _printf_(setw(15)<<" DoubleArrayInput "<enum_type)<<" Size: " << m << "\n"); + for (int i=0;ivalues[i]<<"\n"); + +} +/*}}}*/ +int DoubleArrayInput::Id(void){ return -1; }/*{{{*/ +/*}}}*/ +int DoubleArrayInput::ObjectEnum(void){/*{{{*/ + + return DoubleArrayInputEnum; + +} +/*}}}*/ +Object* DoubleArrayInput::copy() {/*{{{*/ + + return new DoubleArrayInput(this->enum_type,this->values,this->m); + +} +/*}}}*/ +void DoubleArrayInput::Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){ /*{{{*/ + + MARSHALLING_ENUM(DoubleArrayInputEnum); + + MARSHALLING(enum_type); + MARSHALLING(m); + MARSHALLING_DYNAMIC(this->values,IssmDouble,m); +} +/*}}}*/ + +/*DoubleArrayInput management*/ +int DoubleArrayInput::InstanceEnum(void){/*{{{*/ + + return this->enum_type; + +} +/*}}}*/ +void DoubleArrayInput::ResultToMatrix(IssmDouble* values,int ncols,int sid){/*{{{*/ + + int ncols_local = this->GetResultArraySize(); + + /*Some checks*/ + _assert_(values); + _assert_(ncols_local<=ncols); + + /*Fill in arrays*/ + for(int i=0;ivalues[i]; +} +/*}}}*/ +void DoubleArrayInput::GetValues(IssmDouble** pvalues, int *pm){ /*{{{*/ + + /*output: */ + IssmDouble* outvalues= NULL; + + outvalues=xNew(m); + + xMemCpy(outvalues,values,m); + + /*assign output pointers: */ + *pm=m; + *pvalues=outvalues; +} +/*}}}*/ + +/*Object functions*/ +void DoubleArrayInput::ChangeEnum(int newenumtype){/*{{{*/ + this->enum_type=newenumtype; +} +/*}}}*/ +void DoubleArrayInput::Configure(Parameters* parameters){/*{{{*/ + /*do nothing: */ +} +/*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Inputs/DoubleArrayInput.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/DoubleArrayInput.h (revision 19554) +++ /issm/trunk-jpl/src/c/classes/Inputs/DoubleArrayInput.h (revision 19554) @@ -0,0 +1,76 @@ +/*! \file DoubleArrayInput.h + * \brief: header file for vector type input object + */ + +#ifndef _DOUBLE_ARRAY_INPUT_H_ +#define _DOUBLE_ARRAY_INPUT_H_ + +/*Headers:*/ +/*{{{*/ +#include "./Input.h" +/*}}}*/ + +class DoubleArrayInput: public Input{ + + public: + int enum_type; + IssmDouble* values; /*vector*/ + int m; /*size of vector*/ + + /*DoubleArrayInput constructors, destructors: {{{*/ + DoubleArrayInput(); + DoubleArrayInput(int enum_type,IssmDouble* values, int m); + ~DoubleArrayInput(); + /*}}}*/ + /*Object virtual functions definitions:{{{ */ + void Echo(); + void DeepEcho(); + int Id(); + int ObjectEnum(); + Object* copy(); + void Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction); + /*}}}*/ + /*DoubleArrayInput management: {{{*/ + int InstanceEnum(); + Input* SpawnTriaInput(int index1,int index2,int index3){_error_("not implemented yet");}; + Input* SpawnSegInput(int index1,int index2){_error_("not implemented yet");}; + Input* PointwiseDivide(Input* inputB){_error_("not implemented yet");}; + Input* PointwiseMin(Input* inputB){_error_("not implemented yet");}; + Input* PointwiseMax(Input* inputB){_error_("not implemented yet");}; + int GetResultInterpolation(void){return P0ArrayEnum;}; + int GetResultNumberOfNodes(void){return 1;}; + int GetResultArraySize(void){return m;}; + void ResultToMatrix(IssmDouble* values,int ncols,int sid); + void Configure(Parameters* parameters); + void GetValues(IssmDouble** pvalues,int* pm); + /*}}}*/ + /*numerics: {{{*/ + void GetInputValue(bool* pvalue){_error_("not implemented yet");}; + void GetInputValue(int* pvalue){_error_("not implemented yet");}; + void GetInputValue(IssmDouble* pvalue){_error_("not implemented yet");}; + void GetInputValue(IssmDouble* pvalue,Gauss* gauss){_error_("not implemented yet");}; + void GetInputValue(IssmDouble* pvalue,Gauss* gauss,IssmDouble time){_error_("not implemented yet");}; + void GetInputValue(IssmDouble* pvalue,Gauss* gauss ,int index){_error_("not implemented yet");}; + void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");}; + void GetInputAverage(IssmDouble* pvalue){_error_("not implemented yet");}; + void GetInputAllTimeAverages(IssmDouble** pvalues,IssmDouble** ptimes, int* pnumtimes){_error_("not implemented yet");}; + void GetInputUpToCurrentTimeAverages(IssmDouble** pvalues, IssmDouble** ptimes, int* pnumtimes, IssmDouble currenttime){_error_("not implemented yet");}; + void SquareMin(IssmDouble* psquaremin,Parameters* parameters){_error_("not implemented yet");}; + void ConstrainMin(IssmDouble minimum){_error_("not implemented yet");}; + void Set(IssmDouble setvalue){_error_("Set not implemented yet");}; + void Scale(IssmDouble scale_factor){_error_("not implemented yet");}; + void AXPY(Input* xinput,IssmDouble scalar){_error_("not implemented yet");}; + void Constrain(IssmDouble cm_min, IssmDouble cm_max){_error_("not implemented yet");}; + void ChangeEnum(int newenumtype); + IssmDouble InfinityNorm(void){_error_("not implemented yet");}; + IssmDouble Max(void){_error_("not implemented yet");}; + IssmDouble MaxAbs(void){_error_("not implemented yet");}; + IssmDouble Min(void){_error_("not implemented yet");}; + IssmDouble MinAbs(void){_error_("not implemented yet");}; + void Extrude(int start){_error_("not supported yet");}; + void VerticallyIntegrate(Input* thickness_input){_error_("not implemented yet");}; + void GetVectorFromInputs(Vector* vector,int* doflist){_error_("not implemented yet");}; + /*}}}*/ + +}; +#endif /* _DOUBLE_ARRAY_INPUT_H */ Index: /issm/trunk-jpl/src/c/classes/Inputs/DoubleInput.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/DoubleInput.h (revision 19553) +++ /issm/trunk-jpl/src/c/classes/Inputs/DoubleInput.h (revision 19554) @@ -41,4 +41,5 @@ int GetResultInterpolation(void){return P0Enum;}; int GetResultNumberOfNodes(void){return 1;}; + int GetResultArraySize(void){return 1;}; void ResultToPatch(IssmDouble* values,int nodesperelement,int sid){_error_("not supported yet");}; void AddTimeValues(IssmDouble* values,int step,IssmDouble time){_error_("not supported yet");}; Index: /issm/trunk-jpl/src/c/classes/Inputs/Input.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/Input.h (revision 19553) +++ /issm/trunk-jpl/src/c/classes/Inputs/Input.h (revision 19554) @@ -61,5 +61,7 @@ virtual int GetResultInterpolation(void)=0; virtual int GetResultNumberOfNodes(void)=0; - virtual void ResultToPatch(IssmDouble* values,int nodesperelement,int sid){_error_("not supported yet");}; + virtual int GetResultArraySize(void)=0; + virtual void ResultToPatch(IssmDouble* values,int nodesperelement,int sid){_error_("not supported yet");}; + virtual void ResultToMatrix(IssmDouble* values,int ncols,int sid){_error_("not supported yet");}; }; #endif Index: /issm/trunk-jpl/src/c/classes/Inputs/IntInput.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/IntInput.h (revision 19553) +++ /issm/trunk-jpl/src/c/classes/Inputs/IntInput.h (revision 19554) @@ -42,4 +42,5 @@ int GetResultInterpolation(void){return P0Enum;}; int GetResultNumberOfNodes(void){return 1;}; + int GetResultArraySize(void){return 1;}; void ResultToPatch(IssmDouble* values,int nodesperelement,int sid){_error_("not supported yet");}; void AddTimeValues(IssmDouble* values,int step,IssmDouble time){_error_("not supported yet");}; Index: /issm/trunk-jpl/src/c/classes/Inputs/PentaInput.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/PentaInput.h (revision 19553) +++ /issm/trunk-jpl/src/c/classes/Inputs/PentaInput.h (revision 19554) @@ -43,4 +43,5 @@ int GetResultInterpolation(void); int GetResultNumberOfNodes(void); + int GetResultArraySize(void){return 1;}; void ResultToPatch(IssmDouble* values,int nodesperelement,int sid); void AddTimeValues(IssmDouble* values,int step,IssmDouble time){_error_("not supported yet");}; Index: /issm/trunk-jpl/src/c/classes/Inputs/SegInput.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/SegInput.h (revision 19553) +++ /issm/trunk-jpl/src/c/classes/Inputs/SegInput.h (revision 19554) @@ -43,4 +43,5 @@ int GetResultInterpolation(void){_error_("not implemented");}; int GetResultNumberOfNodes(void){_error_("not implemented");}; + int GetResultArraySize(void){_error_("not implemented");}; void ResultToPatch(IssmDouble* values,int nodesperelement,int sid){_error_("not supported yet");}; void AddTimeValues(IssmDouble* values,int step,IssmDouble time){_error_("not supported yet");}; Index: /issm/trunk-jpl/src/c/classes/Inputs/TetraInput.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/TetraInput.h (revision 19553) +++ /issm/trunk-jpl/src/c/classes/Inputs/TetraInput.h (revision 19554) @@ -43,4 +43,5 @@ int GetResultInterpolation(void); int GetResultNumberOfNodes(void); + int GetResultArraySize(void){return 1;}; void ResultToPatch(IssmDouble* values,int nodesperelement,int sid); void AddTimeValues(IssmDouble* values,int step,IssmDouble time){_error_("not supported yet");}; Index: /issm/trunk-jpl/src/c/classes/Inputs/TransientInput.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/TransientInput.cpp (revision 19553) +++ /issm/trunk-jpl/src/c/classes/Inputs/TransientInput.cpp (revision 19554) @@ -357,4 +357,9 @@ } /*}}}*/ +int TransientInput::GetResultArraySize(void){/*{{{*/ + + return 1; +} +/*}}}*/ void TransientInput::Extrude(int start){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/Inputs/TransientInput.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/TransientInput.h (revision 19553) +++ /issm/trunk-jpl/src/c/classes/Inputs/TransientInput.h (revision 19554) @@ -48,4 +48,5 @@ int GetResultInterpolation(void); int GetResultNumberOfNodes(void); + int GetResultArraySize(void); void ResultToPatch(IssmDouble* values,int nodesperelement,int sid){_error_("not supported yet");}; void Configure(Parameters* parameters); Index: /issm/trunk-jpl/src/c/classes/Inputs/TriaInput.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/TriaInput.cpp (revision 19553) +++ /issm/trunk-jpl/src/c/classes/Inputs/TriaInput.cpp (revision 19554) @@ -141,4 +141,10 @@ } /*}}}*/ +int TriaInput::GetResultArraySize(void){/*{{{*/ + + return 1; + +} +/*}}}*/ void TriaInput::ResultToPatch(IssmDouble* values,int nodesperelement,int sid){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/Inputs/TriaInput.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/TriaInput.h (revision 19553) +++ /issm/trunk-jpl/src/c/classes/Inputs/TriaInput.h (revision 19554) @@ -43,4 +43,5 @@ int GetResultInterpolation(void); int GetResultNumberOfNodes(void); + int GetResultArraySize(void); void ResultToPatch(IssmDouble* values,int nodesperelement,int sid); void AddTimeValues(IssmDouble* values,int step,IssmDouble time){_error_("not supported yet");}; Index: /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp (revision 19553) +++ /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp (revision 19554) @@ -39,4 +39,7 @@ dpermil=0; + albedo_snow=0; + albedo_ice=0; + sediment_compressibility=0; sediment_porosity=0; @@ -95,4 +98,8 @@ case SMBforcingEnum: /*Nothing to add*/ + break; + case SMBgembEnum: + iomodel->Constant(&this->albedo_ice,SmbAIceEnum); + iomodel->Constant(&this->albedo_snow,SmbASnowEnum); break; case SMBpddEnum: Index: /issm/trunk-jpl/src/c/classes/Materials/Matpar.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matpar.h (revision 19553) +++ /issm/trunk-jpl/src/c/classes/Materials/Matpar.h (revision 19554) @@ -35,4 +35,8 @@ IssmDouble rlapslgm; IssmDouble dpermil; + + /*albedo: */ + IssmDouble albedo_ice; + IssmDouble albedo_snow; /*hydrology Dual Porous Continuum: */ Index: /issm/trunk-jpl/src/c/classes/classes.h =================================================================== --- /issm/trunk-jpl/src/c/classes/classes.h (revision 19553) +++ /issm/trunk-jpl/src/c/classes/classes.h (revision 19554) @@ -62,4 +62,5 @@ #include "./Inputs/BoolInput.h" #include "./Inputs/DoubleInput.h" +#include "./Inputs/DoubleArrayInput.h" #include "./Inputs/IntInput.h" #include "./Inputs/TetraInput.h" Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp (revision 19554) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp (revision 19554) @@ -0,0 +1,1806 @@ +/*!\file GEMB module from Alex Gardner. + * \brief: calculates SMB + */ + +#include "./SurfaceMassBalancex.h" +#include "../../shared/shared.h" +#include "../../toolkits/toolkits.h" + +const double Pi =3.141592653589793238462643383279502884197169399375105820974944592308; + +void Gembx(FemModel* femmodel){ /*{{{*/ + + for(int i=0;ielements->Size();i++){ + Element* element=xDynamicCast(femmodel->elements->GetObjectByOffset(i)); + element->SmbGemb(); + } + +} /*}}}*/ +void GembgridInitialize(IssmDouble** pdz, int* psize, IssmDouble zTop, IssmDouble dzTop, IssmDouble zMax, IssmDouble zY){ /*{{{*/ + + /* This file sets up the initial grid spacing and total grid depth. The + grid structure is set as constant grid length 'dzTop' for the top + 'zTop' meters of the model grid. Bellow 'zTop' the gid length increases + linearly with depth */ + + /*intermediary:*/ + IssmDouble dgpTop; + int gpTop, gpBottom; + int i; + IssmDouble gp0,z0; + IssmDouble* dzT=NULL; + IssmDouble* dzB=NULL; + + /*output: */ + int m; + IssmDouble* dz=NULL; + + //----------------------Calculate Grid Lengths------------------------------ + //calculate number of top grid points + dgpTop = zTop/dzTop; + + //check to see if the top grid cell structure length (dzTop) goes evenly + //into specified top structure depth (zTop). Also make sure top grid cell + //structure length (dzTop) is greater than 5 cm + if (dgpTop != round(dgpTop)){ + _error_("top grid cell structure length does not go evenly into specified top structure depth, adjust dzTop or zTop"); + } + else if(dzTop < 0.05){ + _printf_("initial top grid cell length (dzTop) is < 0.05 m"); + } + gpTop=round(dgpTop); + + //initialize top grid depth vector + dzT = xNew(gpTop); + for (i=0;i z0){ + gp0= gp0 * zY; + z0 = z0 + gp0; + gpBottom++; + } + //initialize bottom vectors + dzB = xNewZeroInit(gpBottom); + gp0 = dzTop; + z0 = zTop; + for(i=0;i(gpTop+gpBottom); + for(i=0;i 400 kg m-3 because (H is set to zero)------ + // ---------------this is a major limitation of the model------------------- + + // initialize + IssmDouble F = 0, H=0, G=0, E; + + // convert T from K to ºC + T = T - 273.15; + + // temperature coefficient F + if(T>-6.0) F = 0.7 + ((T/-6) * 0.3); + if(T<=-6.0 && T>-22.0) F = 1 - ((T+6)/-16 * 0.8); + if(T<=-22.0 && T>=-40.0) F = 0.2 - ((T+22)/-18 * 0.2); + + // density coefficient F + if(d<150) H=1.0; + + if(d>=150 & d<400) H = 1 - ((d-150)/250); + + // temperature gradient coefficient G + if(dT >= 0.16 && dT < 0.25) G = ((dT - 0.16)/0.09) * 0.1; + if(dT >= 0.25 && dT < 0.4) G = 0.1 + (((dT - 0.25)/0.15) * 0.57); + if(dT >= 0.4 && dT < 0.5) G = 0.67 + (((dT - 0.4)/0.1) * 0.23); + if(dT >= 0.5 && dT < 0.7) G = 0.9 + (((dT - 0.5)/0.2) * 0.1); + if(dT >= .7 ) G = 1.0; + + // model time growth constant E + E = 0.09; //[mm d-1] + + // grouped coefficient Q + return F*H*G*E; + +} /*}}}*/ +void grainGrowth(IssmDouble* re, IssmDouble* gdn, IssmDouble* gsp, IssmDouble* T,IssmDouble* dz,IssmDouble* d, IssmDouble* W,IssmDouble smb_dt,int m,int aIdx){ /*{{{*/ + + /*Created by: Alex S. Gardner, University of Alberta + + *Description*: models the effective snow grain size + + *Reference:* + DENDRITIC SNOW METAMORPHISM: + Brun, E., P. David, M. Sudul, and G. Brunot, 1992: A numerical model to + simulate snow-cover stratigraphy for operational avalanche forecasting. + Journal of Glaciology, 38, 13-22. + + NONDENDRITIC SNOW METAMORPHISM: + Dry snow metamorphism: + Marbouty, D., 1980: An experimental study of temperature-gradient + metamorphism. Journal of Glaciology, 26, 303-312. + + WET SNOW METAMORPHISM: + Brun, E., 1989: Investigation on wet-snow metamorphism in respect of + liquid-water content. Annals of Glaciology, 13, 22-26. + + INPUTS + * T: grid cell temperature [k] + * dz: grid cell depth [m] + * d: grid cell density [kg m-3] + * W: water content [kg/m^2] + * re: effective grain size [mm] + * gdn: grain dentricity + * gsp: grain sphericity + * dt: time step of input data [s] + + OUTPUTS + * re: effective grain size [mm] + * gdn: grain dentricity + * gsp: grain sphericity*/ + + /*intermediary: */ + IssmDouble dt; + IssmDouble* gsz=NULL; + IssmDouble* dT=NULL; + IssmDouble* zGPC=NULL; + IssmDouble* lwc=NULL; + IssmDouble Q=0; + + /*only when aIdx = 1 or 2 do we run grainGrowth: */ + if(aIdx!=1 & aIdx!=2){ + /*come out as we came in: */ + return; + } + + /*Figure out grain size from effective grain radius: */ + gsz=xNew(m); for(int i=0;i(m); for(int i=0;i9) lwc[i]=9.0; + + + /* Calculate temperature gradiant across grid cells + * Returns the avereage gradient across the upper and lower grid cell */ + + //initialize + dT=xNewZeroInit(m); + + //depth of grid point center from surface + zGPC=xNewZeroInit(m); + for(int i=0;i0){ + + //_printf_("Dendritic dry snow metamorphism\n"); + + //index for dentricity > 0 and T gradients < 5 degC m-1 and >= 5 degC m-1 + if(dT[i]<5){ + //determine coefficients + IssmDouble A = - 2e8 * exp(-6e3 / T[i]) * dt; + IssmDouble B = 1e9 * exp(-6e3 / T[i]) * dt; + //new dentricity and sphericity for dT < 5 degC m-1 + gdn[i] += A; + gsp[i] += B; + } + else{ + // new dendricity and sphericity for dT >= 5 degC m-1 + + //determine coefficients + IssmDouble C = (-2e8 * exp(-6e3 / T[i]) * dt) * pow(dT[i],.4); + gdn[i] +=C; + gsp[i] +=C; + } + + // wet snow metamorphism + if(W[i]>0){ + + //_printf_("D}ritic wet snow metamorphism\n"); + + //determine coefficient + IssmDouble D = (1/16) * pow(lwc[i],3.0) * dt; + + // new dendricity and sphericity for wet snow + gdn[i] -= D; + gsp[i] += D; + } + + // dendricity and sphericity can not be > 1 or < 0 + if (gdn[i]<0)gdn[i]=0.0; + if (gsp[i]<0)gsp[i]=0.0; + if (gdn[i]>1)gdn[i]=1.0; + if (gsp[i]>1)gsp[i]=1.0; + + // determine new grain size (mm) + gsz[i] = 0.1 + (1-gdn[i])*0.25 + (0.5-gsp[i])*0.1; + + } + else{ + + /*Dry snow metamorphism (Marbouty, 1980) grouped model coefficinets + *from Marbouty, 1980: Figure 9*/ + + //_printf_("Nond}ritic snow metamorphism\n"); + Q = Marbouty(T[i],d[i],dT[i]); + + // calculate grain growth + gsz[i] += Q* dt; + + //Wet snow metamorphism (Brun, 1989) + if(W[i]>0){ + //_printf_("Nond}ritic wet snow metamorphism\n"); + //wet rate of change coefficient + IssmDouble E = 1.28e-8 + (4.22e-10 * pow(lwc[i],3.0))* (dt *86400); // [mm^3 s^-1] + + // calculate change in grain volume and convert to grain size + gsz[i] = 2 * pow(3/(Pi * 4)*((4 / 3)*Pi*pow(gsz[i]/2,3.0) + E),1.0/3.0); + + } + + // grains with sphericity == 1 can not have grain sizes > 2 mm (Brun, 1992) + if (gsp[i]==1.0 && gsz[i]>2) gsz[i]=2.0; + + // grains with sphericity == 0 can not have grain sizes > 5 mm (Brun, 1992) + if (gsp[i]!=1.0 && gsz[i]>5.0) gsz[i]=5.0; + } + + //convert grain size back to effective grain radius: + re[i]=gsz[i]/2; + } + + /*Free ressources:*/ + xDelete(gsz); + xDelete(dT); + xDelete(zGPC); + xDelete(lwc); + +} /*}}}*/ +void albedo(IssmDouble* a, int aIdx, IssmDouble* re, IssmDouble* d, IssmDouble cldFrac, IssmDouble aIce, IssmDouble aSnow, IssmDouble* TK, IssmDouble* W, IssmDouble P, IssmDouble EC, IssmDouble t0wet, IssmDouble t0dry, IssmDouble K, IssmDouble dt, int m) { /*{{{*/ + + //// Calculates Snow, firn and ice albedo as a function of: + // 1 : effective grain radius (Gardner & Sharp, 2009) + // 2 : effective grain radius (Brun et al., 2009) + // 3 : density and cloud amount (Greuell & Konzelmann, 1994) + // 4 : exponential time decay & wetness (Bougamont & Bamber, 2005) + + //// Inputs + // aIdx = albedo method to use + + // Methods 1 & 2 + // re = surface effective grain radius [mm] + + // Methods 3 + // d = snow surface density [kg m-3] + // n = cloud amount + // aIce = albedo of ice + // aSnow = albedo of fresh snow + + // Methods 4 + // aIce = albedo of ice + // aSnow = albedo of fresh snow + // a = grid cell albedo from prevous time step; + // T = grid cell temperature [k] + // W = pore water [kg] + // P = precipitation [mm w.e.] or [kg m-3] + // EC = surface evaporation (-) condensation (+) [kg m-2] + // t0wet = time scale for wet snow (15-21.9) [d] + // t0dry = warm snow timescale [15] [d] + // K = time scale temperature coef. (7) [d] + // dt = time step of input data [s] + + //// Output + // a = grid cell albedo + + //// Usage + // Method 1 + // a = albedo(1, 0.1); + + // Method 4 + // a = albedo(4, [], [], [], 0.48, 0.85, [0.8 0.5 ... 0.48], ... + // [273 272.5 ... 265], [0 0.001 ... 0], 0, 0.01, 15, 15, 7, 3600) + + //some constants: + const IssmDouble dSnow = 300; // density of fresh snow [kg m-3] + const IssmDouble dIce = 910; // density of ice [kg m-3] + + if(aIdx==1){ + //function of effective grain radius + + //convert effective radius to specific surface area [cm2 g-1] + IssmDouble S = 3.0 / (.091 * re[0]); + + //determine broadband albedo + a[0]= 1.48 - pow(S,-.07); + } + else if(aIdx==2){ + + // Spectral fractions (Lefebre et al., 2003) + // [0.3-0.8um 0.8-1.5um 1.5-2.8um] + + IssmDouble sF[3] = {0.606, 0.301, 0.093}; + + // convert effective radius to grain size in meters + IssmDouble gsz = (re[0] * 2) / 1000.0; + + // spectral range: + // 0.3 - 0.8um + IssmDouble a0 = fmin(0.98, 1 - 1.58 *pow(gsz,0.5)); + // 0.8 - 1.5um + IssmDouble a1 = fmax(0, 0.95 - 15.4 *pow(gsz,0.5)); + // 1.5 - 2.8um + IssmDouble a2 = fmax(0.127, 0.88 + 346.3*gsz - 32.31*pow(gsz,0.5)); + + // broadband surface albedo + a[0] = sF[0]*a0 + sF[1]*a1 + sF[2]*a2; + + } + else if(aIdx==3){ + + // a as a function of density + + // calculate albedo + a[0] = aIce + (d[0] - dIce)*(aSnow - aIce) / (dSnow - dIce) + (0.05 * (cldFrac - 0.5)); + } + else if(aIdx==4){ + + // exponential time decay & wetness + + // change in albedo with time: + // (d_a) = (a - a_old)/(t0) + // where: t0 = timescale for albedo decay + + dt = dt / 86400; // convert from [s] to [d] + + // initialize variables + IssmDouble* t0=xNew(m); + IssmDouble* T=xNew(m); + IssmDouble* t0warm=xNew(m); + IssmDouble* d_a=xNew(m); + + // specify constants + // a_wet = 0.15; // water albedo (0.15) + // a_new = aSnow // new snow albedo (0.64 - 0.89) + // a_old = aIce; // old snow/ice albedo (0.27-0.53) + // t0_wet = t0wet; // time scale for wet snow (15-21.9) [d] + // t0_dry = t0dry; // warm snow timescale [15] [d] + // K = 7 // time scale temperature coef. (7) [d] + // W0 = 300; // 200 - 600 [mm] + const IssmDouble z_snow = 15; // 16 - 32 [mm] + + // determine timescale for albedo decay + for(int i=0;i0)t0[i]=t0wet; // wet snow timescale + for(int i=0;i=-10)t0[i]= t0warm[i]; + for(int i=0;i 0 && T[0] < 0) P = P + (EC/dSnow) * 1000; // add cond to precip [mm] + + a[0] = aSnow - (aSnow - a[0]) * exp(-P/z_snow); + + //----------THIS NEEDS TO BE IMPLEMENTED AT A LATER DATE------------ + // modification of albedo due to thin layer of water on the surface + // a_surf = a_wet - (a_wet - a_surf) * exp(-W_surf/W0); + + /*Free ressources:*/ + xDelete(t0); + xDelete(T); + xDelete(t0warm); + xDelete(d_a); + + } + else _error_("albedo method switch should range from 1 to 4!"); + + // Check for erroneous values + if (a[0] > 1) _printf_("albedo > 1.0\n"); + else if (a[0] < 0) _printf_("albedo is negative\n"); + else if (isnan(a[0])) _error_("albedo == NAN\n"); +} /*}}}*/ +void thermo(IssmDouble* pEC, IssmDouble* T, IssmDouble* dz, IssmDouble* d, IssmDouble* swf, IssmDouble dlwrf, IssmDouble Ta, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble Ws, IssmDouble dt0, int m, IssmDouble Vz, IssmDouble Tz) { /*{{{*/ + + /* ENGLACIAL THERMODYNAMICS*/ + + /* Description: + computes new temperature profile accounting for energy absorption and + thermal diffusion.*/ + + // INPUTS + // T: grid cell temperature [k] + // dz: grid cell depth [m] + // d: grid cell density [kg m-3] + // swf: shortwave radiation fluxes [W m-2] + // dlwrf: downward longwave radiation fluxes [W m-2] + // Ta: 2 m air temperature + // V: wind velocity [m s-1] + // eAir: screen level vapor pressure [Pa] + // Ws: surface water content [kg] + // dt0: time step of input data [s] + // elev: surface elevation [m a.s.l.] + // Vz: air temperature height above surface [m] + // Tz: wind height above surface [m] + + // OUTPUTS + // T: grid cell temperature [k] + // EC: evaporation/condensation [kg] + + /*intermediary: */ + IssmDouble* K = NULL; + IssmDouble* KU = NULL; + IssmDouble* KD = NULL; + IssmDouble* KP = NULL; + IssmDouble* Au = NULL; + IssmDouble* Ad = NULL; + IssmDouble* Ap = NULL; + IssmDouble* Nu = NULL; + IssmDouble* Nd = NULL; + IssmDouble* Np = NULL; + IssmDouble* dzU = NULL; + IssmDouble* dzD = NULL; + IssmDouble* sw = NULL; + IssmDouble* dT_sw = NULL; + IssmDouble* lw = NULL; + IssmDouble* T0 = NULL; + IssmDouble* Tu = NULL; + IssmDouble* Td = NULL; + + IssmDouble z0; + IssmDouble dt; + IssmDouble max_fdt=0; + IssmDouble Ts; + IssmDouble L; + IssmDouble eS; + IssmDouble Ri; + IssmDouble coefM; + IssmDouble coefH; + IssmDouble An; + IssmDouble C; + IssmDouble shf; + IssmDouble SB; + IssmDouble CI; + IssmDouble ds; + IssmDouble dAir; + IssmDouble TCs; + IssmDouble lhf; + IssmDouble EC_day; + IssmDouble dT_turb; + IssmDouble turb; + IssmDouble ulw; + IssmDouble dT_ulw; + IssmDouble dlw; + IssmDouble dT_dlw; + + /*outputs:*/ + IssmDouble EC; + + // INITIALIZE + CI = 2102; // heat capacity of snow/ice (J kg-1 k-1) + // CA = 1005; // heat capacity of air (J kg-1 k-1) + // LF = 0.3345E6; // latent heat of fusion(J kg-1) + // LV = 2.495E6; // latent heat of vaporization(J kg-1) + // dIce = 910; // density of ice [kg m-3] + // dSnow = 300; // density of snow [kg m-3] + SB = 5.67E-8; // Stefan-Boltzmann constant [W m-2 K-4] + + ds = d[0]; // density of top grid cell + + // calculated air density [kg/m3] + dAir = 0.029 * pAir /(8.314 * Ta); + + // thermal capacity of top grid cell [J/k] + TCs = d[0]*dz[0]*CI; + + //initialize Evaporation - Condenstation + EC = 0; + + // check if all SW applied to surface or distributed throught subsurface + // swIdx = length(swf) > 1 + + // SURFACE ROUGHNESS (Bougamont, 2006) + // wind/temperature surface roughness height [m] + if (ds < 910 && Ws == 0) z0 = 0.00012; // 0.12 mm for dry snow + else if (ds >= 910) z0 = 0.0032; // 3.2 mm for ice + else z0 = 0.0013; // 1.3 mm for wet snow + + // if V = 0 goes to infinity therfore if V = 0 change + if(V<.01)V=.01; + + // Bulk-transfer coefficient for turbulent fluxes + An = pow(0.4,2) / pow(log(Tz/z0),2); // Bulk-transfer coefficient + C = An * dAir * V; // shf & lhf common coefficient + + // THERMAL CONDUCTIVITY (Sturm, 1997: J. Glaciology) + // calculate new K profile [W m-1 K-1] + + // initialize conductivity + K= xNewZeroInit(m); + + // for snow and firn (density < 910 kg m-3) (Sturn et al, 1997) + for(int i=0;i= 910 kg m-3) + for(int i=0;i=910) K[i] = 9.828 * exp(-5.7E-3*T[i]); + + // THERMAL DIFFUSION COEFFICIENTS + + // A discretization scheme which truncates the Taylor-Series expansion + // after the 3rd term is used. See Patankar 1980, Ch. 3&4 + + // discretized heat equation: + + // Tp = (Au*Tu° + Ad*Td° + (Ap-Au-Ad)Tp° + S) / Ap + + // where neighbor coefficients Au, Ap, & Ad are + + // Au = [dz_u/2KP + dz_p/2KE]^-1 + // Ad = [dz_d/2KP + dz_d/2KD]^-1 + // Ap = d*CI*dz/Dt + + // and u & d represent grid points up and down from the center grid point + // p and // u & d represent grid points up and down from the center grid + // point p and ° identifies previous time step values. S is a source term. + + // u, d, and p conductivities + KU = xNew(m); + KD = xNew(m); + KP = xNew(m); + + KU[0] = UNDEF; + KD[m-1] = UNDEF; + for(int i=1;i(m); + dzD = xNew(m); + dzU[0]=UNDEF; + dzD[m-1]=UNDEF; + + for(int i=1;i 1/2) [s] + dt=1e12; + for(int i=0;i 1/2 + // must go evenly into one hour or the data frequency if it is smaller + + // all integer factors of the number of second in a day (8600 [s]) + int f[45] = {1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24, 25, 30, 36, 40, 45, 48, 50, 60, + 72, 75, 80, 90, 100, 120, 144, 150, 180, 200, 225, 240, 300, 360, 400, 450, 600, 720, 900, 1200, 1800, 3600}; + + // return the min integer factor that is < dt + max_fdt=0; + for(int i=0;i<45;i++){ + if (f[i]=max_fdt)max_fdt=f[i]; + } + dt=max_fdt; + + // determine mean (harmonic mean) of K/dz for u, d, & p + Au = xNew(m); + Ad = xNew(m); + Ap = xNew(m); + for(int i=0;i(m); + Nd = xNew(m); + Np = xNew(m); + for(int i=0;i(m); + for(int i=0;i(m); + for(int i=0;i(m+2); + Tu=xNew(m); + Td=xNew(m); + + /* CALCULATE ENERGY SOURCES AND DIFFUSION FOR EVERY TIME STEP [dt]*/ + for (int i=1;i<=dt0;i+=dt){ + + // PART OF ENERGY CONSERVATION CHECK + // store initial temperature + //T_init = T; + + // calculate temperature of snow surface (Ts) + // when incoming SW radition is allowed to penetrate the surface, + // the modeled energy balance becomes very sensitive to how Ts is + // calculated. The estimated enegy balance & melt are significanly + // less when Ts is taken as the mean of the x top grid cells. + Ts = (T[0] + T[1])/2; + Ts = fmin(273.15,Ts); // don't allow Ts to exceed 273.15 K (0°C) + + //TURBULENT HEAT FLUX + + // Monin–Obukhov Stability Correction + // Reference: + // Ohmura, A., 1982: Climate and Energy-Balance on the Arctic Tundra. + // Journal of Climatology, 2, 65-84. + + // calculate the Bulk Richardson Number (Ri) + Ri = (2*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 + Ri = fmin(Ri, 0.19); + + // calculate momentum 'coefM' stability factor + if (Ri > 0){ + // if stable + coefM = 1/(1-5.2*Ri); + } + else { + coefM =pow (1-18*Ri,-0.25); + } + + // calculate heat/wind 'coef_H' stability factor + if (Ri < -0.03) coefH = 1.3 * coefM; + else coefH = coefM; + + //// Sensible Heat + // calculate the sensible heat flux [W m-2](Patterson, 1998) + shf = C * 1005 * (Ta - Ts); + + // adjust using Monin–Obukhov stability theory + shf = shf / (coefM * coefH); + + //// Latent Heat + // determine if snow pack is melting & calcualte surface vapour pressure over ice or liquid water + if (Ts >= 273.15){ + L = 2.495E6; + + // for an ice surface Murphy and Koop, 2005 [Equation 7] + eS = exp(9.550426 - 5723.265/Ts + 3.53068 * log(Ts) - 0.00728332 * Ts); + } + else{ + L = 2.8295E6; // latent heat of sublimation for liquid surface (assume liquid on surface when Ts == 0 deg C) + // Wright (1997), US Meteorological Handbook from Murphy and Koop, 2005 Appendix A + eS = 611.21 * exp(17.502 * (Ts - 273.15) / (240.97 + Ts - 273.15)); + } + + // Latent heat flux [W m-2] + lhf = C * L * (eAir - eS) * 0.622 / pAir; + + // adjust using Monin–Obukhov stability theory (if lhf '+' then there is energy and mass gained at the surface, + // if '-' then there is mass and energy loss at the surface. + lhf = lhf / (coefM * coefH); + + //mass loss (-)/acreation(+) due to evaporation/condensation [kg] + EC_day = lhf * 86400 / L; + + // temperature change due turbulent fluxes + turb = (shf + lhf)* dt; + dT_turb = turb / TCs; + + // upward longwave contribution + ulw = - SB * pow(Ts,4.0) * dt; + dT_ulw = ulw / TCs; + + // new grid point temperature + + //SW penetrates surface + for(int i=0;i 1E-6 || isnan(E_diff) + disp(T(1)) + _error_("energy not conserved in thermodynamics equations"); + */ + } + + /*Free ressources:*/ + xDelete(K); + xDelete(KU); + xDelete(KD); + xDelete(KP); + xDelete(Au); + xDelete(Ad); + xDelete(Ap); + xDelete(Nu); + xDelete(Nd); + xDelete(Np); + xDelete(dzU); + xDelete(dzD); + xDelete(sw); + xDelete(dT_sw); + xDelete(lw); + xDelete(T0); + xDelete(Tu); + xDelete(Td); + + + /*Assign output pointers:*/ + *pEC=EC; + +} /*}}}*/ +void shortwave(IssmDouble** pswf, int swIdx, int aIdx, IssmDouble dsw, IssmDouble as, IssmDouble* d, IssmDouble* dz, IssmDouble* re, int m){ /*{{{*/ + + // DISTRIBUTES ABSORBED SHORTWAVE RADIATION WITHIN SNOW/ICE + + // swIdx = 0 : all absorbed SW energy is assigned to the top grid cell + + // swIdx = 1 : absorbed SW is distributed with depth as a function of: + // 1 : snow density (taken from Bassford, 2004) + // 2 : grain size in 3 spectral bands (Brun et al., 1992) + + // Inputs + // swIdx = shortwave allowed to penetrate surface (0 = No, 1 = Yes) + // aIdx = method for calculating albedo (1-4) + // dsw = downward shortwave radiative flux [w m-2] + // as = surface albedo + // d = grid cell density [kg m-3] + // dz = grid cell depth [m] + // re = grid cell effective grain radius [mm] + + // Outputs + // swf = absorbed shortwave radiation [W m-2] + + /*outputs:*/ + IssmDouble* swf = NULL; + + // SHORTWAVE FUNCTION + // initialize variables + swf = xNewZeroInit(m); + + if (swIdx == 0) {// all sw radation is absorbed in by the top grid cell + + // calculate surface shortwave radiation fluxes [W m-2] + swf[0] = (1 - as) * dsw; + } + else{ // sw radation is absorbed at depth within the glacier + + if (aIdx == 2){ // function of effective radius (3 spectral bands) + + IssmDouble * gsz=NULL; + IssmDouble * B1_cum=NULL; + IssmDouble * B2_cum=NULL; + IssmDouble* h =NULL; + IssmDouble* B1 =NULL; + IssmDouble* B2 =NULL; + IssmDouble* exp1 = NULL; + IssmDouble* exp2 = NULL; + IssmDouble* Qs1 = NULL; + IssmDouble* Qs2 = NULL; + + // convert effective radius [mm] to grain size [m] + gsz=xNew(m); + for(int i=0;i(m+1); + B2_cum=xNew(m+1); + for(int i=0;i(m); + B1 =xNew(m); + B2 =xNew(m); + for(int i=0;i(m); + exp2 = xNew(m); + for(int i=0;i(m+1); + Qs2 = xNew(m+1); + for(int i=0;i(gsz); + xDelete(B1_cum); + xDelete(B2_cum); + xDelete(h); + xDelete(B1); + xDelete(B2); + xDelete(exp1); + xDelete(exp2); + xDelete(Qs1); + xDelete(Qs2); + } + else{ //function of grid cell density + + /*intermediary: */ + IssmDouble* B_cum = NULL; + IssmDouble* exp_B = NULL; + IssmDouble* Qs = NULL; + IssmDouble* B = NULL; + + // fraction of sw radiation absorbed in top grid cell (wavelength > 0.8um) + IssmDouble SWs = 0.36; + + // 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 + // respectively. This, however, may not account for the fact that + // the albedo of wavelengths > 800 nm has a much lower albedo. + + // calculate surface shortwave radiation fluxes [W m-2] + IssmDouble swf_s = SWs * (1 - as) * dsw; + + // calculate surface shortwave radiation fluxes [W m-2] + IssmDouble swf_ss = (1-SWs) * (1 - as) * dsw; + + // SW allowed to penetrate into snowpack + IssmDouble Bs = 10; // snow SW extinction coefficient [m-1] (Bassford,2006) + IssmDouble Bi = 1.3; // ice SW extinction coefficient [m-1] (Bassford,2006) + + // calculate extinction coefficient B [m-1] vector + B=xNew(m); + for(int i=0;i(m+1); + exp_B = xNew(m); + for(int i=0;i(m+1); + for(int i=0;i(B_cum); + xDelete(exp_B); + xDelete(Qs); + xDelete(B); + } + } + /*Assign output pointers: */ + *pswf=swf; +} /*}}}*/ +void accumulation(IssmDouble** pT, IssmDouble** pdz, IssmDouble** pd, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pm, IssmDouble T_air, IssmDouble P, IssmDouble dzMin, IssmDouble aSnow){ /*{{{*/ + + // Adds precipitation and deposition to the model grid + + // Author: Alex Gardner, University of Alberta + // Date last modified: JAN, 2008 + + /* Description: + adjusts the properties of the top grid cell to account for accumulation + T_air & T = Air and top grid cell temperatures [K] + dz = topgrid cell length [m] + d = density of top grid gell [kg m-3] + P = precipitation [mm w.e.] or [kg m-3] + re = effective grain radius [mm] + gdn = grain dentricity + gsp = grain sphericity*/ + + // MAIN FUNCTION + // specify constants + const IssmDouble dIce = 910; // density of ice [kg m-3] + const IssmDouble dSnow = 150; // density of snow [kg m-3] + const IssmDouble reNew = 0.1; // new snow grain size [mm] + const IssmDouble gdnNew = 1; // new snow dendricity + const IssmDouble gspNew = 0.5; // new snow sphericity + + /*intermediary: */ + IssmDouble* mInit=NULL; + bool top=true; + IssmDouble mass, massinit, mass_diff; + + /*output: */ + IssmDouble* T=NULL; + IssmDouble* dz=NULL; + IssmDouble* d=NULL; + IssmDouble* W=NULL; + IssmDouble* a=NULL; + IssmDouble* re=NULL; + IssmDouble* gdn=NULL; + IssmDouble* gsp=NULL; + int m; + + /*Recover pointers: */ + T=*pT; + dz=*pdz; + d=*pd; + W=*pW; + a=*pa; + re=*pre; + gdn=*pgdn; + gsp=*pgsp; + m=*pm; + + /*Allocate: */ + mInit=xNew(m); + + if (P > 0){ + // determine initial mass + for(int i=0;i dzMin){ + + newcell(&T,T_air,top,m); //new cell T + newcell(&dz,z_snow,top,m); //new cell dz + newcell(&d,dSnow,top,m); //new cell d + newcell(&W,0,top,m); //new cell W + newcell(&a,aSnow,top,m); //new cell a + newcell(&re,reNew,top,m); //new cell grain size + newcell(&gdn,gdnNew,top,m); //new cell grain dendricity + newcell(&gsp,gspNew,top,m); //new cell grain sphericity + m=m+1; + } + else { // if snow depth is less than specified minimum dz snow + + IssmDouble mass = mInit[0] + P; // grid cell adjust mass + + dz[0] = dz[0] + P/dSnow; // adjust grid cell depth + d[0] = mass / dz[0]; // adjust grid cell density + + // adjust variables as a linearly weighted function of mass + // adjust temperature (assume P is same temp as air) + T[0] = (T_air * P + T[0] * mInit[0])/mass; + + // adjust a, re, gdn & gsp + a[0] = (aSnow * P + a[0] * mInit[0])/mass; + re[0] = (reNew * P + re[0] * mInit[0])/mass; + gdn[0] = (gdnNew * P + gdn[0] * mInit[0])/mass; + gsp[0] = (gspNew * P + gsp[0] * mInit[0])/mass; + } + } + else{ // if rain + + /*rain is added by increasing the mass and temperature of the ice + of the top grid cell. Temperatures are set artifically high to + account for the latent heat of fusion. This is the same as + directly adding liquid water to the the snow pack surface but + makes the numerics easier.*/ + + IssmDouble LF = 0.3345E6; // latent heat of fusion(J kg-1) + IssmDouble CI = 2102; // specific heat capacity of snow/ice (J kg-1 k-1) + + // grid cell adjust mass + mass = mInit[0] + P; + + // adjust temperature + // liquid: must account for latent heat of fusion + T[0] = (P *(T_air + LF/CI) + T[0] * mInit[0]) / mass; + + // adjust grid cell density + d[0] = mass / dz[0]; + + // if d > the density of ice, d = dIce + if (d[0] > dIce){ + d[0] = dIce; // adjust d + dz[0] = mass / d[0]; // dz is adjusted to conserve mass + } + } + + // check for conservation of mass + mass=0; for(int i=0;i 0) _error_("mass not conserved in accumulation function"); + + } + /*Free ressources:*/ + if(mInit)xDelete(mInit); + + /*Assign output pointers:*/ + *pT=T; + *pdz=dz; + *pd=d; + *pW=W; + *pa=a; + *pre=re; + *pgdn=gdn; + *pgsp=gsp; + *pm=m; +} /*}}}*/ +void melt(IssmDouble* pM, IssmDouble* pR, IssmDouble* pmAdd, IssmDouble** pT, IssmDouble** pd, IssmDouble** pdz, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pn, IssmDouble dzMin, IssmDouble zMax, IssmDouble zMin){ /*{{{*/ + + //// MELT ROUTINE + + // Description: + // computes the quantity of meltwater due to snow temperature in excess of + // 0 deg C, determines pore water content and adjusts grid spacing + + /*intermediary:*/ + IssmDouble* m=NULL; + IssmDouble* maxF=NULL; + IssmDouble* dW=NULL; + IssmDouble* exsW=NULL; + IssmDouble* exsT=NULL; + IssmDouble* surpT=NULL; + IssmDouble* surpE=NULL; + IssmDouble* F=NULL; + IssmDouble* flxDn=NULL; + IssmDouble* R=NULL; + IssmDouble* ER=NULL; + IssmDouble* EI=NULL; + IssmDouble* EW=NULL; + IssmDouble* M=NULL; + int* D=NULL; + + IssmDouble sumM; + IssmDouble sumER; + IssmDouble addE; + IssmDouble mSum0; + IssmDouble sumE0; + IssmDouble mSum1; + IssmDouble sumE1; + IssmDouble dE; + IssmDouble dm; + IssmDouble X; + IssmDouble Wi; + int D_size; + + /*outputs:*/ + IssmDouble mAdd; + IssmDouble Rsum; + IssmDouble* T=*pT; + IssmDouble* d=*pd; + IssmDouble* dz=*pdz; + IssmDouble* W=*pW; + IssmDouble* a=*pa; + IssmDouble* re=*pre; + IssmDouble* gdn=*pgdn; + IssmDouble* gsp=*pgsp; + int n=*pn; + + //// INITIALIZATION + + /*Allocations: */ + M=xNewZeroInit(n); + maxF=xNew(n); + dW=xNew(n); + + // specify constants + const IssmDouble CtoK = 273.15; // clecius to Kelvin conversion + const IssmDouble CI = 2102; // specific heat capacity of snow/ice (J kg-1 k-1) + const IssmDouble LF = 0.3345E6; // latent heat of fusion(J kg-1) + const IssmDouble dPHC = 830; // pore hole close off density[kg m-3] + const IssmDouble dIce = 910; // density of ice [kg m-3] + + // store initial mass [kg] and energy [J] + m=xNew(n); for(int i=0;i(n); for(int i=0;i(n); for(int i=0;i(n); + for(int i=0;i 0){ + // _printf_("PORE WATER REFREEZE"); + // calculate maximum freeze amount, maxF [kg] + for(int i=0;idIce)d[i]=m[i]/d[i]; + } + + // squeeze water from snow pack + exsW=xNew(n); for(int i=0;i 0) || (cellsum(exsW,n) > 0)){ + + // _printf_(""MELT OCCURS"); + // check to see if thermal energy exceeds energy to melt entire cell + // if so redistribute temperature to lower cells (temperature surplus) + // (maximum T of snow before entire grid cell melts is a constant + // LF/CI = 159.1342) + surpT=xNew(n); for(int i=0;i 0 ){ + // _printf_("T Surplus"); + // calculate surplus energy + surpE=xNew(n); for(int i=0;i 0){ + // use surplus energy to increase the temperature of lower cell + T[i+i] = surpE[i] * m[i+1]/CI + T[i+i]; + surpT[i+1] = fmax(0, (T[i+1] - CtoK - 159.1342)); + surpE[i+1] = surpT[i+1] * CI / m[i+1]; + + // adjust current cell properties (again 159.1342 is the max T) + T[i] = CtoK + 159.1342; + surpE[i] = 0; + i = i + 1; + } + // recalculate temperature excess above 0 deg C + for(int i=0;i(n); + R = xNewZeroInit(n); + for(int i=0;i(n+1); for(int i=0;i=0;i--){ + if(M[i]>0 || exsW[i]){ + X=i; + break; + } + } + + //// meltwater percolation + for(int i=0;i mw_depth + if (inM == 0 && i > X){ + break; + } + + // if reaches impermeable ice layer all liquid water runs off (R) + else if (d[i] >= dIce){ // dPHC = pore hole close off [kg m-3] + // _printf_("ICE LAYER"); + // no water freezes in this cell + // no water percolates to lower cell + // cell ice temperature & density do not change + + m[i] = m[i] - M[i]; // mass after melt + Wi = (910-d[i]) * Swi * (m[i]/d[i]); // irreducible water + dW[i] = fmin(inM, Wi - W[i]); // change in pore water + R[i] = fmax(0, inM - dW[i]); // runoff + } + // check if no energy to refreeze meltwater + else if (maxF[i] == 0){ + // _printf_("REFREEZE == 0"); + // no water freezes in this cell + // cell ice temperature & density do not change + + m[i] = m[i] - M[i]; // mass after melt + Wi = (910-d[i]) * Swi * (m[i]/d[i]); // irreducible water + dW[i] = fmin(inM, Wi-W[i]); // change in pore water + flxDn[i+1] = fmax(0, inM-dW[i]); // meltwater out + F[i] = 0; // no freeze + } + // some or all meltwater refreezes + else{ + // change in density density and temperature + // _printf_("MELT REFREEZE"); + //-----------------------melt water----------------------------- + IssmDouble dz_0 = m[i]/d[i]; + IssmDouble dMax = (dIce - d[i])*dz_0; // d max = dIce + IssmDouble F1 = fmin(fmin(inM,dMax),maxF[i]); // maximum refreeze + m[i] = m[i] + F1; // mass after refreeze + d[i] = m[i]/dz_0; + + //-----------------------pore water----------------------------- + Wi = (910-d[i])* Swi * dz_0; // irreducible water + dW[i] = fmin(inM - F1, Wi-W[i]); // change in pore water + if (-dW[i]>W[i] ){ + dW[i]= W[i]; + } + IssmDouble F2 = 0; + + if (dW[i] < 0){ // excess pore water + dMax = (dIce - d[i])*dz_0; // maximum refreeze + IssmDouble maxF2 = fmin(dMax, maxF[i]-F1); // maximum refreeze + F2 = fmin(-dW[i], maxF2); // pore water refreeze + m[i] = m[i] + F2; // mass after refreeze + d[i] = m[i]/dz_0; + } + + flxDn[i+1] = inM - F1 - dW[i] - F2; // meltwater out + T[i] = T[i] + ((F1+F2)*(LF+(CtoK - T[i])*CI)/(m[i]*CI));// change in temperature + + + // check if an ice layer forms + if (d[i] == dIce){ + // _printf_("ICE LAYER FORMS"); + // excess water runs off + R[i] = flxDn[i+1]; + // no water percolates to lower cell + flxDn[i+1] = 0; + } + } + } + + + //// GRID CELL SPACING AND MODEL DEPTH + for(int i=0;i(D_size); + D_size=0; for(int i=0;i=0;i--){ + if(dz[i] dzMin * 2') + // adjust variables as a linearly weighted function of mass + IssmDouble m_new = m[i] + m[i+1]; + T[i+1] = (T[i]*m[i] + T[i+1]*m[i+1]) / m_new; + a[i+1] = (a[i]*m[i] + a[i+1]*m[i+1]) / m_new; + re[i+1] = (re[i]*m[i] + re[i+1]*m[i+1]) / m_new; + gdn[i+1] = (gdn[i]*m[i] + gdn[i+1]*m[i+1]) / m_new; + gsp[i+1] = (gsp[i]*m[i] + gsp[i+1]*m[i+1]) / m_new; + + // merge with underlying grid cell and delete old cell + dz [i+1] = dz[i] + dz[i+1]; // combine cell depths + d[i+1] = m_new / dz[i+1]; // combine top densities + W[i+1] = W[i+1] + W[i]; // combine liquid water + m[i+1] = m_new; // combine top masses + + // set cell to 99999 for deletion + m[i] = 99999; + } + } + + // delete combined cells + xDelete(D); D_size=0; for(int i=0;i(D_size); + D_size=0; for(int i=0;i zMax){ // check if model is too deep + // _printf_("z > zMax') + // // mass and energy loss + // mAdd = -(m(end) + W(end)); + // addE = -(T(end) * m(end) * CI); + // + // // add a grid cell of the same size and temperature to the bottom + // dz(end) = []; T(end) = []; W(end) = []; m(end) = []; + // d(end) = []; a(end) = []; re(end) = []; gdn(end) = []; + // gsp(end) = []; + // } + + //// CHECK FOR MASS AND ENERGY CONSERVATION + + // calculate final mass [kg] and energy [J] + ER=xNew(n); + Rsum = cellsum(R,n); + for(int i=0;i(m); + if(EI)xDelete(EI); + if(maxF)xDelete(maxF); + if(dW)xDelete(dW); + if(exsW)xDelete(exsW); + if(exsT)xDelete(exsT); + if(surpT)xDelete(surpT); + if(surpE)xDelete(surpE); + if(F)xDelete(F); + if(flxDn)xDelete(flxDn); + if(D)xDelete(D); + if(R)xDelete(R); + if(M)xDelete(M); + + /*Assign output pointers:*/ + *pM=sumM; + *pR=Rsum; + *pmAdd=mAdd; + + *pT=T; + *pd=d; + *pdz=dz; + *pW=W; + *pa=a; + *pre=re; + *pgdn=gdn; + *pgsp=gsp; + *pn=n; + +} /*}}}*/ +void densification(IssmDouble* d,IssmDouble* dz, IssmDouble* T, IssmDouble* re, int denIdx, IssmDouble C, IssmDouble dt, IssmDouble Tmean,int m){ /*{{{*/ + + //// THIS NEEDS TO BE DOUBLE CHECKED AS THERE SEAMS TO BE LITTLE DENSIFICATION IN THE MODEL OUTOUT [MAYBE COMPATION IS COMPNSATED FOR BY TRACES OF SNOW???] + + //// FUNCTION INFO + + // Author: Alex Gardner, University of Alberta + // Date last modified: FEB, 2008 + + // Description: + // computes the densification of snow/firn using the emperical model of + // Herron and Langway (1980) or the semi-emperical model of Anthern et al. + // (2010) + + // Inputs: + // denIdx = densification model to use: + // 1 = emperical model of Herron and Langway (1980) + // 2 = semi-imerical model of Anthern et al. (2010) + // 3 = physical model from Appendix B of Anthern et al. (2010) + // d = initial snow/firn density [kg m-3] + // T = temperature [K] + // dz = grid cell size [m] + // C = average accumulation rate [kg m-2 yr-1] + // dt = time lapsed [s] + // re = effective grain radius [mm]; + // Ta = mean annual temperature + + // Reference: + // Herron and Langway (1980), Anthern et al. (2010) + + //// FOR TESTING + // denIdx = 2; + // d = 800; + // T = 270; + // dz = 0.005; + // C = 200; + // dt = 60*60; + // re = 0.7; + // Tmean = 273.15-18; + + //// MAIN FUNCTION + // specify constants + const IssmDouble dIce = 910; // density of ice [kg m-3] + dt = dt / 86400; // convert from [s] to [d] + // R = 8.314 // gas constant [mol-1 K-1] + // Ec = 60 // activation energy for self-diffusion of water + // // molecules through the ice tattice [kJ mol-1] + // Eg = 42.4 // activation energy for grain growth [kJ mol-1] + + /*intermediary: */ + IssmDouble c0,c1,H; + + // initial mass + IssmDouble* mass_init = xNew(m);for(int i=0;i(m-1); + cumdz[0]=dz[0]; + for(int i=1;i(m); + obp[0]=0; + for(int i=1;i 550 [kg m-3] + + + for(int i=0;idIce)d[i]=dIce; + + // calculate new grid cell length + dz[i] = mass_init[i] / d[i]; + } + /*Free ressources:*/ + xDelete(mass_init); + xDelete(cumdz); + xDelete(obp); + +} /*}}}*/ +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){ /*{{{*/ + + //// TURBULENT HEAT FLUX + + // Description: + // computed the surface sensible and latent heat fluxes [W m-2], this + // function also calculated the mass loss/acreation due to + // condensation/evaporation [kg] + + // Reference: + // Dingman, 2002. + + //// INPUTS: + // Ta: 2m air temperature [K] + // Ts: snow/firn/ice surface temperature [K] + // V: wind speed [m s^-^1] + // eAir: screen level vapor pressure [Pa] + // pAir: surface pressure [Pa] + // ds: surface density [kg/m^3] + // Ws: surface liquid water content [kg/m^2] + // Vz: height above ground at which wind (V) eas sampled [m] + // Tz: height above ground at which temperature (T) was sampled [m] + + //// FUNCTION INITILIZATION + + // CA = 1005; // heat capacity of air (J kg-1 k-1) + // LF = 0.3345E6; // latent heat of fusion(J kg-1) + // LV = 2.495E6; // latent heat of vaporization(J kg-1) + // dIce = 910; // density of ice [kg m-3] + + /*intermediary:*/ + IssmDouble d_air; + IssmDouble Ri; + IssmDouble z0; + IssmDouble coef_M,coef_H; + IssmDouble An, C; + IssmDouble L, eS; + + /*output: */ + IssmDouble shf, lhf, EC; + + // calculated air density [kg/m3] + d_air = 0.029 * pAir /(8.314 * Ta); + + //// Determine Surface Roughness + // Bougamont, 2006 + // wind/temperature surface roughness height [m] + if (ds < 910 && Ws == 0) z0 = 0.00012; // 0.12 mm for dry snow + else if (ds >= 910) z0 = 0.0032; // 3.2 mm for ice + else z0 = 0.0013; // 1.3 mm for wet snow + + //// Monin–Obukhov Stability Correction + // Reference: + // Ohmura, A., 1982: Climate and Energy-Balance on the Arctic Tundra. + // Journal of Climatology, 2, 65-84. + + // if V = 0 goes to infinity therfore if V = 0 change + if(V< .01) V=.01; + + // calculate the Bulk Richardson Number (Ri) + Ri = (2*9.81* (Vz - z0) * (Ta - Ts)) / ((Ta + Ts)* pow(V,2)); + + // calculate Monin–Obukhov stability factors 'coef_M' and 'coef_H' + + // do not allow Ri to exceed 0.19 + if(Ri>.19)Ri= 0.19; + + // calculate momentum 'coef_M' stability factor + if (Ri > 0) coef_M = pow(1-5.2*Ri,-1); // if stable + else coef_M = pow(1-18*Ri,-0.25); + + // calculate heat/wind 'coef_H' stability factor + if (Ri < -0.03) coef_H = 1.3 * coef_M; + else coef_H = coef_M; + + //// Bulk-transfer coefficient + An = pow(0.4,2) / pow(log(Tz/z0),2); // Bulk-transfer coefficient + C = An * d_air * V; // shf & lhf common coefficient + + //// Sensible Heat + // calculate the sensible heat flux [W m-2](Patterson, 1998) + shf = C * 1005 * (Ta - Ts); + + // adjust using Monin–Obukhov stability theory + shf = shf / (coef_M * coef_H); + + //// Latent Heat + // determine if snow pack is melting & calcualte surface vapour pressure + // over ice or liquid water + if (Ts >= 273.15){ + L = 2.495E6; + + // for an ice surface Murphy and Koop, 2005 [Equation 7] + eS = exp(9.550426 - 5723.265/Ts + 3.53068 * log(Ts)- 0.00728332 * Ts); + } + else{ + L = 2.8295E6; // latent heat of sublimation + // for liquid surface (assume liquid on surface when Ts == 0 deg C) + // Wright (1997), US Meteorological Handbook from Murphy and Koop, + // 2005 Apendix A + eS = 611.21 * exp(17.502 * (Ts - 273.15) / (240.97 + Ts - 273.15)); + } + + // Latent heat flux [W m-2] + lhf = C * L * (eAir - eS) * 0.622 / pAir; + + // adjust using Monin–Obukhov stability theory (if lhf '+' then there is + // energy and mass gained at the surface, if '-' then there is mass and + // energy loss at the surface. + lhf = lhf / (coef_M * coef_H); + + // mass loss (-)/acreation(+) due to evaporation/condensation [kg] + EC = lhf * 86400 / L; + + /*assign output poitners: */ + *pshf=shf; + *plhf=lhf; + *pEC=EC; + +} /*}}}*/ Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h (revision 19553) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h (revision 19554) @@ -17,5 +17,17 @@ void SmbHenningx(FemModel* femmodel); void SmbComponentsx(FemModel* femmodel); -void SmbMeltComponentsx(FemModel* femmodel); +void SmbMeltComponentsx(FemModel* femmodel); +/*GEMB: */ +void Gembx(FemModel* femmodel); +void GembgridInitialize(IssmDouble** pdz, int* psize, IssmDouble zTop, IssmDouble dzTop, IssmDouble zMax, IssmDouble zY); +IssmDouble Marbouty(IssmDouble T, IssmDouble d, IssmDouble dT); +void grainGrowth(IssmDouble* pre, IssmDouble* pgdn, IssmDouble* pgsp, IssmDouble* T,IssmDouble* dz,IssmDouble* d, IssmDouble* W,IssmDouble smb_dt,int m,int aIdx); +void albedo(IssmDouble* a,int aIdx, IssmDouble* re, IssmDouble* d, IssmDouble cldFrac, IssmDouble aIce, IssmDouble aSnow, IssmDouble* T, IssmDouble* W, IssmDouble P, IssmDouble EC, IssmDouble t0wet, IssmDouble t0dry, IssmDouble K, IssmDouble dt,int m); +void shortwave(IssmDouble** pswf, int swIdx, int aIdx, IssmDouble dsw, IssmDouble as, IssmDouble* d, IssmDouble* dz, IssmDouble* re, int m); +void thermo(IssmDouble* pEC, IssmDouble* T, IssmDouble* dz, IssmDouble* d, IssmDouble* swf, IssmDouble dlw, IssmDouble Ta, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble Ws, IssmDouble dt0, int m, IssmDouble Vz, IssmDouble Tz); +void accumulation(IssmDouble** pT, IssmDouble** pdz, IssmDouble** pd, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pm, IssmDouble Ta, IssmDouble P, IssmDouble dzMin, IssmDouble aSnow); +void melt(IssmDouble* pM, IssmDouble* pR, IssmDouble* pmAdd, IssmDouble** pT, IssmDouble** pd, IssmDouble** pdz, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pn, IssmDouble dzMin, IssmDouble zMax, IssmDouble zMin); +void densification(IssmDouble* d,IssmDouble* dz, IssmDouble* T, IssmDouble* re, int denIdx, IssmDouble C, IssmDouble dt, IssmDouble Tmean,int m); +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); #endif /* _SurfaceMassBalancex_H*/ Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 19553) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 19554) @@ -354,4 +354,5 @@ SmbNumRequestedOutputsEnum, SmbRequestedOutputsEnum, + SmbIsInitializedEnum, /*SMBforcing*/ SMBforcingEnum, @@ -369,6 +370,7 @@ SmbCEnum, SmbTzEnum, - SmbVzEnum, - SmbSpinUpEnum, + SmbVzEnum, + SmbDtEnum, + SmbDzEnum, SmbAIdxEnum, SmbSwIdxEnum, @@ -387,4 +389,20 @@ SmbT0dryEnum, SmbKEnum, + SmbDEnum, + SmbReEnum, + SmbGdnEnum, + SmbGspEnum, + SmbECEnum, + SmbWEnum, + SmbAEnum, + SmbTEnum, + SmbIsgraingrowthEnum, + SmbIsalbedoEnum, + SmbIsshortwaveEnum, + SmbIsthermalEnum, + SmbIsaccumulationEnum, + SmbIsmeltEnum, + SmbIsdensificationEnum, + SmbIsturbulentfluxEnum, /*SMBpdd*/ SMBpddEnum, @@ -520,4 +538,5 @@ DatasetInputEnum, DoubleInputEnum, + DoubleArrayInputEnum, DataSetParamEnum, DoubleMatArrayParamEnum, @@ -696,4 +715,5 @@ /*Element Interpolations{{{*/ P0Enum, + P0ArrayEnum, P1Enum, P1DGEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 19553) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 19554) @@ -358,4 +358,5 @@ case SmbNumRequestedOutputsEnum : return "SmbNumRequestedOutputs"; case SmbRequestedOutputsEnum : return "SmbRequestedOutputs"; + case SmbIsInitializedEnum : return "SmbIsInitialized"; case SMBforcingEnum : return "SMBforcing"; case SmbMassBalanceEnum : return "SmbMassBalance"; @@ -372,5 +373,6 @@ case SmbTzEnum : return "SmbTz"; case SmbVzEnum : return "SmbVz"; - case SmbSpinUpEnum : return "SmbSpinUp"; + case SmbDtEnum : return "SmbDt"; + case SmbDzEnum : return "SmbDz"; case SmbAIdxEnum : return "SmbAIdx"; case SmbSwIdxEnum : return "SmbSwIdx"; @@ -389,4 +391,20 @@ case SmbT0dryEnum : return "SmbT0dry"; case SmbKEnum : return "SmbK"; + case SmbDEnum : return "SmbD"; + case SmbReEnum : return "SmbRe"; + case SmbGdnEnum : return "SmbGdn"; + case SmbGspEnum : return "SmbGsp"; + case SmbECEnum : return "SmbEC"; + case SmbWEnum : return "SmbW"; + case SmbAEnum : return "SmbA"; + case SmbTEnum : return "SmbT"; + case SmbIsgraingrowthEnum : return "SmbIsgraingrowth"; + case SmbIsalbedoEnum : return "SmbIsalbedo"; + case SmbIsshortwaveEnum : return "SmbIsshortwave"; + case SmbIsthermalEnum : return "SmbIsthermal"; + case SmbIsaccumulationEnum : return "SmbIsaccumulation"; + case SmbIsmeltEnum : return "SmbIsmelt"; + case SmbIsdensificationEnum : return "SmbIsdensification"; + case SmbIsturbulentfluxEnum : return "SmbIsturbulentflux"; case SMBpddEnum : return "SMBpdd"; case SmbDelta18oEnum : return "SmbDelta18o"; @@ -512,4 +530,5 @@ case DatasetInputEnum : return "DatasetInput"; case DoubleInputEnum : return "DoubleInput"; + case DoubleArrayInputEnum : return "DoubleArrayInput"; case DataSetParamEnum : return "DataSetParam"; case DoubleMatArrayParamEnum : return "DoubleMatArrayParam"; @@ -680,4 +699,5 @@ case GiaWEnum : return "GiaW"; case P0Enum : return "P0"; + case P0ArrayEnum : return "P0Array"; case P1Enum : return "P1"; case P1DGEnum : return "P1DG"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 19553) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 19554) @@ -364,4 +364,5 @@ else if (strcmp(name,"SmbNumRequestedOutputs")==0) return SmbNumRequestedOutputsEnum; else if (strcmp(name,"SmbRequestedOutputs")==0) return SmbRequestedOutputsEnum; + else if (strcmp(name,"SmbIsInitialized")==0) return SmbIsInitializedEnum; else if (strcmp(name,"SMBforcing")==0) return SMBforcingEnum; else if (strcmp(name,"SmbMassBalance")==0) return SmbMassBalanceEnum; @@ -378,13 +379,14 @@ else if (strcmp(name,"SmbTz")==0) return SmbTzEnum; else if (strcmp(name,"SmbVz")==0) return SmbVzEnum; - else if (strcmp(name,"SmbSpinUp")==0) return SmbSpinUpEnum; + else if (strcmp(name,"SmbDt")==0) return SmbDtEnum; + else if (strcmp(name,"SmbDz")==0) return SmbDzEnum; else if (strcmp(name,"SmbAIdx")==0) return SmbAIdxEnum; else if (strcmp(name,"SmbSwIdx")==0) return SmbSwIdxEnum; - else if (strcmp(name,"SmbDenIdx")==0) return SmbDenIdxEnum; - else if (strcmp(name,"SmbZTop")==0) return SmbZTopEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"SmbDzTop")==0) return SmbDzTopEnum; + if (strcmp(name,"SmbDenIdx")==0) return SmbDenIdxEnum; + else if (strcmp(name,"SmbZTop")==0) return SmbZTopEnum; + else if (strcmp(name,"SmbDzTop")==0) return SmbDzTopEnum; else if (strcmp(name,"SmbDzMin")==0) return SmbDzMinEnum; else if (strcmp(name,"SmbZY")==0) return SmbZYEnum; @@ -398,4 +400,20 @@ else if (strcmp(name,"SmbT0dry")==0) return SmbT0dryEnum; else if (strcmp(name,"SmbK")==0) return SmbKEnum; + else if (strcmp(name,"SmbD")==0) return SmbDEnum; + else if (strcmp(name,"SmbRe")==0) return SmbReEnum; + else if (strcmp(name,"SmbGdn")==0) return SmbGdnEnum; + else if (strcmp(name,"SmbGsp")==0) return SmbGspEnum; + else if (strcmp(name,"SmbEC")==0) return SmbECEnum; + else if (strcmp(name,"SmbW")==0) return SmbWEnum; + else if (strcmp(name,"SmbA")==0) return SmbAEnum; + else if (strcmp(name,"SmbT")==0) return SmbTEnum; + else if (strcmp(name,"SmbIsgraingrowth")==0) return SmbIsgraingrowthEnum; + else if (strcmp(name,"SmbIsalbedo")==0) return SmbIsalbedoEnum; + else if (strcmp(name,"SmbIsshortwave")==0) return SmbIsshortwaveEnum; + else if (strcmp(name,"SmbIsthermal")==0) return SmbIsthermalEnum; + else if (strcmp(name,"SmbIsaccumulation")==0) return SmbIsaccumulationEnum; + else if (strcmp(name,"SmbIsmelt")==0) return SmbIsmeltEnum; + else if (strcmp(name,"SmbIsdensification")==0) return SmbIsdensificationEnum; + else if (strcmp(name,"SmbIsturbulentflux")==0) return SmbIsturbulentfluxEnum; else if (strcmp(name,"SMBpdd")==0) return SMBpddEnum; else if (strcmp(name,"SmbDelta18o")==0) return SmbDelta18oEnum; @@ -488,5 +506,8 @@ else if (strcmp(name,"MeshdeformationSolution")==0) return MeshdeformationSolutionEnum; else if (strcmp(name,"MeshdeformationAnalysis")==0) return MeshdeformationAnalysisEnum; - else if (strcmp(name,"LevelsetAnalysis")==0) return LevelsetAnalysisEnum; + else stage=5; + } + if(stage==5){ + if (strcmp(name,"LevelsetAnalysis")==0) return LevelsetAnalysisEnum; else if (strcmp(name,"LevelsetStabilization")==0) return LevelsetStabilizationEnum; else if (strcmp(name,"ExtrapolationAnalysis")==0) return ExtrapolationAnalysisEnum; @@ -506,8 +527,5 @@ else if (strcmp(name,"DataSet")==0) return DataSetEnum; else if (strcmp(name,"Constraints")==0) return ConstraintsEnum; - else stage=5; - } - if(stage==5){ - if (strcmp(name,"Loads")==0) return LoadsEnum; + else if (strcmp(name,"Loads")==0) return LoadsEnum; else if (strcmp(name,"Materials")==0) return MaterialsEnum; else if (strcmp(name,"Nodes")==0) return NodesEnum; @@ -524,4 +542,5 @@ else if (strcmp(name,"DatasetInput")==0) return DatasetInputEnum; else if (strcmp(name,"DoubleInput")==0) return DoubleInputEnum; + else if (strcmp(name,"DoubleArrayInput")==0) return DoubleArrayInputEnum; else if (strcmp(name,"DataSetParam")==0) return DataSetParamEnum; else if (strcmp(name,"DoubleMatArrayParam")==0) return DoubleMatArrayParamEnum; @@ -610,5 +629,8 @@ else if (strcmp(name,"MassFlux")==0) return MassFluxEnum; else if (strcmp(name,"MeltingOffset")==0) return MeltingOffsetEnum; - else if (strcmp(name,"Misfit")==0) return MisfitEnum; + else stage=6; + } + if(stage==6){ + if (strcmp(name,"Misfit")==0) return MisfitEnum; else if (strcmp(name,"Pressure")==0) return PressureEnum; else if (strcmp(name,"PressurePicard")==0) return PressurePicardEnum; @@ -629,8 +651,5 @@ else if (strcmp(name,"ThicknessAbsMisfit")==0) return ThicknessAbsMisfitEnum; else if (strcmp(name,"SurfaceAbsMisfit")==0) return SurfaceAbsMisfitEnum; - else stage=6; - } - if(stage==6){ - if (strcmp(name,"Vel")==0) return VelEnum; + else if (strcmp(name,"Vel")==0) return VelEnum; else if (strcmp(name,"Velocity")==0) return VelocityEnum; else if (strcmp(name,"VxAverage")==0) return VxAverageEnum; @@ -695,4 +714,5 @@ else if (strcmp(name,"GiaW")==0) return GiaWEnum; else if (strcmp(name,"P0")==0) return P0Enum; + else if (strcmp(name,"P0Array")==0) return P0ArrayEnum; else if (strcmp(name,"P1")==0) return P1Enum; else if (strcmp(name,"P1DG")==0) return P1DGEnum; @@ -732,5 +752,8 @@ else if (strcmp(name,"Outputdefinition4")==0) return Outputdefinition4Enum; else if (strcmp(name,"Outputdefinition5")==0) return Outputdefinition5Enum; - else if (strcmp(name,"Outputdefinition6")==0) return Outputdefinition6Enum; + else stage=7; + } + if(stage==7){ + if (strcmp(name,"Outputdefinition6")==0) return Outputdefinition6Enum; else if (strcmp(name,"Outputdefinition7")==0) return Outputdefinition7Enum; else if (strcmp(name,"Outputdefinition8")==0) return Outputdefinition8Enum; @@ -752,8 +775,5 @@ else if (strcmp(name,"Outputdefinition24")==0) return Outputdefinition24Enum; else if (strcmp(name,"Outputdefinition25")==0) return Outputdefinition25Enum; - else stage=7; - } - if(stage==7){ - if (strcmp(name,"Outputdefinition26")==0) return Outputdefinition26Enum; + else if (strcmp(name,"Outputdefinition26")==0) return Outputdefinition26Enum; else if (strcmp(name,"Outputdefinition27")==0) return Outputdefinition27Enum; else if (strcmp(name,"Outputdefinition28")==0) return Outputdefinition28Enum; @@ -855,5 +875,8 @@ else if (strcmp(name,"MaxVx")==0) return MaxVxEnum; else if (strcmp(name,"MaxAbsVx")==0) return MaxAbsVxEnum; - else if (strcmp(name,"MinVy")==0) return MinVyEnum; + else stage=8; + } + if(stage==8){ + if (strcmp(name,"MinVy")==0) return MinVyEnum; else if (strcmp(name,"MaxVy")==0) return MaxVyEnum; else if (strcmp(name,"MaxAbsVy")==0) return MaxAbsVyEnum; @@ -875,8 +898,5 @@ else if (strcmp(name,"None")==0) return NoneEnum; else if (strcmp(name,"AggressiveMigration")==0) return AggressiveMigrationEnum; - else stage=8; - } - if(stage==8){ - if (strcmp(name,"SoftMigration")==0) return SoftMigrationEnum; + else if (strcmp(name,"SoftMigration")==0) return SoftMigrationEnum; else if (strcmp(name,"SubelementMigration")==0) return SubelementMigrationEnum; else if (strcmp(name,"SubelementMigration2")==0) return SubelementMigration2Enum; Index: /issm/trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp (revision 19553) +++ /issm/trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp (revision 19554) @@ -557,2 +557,72 @@ for(int i=0;i<4;i++) X[i]=Ainv[i][0]*B[0] + Ainv[i][1]*B[1] + Ainv[i][2]*B[2] + Ainv[i][3]*B[3]; }/*}}}*/ + +void newcell(IssmDouble** pcell, IssmDouble newvalue, bool top, int m){ /*{{{*/ + + IssmDouble* cell=NULL; + IssmDouble* dummy=NULL; + + /*recover pointer: */ + cell=*pcell; + + /*reallocate:*/ + dummy=xNew(m+1); + + /*copy data:*/ + for(int i=0;i(cell); cell=dummy; + + /*assign output pointer:*/ + *pcell=cell; +} /*}}}*/ +IssmDouble cellsum(IssmDouble* cell, int m){ /*{{{*/ + + IssmDouble sum=0; + + for(int i=0;i(nind); + + for(int i=0;i(cell); + + /*assign output pointers: */ + *pcell=newcell; +} /*}}}*/ +void cellsplit(IssmDouble** pcell, int m, int i,IssmDouble scale) { /*{{{*/ + + /*input: */ + IssmDouble* cell=*pcell; + + /*output: */ + IssmDouble* newcell=xNew(m+1); + + for(int j=0;j(cell); + + /*assign output pointers: */ + *pcell=newcell; +} /*}}}*/ Index: /issm/trunk-jpl/src/c/shared/Matrix/matrix.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Matrix/matrix.h (revision 19553) +++ /issm/trunk-jpl/src/c/shared/Matrix/matrix.h (revision 19554) @@ -25,4 +25,8 @@ void Matrix4x4Determinant(IssmDouble* Adet,IssmDouble* A); void Matrix4x4Solve(IssmDouble* X,IssmDouble* A,IssmDouble* B); - + +void newcell(IssmDouble** pcell, IssmDouble newvalue, bool top, int m); +IssmDouble cellsum(IssmDouble* cell, int m); +void celldelete(IssmDouble** pcell, int m, int* indices, int nind); +void cellsplit(IssmDouble** pcell, int m, int i,IssmDouble scale) ; #endif //ifndef _MATRIXUTILS_H_ Index: /issm/trunk-jpl/src/m/classes/SMBgemb.m =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBgemb.m (revision 19553) +++ /issm/trunk-jpl/src/m/classes/SMBgemb.m (revision 19554) @@ -12,4 +12,14 @@ %from an automatic weather stations (AWS). Each property is therefore a matrix, of size (numberofvertices x number %of time steps. ) + + %solution choices + isgraingrowth; + isalbedo; + isshortwave; + isthermal; + isaccumulation; + ismelt; + isdensification; + isturbulentflux; %inputs: @@ -28,5 +38,4 @@ %settings: - spinUp = NaN; %number of cycles of met data run before output is calculated (default is 0) aIdx = NaN; %method for calculating albedo and subsurface absorption (default is 1) % 1: effective grain radius [Gardner & Sharp, 2009] @@ -97,6 +106,13 @@ function self = setdefaultparameters(self,mesh,geometry) % {{{ - - self.spinUp=0; + self.isgraingrowth=1; + self.isalbedo=1; + self.isshortwave=1; + self.isthermal=1; + self.isaccumulation=1; + self.ismelt=1; + self.isdensification=1; + self.isturbulentflux=1; + self.aIdx = 1; self.swIdx = 1; @@ -123,4 +139,14 @@ function md = checkconsistency(self,md,solution,analyses) % {{{ + + md = checkfield(md,'fieldname','smb.isgraingrowth','values',[0 1]); + md = checkfield(md,'fieldname','smb.isalbedo','values',[0 1]); + md = checkfield(md,'fieldname','smb.isshortwave','values',[0 1]); + md = checkfield(md,'fieldname','smb.isthermal','values',[0 1]); + md = checkfield(md,'fieldname','smb.isaccumulation','values',[0 1]); + md = checkfield(md,'fieldname','smb.ismelt','values',[0 1]); + md = checkfield(md,'fieldname','smb.isdensification','values',[0 1]); + md = checkfield(md,'fieldname','smb.isturbulentflux','values',[0 1]); + md = checkfield(md,'fieldname','smb.Ta','timeseries',1,'NaN',1,'>',273-60,'<',273+60); %60 celsius max value md = checkfield(md,'fieldname','smb.V','timeseries',1,'NaN',1,'>=',0,'<',45); %max 500 km/h @@ -130,10 +156,9 @@ md = checkfield(md,'fieldname','smb.eAir','timeseries',1,'NaN',1); - md = checkfield(md,'fieldname','smb.Tmean','NaN',1,'>',273-60,'<',273+60); %60 celsius max value - md = checkfield(md,'fieldname','smb.C','NaN',1,'>=',0); - md = checkfield(md,'fieldname','smb.Tz','NaN',1,'>=',0,'<=',5000); - md = checkfield(md,'fieldname','smb.Vz','NaN',1,'>=',0,'<=',5000); - - md = checkfield(md,'fieldname','smb.spinUp','NaN',1,'>=',0); + md = checkfield(md,'fieldname','smb.Tmean','size',[md.mesh.numberofelements 1],'NaN',1,'>',273-60,'<',273+60); %60 celsius max value + md = checkfield(md,'fieldname','smb.C','size',[md.mesh.numberofelements 1],'NaN',1,'>=',0); + md = checkfield(md,'fieldname','smb.Tz','size',[md.mesh.numberofelements 1],'NaN',1,'>=',0,'<=',5000); + md = checkfield(md,'fieldname','smb.Vz','size',[md.mesh.numberofelements 1],'NaN',1,'>=',0,'<=',5000); + md = checkfield(md,'fieldname','smb.aIdx','NaN',1,'values',[1,2,3,4]); md = checkfield(md,'fieldname','smb.swIdx','NaN',1,'values',[0,1]); @@ -170,4 +195,12 @@ disp(sprintf(' surface forcings for SMB GEMB model :')); + fielddisplay(self,'isgraingrowth','run grain growth module (default true)'); + fielddisplay(self,'isalbedo','run albedo module (default true)'); + fielddisplay(self,'isshortwave','run short wave module (default true)'); + fielddisplay(self,'isthermal','run thermal module (default true)'); + fielddisplay(self,'isaccumulation','run accumulation module (default true)'); + fielddisplay(self,'ismelt','run melting module (default true)'); + fielddisplay(self,'isdensification','run densification module (default true)'); + fielddisplay(self,'isturbulentflux','run turbulant heat fluxes module (default true)'); fielddisplay(self,'Ta','2 m air temperature, in Kelvin'); fielddisplay(self,'V','wind speed (m/s-1)'); @@ -188,5 +221,4 @@ fielddisplay(self,'zY','strech grid cells bellow top_z by a [top_dz * y ^ (cells bellow top_z)]'); fielddisplay(self,'outputFreq','output frequency in days (default is monthly, 30)'); - fielddisplay(self,'spinUp','number of cycles of met data run before output is calcualted (default is 0)'); fielddisplay(self,'aIdx',{'method for calculating albedo and subsurface absorption (default is 1)',... '1: effective grain radius [Gardner & Sharp, 2009]',... @@ -198,5 +230,5 @@ case {1 2} fielddisplay(self,'aSnow','new snow albedo (0.64 - 0.89)'); - fielddisplay(self,'aIce','range 0.27-0.58 for old snow'); + fielddisplay(self,'aIce','albedo of ice (0.27-0.58)'); case 3 fielddisplay(self,'cldFrac','average cloud amount'); @@ -224,27 +256,38 @@ WriteData(fid,'enum',SmbEnum(),'data',SMBgembEnum(),'format','Integer'); - WriteData(fid,'object',self,'class','smb','fieldname','Ta','format','DoubleMat','mattype',1,'scale',1,'timeserieslength',md.mesh.numberofelements+1); - WriteData(fid,'object',self,'class','smb','fieldname','V','format','DoubleMat','mattype',1,'scale',1,'timeserieslength',md.mesh.numberofelements+1); - WriteData(fid,'object',self,'class','smb','fieldname','dswrf','format','DoubleMat','mattype',1,'scale',1,'timeserieslength',md.mesh.numberofelements+1); - WriteData(fid,'object',self,'class','smb','fieldname','dlwrf','format','DoubleMat','mattype',1,'scale',1,'timeserieslength',md.mesh.numberofelements+1); - WriteData(fid,'object',self,'class','smb','fieldname','P','format','DoubleMat','mattype',1,'scale',1,'timeserieslength',md.mesh.numberofelements+1); - WriteData(fid,'object',self,'class','smb','fieldname','eAir','format','DoubleMat','mattype',1,'scale',1,'timeserieslength',md.mesh.numberofelements+1); - WriteData(fid,'object',self,'class','smb','fieldname','pAir','format','DoubleMat','mattype',1,'scale',1,'timeserieslength',md.mesh.numberofelements+1); - - WriteData(fid,'object',self,'class','smb','fieldname','Tmean','format','Double','scale',1); - WriteData(fid,'object',self,'class','smb','fieldname','C','format','Double','scale',1); - WriteData(fid,'object',self,'class','smb','fieldname','Tz','format','Double','scale',1); - WriteData(fid,'object',self,'class','smb','fieldname','Vz','format','Double','scale',1); - WriteData(fid,'object',self,'class','smb','fieldname','spinUp','format','Integer','scale',1); + WriteData(fid,'object',self,'class','smb','fieldname','isgraingrowth','format','Boolean'); + WriteData(fid,'object',self,'class','smb','fieldname','isalbedo','format','Boolean'); + WriteData(fid,'object',self,'class','smb','fieldname','isshortwave','format','Boolean'); + WriteData(fid,'object',self,'class','smb','fieldname','isthermal','format','Boolean'); + WriteData(fid,'object',self,'class','smb','fieldname','isaccumulation','format','Boolean'); + WriteData(fid,'object',self,'class','smb','fieldname','ismelt','format','Boolean'); + WriteData(fid,'object',self,'class','smb','fieldname','isdensification','format','Boolean'); + WriteData(fid,'object',self,'class','smb','fieldname','isturbulentflux','format','Boolean'); + WriteData(fid,'object',self,'class','smb','fieldname','isgraingrowth','format','Boolean'); + WriteData(fid,'object',self,'class','smb','fieldname','isgraingrowth','format','Boolean'); + + WriteData(fid,'object',self,'class','smb','fieldname','Ta','format','DoubleMat','mattype',2,'scale',1,'timeserieslength',md.mesh.numberofelements+1); + WriteData(fid,'object',self,'class','smb','fieldname','V','format','DoubleMat','mattype',2,'scale',1,'timeserieslength',md.mesh.numberofelements+1); + WriteData(fid,'object',self,'class','smb','fieldname','dswrf','format','DoubleMat','mattype',2,'scale',1,'timeserieslength',md.mesh.numberofelements+1); + WriteData(fid,'object',self,'class','smb','fieldname','dlwrf','format','DoubleMat','mattype',2,'scale',1,'timeserieslength',md.mesh.numberofelements+1); + WriteData(fid,'object',self,'class','smb','fieldname','P','format','DoubleMat','mattype',2,'scale',1,'timeserieslength',md.mesh.numberofelements+1); + WriteData(fid,'object',self,'class','smb','fieldname','eAir','format','DoubleMat','mattype',2,'scale',1,'timeserieslength',md.mesh.numberofelements+1); + WriteData(fid,'object',self,'class','smb','fieldname','pAir','format','DoubleMat','mattype',2,'scale',1,'timeserieslength',md.mesh.numberofelements+1); + + WriteData(fid,'object',self,'class','smb','fieldname','Tmean','format','DoubleMat','mattype',2,'scale',1); + WriteData(fid,'object',self,'class','smb','fieldname','C','format','DoubleMat','mattype',2,'scale',1); + WriteData(fid,'object',self,'class','smb','fieldname','Tz','format','DoubleMat','mattype',2,'scale',1); + WriteData(fid,'object',self,'class','smb','fieldname','Vz','format','DoubleMat','mattype',2,'scale',1); + WriteData(fid,'object',self,'class','smb','fieldname','zTop','format','DoubleMat','mattype',2,'scale',1); + WriteData(fid,'object',self,'class','smb','fieldname','dzTop','format','DoubleMat','mattype',2,'scale',1); + WriteData(fid,'object',self,'class','smb','fieldname','dzMin','format','DoubleMat','mattype',2,'scale',1); + WriteData(fid,'object',self,'class','smb','fieldname','zY','format','DoubleMat','mattype',2,'scale',1); + WriteData(fid,'object',self,'class','smb','fieldname','zMax','format','DoubleMat','mattype',2,'scale',1); + WriteData(fid,'object',self,'class','smb','fieldname','zMin','format','DoubleMat','mattype',2,'scale',1); + WriteData(fid,'object',self,'class','smb','fieldname','aIdx','format','Integer','scale',1); WriteData(fid,'object',self,'class','smb','fieldname','swIdx','format','Integer','scale',1); WriteData(fid,'object',self,'class','smb','fieldname','denIdx','format','Integer','scale',1); - WriteData(fid,'object',self,'class','smb','fieldname','zTop','format','DoubleMat','mattype',1,'scale',1); - WriteData(fid,'object',self,'class','smb','fieldname','dzTop','format','DoubleMat','mattype',1,'scale',1); - WriteData(fid,'object',self,'class','smb','fieldname','dzMin','format','DoubleMat','mattype',1,'scale',1); - WriteData(fid,'object',self,'class','smb','fieldname','zY','format','DoubleMat','mattype',1,'scale',1); - WriteData(fid,'object',self,'class','smb','fieldname','zMax','format','DoubleMat','mattype',1,'scale',1); - WriteData(fid,'object',self,'class','smb','fieldname','zMin','format','DoubleMat','mattype',1,'scale',1); WriteData(fid,'object',self,'class','smb','fieldname','outputFreq','format','Double','scale',1); @@ -255,4 +298,10 @@ WriteData(fid,'object',self,'class','smb','fieldname','t0dry','format','Double','scale',1); WriteData(fid,'object',self,'class','smb','fieldname','K','format','Double','scale',1); + + %figure out dt from forcings: + time=self.Ta(end,:); %assume all forcings are one the same time step + dtime=diff(time,1); + dt=min(dtime); + WriteData(fid,'data',dt,'enum',SmbDtEnum,'format','Double','scale',yts); %process requested outputs @@ -264,5 +313,4 @@ end WriteData(fid,'data',outputs,'enum',SmbRequestedOutputsEnum,'format','StringArray'); - end % }}} end Index: /issm/trunk-jpl/src/m/enum/DoubleArrayInputEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/DoubleArrayInputEnum.m (revision 19554) +++ /issm/trunk-jpl/src/m/enum/DoubleArrayInputEnum.m (revision 19554) @@ -0,0 +1,11 @@ +function macro=DoubleArrayInputEnum() +%DOUBLEARRAYINPUTENUM - Enum of DoubleArrayInput +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=DoubleArrayInputEnum() + +macro=StringToEnum('DoubleArrayInput'); Index: /issm/trunk-jpl/src/m/enum/EnumDefinitions.py =================================================================== --- /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 19553) +++ /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 19554) @@ -350,4 +350,5 @@ def SmbNumRequestedOutputsEnum(): return StringToEnum("SmbNumRequestedOutputs")[0] def SmbRequestedOutputsEnum(): return StringToEnum("SmbRequestedOutputs")[0] +def SmbIsInitializedEnum(): return StringToEnum("SmbIsInitialized")[0] def SMBforcingEnum(): return StringToEnum("SMBforcing")[0] def SmbMassBalanceEnum(): return StringToEnum("SmbMassBalance")[0] @@ -364,5 +365,6 @@ def SmbTzEnum(): return StringToEnum("SmbTz")[0] def SmbVzEnum(): return StringToEnum("SmbVz")[0] -def SmbSpinUpEnum(): return StringToEnum("SmbSpinUp")[0] +def SmbDtEnum(): return StringToEnum("SmbDt")[0] +def SmbDzEnum(): return StringToEnum("SmbDz")[0] def SmbAIdxEnum(): return StringToEnum("SmbAIdx")[0] def SmbSwIdxEnum(): return StringToEnum("SmbSwIdx")[0] @@ -381,4 +383,20 @@ def SmbT0dryEnum(): return StringToEnum("SmbT0dry")[0] def SmbKEnum(): return StringToEnum("SmbK")[0] +def SmbDEnum(): return StringToEnum("SmbD")[0] +def SmbReEnum(): return StringToEnum("SmbRe")[0] +def SmbGdnEnum(): return StringToEnum("SmbGdn")[0] +def SmbGspEnum(): return StringToEnum("SmbGsp")[0] +def SmbECEnum(): return StringToEnum("SmbEC")[0] +def SmbWEnum(): return StringToEnum("SmbW")[0] +def SmbAEnum(): return StringToEnum("SmbA")[0] +def SmbTEnum(): return StringToEnum("SmbT")[0] +def SmbIsgraingrowthEnum(): return StringToEnum("SmbIsgraingrowth")[0] +def SmbIsalbedoEnum(): return StringToEnum("SmbIsalbedo")[0] +def SmbIsshortwaveEnum(): return StringToEnum("SmbIsshortwave")[0] +def SmbIsthermalEnum(): return StringToEnum("SmbIsthermal")[0] +def SmbIsaccumulationEnum(): return StringToEnum("SmbIsaccumulation")[0] +def SmbIsmeltEnum(): return StringToEnum("SmbIsmelt")[0] +def SmbIsdensificationEnum(): return StringToEnum("SmbIsdensification")[0] +def SmbIsturbulentfluxEnum(): return StringToEnum("SmbIsturbulentflux")[0] def SMBpddEnum(): return StringToEnum("SMBpdd")[0] def SmbDelta18oEnum(): return StringToEnum("SmbDelta18o")[0] @@ -504,4 +522,5 @@ def DatasetInputEnum(): return StringToEnum("DatasetInput")[0] def DoubleInputEnum(): return StringToEnum("DoubleInput")[0] +def DoubleArrayInputEnum(): return StringToEnum("DoubleArrayInput")[0] def DataSetParamEnum(): return StringToEnum("DataSetParam")[0] def DoubleMatArrayParamEnum(): return StringToEnum("DoubleMatArrayParam")[0] @@ -672,4 +691,5 @@ def GiaWEnum(): return StringToEnum("GiaW")[0] def P0Enum(): return StringToEnum("P0")[0] +def P0ArrayEnum(): return StringToEnum("P0Array")[0] def P1Enum(): return StringToEnum("P1")[0] def P1DGEnum(): return StringToEnum("P1DG")[0] Index: /issm/trunk-jpl/src/m/enum/P0ArrayEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/P0ArrayEnum.m (revision 19554) +++ /issm/trunk-jpl/src/m/enum/P0ArrayEnum.m (revision 19554) @@ -0,0 +1,11 @@ +function macro=P0ArrayEnum() +%P0ARRAYENUM - Enum of P0Array +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=P0ArrayEnum() + +macro=StringToEnum('P0Array'); Index: /issm/trunk-jpl/src/m/enum/SmbAEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/SmbAEnum.m (revision 19554) +++ /issm/trunk-jpl/src/m/enum/SmbAEnum.m (revision 19554) @@ -0,0 +1,11 @@ +function macro=SmbAEnum() +%SMBAENUM - Enum of SmbA +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=SmbAEnum() + +macro=StringToEnum('SmbA'); Index: /issm/trunk-jpl/src/m/enum/SmbDEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/SmbDEnum.m (revision 19554) +++ /issm/trunk-jpl/src/m/enum/SmbDEnum.m (revision 19554) @@ -0,0 +1,11 @@ +function macro=SmbDEnum() +%SMBDENUM - Enum of SmbD +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=SmbDEnum() + +macro=StringToEnum('SmbD'); Index: /issm/trunk-jpl/src/m/enum/SmbDtEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/SmbDtEnum.m (revision 19554) +++ /issm/trunk-jpl/src/m/enum/SmbDtEnum.m (revision 19554) @@ -0,0 +1,11 @@ +function macro=SmbDtEnum() +%SMBDTENUM - Enum of SmbDt +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=SmbDtEnum() + +macro=StringToEnum('SmbDt'); Index: /issm/trunk-jpl/src/m/enum/SmbDzEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/SmbDzEnum.m (revision 19554) +++ /issm/trunk-jpl/src/m/enum/SmbDzEnum.m (revision 19554) @@ -0,0 +1,11 @@ +function macro=SmbDzEnum() +%SMBDZENUM - Enum of SmbDz +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=SmbDzEnum() + +macro=StringToEnum('SmbDz'); Index: /issm/trunk-jpl/src/m/enum/SmbECEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/SmbECEnum.m (revision 19554) +++ /issm/trunk-jpl/src/m/enum/SmbECEnum.m (revision 19554) @@ -0,0 +1,11 @@ +function macro=SmbECEnum() +%SMBECENUM - Enum of SmbEC +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=SmbECEnum() + +macro=StringToEnum('SmbEC'); Index: /issm/trunk-jpl/src/m/enum/SmbGdnEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/SmbGdnEnum.m (revision 19554) +++ /issm/trunk-jpl/src/m/enum/SmbGdnEnum.m (revision 19554) @@ -0,0 +1,11 @@ +function macro=SmbGdnEnum() +%SMBGDNENUM - Enum of SmbGdn +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=SmbGdnEnum() + +macro=StringToEnum('SmbGdn'); Index: /issm/trunk-jpl/src/m/enum/SmbGspEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/SmbGspEnum.m (revision 19554) +++ /issm/trunk-jpl/src/m/enum/SmbGspEnum.m (revision 19554) @@ -0,0 +1,11 @@ +function macro=SmbGspEnum() +%SMBGSPENUM - Enum of SmbGsp +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=SmbGspEnum() + +macro=StringToEnum('SmbGsp'); Index: /issm/trunk-jpl/src/m/enum/SmbIsInitializedEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/SmbIsInitializedEnum.m (revision 19554) +++ /issm/trunk-jpl/src/m/enum/SmbIsInitializedEnum.m (revision 19554) @@ -0,0 +1,11 @@ +function macro=SmbIsInitializedEnum() +%SMBISINITIALIZEDENUM - Enum of SmbIsInitialized +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=SmbIsInitializedEnum() + +macro=StringToEnum('SmbIsInitialized'); Index: /issm/trunk-jpl/src/m/enum/SmbIsaccumulationEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/SmbIsaccumulationEnum.m (revision 19554) +++ /issm/trunk-jpl/src/m/enum/SmbIsaccumulationEnum.m (revision 19554) @@ -0,0 +1,11 @@ +function macro=SmbIsaccumulationEnum() +%SMBISACCUMULATIONENUM - Enum of SmbIsaccumulation +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=SmbIsaccumulationEnum() + +macro=StringToEnum('SmbIsaccumulation'); Index: /issm/trunk-jpl/src/m/enum/SmbIsalbedoEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/SmbIsalbedoEnum.m (revision 19554) +++ /issm/trunk-jpl/src/m/enum/SmbIsalbedoEnum.m (revision 19554) @@ -0,0 +1,11 @@ +function macro=SmbIsalbedoEnum() +%SMBISALBEDOENUM - Enum of SmbIsalbedo +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=SmbIsalbedoEnum() + +macro=StringToEnum('SmbIsalbedo'); Index: /issm/trunk-jpl/src/m/enum/SmbIsdensificationEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/SmbIsdensificationEnum.m (revision 19554) +++ /issm/trunk-jpl/src/m/enum/SmbIsdensificationEnum.m (revision 19554) @@ -0,0 +1,11 @@ +function macro=SmbIsdensificationEnum() +%SMBISDENSIFICATIONENUM - Enum of SmbIsdensification +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=SmbIsdensificationEnum() + +macro=StringToEnum('SmbIsdensification'); Index: /issm/trunk-jpl/src/m/enum/SmbIsgraingrowthEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/SmbIsgraingrowthEnum.m (revision 19554) +++ /issm/trunk-jpl/src/m/enum/SmbIsgraingrowthEnum.m (revision 19554) @@ -0,0 +1,11 @@ +function macro=SmbIsgraingrowthEnum() +%SMBISGRAINGROWTHENUM - Enum of SmbIsgraingrowth +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=SmbIsgraingrowthEnum() + +macro=StringToEnum('SmbIsgraingrowth'); Index: /issm/trunk-jpl/src/m/enum/SmbIsmeltEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/SmbIsmeltEnum.m (revision 19554) +++ /issm/trunk-jpl/src/m/enum/SmbIsmeltEnum.m (revision 19554) @@ -0,0 +1,11 @@ +function macro=SmbIsmeltEnum() +%SMBISMELTENUM - Enum of SmbIsmelt +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=SmbIsmeltEnum() + +macro=StringToEnum('SmbIsmelt'); Index: /issm/trunk-jpl/src/m/enum/SmbIsshortwaveEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/SmbIsshortwaveEnum.m (revision 19554) +++ /issm/trunk-jpl/src/m/enum/SmbIsshortwaveEnum.m (revision 19554) @@ -0,0 +1,11 @@ +function macro=SmbIsshortwaveEnum() +%SMBISSHORTWAVEENUM - Enum of SmbIsshortwave +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=SmbIsshortwaveEnum() + +macro=StringToEnum('SmbIsshortwave'); Index: /issm/trunk-jpl/src/m/enum/SmbIsthermalEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/SmbIsthermalEnum.m (revision 19554) +++ /issm/trunk-jpl/src/m/enum/SmbIsthermalEnum.m (revision 19554) @@ -0,0 +1,11 @@ +function macro=SmbIsthermalEnum() +%SMBISTHERMALENUM - Enum of SmbIsthermal +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=SmbIsthermalEnum() + +macro=StringToEnum('SmbIsthermal'); Index: /issm/trunk-jpl/src/m/enum/SmbIsturbulentfluxEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/SmbIsturbulentfluxEnum.m (revision 19554) +++ /issm/trunk-jpl/src/m/enum/SmbIsturbulentfluxEnum.m (revision 19554) @@ -0,0 +1,11 @@ +function macro=SmbIsturbulentfluxEnum() +%SMBISTURBULENTFLUXENUM - Enum of SmbIsturbulentflux +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=SmbIsturbulentfluxEnum() + +macro=StringToEnum('SmbIsturbulentflux'); Index: /issm/trunk-jpl/src/m/enum/SmbReEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/SmbReEnum.m (revision 19554) +++ /issm/trunk-jpl/src/m/enum/SmbReEnum.m (revision 19554) @@ -0,0 +1,11 @@ +function macro=SmbReEnum() +%SMBREENUM - Enum of SmbRe +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=SmbReEnum() + +macro=StringToEnum('SmbRe'); Index: /issm/trunk-jpl/src/m/enum/SmbTEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/SmbTEnum.m (revision 19554) +++ /issm/trunk-jpl/src/m/enum/SmbTEnum.m (revision 19554) @@ -0,0 +1,11 @@ +function macro=SmbTEnum() +%SMBTENUM - Enum of SmbT +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=SmbTEnum() + +macro=StringToEnum('SmbT'); Index: /issm/trunk-jpl/src/m/enum/SmbWEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/SmbWEnum.m (revision 19554) +++ /issm/trunk-jpl/src/m/enum/SmbWEnum.m (revision 19554) @@ -0,0 +1,11 @@ +function macro=SmbWEnum() +%SMBWENUM - Enum of SmbW +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=SmbWEnum() + +macro=StringToEnum('SmbW'); Index: /issm/trunk-jpl/test/NightlyRun/gemb.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/gemb.m (revision 19553) +++ /issm/trunk-jpl/test/NightlyRun/gemb.m (revision 19554) @@ -23,11 +23,21 @@ md.smb.pAir=[repmat(inputs.pAir0',md.mesh.numberofelements,1);dateN']; -md.smb.Vz=inputs.LP.Vz; -md.smb.Tz=inputs.LP.Tz; -md.smb.Tmean=inputs.LP.Tmean; -md.smb.C=inputs.LP.C; +md.smb.Vz=repmat(inputs.LP.Vz,md.mesh.numberofelements,1); +md.smb.Tz=repmat(inputs.LP.Tz,md.mesh.numberofelements,1); +md.smb.Tmean=repmat(inputs.LP.Tmean,md.mesh.numberofelements,1); +md.smb.C=repmat(inputs.LP.C,md.mesh.numberofelements,1); %settings -md.smb.spinUp=2; +md.smb.ismelt=0; +md.smb.isthermal=0; +md.smb.isalbedo=0; +md.transient.isstressbalance=0; +md.transient.ismasstransport=0; +md.transient.isthermal=0; +md.timestepping.start_time=2000; +md.timestepping.final_time=2001; +md.timestepping.time_step=1/12; %monthly + +md.verbose=verbose('11111111'); %Run transient