Ice Sheet System Model  4.18
Code documentation
Functions
FloatingiceMeltingRatePicox.h File Reference

header file for Floatingice melting rate More...

#include "../../classes/classes.h"

Go to the source code of this file.

Functions

void FloatingiceMeltingRatePicox (FemModel *femmodel)
 
void UpdateBoxIdsPico (FemModel *femmodel)
 
void ComputeBoxAreasPico (FemModel *femmodel)
 
void UpdateBoxPico (FemModel *femmodel, int loopboxid)
 
void ComputeAverageOceanvarsPico (FemModel *femmodel, int boxid)
 
void ComputeBasalMeltPlume (FemModel *femmodel)
 

Detailed Description

header file for Floatingice melting rate

Definition in file FloatingiceMeltingRatePicox.h.

Function Documentation

◆ FloatingiceMeltingRatePicox()

void FloatingiceMeltingRatePicox ( FemModel *  femmodel)

Definition at line 10 of file FloatingiceMeltingRatePicox.cpp.

10  {/*{{{*/
11 
12  int maxbox;
13  bool isplume;
14 
15  /*First, reset all melt to 0 */
16  for(int i=0;i<femmodel->elements->Size();i++){
17  Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
18  int numvertices = element->GetNumberOfVertices();
19  IssmDouble* values = xNewZeroInit<IssmDouble>(numvertices);
20  element->AddInput2(BasalforcingsFloatingiceMeltingRateEnum,values,P1Enum);
21  xDelete<IssmDouble>(values);
22  }
23 
24  /*PICO melt rate parameterization (Reese et al., 2018)*/
25  femmodel->parameters->FindParam(&maxbox,BasalforcingsPicoMaxboxcountEnum);
26  UpdateBoxIdsPico(femmodel);
27  ComputeBoxAreasPico(femmodel);
28  for(int i=0;i<maxbox;i++){
29  UpdateBoxPico(femmodel,i);
30  ComputeAverageOceanvarsPico(femmodel,i);
31  }
32 
33  /*Optional buoyant plume melt rate parameterization (Lazeroms et al., 2018) */
34  femmodel->parameters->FindParam(&isplume,BasalforcingsPicoIsplumeEnum);
35  if(isplume) ComputeBasalMeltPlume(femmodel);
36 }/*}}}*/

◆ UpdateBoxIdsPico()

void UpdateBoxIdsPico ( FemModel *  femmodel)

Definition at line 38 of file FloatingiceMeltingRatePicox.cpp.

38  {/*{{{*/
39 
40  int numvertices,num_basins,maxbox,basinid;
41  IssmDouble dist_max;
42  IssmDouble* distances=NULL;
43 
44  femmodel->parameters->FindParam(&num_basins,BasalforcingsPicoNumBasinsEnum);
45  femmodel->parameters->FindParam(&maxbox,BasalforcingsPicoMaxboxcountEnum);
46  IssmDouble* dmax_basin_cpu=xNew<IssmDouble>(num_basins);
47 
48  InputDuplicatex(femmodel,MaskOceanLevelsetEnum,DistanceToGroundinglineEnum);
49  femmodel->DistanceToFieldValue(MaskOceanLevelsetEnum,0.,DistanceToGroundinglineEnum);
50 
51  InputDuplicatex(femmodel,MaskIceLevelsetEnum,DistanceToCalvingfrontEnum);
52  femmodel->DistanceToFieldValue(MaskIceLevelsetEnum,0.,DistanceToCalvingfrontEnum);
53 
54  /*find maximum distance to grounding line per domain and per basin*/
55  IssmDouble maxdist_cpu=-1.;
56  for(int i=0;i<num_basins;i++){dmax_basin_cpu[i]=-1;}
57  for(int i=0;i<femmodel->elements->Size();i++){
58  Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
59  if(!element->IsIceInElement() || !element->IsFloating()) continue;
60  numvertices = element->GetNumberOfVertices();
61  distances=xNew<IssmDouble>(numvertices);
62  element->GetInputListOnVertices(&distances[0],DistanceToGroundinglineEnum);
63  element->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum);
64  for(int k=0; k<numvertices; k++){
65  if(fabs(distances[k])>maxdist_cpu){maxdist_cpu=fabs(distances[k]);}
66  if(fabs(distances[k])>dmax_basin_cpu[basinid]){dmax_basin_cpu[basinid]=fabs(distances[k]);}
67  }
68  xDelete<IssmDouble>(distances);
69  }
70 
71  /*Synchronize across cpus*/
72  IssmDouble* dmax_basin=xNew<IssmDouble>(num_basins);
73  ISSM_MPI_Allreduce((void*)&maxdist_cpu,(void*)&dist_max,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,IssmComm::GetComm());
74  ISSM_MPI_Allreduce(dmax_basin_cpu,dmax_basin,num_basins,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,IssmComm::GetComm());
75 
76  /*Define maximum number of boxes per basin*/
77  int* nd=xNew<int>(num_basins);
78  for(int i=0; i<num_basins;i++){
79  IssmDouble val=sqrt(dmax_basin[i]/dist_max)*(maxbox-1);
80 
81  #ifdef _HAVE_AD_
82  _error_("Check the implementation of floor below");
83  /*Do not use floor when AD is on*/
84  int k=0; while(k<val+.5){k++;}
85  nd[i]=k;
86 
87  #else
88  nd[i]= reCast<int>(floor(val));
89  #endif
90  }
91 
92  /*Assign box numbers*/
93  for(int i=0;i<femmodel->elements->Size();i++){
94  Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
95  element->PicoUpdateBoxid(nd);
96  }
97 
98  /*Cleanup and return */
99  xDelete<int>(nd);
100  xDelete<IssmDouble>(dmax_basin);
101  xDelete<IssmDouble>(dmax_basin_cpu);
102 
103 }/*}}}*/

