Ice Sheet System Model  4.18
Code documentation
Functions
OceanExchangeDatax.cpp File Reference
#include "../../shared/shared.h"
#include "../../toolkits/toolkits.h"
#include "../../classes/classes.h"
#include "../modules.h"
#include "./OceanExchangeDatax.h"

Go to the source code of this file.

Functions

void OceanExchangeDatax (FemModel *femmodel, bool init_stage)
 

Function Documentation

◆ OceanExchangeDatax()

void OceanExchangeDatax ( FemModel *  femmodel,
bool  init_stage 
)

Definition at line 11 of file OceanExchangeDatax.cpp.

11  {
12 
13  #ifndef _HAVE_AD_
14  if(VerboseSolution()) _printf0_(" ocean coupling: exchanging information\n");
15  int my_rank;
16  ISSM_MPI_Comm tomitgcmcomm;
17  ISSM_MPI_Status status;
18 
19  my_rank=IssmComm::GetRank();
20  GenericParam<ISSM_MPI_Comm>* parcom = dynamic_cast<GenericParam<ISSM_MPI_Comm>*>(femmodel->parameters->FindParamObject(ToMITgcmCommEnum));
21  if(!parcom)_error_("TransferForcing error message: could not find ToMITgcmCommEnum communicator");
22  tomitgcmcomm=parcom->GetParameterValue();
23 
24  int oceangridnxsize,oceangridnysize,ngrids_ocean,nels_ocean;
25  IssmDouble oceantime,coupling_time,time,yts;
26  IssmDouble rho_ice;
27  IssmDouble *oceanmelt = NULL;
28  IssmDouble *oceangridx;
29  IssmDouble *oceangridy;
30  IssmDouble *icebase_oceangrid = NULL;
31  IssmDouble *icemask_oceangrid = NULL;
32  IssmDouble* x_ice = NULL;
33  IssmDouble* y_ice = NULL;
34  IssmDouble* lat_ice = NULL;
35  IssmDouble* lon_ice = NULL;
36  IssmDouble* icebase = NULL;
37  IssmDouble* icemask = NULL;
38  IssmDouble* melt_mesh = NULL;
39  int* index_ice = NULL;
40  int* index_ocean = NULL;
41  int ngrids_ice=femmodel->vertices->NumberOfVertices();
42  int nels_ice=femmodel->elements->NumberOfElements();
43 
44  /*Recover fixed parameters and store them*/
45  femmodel->parameters->FindParam(&coupling_time,TimesteppingCouplingTimeEnum);
46  femmodel->parameters->FindParam(&time,TimeEnum);
47 
48  /*Exchange or recover mesh and inputs needed*/
49  if(init_stage==true){
50  if(my_rank==0){
51  ISSM_MPI_Send(&coupling_time,1,ISSM_MPI_DOUBLE,0,10001000,tomitgcmcomm);
52  ISSM_MPI_Recv(&oceangridnxsize,1,ISSM_MPI_INT,0,10001003,tomitgcmcomm,&status);
53  ISSM_MPI_Recv(&oceangridnysize,1,ISSM_MPI_INT,0,10001004,tomitgcmcomm,&status);
54  }
55  ngrids_ocean=oceangridnxsize*oceangridnysize;
56  ISSM_MPI_Bcast(&oceangridnxsize,1,ISSM_MPI_INT,0,IssmComm::GetComm());
57  ISSM_MPI_Bcast(&oceangridnysize,1,ISSM_MPI_INT,0,IssmComm::GetComm());
58  ISSM_MPI_Bcast(&ngrids_ocean,1,ISSM_MPI_INT,0,IssmComm::GetComm());
59  ISSM_MPI_Bcast(&oceantime,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
60  femmodel->parameters->SetParam(oceangridnxsize,OceanGridNxEnum);
61  femmodel->parameters->SetParam(oceangridnysize,OceanGridNyEnum);
62  if(my_rank==0){
63  oceangridx = xNew<IssmDouble>(ngrids_ocean);
64  ISSM_MPI_Recv(oceangridx,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001005,tomitgcmcomm,&status);
65  oceangridy = xNew<IssmDouble>(ngrids_ocean);
66  ISSM_MPI_Recv(oceangridy,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001006,tomitgcmcomm,&status);
67 
68  /*Exchange varying parameters for the initialization*/
69  ISSM_MPI_Send(&time,1,ISSM_MPI_DOUBLE,0,10001001,tomitgcmcomm);
70  ISSM_MPI_Recv(&oceantime,1,ISSM_MPI_DOUBLE,0,10001002,tomitgcmcomm,&status);
71  }
72  if(my_rank!=0){
73  oceangridx=xNew<IssmDouble>(ngrids_ocean);
74  oceangridy=xNew<IssmDouble>(ngrids_ocean);
75  }
76  ISSM_MPI_Bcast(oceangridx,ngrids_ocean,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
77  ISSM_MPI_Bcast(oceangridy,ngrids_ocean,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
78  femmodel->parameters->SetParam(oceangridx,ngrids_ocean,OceanGridXEnum);
79  femmodel->parameters->SetParam(oceangridy,ngrids_ocean,OceanGridYEnum);
80  }
81  else{
82  femmodel->parameters->FindParam(&oceangridx,&ngrids_ocean,OceanGridXEnum);
83  femmodel->parameters->FindParam(&oceangridy,&ngrids_ocean,OceanGridYEnum);
84  }
85 
86  /*Interpolate ice base and mask onto ocean grid*/
87  femmodel->GetMesh(femmodel->vertices,femmodel->elements,&x_ice,&y_ice,&index_ice);
88  BamgTriangulatex(&index_ocean,&nels_ocean,oceangridx,oceangridy,ngrids_ocean);
89  femmodel->vertices->LatLonList(&lat_ice,&lon_ice);
90  GetVectorFromInputsx(&icebase,femmodel,BaseEnum,VertexSIdEnum);
91  Options* options = new Options();
92  GenericOption<double> *odouble = new GenericOption<double>();
93  const char* name = "default";
94  odouble->name =xNew<char>(strlen(name)+1);
95  memcpy(odouble->name,name,(strlen(name)+1)*sizeof(char));
96  odouble->value=+9999.;
97  odouble->size[0]=1;
98  odouble->size[1]=1;
99  options->AddOption(odouble);
100  InterpFromMeshToMesh2dx(&icebase_oceangrid,index_ice,lon_ice,lat_ice,ngrids_ice,nels_ice,
101  icebase,ngrids_ice,1,oceangridx,oceangridy,ngrids_ocean,options);
102  delete options;
103  xDelete<IssmDouble>(icebase);
104 
105  GetVectorFromInputsx(&icemask,femmodel,MaskIceLevelsetEnum,VertexSIdEnum);
106  Options* options2 = new Options();
107  GenericOption<double> *odouble2 = new GenericOption<double>();
108  const char* name2 = "default";
109  odouble2->name =xNew<char>(strlen(name2)+1);
110  memcpy(odouble2->name,name2,(strlen(name2)+1)*sizeof(char));
111  odouble2->value=+1.;
112  odouble2->size[0]=1;
113  odouble2->size[1]=1;
114  options2->AddOption(odouble2);
115  InterpFromMeshToMesh2dx(&icemask_oceangrid,index_ice,lon_ice,lat_ice,ngrids_ice,nels_ice,
116  icemask,ngrids_ice,1,oceangridx,oceangridy,ngrids_ocean,options2);
117  delete options2;
118  xDelete<IssmDouble>(icemask);
119 
120  /*Put +9999 for places where there is no ice!*/
121  for(int i=0;i<ngrids_ocean;i++) if(icemask_oceangrid[i]>0.) icebase_oceangrid[i]=+9999.;
122  xDelete<IssmDouble>(icemask_oceangrid);
123 
124  if(init_stage==true){ //just send icebase
125  if(my_rank==0){
126  ISSM_MPI_Send(icebase_oceangrid,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001008,tomitgcmcomm);
127  }
128  }
129  else{ //send and receive exchanged data
130  femmodel->parameters->FindParam(&rho_ice,MaterialsRhoIceEnum);
131  femmodel->parameters->FindParam(&yts,ConstantsYtsEnum);
132  if(my_rank==0){
133  ISSM_MPI_Send(&time,1,ISSM_MPI_DOUBLE,0,10001001,tomitgcmcomm);
134  ISSM_MPI_Recv(&oceantime,1,ISSM_MPI_DOUBLE,0,10001002,tomitgcmcomm,&status);
135  if((oceantime - time > 0.1*yts) & (oceantime - time < -0.1*yts)) _error_("Ocean and ice time are starting to diverge");
136  oceanmelt = xNew<IssmDouble>(ngrids_ocean);
137  ISSM_MPI_Recv(oceanmelt,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001007,tomitgcmcomm,&status);
138  ISSM_MPI_Send(icebase_oceangrid,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001008,tomitgcmcomm);
139  }
140  ISSM_MPI_Bcast(&oceantime,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
141  if(my_rank!=0) oceanmelt=xNew<IssmDouble>(ngrids_ocean);
142  ISSM_MPI_Bcast(oceanmelt,ngrids_ocean,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
143 
144  /*Interp melt onto ice grid*/
145  InterpFromMeshToMesh2dx(&melt_mesh,index_ocean,oceangridx,oceangridy,ngrids_ocean,nels_ocean,
146  oceanmelt,ngrids_ocean,1,
147  lon_ice,lat_ice,ngrids_ice,NULL);
148 
149  for(int i=0;i<ngrids_ice;i++) melt_mesh[i]=-melt_mesh[i]/rho_ice; //heat flux provided by ocean is in kg/m^2/s
150  InputUpdateFromVectorx(femmodel,melt_mesh,BasalforcingsFloatingiceMeltingRateEnum,VertexSIdEnum);
151  }
152 
153  /*Delete*/
154  xDelete<int>(index_ice);
155  xDelete<int>(index_ocean);
156  xDelete<IssmDouble>(lat_ice);
157  xDelete<IssmDouble>(lon_ice);
158  xDelete<IssmDouble>(x_ice);
159  xDelete<IssmDouble>(y_ice);
160  xDelete<IssmDouble>(icebase_oceangrid);
161  xDelete<IssmDouble>(oceangridx);
162  xDelete<IssmDouble>(oceangridy);
163  xDelete<IssmDouble>(melt_mesh);
164  xDelete<IssmDouble>(oceanmelt);
165  #else
166  _error_("not supported");
167  #endif
168 }
OceanGridXEnum
@ OceanGridXEnum
Definition: EnumDefinitions.h:279
GenericOption::size
int size[2]
Definition: GenericOption.h:27
BaseEnum
@ BaseEnum
Definition: EnumDefinitions.h:495
Vertices::LatLonList
void LatLonList(IssmDouble **lat, IssmDouble **lon)
Definition: Vertices.cpp:135
Options
Definition: Options.h:9
FemModel::vertices
Vertices * vertices
Definition: FemModel.h:49
IssmDouble
double IssmDouble
Definition: types.h:37
_printf0_
#define _printf0_(StreamArgs)
Definition: Print.h:29
GenericOption::name
char * name
Definition: GenericOption.h:25
MaskIceLevelsetEnum
@ MaskIceLevelsetEnum
Definition: EnumDefinitions.h:641
FemModel::parameters
Parameters * parameters
Definition: FemModel.h:46
TimeEnum
@ TimeEnum
Definition: EnumDefinitions.h:427
IssmComm::GetComm
static ISSM_MPI_Comm GetComm(void)
Definition: IssmComm.cpp:30
GetVectorFromInputsx
void GetVectorFromInputsx(IssmDouble **pvector, int *pvector_size, FemModel *femmodel, int name)
Definition: GetVectorFromInputsx.cpp:81
Vertices::NumberOfVertices
int NumberOfVertices(void)
Definition: Vertices.cpp:255
MaterialsRhoIceEnum
@ MaterialsRhoIceEnum
Definition: EnumDefinitions.h:264
ConstantsYtsEnum
@ ConstantsYtsEnum
Definition: EnumDefinitions.h:104
GenericOption::value
OptionType value
Definition: GenericOption.h:26
OceanGridNxEnum
@ OceanGridNxEnum
Definition: EnumDefinitions.h:277
BamgTriangulatex
int BamgTriangulatex(int **pindex, int *pnels, double *x, double *y, int nods)
Definition: BamgTriangulatex.cpp:14
ISSM_MPI_DOUBLE
#define ISSM_MPI_DOUBLE
Definition: issmmpi.h:125
Parameters::SetParam
void SetParam(bool boolean, int enum_type)
Definition: Parameters.cpp:441
ISSM_MPI_INT
#define ISSM_MPI_INT
Definition: issmmpi.h:127
GenericParam::GetParameterValue
P & GetParameterValue()
Definition: GenericParam.h:62
VertexSIdEnum
@ VertexSIdEnum
Definition: EnumDefinitions.h:1325
InputUpdateFromVectorx
void InputUpdateFromVectorx(FemModel *femmodel, Vector< IssmDouble > *vector, int name, int type)
Definition: InputUpdateFromVectorx.cpp:9
BasalforcingsFloatingiceMeltingRateEnum
@ BasalforcingsFloatingiceMeltingRateEnum
Definition: EnumDefinitions.h:476
ISSM_MPI_Send
int ISSM_MPI_Send(void *buf, int count, ISSM_MPI_Datatype datatype, int dest, int tag, ISSM_MPI_Comm comm)
Definition: issmmpi.cpp:484
ISSM_MPI_Status
int ISSM_MPI_Status
Definition: issmmpi.h:121
FemModel::elements
Elements * elements
Definition: FemModel.h:44
ISSM_MPI_Bcast
int ISSM_MPI_Bcast(void *buffer, int count, ISSM_MPI_Datatype datatype, int root, ISSM_MPI_Comm comm)
Definition: issmmpi.cpp:162
TimesteppingCouplingTimeEnum
@ TimesteppingCouplingTimeEnum
Definition: EnumDefinitions.h:429
ToMITgcmCommEnum
@ ToMITgcmCommEnum
Definition: EnumDefinitions.h:437
OceanGridNyEnum
@ OceanGridNyEnum
Definition: EnumDefinitions.h:278
IssmComm::GetRank
static int GetRank(void)
Definition: IssmComm.cpp:34
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
VerboseSolution
bool VerboseSolution(void)
Definition: Verbosity.cpp:24
ISSM_MPI_Comm
int ISSM_MPI_Comm
Definition: issmmpi.h:118
Parameters::FindParam
void FindParam(bool *pinteger, int enum_type)
Definition: Parameters.cpp:262
Elements::NumberOfElements
int NumberOfElements(void)
Definition: Elements.cpp:67
ISSM_MPI_Recv
int ISSM_MPI_Recv(void *buf, int count, ISSM_MPI_Datatype datatype, int source, int tag, ISSM_MPI_Comm comm, ISSM_MPI_Status *status)
Definition: issmmpi.cpp:342
Parameters::FindParamObject
Param * FindParamObject(int enum_type)
Definition: Parameters.cpp:588
Options::AddOption
int AddOption(Option *in_oobject)
Definition: Options.cpp:33
InterpFromMeshToMesh2dx
int InterpFromMeshToMesh2dx(double **pdata_interp, int *index_data, double *x_data, double *y_data, int nods_data, int nels_data, double *data, int M_data, int N_data, double *x_interp, double *y_interp, int N_interp, Options *options)
Definition: InterpFromMeshToMesh2dx.cpp:16
FemModel::GetMesh
void GetMesh(Vertices *femmodel_vertices, Elements *femmodel_elements, IssmDouble **px, IssmDouble **py, int **pelementslist)
Definition: FemModel.cpp:3731
OceanGridYEnum
@ OceanGridYEnum
Definition: EnumDefinitions.h:280
GenericParam
Definition: GenericParam.h:26
GenericOption
Definition: GenericOption.h:22
femmodel
FemModel * femmodel
Definition: esmfbinders.cpp:16
Generated on Thu Jul 2 2020 08:09:24 for Ice Sheet System Model by   doxygen 1.8.19