Ice Sheet System Model  4.18
Code documentation
Radar.cpp
Go to the documentation of this file.
1 
5 /*Headers:*/
6 #ifdef HAVE_CONFIG_H
7  #include <config.h>
8 #else
9 #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
10 #endif
11 #include "./classes.h"
12 #include "./Elements/Element.h"
13 #include "./Elements/Elements.h"
16 #include "../datastructures/datastructures.h"
17 #include "./FemModel.h"
18 #include "../classes/Params/Parameters.h"
19 #include "../classes/gauss/Gauss.h"
20 #include "./Radar.h"
21 
22 /*Element macros*/
23 #define NUMVERTICES 6
24 #define NUMVERTICES2D 3
25 
26 /*Radar constructors, destructors :*/
27 Radar::Radar(){/*{{{*/
28  this->definitionenum = -1;
29  this->name = NULL;
30 }
31 /*}}}*/
32 Radar::Radar(char* in_name, int in_definitionenum){/*{{{*/
33  this->definitionenum=in_definitionenum;
34  this->name = xNew<char>(strlen(in_name)+1);
35  xMemCpy<char>(this->name,in_name,strlen(in_name)+1);
36 }
37 /*}}}*/
38 Radar::~Radar(){/*{{{*/
39  if(this->name)xDelete(this->name);
40 }
41 /*}}}*/
42 /*Object virtual function resolution: */
43 Object* Radar::copy() {/*{{{*/
44  Radar* out =new Radar(this->name,this->definitionenum);
45  return (Object*)out;
46 }
47 /*}}}*/
48 void Radar::DeepEcho(void){/*{{{*/
49  this->Echo();
50 }
51 /*}}}*/
52 void Radar::Echo(void){/*{{{*/
53  _printf_(" Radar: " << name << " " << this->definitionenum << "\n");
54 }
55 /*}}}*/
56 int Radar::Id(void){/*{{{*/
57  return -1;
58 }
59 /*}}}*/
60 void Radar::Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){/*{{{*/
61  _error_("not implemented yet!");
62 }
63 /*}}}*/
64 int Radar::ObjectEnum(void){/*{{{*/
65  return RadarEnum;
66 }
67 /*}}}*/
68 /*Definition virtual function resolutoin: */
69 int Radar::DefinitionEnum(){/*{{{*/
70  return this->definitionenum;
71 }
72 /*}}}*/
73 char* Radar::Name(){/*{{{*/
74  char* name2=xNew<char>(strlen(this->name)+1);
75  xMemCpy(name2,this->name,strlen(this->name)+1);
76  return name2;
77 }
78 /*}}}*/
80  int i;
81  for(i=0;i<femmodel->elements->Size();i++){
83  this->ComputeRadarAttenuation(element);
84  this->ComputeRadarPower(element);
85  }
86  return 0.;
87 }
88  /*}}}*/
90  //int numvertices = element->GetNumberOfVertices();
91  IssmDouble eps0=8.85418782e-12;
92  IssmDouble eps_ice=3.17;
93  IssmDouble e=2.7183;
94  IssmDouble c=299792458;
95  IssmDouble k=1.3806488e-23;
96  IssmDouble eVtoJ=1.6e-19;
97  IssmDouble Tr_M07=252.1500;
98  IssmDouble Tr_W97=258.1500;
99  IssmDouble sig_ice_M07=9.2;
100  IssmDouble sig_ice_W97=9;
101  IssmDouble cond_H_M07=3.2;
102  IssmDouble cond_H_W97=4;
103  IssmDouble cond_Cl_M07=0.43;
104  IssmDouble cond_Cl_W97=0.55;
105  IssmDouble cond_NH_M07=0.8;
106  IssmDouble cond_NH_W97=1;
107  IssmDouble E_M07=8.1600e-20;
108  IssmDouble E_W97=9.2800e-20;
109  IssmDouble E_H_M07=3.2000e-20;
110  IssmDouble E_H_W97=3.3600e-20;
111  IssmDouble E_Cl_M07=3.0400e-20;
112  IssmDouble E_Cl_W97=3.6800e-20;
113  IssmDouble E_NH=3.6800e-20;
114  IssmDouble mol_H_hol=1.6;
115  IssmDouble mol_H_lgp=0.2;
116  IssmDouble mol_Cl_hol=0.4;
117  IssmDouble mol_Cl_lgp=1.8;
118  IssmDouble mol_NH_hol=0.5;
119  IssmDouble mol_NH_lgp=0.4;
120  IssmDouble mol_H, mol_Cl, mol_NH;
121  IssmDouble attenuation_rate_macgregor[NUMVERTICES];
122  IssmDouble attenuation_rate_wolff[NUMVERTICES];
123  IssmDouble temperature, ice_period, attenuation_rate_M07_pureice, attenuation_rate_M07_H, attenuation_rate_M07_Cl, attenuation_rate_M07_NH;
124  IssmDouble attenuation_rate_W97_pureice, attenuation_rate_W97_H, attenuation_rate_W97_Cl, attenuation_rate_W97_NH;
125  IssmDouble m1, m2, m3, m4, m5, m6, m7, m8;
126  IssmDouble w2, w3, w4, w5, w6, w7, w8;
127  GaussPenta* gauss=NULL;
128 
129  /*Retrieve all inputs we will be needing: */
130  Input2* temp_input=element->GetInput2(TemperatureEnum); _assert_(temp_input);
131  Input2* ice_period_input=element->GetInput2(RadarIcePeriodEnum); _assert_(ice_period_input);
132 
133  /* Start looping on the number of vertices: */
134  gauss=new GaussPenta();
135 
136  for (int iv=0;iv<NUMVERTICES;iv++){
137  gauss->GaussVertex(iv);
138 
139  /*Get ice temperature: */
140  temp_input->GetInputValue(&temperature,gauss);
141  ice_period_input->GetInputValue(&ice_period,gauss);
142 
143  if(ice_period>0){;
144  mol_H=mol_H_hol;
145  mol_Cl=mol_Cl_hol;
146  mol_NH=mol_NH_hol;
147  }
148  else{
149  mol_H=mol_H_lgp;
150  mol_Cl=mol_Cl_lgp;
151  mol_NH=mol_NH_lgp;
152  }
153 
154  /*Compute M07 radar conductivity constant: */
155  m1=(10*log10(e))/(1000*eps0*sqrt(eps_ice)*c);
156  m2=E_M07/k;
157  m3=cond_H_M07*mol_H;
158  m4=E_H_M07/k;
159  m5=cond_Cl_M07*mol_Cl;
160  m6=E_Cl_M07/k;
161  m7=cond_NH_M07*mol_NH;
162  m8=E_NH/k;
163 
164  /*Compute MacGregor (M07) attenuation rate: */
165  attenuation_rate_M07_pureice=m1*sig_ice_M07*exp(m2*((1/Tr_M07)-(1/temperature)));
166  attenuation_rate_M07_H=m3*exp(m4*((1/Tr_M07)-(1/temperature)));
167  attenuation_rate_M07_Cl=m5*exp(m6*((1/Tr_M07)-(1/temperature)));
168  attenuation_rate_M07_NH=m7*exp(m8*((1/Tr_M07)-(1/temperature)));
169  attenuation_rate_macgregor[iv]=attenuation_rate_M07_pureice+attenuation_rate_M07_H+attenuation_rate_M07_Cl+attenuation_rate_M07_NH;
170 
171  /*Compute Wolff (W97) radar conductivity constant: */
172  w2=E_W97/k;
173  w3=cond_H_W97*mol_H;
174  w4=E_H_W97/k;
175  w5=cond_Cl_W97*mol_Cl;
176  w6=E_Cl_W97/k;
177  w7=cond_NH_W97*mol_NH;
178  w8=E_NH/k;
179 
180  /*Compute Wolff attenuation rate: */
181  attenuation_rate_W97_pureice=m1*sig_ice_W97*exp(w2*((1/Tr_W97)-(1/temperature)));
182  attenuation_rate_W97_H=w3*exp(w4*((1/Tr_W97)-(1/temperature)));
183  attenuation_rate_W97_Cl=w5*exp(w6*((1/Tr_W97)-(1/temperature)));
184  attenuation_rate_W97_NH=w7*exp(w8*((1/Tr_W97)-(1/temperature)));
185  attenuation_rate_wolff[iv]=attenuation_rate_W97_pureice+attenuation_rate_W97_H+attenuation_rate_W97_Cl+attenuation_rate_W97_NH;
186  }
187 
188  /*Add Attenuation rate results into inputs*/
189  element->AddInput2(RadarAttenuationMacGregorEnum,&attenuation_rate_macgregor[0],P1Enum);
190  element->AddInput2(RadarAttenuationWolffEnum,&attenuation_rate_wolff[0],P1Enum);
191 
192  /*Clean up*/
193  delete gauss;
194 
195 }/*}}}*/
196 void Radar::ComputeRadarPower(Element* element){/*{{{*/
197 
198  IssmDouble *xyz_list=NULL;
199  IssmDouble power_M07[NUMVERTICES];
200  IssmDouble power_W97[NUMVERTICES];
201  IssmDouble depth[NUMVERTICES];
202  IssmDouble aircraft_elev=0.5;
203  IssmDouble eps_ice=3.15;
204  IssmDouble t_tp=273.15; /* triple point temperature [K] */
205  IssmDouble p_tp=611.73; /* water pressure [Pa] */
206  IssmDouble gamma=7.4200e-07; /* Clausius-Clapeyron constant [K/kPa] */
207  IssmDouble attenuation_rate_macgregor, attenuation_rate_wolff, attenuation_total_M07, attenuation_total_W97;
208  IssmDouble thickness, surface, z, temperature, geometric_loss, reflectivity;
209  IssmDouble rho_ice, gravity, pressure, pressure_melting_pt, frozen_temp, basal_temp, basal_pmp;
210  GaussPenta* gauss=NULL;
211 
212  /* Get node coordinates*/
213  element->GetVerticesCoordinates(&xyz_list);
214  Input2 *atten_input_M07 = element->GetInput2(RadarAttenuationMacGregorEnum); _assert_(atten_input_M07);
215  Input2 *atten_input_W97 = element->GetInput2(RadarAttenuationWolffEnum); _assert_(atten_input_W97);
216  Input2 *surf_input = element->GetInput2(SurfaceEnum); _assert_(surf_input);
217  Input2 *thick_input = element->GetInput2(ThicknessEnum); _assert_(thick_input);
218  Input2 *temp_input = element->GetInput2(TemperatureEnum); _assert_(temp_input);
219 
220  /* Start looping on the number of vertices: */
221  gauss=new GaussPenta();
222  for (int iv=0;iv<NUMVERTICES;iv++){
223  gauss->GaussVertex(iv);
224 
225  /*Get all the inputs: */
226  atten_input_M07->GetInputValue(&attenuation_rate_macgregor,gauss);
227  atten_input_W97->GetInputValue(&attenuation_rate_wolff,gauss);
228  thick_input->GetInputValue(&thickness,gauss);
229  temp_input->GetInputValue(&temperature,gauss);
230  surf_input->GetInputValue(&surface,gauss);
231 
232  /*Compute depth below the ice surface: */
233  z=xyz_list[3*iv+2];
234  depth[iv]=(surface-z)/1e3;
235 
236  /*Compute total attenuation: */
237  attenuation_total_M07=attenuation_rate_macgregor*depth[iv];
238  attenuation_total_W97=attenuation_rate_wolff*depth[iv];
239 
240  /*Compute geometric loss: */
241  geometric_loss=10*log10((depth[iv]+aircraft_elev)/sqrt(eps_ice));
242 
243  /*Compute radar power: */
244  power_M07[iv]=-geometric_loss-attenuation_total_M07;
245  power_W97[iv]=-geometric_loss-attenuation_total_W97;
246 
247  /*Identify basal elements: */
248  if(element->IsOnBase() && iv<NUMVERTICES2D){
249 
250  /*Compute pressure melting point: */
251  rho_ice=element->FindParam(MaterialsLatentheatEnum);
252  gravity=element->FindParam(ConstantsGEnum);
253  pressure=rho_ice*gravity*thickness;
254  pressure_melting_pt=t_tp-gamma*(pressure-p_tp);
255 
256  if((temperature-pressure_melting_pt)<=-1){
257  reflectivity=-40;
258  }
259  else if((temperature-pressure_melting_pt)>-1 && (temperature-pressure_melting_pt)<0){
260  reflectivity=0;
261  }
262  else{
263  reflectivity=70;
264  }
265  power_M07[iv]=power_M07[iv]+reflectivity;
266  power_W97[iv]=power_W97[iv]+reflectivity;
267  }
268  }
269 
270  /*Add power results into inputs*/
271  element->AddInput2(RadarPowerMacGregorEnum,&power_M07[0],P1Enum);
272  element->AddInput2(RadarPowerWolffEnum,&power_W97[0],P1Enum);
273 
274  /*Clean up and return*/
275  delete gauss;
276 }/*}}}*/
RadarPowerWolffEnum
@ RadarPowerWolffEnum
Definition: EnumDefinitions.h:670
DataSet::Size
int Size()
Definition: DataSet.cpp:399
Radar::Radar
Radar()
Definition: Radar.cpp:27
Radar::name
char * name
Definition: Radar.h:18
_assert_
#define _assert_(ignore)
Definition: exceptions.h:37
GaussPenta::GaussVertex
void GaussVertex(int iv)
Definition: GaussPenta.cpp:785
IssmDouble
double IssmDouble
Definition: types.h:37
Element::IsOnBase
bool IsOnBase()
Definition: Element.cpp:1984
Element::FindParam
void FindParam(bool *pvalue, int paramenum)
Definition: Element.cpp:933
_printf_
#define _printf_(StreamArgs)
Definition: Print.h:22
RadarAttenuationMacGregorEnum
@ RadarAttenuationMacGregorEnum
Definition: EnumDefinitions.h:666
Radar::definitionenum
int definitionenum
Definition: Radar.h:19
MaterialsLatentheatEnum
@ MaterialsLatentheatEnum
Definition: EnumDefinitions.h:254
Radar
Definition: Radar.h:15
Radar.h
: header file for Radar object
Radar::ObjectEnum
int ObjectEnum(void)
Definition: Radar.cpp:64
Element::GetInput2
virtual Input2 * GetInput2(int inputenum)=0
RadarAttenuationWolffEnum
@ RadarAttenuationWolffEnum
Definition: EnumDefinitions.h:667
Element
Definition: Element.h:41
Elements.h
P1Enum
@ P1Enum
Definition: EnumDefinitions.h:662
Radar::Echo
void Echo(void)
Definition: Radar.cpp:52
Element.h
abstract class for Element object This class is a place holder for the Tria and the Penta elements....
NUMVERTICES
#define NUMVERTICES
Definition: Radar.cpp:23
RadarPowerMacGregorEnum
@ RadarPowerMacGregorEnum
Definition: EnumDefinitions.h:669
Object
Definition: Object.h:13
ConstantsGEnum
@ ConstantsGEnum
Definition: EnumDefinitions.h:102
Radar::copy
Object * copy()
Definition: Radar.cpp:43
Element::AddInput2
virtual void AddInput2(int input_enum, IssmDouble *values, int interpolation_enum)
Definition: Element.h:216
xDelete
void xDelete(T **&aT_pp)
Definition: MemOps.h:97
GaussPenta
Definition: GaussPenta.h:13
ExternalResult.h
abstract class for ExternalResult object
NUMVERTICES2D
#define NUMVERTICES2D
Definition: Radar.cpp:24
Element::GetVerticesCoordinates
void GetVerticesCoordinates(IssmDouble **xyz_list)
Definition: Element.cpp:1446
Radar::Name
char * Name()
Definition: Radar.cpp:73
SurfaceEnum
@ SurfaceEnum
Definition: EnumDefinitions.h:823
Radar::Marshall
void Marshall(char **pmarshalled_data, int *pmarshalled_data_size, int marshall_direction)
Definition: Radar.cpp:60
FemModel::elements
Elements * elements
Definition: FemModel.h:44
Input2
Definition: Input2.h:18
FemModel
Definition: FemModel.h:31
Radar::Response
IssmDouble Response(FemModel *femmodel)
Definition: Radar.cpp:79
TemperatureEnum
@ TemperatureEnum
Definition: EnumDefinitions.h:831
Radar::ComputeRadarAttenuation
void ComputeRadarAttenuation(Element *element)
Definition: Radar.cpp:89
RadarEnum
@ RadarEnum
Definition: EnumDefinitions.h:665
Radar::~Radar
~Radar()
Definition: Radar.cpp:38
Radar::DeepEcho
void DeepEcho(void)
Definition: Radar.cpp:48
RadarIcePeriodEnum
@ RadarIcePeriodEnum
Definition: EnumDefinitions.h:668
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
DataSet::GetObjectByOffset
Object * GetObjectByOffset(int offset)
Definition: DataSet.cpp:334
FemModel.h
ThicknessEnum
@ ThicknessEnum
Definition: EnumDefinitions.h:840
Input2::GetInputValue
virtual void GetInputValue(IssmDouble *pvalue, Gauss *gauss)
Definition: Input2.h:38
Radar::Id
int Id(void)
Definition: Radar.cpp:56
xMemCpy
T * xMemCpy(T *dest, const T *src, unsigned int size)
Definition: MemOps.h:152
Radar::ComputeRadarPower
void ComputeRadarPower(Element *element)
Definition: Radar.cpp:196
Results.h
classes.h
Radar::DefinitionEnum
int DefinitionEnum()
Definition: Radar.cpp:69
femmodel
FemModel * femmodel
Definition: esmfbinders.cpp:16