◆ ComputeBoxAreasPico()

void ComputeBoxAreasPico ( FemModel *  femmodel)

Definition at line 104 of file FloatingiceMeltingRatePicox.cpp.

104  {/*{{{*/
105 
106  int num_basins,maxbox,basinid,boxid,domaintype;
107  IssmDouble dist_max,area;
108 
109  femmodel->parameters->FindParam(&num_basins,BasalforcingsPicoNumBasinsEnum);
110  femmodel->parameters->FindParam(&maxbox,BasalforcingsPicoMaxboxcountEnum);
111 
112  IssmDouble* boxareas=xNewZeroInit<IssmDouble>(num_basins*maxbox);
113 
114  for(int i=0;i<femmodel->elements->Size();i++){
115  Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
116  if(!element->IsOnBase()) continue;
117  Element* basalelement = element->SpawnBasalElement();
118  if(!basalelement->IsIceInElement() || !basalelement->IsFloating()) continue;
119  basalelement->GetInputValue(&boxid,BasalforcingsPicoBoxIdEnum);
120  basalelement->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum);
121  boxareas[basinid*maxbox+boxid]+=basalelement->GetHorizontalSurfaceArea();
122  basalelement->FindParam(&domaintype,DomainTypeEnum);
123  if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
124  }
125 
126  /*Synchronize across cpus*/
127  IssmDouble* sumareas =xNew<IssmDouble>(num_basins*maxbox);
128  ISSM_MPI_Allreduce(boxareas,sumareas,num_basins*maxbox,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
129  //if(sumareas[0]==0){_error_("No elements in box 0, basal meltrates will be 0. Consider decreasing md.basalforcings.maxboxcount or refining your mesh!");}
130 
131  /*Update parameters to keep track of the new areas in future calculations*/
132  femmodel->parameters->AddObject(new DoubleVecParam(BasalforcingsPicoBoxAreaEnum,sumareas,maxbox*num_basins));
133 
134  /*Cleanup and return */
135  xDelete<IssmDouble>(boxareas);
136  xDelete<IssmDouble>(sumareas);
137 
138 }/*}}}*/

◆ UpdateBoxPico()

void UpdateBoxPico ( FemModel *  femmodel,
int  loopboxid 
)

Definition at line 139 of file FloatingiceMeltingRatePicox.cpp.

139  {/*{{{*/
140  for(int i=0;i<femmodel->elements->Size();i++){
141  Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
142  element->PicoUpdateBox(loopboxid);
143  }
144 }/*}}}*/

◆ ComputeAverageOceanvarsPico()

void ComputeAverageOceanvarsPico ( FemModel *  femmodel,
int  boxid 
)

Definition at line 145 of file FloatingiceMeltingRatePicox.cpp.

