Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 23713) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 23714) @@ -520,4 +520,9 @@ endif #}}} +#Oceansources {{{ +if OCEAN +issm_sources += ./modules/OceanExchangeDatax/OceanExchangeDatax.cpp +endif +#}}} #Slr sources {{{ if SEALEVELRISE Index: /issm/trunk-jpl/src/c/cores/transient_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/transient_core.cpp (revision 23713) +++ /issm/trunk-jpl/src/c/cores/transient_core.cpp (revision 23714) @@ -79,130 +79,11 @@ #endif - if(isoceancoupling){ /*{{{*/ - #ifndef _HAVE_AD_ - if(VerboseSolution()) _printf0_(" ocean coupling: initialization \n"); - int my_rank; - ISSM_MPI_Comm tomitgcmcomm; - ISSM_MPI_Status status; - - my_rank=IssmComm::GetRank(); - GenericParam* parcom = dynamic_cast*>(femmodel->parameters->FindParamObject(ToMITgcmCommEnum)); - if(!parcom)_error_("TransferForcing error message: could not find ToMITgcmCommEnum communicator"); - tomitgcmcomm=parcom->GetParameterValue(); - - int oceangridnxsize,oceangridnysize,ngrids_ocean,nels_ocean; - IssmDouble oceantime,coupling_time; - IssmDouble *oceangridx; - IssmDouble *oceangridy; - IssmDouble *icebase_oceangrid = NULL; - IssmDouble *icemask_oceangrid = NULL; - IssmDouble* x_ice = NULL; - IssmDouble* y_ice = NULL; - IssmDouble* lat_ice = NULL; - IssmDouble* lon_ice = NULL; - IssmDouble* icebase = NULL; - IssmDouble* icemask = NULL; - int* index_ice = NULL; - int* index_ocean = NULL; - int ngrids_ice=femmodel->vertices->NumberOfVertices(); - int nels_ice=femmodel->elements->NumberOfElements(); - - /*Recover fixed parameters and store them*/ - femmodel->parameters->FindParam(&coupling_time,TimesteppingCouplingTimeEnum); - if(my_rank==0){ - ISSM_MPI_Send(&coupling_time,1,ISSM_MPI_DOUBLE,0,10001000,tomitgcmcomm); - ISSM_MPI_Recv(&oceangridnxsize,1,ISSM_MPI_INT,0,10001003,tomitgcmcomm,&status); - ISSM_MPI_Recv(&oceangridnysize,1,ISSM_MPI_INT,0,10001004,tomitgcmcomm,&status); - } - ngrids_ocean=oceangridnxsize*oceangridnysize; - ISSM_MPI_Bcast(&oceangridnxsize,1,ISSM_MPI_INT,0,IssmComm::GetComm()); - ISSM_MPI_Bcast(&oceangridnysize,1,ISSM_MPI_INT,0,IssmComm::GetComm()); - ISSM_MPI_Bcast(&ngrids_ocean,1,ISSM_MPI_INT,0,IssmComm::GetComm()); - ISSM_MPI_Bcast(&oceantime,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - femmodel->parameters->SetParam(oceangridnxsize,OceanGridNxEnum); - femmodel->parameters->SetParam(oceangridnysize,OceanGridNyEnum); - if(my_rank==0){ - oceangridx = xNew(ngrids_ocean); - ISSM_MPI_Recv(oceangridx,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001005,tomitgcmcomm,&status); - oceangridy = xNew(ngrids_ocean); - ISSM_MPI_Recv(oceangridy,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001006,tomitgcmcomm,&status); - - /*Exchange varying parameters for the initialization*/ - ISSM_MPI_Send(&time,1,ISSM_MPI_DOUBLE,0,10001001,tomitgcmcomm); - ISSM_MPI_Recv(&oceantime,1,ISSM_MPI_DOUBLE,0,10001002,tomitgcmcomm,&status); - } - if(my_rank!=0){ - oceangridx=xNew(ngrids_ocean); - oceangridy=xNew(ngrids_ocean); - } - ISSM_MPI_Bcast(oceangridx,ngrids_ocean,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - ISSM_MPI_Bcast(oceangridy,ngrids_ocean,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - femmodel->parameters->SetParam(oceangridx,ngrids_ocean,OceanGridXEnum); - femmodel->parameters->SetParam(oceangridy,ngrids_ocean,OceanGridYEnum); - - /*Interpolate ice base and mask onto ocean grid*/ - femmodel->GetMesh(femmodel->vertices,femmodel->elements,&x_ice,&y_ice,&index_ice); - BamgTriangulatex(&index_ocean,&nels_ocean,oceangridx,oceangridy,ngrids_ocean); - femmodel->vertices->LatLonList(&lat_ice,&lon_ice); - GetVectorFromInputsx(&icebase,femmodel,BaseEnum,VertexSIdEnum); - Options* options = new Options(); - GenericOption *odouble = new GenericOption(); - const char* name = "default"; - odouble->name =xNew(strlen(name)+1); - memcpy(odouble->name,name,(strlen(name)+1)*sizeof(char)); - odouble->value=+9999.; - odouble->size[0]=1; - odouble->size[1]=1; - options->AddOption(odouble); - InterpFromMeshToMesh2dx(&icebase_oceangrid,index_ice,lon_ice,lat_ice,ngrids_ice,nels_ice, - icebase,ngrids_ice,1, - oceangridx,oceangridy,ngrids_ocean,options); - delete options; - xDelete(icebase); - - GetVectorFromInputsx(&icemask,femmodel,MaskIceLevelsetEnum,VertexSIdEnum); - Options* options2 = new Options(); - GenericOption *odouble2 = new GenericOption(); - const char* name2 = "default"; - odouble2->name =xNew(strlen(name2)+1); - memcpy(odouble2->name,name2,(strlen(name2)+1)*sizeof(char)); - odouble2->value=+1.; - odouble2->size[0]=1; - odouble2->size[1]=1; - options2->AddOption(odouble2); - InterpFromMeshToMesh2dx(&icemask_oceangrid,index_ice,lon_ice,lat_ice,ngrids_ice,nels_ice, - icemask,ngrids_ice,1,oceangridx,oceangridy,ngrids_ocean,options2); - delete options2; - xDelete(icemask); - - /*Put +9999 for places where there is no ice!*/ - for(int i=0;i0.) icebase_oceangrid[i]=+9999.; - xDelete(icemask_oceangrid); - - if(my_rank==0){ - ISSM_MPI_Send(icebase_oceangrid,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001008,tomitgcmcomm); - } - - /*Delete*/ - xDelete(index_ice); - xDelete(index_ocean); - xDelete(lat_ice); - xDelete(lon_ice); - xDelete(x_ice); - xDelete(y_ice); - xDelete(icebase_oceangrid); - xDelete(oceangridx); - xDelete(oceangridy); - #else - _error_("not supported"); - #endif - } - /*}}}*/ - - IssmDouble output_value; - int num_dependents; - IssmPDouble *dependents; - DataSet* dependent_objects=NULL; - IssmDouble J=0.; + if(isoceancoupling) OceanExchangeDatax(femmodel); + + IssmDouble output_value; + int num_dependents; + IssmPDouble *dependents; + DataSet* dependent_objects=NULL; + IssmDouble J=0.; if(iscontrol && isautodiff){ Index: /issm/trunk-jpl/src/c/modules/OceanExchangeDatax/OceanExchangeDatax.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/OceanExchangeDatax/OceanExchangeDatax.cpp (revision 23714) +++ /issm/trunk-jpl/src/c/modules/OceanExchangeDatax/OceanExchangeDatax.cpp (revision 23714) @@ -0,0 +1,130 @@ +/*!\file OceanExchangeDatax + * \brief: exchange of data with ocean model + */ + +#include "../../shared/shared.h" +#include "../../toolkits/toolkits.h" +#include "../../classes/classes.h" +#include "../modules.h" +#include "./OceanExchangeDatax.h" + +void OceanExchangeDatax(FemModel* femmodel){ + + #ifndef _HAVE_AD_ + if(VerboseSolution()) _printf0_(" ocean coupling: initialization \n"); + int my_rank; + ISSM_MPI_Comm tomitgcmcomm; + ISSM_MPI_Status status; + + my_rank=IssmComm::GetRank(); + GenericParam* parcom = dynamic_cast*>(femmodel->parameters->FindParamObject(ToMITgcmCommEnum)); + if(!parcom)_error_("TransferForcing error message: could not find ToMITgcmCommEnum communicator"); + tomitgcmcomm=parcom->GetParameterValue(); + + int oceangridnxsize,oceangridnysize,ngrids_ocean,nels_ocean; + IssmDouble oceantime,coupling_time,time; + IssmDouble *oceangridx; + IssmDouble *oceangridy; + IssmDouble *icebase_oceangrid = NULL; + IssmDouble *icemask_oceangrid = NULL; + IssmDouble* x_ice = NULL; + IssmDouble* y_ice = NULL; + IssmDouble* lat_ice = NULL; + IssmDouble* lon_ice = NULL; + IssmDouble* icebase = NULL; + IssmDouble* icemask = NULL; + int* index_ice = NULL; + int* index_ocean = NULL; + int ngrids_ice=femmodel->vertices->NumberOfVertices(); + int nels_ice=femmodel->elements->NumberOfElements(); + + /*Recover fixed parameters and store them*/ + femmodel->parameters->FindParam(&coupling_time,TimesteppingCouplingTimeEnum); + if(my_rank==0){ + ISSM_MPI_Send(&coupling_time,1,ISSM_MPI_DOUBLE,0,10001000,tomitgcmcomm); + ISSM_MPI_Recv(&oceangridnxsize,1,ISSM_MPI_INT,0,10001003,tomitgcmcomm,&status); + ISSM_MPI_Recv(&oceangridnysize,1,ISSM_MPI_INT,0,10001004,tomitgcmcomm,&status); + } + ngrids_ocean=oceangridnxsize*oceangridnysize; + ISSM_MPI_Bcast(&oceangridnxsize,1,ISSM_MPI_INT,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(&oceangridnysize,1,ISSM_MPI_INT,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(&ngrids_ocean,1,ISSM_MPI_INT,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(&oceantime,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + femmodel->parameters->SetParam(oceangridnxsize,OceanGridNxEnum); + femmodel->parameters->SetParam(oceangridnysize,OceanGridNyEnum); + if(my_rank==0){ + oceangridx = xNew(ngrids_ocean); + ISSM_MPI_Recv(oceangridx,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001005,tomitgcmcomm,&status); + oceangridy = xNew(ngrids_ocean); + ISSM_MPI_Recv(oceangridy,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001006,tomitgcmcomm,&status); + + /*Exchange varying parameters for the initialization*/ + ISSM_MPI_Send(&time,1,ISSM_MPI_DOUBLE,0,10001001,tomitgcmcomm); + ISSM_MPI_Recv(&oceantime,1,ISSM_MPI_DOUBLE,0,10001002,tomitgcmcomm,&status); + } + if(my_rank!=0){ + oceangridx=xNew(ngrids_ocean); + oceangridy=xNew(ngrids_ocean); + } + ISSM_MPI_Bcast(oceangridx,ngrids_ocean,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(oceangridy,ngrids_ocean,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + femmodel->parameters->SetParam(oceangridx,ngrids_ocean,OceanGridXEnum); + femmodel->parameters->SetParam(oceangridy,ngrids_ocean,OceanGridYEnum); + + /*Interpolate ice base and mask onto ocean grid*/ + femmodel->GetMesh(femmodel->vertices,femmodel->elements,&x_ice,&y_ice,&index_ice); + BamgTriangulatex(&index_ocean,&nels_ocean,oceangridx,oceangridy,ngrids_ocean); + femmodel->vertices->LatLonList(&lat_ice,&lon_ice); + GetVectorFromInputsx(&icebase,femmodel,BaseEnum,VertexSIdEnum); + Options* options = new Options(); + GenericOption *odouble = new GenericOption(); + const char* name = "default"; + odouble->name =xNew(strlen(name)+1); + memcpy(odouble->name,name,(strlen(name)+1)*sizeof(char)); + odouble->value=+9999.; + odouble->size[0]=1; + odouble->size[1]=1; + options->AddOption(odouble); + InterpFromMeshToMesh2dx(&icebase_oceangrid,index_ice,lon_ice,lat_ice,ngrids_ice,nels_ice, + icebase,ngrids_ice,1, + oceangridx,oceangridy,ngrids_ocean,options); + delete options; + xDelete(icebase); + + GetVectorFromInputsx(&icemask,femmodel,MaskIceLevelsetEnum,VertexSIdEnum); + Options* options2 = new Options(); + GenericOption *odouble2 = new GenericOption(); + const char* name2 = "default"; + odouble2->name =xNew(strlen(name2)+1); + memcpy(odouble2->name,name2,(strlen(name2)+1)*sizeof(char)); + odouble2->value=+1.; + odouble2->size[0]=1; + odouble2->size[1]=1; + options2->AddOption(odouble2); + InterpFromMeshToMesh2dx(&icemask_oceangrid,index_ice,lon_ice,lat_ice,ngrids_ice,nels_ice, + icemask,ngrids_ice,1,oceangridx,oceangridy,ngrids_ocean,options2); + delete options2; + xDelete(icemask); + + /*Put +9999 for places where there is no ice!*/ + for(int i=0;i0.) icebase_oceangrid[i]=+9999.; + xDelete(icemask_oceangrid); + + if(my_rank==0){ + ISSM_MPI_Send(icebase_oceangrid,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001008,tomitgcmcomm); + } + + /*Delete*/ + xDelete(index_ice); + xDelete(index_ocean); + xDelete(lat_ice); + xDelete(lon_ice); + xDelete(x_ice); + xDelete(y_ice); + xDelete(icebase_oceangrid); + xDelete(oceangridx); + xDelete(oceangridy); + #else + _error_("not supported"); + #endif +} Index: /issm/trunk-jpl/src/c/modules/OceanExchangeDatax/OceanExchangeDatax.h =================================================================== --- /issm/trunk-jpl/src/c/modules/OceanExchangeDatax/OceanExchangeDatax.h (revision 23714) +++ /issm/trunk-jpl/src/c/modules/OceanExchangeDatax/OceanExchangeDatax.h (revision 23714) @@ -0,0 +1,12 @@ +/*!\file: OceanExchangeDatax.h + * \brief header file for Exchange of ocean data + */ + +#ifndef _OCEANEXCHANGEDATAX_H +#define _OCEANEXCHANGEDATAX_H + +#include "../../classes/classes.h" + +/* local prototypes: */ +void OceanExchangeDatax(FemModel* femmodel); +#endif /* _OCEANEXCHANGEDATAX_H */ Index: /issm/trunk-jpl/src/c/modules/modules.h =================================================================== --- /issm/trunk-jpl/src/c/modules/modules.h (revision 23713) +++ /issm/trunk-jpl/src/c/modules/modules.h (revision 23714) @@ -100,3 +100,8 @@ #include "./UpdateDynamicConstraintsx/UpdateDynamicConstraintsx.h" #include "./VertexCoordinatesx/VertexCoordinatesx.h" + +#ifdef _HAVE_OCEAN_ +#include "./OceanExchangeDatax/OceanExchangeDatax.h" #endif + +#endif