Index: /issm/trunk-jpl/src/c/cores/transient_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/transient_core.cpp (revision 23716) +++ /issm/trunk-jpl/src/c/cores/transient_core.cpp (revision 23717) @@ -80,5 +80,5 @@ #if defined(_HAVE_OCEAN_ ) - if(isoceancoupling) OceanExchangeDatax(femmodel); + if(isoceancoupling) OceanExchangeDatax(femmodel,true); #endif @@ -126,119 +126,7 @@ femmodel->parameters->SetParam(save_results,SaveResultsEnum); - if(isoceancoupling){ /*{{{*/ - - #ifndef _HAVE_AD_ - if(VerboseSolution()) _printf0_(" ocean coupling: exchanging information\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 ngrids_ocean, nels_ocean; - IssmDouble oceantime; - IssmDouble rho_ice; - IssmDouble *oceanmelt = NULL; - IssmDouble *icebase_oceangrid = NULL; - IssmDouble *icemask_oceangrid = NULL; - IssmDouble *oceangridx = NULL; - IssmDouble *oceangridy = NULL; - IssmDouble *x_ice = NULL; - IssmDouble *y_ice = NULL; - IssmDouble *lat_ice = NULL; - IssmDouble *lon_ice = NULL; - IssmDouble *icebase = NULL; - IssmDouble *icemask = NULL; - IssmDouble *melt_mesh = NULL; - int *index_ice = NULL; - int *index_ocean = NULL; - int ngrids_ice=femmodel->vertices->NumberOfVertices(); - int nels_ice=femmodel->elements->NumberOfElements(); - - /*Recover mesh and inputs needed*/ - femmodel->parameters->FindParam(&rho_ice,MaterialsRhoIceEnum); - femmodel->GetMesh(femmodel->vertices,femmodel->elements,&x_ice,&y_ice,&index_ice); - femmodel->parameters->FindParam(&oceangridx,&ngrids_ocean,OceanGridXEnum); - femmodel->parameters->FindParam(&oceangridy,&ngrids_ocean,OceanGridYEnum); - BamgTriangulatex(&index_ocean,&nels_ocean,oceangridx,oceangridy,ngrids_ocean); - - femmodel->vertices->LatLonList(&lat_ice,&lon_ice); - - /*Interpolate ice base and mask onto ocean grid*/ - 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); - - /*Send and receive data*/ - if(my_rank==0){ - ISSM_MPI_Send(&time,1,ISSM_MPI_DOUBLE,0,10001001,tomitgcmcomm); - ISSM_MPI_Recv(&oceantime,1,ISSM_MPI_DOUBLE,0,10001002,tomitgcmcomm,&status); - if((oceantime - time > 0.1*yts) & (oceantime - time < -0.1*yts)) _error_("Ocean and ice time are starting to diverge"); - oceanmelt = xNew(ngrids_ocean); - ISSM_MPI_Recv(oceanmelt,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001007,tomitgcmcomm,&status); - ISSM_MPI_Send(icebase_oceangrid,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001008,tomitgcmcomm); - } - ISSM_MPI_Bcast(&oceantime,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - if(my_rank!=0) oceanmelt=xNew(ngrids_ocean); - ISSM_MPI_Bcast(oceanmelt,ngrids_ocean,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - - /*Interp melt onto ice grid*/ - InterpFromMeshToMesh2dx(&melt_mesh,index_ocean,oceangridx,oceangridy,ngrids_ocean,nels_ocean, - oceanmelt,ngrids_ocean,1, - lon_ice,lat_ice,ngrids_ice,NULL); - - for(int i=0;i(index_ice); - xDelete(index_ocean); - xDelete(lat_ice); - xDelete(lon_ice); - xDelete(x_ice); - xDelete(y_ice); - xDelete(melt_mesh); - xDelete(oceangridx); - xDelete(oceangridy); - xDelete(oceanmelt); - xDelete(icebase_oceangrid); - - #else - _error_("not supported"); - #endif - } - /*}}}*/ + #if defined(_HAVE_OCEAN_ ) + if(isoceancoupling) OceanExchangeDatax(femmodel,false); + #endif if(isthermal && domaintype==Domain3DEnum){ Index: /issm/trunk-jpl/src/c/modules/OceanExchangeDatax/OceanExchangeDatax.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/OceanExchangeDatax/OceanExchangeDatax.cpp (revision 23716) +++ /issm/trunk-jpl/src/c/modules/OceanExchangeDatax/OceanExchangeDatax.cpp (revision 23717) @@ -9,8 +9,8 @@ #include "./OceanExchangeDatax.h" -void OceanExchangeDatax(FemModel* femmodel){ +void OceanExchangeDatax(FemModel* femmodel, bool init_stage){ #ifndef _HAVE_AD_ - if(VerboseSolution()) _printf0_(" ocean coupling: initialization \n"); + if(VerboseSolution()) _printf0_(" ocean coupling: exchanging information\n"); int my_rank; ISSM_MPI_Comm tomitgcmcomm; @@ -23,5 +23,7 @@ int oceangridnxsize,oceangridnysize,ngrids_ocean,nels_ocean; - IssmDouble oceantime,coupling_time,time; + IssmDouble oceantime,coupling_time,time,yts; + IssmDouble rho_ice; + IssmDouble *oceanmelt = NULL; IssmDouble *oceangridx; IssmDouble *oceangridy; @@ -34,4 +36,5 @@ IssmDouble* icebase = NULL; IssmDouble* icemask = NULL; + IssmDouble* melt_mesh = NULL; int* index_ice = NULL; int* index_ocean = NULL; @@ -41,34 +44,42 @@ /*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); + + /*Exchange or recover mesh and inputs needed*/ + if(init_stage==true){ + 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); } - 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); + else{ + femmodel->parameters->FindParam(&oceangridx,&ngrids_ocean,OceanGridXEnum); + femmodel->parameters->FindParam(&oceangridy,&ngrids_ocean,OceanGridYEnum); } - 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*/ @@ -87,6 +98,5 @@ 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); + icebase,ngrids_ice,1,oceangridx,oceangridy,ngrids_ocean,options); delete options; xDelete(icebase); @@ -111,6 +121,31 @@ xDelete(icemask_oceangrid); - if(my_rank==0){ - ISSM_MPI_Send(icebase_oceangrid,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001008,tomitgcmcomm); + if(init_stage==true){ //just send icebase + if(my_rank==0){ + ISSM_MPI_Send(icebase_oceangrid,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001008,tomitgcmcomm); + } + } + else{ //send and receive exchanged data + femmodel->parameters->FindParam(&rho_ice,MaterialsRhoIceEnum); + femmodel->parameters->FindParam(&yts,ConstantsYtsEnum); + if(my_rank==0){ + ISSM_MPI_Send(&time,1,ISSM_MPI_DOUBLE,0,10001001,tomitgcmcomm); + ISSM_MPI_Recv(&oceantime,1,ISSM_MPI_DOUBLE,0,10001002,tomitgcmcomm,&status); + if((oceantime - time > 0.1*yts) & (oceantime - time < -0.1*yts)) _error_("Ocean and ice time are starting to diverge"); + oceanmelt = xNew(ngrids_ocean); + ISSM_MPI_Recv(oceanmelt,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001007,tomitgcmcomm,&status); + ISSM_MPI_Send(icebase_oceangrid,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001008,tomitgcmcomm); + } + ISSM_MPI_Bcast(&oceantime,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + if(my_rank!=0) oceanmelt=xNew(ngrids_ocean); + ISSM_MPI_Bcast(oceanmelt,ngrids_ocean,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + + /*Interp melt onto ice grid*/ + InterpFromMeshToMesh2dx(&melt_mesh,index_ocean,oceangridx,oceangridy,ngrids_ocean,nels_ocean, + oceanmelt,ngrids_ocean,1, + lon_ice,lat_ice,ngrids_ice,NULL); + + for(int i=0;i(oceangridx); xDelete(oceangridy); + xDelete(melt_mesh); + xDelete(oceanmelt); #else _error_("not supported"); Index: /issm/trunk-jpl/src/c/modules/OceanExchangeDatax/OceanExchangeDatax.h =================================================================== --- /issm/trunk-jpl/src/c/modules/OceanExchangeDatax/OceanExchangeDatax.h (revision 23716) +++ /issm/trunk-jpl/src/c/modules/OceanExchangeDatax/OceanExchangeDatax.h (revision 23717) @@ -9,4 +9,4 @@ /* local prototypes: */ -void OceanExchangeDatax(FemModel* femmodel); +void OceanExchangeDatax(FemModel* femmodel, bool init_stage); #endif /* _OCEANEXCHANGEDATAX_H */