145  {/*{{{*/
146 
147  int num_basins, basinid, maxbox, M, domaintype;
148  IssmDouble area, toc, soc, overturning;
149  IssmDouble* boxareas=NULL;
150  IssmDouble* overturning_weighted_avg=NULL;
151 
152  femmodel->parameters->FindParam(&num_basins,BasalforcingsPicoNumBasinsEnum);
153  femmodel->parameters->FindParam(&maxbox,BasalforcingsPicoMaxboxcountEnum);
154  femmodel->parameters->FindParam(&boxareas,&M, BasalforcingsPicoBoxAreaEnum);
155  IssmDouble* toc_weighted_avg = xNewZeroInit<IssmDouble>(num_basins);
156  IssmDouble* soc_weighted_avg = xNewZeroInit<IssmDouble>(num_basins);
157  IssmDouble* toc_sumweightedavg = xNewZeroInit<IssmDouble>(num_basins);
158  IssmDouble* soc_sumweightedavg = xNewZeroInit<IssmDouble>(num_basins);
159  IssmDouble* overturning_sumweightedavg = xNewZeroInit<IssmDouble>(num_basins);
160 
161  /* Compute Toc and Soc weighted avg (boxes 0 to n-1) */
162  for(int i=0;i<femmodel->elements->Size();i++){
163  Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
164  if(!element->IsOnBase()) continue;
165  Element* basalelement = element->SpawnBasalElement();
166  if(!basalelement->IsIceInElement() || !basalelement->IsFloating()) continue;
167  int el_boxid;
168  basalelement->GetInputValue(&el_boxid,BasalforcingsPicoBoxIdEnum);
169  if(el_boxid!=boxid) continue;
170 
171  Input2* tocs_input=basalelement->GetInput2(BasalforcingsPicoSubShelfOceanTempEnum); _assert_(tocs_input);
172  Input2* socs_input=basalelement->GetInput2(BasalforcingsPicoSubShelfOceanSalinityEnum); _assert_(socs_input);
173 
174  basalelement->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum);
175  Gauss* gauss=basalelement->NewGauss(1); gauss->GaussPoint(0);
176  tocs_input->GetInputValue(&toc,gauss);
177  socs_input->GetInputValue(&soc,gauss);
178  delete gauss;
179  area=basalelement->GetHorizontalSurfaceArea();
180  toc_weighted_avg[basinid]+=toc*area;
181  soc_weighted_avg[basinid]+=soc*area;
182  basalelement->FindParam(&domaintype,DomainTypeEnum);
183  if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
184  }
185 
186  /*Syncronize across cpus*/
187  ISSM_MPI_Allreduce(toc_weighted_avg,toc_sumweightedavg,num_basins,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
188  ISSM_MPI_Allreduce(soc_weighted_avg,soc_sumweightedavg,num_basins,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
189 
190  for(int k=0;k<num_basins;k++){
191  int p=k*maxbox+boxid;
192  if(boxareas[p]==0) continue;
193  toc_sumweightedavg[k] = toc_sumweightedavg[k]/boxareas[p];
194  soc_sumweightedavg[k] = soc_sumweightedavg[k]/boxareas[p];
195  }
196 
197  femmodel->parameters->AddObject(new DoubleVecParam(BasalforcingsPicoAverageTemperatureEnum,toc_sumweightedavg,num_basins));
198  femmodel->parameters->AddObject(new DoubleVecParam(BasalforcingsPicoAverageSalinityEnum,soc_sumweightedavg,num_basins));
199 
200  /* Compute overturning weighted avg (box 0 only) */
201  if(boxid==0){
202  overturning_weighted_avg=xNewZeroInit<IssmDouble>(num_basins);
203  for(int i=0;i<femmodel->elements->Size();i++){
204  Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
205  if(!element->IsOnBase()) continue;
206  Element* basalelement = element->SpawnBasalElement();
207  if(!basalelement->IsIceInElement() || !basalelement->IsFloating()) continue;
208  int el_boxid;
209  basalelement->GetInputValue(&el_boxid,BasalforcingsPicoBoxIdEnum);
210  if(el_boxid!=boxid) continue;
211 
212  Input2* overturnings_input=basalelement->GetInput2(BasalforcingsPicoSubShelfOceanOverturningEnum); _assert_(overturnings_input);
213 
214  basalelement->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum);
215  Gauss* gauss=basalelement->NewGauss(1); gauss->GaussPoint(0);
216  overturnings_input->GetInputValue(&overturning,gauss);
217  delete gauss;
218  area=basalelement->GetHorizontalSurfaceArea();
219  overturning_weighted_avg[basinid]+=overturning*area;
220  basalelement->FindParam(&domaintype,DomainTypeEnum);
221  if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
222  }
223 
224  /*Syncronize across cpus*/
225  ISSM_MPI_Allreduce(overturning_weighted_avg,overturning_sumweightedavg,num_basins,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
226 
227  for(int k=0;k<num_basins;k++){
228  int p=k*maxbox+boxid;
229  if(boxareas[p]==0.) continue;
230  overturning_sumweightedavg[k] = overturning_sumweightedavg[k]/boxareas[p];
231  }
232  femmodel->parameters->AddObject(new DoubleVecParam(BasalforcingsPicoAverageOverturningEnum,overturning_sumweightedavg,num_basins));
233  }
234 
235  /*Cleanup and return */
236  xDelete<IssmDouble>(overturning_sumweightedavg);
237  xDelete<IssmDouble>(toc_sumweightedavg);
238  xDelete<IssmDouble>(soc_sumweightedavg);
239  xDelete<IssmDouble>(overturning_weighted_avg);
240  xDelete<IssmDouble>(toc_weighted_avg);
241  xDelete<IssmDouble>(soc_weighted_avg);
242  xDelete<IssmDouble>(boxareas);
243 }/*}}}*/

◆ ComputeBasalMeltPlume()

void ComputeBasalMeltPlume ( FemModel *  femmodel)

Definition at line 244 of file FloatingiceMeltingRatePicox.cpp.

244  {/*{{{*/
245  for(int i=0;i<femmodel->elements->Size();i++){
246  Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
247  element->PicoComputeBasalMelt();
248  }
249 }/*}}}*/
DataSet::Size
int Size()
Definition: DataSet.cpp:399
BasalforcingsPicoBasinIdEnum
@ BasalforcingsPicoBasinIdEnum
Definition: EnumDefinitions.h:486
_assert_
#define _assert_(ignore)
Definition: exceptions.h:37
ComputeBasalMeltPlume
void ComputeBasalMeltPlume(FemModel *femmodel)
Definition: FloatingiceMeltingRatePicox.cpp:244
IssmDouble
double IssmDouble
Definition: types.h:37
Element::IsOnBase
bool IsOnBase()
Definition: Element.cpp:1984
ISSM_MPI_Allreduce
int ISSM_MPI_Allreduce(void *sendbuf, void *recvbuf, int count, ISSM_MPI_Datatype datatype, ISSM_MPI_Op op, ISSM_MPI_Comm comm)
Definition: issmmpi.cpp:116
Element::FindParam
void FindParam(bool *pvalue, int paramenum)
Definition: Element.cpp:933
BasalforcingsPicoAverageSalinityEnum
@ BasalforcingsPicoAverageSalinityEnum
Definition: EnumDefinitions.h:77
Element::PicoUpdateBoxid
void PicoUpdateBoxid(int *pmax_boxid_basin)
Definition: Element.cpp:2509
BasalforcingsPicoSubShelfOceanOverturningEnum
@ BasalforcingsPicoSubShelfOceanOverturningEnum
Definition: EnumDefinitions.h:489
MaskOceanLevelsetEnum
@ MaskOceanLevelsetEnum
Definition: EnumDefinitions.h:640
UpdateBoxPico
void UpdateBoxPico(FemModel *femmodel, int loopboxid)
Definition: FloatingiceMeltingRatePicox.cpp:139
MaskIceLevelsetEnum
@ MaskIceLevelsetEnum
Definition: EnumDefinitions.h:641
ISSM_MPI_SUM
#define ISSM_MPI_SUM
Definition: issmmpi.h:134
FemModel::parameters
Parameters * parameters
Definition: FemModel.h:46
IssmComm::GetComm
static ISSM_MPI_Comm GetComm(void)
Definition: IssmComm.cpp:30
Element::SpawnBasalElement
virtual Element * SpawnBasalElement(void)=0
ComputeBoxAreasPico
void ComputeBoxAreasPico(FemModel *femmodel)
Definition: FloatingiceMeltingRatePicox.cpp:104
Element::PicoUpdateBox
void PicoUpdateBox(int loopboxid)
Definition: Element.cpp:2550
ISSM_MPI_MAX
#define ISSM_MPI_MAX
Definition: issmmpi.h:131
Element::IsFloating
bool IsFloating()
Definition: Element.cpp:1987
UpdateBoxIdsPico
void UpdateBoxIdsPico(FemModel *femmodel)
Definition: FloatingiceMeltingRatePicox.cpp:38
Parameters::AddObject
void AddObject(Param *newparam)
Definition: Parameters.cpp:67
BasalforcingsPicoAverageTemperatureEnum
@ BasalforcingsPicoAverageTemperatureEnum
Definition: EnumDefinitions.h:78
Element::GetInput2
virtual Input2 * GetInput2(int inputenum)=0
Element::DeleteMaterials
void DeleteMaterials(void)
Definition: Element.cpp:429
Element::GetHorizontalSurfaceArea
virtual IssmDouble GetHorizontalSurfaceArea(void)
Definition: Element.h:245
Element
Definition: Element.h:41
BasalforcingsPicoNumBasinsEnum
@ BasalforcingsPicoNumBasinsEnum
Definition: EnumDefinitions.h:85
P1Enum
@ P1Enum
Definition: EnumDefinitions.h:662
Domain2DhorizontalEnum
@ Domain2DhorizontalEnum
Definition: EnumDefinitions.h:534
ISSM_MPI_DOUBLE
#define ISSM_MPI_DOUBLE
Definition: issmmpi.h:125
ComputeAverageOceanvarsPico
void ComputeAverageOceanvarsPico(FemModel *femmodel, int boxid)
Definition: FloatingiceMeltingRatePicox.cpp:145
DomainTypeEnum
@ DomainTypeEnum
Definition: EnumDefinitions.h:124
Element::AddInput2
virtual void AddInput2(int input_enum, IssmDouble *values, int interpolation_enum)
Definition: Element.h:216
Element::NewGauss
virtual Gauss * NewGauss(void)=0
BasalforcingsPicoSubShelfOceanSalinityEnum
@ BasalforcingsPicoSubShelfOceanSalinityEnum
Definition: EnumDefinitions.h:490
DistanceToCalvingfrontEnum
@ DistanceToCalvingfrontEnum
Definition: EnumDefinitions.h:532
BasalforcingsFloatingiceMeltingRateEnum
@ BasalforcingsFloatingiceMeltingRateEnum
Definition: EnumDefinitions.h:476
BasalforcingsPicoBoxAreaEnum
@ BasalforcingsPicoBoxAreaEnum
Definition: EnumDefinitions.h:79
BasalforcingsPicoBoxIdEnum
@ BasalforcingsPicoBoxIdEnum
Definition: EnumDefinitions.h:487
InputDuplicatex
void InputDuplicatex(FemModel *femmodel, int original_enum, int new_enum)
Definition: InputDuplicatex.cpp:10
FemModel::elements
Elements * elements
Definition: FemModel.h:44
Input2
Definition: Input2.h:18
BasalforcingsPicoMaxboxcountEnum
@ BasalforcingsPicoMaxboxcountEnum
Definition: EnumDefinitions.h:84
Element::IsIceInElement
bool IsIceInElement()
Definition: Element.cpp:2021
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
DataSet::GetObjectByOffset
Object * GetObjectByOffset(int offset)
Definition: DataSet.cpp:334
DoubleVecParam
Definition: DoubleVecParam.h:20
DistanceToGroundinglineEnum
@ DistanceToGroundinglineEnum
Definition: EnumDefinitions.h:533
Gauss::GaussPoint
virtual void GaussPoint(int ig)=0
Element::GetInputValue
void GetInputValue(bool *pvalue, int enum_type)
Definition: Element.cpp:1177
Element::PicoComputeBasalMelt
void PicoComputeBasalMelt()
Definition: Element.cpp:2684
Parameters::FindParam
void FindParam(bool *pinteger, int enum_type)
Definition: Parameters.cpp:262
Input2::GetInputValue
virtual void GetInputValue(IssmDouble *pvalue, Gauss *gauss)
Definition: Input2.h:38
Element::GetNumberOfVertices
virtual int GetNumberOfVertices(void)=0
Element::GetInputListOnVertices
void GetInputListOnVertices(IssmDouble *pvalue, int enumtype)
Definition: Element.cpp:1131
FemModel::DistanceToFieldValue
void DistanceToFieldValue(int fieldenum, IssmDouble fieldvalue, int distanceenum)
Definition: FemModel.cpp:1103
Gauss
Definition: Gauss.h:8
BasalforcingsPicoAverageOverturningEnum
@ BasalforcingsPicoAverageOverturningEnum
Definition: EnumDefinitions.h:76
BasalforcingsPicoIsplumeEnum
@ BasalforcingsPicoIsplumeEnum
Definition: EnumDefinitions.h:83
BasalforcingsPicoSubShelfOceanTempEnum
@ BasalforcingsPicoSubShelfOceanTempEnum
Definition: EnumDefinitions.h:491
femmodel
FemModel * femmodel
Definition: esmfbinders.cpp:16
Generated on Thu Jul 2 2020 08:09:23 for Ice Sheet System Model by   doxygen 1.8.19